orthonormal version of Eigensystem[]?

[John Sidles]

Dear Mathematica users

I admire the version of Eigensystem[] on Mathematica very much. It is
fast and accurate, even for the 1000 x 1000 symmetric real matrices
than I am looking at.

But Eigensystem[] has one very important limitation when solving a
symmetric real matrix, or more generally when solving a complex
Hermitian matrix. This limitation is poorly documented and can "bite"
the unwary user, as I found from personal experience.

Namely, although Eigensystem[] usually returns a set of orthonormal
eigenvectors when given a real symmetric matrix, it is not *guaranteed*
do so. Specifically in the case where two eigenvalues are degenerate,
the returned eigenvectors are (occasionally but not always) grossly
nonorthogonal.

Now, this behavior is perfectly consistent with the documentation of
Eigensystem[], which is completely silent on the issue of eigenvector
orthonormality. The reason the behavior is a problem for a typical
Mathematica user is that the eigenvectors of a symmetric real matrix
can always be chosen to be orthonormal, and *usually* the vectors
returned by Eigensystem[] *are* orthonormal, and unwary users are
therefore be led to assume that this is always the case. Not so!

Now, there is a fix for this within Mathematica, but the fix itself has
a problem. Namely, we can always compute a set of orthonormal
eigenvectors by the following simple commands:

<< LinearAlgebra`Orthogonalization`;
m := (* some real symmetric matrix *);
vecs = Normalize/@(Eigenvectors[m]//GramSchmidt);

This fix is guaranteed to work, but has the problem that it is hideously
inefficient. For big matrices, the Gram-Schmidt orthogonalization
takes orders of magnitude more time than Eigenvectors[] (even though it
shouldn't from a technical point of view).

So as far as I can tell, therefore, the problem is unresolved: how can
the eigenvalues and orthonormal eigenvectors of large Hermition
matrices be efficiently computed using Mathematica? And where on
MathSource can I find better information?

Mathematica employees: please mention this topic on the Mathematica web
site FAQ in the area on matrix computations (presently "under
construction").

Thanks ... JAS

PS: Here is a function called "orthonormalEigensystem[]" which takes
about 25% longer to evaluate than "Eigensystem[]" for large matrices,
but which is guaranteed to return sorted eigenvalues together with
orthonormal eigenvectors. Note that it accepts only real numerical
matrices as input, and keeps only the symmetric part of such matrices,
and is in this respect more specialized than "Eigensystem[]".

<< LinearAlgebra`Orthogonalization`;

realQ[x_] := NumberQ[x] && (Im[x] == 0.0);

sortEigs[x_] := x//Transpose//Sort//Transpose ;

quickerGramSchmidt[x_] := Module[
{err,eout,len,istart,iend},

(* eigenvalues closer than "err" are treated as degenerate *)
err = 10^6*$MachineEpsilon*(x[[1]]//Abs//Max);

(* create an empty list to store the results *)
eout = {{},{}};
len = x[[1]]//Length;

(* loop over starting indices of blocks of degenerate eigenvalues *)
istart=1;
While[istart<=len,

(* locate the end of the current degenerate block *)
iend = istart;
While[ iend<len &&
(x[[1,istart]]-x[[1,iend+1]]//Abs) < err,iend++];

(* end of degenerate block has been found, so process it *)
If[iend == istart,

(* normalize the eigenvectors of isolated eigenvalues *)
AppendTo[eout[[1]],{x[[1,istart]]}];
AppendTo[eout[[2]],{x[[2,istart]]//Normalize}];,

(* orthogonalize the eigenvectors of degenerate eigenvalues
*)
AppendTo[eout[[1]],Take[x[[1]],{istart,iend}]];
AppendTo[eout[[2]],
Normalize/@Take[x[[2]],{istart,iend}]//GramSchmidt]];

(* look for the next block of eigenvectors *)
istart = iend+1;];

(* flatten the created list *)
{Flatten[eout[[1]],1],Flatten[eout[[2]],1]} ];

orthonormalEigensystem[x_] :=(Eigensystem[0.5*(x+(x//Transpose))]//
sortEigs// quickerGramSchmidt)/;MatrixQ[x,realQ];

orthonormalEigensystem::usage :=
"orthonormalEigensystem[] is similar to the Mathematica function \
Eigensystem[] except:(1) only real numerical matrices are accepted, (2)
only \ the symmetric part of the input matrix is kept for computation,
and (3) a \ sorted list of eigenvalues and their orthonormal
eigenvectors is returned. \ Smallest eigenvalues are listed first.
Note: the Mathematica function \ Eigensystem[] sometimes returns
non-orthonormal eigenvectors. In contrast, \ the eigenvectors returned
by orthonormalEigensystem[] are guaranteed to be \ orthonormal.";