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>