       Re: Eigensystem sometimes returns non-orthonormal

• To: mathgroup at smc.vnet.net
• Subject: [mg96874] Re: [mg96805] Eigensystem sometimes returns non-orthonormal
• From: Yen Lee Loh <yloh at mps.ohio-state.edu>
• Date: Thu, 26 Feb 2009 07:59:03 -0500 (EST)
• References: <200902250904.EAA15732@smc.vnet.net> <49A55D8A.6080004@wolfram.com>

Dear Daniel,

Thanks for your suggestion.  Unfortunately it doesn't fix the problem:

h = [a certain 3x3 hermitian matrix]
{u, d} = Eigensystem[h];   d=DiagonalMatrix[d]

The eigenvalue equation is satisfied, as per the specifications:
h.Transpose[u] - Transpose[u].d == 0.
The eigenvectors are normalized, but not orthonormal:
u.ConjugateTranspose[u] != IdentityMatrix.
If we compute the q part of the QRDecomposition,

q = First@QRDecomposition[u]  ,

then q.ConjugateTranspose[q] == IdentityMatrix, but h.Transpose[q] -
Transpose[q].d ==0,
so the eigenvalue equation is no longer satisfied.   q =Transpose@First
@QRDecomposition[Transpose@u] doesn't work either.

I haven't analysed things in detail, but since the task at hand is to
"orthogonalize only the sets of eigenvectors u corresponding to degenerate
eigenvalues", it may not be safe to use an approach
(Orthogonalize/QRDecomposition) that doesn't "know" about the pattern of
eigenvalues.  Of course one can write a program beginning with
degenEigvalIndices=Split[Sort[d]]
and then calling Orthogonalize within each set, but this is kinda messy.
I'm wondering if it is actually easier to do the calculation in Fortran/C++
and use MathLink, somehow.  Does Wolfram have any plans to add support for
orthonormal eigenvectors?

Out of interest, does anyone know WHY Eigensystem can returns
non-orthonormal eigenvectors?  I was under the impression that Eigensystem[]
first checks to see if the matrix is Hermitian, and if it is, it uses
LAPACK-based algorithms for Hermitian matrices.  I thought LAPACK's
algorithms involved tridiagonalization followed by some form of QR
factorization/divide-and-conquer Givens rotations, in which the eigenvector
matrix started off as the identity matrix and accumulated successive unitary
transformations, thus remaining unitary at every stage;  I'm curious how
Eigensystem[] even manages to break orthonormality at all.

p.s.  DrMajorBob: Thanks for your suggestion, but my eigenvectors are
already normalized --- the problem is that they're not mutually orthogonal,
which seems to be trickier to fix.
p.p.s.  Sorry this post is rather badly written.  Maybe when the final
solution appears I will summarize things properly for the benefit of future

On Wed, Feb 25, 2009 at 10:02 AM, Daniel Lichtblau <danl at wolfram.com> wrote:

> Yen Lee Loh wrote:
>
>> Dear Mathematica users,
>>
>> In Mathematica 7.0.0 (for Linux), calling Eigenvectors[H] or
>> Eigensystem[H]
>> for a numerical Hermitian matrix H sometimes returns eigenvectors that are
>> not orthonormal.
>> This happens when some eigenvalues are degenerate.  (I can supply example
>> code that illustrates the problem, if necessary.)
>>
>> This is not really a bug -- the documentation for Eigensystem[] doesn't
>> make
>> any guarantees of orthonormality -- but nevertheless it is an annoying
>> part
>> of Mathematica's design.
>> This issue has been raised 11 years ago (
>> http://forums.wolfram.com/mathgroup/archive/1998/Mar/msg00418.html ), but
>> that post is corrupted!  (Surely Wolfram isn't resorting to censorship?)
>>
>
> A non-corrupted version is located at the URL below.
>
> http://forums.wolfram.com/mathgroup/archive/1998/Mar/msg00425.html
>
> Note that the behavior in that report was of a more serious nature than
> that which you describe (eigenvectors sometimes had complex values).
>
>
>  I
>> was hoping that in Mathematica 7 I would be able to write something like
>>
>>    Eigensystem[H, Method->"LAPACK-ZHEEVR"]
>>
>> or
>>
>>    Eigensystem[H, OrthonormalizeEigenvectors->True]
>>
>> but no such options seem to exist.  One workaround is to apply
>> Orthogonalize[] to the matrix of eigenvectors, but the documentation for
>> Orthogonalize[] doesn't guarantee that the orthonormalization will only
>> occur within the "degenerate subspace".  So one has to resort to
>> complicated
>> fixes (e.g., http://arxiv.org/pdf/hep-ph/9607313 ).
>> Does anyone have a simpler solution?  (For example, is there an easy way
>> to
>> call LAPACK'S ZHEEVR routine, which guarantees orthornormal eigenvectors,
>> from Mathematica?)
>>
>> Thanks a lot.
>> Yen Lee Loh
>>
>
> If you apply QRDecomposition to the eigenvectors, then I believe the Q part
> will provide what you are looking for.
>
>
> Daniel Lichtblau
> Wolfram Research
>
>
>

--
Yen Lee Loh
Postdoctoral Associate, The Ohio State University
Home: 544 Stinchcomb Dr Apt 10, Columbus OH 43202-1728, USA
Office: 2043 Physics Research Building, 191 W Woodruff Ave, Columbus OH
43210-1117, USA
Office phone: +1 614 247 4772
Mobile phone: +1 765 532 9457
Email: yloh at mps.ohio-state.edu
Web: http://www.physics.ohio-state.edu/~yloh/<http://www.physics.ohio-state.edu/%7Eyloh/>

• Prev by Date: Re: Graphics export for high quality documents
• Next by Date: Re: Re: listing user defined, independent variables
• Previous by thread: Re: Eigensystem sometimes returns non-orthonormal
• Next by thread: Re: Eigensystem[hermitianMatrix]sometimes returns