# Unitary matrix with non-orthogonal eigenvectors?

6 views (last 30 days)
Joshua Feis on 2 Apr 2022
Commented: Joshua Feis on 11 Apr 2022
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?

Christine Tobler on 7 Apr 2022
Edited: Christine Tobler on 11 Apr 2022
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.
Joshua Feis on 11 Apr 2022
Thank you - this does exactly what I want it to do!

David Goodmanson on 2 Apr 2022
Edited: David Goodmanson on 4 Apr 2022
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.

R2020b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!