MATLAB Answers

Diagonalising a Skew-Symmetric Matrix

조회 수: 19(최근 30일)
Jonah Berean-Dutcher
Jonah Berean-Dutcher 2020년 9월 22일
댓글: Christine Tobler 2021년 4월 15일
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?

채택된 답변

Christine Tobler
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
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
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.
  댓글 수: 1
Jonah Berean-Dutcher
Jonah Berean-Dutcher 2020년 9월 22일
It fails for my original matrix:
B = [ ...
0 0 -1 0 0 0 0 0 0 0 0 0; ...
0 0 0 -1 0 0 0 0 0 0 0 0; ...
1 0 0 0 0 0 0 0 0 0 0 0; ...
0 1 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 -1 0 -1 0 0; ...
0 0 0 0 0 0 1 0 0 0 -1 0; ...
0 0 0 0 0 0 0 0 0 0 0 -1; ...
0 0 0 0 0 0 1 0 0 0 -1 0; ...
0 0 0 0 0 0 0 1 0 1 0 0; ...
0 0 0 0 0 0 0 0 1 0 0 0];

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


Paul
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
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
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.

제품


릴리스

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by