[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: |
---------------------------------+------------------------------------------
Hi,
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''
library.
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:
[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 for X, and ''lobpcg'' returned the correct
values.
I am wondering if anyone else has encountered problems using matrices of
this size with ''lobpcg''.
Thanks for any suggestions on this,
Andrew
--
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