Eigenvectors are not orthogonal for some skew-symmetric matrices, why?

             0   -0.5000         0         0         0    0.5000
        0.5000         0   -0.5000         0         0         0
             0    0.5000         0   -0.5000         0         0
A =          0         0    0.5000         0   -0.5000         0
             0         0         0    0.5000         0   -0.5000
       -0.5000         0         0         0    0.5000         0

The above matrix is skew-symmetric. When I use [U E] = eig(A), to find the eigenvectors of the matrix. These eigenvectors must be orthogonal, i.e., U*U' matix must be Identity matrix. However, I am getting U*U' as

    0.9855   -0.0000    0.0410   -0.0000   -0.0265    0.0000
   -0.0000    0.9590    0.0000    0.0265   -0.0000    0.0145
    0.0410    0.0000    0.9735   -0.0000   -0.0145    0.0000
   -0.0000    0.0265   -0.0000    1.0145    0.0000   -0.0410
   -0.0265   -0.0000   -0.0145    0.0000    1.0410   -0.0000
    0.0000    0.0145    0.0000   -0.0410   -0.0000    1.0265

Here we can observe a substantial error. This happens for some other skew-symmetric matrices also. Why this large error is being observed and how do I get correct eigen-decomposition for all skew-symmetric matrices?

Answer 1

[Roger Stafford]

Your matrix A is "defective" , meaning that its eigenvalues are not all distinct. In fact, it has only three distinct eigenvalues. Consequently the space of eigenvectors does not fully span six-dimensional vector space. See the Wikipedia article:

http://en.wikipedia.org/wiki/Defective_matrix

What you are seeing is not an error on Matlab's part. It is a mathematical property of such matrices. You cannot achieve what you call "correct eigen-decomposition" for such matrices.

[Roger Stafford]

@Rahul: I retract my statement of yesterday that your skew-symmetric matrix is defective. Actually its eigenvectors do in fact span the entire six-dimensional space, as can be seen using the 'null' function. Moreover, according to the article at

http://en.wikipedia.org/wiki/Skew-symmetric_matrix

any real skew-symmetric matrix should always be diagonalizable by a unitary matrix, which I interpret to mean that its eigenvectors should be expressible as an orthonormal set of vectors. I don't know why Matlab doesn't produce such a set with its 'eig' function, but it also failed to do so on my own version of Matlab.

It would seem that I haven't earned the four points you kindly gave to me with your acceptance.

[Star Strider]

@Roger — That sounds like a bug report or enhancement request. TMW needs to know if MATLAB functions are producing anomalous results.

[Roger Stafford]

@ Rahul & Star: I have just verified that the eigenvectors produced (on my computer) for Rahul's matrix can easily be adjusted so as to constitute an orthonormal set and still be valid eigenvectors for A. Presumably this would also be true for Rahul's version of 'eig'. So, this leaves the question of why 'eig' is seemingly not handling these cases of repeated eigenvalues properly. (In Rahul's A there were three distinct eigenvalues, each repeated once.)

Added: I see in the documentation for 'eig' at

http://www.mathworks.com/help/matlab/ref/eig.html

that Mathworks only claims that the eigenvectors will be orthonormal for a real symmetric matrix, so the above complaint would probably have to be in the form of a request for enhancement.

[Rahul Singh]

@Roger: How did you adjust the eigenvectors to constitute an orthonormal set, and can it be generalized for this error?

How this error is going to be resolved? How to request for enhancement?

[Star Strider]

Requesting the enhancement is straightforward. Go here, click on ‘Create Service Request’, and then choose the appropriate option when the opportunity presents itself. Include the URL of this thread.

[Roger Stafford]

@Rahul: Once the repetition of eigenvalues is recognized, the adjustment can easily be done using matlab's 'orth' function. In the case of your particular A, the first and fourth eigenvalues in diag(E) are equal (on my computer), and similarly for the second and fifth, and for the third and sixth, so the adjustment (in my case) would simply be:

   U(:,[1,4]) = orth(U(:,[1,4]));
   U(:,[2,5]) = orth(U(:,[2,5]));
   U(:,[3,6]) = orth(U(:,[3,6]));

Note that the resulting U, besides being orthonormal, will still constitute a valid set of eigenvectors, even though there is an arbitrary aspect to the result of 'orth'.

The hard part is recognizing such eigenvalue equalities and in fact recognizing matrices that are diagonalizable. There really are matrices that cannot be diagonalized and which are therefore designated as "defective". For these reasons it ought to be Mathworks that carries out such a revision.

However, you can experiment on your own using 'orth' to see how it works. Remember, both the eigenvalues and the eigenvectors will be complex-valued for your skew-symmetric matrices, and in testing the adjusted U'*U you will get tiny imaginary components due to rounding errors.

By the way, in requesting a change, you should probably not refer to the current version as being in "error", since Mathworks up to this point hasn't promised orthonormal eigenvectors for matrices that are other than real and symmetric.

[Rahul Singh]

@Roger and Star: Thanks a lot.

[Lorenzo]

That matrix is not defective (1i times the matrix is hermitian and so it has a complete set of eigenvectors), it has however degenerate eigenvalues and this is the reason why U fails to be unitary. You should use schur in this case, which always return a unitary matrix

Answer 2

[Rahul Singh]

Is there any other way (other than Matlab) of computing orthogonal eigenvectors for this particular skew-symmetric matrix ? I tried NumPy package also, which gave me same results as Matlab.

[Rahul Singh]

@Roger: I tried your method using ''orth,'' but U'*U and U*U' are still entirely different.

Also, as the eigenvector pairs must be complex conjugate of each other in this case; using ''orth'' is spoiling that too.

[Roger Stafford]

The adjusted U'*U and U*U' values using 'orth' are not significantly different on my computer. They differ only from the 6 x 6 identity matrix and from each other out at the 15-th or 16-th decimal place which is what can be expected from round-off errors. Correspondingly tiny imaginary parts are also present, again due to round-off errors. Are you sure you have like eigenvalues matched? Yours may have been given in a different order from mine.

[Rahul Singh]

Are the eigenvector pairs still orthogonal?

[Roger Stafford]

Yes, all the eigenvectors come out orthogonal after that adjustment I described. The fact that U'*U gives the identity matrix implies that. You should be able to check that for yourself.

Answer 3

[Christine Tobler]

Since, as Lorenzo points out in a comment above, 1i*A is hermitian, you could apply eig to that matrix:

>> [U, D] = eig(1i*A);
>> D = D/1i;
>> norm(U'*U - eye(6))
ans =
   1.4373e-15
>> norm(A*U - U*D)
ans =
   7.8098e-16

[Paul]

"As of R2021a, eig will now detect if a matrix is skew-symmetric or skew-hermitian, and will return purely imaginary eigenvalues and orthogonal eigenvectors in that case. See the Release Notes. This does require the input to be exactly skew-hermitian, not just up to round-off (you can use ishermitian(A, 'skew') to check this)." Source

A = [        0   -0.5000         0         0         0    0.5000
        0.5000         0   -0.5000         0         0         0
             0    0.5000         0   -0.5000         0         0
             0         0    0.5000         0   -0.5000         0
             0         0         0    0.5000         0   -0.5000
       -0.5000         0         0         0    0.5000         0];

[U,D] = eig(A,'vector');
any(real(D))

ans = logical
   0

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

ans = 3.6854e-15

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

ans = 3.6670e-15