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?
0 个评论
采纳的回答
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)
% Compute eigenvalues using eig:
[U, D] = eig(Q);
res = norm(Q*U - U*D) % residual shows these are correct eigenvalues and eigenvectors
orthogonality = norm(U'*U - eye(10)) % but the eigenvectors aren't orthogonal
% Instead, compute the Schur decomposition:
[U, S] = schur(Q, 'complex');
S
% 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
orthogonality = norm(U'*U - eye(10)) % and the eigenvectors are also orthogonal
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 个评论
更多回答(1 个)
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.
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!