Eigensolver for Large Sparse Matrices in Python

[Javier]

Hi,

I think you can also use scipy.sparse.linalg.eigen.arpack in addition to 
scipy.sparse.linalg.eigen.lobpcg

Also, from my experience with this routines I can tell you that they
don't like to be asked a small number of eigenvalues.
Contrary to common sense I have found these routines to prefer to
calculate 30 eigenvalues than 5 eigenvalues.  Try to ask it for more
eigenvalues.

Does anybody know why the routines don't work well when they are aked
for  small number of eigenvalues?

Andrew MacLean <andrew.maclean at gmail.com> wrote:
> 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