# Complex eigenvalues from a sparse Hermitian matrix

Bug introduced in 9.0 or earlier and persisting through 12.2.

I notice in the following example that wrong complex eigenvalues are resulted if calculating from a Hermitian sparse matrix, which should by no means have unreal eigenvalues. However, it gives correct result if we

• calculate from the corresponding normal matrix

I found many cases with this behavior. Actually, it becomes complex once n>24.

I construct the big matrix from 2by2 building blocks A and B. It's just a piece of concise code to make the big matrix of the form like $$\begin{bmatrix} B & A[1] & 0 & 0 \\ A[-1] & B & A[1] & 0 \\ 0 & A[-1] & B & A[1] \\ 0 & 0 & A[-1] & B \end{bmatrix}$$.

n = 25; Nless = 10;
BlockDiag[tt_, offset_] :=
DiagonalMatrix[Hold /@ tt, offset] // ReleaseHold // ArrayFlatten;
A[pm_] := {{-1, pm I}, {pm I, 1}};
B = (2.0 - Sqrt[2]) {{1, 0}, {0, -1}};
M = SparseArray[
BlockDiag[Table[B, {j, 1, n}], 0] +
BlockDiag[Table[A[1], n - 1], 2] +
BlockDiag[Table[A[-1], n - 1], -2]];
Reverse[Eigenvalues[M, -Nless]]


The result is

{-8.50551*10^-14, 8.50983*10^-14,
1.42089 + 0.0000210718 I, -1.42139 + 0.0000878787 I, -1.43983 -
0.0000374648 I, 1.44086 - 0.0000496205 I,
1.47277 + 0.000192942 I, -1.47298 - 0.0000263272 I, -1.516 -
0.000610613 I, 1.5161 - 0.0000292111 I}


The imaginary parts are not negligibly small. Knowing other strange behavior of sparse matrix, it looks not quite safe to calculate eigensystem of sparse matrix by default in Mathematica?

• You should report this to support@wolfram.com. – user21 Nov 19 '18 at 11:45
• The eigenvalues output by the example code are similar on Mathematica 9.0 and 10.1 (after changing the iterator format in Table to one compatible with the earlier versions), so this bug is definitely older than 11.3. The results are not precisely the same, but the spurious imaginary values persist. – eyorble Nov 20 '18 at 8:23
• Bug report submitted to Wolfram.com. – xiaohuamao Nov 20 '18 at 17:25
• Please use the bug header verbatim to make it easier to extract information automatically. mathematica.stackexchange.com/tags/bugs/info – Szabolcs Nov 20 '18 at 18:18

Very good observation. Indeed, this issue is really frustrating. To single out the issue: It seems that Arnoldi's method is to blame:

Max@Abs@Im@Eigenvalues[M, -Nless]
Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "Arnoldi"]
Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "Direct"]
Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "FEAST"]


0.000610613

0.000610613

0

0

IRRC, Arnoldi's method may have problems when eigenvalues cluster around 0. Sometimes, introducing a shift into Arnoldi's method can help. Usually one shifts by a positive real number in order to make the matrix positive-definite. However, this changes also the ordering of the eigenvalues if the matrix is Hermiatian but not indefinite (an issue that has been also observed here several times). In my desperation, I tried to shift by I, and in this case, the imaginary parts are much smaller:

Max@Abs@Im@Eigenvalues[M, -Nless, Method -> {"Arnoldi", "Shift" -> I}]
Eigenvalues[M, -Nless, Method -> {"Arnoldi", "Shift" -> I}]


1.11022*10^-15

{1.51607, -1.51607, -1.473, 1.473, -1.44084, 1.44084, -1.42095, 1.42095, -8.51975*10^-14 + 8.88178*10^-16 I, 8.49552*10^-14 - 1.11022*10^-15 I}

To my surprise, the ordering of the eigenvalues seems to be more or less consistent with the outputs of the other methods (M has pairs of eigenvalues of same magnitude but with opposite signs. Since each numerical method may induce small erros, this can effect the default ordering which is by magnitude.)

No guarantees for the correctness of the results, though. I think a bug report is a good idea anyways.

# Edit

While I first wondered why shifting by I works, it just came to my mind that the function $$x \mapsto |x+I|$$ is monotonically increasing on the positive real axis:

ParametricPlot[{Abs[x], Abs[x + I]}, {x, -4, 4},
AxesLabel -> {"Abs[x]", "Abs[x+I]"}]


So for a Hermitian matrix M, the corresponding eigenvalues of M and M + I are in consistent ordering (up to numerical errors). And of course, this shift guarantees that 0 is not an eigenvalue of M + I. So, now I am more confident to suggest this hack.

# Edit 2

Another curiosum: Also shifting by 0 or None seems to force the implementation of Arnoldi's method to branch to a more stable subroutine:

Max@Abs@Im@Eigenvalues[M, -Nless, Method -> {"Arnoldi", "Shift" -> 0}]
Max@Abs@Im@Eigenvalues[M, -Nless, Method -> {"Arnoldi", "Shift" -> None}]


9.76729*10^-6

4.84089*10^-17

I would not recommend to use this for larger matrices, though. Here is a significantly faster way to build larger versions of your matrix:

n = 250;
Nless = 10;
A[pm_] := N[{{-1, pm I}, {pm I, 1}}];
B = (2.0 - Sqrt[2]) {{1, 0}, {0, -1}};
M = Plus[
SparseArray[
{
Band[{1, 1}] -> Table[B, {j, 1, n}],
Band[{1, 3}] -> Table[A[1], n - 1],
Band[{3, 1}] -> Table[A[-1], n - 1]
},
{n 2, n 2}], 0.
];


Running Arnoldi's method with "Shift" -> 0 for this bigger matrix returns an error:

Eigenvalues[M, -Nless,
Method -> {"Arnoldi", "Shift" -> 0, MaxIterations -> 10000}]


But it still seems to produce plausible results with "Shift" -> None without any complaints.

• @ΑλέξανδροςΖεγγ I did not mean to bite you away. Your example with SetPrecision is great because it shows that there is indeed a precision issue with Arnoldi's method. – Henrik Schumacher Nov 18 '18 at 12:29
• Never mind, it is more important to figure out what really is wrong :). – Αλέξανδρος Ζεγγ Nov 18 '18 at 14:24
• Actually I wanted to post this question: why sometimes eigenvalue does not bring the eigenvalue symmetrically around zero. So I realize now, it can also be related to the shift. – Rasoul-Ghadimi Jul 13 '20 at 2:52