Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #7244 > unrolled thread

Eigensolver for Large Sparse Matrices in Python

Back to article view | Back to comp.lang.python

Contents

  Eigensolver for Large Sparse Matrices in Python Andrew MacLean <andrew.maclean@gmail.com> - 2011-06-08 08:44 -0700
    Re: Eigensolver for Large Sparse Matrices in Python Javier <nospam@nospam.com> - 2011-06-09 19:59 +0000

#7244 — Eigensolver for Large Sparse Matrices in Python

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

[toc] | [next] | [standalone]

#7315

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@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

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web