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
