Re: Mathematica 6 obtains imaginary eigenvalues for a Hermitian
- To: mathgroup at smc.vnet.net
- Subject: [mg86174] Re: [mg86010] Mathematica 6 obtains imaginary eigenvalues for a Hermitian
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Tue, 4 Mar 2008 02:10:30 -0500 (EST)
- References: <200803010947.EAA23583@smc.vnet.net>
Sebastian Meznaric wrote:
> I have a 14x14 Hermitian matrix, posted at the bottom of this message.
> The eigenvalues that Mathematica obtains using the
> N[Eigenvalues[matrix]] include non-real numbers:
> {-9.41358 + 0.88758 \[ImaginaryI], -9.41358 -
> 0.88758 \[ImaginaryI], -7.37965 + 2.32729 \[ImaginaryI], -7.37965 -
> 2.32729 \[ImaginaryI], -4.46655 + 2.59738 \[ImaginaryI], -4.46655 -
> 2.59738 \[ImaginaryI], 4.36971, 3.21081, -2.32456 +
> 2.10914 \[ImaginaryI], -2.32456 - 2.10914 \[ImaginaryI],
> 2.04366+ 0.552265 \[ImaginaryI],
> 2.04366- 0.552265 \[ImaginaryI], -0.249588 +
> 1.29034 \[ImaginaryI], -0.249588 - 1.29034 \[ImaginaryI]}.
> However, if you do Eigenvalues[N[matrix]] it obtains different results
> {-9.09122, -7.41855, -7.41855, -7.2915, 4.33734, -4., -4., 3.2915, \
> -3.24612, -2.38787, -2.38787, 1.80642, 1.80642, 0}.
>
> These results agree with Solve[CharacteristicPolynomial[matrix,x],x].
> Therefore I assume that the latter are correct. Has anyone seen this?
> I am using 6.0.0.
>
>
> Here is the matrix:
> {{-6, 0, -Sqrt[3], 0, 0, Sqrt[3], 0, 0, 0, 0, 0, 0, 0, 0}, {0, -6,
> 0, -Sqrt[3], 0, 0, Sqrt[3], 0, 0, 0, 0, 0, 0, 0}, {-Sqrt[3], 0, -4,
> 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], 0, 0, Sqrt[3], 0,
> 0, 0, 0, 0, 0}, {0, -Sqrt[3],
> 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -4/3, -(2 Sqrt[2])/3, 0, 0, 0,
> Sqrt[3], 0, 0, 0, 0, 0}, {0, 0, 2 Sqrt[2/3], -(2 Sqrt[2])/3, 7/3, 0,
> 0, 0, 0, Sqrt[3], 0, 0, 0, 0}, {Sqrt[3], 0, 0, 0, 0, -4, 0,
> 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, 0, 2 Sqrt[2/3], 0, 0, 0}, {0,
> Sqrt[3], 0, 0, 0, 0, -4, 0, 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, 0,
> 2 Sqrt[2/3], 0, 0}, {0, 0, Sqrt[3], 0, 0,
> 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, -14/3,
> 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], (2 Sqrt[2])/3, 0,
> 0, 0}, {0, 0, 0, Sqrt[3], 0, 0, 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4),
> 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -2, -(2 Sqrt[2])/3, 0, (
> 2 Sqrt[2])/3, 0, 0}, {0, 0, 0, 0, Sqrt[3], 0, 0,
> 2 Sqrt[2/3], -(2 Sqrt[2])/3, -7/3, 0, 0,
> 2 (1/(3 Sqrt[2]) + (2 Sqrt[2])/3), Sqrt[10/3]}, {0, 0, 0, 0, 0,
> 2 Sqrt[2/3], 0, (2 Sqrt[2])/3, 0, 0, -16/3,
> 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], 0}, {0, 0, 0, 0, 0,
> 0, 2 Sqrt[2/3], 0, (2 Sqrt[2])/3, 0,
> 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -8/3, -(2 Sqrt[2])/3, 0}, {0, 0,
> 0, 0, 0, 0, 0, 0, 0, 2 (1/(3 Sqrt[2]) + (2 Sqrt[2])/3),
> 2 Sqrt[2/3], -(2 Sqrt[2])/3, 1/2,
> 2 (-Sqrt[5/3]/16 - Sqrt[15]/16)}, {0, 0, 0, 0, 0, 0, 0, 0, 0, Sqrt[
> 10/3], 0, 0, 2 (-Sqrt[5/3]/16 - Sqrt[15]/16), 7/2}}
I wanted to confirm that this is a bug in the symbolic eigenvalue code.
It was fixed in our development version back in December; the fix was
deemed too risky to try to backport to version 6.0.2 (which was very
near to completion at that time).
The problem is a bit obscure and has to do with some details of row
reduction when a matrix contains algebraics. It is not obvious that the
problem lies there, because CharacteristicPolynomial gets the correct
polynomial. What happens is that Eigenvalues first shuffles rows and
columns (symmetrically) according to some heuristic, and then attempts
to compute the characteristic polynomial of this modified matrix. It
should agree with CharacteristicPolynomial of the original, up to a
factor of +-1. But it does not, because with the shuffled matrix the row
reduction used to find that char poly runs afoul of said bug.
Here we have IdentityMatrix[14]*x - shuffledmatrix.
m2 = {{6 + x, Sqrt[3], 0, 0, 0, 0, 0, -Sqrt[3], 0, 0, 0, 0, 0, 0},
{Sqrt[3], 4 + x, -2*(-1/(4*Sqrt[3]) + (3*Sqrt[3])/4), 0, 0,
0, -Sqrt[3], 0, 0, 0, 0, 0, 0, -2*Sqrt[2/3]},
{0, -2*(-1/(4*Sqrt[3]) + (3*Sqrt[3])/4), 4/3 + x, Sqrt[3], 0,
-Sqrt[3], 0, 0, 0, 0, 0, 0, 0, (2*Sqrt[2])/3},
{0, 0, Sqrt[3], 6 + x, -Sqrt[3], 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, -Sqrt[3], 4 + x, -2*(-1/(4*Sqrt[3]) + Sqrt[3]/4), 0, 0,
0, -2*Sqrt[2/3], 0, 0, 0, 0},
{0, 0, -Sqrt[3], 0, -2*(-1/(4*Sqrt[3]) + Sqrt[3]/4), 2 + x,
-2*(-1/(4*Sqrt[3]) + (3*Sqrt[3])/4), 0, 0, (-2*Sqrt[2])/3, 0,
(2*Sqrt[2])/3, 0, 0},
{0, -Sqrt[3], 0, 0, 0, -2*(-1/(4*Sqrt[3]) + (3*Sqrt[3])/4), 14/3 + x,
-2*(-1/(4*Sqrt[3]) + Sqrt[3]/4), (-2*Sqrt[2])/3, 0, 0,
-2*Sqrt[2/3], 0, 0},
{-Sqrt[3], 0, 0, 0, 0, 0, -2*(-1/(4*Sqrt[3]) + Sqrt[3]/4), 4 + x,
-2*Sqrt[2/3], 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, (-2*Sqrt[2])/3, -2*Sqrt[2/3], 16/3 + x,
-2*(-1/(4*Sqrt[3]) + (3*Sqrt[3])/4), -2*Sqrt[2/3], 0, 0, 0},
{0, 0, 0, 0, -2*Sqrt[2/3], (-2*Sqrt[2])/3, 0, 0,
-2*(-1/(4*Sqrt[3]) + (3*Sqrt[3])/4), 8/3 + x, (2*Sqrt[2])/3, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, -2*Sqrt[2/3], (2*Sqrt[2])/3, -1/2 + x,
-2*(1/(3*Sqrt[2]) + (2*Sqrt[2])/3), -2*(-Sqrt[5/3]/16 -
Sqrt[15]/16), 0},
{0, 0, 0, 0, 0, (2*Sqrt[2])/3, -2*Sqrt[2/3], 0, 0, 0,
-2*(1/(3*Sqrt[2]) + (2*Sqrt[2])/3), 7/3 + x, -Sqrt[10/3], -Sqrt[3]},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2*(-Sqrt[5/3]/16 - Sqrt[15]/16),
-Sqrt[10/3],
-7/2 + x, 0},
{0, -2*Sqrt[2/3], (2*Sqrt[2])/3, 0, 0, 0, 0, 0, 0, 0,
0, -Sqrt[3], 0, -7/3 + x}};
This is the computation that gets the characteristic polynomial wrong.
Expand[Det[m2]]
This has been wrong since version 5.1, although I suspect the underlying
bug predates that and some small change (pivot selection, perhaps)
brought it out in this example at that time.
For the record, this was one of the longest bug hunts I ever set out on.
It involved a subtlety in a method known as "one-step" row reduction.
Subduing the creature required a wooden stake, a silver bullet, and
several incantations from bug-killing lore that are not repeatable.
Daniel Lichtblau
Wolfram Research
- References:
- Mathematica 6 obtains imaginary eigenvalues for a Hermitian matrix
- From: Sebastian Meznaric <meznaric@gmail.com>
- Mathematica 6 obtains imaginary eigenvalues for a Hermitian matrix