Should I use EIG or SVD to compute the eigenvalues of a symmetric matrix?

I saw ths question posed as an answer to another question. So I'll pose the question myself, then post my own thoughts as my answer. But feel free to add your own ideas.

The question there was given a matrix like

A = [6 2 1;2 5 2;1 2 3]

A = 3×3

6 2 1 2 5 2 1 2 3

we see that both eig and svd can be used to compute the eigenvalues and eigenvectors. Thus:

[W,D] = eig(A)

W = 3×3

0.0637 -0.7224 -0.6885 -0.5513 0.5496 -0.6277 0.8319 0.4196 -0.3633

D = 3×3

1.7511 0 0 0 3.8978 0 0 0 8.3511

[U,S,V] = svd(A)

U = 3×3

-0.6885 0.7224 0.0637 -0.6277 -0.5496 -0.5513 -0.3633 -0.4196 0.8319

S = 3×3

8.3511 0 0 0 3.8978 0 0 0 1.7511

V = 3×3

-0.6885 0.7224 0.0637 -0.6277 -0.5496 -0.5513 -0.3633 -0.4196 0.8319

And the two seem to return subtly different results. Is one better? Is one more correct? Is one faster, more efficient?

Answer 1

[John D'Errico]

Again, these are just my thoughts as I write them down. I'll possibly fail to get everything answered in here. Such is life.

What happens when we use EIG versus SVD? In the example given, it points out a real issue to remember.

EIG and SVD do not necessarily return their results in the same order. While SVD I recall always returns the singular values sorted in decreasing order, that is not true of eig.

A = [6 2 1;2 5 2;1 2 3]

A = 3×3

6 2 1 2 5 2 1 2 3

[W,D] = eig(A)

W = 3×3

0.0637 -0.7224 -0.6885 -0.5513 0.5496 -0.6277 0.8319 0.4196 -0.3633

D = 3×3

1.7511 0 0 0 3.8978 0 0 0 8.3511

[U,S,V] = svd(A)

U = 3×3

-0.6885 0.7224 0.0637 -0.6277 -0.5496 -0.5513 -0.3633 -0.4196 0.8319

S = 3×3

8.3511 0 0 0 3.8978 0 0 0 1.7511

V = 3×3

-0.6885 0.7224 0.0637 -0.6277 -0.5496 -0.5513 -0.3633 -0.4196 0.8319

The first thing to recognize is that EIG here has returned eigenvalues sorted in INCREASING order, while SVD used decreasing order. In fact, EIG explicitly does not insure a sorted sequence. Anyway, we might choose to shuffle the results from EIG to be consistent here with SVD.

ind = [3 2 1];
D = D(ind,ind)

D = 3×3

8.3511 0 0 0 3.8978 0 0 0 1.7511

W = W(:,ind)

W = 3×3

-0.6885 -0.7224 0.0637 -0.6277 0.5496 -0.5513 -0.3633 0.4196 0.8319

The next point of importance, is that even if we re-sorted the eigenvalues and eigenvectors, we can still see a difference. That is because eigenvectors are NOT unique. We can arbitrarily multiply an eigenvector by any scalar constant, and it will still work as an eigenvector. But since the eigenvectors and singular vectors are normalized by these codes to have norm of 1, that allows the codes to arbitrarily sometimes return some of those vectors with a flipped set of signs. This is the case if we now compare

[W(:,2),U(:,2)]

ans = 3×2

-0.7224 0.7224 0.5496 -0.5496 0.4196 -0.4196

Same vectors, but flipped signs. Neither one is preferred, but you need to recognize this as a non-significant "difference", and not worry about it.

Next, for a symmetric positive semi-definite matrix A, we will have (to within floating point trash) for EIG

W' * D * W == A

and for the SVD we will have

U * S * V' == A

So the columns of W are eigenvectors, as are the columns of U, but also the columns of V. Again, if A is SPD, then U==V should be the case. Here SVD comes close, but is not perfect.

U - V

ans = 3×3

1.0e-14 * 0.0222 -0.0444 0.1929 0 0.1221 -0.1443 0.0056 -0.2331 -0.1332

There are differences that are close to machine epsilon. But SVD is not explicitly designed to solve the EIGENVALUE problem. The SVD is computing a singular value decomposition, that in the case of a SPD matrix, just happens to be mathematically the same as what EIG will generate. But remember that SVD returned subtly different results for U and V here. If we ignore that, will SVD be as accurate? I'll use a larger matrix here to stress the linear algebra just a bit more.

A = cov(randn(2000,1000));
[W,D] = eig(A);
[U,S,V] = svd(A);

norm(W*D*W' - A)

ans = 1.0306e-14

norm(U*S*U' - A)

ans = 1.3701e-14

So even for this much larger problem, we are still doing impressively well just using U as the eigenvectors of A. SVD was slightly worse, but not by a lot.

Is one faster than the other? Well, yes. This is one main reason why using eig is surely preferred.

timeit(@() eig(A))

ans = 0.0899

timeit(@() svd(A))

ans = 0.1928

As you can see, eig was roughly twice as fast. For much smaller matrices, really, you don't care about time.

Finally, there are subtle issues when your matrix happens to be numerically not SPD. Then eig can generate negative eigenvalues, or in some cases, it might generate complex eigenvalues, depending on how bad things are. SVD will not do that, but now the presumption that U==V will begin to fail.

[Paul]

Why is only SPD being considered here? Isn't it true that theoretically U = W for any symmetric matrix?

Note that the original thread that motivated this Question was discussing eigenvectors, not eigenvalues, though I realize now that this Question is about eigenvalues, though this answer is also discussing eigenvectors.

[Christine Tobler]

It's a tricky question - there's an advantage in having a symmetric form U*D*U' from calling EIG, but there's also an advantage in SVD guaranteeing all nonnegative singular values, while EIG will return negative eigenvalues sometimes.

I think this question often comes up for the case where A is ill-conditioned SPD, where EIG may return some negative eigenvalues on the round-off level. Using the SVD there seems helpful, but the assumption of U and V being largely the same isn't true there anymore:

rng default;
A = randn(7, 10); A = A'*A;

[U, D] = eig(A);
norm(A - U*D*U')

ans = 2.9031e-14

norm(U'*U - eye(10))

ans = 9.1545e-16

[U, S, V] = svd(A);
norm(A - U*S*V')

ans = 2.8517e-14

norm(U'*V - eye(10))

ans = 2.0000

The large error at the end isn't down to scaling of the columns of U and V, instead those columns are in a different order in U and in V. Basically the last three columns of U and V each are a basis of the null-space of A, but they don't have to be the same basis.

No matter if EIG or SVD is used, one usually has to add some separate code to deal with round-off level eigenvalues in some way (setting them all to exact zero for example).

Answer 2

[Bruno Luong]

I always use EIG on SDP matrix, for the "same" reason that I use

A\b

preferable to

(A'*A) \ (A'*b) 

to solve least square solution of A*x = b for tall A.

I simply have intuition that SVD works on A'*A or A*A' behind the scene.