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.macl...@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 -- http://mail.python.org/mailman/listinfo/python-list