python scipy.sparse.linalg.eigs giving different results for consecutive calls

Question

I am trying to compute the spectral radius of a sparse matrix in python. This is what I have:

from scipy.sparse.linalg import eigs
from scipy import sparse

w = sparse.rand(10, 10, 0.1)
spec_radius = max(abs(eigs(w)[0]))

where the values of w get scaled to be in the range of [-1,1]. However, running that command gives a different result every time:

>>> print max(abs(eigs(w)[0]))
4.51859016293e-05

>>> print max(abs(eigs(w)[0]))
4.02309443625e-06

>>> print max(abs(eigs(w)[0]))
3.7611221426e-05

What gives? I would have thought it would be the same every time. Am I misunderstanding how these commands work?

Answer 1

You apply the methods correctly and they will give you the same results if the absolute value of the largest eigenvalue is significantly larger than 0. The reason for the outcome you observe is based on the iterative nature of the algorithm that is used to determine the eigenvalues. From the documentation:

"This function is a wrapper to the ARPACK [R209] SNEUPD, DNEUPD, CNEUPD, ZNEUPD, functions which use the Implicitly Restarted Arnoldi Method to find the eigenvalues and eigenvectors [R210]." In case you are interested in details of this algorithm, you can find an explanation e.g. here.

As in all numerical methods you can determine the desired value only with a certain precision. For eigenvalues that are significantly unequal to 0, you will always obtain the same output; for values that are close to 0 you might obtain different values as you observed in the above example.

You can try to vary the parameters "maxiter" and "tol" (check the above cited documentation for details) which you can pass to "eigs". Maxiter is the maximum number of Arnoldi update iterations allowed - by increasing the number you should get more accurate results. "Tol" is the relative accuracy for eigenvalues and stopping criterion of the algorithm.

Answer 2

Sorry for responding to an old question here, but the other answer is not quite satisfactory.

The randomness is not part of the algorithm bundled with ARPACK, but rather the initialization of the algorithm. From the scipy documentation the initialization, v0, is random, unless specified by the user. Sure enough, we see this (note the setup is slightly different--entries of w are scaled to be in [0,1]):

import numpy
from scipy.sparse.linalg import eigs
from scipy import sparse
w = sparse.rand(10, 10, 0.1)
w = w/w.max()

if we do not specify v0, we get some (slight) randomness:

>>> print max(abs(eigs(w)[0]))
0.00024188777676476916
>>> print max(abs(eigs(w)[0]))
0.00028073646868200566
>>> print max(abs(eigs(w)[0]))
0.00025250058038424729
>>> print max(abs(eigs(w)[0]))
0.00018183677959035711

but, if we specify, the initialization, we always get the same answer:

>>> print numpy.all([max(abs(eigs(w, v0 = numpy.ones(10))[0])) == 0.00026363015600771211 for k in range(1000)])
True