# how to diagonalize a large sparse symmetric matrix, to get the eigenvalues and eigenvectors

How does one diagonalize a large sparse symmetric matrix to get the eigenvalues and the eigenvectors?

The problem is the matrix could be very large (though it is sparse), at most 2500*2500. Is there a good algorithm to do that and most importantly, I can implement it into my own code. Thanks a lot!

-
Why would you want to implement it on your own? Matlab can do it just fine. –  Yuval Filmus Nov 28 '10 at 5:15
I am somewhat interested in this question because I know nothing about algorithmic efficiency. My naive thoughts here are that the usual diagonalization algorithm (i.e., performing simultaneous row and column operations) should go faster the sparser the matrix is. From a practical standpoint, it would be useful to have a "sparse matrix" datatype, so that the computer knows from the beginning that most of the row-column operations do not need to be performed. But is there more to it than this? I.e., does one actually use a different algorithm rather than just different implementation? –  Pete L. Clark Nov 28 '10 at 8:05
1. Is there any structure in your sparse matrix? The efficiency of a sparse eigensolver ultimately rests on how well you wrote your matrix-vector multiplication routine, since Lanczos/Arnoldi requires at its core the multiplication of your sparse matrix with a number of vectors per iteration. –  J. M. is back. Nov 28 '10 at 12:33
2. Do you really need all the eigenvalues and eigenvectors? For most applications of sparse eigenproblems, one only needs the largest few and/or the smallest few. In addition to that, your 2500-by-2500 matrix of eigenvectors is guaranteed not to be sparse; so unless you have provisions for storing all 2500 vectors, as well as the time needed to wait for them (Lanczos/Arnoldi is fast for one vector, but for 2500 vectors...), you might want to reconsider your wants/needs. –  J. M. is back. Nov 28 '10 at 12:37
@Pete: It might interest you to know that sparse matrix storage techniques borrow a lot from graph theory, but not being a graph theory expert, that's all I can say about it. For tridiagonalizing a symmetric matrix with a neat pattern, I suppose judicious use of rotations might work, but if one tries Jacobi's scheme for diagonalization on a sparse matrix, you get a lot of fill-in in the first few iterations (even though theoretically the off-diagonal elements should converge to zero in the limit). –  J. M. is back. Nov 28 '10 at 15:22

$2500 \times 2500$ is a small matrix by current standards. The standard eig command of matlab should be able to handle this size with ease. Iterative sparse matrix eigensolvers like those implemented in ARPACK, or SLEPc will become more preferable if the matrix is much larger.

Also, if you want to implement an eigensolver into your own code, just use the LAPACK library that comes with very well developed routines for such purpose. Matlab also ultimately invokes LAPACK routines for doing most of its numerical linear algebra.

Semi-related note: the matrix need not be explicitly available for the large sparse solvers, because they usually just depend on being able to compute $A*x$ and $A'*x$.

-
Sure, it's relatively small, but sometimes the need to exploit sparsity is due to time (QR or MRRR or divide-and-conquer take longer than Lanczos/Arnoldi if you need only a few eigenvalues/eigenvectors.) –  J. M. is back. Nov 28 '10 at 12:39
yes, of course: but the OP seems to want all eigenvectors and eigenvalues, then probably viewing it as a sparse matrix might not be that great! Also, for this size of a matrix, having multicore or GPU implementations of dense linear algebra are significantly easier (or at least have better scaling) than the sparse case. –  user1709 Nov 28 '10 at 13:28
Hence my questions to him/her; I am rather doubtful that s/he needs all of them! Still, I would have to agree that LAPACK/BLAS is pretty well-tuned as it stands! –  J. M. is back. Nov 28 '10 at 14:35
yes, I think your question to him/her regarding whether he/she really needs all eigenvals+vecs is a very important one. one setting where one may require (in the worst case) all might be when having to repeatedly ensure positive-definiteness of a matrix during an iterative process. –  user1709 Nov 28 '10 at 15:38