LAPACK fails calculating eigenvectors

Question

I wrote some code to solve the general eigenvalue problem and now I am comparing my results against LAPACK's DSPGVX function. I just worked with this example.

So I obtained the 4 auto vectors

{
 {-0.0319133, -0.265466, -0.713483,  0.64765},
 {-0.425628,  -0.520961, -0.714215,  0.193227},
 { 0.32702,    0.565845, -0.37129,  -0.659561},
 {-0.682699,  -0.056645,  0.0771025, 0.724409}
}

and auto values

{-2.22545, 1.12704, -0.454756, 0.100076}

both with my code and with Mathematica and results agree.

But in the previous link, auto vectors reported from LAPACK are completely different.

 Eigenvalues
    -0.4548  0.1001
 Selected eigenvectors
          1       2
 1   0.3080  0.4469
 2   0.5329  0.0371
 3  -0.3496 -0.0505
 4  -0.6211 -0.4743

Whom should I trust?

P.S. I also checked that my auto values/autovectors are correct since they yield A*x-lambda*B*x=0, while the values from LAPACK do not.

Answer 1

It looks like DSGPVX is solving A*lambda = B*x*lambda; Matlab gives the DSGPVX solution to your problem using "eig", though Matlab's documentation is correct. My guess is this a bug in the DSGPVX documentation.

>> a=[0.24 0.39 0.42 -0.16;0.39 -0.11 0.79 0.63;0.42 0.79 -0.25 0.48;-0.16 0.63 0.48 -0.03];
>> b=[4.16 -3.12 0.56 -0.1;-3.12 5.03 -0.83 1.09;0.56 -0.83 0.76 0.34;-0.1 1.09 0.34 1.18];
>> [v,d]=eig(a,b)

v =

   -0.0690    0.3080   -0.4469   -0.5528
   -0.5740    0.5329   -0.0371   -0.6766
   -1.5428   -0.3496    0.0505   -0.9276
    1.4004   -0.6211    0.4743    0.2510


d =

   -2.2254         0         0         0
         0   -0.4548         0         0
         0         0    0.1001         0
         0         0         0    1.1270

>> norm(a*v-b*v*d)

ans =

   1.5001e-15

Answer 2

I don't know why you think that LAPACK is giving incorrect answers, they seem fine to me. Using the four figure decimals that you quote I get residuals (r = A*x - lambda*B*x) such that

norm(r1) = 1.5921e-04, norm(r2) = 6.0842e-05.

Since norm(A) = 1.2994 and norm(B) = 7.9874, these residuals seem very satisfactory.

The eigenvectors produced by DSPGVX are normalized so that

norm(x'*B*x) = 1.

Answer 3

It looks like Lapack's results do actually correspond to the last two eigenvalues generated by your code and Mathematica, though with lower-order bits coming out quite different. The corresponding vectors are quite close, just scaled differently.

Clearly, if your/Mathematica's values check out and Lapack's do not, you should trust the one that produces the correct answers. It might be valuable to investigate what it is about your problem and and Lapack's algorithms that make it provide very imprecise answers.