How should I compute the eigenvectors of a sparse, real, symmetric matrix?
34 views (last 30 days)
Show older comments
I need to compute all the eigenvectors of a sparse, real, symmetric matrix, A.
I've tried the following:
>> [V, D] = eig(A);
Error using eig
For standard eigenproblem EIG(A) when A is sparse, EIG does not support computing eigenvectors. Use EIGS instead.
>> [V, D] = eigs(A, size(A, 1));
Warning: For real symmetric problems, must have number of eigenvalues k < n.
Using eig instead.
So MATLAB gives me conflicting advice: first, to use eigs, and then to use eig! What it does internally is convert my sparse matrix to a full one. This seems inefficient if I have lots of sparsity. Any better suggestions?
2 Comments
Accepted Answer
More Answers (2)
Matt J
on 1 Nov 2013
Edited: Matt J
on 1 Nov 2013
For the same reasons it has an advantage in other computations.
Sparsity just isn't always helpful. You have to know if/why it will be useful in a given situation.
In the case of eigenvalue decomposition, it's hard to see how sparsity could be exploited. Eigendecomposition is based on QR decomposition and the QR decompositions of sparse matrices are not sparse. That could be the reason why qr() is so much slower for sparse matrices, e.g.,
>> As=sprand(3e3,3e3,1e-3); Af=full(As);
>> tic; [Q,R]=qr(Af);toc
Elapsed time is 2.506753 seconds.
>> tic; [Q,R]=qr(As);toc
Elapsed time is 15.731333 seconds.
There are also operations that are slower for MATLAB's sparse type because of software design compromises, even simple transposition:
>> Af=zeros(1e7,1); As=sparse(Af);
>> tic; Af.'; toc
Elapsed time is 0.000012 seconds.
>> tic; As.'; toc
Elapsed time is 0.054525 seconds.
5 Comments
Matt J
on 1 Nov 2013
Edited: Matt J
on 1 Nov 2013
But do you not get my perspective here?
It's an understandable perspective from someone convinced that there was an "implication of fallacy" in my initial comment. However, you could also consider the possibility that I in fact thought you might be right, but that your reasoning was not obvious, and that I was just looking for elaboration...
Why my answer is considered off point, I still don't get. I explained what you eventually concluded yourself in your Answer, that "sparse isn't always faster".
I also don't understand why you decided your question was resolved and closed the thread. Your original question was if/why it's not possible to do better than eig(full(S)). Do you now believe you know why? If so, why is my theory about QR so outlandish? I think we agree that eig() uses QR. If you're now convinced that eig(), and hence QR, is the optimal eigendecomposition method, isn't it relevant that QR won't be faster in sparse form?
Matt J
on 3 Nov 2013
Edited: Matt J
on 3 Nov 2013
It does appear that the Lanczos algorithm, if that's what is being used by EIGS, is not being used optimally for the case when only the eigenvalues are requested as output. According to here, Lanczos should be able to derive the eigenvalues in O(N^2) for a sparse matrix of density 1/N. However, in my tests below, computation time for the eigenvalues does seem to go cubically with N,
for N=[1:4]*1e3;
A = sprand(N, N, 1/N); A=A+A.';
tic; E = eigs(A, N-1); toc;
end
Elapsed time is 1.124703 seconds.
Elapsed time is 8.516908 seconds.
Elapsed time is 27.762566 seconds.
Elapsed time is 65.062230 seconds.
It might therefore be worth trying some of the external MATLAB Lanczos implementations, also at the link above.
When eigenvectors are requested as well, all bets are off. There is nothing I can see in the Lanczos algorithm that gaurantees less than O(N^3) computation. In particular, the final conversion from Lanczos vectors to eigenvectors would require an NxN matrix multiplication. Also, non-sparse memory allocation is required just to hold the N eigenvectors, so who kows what overhead that will bring. It doesn't appear that Lanczos was meant for doing full eigenvector decomposition, only partial ones.
0 Comments
See Also
Categories
Find more on Sparse Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!