Matlab Not Returning Orthonormal Matrix of Eigenvectors

Question

When I try to find the eigen-decomposition of a matrix in Matlab that has a repeated eigenvalue but is NOT defective, it is not returning an orthonormal matrix of eignevectors. For example:

k = 5;
repeats = 1;

% First generate a random matrix of eignevectors that is orthonormal
V = orth(rand(k));

% Now generate a vector of eigenvalues with the given number of repeats
D = rand(k,1);
for i = 1:repeats
    % Put one random value into another (note this sometimes will result in
    % less than the given number of repeats if we ever input the same
    % number)
    D(ceil(k*rand())) = D(ceil(k*rand()));
end

A = V'*diag(D)*V;

% Now test the eignevector matrix of A
[V_A, D_A] = eig(A);

disp(V_A*V_A' - eye(k))

I am finding that my matrix of eigenvectors V_A is not orthogonal i.e. V_A*V_A' is not equalling the identity matrix (taking into account rounding errors).

I was under the impression that if my matrix was real and symmetric then Matlab would return an orthogonal matrix of eigenvectors, so what is the issue here?

Answer 1

This seems to be a numerical precision issue.

The eigenvectors of a real symmetric matrix are orthogonal. But your input matrix A is not exactly symmetric. The differences are on the order of eps, as expected from numerical errors.

>> A-A.'
ans =
   1.0e-16 *
         0   -0.2082   -0.2776         0    0.1388
    0.2082         0         0   -0.1388         0
    0.2776         0         0   -0.2776         0
         0    0.1388    0.2776         0   -0.5551
   -0.1388         0         0    0.5551         0

If you force A to be exactly symmetric you'll get an orthogonal V_A, up to numerical errrors on the order of eps:

>> A = (A+A.')/2;
>> A-A.'
ans =
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
>> [V_A, D_A] = eig(A);
>> disp(V_A*V_A' - eye(k))
   1.0e-15 *
   -0.3331    0.2220    0.0755    0.1804         0
    0.2220   -0.2220    0.0572   -0.1665    0.1110
    0.0755    0.0572   -0.8882   -0.0590   -0.0763
    0.1804   -0.1665   -0.0590         0   -0.0555
         0    0.1110   -0.0763   -0.0555         0

Still, it's surprising that so wildly different results are obtained for V_A when A is symmetric and when A is nearly symmetric. This is my bet as to what's happening: as noted by @ArturoMagidin,

(1) Eigenvectors corresponding to distinct eigenvalues of a symmetric matrix must be orthogonal to each other. Eigenvectors corresponding to the same eigenvalue need not be orthogonal to each other.

(2) However, since every subspace has an orthonormal basis,you can find orthonormal bases for each eigenspace, so you can find an orthonormal basis of eigenvectors.

Matlab is probably taking route (2) (thus forcing V_a to be orthogonal) only if A is symmetric. For A not exactly symmetric it probably takes route (1) and gives you a basis of each subspace, but not necessarily with orthogonal vectors.