generalized eigenvalue problem for sparse matrices

classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

generalized eigenvalue problem for sparse matrices

Abhinav Sarkar
Hi
  I am trying to solve a generalized eigenvalue problem for sparse
 matrices. The problem is in form
    A*x = 𝜎*M*x
 where A and M are sparse matrices and x is a vector and sigma is a
 scalar whose value is to be found.

 For this I am using the Arpack function ARPACK_gen_eigs provided in
 the module scipy.sparse.linalg.eigen.arpack.speigs . However the
 solution which I am getting does not match with the solution obtained
 from the eig function in MATLAB. For example:

 A = [1 0 0 -33 16 0; 0 1 0 16 -33 16; 0 0 1 0 16 -33; 1601 -96 256
 -1800 0 0; -96 1601 -96 0 -1800 0; 256 -96 1601 0 0 -1800]
 M = [0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; -33 16 0 0 0 0; 16 -33 16
 0 0 0; 0 16 -33 0 0 0]

 The matrices are written in MATLAB format, spaces separating the row
 elements and ";" separating the rows.
 When I solve this problem in MATLAB using the function eig, I get the
 following eigenvalues:
 -155.0345, -56.9898, -45.2209, -31.9376, -28.5367, -9.1611

 However when I solve it using ARPACK_gen_eigs to find the two
 eigenvalues near 10, I get the following solution:
 13.83675153-1.71439075j, 13.83675153+1.71439075j
 which is certainly not correct.

 I am using the following code to solve the problem:


 from scipy import sparse
 from scipy.sparse.linalg import dsolve
 from scipy.sparse.linalg import eigen
 from scipy.sparse.linalg.eigen.arpack.speigs import ARPACK_gen_eigs

 n = 3
 h = 1.0/(n+1)

 a = 1
 Pr = 7.0
 Ra = 1000

 A = get_A_mat(n, h, a, Pr, Ra)
 M = get_M_mat(n, h, a, Pr, Ra)

 sigma = 10.0
 B = A - sigma*M
 s = dsolve.splu(B)
 e, v = ARPACK_gen_eigs(M.matvec, s.solve, 2*n, sigma, 2, 'LR')
 print e


 which also seems to be correct to me. Please tell if the method I am
 using is correct or not and why am I not getting correct solutions. Is
 the ARPACK_gen_eigs broken? Or is there a problem in my code?

 Regards
 --
 Abhinav Sarkar
 4th year Undergraduate Student
 Deptt. of Mechanical Engg.
 Indian Institute of Technology, Kharagpur
 India

 Web: http://claimid.com/abhin4v
 Twitter: http://twitter.com/abhin4v
 ---------
 The world is a book, those who do not travel read only one page.
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Abhinav Sarkar
abhinav sarkar <abhinav.sarkar <at> gmail.com> writes:

>
> Hi
>   I am trying to solve a generalized eigenvalue problem for sparse
>  matrices. The problem is in form
>     A*x = 𝜎*M*x
>  where A and M are sparse matrices and x is a vector and sigma is a
>  scalar whose value is to be found.
>
>  For this I am using the Arpack function ARPACK_gen_eigs provided in
>  the module scipy.sparse.linalg.eigen.arpack.speigs . However the
>  solution which I am getting does not match with the solution obtained
>  from the eig function in MATLAB. For example:
>
>  A = [1 0 0 -33 16 0; 0 1 0 16 -33 16; 0 0 1 0 16 -33; 1601 -96 256
>  -1800 0 0; -96 1601 -96 0 -1800 0; 256 -96 1601 0 0 -1800]
>  M = [0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; -33 16 0 0 0 0; 16 -33 16
>  0 0 0; 0 16 -33 0 0 0]
>
>  The matrices are written in MATLAB format, spaces separating the row
>  elements and ";" separating the rows.
>  When I solve this problem in MATLAB using the function eig, I get the
>  following eigenvalues:
>  -155.0345, -56.9898, -45.2209, -31.9376, -28.5367, -9.1611
>
>  However when I solve it using ARPACK_gen_eigs to find the two
>  eigenvalues near 10, I get the following solution:
>  13.83675153-1.71439075j, 13.83675153+1.71439075j
>  which is certainly not correct.
>  --
>  Abhinav Sarkar


I tried the eigenvalue problem solver provided for the dense matrices at
scipy.linalg.eig and it gives the same result as the MATLAB functin eig. Then
why is scipy.sparse.linalg.eigen.arpack.speigs.ARPACK_gen_eigs giving different
result? Please help me out.

Regards
---
Abhinav Sarkar

_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Neilen Marais-4
In reply to this post by Abhinav Sarkar
Abinav,

On Tue, 15 Apr 2008 01:50:20 +0530, abhinav sarkar wrote:

> Hi
>  n = 3
>  h = 1.0/(n+1)
>
>  a = 1
>  Pr = 7.0
>  Ra = 1000
>
>  A = get_A_mat(n, h, a, Pr, Ra)
>  M = get_M_mat(n, h, a, Pr, Ra)
>
>  sigma = 10.0
^^^
I think you're looking for -10.

>  B = A - sigma*M
>  s = dsolve.splu(B)
>  e, v = ARPACK_gen_eigs(M.matvec, s.solve, 2*n, sigma, 2, 'LR')
                                                             ^^
You may try futzing around with this parameter.

Anyway, I did play around a bit with your example and never really got a
good answer. But my experience has been that ARPACK doesn't really work
well with small matrices. Try a bigger problem and see if things turn out
better, I mean, obviously you don't need a sparse solver to solve a 6x6
matrix system ;). You may also find that if you force Matlab to do this
problem using a sparse solver you'll run into the same issues.
Unfortunately I don't know enough about the ARPACK routines to give more
insight. They do work quite well on my problems, where the eigenvalues
are always positive.

>  which also seems to be correct to me. Please tell if the method I am
>  using is correct or not and why am I not getting correct solutions. Is
>  the ARPACK_gen_eigs broken? Or is there a problem in my code?

Seems OK, but I'm no expert ;)

Regards
Neilen

>
>  Regards
>  --
>  Abhinav Sarkar
>  4th year Undergraduate Student
>  Deptt. of Mechanical Engg.
>  Indian Institute of Technology, Kharagpur India
>
>  Web: http://claimid.com/abhin4v
>  Twitter: http://twitter.com/abhin4v
>  ---------
>  The world is a book, those who do not travel read only one page.
> _______________________________________________ SciPy-user mailing list
> [hidden email]
> http://projects.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Nathan Bell-4
In reply to this post by Abhinav Sarkar
On Tue, Apr 15, 2008 at 1:48 AM, Abhinav Sarkar
<[hidden email]> wrote:
>  I tried the eigenvalue problem solver provided for the dense matrices at
>  scipy.linalg.eig and it gives the same result as the MATLAB functin eig. Then
>  why is scipy.sparse.linalg.eigen.arpack.speigs.ARPACK_gen_eigs giving different
>  result? Please help me out.

Have you also tried eigs() in MATLAB?  I'm not sure what happens when
you call eig() with a sparse matrix.

As Neilen said, this isn't a good test for ARPACK.  Try setting the
'tol' parameter to something smaller.

--
Nathan Bell [hidden email]
http://graphics.cs.uiuc.edu/~wnbell/
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Matthieu Brucher-2
Hi,

I've tested this function with a 2000x2000 array, and like Abhinav, I didn't manage to find correct eigenvalues. Everything worked as expected with a dense array with the same element with numpy.linalg.eig, but not with Arpack...

Matthieu

2008/4/15, Nathan Bell <[hidden email]>:
On Tue, Apr 15, 2008 at 1:48 AM, Abhinav Sarkar
<[hidden email]> wrote:
>  I tried the eigenvalue problem solver provided for the dense matrices at
>  scipy.linalg.eig and it gives the same result as the MATLAB functin eig. Then
>  why is scipy.sparse.linalg.eigen.arpack.speigs.ARPACK_gen_eigs giving different
>  result? Please help me out.


Have you also tried eigs() in MATLAB?  I'm not sure what happens when
you call eig() with a sparse matrix.

As Neilen said, this isn't a good test for ARPACK.  Try setting the
'tol' parameter to something smaller.


--
Nathan Bell [hidden email]
http://graphics.cs.uiuc.edu/~wnbell/

_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user



--
French PhD student
Website : http://matthieu-brucher.developpez.com/
Blogs : http://matt.eifelle.com and http://blog.developpez.com/?blog=92
LinkedIn : http://www.linkedin.com/in/matthieubrucher
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Abhinav Sarkar
In reply to this post by Neilen Marais-4
On Tue, Apr 15, 2008 at 6:24 PM, Neilen Marais <[hidden email]> wrote:

> Abinav,
>
>  On Tue, 15 Apr 2008 01:50:20 +0530, abhinav sarkar wrote:
>
>  > Hi
>
> >  n = 3
>  >  h = 1.0/(n+1)
>  >
>  >  a = 1
>  >  Pr = 7.0
>  >  Ra = 1000
>  >
>  >  A = get_A_mat(n, h, a, Pr, Ra)
>  >  M = get_M_mat(n, h, a, Pr, Ra)
>  >
>  >  sigma = 10.0
>  ^^^
>  I think you're looking for -10.
i am looking for the largest eigenvalue

>
>
>  >  B = A - sigma*M
>  >  s = dsolve.splu(B)
>  >  e, v = ARPACK_gen_eigs(M.matvec, s.solve, 2*n, sigma, 2, 'LR')
>                                                              ^^
>  You may try futzing around with this parameter.
>
>  Anyway, I did play around a bit with your example and never really got a
>  good answer. But my experience has been that ARPACK doesn't really work
>  well with small matrices. Try a bigger problem and see if things turn out
>  better, I mean, obviously you don't need a sparse solver to solve a 6x6
>  matrix system ;). You may also find that if you force Matlab to do this
>  problem using a sparse solver you'll run into the same issues.
>  Unfortunately I don't know enough about the ARPACK routines to give more
>  insight. They do work quite well on my problems, where the eigenvalues
>  are always positive.
>
>
>  >  which also seems to be correct to me. Please tell if the method I am
>  >  using is correct or not and why am I not getting correct solutions. Is
>  >  the ARPACK_gen_eigs broken? Or is there a problem in my code?
>
>  Seems OK, but I'm no expert ;)
>
>  Regards
>  Neilen
>
>
>  >
>  >  Regards
>  >  --
>  >  Abhinav Sarkar
>  >  4th year Undergraduate Student
>  >  Deptt. of Mechanical Engg.
>  >  Indian Institute of Technology, Kharagpur India
>  >
>  >  Web: http://claimid.com/abhin4v
>  >  Twitter: http://twitter.com/abhin4v
>  >  ---------
>  >  The world is a book, those who do not travel read only one page.
>  > _______________________________________________ SciPy-user mailing list
>
>
> > [hidden email]
>  > http://projects.scipy.org/mailman/listinfo/scipy-user
>
>
>  _______________________________________________
>  SciPy-user mailing list
>  [hidden email]
>  http://projects.scipy.org/mailman/listinfo/scipy-user
>

 I have formulated the problem in a slightly different way and now I
am able to obtain the correct eigenvalues. Here is the new code:

class Solver:
        def __init__(self, A, M, sigma):
                self.A = A
                self.M = M
                self.sigma = sigma
               
        def solve(self, b):
                #(M-sigma*A)*x=A*b
                #return x
                return dsolve.spsolve((self.M-self.sigma*self.A), self.A*b)

A = get_A_mat(n, h, a, Pr, Ra)
M = get_M_mat(n, h, a, Pr, Ra)
I = sparse.eye(2*n,2*n)
sigma = 0.0
sigma_solve = Solver(A, M, sigma).solve
e = 1.0/ARPACK_gen_eigs(I.matvec, sigma_solve, I.shape[0], sigma, 5, 'LR')[0]

Here I have converted the A*x = sigma*M*x to inv(A)*M*x = nu*I*x where
nu = 1/sigma. And replaced s.solve by my own Solver class. Now
solutions obtained are correct.
--
Abhinav Sarkar
4th year Undergraduate Student
Deptt. of Mechanical Engg.
Indian Institute of Technology, Kharagpur
India

Mobile:+91-99327-41517
Residence:+91-631-2429887
Web: http://claimid.com/abhin4v
Twitter: http://twitter.com/abhin4v
---------
The world is a book, those who do not travel read only one page.
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Ondrej Certik
In reply to this post by Matthieu Brucher-2
On Tue, Apr 15, 2008 at 4:29 PM, Matthieu Brucher
<[hidden email]> wrote:
> Hi,
>
> I've tested this function with a 2000x2000 array, and like Abhinav, I didn't
> manage to find correct eigenvalues. Everything worked as expected with a
> dense array with the same element with numpy.linalg.eig, but not with
> Arpack...


When I compiled arpack myself and used it myself, I also didn't manage
to get correct results. Anyway, I wouldn't recommend to use arpack.

Ondrej
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Matthieu Brucher-2


2008/4/15, Ondrej Certik <[hidden email]>:
On Tue, Apr 15, 2008 at 4:29 PM, Matthieu Brucher
<[hidden email]> wrote:
> Hi,
>
> I've tested this function with a 2000x2000 array, and like Abhinav, I didn't
> manage to find correct eigenvalues. Everything worked as expected with a
> dense array with the same element with numpy.linalg.eig, but not with
> Arpack...



When I compiled arpack myself and used it myself, I also didn't manage
to get correct results. Anyway, I wouldn't recommend to use arpack.


Ondrej

What do you recommend for generalized eigenproblems with sparse arrays ?

Matthieu
--
French PhD student
Website : http://matthieu-brucher.developpez.com/
Blogs : http://matt.eifelle.com and http://blog.developpez.com/?blog=92
LinkedIn : http://www.linkedin.com/in/matthieubrucher
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Ondrej Certik
On Tue, Apr 15, 2008 at 4:55 PM, Matthieu Brucher
<[hidden email]> wrote:

>
>
> 2008/4/15, Ondrej Certik <[hidden email]>:
>
> > On Tue, Apr 15, 2008 at 4:29 PM, Matthieu Brucher
> > <[hidden email]> wrote:
> > > Hi,
> > >
> > > I've tested this function with a 2000x2000 array, and like Abhinav, I
> didn't
> > > manage to find correct eigenvalues. Everything worked as expected with a
> > > dense array with the same element with numpy.linalg.eig, but not with
> > > Arpack...
> >
> >
> >
> > When I compiled arpack myself and used it myself, I also didn't manage
> > to get correct results. Anyway, I wouldn't recommend to use arpack.
> >
> >
> > Ondrej
> >
>
> What do you recommend for generalized eigenproblems with sparse arrays ?

I have a very good experience with pysparse. That's why I was
objecting to it's removal from scipy. But I think I got ok to put it
back in, I just didn't find time for that.

Plus there are other good solvers like bzlpack or primme, both should
be able to do generalized eigenvalues and both are more modern than
arpack, but I haven't yet tried them extensively, but it's on my todo
list to create a scikit with easy to use interface to all of them. Any
help appreciated. :)

Ondrej
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: generalized eigenvalue problem for sparse matrices

Neilen Marais-4
In reply to this post by Abhinav Sarkar
Abihnav,

On Tue, 15 Apr 2008 20:13:38 +0530, abhinav sarkar wrote:

> Here I have converted the A*x = sigma*M*x to inv(A)*M*x = nu*I*x where
> nu = 1/sigma. And replaced s.solve by my own Solver class. Now solutions
> obtained are correct.

I'm glad you found a method to make it work. Not sure why the obvious
solution didn't work, but at least anyone else who needs to solve a
similar problem can use your method for now :)

Regards
Neilen

_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user