[Scipy-tickets] [SciPy] #1456: scipy.sparse.linalg.lobpcg fails for large sparse matrices

SciPy Trac scipy-tickets@scipy....
Wed Jun 8 09:52:24 CDT 2011

#1456: scipy.sparse.linalg.lobpcg fails for large sparse matrices
 Reporter:  AndrewMacLean        |       Owner:  wnbell
     Type:  defect               |      Status:  new   
 Priority:  normal               |   Milestone:        
Component:  scipy.sparse.linalg  |     Version:  0.7.2 
 Keywords:                       |  

 I have been using ''lobpcg'' (scipy.sparse.linalg.lobpcg) to solve
 symmetric generalized eigenvalue problems with large, sparse stiffness and
 mass matrices, say '''A''' and '''B'''. The problem is of the form
 '''Av''' =lambda'''BV'''.

 They are both symmetric, and '''B''' is positive definite, and both are 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''

 The part of my code using lobpcg to find 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:


 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.


 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 for X, and ''lobpcg'' returned the correct

 I am wondering if anyone else has encountered problems using matrices of
 this size with ''lobpcg''.

 Thanks for any suggestions on this,

Ticket URL: <http://projects.scipy.org/scipy/ticket/1456>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.

More information about the Scipy-tickets mailing list