scipy generalized eigenproblem with positive semidefinite

Question

Hi, guys!!!

I want to compute generalized eigendecomposition of the form:

Lf = lambda Af

by using scipy.sparse.linalg.eigs function, but get this error:

/usr/local/lib/python2.7/dist-packages/scipy/linalg/decomp_lu.py:61: RuntimeWarning: Diagonal number 65 is exactly zero. Singular matrix. RuntimeWarning) ** On entry to DLASCL parameter number 4 had an illegal value

I am passing three arguments, a diagonal matrix, a positive semi-definite (PSD) matrix and numeric value K (first K eigenvalues). Matlab's eigs function performs well using the same input parameters, but in SciPy as I have understood, in order to compute with PSD I need to specify sigma parameter as well.

So, my question is: is there a way to avoid setting sigma parameter, as it is in MatLab, or if not, how to pick up sigma value?

Looking forward to getting advices or hints... Thank you in advance!

Answer 1

The error appears to mean that in your generalized eigenproblem

L x = lambda A x

the matrix A is not positive definite (check the eigs docstring -- in your case the matrix is probably singular). This is a requirement for ARPACK mode 2. However, you can try specifying sigma=0 to switch to ARPACK mode 3 (but note that the meaning of the which parameter is inverted in this case!).

Now, I'm not sure what Matlab does, but a possibility is that it's calculating the pseudoinverse rather than the inverse of A. To emulate this, do

from scipy.sparse.linalg import LinearOperator
from scipy.linalg import lstsq

Ainv = LinearOperator(matvec=lambda x: lstsq(A, x)[0], shape=A.shape)
w, v = eigs(L, M=A, Minv=Ainv)

Check the results --- I don't know what will happen in this case.

Alternatively, you may try to specify a nonzero sigma. What you should select depends on the matrices involved. It affects the eigenvalues that are picked --- for instance with which='LM' are those for which lambda' = 1/(lambda - sigma) is large. Otherwise, it can probably be chosen arbitrarily, of course it's probably better for the Krylov progress if the transformed eigenvalues lambda' which you are interested in become well separated from the other eigenvalues.