RE: Eigensystem Error on Machine Precision Matrices
- To: mathgroup at smc.vnet.net
- Subject: [mg22888] RE: [mg22850] Eigensystem Error on Machine Precision Matrices
- From: "Ingolf Dahl" <f9aid at fy.chalmers.se>
- Date: Tue, 4 Apr 2000 01:25:25 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Dear David and Mathgroup,
I am resending this letter, because my computer consumed my letter
instead of sending it. Have you encountered the same error in
Mathematica 4 that I addressed in my message "[mg21057] Eigensystem
applied to a unitary matrix crashes Mathematica 4" from December 1999.
I have had some communication with Wolfram about it, but I have still
not got any solution from them. However, I found a simple fix, that I
posted in [mg21556] in January. The fix seems to work good to me, and
I have tested it for some thousand random samples.
After my letter to the Mathgroup, I read in my old textbook "Numerical
Methods" by Bjork and Dahlquist, and found in section 5.98 a discussion
about the QR-method to solve the eigenvalue problem. They say that it
is a very promising method, but also say that eigenvalues with the same
modulus (absolute value) could give bad convergence. The eigenvalues to
the orthogonal matrices all have modulus one. They also suggest the
use of shift parameters, similar to my fix. They suggest the use of
different shift parameters in each iteration, and maybe that is a good
idea. So someone has thought about these problems a long time ago.
By the way, the way you generate random orthogonal matrices (with
random integers) is a bit rude, but probably quite Ok for generating
some test matrices. The use of the vector {0,1,0} as second vector in
the matrix will after the GramSchmidt procedure still give a bias for
this vector, so the probability distribution of the orthogonal matrices
will not be uniform. A more sophisticated and mathematical sound way to
generate the random orthogonal matrices is to use coordinates taken as
random numbers from the Normal Distribution. That will result in
vectors with a distribution that is invariant to rotations of the
coordinate system, and thus uniform over the unit sphere, in three or N
dimensions, whatever you want. Here is the code, n is the dimension:
RandomOrthogonalMatrix[n_] := Module[{deter, a},
Needs["LinearAlgebra`Orthogonalization`"];
Needs["Statistics`NormalDistribution`"];
deter = 0.;
While[Abs[deter] < 0.5,
a =
LinearAlgebra`Orthogonalization`GramSchmidt[
Table[Random[
Statistics`NormalDistribution`NormalDistribution[0,1]], {n},
{n}],
LinearAlgebra`Orthogonalization`Normalized -> True];
deter = Det[a]];
a]
I hope that someone will find it useful.
Sincerely yours
Ingolf Dahl
Chalmers University
f9aid at fy.chalmers.se