MATLAB Answers

Diagonalising a Skew-Symmetric Matrix

33 views (last 30 days)
Jonah Berean-Dutcher
Jonah Berean-Dutcher on 22 Sep 2020 at 21:42
Answered: Christine Tobler on 28 Sep 2020 at 16:00
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 Comments

Sign in to comment.

Accepted Answer

Christine Tobler
Christine Tobler on 28 Sep 2020 at 16:00
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

  0 Comments

Sign in to comment.

More Answers (3)

sushanth govinahallisathyanarayana
sushanth govinahallisathyanarayana on 22 Sep 2020 at 21:54
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
Jonah Berean-Dutcher on 22 Sep 2020 at 22:14
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];

Sign in to comment.


Paul
Paul on 22 Sep 2020 at 22:34
Edited: Paul on 23 Sep 2020 at 1:42
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
Jonah Berean-Dutcher on 22 Sep 2020 at 23:36
Thanks, are you able to explain why eig(B) doesn't produce eigenvectors that are orthonormal for this skew-symmetric case? And why jordan(B) does?
Paul
Paul on 23 Sep 2020 at 1:57
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

Sign in to comment.


John D'Errico
John D'Errico on 23 Sep 2020 at 0:24
Edited: John D'Errico on 23 Sep 2020 at 0:30
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 Comments

Sign in to comment.

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!