Unitary matrix with non-orthogonal eigenvectors?

I have the following issue: A matrix I have come across (attached) is pretty clearly unitary (as tested below) but its eigenvectors as calculated by eig() are not orthonormal (as seen by the matrix of eigenvectors not being unitary). Unless I am missing something mathematically, that shouldn't happen. So what went wrong here?

Here is some minimal code to calculate the eigenvectors/values and check for unitarity:

[M,ei] = eig(full(W));

norm(full(W*W')-eye(size(W))) %spits out 5.8622e-16 - pretty unitary.
norm(full(M*M')-eye(size(M))) %spits out 6.2597 - clearly not unitary.

A quick

plot(abs(diag(ei)),'o')

confirms that at least the eigenvalues of this matrix are behaving as expected - their abs() is 1.

Any ideas?

Answer 1

[Christine Tobler]

The short answer is that it's possible to compute an orthonormal basis of eigenvectors for an orthogonal matrix, but that MATLAB doesn't check for orthogonal matrices in EIG and so just provides an answer for a generic nonsymmetric matrix. For general nonsymmetric matrices, there typically isn't an orthogonal basis of eigenvectors.

However, if you know that your matrix is orthogonal, you can compute an orthogonal eigenbasis as follows:

% Construct an orthogonal matrix with a repeated eigenvalue 1 (if they are
% all different, eig will already return an orthonormal basis of
% eigenvectors)
rng default;
Q = blkdiag(eye(5), orth(randn(5)));
V = orth(randn(10));
Q = V*Q*V';
eig(Q)

ans = 
   0.0621 + 0.9981i
   0.0621 - 0.9981i
   0.9427 + 0.3335i
   0.9427 - 0.3335i
   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0000i

% Compute eigenvalues using eig:
[U, D] = eig(Q);
res = norm(Q*U - U*D) % residual shows these are correct eigenvalues and eigenvectors

res = 1.3813e-15

orthogonality = norm(U'*U - eye(10)) % but the eigenvectors aren't orthogonal

orthogonality = 0.4988

% Instead, compute the Schur decomposition:
[U, S] = schur(Q, 'complex');
S

S = 
   0.0621 + 0.9981i  -0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i
   0.0000 + 0.0000i   0.0621 - 0.9981i   0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.9427 + 0.3335i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.9427 - 0.3335i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.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   0.0000 + 0.0000i  -0.0000 + 0.0000i
   0.0000 + 0.0000i   0.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  -0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.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
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.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   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.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   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   1.0000 + 0.0000i

 % since S is a triangular matrix which we know to also be orthogonal up to round-off,
 % it must be diagonal up to round-off. Set D to be a diagonal matrix
 % containing the same diagonal as S.
D = diag(diag(S));

res = norm(Q*U - U*D) % residual shows these are correct eigenvalues and eigenvectors

res = 2.6971e-15

orthogonality = norm(U'*U - eye(10)) % and the eigenvectors are also orthogonal

orthogonality = 1.2572e-15

We're unlikely to make EIG try to detect if an input matrix is orthogonal and do this kind of operation, because we can only determine if a matrix is orthogonal up to a tolerance, and would have to make a choice of a cut-off at which we don't count a matrix as orthogonal anymore. These kinds of hidden behaviors have usually been more trouble than they are worth: imagine if you were passing a bunch of orthogonal matrices to EIG in some function, and relying on the eigenvectors returned being orthogonal. If one of them was slightly further from being orthogonal, the eigenvectors could suddenly be returned as not orthogonal at all for this matrix.

[David Goodmanson]

Hi Christine,

good one. I had not realized before that a triangular unitary matrix has to be diagonal. I decide to check the whole process for matrices with complex entries and with complex repeated eigenvalues, and naturally it worked.

[Bruno Luong]

+1, it's the first time I see Schur's decomposition does something neat.

[Joshua Feis]

Thank you - this does exactly what I want it to do!

Answer 2

[David Goodmanson]

Hi Joshua,

I believe this happens because W has at least one set of degenerate eigenvalues, and does not happen if all the eigenvalues are distinct.

mindiff = min(abs(diff(sort(eig(full(W))))))    % differences between sorted eigenvalues
   mindiff = 2.8189e-18                         % numerically consistent with zero.

and W appears to have some degenerate eigenvalues. Although It might not be all that easy to do, the nonumitarity of M appears to be fixable. To take the simplest possible case, suppose

W = [1 0;
     0 1]              % the identity
   
[M lam] = eig(W,'vec')
   M =
       1     0
       0     1
  lam =
         1             % degenerate eigenvalues
         1

Obviously W is unitary, and so is M. But the solution M for the eigenvalue problem is not the only one, Another is, e.g.

   M1   = [1  .8;
           0  .6]

which has two normalized column vectors that span the same space as M. This solution theoretically could have been chosen by eig. However, M1 is not even close to being unitary.

Any two linearly independent vectors that span the same space as M or M1 are also a good solution, so a fix is

M2 = orth(M1)
   M2 =  0.9487   -0.3162
         0.3162    0.9487

which is unitary. So a way out would be to find all sets of degenerate eigenvalues (remembering that eig does not sort eigenvalues) within some appropriate tolerance, and then use orth on whatever set of eigenvectors correspond to each set.