Eigensolver for Large Sparse Matrices in Python

[Andrew MacLean]

Hi,

I need to solve symmetric generalized eigenvalue problems with large,
sparse stiffness
and mass matrices, say 'A' and 'B'. The problem is of the form Av =
lambdaBV. I have been using lobpcg (scipy.sparse.linalg.lobpcg), in
Scipy 0.7.2, although this has been giving me incorrect values that
are also about 1000 times too large.

They are both Hermitian, and 'B' is positive definite, and both are of
approximately of size 70 000 x 70 000. 'A' has about 5 million non-
zero values and 'B'
has about 1.6 million. 'A' has condition number on the order of 10^9
and 'B' has a condition number on the order of 10^6. I have stored
them both as "csc" type sparse matrices from the scipy.sparse library.

The part of my code using lobpcg is fairly simple (for the 20 smallest
eigenvalues):
--------------------------------------------------------------------------------------------------------
from scipy.sparse.linalg import lobpcg
from scipy import rand

X = rand(A.shape[0], 20)

W, V = lobpcg (A, X, B = B, largest = False, maxiter = 40)
-------------------------------------------------------------------------------------------------------

I tested lobpcg on a "scaled down" version of my problem, with 'A' and
'B' matrices of size 10 000 x 10 000, and got the correct values
(using maxiter = 20), as confirmed by using the "eigs" function in
Matlab. I used it here to find the smallest 10 eigenvalues, and here
is a table of my results, showing that the eigenvalues computed from
lobpcg in Python are very close to those using eigs in Matlab:

https://docs.google.com/leaf?id=0Bz-X2kbPhoUFMTQ0MzM2MGMtNjgwZi00N2U0...

With full sized 'A' and 'B' matrices, I could not get the correct
values, and it became clear that increasing the maximum number of
iterations indefinitely would not work (see graph below). I made a
graph for the 20th smallest eigenvalue vs. the number of iterations. I
compared 4 different guesses for X, 3 of which were just random
matrices (as in the code above), and a 4th orthonormalized one.

https://docs.google.com/leaf?id=0Bz-X2kbPhoUFYTM4OTIxZGQtZmE0Yi00MTMy...

It appears that it will take a very large number of iterations to get
the correct eigenvalues.  As well, I tested lobpcg by using
eigenvectors generated by eigs in
Matlab as X, and lobpcg returned the correct values.

I don't believe it is a bug that was solved for lobpcg in newer
versions of Scipy, as I have also gotten the same problem using the
most recent version (4.12) of lobpcg for Matlab.

If anyone has any suggestions for how to improve the results of
lobpcg, or has experience with an alternate solver (such as JDSYM from
Pysparse or eigsh in newer versions of Scipy) with matrices of this
size, any recommendations would be grealty appreciated.

Thanks,
Andrew