Diagonalising a Skew-Symmetric Matrix

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?

Answer 1

[Christine Tobler]

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

[Christine Tobler]

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.

Answer 2

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

[Jonah Berean-Dutcher]

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];

Answer 3

[Paul]

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

[Jonah Berean-Dutcher]

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]

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

Answer 4

[John D'Errico]

Matrices that are defective will not have a complete set of eigenvectors.

https://en.wikipedia.org/wiki/Defective_matrix

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