evecs not correct with jacobi, tred2, tqli


Ward Ciac
03-20-2007, 01:27 PM
I am trying to calculate evals and evecs of a real symmetric matrix - thought I should use either jacobi or tred2-tqli with Fortran 77. Originally I used NR v. 2.01 then upgraded the subroutines to v. 2.08 in case there were any errors on the earlier version. Running on my IBM X60 laptop (Centrino Duo) with XP Pro and Intel® Visual Fortran Compiler 9.0 for Windows, IA-32 processor. I am just doing 3 x 3 matrices with known answers. I obtained these from Schaum's Outline Series with answers and confirmed the answers at link http://www.arndt-bruenner.de/mathe/scripts/engl_eigenwert.htm. The evals come out correctly but not always the evecs. In the NR book, I notice a caution about 'underflow' - not sure what this is - might this be a problem? I did change all reals to double precision.

Jacobi and tred2-tqli calculate some 3x3's correctly, for example:

(3 2 -1)
(2 3 -1)
(-1 -1 4)

has correct evals (1 3 6) and evecs (1,-1,0), (1 1 2), and (1 1 -1), which with jacobi and tred2-tqli comes out as:
evals (1 6 3) and evecs (.707,-.707,0), (.577,.577,-.577), and (-.408,-.408,-.816)
(just a difference in order and scaling for evals/evecs).

However, for example for the following matrix the evecs are not correct with jacobi or tred2-tqli:

(3,-1,1)
(-1,3,-1)
(1,-1,3)

has correct evals (2 2 5) and correct evecs (1 1 0), (-1 0 1), (1 -1 1) but with tred2-tqli:

evals (2 5 2), evecs (.707,.707,0), (-.577,.577,-.577), (.408,-.408,-.816) or scaled this is evecs (1,1,0), (-1,1,-1), (,1,-1,-2)
which is not the same.
With jacobi:
evals (5 2 2 ), evecs: (.577,-.577,.577), (.789,.211,-.577), (.211,.788,.577) which is not the same.

Another example:

(1,1,-1)
(1,1,-1)
(-1,-1,1)

has correct evals (0 0 3) and correct evecs (-1 1 0), (1 0 1), (-1 -1 1)

but with tred2-tqli:
evals (0 3 0), evecs (.707 -.707 0), (.577 .577 -.577), (-.408 -.408 -.816)

with jacobi:
evals (0 3 0), evecs (.789 -.211 .577), (.577 .577 -.577), (-.211 .789 .577)

Thanks for listening!
Ward Ciac

Saul Teukolsky
03-20-2007, 07:11 PM
Hi Ward,

If two or more eigenvalues are the same, then any linear combination of the corresponding eigenvectors is also an eigenvector. You can check that the vectors you claimed are wrong are in fact correct.

Saul Teukolsky

Ward Ciac
03-21-2007, 03:38 PM
Dear Dr. Teukolsky,

I think I see my mistake. Two eigenvectors corresponding to two different eigenvalues are mutually orthogonal. However, the eigenvectors x of degenerate eigenvalues lambda are not necessarily mutually orthogonal

(ref: http://farside.ph.utexas.edu/teaching/336k/lectures/node72.html)

although they satisfy the matrix eigenvalue equation:

A x = lambda x

For the degenerate eigenvalues, the eigenvectors given in the other references were not mutually orthogonal – as can be observed visually by plotting their Cartesian coordinates in 3D, such as with map3d:

http://www.sci.utah.edu/cibc/software/map3d.html

The NR subroutines transform the eigenvectors of degenerate eigenvalues so that they are mutually orthogonal. That is very nice, because I need orthogonal eigenvectors for two of our projects which use Karhunen-Loève transformation: face recognition, reference:

M. Turk and A. Pentland, "Eigenfaces for recognition", J Cognitive Neuroscience 1991;3:71–86

and Fukunaga-Koontz transformation, reference:

K. Fukunaga and W.L.G. Koontz, "Application of the Karhunen-Loève Expansion to Feature Selection and Ordering", IEEE Trans Comp C-19, pp 311-318, 1970

I think it also answers the question posed in a thread in this forum by Denis Davidoff ‘Eigenproblem (Jacobi) multiplicity>1’.

Thank you!
Ward