Unitary matrix with non-orthogonal eigenvectors?

조회 수: 14(최근 30일)
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
Christine Tobler 2022년 4월 7일
편집: Christine Tobler 2022년 4월 11일
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.
  댓글 수: 3
Joshua Feis
Joshua Feis 2022년 4월 11일
Thank you - this does exactly what I want it to do!

댓글을 달려면 로그인하십시오.

추가 답변(1개)

David Goodmanson
David Goodmanson 2022년 4월 2일
편집: David Goodmanson 2022년 4월 4일
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.

Community Treasure Hunt

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

Start Hunting!

Translated by