LAPACK fails calculating eigenvectors

Open discussion regarding features, bugs, issues, vendors, etc.

LAPACK fails calculating eigenvectors

Postby emilio » Mon Mar 12, 2012 11:49 am

I wrote some code to solve the general eigenvalue problem and now I am comparing my results agains LAPACK, 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.
emilio
 
Posts: 2
Joined: Fri Feb 24, 2012 7:18 am

Re: LAPACK fails calculating eigenvectors

Postby cottrell » Mon Mar 19, 2012 8:30 pm

It should be obvious that you have not provided enough information
for anyone to assess your claim or offer any solution or advice.

What is the original matrix that is the subject of your eigen-analysis?
cottrell
 
Posts: 68
Joined: Thu Jan 15, 2009 1:40 pm

Re: LAPACK fails calculating eigenvectors

Postby Julien Langou » Tue Mar 20, 2012 12:19 am

A link is missing from Horacio's post. The example Horacio was referring to is: http://www.nag.com/lapack-ex/node102.html
The same post is also at: http://stackoverflow.com/questions/9669736/lapack-fails-calculating-eigenvectors

Most of the conversation happened by email between Sven and Horacio. I copy-paste for further references.

Cheers, Julien.

Dear Horacio,

[... ]
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.

[ ... ]

DSPGVX will have computed the solution in double precision,
but the web site has rounded the solution to four decimal figures for
convenience of printed output. If I take the solution direct from
DSPGVX then on my machine I get residuals

norm(r1) = 5.4177e-16, norm(r2) = 5.3663e-16,

so DSPGVX is indeed computing satisfactory results.

Best wishes,

Sven.
Julien Langou
 
Posts: 734
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: LAPACK fails calculating eigenvectors

Postby error5772 » Thu Mar 29, 2012 7:43 am

Hi emilio!
You can not trust both and my comment to the generalized
eigensystem problem is: Avoid it!

Nevertheless you should construct your own test-matrices from
matrices with known eigenvalues and eigenvectors! In the case of
a known eigensystem of a symmetric real matrix C you can write
it as: D = U^T C U <=> C = U D U^T, where U contains the
eigenvectors vec_k as columns in the same order the diagonal
matrix D contains the eigenvalues lambda_k as diagonal elements
and U^TU = 1 for normalized eigenvectors! Now look at:
C vec_k = lambda_k vec_k <=> U D U^T vec_k = lambda_k vec_k,
which can be converted to: D U^T vec_k = lambda_k U^T vec_k
and identified with: A vec_k = lambda_k B vec_k,
where: A = D U^T and: B = U^T. So the same eigenvectors and
eigenvalues fullfill this generalized eigensystem problem and you
can test the routine! You can try the known eigensystem of the real
tridiogonal matrix, where all diagonal elements are a and the direct
neighbor elements are b and all else zero:
ab000
bab00
0bab0
.........
of dimension N with eigenvalues: lambda_k = a + 2*b*Cos(k*Pi/(N+1))
and eigenvectors consisting of the components: (vec_k)_i = Sqrt(2/(N+1))*Sin(i*k*Pi/(N+1))
for i=1,2,...,N where each pair has a fixed k=1,2,...,N.

This is a very simple example for testing the normal eigensystem
routine (dsyev() or dsyevx()), but if you want to test the generalized
eigensystem routine I am not sure if the routine gets problems with
zero eigenvalues or eigenvalues getting close to zero!
You should also check the quality of the LU decomposition
of your matrix B, because this is needed for the algorithm
(Or not?).

Please report the results if you read and try it!

Good luck!

Michael
error5772
 
Posts: 19
Joined: Thu Sep 01, 2011 5:02 am


Return to User Discussion

Who is online

Users browsing this forum: Yahoo [Bot] and 2 guests

cron