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[3].
If we compute the q part of the QRDecomposition,
q = First@QRDecomposition[u] ,
then q.ConjugateTranspose[q] == IdentityMatrix[3], 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
readers.
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/>
- References:
- Eigensystem[hermitianMatrix] sometimes returns non-orthonormal
- From: Yen Lee Loh <yloh@mps.ohio-state.edu>
- Eigensystem[hermitianMatrix] sometimes returns non-orthonormal