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 |
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 |
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 |
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 |
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 -- 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 |
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. > > > > 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 |
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 |
2008/4/15, Ondrej Certik <[hidden email]>: On Tue, Apr 15, 2008 at 4:29 PM, Matthieu Brucher 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 |
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 |
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 |
Free forum by Nabble | Edit this page |