Diagonalising a Skew-Symmetric Matrix
조회 수: 15 (최근 30일)
이전 댓글 표시
I have a skew-symmetric matrix B, and when I run:
[V,D] = eig(B)
the V matrix returned is not unitary, as I desire it to be. I think this is becuase I should not be using 'eig' for this purpose. What is the correct function, or algorithm to be using for diagonalizing a skew-symmetric matrix, in a complex vector space?
댓글 수: 0
채택된 답변
Christine Tobler
2020년 9월 28일
When EIG is called with an exactly symmetric/hermitian matrix, MATLAB falls back to a specialized algorithm that guarantees that U is orthogonal/unitary and that D is real. There is no such special algorithm choice for skew-symmetric matrices, so there is no guarantee here, even though if the problem is nicely conditioned, the result will be close to that:
>> rng default; A = randn(10); A = A - A';
>> [U, D] = eig(A);
>> max(abs(real(diag(D))))
ans =
2.1034e-16
>> norm(U'*U - eye(10))
ans =
4.7239e-15
However, if matrix B is (exactly) skew-symmetric, it implies that matrix A = 1i*B is hermitian, and passing this matrix to EIG will result in unitary eigenvectors and all-real eigenvalues, which you can then transform back:
[U, D] = eig(1i*A);
D = -1i*D;
>> max(abs(real(diag(D))))
ans =
0
>> norm(U'*U - eye(10))
ans =
1.5001e-15
댓글 수: 1
Christine Tobler
2021년 4월 15일
As of R2021a, eig will now detect if a matrix is skew-symmetric or skew-hermitian, and will return purely imaginary eigenvalues and orthogonal eigenvectors in that case. See the Release Notes. This does require the input to be exactly skew-hermitian, not just up to round-off (you can use ishermitian(A, 'skew') to check this).
Thank you for pointing out this case as a useful special case treatment to add to EIG, it's helped us improve the functionality.
추가 답변 (3개)
sushanth govinahallisathyanarayana
2020년 9월 22일
I tested this out in Matlab R2018 A
m=[0,-1;1,0];
[V,D]=eig(m)
V'*V seems to be identity for me, m is skew symmetric here. I may be missing something about your original matrix, or you could check to see if it really is skew symmetric.
Paul
2020년 9월 22일
편집: Paul
2020년 9월 23일
For your matrix B, you can diagonalize it and get the associated trasnsformation matrix as follows:
[T,J]=jordan(B);
any(any(J-diag(diag(J)))) % prove J is diagonal
ans =
logical
0
You can inspect J and see that the diagonal elements are equal to eig(B).
max(abs([cplxpair(diag(J))-cplxpair(eig(B))]))
ans =
4.4409e-16
댓글 수: 2
Paul
2020년 9월 23일
편집: Paul
2020년 9월 23일
It certainly appears that a) the columns of T are linearly independent and b) T contains eigenvectors of B that, as we've seen, diagonalizes B:
>> rank(T)
ans =
12
>> max(abs(B*T-T*J),[],'all')
ans =
0
Furthermore, now that I think about it some more, if B is defective then it should have at least one entry of 1 on the superdiagonal of J, But it doesn't. Which again indicates that B is not defective. Let's try the symbolic approach:
[Vs,Ds]=eig(sym(B));
Vs\B*Vs
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, 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, 0, 0, -2i, 0, 0, 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 2i, 0, 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, -1i, 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, -1i, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, -1i, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1i, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1i, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1i]
>> [Vs-T]
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, 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, 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, 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, 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, 0, 0, 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
So it certainly seems like B should be diagonalizable, and it is either via the transformation returned from Jordan, or via the eigenvectors returned from the Symbolic toolbox.
Try a different option for eig, which seems to work better for your case.
[V1,D1]=eig(B,'nobalance');
rank(V1)
ans =
12
max(max(abs(V1\B*V1-D1)))
ans =
1.7114e-08
John D'Errico
2020년 9월 23일
편집: John D'Errico
2020년 9월 23일
Matrices that are defective will not have a complete set of eigenvectors.
" In particular, an n × n matrix is defective if and only if it does not have n linearly independent eigenvectors."
The classic example that I know of is
>> A = triu(ones(2));
>> [V,D] = eig(A);
V =
1 -1
0 2.22044604925031e-16
D =
1 0
0 1
As you can see, there is a duplicate eigenvalue. But you can also see the two eigenvectors (columns of V) are not orthogonal.
In this case, the output from eig still satisfies the relation that A*V == V*D, at least within floating point trash.
>> norm(A*V - V*D)
ans =
2.22044604925031e-16
How about the case of your matrix B?
>> [V,D] = eig(B);
>> norm(B*V - V*D)
ans =
8.74077768365071e-16
So again, at best we can see this norm is zero. But V is not a complete set of eigenvectors. Why not? Because B is defective.
What can I say? If you read the help for eig, all it can promise is
"[V,D] = eig(A) produces a diagonal matrix D of eigenvalues and
a full matrix V whose columns are the corresponding eigenvectors
so that A*V = V*D."
When the matrix is defective, it can do no better than that, since the set of eigenvectors you are asking it to produce apparently don't exist.
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Matrix Indexing에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!