Is there a bug in Eigensystem[]?

Does Eigensystem[] produce incorrect output for symmetric matrices with integer components?

The following eigensystem decomposition of a 12x12 matrix and its subsequent reconstruction depends on whether the matrix components are integer or real values:

dim=12;
m = {{12, 0, 0, -4, 0, 0, -4, 0, 0, -4, 0, 0},
{0, 12, 0, 0, -4, 0,  0, -4, 0, 0, -4, 0},
{0, 0, 12, 0, 0, -4, 0, 0, -4, 0, 0, -4},
{-4, 0, 0, 12, 0, 0, -4, 0, 0, -4, 0, 0},
{0, -4, 0, 0, 12, 0, 0, -4, 0, 0, -4, 0},
{0, 0, -4, 0, 0, 12, 0, 0, -4, 0, 0, -4},
{-4, 0, 0, -4, 0, 0, 12, 0, 0, -4, 0, 0},
{0, -4, 0, 0, -4, 0, 0, 12, 0, 0, -4, 0},
{0, 0, -4, 0, 0, -4, 0, 0, 12, 0, 0, -4},
{-4, 0, 0, -4, 0, 0, -4, 0, 0, 12, 0, 0},
{0, -4, 0, 0, -4, 0, 0, -4, 0, 0, 12, 0},
{0, 0, -4, 0, 0, -4, 0, 0, -4, 0, 0, 12}};
es = Eigensystem[m];
eveks = Table[Normalize[es[[2, i]]], {i, 1, dim}];
Print["Check eigenvector orthogonality"]
Table[eveks[[i]] . eveks[[j]], {i, 1, dim}, {j, 1, dim}] // Chop // MatrixForm
Print["Check reconstruction of matrix from eigensystem"]
(Sum[es[[1, i]] Outer[Times, eveks[[i]], eveks[[i]]], {i, 1, dim}] - m) // Chop


In MMA 13.2 the above produces the eigenvector metric

{{1, 0, 0, 1/2, 0, 0, 1/2, 0, 0, 0, 0, 0},
{0, 1, 0, 0, 1/2, 0, 0, 1/2, 0, 0, 0, 0},
{0, 0, 1, 0, 0, 1/2, 0, 0, 1/2, 0, 0, 0},
{1/2, 0, 0, 1, 0, 0, 1/2, 0, 0, 0, 0, 0},
{0, 1/2, 0, 0, 1, 0, 0, 1/2, 0, 0, 0, 0},
{0, 0, 1/2, 0, 0, 1, 0, 0, 1/2, 0, 0, 0},
{1/2, 0, 0, 1/2, 0, 0, 1, 0, 0, 0, 0, 0},
{0, 1/2, 0, 0, 1/2, 0, 0, 1, 0, 0, 0, 0},
{0, 0, 1/2, 0, 0, 1/2, 0, 0, 1, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}


and the difference between the matrix and its reconstruction in spectral form is

{{12, 0, 0, -4, 0, 0, -4, 0, 0, -4, 0, 0},
{0, 12, 0, 0, -4, 0, 0, -4, 0, 0, -4, 0},
{0, 0, 12, 0, 0, -4, 0, 0, -4, 0, 0, -4},
{-4, 0, 0, -4, 0, 0, 4, 0, 0, 4, 0, 0},
{0, -4, 0, 0, -4, 0, 0, 4, 0, 0, 4, 0},
{0, 0, -4, 0, 0, -4, 0, 0, 4, 0, 0, 4},
{-4, 0, 0, 4, 0, 0, -4, 0, 0, 4, 0, 0},
{0, -4, 0, 0, 4, 0, 0, -4, 0, 0, 4, 0},
{0, 0, -4, 0, 0, 4, 0, 0, -4, 0, 0, 4},
{-4, 0, 0, 4, 0, 0, 4, 0, 0, -4, 0, 0},
{0, -4, 0, 0, 4, 0, 0, 4, 0, 0, -4, 0},
{0, 0, -4, 0, 0, 4, 0, 0, 4, 0, 0, -4}}


Clearly, the eigenvectors are not orthogonal and the reconstruction fails. If I convert the matrix components to floats, e.g. by changing the Eigensystem argument to es = Eigensystem[N[m]], everything works as expected.

1 Answer

Notice that there are two eigenvalues, both degenerate:

In[23]:= {eval, evec} = Eigensystem[m];

In[24]:= eval
Out[24]= {16, 16, 16, 16, 16, 16, 16, 16, 16, 0, 0, 0}


The corresponding two eigenspaces will be orthogonal, but within an eigenspace any vector is an eigenvector, and of course they are not all orthogonal to each other.

There is no bug here.

Since the computation is exact, we can easily verify that all eigenvectors are correct:

In[25]:= Table[
m . evec[[i]] == eval[[i]] evec[[i]],
{i, Length[m]}
]

Out[25]= {True, True, True, True, True, True, True, True, True, True, True, True}


We can also check that they are independent of each other:

In[26]:= MatrixRank[evec]

Out[26]= 12

• The question still remains why MMA gives different solutions for Integers and Reals? Why does MMA orthogonalize the Eigenvectors for Reals? Mar 14 at 20:46
• @granular I suspect they emerge orthogonal from underlying Lapack/BLAS routines in the machine number case. Mar 15 at 17:23
• @granular Mathematica uses very different algorithms for the exact case than the machine precision one, so it's no surprise that the results are different. Mar 15 at 19:38