Is there a bug in Eigensystem[]?

Question

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.

Answer 1

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