generalized eigen-problom, complex non-symmetric matrix, how to find orthogonal eigenvectors

I want to solve eigen problem

with normalization condition

My code is

[U,D] = eig(B\A)
[W,D2] = eig(B.'\A.')

I use some code to rank the eigenvalue and shift the eigenvector in U and W accrodingly.

say now we have U_new, W_new

is not diagonal, in fact it's almost diagonal with some non-diagonal elements when two eigenvector share same eigenvalue.

What should i do to get

Answer 1

[Christine Tobler]

There's an easier way to compute both U and W:

[U, D, W] = eig(A, B);

This will make sure that the eigenvectors U(:, i) and V(:, i) belong together. It's also going to be more accurate because B does not need to be inverted.

I'm not sure that this will satisfy your requirement that W'*A*U is diagonal, though. I tried to construct an example matrix, but couldn't find one where it wasn't diagonal, but I couldn't find a clear statement that this will always be the case, either.

[Qingqing Zhao]

Hey Christine, thanks for your notes and help. I did try [U, D, W] = eig(A, B) before. The problem of this is when i do the senitary checking:

>> imagesc(real((A*U-B*U*D)))

this is good, i get mostly zeros ~10e-9

>> imagesc(real((W.'*A-D*W.'*B)))

ohhh i get some big values ~100

Do you know where might go wrong?

My matrix A and B are of size 2000*2000 and can go up to 20000*20000, and A is complex non-symmetry. The main issue is that there are lots of eigenvectors with same eigenvalue, over those states, it seems the algorithm didn't pick the eigenvectors that satisfy the desired orthogonality condition, i.e. W'*A*U is diagonal.

[David Goodmanson]

Hi QZ,

If there are complex values running around, should your agreement check be for

W'*A-D*W'*B, not W.'*A-D*W.'*B   ?

[Christine Tobler]

I agree with David, the issue is that the 3-output syntax of EIG is using the conjugate transpose, not the non-conjugate transposed. I had misinterpreted your A^T as being the conjugate transpose, not just the transpose.

I'm not aware of any direct algorithms for computing both U and W when this is defined using W.' instead of W'.

[David Goodmanson]

Hi Christine & QZ,

It may not help for practical computation, but if you want

A*U = B*U*D        W'*A = D*W'*B    then   
W' = (B*U)^-1     works and ensures that W'*A*U = D

and if you want

A*U = B*U*D        Wt*A = D*Wt*B

(where Wt is non-optimal notation for the transpose of W but maybe better than .' for the moment) then

Wt = (B*U)^-1     works and ensures that Wt*A*U = D

<comment moved down>

[Qingqing Zhao]

I see. Thanks, David and Christine!

It's really helpful!

[David Goodmanson]

Hi Christine & QZ,

As for computing Wt instead of W' in eig, let

[U1 D1 W1] = eig(conj(A), conj(B))

and take

U = conj(U1)
D = conj(D1)
W = W1

Then if you do the algebra I believe you get

A*U = B*U*D
Wt*A = D*Wt*B