Mathematica 6 obtains imaginary eigenvalues for a Hermitian matrix

[Sebastian Meznaric]

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}}

[Andrzej Kozlowski]

In fact, for your matrix exact computation with Eigenvalues gives a
wrong answer. Thus this is not a numerical problem, as often happens
in such cases, but something much more basic.

mat = {{-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}};

Take the "largest" eigenvalue:

eig = First[Eigenvalues[mat, 1]];

f[x_] = CharacteristicPolynomial[mat, x]

N[f[eig], 10]
-1.8696095421628027977`10.111112547559953*^10 -
8.339392668664626899`9.760496055115645*^9*I

which is far from zero.

Another related point is that setting the Options Cubics->True should
give a radical expression, in cases when the Characteristic Polynomial
can be factored into linear factors, quadratics and quadratics as is
the case here:

Factor[f[x]]
x*(x + 4)^2*(x^2 + 4*x - 24)*(x^3 + 8*x^2 - 32)^2*(x^3 + 8*x^2 - 24*x
- 128)

But actually using this options makes no difference. It looks very
much like Eigenvalues is solving the wrong polynomial of the same
degree since

g = MinimalPolynomial[eig, x]

17006112*x^14 + 612220032*x^13 + 8162933760*x^12 + 40814668800*x^11 -
81629337600*x^10 - 1671768834048*x^9 -
4667021328384*x^8 + 11702381838336*x^7 + 76901366366208*x^6 +
42669909636920*x^5 - 337697299646488*x^4 -
551524516377720*x^3 + 576589503454408*x^2 +
260729243794329*x + 2600447009874406

This is a prime polynomial and it is not clear how, if at all, it is
related to the characteristic polynomial.

Andrzej Kozlowski

[David Reiss]

Note that your matrix is not invertible (its determinant is zero). So
this is the source of your problem...

Hope that this helps...

-David
A WorkLife FrameWork
E x t e n d i n g MATHEMATICA's Reach...
http://scientificarts.com/worklife/

On Mar 1, 4:57=A0am, Sebastian Meznaric <mezna...@gmail.com> 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 -

> =A0 0.88758 \[ImaginaryI], -7.37965 + 2.32729 \[ImaginaryI], -7.37965 -
> =A0 2.32729 \[ImaginaryI], -4.46655 + 2.59738 \[ImaginaryI], -4.46655 -
> =A0 2.59738 \[ImaginaryI], 4.36971, 3.21081, -2.32456 +
> =A0 2.10914 \[ImaginaryI], -2.32456 - 2.10914 \[ImaginaryI],
> =A02.04366+ 0.552265 \[ImaginaryI],
> =A02.04366- 0.552265 \[ImaginaryI], -0.249588 +
> =A0 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,

> =A0 0, -Sqrt[3], 0, 0, Sqrt[3], 0, 0, 0, 0, 0, 0, 0}, {-Sqrt[3], 0, -4,
> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], 0, 0, Sqrt[3], 0,
> =A0 0, 0, 0, 0, 0}, {0, -Sqrt[3],
> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -4/3, -(2 Sqrt[2])/3, 0, 0, 0,
> =A0 Sqrt[3], 0, 0, 0, 0, 0}, {0, 0, 2 Sqrt[2/3], -(2 Sqrt[2])/3, 7/3, 0,
> =A0 =A00, 0, 0, Sqrt[3], 0, 0, 0, 0}, {Sqrt[3], 0, 0, 0, 0, -4, 0,
> =A0 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, 0, 2 Sqrt[2/3], 0, 0, 0}, {0,
> =A0 Sqrt[3], 0, 0, 0, 0, -4, 0, 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, 0,
> =A0 2 Sqrt[2/3], 0, 0}, {0, 0, Sqrt[3], 0, 0,
> =A0 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, -14/3,
> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], (2 Sqrt[2])/3, 0,
> =A0 0, 0}, {0, 0, 0, Sqrt[3], 0, 0, 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4),
> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -2, -(2 Sqrt[2])/3, 0, (
> =A0 2 Sqrt[2])/3, 0, 0}, {0, 0, 0, 0, Sqrt[3], 0, 0,
> =A0 2 Sqrt[2/3], -(2 Sqrt[2])/3, -7/3, 0, 0,
> =A0 2 (1/(3 Sqrt[2]) + (2 Sqrt[2])/3), Sqrt[10/3]}, {0, 0, 0, 0, 0,
> =A0 2 Sqrt[2/3], 0, (2 Sqrt[2])/3, 0, 0, -16/3,
> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], 0}, {0, 0, 0, 0, 0,
> =A0 =A00, 2 Sqrt[2/3], 0, (2 Sqrt[2])/3, 0,
> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -8/3, -(2 Sqrt[2])/3, 0}, {0, 0,
> =A0 =A00, 0, 0, 0, 0, 0, 0, 2 (1/(3 Sqrt[2]) + (2 Sqrt[2])/3),
> =A0 2 Sqrt[2/3], -(2 Sqrt[2])/3, 1/2,
> =A0 2 (-Sqrt[5/3]/16 - Sqrt[15]/16)}, {0, 0, 0, 0, 0, 0, 0, 0, 0, Sqrt[
> =A0 10/3], 0, 0, 2 (-Sqrt[5/3]/16 - Sqrt[15]/16), 7/2}}

[David Bailey]

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.
>

I am guessing here, but if you display the matrix it contains terms like
Sqrt[3] which has two possible answers +1.73.. and -1.73.... Although
this will not matter if you apply N to the matrix first, and then take
the eigenvalues, if you start by taking eigenvalues of the matrix you
get a horrendous expression involving Root objects (try it) and perhaps
the above confusion has arisen in that process.

In any case, for both efficiency and accuracy reasons it is preferable
to use use N as early as possible - EigenValues[N[matrix]]

David Bailey
http://www.dbaileyconsultancy.co.uk

[Andrzej Kozlowski]

Eigenvalues or characteristic values of a are defined (or rather, can
be defined) as the roots of the characteristic polynomial - and it
does not matter is the matrix is invertible or not. Indeed, for a
nilpotent matrix, such as

M = {{-1, I}, {I, 1}}

we have

In[39]:= Eigenvalues[M]
Out[39]= {0, 0}

and

In[40]:= CharacteristicPolynomial[M, x]
Out[40]= x^2

Moreover, the problem has nothing to do with numerical precision
because in this case the exact eigenvalues do not satisfy the
characteristic polynomial and in fact are the exact roots of a
different polynomial of the same degree (as shown in my first post in
this thread). Very weird.

Andrzej Kozlowski

On 2 Mar 2008, at 19:56, David Reiss wrote:

> Note that your matrix is not invertible (its determinant is zero). So
> this is the source of your problem...
>
> Hope that this helps...
>
> -David
> A WorkLife FrameWork
> E x t e n d i n g MATHEMATICA's Reach...
> http://scientificarts.com/worklife/
>
>
>
>
> On Mar 1, 4:57=A0am, Sebastian Meznaric <mezna...@gmail.com> 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 -

>> =A0 0.88758 \[ImaginaryI], -7.37965 + 2.32729 \[ImaginaryI],
>> -7.37965 -
>> =A0 2.32729 \[ImaginaryI], -4.46655 + 2.59738 \[ImaginaryI],
>> -4.46655 -

>> =A0 2.59738 \[ImaginaryI], 4.36971, 3.21081, -2.32456 +
>> =A0 2.10914 \[ImaginaryI], -2.32456 - 2.10914 \[ImaginaryI],
>> =A02.04366+ 0.552265 \[ImaginaryI],

>> =A02.04366- 0.552265 \[ImaginaryI], -0.249588 +
>> =A0 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,

>> =A0 0, -Sqrt[3], 0, 0, Sqrt[3], 0, 0, 0, 0, 0, 0, 0}, {-Sqrt[3], 0,
>> -4,
>> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], 0, 0, Sqrt[3],
>> 0,
>> =A0 0, 0, 0, 0, 0}, {0, -Sqrt[3],
>> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -4/3, -(2 Sqrt[2])/3, 0, 0,
>> 0,
>> =A0 Sqrt[3], 0, 0, 0, 0, 0}, {0, 0, 2 Sqrt[2/3], -(2 Sqrt[2])/3,
>> 7/3, 0,
>> =A0 =A00, 0, 0, Sqrt[3], 0, 0, 0, 0}, {Sqrt[3], 0, 0, 0, 0, -4, 0,
>> =A0 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, 0, 2 Sqrt[2/3], 0, 0, 0}, {0,
>> =A0 Sqrt[3], 0, 0, 0, 0, -4, 0, 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, 0,
>> =A0 2 Sqrt[2/3], 0, 0}, {0, 0, Sqrt[3], 0, 0,
>> =A0 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4), 0, -14/3,
>> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], (2 Sqrt[2])/3,
>> 0,
>> =A0 0, 0}, {0, 0, 0, Sqrt[3], 0, 0, 2 (-1/(4 Sqrt[3]) + Sqrt[3]/4),
>> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -2, -(2 Sqrt[2])/3, 0, (
>> =A0 2 Sqrt[2])/3, 0, 0}, {0, 0, 0, 0, Sqrt[3], 0, 0,
>> =A0 2 Sqrt[2/3], -(2 Sqrt[2])/3, -7/3, 0, 0,
>> =A0 2 (1/(3 Sqrt[2]) + (2 Sqrt[2])/3), Sqrt[10/3]}, {0, 0, 0, 0, 0,
>> =A0 2 Sqrt[2/3], 0, (2 Sqrt[2])/3, 0, 0, -16/3,
>> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), 2 Sqrt[2/3], 0}, {0, 0, 0,
>> 0, 0,
>> =A0 =A00, 2 Sqrt[2/3], 0, (2 Sqrt[2])/3, 0,
>> =A0 2 (-1/(4 Sqrt[3]) + (3 Sqrt[3])/4), -8/3, -(2 Sqrt[2])/3, 0},
>> {0, 0,
>> =A0 =A00, 0, 0, 0, 0, 0, 0, 2 (1/(3 Sqrt[2]) + (2 Sqrt[2])/3),
>> =A0 2 Sqrt[2/3], -(2 Sqrt[2])/3, 1/2,
>> =A0 2 (-Sqrt[5/3]/16 - Sqrt[15]/16)}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
>> Sqrt[
>> =A0 10/3], 0, 0, 2 (-Sqrt[5/3]/16 - Sqrt[15]/16), 7/2}}
>
>

[Daniel Lichtblau]

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

[David Reiss]

Thanks.... I stand (or sit) corrected.

I note that although N[Eigenvalues[matrix]] does not give correct
results (with the bug that Daniel Lichtblau explains),

Eigenvalues[N[matrix]]

does indeed work as expected...

--David

On Mar 3, 4:48 am, Andrzej Kozlowski <a...@mimuw.edu.pl> wrote:
> Eigenvalues or characteristic values of a are defined (or rather, can =

> be defined) as the roots of the characteristic polynomial - and it
> does not matter is the matrix is invertible or not. Indeed, for a
> nilpotent matrix, such as
>
> M = {{-1, I}, {I, 1}}
>
> we have
>
> In[39]:= Eigenvalues[M]
> Out[39]= {0, 0}
>
> and
>
> In[40]:= CharacteristicPolynomial[M, x]
> Out[40]= x^2
>
> Moreover, the problem has nothing to do with numerical precision
> because in this case the exact eigenvalues do not satisfy the
> characteristic polynomial and in fact are the exact roots of a
> different polynomial of the same degree (as shown in my first post in
> this thread). Very weird.
>
> Andrzej Kozlowski
>
> On 2 Mar 2008, at 19:56, David Reiss wrote:
>

> > Note that your matrix is not invertible (its determinant is zero). So=

> > this is the source of your problem...
>
> > Hope that this helps...
>
> > -David
> > A WorkLife FrameWork
> > E x t e n d i n g MATHEMATICA's Reach...
> >http://scientificarts.com/worklife/
>

> > On Mar 1, 4:57 am, Sebastian Meznaric <mezna...@gmail.com> 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, =

[Andrzej Kozlowski]

Perhaps I should have added that there are good reasons for having
both kinds of algorithms: symbolic and numeric. This example
illustrates a situation when there is a bug in a symbolic algorithm,
but numeric algorithms can also be problematic, particularly when
machine precision is used in the presence of numerical instability
(ill-conditioning). So, while in this case Eigenvalues[N[matrix]]
(numeric) is of course preferable to N[Eigenvalues[matrix]]
("essentially" symbolic) I would not recommend adopting this as a
general principle.

Andrzej Kozlowski

On 5 Mar 2008, at 11:23, Andrzej Kozlowski wrote:

> Yes, since it uses quite different algorithm. Generally algorithms
> used for numerical solution of equations or various matrix
> computations are quite different from (and much faster than) those
> used in the symbolic case.
>
> Andrzej Kozlowksi

>>>>> 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, =

[Andrzej Kozlowski]

Andrzej Kozlowksi

>>>> 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, =

[David Reiss]

At the expense of stating the obvious, I suspect that most readers of
this newsgroup understand this point...

-_David

On Mar 6, 3:09 am, Andrzej Kozlowski <a...@mimuw.edu.pl> wrote:
> Yes, since it uses quite different algorithm. Generally algorithms
> used for numerical solution of equations or various matrix
> computations are quite different from (and much faster than) those
> used in the symbolic case.
>
> Andrzej Kozlowksi
>
> On 5 Mar 2008, at 09:38, David Reiss wrote:
>
> > Thanks.... I stand (or sit) corrected.
>
> > I note that although N[Eigenvalues[matrix]] does not give correct
> > results (with the bug that Daniel Lichtblau explains),
>
> > Eigenvalues[N[matrix]]
>
> > does indeed work as expected...
>
> > --David
>
> > On Mar 3, 4:48 am, Andrzej Kozlowski <a...@mimuw.edu.pl> wrote:
> >> Eigenvalues or characteristic values of a are defined (or rather,

> >> can be defined) as the roots of the characteristic polynomial - and it

> >> does not matter is the matrix is invertible or not. Indeed, for a
> >> nilpotent matrix, such as
>
> >> M = {{-1, I}, {I, 1}}
>
> >> we have
>
> >> In[39]:= Eigenvalues[M]
> >> Out[39]= {0, 0}
>
> >> and
>
> >> In[40]:= CharacteristicPolynomial[M, x]
> >> Out[40]= x^2
>
> >> Moreover, the problem has nothing to do with numerical precision
> >> because in this case the exact eigenvalues do not satisfy the
> >> characteristic polynomial and in fact are the exact roots of a
> >> different polynomial of the same degree (as shown in my first post in
> >> this thread). Very weird.
>
> >> Andrzej Kozlowski
>
> >> On 2 Mar 2008, at 19:56, David Reiss wrote:
>
> >>> Note that your matrix is not invertible (its determinant is
> >>> zero). So

> >>> this is the source of your problem...
>
> >>> Hope that this helps...
>
> >>> -David
> >>> A WorkLife FrameWork
> >>> E x t e n d i n g MATHEMATICA's Reach...
> >>>http://scientificarts.com/worklife/
>
> >>> On Mar 1, 4:57 am, Sebastian Meznaric <mezna...@gmail.com> wrote:

[Andrzej Kozlowski]

Yeah, you are probably right. I guess if you teach the same boring
undergraduate classes for too many years, you yourself become a bore...
However, since there may still be *some* people on this forum to whom
these things are not entirely obvious, I feel compelled (sorry) to
observe that your original remark that the problem was caused by the
matrix not being invertible (e.g. determinant being zero) would have
been a reasonable guess had it been the case that
Eigenvalues[N[matrix]] did not work as expected ...

Andrzej Kozlowski