Diagonalising a Skew-Symmetric Matrix
Show older comments
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?
Accepted Answer
Christine Tobler
on 28 Sep 2020
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 Comment
Christine Tobler
on 15 Apr 2021
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.
More Answers (3)
sushanth govinahallisathyanarayana
on 22 Sep 2020
0 votes
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 Comment
Jonah Berean-Dutcher
on 22 Sep 2020
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 Comments
Jonah Berean-Dutcher
on 22 Sep 2020
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
on 23 Sep 2020
Edited: John D'Errico
on 23 Sep 2020
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.
Categories
Find more on Linear Algebra in Help Center and File Exchange
on 22 Sep 2020
on 15 Apr 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
