Eigensystem Error on Machine Precision Matrices
- To: mathgroup at smc.vnet.net
- Subject: [mg22850] Eigensystem Error on Machine Precision Matrices
- From: "David Park" <djmp at earthlink.net>
- Date: Sat, 1 Apr 2000 02:51:10 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
Dear MathGroup,
I believe that the Eigensystem command in Version 4 gives errors about 10%
of the time when used on machine precision orthogonal matrices. It does not
give errors when extended precision is used for the matrices, or if the
calculation is enclosed in SetPrecision, or if the matrix has exact values.
These last methods, however, take 10 to 20 times as long. The machine
precision results have complex eigenvectors that differ from the proper
results (by more that just a complex factor) and they do not satisfy the
eigenvalue equation.
In the following I give one particular test case which gives an incorrect
answer, and statements to evaluate various cases. I also give the code to
generate and test random orthogonal matrices. To keep the posting short, I
do not show the results.
Matrix which gives an error:
omat = {{2/Sqrt[13], 0, 3/Sqrt[13]}, {3/Sqrt[182],
Sqrt[13/14], -Sqrt[2/91]},
{-(3/Sqrt[14]), 1/Sqrt[14], Sqrt[2/7]}};
The exact solution and test:
{evalse, evectse} = Eigensystem[omat] // N // Chop;
Chop[omat.Transpose[evectse]] ==
Chop[Transpose[evectse].DiagonalMatrix[evalse]]
Display of eigenvalues and eigenvectors:
{MatrixForm[Transpose[{evalse}]], MatrixForm[evectse]}
Display of the two sides of the eigenvalue equation:
{omat.Transpose[evectse] // Chop // MatrixForm,
Transpose[evectse].DiagonalMatrix[evalse] // Chop // MatrixForm}
Evaluation with a machine precision matrix (gives incorrect answer):
{evals, evects} = Eigensystem[N[omat]] // Chop;
Chop[omat.Transpose[evects]] ==
Chop[Transpose[evects].DiagonalMatrix[evals]]
Evaluation with extended precision (gives correct answer):
{evalse, evectse} = Eigensystem[N[omat, 20]];
Chop[omat.Transpose[evectse]] ==
Chop[Transpose[evectse].DiagonalMatrix[evalse]]
Evaluation within SetPrecision (gives correct answer):
{evalse, evectse} = SetPrecision[Eigensystem[omat], 5];
Chop[omat.Transpose[evectse]] ==
Chop[Transpose[evectse].DiagonalMatrix[evalse]]
This is code to generate exact random orthogonal matrices and then test
them:
Needs["LinearAlgebra`Orthogonalization`"]
randomortho::usage =
"randomortho[1] with generate a random orientation-preserving orthogonal
\
matrix. randomortho[-1] will generate a random orientation-reversing \
orthogonal matrix.";
randomortho[orientation_?(MatchQ[#, -1 | 1] &)] :=
Module[{v0, mat},
gen := (v0 = {Random[Integer, {-10, 10}], Random[Integer, {-10, 10}],
Random[Integer, {-10, 10}]};
mat = GramSchmidt[{v0, {0, 1, 0}, {0, 0, 1}}]);
gen; While[Det[mat] == 0, gen];
If[Det[mat] == orientation, mat, Reverse[mat]]
]
omat = randomortho[1];
{evals, evects} = Eigensystem[N[omat]] // Chop;
Chop[omat.Transpose[evects]] ==
Chop[Transpose[evects].DiagonalMatrix[evals]]
I would appreciate if some of the greater mathematicians in the group could
confirm if this is indeed a bug, or else inform me if there is something
about orthogonal matrices and eigensystems that I don't understand.
David Park
djmp at earthlink.net
http://home.earthlink.net/~djmp/