Why does "polyeig" returns eigenvectors that have smaller size than number of eigenvalues?

I tried to solve an eigenvalue problem with " polyeig ". It correctly evaluated eigenvalues, however eigenvectors are not same as correct answer, which are obtained using " eig ". Do you know why " polyeig " returns smaller size eigenvectors?

This matter is clearly noted in the explanation of " polyeig ": [X,e] = polyeig(A0,A1,...,Ap) also returns matrix X, of size n-by-n*p, whose columns are the eigenvectors.

[Steven Lord]

Please show us a small sample of code that illustrates the difference.

[Masoud]

In this code, eigenvectors of matrix A is calculated using both "eig" and "polyeig" and size of answer matrices are not same (I mean number of rows of X1 and X2).

Answer 1

[David Goodmanson]

Hi Masoud,

Your matrix A is 402x402. Denote Gamma by G and Omega by W, both 201x201. Let each eigenvector of A be the concatenation [u; v] where u and v are each 201x1. The system of equations is

   Wu + Gv = bu     (1)
 -G'u - Wv = bv     (2)

with eigenvalue b.

Working out polyeig for u gives the formulas you got. Both polyeig(A0,A21,A2) and eig(A) have 402 eigenvalues, and each polyeig eigenvector u is 201x1 as you noted. To obtain the second half of the eigenvector of A, you can invert (1) and use

v = G\(bu-Wu)

[Masoud]

Dear David,

Thank you for your answer. I tried your solution and here is the results that I got.

The blue points are results from "eig" and the red ones are results from the way you suggested. A I doing something wrong? (I will double check and let you if I get the right answer)

I really need to solve this problem.

[David Goodmanson]

Hi Masoud,

I'm taking the point of view that it doesn't really matter what Gamma and Omega are as long as they are complex square matrices of the same size. The code below is basically a copy of yours once the matrices have been created. All the eigenvalues come out the same between the two methods. So it appears to work. But your case appears to come out differently ?!

If in your case Gamma has a large condition number, there could of course be some problems.

n = 201;
Gamma = 2*rand(n,n)-1 +i*(2*rand(n,n)-1);
Omega = 2*rand(n,n)-1 +i*(2*rand(n,n)-1);

A=[Omega Gamma;-Gamma' -Omega];

A0=Gamma'-((Omega/Gamma)*Omega);
A1=(Omega/Gamma)-(Gamma\Omega);
A2=inv(Gamma);

[X1 e1] = polyeig(A0,A1,A2);     % e1 is a vector
[X2 e2] = eig(A);                % e2 is a diagonal matrix 
e2vec = diag(e2);

figure(1)
plot(e1)
hold on
plot(e2vec,'o')
hold off

[Masoud]

Why do you plotted eigenvalues in your code? My problem is eigenvectors, and even in your example, your proposed formula is not working.