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

**NDSolve: Distinguishing equations**

**Re: Problem with evaluation of delayed rules**

**Re: Eigensystem applied to a unitary matrix crashes Mathematica 4.**

**Re: Factor MatrixForm**