32 views (last 30 days)

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

sushanth govinahallisathyanarayana
on 22 Sep 2020

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
on 22 Sep 2020

Edited: Paul
on 23 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

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

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

Start Hunting!
## 0 Comments

Sign in to comment.