Eigenvalue problem which gives incorrect eigenvalues and eigenvectors

[troels_roennow]

Hey,

I am having a pretty strange problem. I am using an DSYGV to solve get the eigenvalues E of of the generalized eigenvalue problem H*v=E*S*v. Letting H be

H = [ -5.01054 -2.35237 -2.35237 -4.02725 -2.18694 -2.18694
-2.35237 9.67769 -0.8082 -2.18694 5.56261 -4.30757
-2.35237 -0.8082 9.67769 -2.18694 -4.30757 5.56261
-4.02725 -2.18694 -2.18694 19.3205 12.5389 12.5389
-2.18694 5.56261 -4.30757 12.5389 15.5401 7.40994
-2.18694 -4.30757 5.56261 12.5389 7.40994 15.5401];

and S as

S = [
0.9069 0.57735 0.57735 0.57735 0.367553 0.367553
0.57735 0.45345 0.367553 0.367553 0.288675 0.216506
0.57735 0.367553 0.45345 0.367553 0.216506 0.288675
0.57735 0.367553 0.367553 0.45345 0.288675 0.288675
0.367553 0.288675 0.216506 0.288675 0.226725 0.163357
0.367553 0.216506 0.288675 0.288675 0.163357 0.226725];

I get the eigenvalues:

878.5320, -5.6977, 261.0238, 122.8077, -394.6588, 131.3400.

This result is reproducable with Matlab, Octave and Intel Lapack Routines. Thus, it should be expected to be correct. Howvever, examining the problem more closely we see

chol(S) =

0.95231 0.60626 0.60626 0.60626 0.38596 0.38596
0.00000 0.29308 0.00000 0.00000 0.18658 -0.05966
0.00000 0.00000 0.29308 0.00000 -0.05966 0.18658
0.00000 0.00000 0.00000 0.29308 0.18658 0.18658
0.00000 0.00000 0.00000 0.00000 0.06766 0.02726
0.00000 0.00000 0.00000 0.00000 0.00000 0.06192

and

A = inv(chol(S)')*H*inv(chol(S)) =

-5.52491 3.00041 3.00038 -3.00050 0.22102 0.14418
3.00041 123.89277 1.81760 -1.82022 -1.41899 -59.34131
3.00038 1.81760 123.89274 -1.82020 -54.88206 22.61392
-3.00050 -1.82022 -1.82020 260.97810 3.38291 2.20684
0.22102 -1.41899 -54.88206 3.38291 499.09770 576.99797
0.14418 -59.34131 22.61392 2.20684 576.99797 -8.98945

which is just the Cholesky factorisation. Now the eigenvalues of A*u=Eu is still

878.5320, -5.6977, 261.0238, 122.8077, -394.6588, 131.3400

as they are suppose to (thus the factorisation were done correct), but we can no see that the eigenvectors are incorrect as

[V,E] = eig(A)

gives

V =
2.7700e-05 9.9940e-01 1.1527e-02 3.2697e-02 1.1019e-17 5.2984e-17
4.4573e-02 -2.3209e-02 1.1754e-02 7.0522e-01 -9.3845e-02 7.0085e-01
4.4573e-02 -2.3209e-02 1.1754e-02 7.0522e-01 9.3846e-02 -7.0086e-01
-6.7902e-03 1.0960e-02 -9.9976e-01 1.7453e-02 -5.0707e-17 -2.7167e-16
-8.3586e-01 -1.7756e-03 6.5768e-03 5.2662e-02 5.4154e-01 7.2513e-02
-5.4527e-01 -1.1583e-03 4.2904e-03 3.4354e-02 -8.3013e-01 -1.1116e-01

E =
878.53198 0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 -5.69773 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 261.02377 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 122.80774 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 -394.65885 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000 131.34004

and by checking one by one we see

(A*V(:,1))./V(:,1) )' =
878.53 878.53 878.53 878.53 878.53 878.53

(A*V(:,2))./V(:,2) )' =
-5.6977 -5.6977 -5.6977 -5.6977 -5.6977 -5.6977

(A*V(:,3))./V(:,3) )'=
261.02 261.02 261.02 261.02 261.02 261.02

(A*V(:,4))./V(:,4) )'=
122.81 122.81 122.81 122.81 122.81 122.81

(A*V(:,5))./V(:,5) )'=
-396.44 -394.66 -394.66 -287.12 -394.66 -394.66

(A*V(:,6))./V(:,6) )'=
130.28 131.34 131.34 264.02 131.34 131.34

While the first 4 eigenvalues and eigenvectors are correct, the last two are obviously not relaiable. As mentioned in the top this results is reproducable with Matlab, Octave as well as the Lapack libraries. My question is:

- Has anyone a solutions such I only would get reliable eigenvalues without having to check the reliability?
- What is the reason for this unstability? My guess is that the S matrix is close to being not positive definite, but that is just a guess with no argument!

Any suggestions that might enligthen me would be greatly appreciated.

Best,
Troels

[Michael_C_Intel4]

Troels,

when computing this sort of ratio for 5th and 6theigenpairyou're dividing the numbers close to zero by each other in 1st and 4th components - this is unstable and not reliable - I even get Inf on my machine with this division.The normal way tounsure the correctness is to check 2 kinds of ratios:
1) eigenproblem equation: || H*v -E*S*v || should be 0 within error boundaries for every eigenpair
2) orthogonality of the eigenvectors: || V*S*V' - I ||should be 0 within error boundaries.

Michael.

[troels_roennow]

Quoting - Michael Chuvelev (Intel)

Troels,

when computing this sort of ratio for 5th and 6theigenpairyou're dividing the numbers close to zero by each other in 1st and 4th components - this is unstable and not reliable - I even get Inf on my machine with this division.The normal way tounsure the correctness is to check 2 kinds of ratios:
1) eigenproblem equation: || H*v -E*S*v || should be 0 within error boundaries for every eigenpair
2) orthogonality of the eigenvectors: || V*S*V' - I ||should be 0 within error boundaries.

Michael.

Agghhr... yeah. You are completely right. I should have seen that myself!!

Thanks