How should I compute the eigenvectors of a sparse, real, symmetric matrix?

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?

[Matt J]

What it does internally is convert my sparse matrix to a full one.This seems inefficient if I have lots of sparsity.

Why do you think sparsity should have an advantage in eigenvalue computation?

[Oliver Woodford]

For the same reasons it has an advantage in other computations.

Answer 1

[Oliver Woodford]

Just convert the matrix to a full one, and use eig. Unless less you want just a few eigenvectors, then the decomposition using the sparse matrix will generally be slower anyway. E.g.:

>> A = sprand(2e3, 2e3, 2e-3);
>> A = A' + A;
>> tic; [V, D] = eigs(A, 2e3-2); toc;
Elapsed time is 21.727999 seconds.
>> tic; [V, D] = eig(full(A)); toc
Elapsed time is 2.488241 seconds.

Moral of the story: sparse isn't always faster!

Answer 2

[Matt J]

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.

[Oliver Woodford]

It would be better to give the timings for eig and eigs - that's what I asked about. I understand the argument you were trying to make, but QR decomposition isn't even used in eigs on symmetric, sparse matrices, the Lanczos algorithm is, and the transposition timings seem irrelevant for a number of reasons (too many to write here!).

A comment on your style here: I'm sure you intended to help with this answer, but I didn't appreciate it much. Here's why: first you asked me to defend my question, then you smacked down my answer with wording like "You have to know if/why it will be useful" and "it's hard to see how sparsity could be exploited" - I found this particularly condescending in tone. It would have been less irritating if the answer were obvious, or the argument even correct, neither of which appear to me to be the case, but I don't think it's ever a good idea.

This is a question and answer site - we should encourage questions, not discourage them by insisting people justify them, and we should give straightforward rather than lecturing answers. I hope I'm being constructive here - we should encourage answers too, after all. I have given what I think is a constructive and helpful answer above. Am I justified here, or am I the only person who is annoyed by such a response to a question posed in good faith?

[Matt J]

Oliver, I'm sure you didn't mean to be condescending either, when you said "For the same reasons it has an advantage in other computations" and nothing further. But it could be construed as if you saw it as too obvious to dignify with further commentary.

The constructive thing would have been to answer: "if eig/eigs runs such and such an algorithm and uses operations x, y, and z, I would reasonably expect sparsity to have an advantage" and that's what I was looking for as a beginning to the discussion.

However, since you didn't answer that way, you didn't seem to understand why sparsity might not be helpful in all cases. Therefore, I don't know where else the discussion could go. And, since I got 2 votes for the response, it looks like others thought it necessary, too.

I would also point out that I did only ever offer my answer as a maybe. I said it's "hard to see" why eigenvalue decomp could benefit from sparsity, but I didn't say I knew it to be impossible and, again, you did nothing to help us see the unobvious...

[Oliver Woodford]

Look, I think I get where you're coming from. I guess you wanted to provide more insight, to understand why someone thinks what they do so you can correct it. I get that - you want to educate (important to actually be right when you do, though). I am prepared to say that I understand and respect why you asked for more info. My response didn't reflect that and for that I apologize. You were asking me for more info so you could tell me exactly where I was wrong, and that annoyed me, regrettably.

But do you not get my perspective here? What I got from you was a loaded question with an implication of fallacy and an answer that was off the point and condescending. By contrast, my example answer requires no more info, doesn't point out where the questioner's thinking is flawed, it just answers the question in a nice way, proving the "unobvious" simply wasn't necessary. Tell me you get it, Matt. If so, we understand each other. If not, oh well.

[Matt J]

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?

[Oliver Woodford]

You wouldn't use QR with a sparse matrix, you'd use Lanczos. I don't actually know what eig uses, and I only assume eigs uses Lanczos from what I read online, and looking at the code, after your answer. My reason for originally believing that sparse would be more efficient is that solving linear problems involving sparse matrices generally is, so I was in no position to either provide the answer you were looking for, or indeed verify your answer. Hence suggesting the more straightforward answer. I don't know if one can do better, but it doesn't seem likely in MATLAB.

Answer 3

[Matt J]

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.