Eig function for symmetric matrices

Two of the properties of symmetric matrices are that their eigenvalues are always real, and that they are always orthogonally diagonalizable. Given any two distinct eigenvalues, the corresponding eigenvectors are orthonormal. Yet if some eigenvalue L with multiplicity greater than 1 appear, it is necessary to use, for example, Gram Schmidt for ensuring orthogonality between the vectors representing the span of Ker(A - Id L). I was so glad to see that computing the MatLab function eig over symmetric matrices having multiple eigenvalues would output an orthogonal eigenvector matrix, meaning that MatLab doesn't only normalize the vectors, but also it make them orthogonal! I tried it on a few matrices, I just need to know if I was "lucky" with them, and if usually eig doesn't do that, or, elsewise, what procedure is implemented in eig for dealing with this situation. Thank you for your help.

Answer 1

[Matt J]

I doubt you were just lucky. I've never seen EIG fail to give orthogonalized eigenvectors for symmetric matrices if it can recognize the matrix as symmetric. If there's a chance you'll have numerical error in the creation of the matrix (such that it is not perfectly symmetric), it would be safest to use SVD instead.

[Walter Roberson]

Or force them to be symmetric by using

A = (A + A.')/2;

Answer 2

[Steven Lord]

According to the documentation page for eig, specifically the section describing the output argument V:

 [V,D] = eig(A) returns matrix V, whose columns are the right eigenvectors of A such 
 that A*V = V*D. The eigenvectors in V are normalized so that the 2-norm
 of each is 1.

 If A is real symmetric, then the right eigenvectors, V, are orthonormal.

[Edited to fix formatting of quoted text.]

[Lorenzo]

Hi Steven,

I noticed that if A is complex and hermitian with repeated eigenvalues, then the right eigenvectors, V, fail to be orthonormal (V is not unitary). Do you know a way to fix this?

[Bruno Luong]

Lorenzo, make sure you input matrix is exactly Hermitian so that MATLAB engine can correctly detect it form. This can accomplish by "symmetrizing" it.

>> [Q,~]=qr(randn(3)+1i*rand(3));
>> A=Q'*diag([1 1 100])*Q

A =

    15.3593 + 0.0000i  -7.9491 - 8.0340i  23.9826 +22.6380i
    -7.9491 + 8.0340i   9.8955 - 0.0000i -25.9424 + 0.8862i
    23.9826 -22.6380i -25.9424 - 0.8862i  76.7452 + 0.0000i

>> A'-A

ans =

     1.0e-14 *

     0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
     0.0000 + 0.0000i   0.0000 + 0.1776i  -0.3553 - 0.3553i
     0.0000 + 0.0000i   0.3553 - 0.3553i   0.0000 + 0.0000i

>> [V D]=eig(A)

V =

     0.9246 + 0.0000i   0.2770 + 0.2614i  -0.4084 + 0.0100i
     0.0868 - 0.0878i  -0.2996 + 0.0102i   0.8103 + 0.0000i
    -0.2620 + 0.2473i   0.8747 + 0.0000i   0.4039 - 0.1157i

D =

     1.0e+02 *

     0.0100 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
     0.0000 + 0.0000i   1.0000 + 0.0000i   0.0000 + 0.0000i
     0.0000 + 0.0000i   0.0000 + 0.0000i   0.0100 - 0.0000i

>> V'*V % not orthonormal

ans =

     1.0000 + 0.0000i   0.0000 - 0.0000i  -0.4417 + 0.0108i
     0.0000 + 0.0000i   1.0000 + 0.0000i  -0.0000 + 0.0000i
    -0.4417 - 0.0108i  -0.0000 - 0.0000i   1.0000 + 0.0000i

>> A=0.5*(A'+A) % symmetrize it

A =

    15.3593 + 0.0000i  -7.9491 - 8.0340i  23.9826 +22.6380i
    -7.9491 + 8.0340i   9.8955 + 0.0000i -25.9424 + 0.8862i
    23.9826 -22.6380i -25.9424 - 0.8862i  76.7452 + 0.0000i

>> A'-A

ans =

       0     0     0
       0     0     0
       0     0     0

>> [V D]=eig(A)

V =

     0.3099 + 0.2813i  -0.5848 - 0.5812i   0.2770 + 0.2614i
    -0.7875 + 0.0165i  -0.5382 - 0.0087i  -0.2996 + 0.0102i
    -0.4521 + 0.0000i   0.1746 + 0.0000i   0.8747 + 0.0000i

D =

      1.0000         0         0
           0    1.0000         0
           0         0  100.0000

>> V'*V % now V is orthonormal

ans =

     1.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i
    -0.0000 + 0.0000i   1.0000 + 0.0000i  -0.0000 + 0.0000i
     0.0000 + 0.0000i  -0.0000 + 0.0000i   1.0000 + 0.0000i

>>

[Bruno Luong]

The second way is calling SVD which return orthonormal basis in any case, and happens to be eigen vectors as well for Hermitian matrix, even slightly off by numerical error:

[Q,~]=qr(randn(3)+1i*rand(3));
A=Q'*diag([1 1 100])*Q

A =

    33.3677 + 0.0000i -30.6395 - 1.8052i  19.3512 +28.9867i
   -30.6395 + 1.8052i  30.1043 - 0.0000i -19.9347 -26.3597i
    19.3512 -28.9867i -19.9347 +26.3597i  38.5280 + 0.0000i

[V D]=eig(A)

V =

     0.8204 + 0.0000i   0.3175 + 0.4756i   0.6104 - 0.1457i
     0.3772 - 0.0222i  -0.3270 - 0.4325i   0.7529 + 0.0000i
    -0.2383 + 0.3569i   0.6157 + 0.0000i   0.1977 + 0.0178i

D =

     1.0e+02 *

     0.0100 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
     0.0000 + 0.0000i   1.0000 + 0.0000i   0.0000 + 0.0000i
     0.0000 + 0.0000i   0.0000 + 0.0000i   0.0100 + 0.0000i

V'*V % V is NOT orthonormal

ans =

     1.0000 + 0.0000i  -0.0000 - 0.0000i   0.7440 - 0.1776i
    -0.0000 + 0.0000i   1.0000 + 0.0000i  -0.0000 + 0.0000i
     0.7440 + 0.1776i  -0.0000 - 0.0000i   1.0000 + 0.0000i

[U D V]=svd(A)

U =

    -0.5718 - 0.0000i  -0.0000 + 0.0000i  -0.8204 + 0.0000i
     0.5413 - 0.0319i   0.4742 - 0.5817i  -0.3772 + 0.0222i
    -0.3418 + 0.5121i  -0.1567 - 0.6421i   0.2383 - 0.3569i

D =

    100.0000         0         0
           0    1.0000         0
           0         0    1.0000

V =

    -0.5718 + 0.0000i   0.0000 + 0.0000i  -0.8204 + 0.0000i
     0.5413 - 0.0319i   0.4742 - 0.5817i  -0.3772 + 0.0222i
    -0.3418 + 0.5121i  -0.1567 - 0.6421i   0.2383 - 0.3569i

V'*A*V

ans =

   100.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
    -0.0000 - 0.0000i   1.0000 - 0.0000i  -0.0000 + 0.0000i
     0.0000 - 0.0000i  -0.0000 - 0.0000i   1.0000 - 0.0000i

V'*V % V is orthonormal

ans =

     1.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i
    -0.0000 - 0.0000i   1.0000 + 0.0000i   0.0000 + 0.0000i
     0.0000 - 0.0000i   0.0000 - 0.0000i   1.0000 + 0.0000i

[Matt J]

The second way is calling SVD which return orthonormal basis in any case, and happens to be eigen vectors as well for Hermitian matrix

That is true as long as U=V' and a Hermitian matrix always has one SVD satisfying this. However, here the SVD is non-unique. Can we be sure Matlab will always return a Hermitian symmetric U,V pair even in non-unique cases?

[Bruno Luong]

If H is input hermitian and G:=sqrtm(H) (G is hermitian as well)

If [U,D,V] is (any) svd(H), then [U,sqrt(D),V] is svd of G. Therefore by simple calculation V'*D*V is svd of G'*G = G*G = H. QED.

BTW you might ask the same question when you propose using ORTH() to fix V.

[Matt J]

Sorry, how is [U,sqrt(D),V] the svd of G for all U,V arising from svd(H)? I can see that it is true when U=V', but not for arbitrary U,V.

[Bruno Luong]

Assuming first H is full-rank, and let define R := V'*U.

From H = H' we have

U*D*V' = V*D*U'

or in other term R*D*R = D;

That means

R(i,:)*D(i,i)*R(:,j) = D(i,i)*delta(i,j)

with delta(.) is kronecker's symbol.

Since D(i,i) ~= 0, we conclude

R(i,:)*R(:,j) = delta(i,j)

Meaning: R^2 = R*R = I;

Thus R = I (let's recall R is orthonormal); meaning

U=V;

If H is not full-rank, U=V is generally not true but still this is true V*D*V' = H. I won't show that.

[Christine Tobler]

If A is symmetric (or complex hermitian), both U and V will be orthogonal matrices. However, they need not be identical if the matrix A is not positive definite.

Example:

>> A = diag([-3, 3]);
>> [U, D] = eig(A)
U =
     1     0
     0     1
D =
    -3     0
     0     3
>> [U, S, V] = svd(A)
U =
     1     0
     0     1
S =
     3     0
     0     3
V =
    -1     0
     0     1

This is because S must have nonnegative values on its diagonal, while D can have any real numbers.

This comes down to the step saying that R*R = I implies R = I; in fact, R can be any diagonal matrix with 1 and -1 on the diagonal.

[Matt J]

This comes down to the step saying that R*R = I implies R = I; in fact, R can be any diagonal matrix with 1 and -1 on the diagonal.

It can be other things as well, e.g., R=[1,1;1,-1]/sqrt(2)

[Bruno Luong]

Yeap Christine and you catch the mistake.

Actually what the proof need to work is one assumption

R (=V'*U) be in SO(n) and not only on U(n).

So in case H definite positive MATLAB SVD always return R in SO(n).

[Lorenzo]

Thank you Bruno and to all in this discussion. In the meantime I'm using schur to diagonalize an hermitian (or hermitian up to numerical errors) matrix as it seems to always return a unitary matrix. Also it seems not slower (actually a tiny bit faster) than eig. Moreover it does not have the problems of svd and I can directly get the eigenvalues (which I need).

This begs the question: shouldn't schur be the method of choice to diagonalize hermitian matrices?

[Bruno Luong]

Great find of using SCHUR

[Christine Tobler]

Calling SCHUR for a hermitian matrix should be slower than calling EIG. I still think the problem is that your matrix is not exactly hermitian.

Here's an example:

>> n = 1000;
>> A = randn(n) + 1i*randn(n);
>> A = A*diag(rand(n, 1))*A'; % nearly hermitian
>> norm(A - A')
ans =
   7.3104e-13
>> tic; [U, D] = eig(A); toc
Elapsed time is 4.432694 seconds.
>> tic; [U, T] = schur(A); toc
Elapsed time is 3.521355 seconds.
>> tic; [U, D] = eig((A+A')/2); toc
Elapsed time is 0.637474 seconds.

Because A is not hermitian, EIG uses the non-hermitian algorithm, which is based on calling SCHUR, and then recovering the eigenvectors from the Schur vectors. Since A is nearly hermitian, its eigenvectors are nearly identical to the Schur vectors.

However, for a hermitian matrix, the Schur decomposition's matrix T is diagonal. There is no need to compute the upper triangle of T, as is done in SCHUR, and the hermitian algorithm is faster because it doesn't do this.