eig() gives wrong results for complex eigenvalues in comparison to eigs()

Hi,

I'm working on a eigenvalue problem of the following two matrices A and B

A = [M zeros(N); zeros(N) -K];  
B = [zeros(N) M; M D];

M, D and K are submatrices of size NxN. Let's say:

N = 1;
K = rand(N); D = rand(N); M = rand(N);

The eigenvalue problem can be formulated as:

A*x = lambda*B*x

or

y.' *A = lambda * y.'*B

with the right and left eigenvectors or modal matrices x, y and X, Y.

Obviously, there are 6 ways to calculate the eigenvalue problem:

[X1, L1, Y1] = eig(A, B);  
[X2, L21] = eigs(A, B);  
[Y2, L22] = eigs(A.', B.');
%
[X3, L3, Y3] = eig(B \ A);  
   Y3 = (Y3.' / B).';  
[X4, L41] = eigs(B \ A);  
[Y4, L42] = eigs((B \ A).');  
   Y4 = (Y4.' / B).';
%
[X5, L5, Y5] = eig(A / B);  
   X5 = B \ X5;  
[X6, L61] = eigs(A / B);  
[Y6, L62] = eigs((A / B).');  
   X6 = B \ X6;

The corrections terms of the modal matrices Y3, Y4, X5 and X6 are required in case of using left and right multiplication of the matrix B.

We will find the eigenvalues lambda1-lambda6 as diagonal entries of the matrices L1-L6 to be equal.

But, the following equation should give us a zero solution vector

          (A - lambda_i * B) * x_i = 0
  y_j.' * (A - lambda_j * B)       = 0

Unfortunately, this is not the case for (X1, Y1), (X3, Y3) and (X5, Y5) for complex eigenvalues. However, for completely real eigenvalues, they give the correct solution.

for qq = 1:length(l1)
%
    right1(:, qq) = (A - l1(qq) * B) * X1(:, qq);
    left1(qq, :) = Y1(:, qq).' * (A - l1(qq) * B);
%
    right2(:, qq) = (A - l2(qq) * B) * X2(:, qq);
    left2(qq, :) = Y2(:, qq).' * (A - l2(qq) * B);
%
    right3(:, qq) = (A - l3(qq) * B) * X3(:, qq);
    left3(qq, :) = Y3(:, qq).' * (A - l3(qq) * B);
%
    right4(:, qq) = (A - l4(qq) * B) * X4(:, qq);
    left4(qq, :) = Y4(:, qq).' * (A - l4(qq) * B);
%
    right5(:, qq) = (A - l5(qq) * B) * X5(:, qq);
    left5(qq, :) = Y5(:, qq).' * (A - l5(qq) * B);
%
    right6(:, qq) = (A - l6(qq) * B) * X6(:, qq);
    left6(qq, :) = Y6(:, qq).' * (A - l6(qq) * B);
end
[left1 right1;
 left2 right2;
 left3 right3;
 left4 right4;
 left5 right5;
 left6 right6]

Can you please tell me, why "eig()" gives "wrong" eigenvectors in this case?

Thank you very much for your help, regards

Lennart

Answer 1

[Nalini Vishnoi]

Hi Lennart,

It seems there are a few things that are not correct with the last part of the code:

1. Instead of li(qq), it should be li(qq,qq).

2. You should not use (.'), instead (') should be used. (.') computes the nonconjugate transpose of a matrix whereas (.) computes complex conjugate transpose. Since the eigenvalues and eigenvectors could be complex in nature even for a real matrix, the latter is the correct form to be used.

3. The code below gives you the desired results:

for qq = 1:length(L1)
  %
  right1(:, qq) = (A - L1(qq,qq) * B) * X1(:, qq);
  left1(qq, :) = Y1(:, qq)' * (A - L1(qq,qq) * B);
  %
  right2(:, qq) = (A - L21(qq,qq) * B) * X2(:, qq);
  left2(qq, :) = Y2(:, qq)' * (A - L22(qq,qq)' * B);
  %
  right3(:, qq) = (A - L3(qq,qq) * B) * X3(:, qq);
  left3(qq, :) = Y3(:, qq)' * (A - L3(qq,qq) * B);
  %
  right4(:, qq) = (A - L41(qq,qq) * B) * X4(:, qq);
  left4(qq, :) = Y4(:, qq)' * (A - L42(qq,qq)' * B);
  %
  right5(:, qq) = (A - L5(qq,qq) * B) * X5(:, qq);
  left5(qq, :) = Y5(:, qq)' * (A - L5(qq,qq) * B);
  %
  right6(:, qq) = (A - L61(qq,qq) * B) * X6(:, qq);
  left6(qq, :) = Y6(:, qq)' * (A - L62(qq,qq)' * B);
end
[left1 right1;
left2 right2;
left3 right3;
left4 right4;
left5 right5;
left6 right6]

Please note that I have changed Li2(qq,qq) to Li2(qq,qq)' for the computation of left2, left4 and left6 in the 2nd, 4th and 6th case. Here is why: As per the documentation of eig,

[Y2, L22] = eigs(A', B');

ensures that the solution satisfies the following:

A'*Y2 = B'*Y2*L22   % which gives  
Y2'*A = L22'*Y2'*B  % (taking transpose of both sides)
Y2'*A - L22'*Y2'*B = 0  
Y2'*(A-L22'*B) = 0

(since values of L22 are taken one by one in your calculations, so it can behave as a scalar and moving its position does not change the result). However, taking its (complex conjugate) transpose is crucial since it is a complex number. Instead of computing the residuals inside a loop, it can be written as one single statement (treating L22 as a matrix):

left2 = Y2'*A - L22'*Y2'*B;

The same reasoning applies for cases 4 and 6. Please note this reasoning does not apply to cases 1, 3 and 5 and the original form of equation stays valid.

I hope the above helps.

Nalini

[lennart]

Thank you Nalini.

I didn't know that the transpose of the left eigenvector and the modal matrix also requires the complex conjugate, since I just found a simple y^T in the literature.

Just one remark: The last step in your solution

 Y2'*A - L22'*Y2'*B = 0  
 Y2'*(A-L22'*B) = 0

is not valid, since Y2' is between L22' and B in the first line and changes position in the second, which is thus not the exact solution anymore.

[Nalini Vishnoi]

Hi Lennart,

You are right, if L22 is a matrix, this is not true. However, I did mention in my previous answer that for this specific case I am treating L22 as a scalar (since you are calculating values for each eigenvalue one by one so shifting a scalar value would not change the answer). If L22 is treated as a matrix, the above won't hold true.