Re: Eigensystem applied to a unitary matrix crashes Mathematica 4.
- To: mathgroup at smc.vnet.net
- Subject: [mg21556] Re: Eigensystem applied to a unitary matrix crashes Mathematica 4.
- From: "Ingolf Dahl" <f9aid at fy.chalmers.se>
- Date: Sat, 15 Jan 2000 02:04:00 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
Some further comments on Eigensystem applied to unitary matrices:
I have written a short program, which tests the Eigensystem on randomly
chosen real unitary 3x3 matrices, and then find a failure in approx 18 % of
the cases. Randomised real Hermitian 3x3 matrices or randomised arbitrary
real 3x3 matrices give zero failures for ten thousand trials. The failure
rate for 5x5 unitary matrices is 44%, for 7x7 unitary matrices 60% and for
9x9 unitary matrices 73% for one thousand trials. If I try with 4x4, 6x6 or
8x8 unitary matrices, the Mathematica kernel crashes after a few trials.
This failure for unitary matrices might be special for Mathematica 4.0.0.0
under Windows 98, I have not checked. Unitary matrices are sometimes also
called orthogonal matrices, and have the property that their inverse is
equal to the transpose of the matrix. The eigenvalues of a unitary matrix is
always ±1 or complex conjugate pairs of modulus one. Unitary matrices have
several uses, and it is a bit of a catastrophe (at least for me) if
Mathematica fails on extracting the correct eigenvectors. As a mathematical
problem, the eigenvalue problem for a unitary matrix should in general be
numerically stable.
However, I have found a quite simple ad hoc fix, that seem to work quite
good: I add the identity matrix times the number 0.5647382, and then apply
Eigensystem to that matrix. The matrix submitted to Eigensystem is then no
longer unitary, and thus not subject to the same failure mechanism. Then I
just subtract the same number 0.5647382 from the obtain eigenvalues. And
voila! It works!
I have tested the fix for 1000 random cases of unitary matrices for 3x3, 4x4
up to 9x9 matrices, but it also seem to work for a 100x100 random unitary
matrix. A 200x200 random matrix is too heavy for my computer. How to
generate these random matrices is another story (please ask if you are
interested!).
Here is the fix as a Mathematica routine, that :
eigensystemfixed[matrix_] :=
Module[{vals, vects}, {vals, vects} =
Eigensystem[matrix + 0.5647382*IdentityMatrix[Length[matrix]]];
{vals - 0.5647382, vects}]
Why the number 0.5647382? You could probably choose some other queer number,
but I just thought it was improbable that anyone will apply this function to
any matrix that becomes unitary when you subtract the identity matrix times
the number 0.5647382 from it. Here I can give a suggestion about a dirty
trick to Wolfram: if this fix is incorporated in the compiled version of
Eigensystem, using a appropriate random number instead of 0.5647382, the
Eigensystem routine will just fail in some rare, non-reproducible cases, and
the users will not be able to complain. A cleaner and nicer method could be
to check the result for failures and then apply this fix as a second trial
whenever that is needed.
Sincerely yours
Ingolf Dahl
Chalmers University
f9aid at fy.chalmers.se