Search code examples
pythonscipysparse-matrixeigenvalue

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

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?

Solution

  • 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