# Inconsistency in eigenvalues of matrices in a specific form (sparse & non-Hermitian)

Suppose one has a non-Hermitian sparse matrix defined as below

Clear[h]
h[a_, n_: 10] := SparseArray[{Band[{2, 1}, {n, n - 1}] -> {1, a},
Band[{1, 2}, {n - 1, n}] -> {1, 0}}, n]


looked harmless and simple, with nonvanishing entries only distributed on the two secondary diagonal positions, which can be seen by running codes

h[a] // MatrixForm


and one gets

Well then if one is interested in its eigenvalues, with a set to various values,

Eigenvalues[h[#]] & /@ {a, 1, 1.}


Lo and behold!

{{-1, -1, -1, -1, -1, 1, 1, 1, 1, 1},
{-1, -1, -1, -1, -1, 1, 1, 1, 1, 1},
{-1.00044 + 0. I, -1.00014 + 0.000415446 I, -1.00014 - 0.000415446 I,
1. + 7.33851*10^-9 I, 1. - 7.33851*10^-9 I, 1. + 1.28894*10^-8 I,
1. - 1.28894*10^-8 I, 1. + 0. I, -0.999647 + 0.000256825 I,
-0.999647 - 0.000256825 I}
}


The first row of eigenvalues say that the eigenvalues are free of a, but with a = 1.0 non-negligible imaginary parts appear, which would get larger with larger ns.

I don't think it a bug, because it is said that similar behavior exists also in other languages, e.g., Python & MATLAB. How to explain this and how to correct it?

• Wow, these are really huge errors. Hm. This matrix as two eigenvalues with high multiplicity. Many numerical algorithms can perform computations only stably if the eigenvalues are well separated... – Henrik Schumacher Sep 1 '18 at 10:16
• This is a precision problem. See, in the following sense, the "wrong" eigensystem is quite correct: A = h[1.]; {\[Lambda], U} = Eigensystem[A]; Max[Abs[Table[ A.U[[i]] - \[Lambda][[i]] U[[i]], {i, 1, Length[\[Lambda]]}]]]. Raising the working precision above \$MachinePrecision helps; have a look at A = h[SetPrecision[1, 100]]; {\[Lambda], U} = Chop@Eigensystem[A]; – Henrik Schumacher Sep 1 '18 at 10:22
• @HenrikSchumacher Thanks, SetPrecision can deal with it. But there is another strange phenomenon that my above defined h, will not suffer this precision problem if transposed. Why? – Αλέξανδρος Ζεγγ Sep 1 '18 at 11:15
• Phew. Honestly, I don't know. That depends on which algorithm is used for computing the matrix. Often, a QR-factorization is performed and that induces a breaking of transposition-symmetry. – Henrik Schumacher Sep 1 '18 at 11:43
• With Julia I don't get large errors (with MATLAB I do). – Szabolcs Sep 1 '18 at 11:55

It naively looks like a precision problem, but increasing the precision does not immediately help. For example, with OP's definitions, I find on my computer following results:

Eigenvalues[h[1]]
Eigenvalues[h[1.5]]
Eigenvalues[h[1.10]]
Eigenvalues[h[1.100]]


{-1, -1, -1, -1, -1, 1, 1, 1, 1, 1}

{-1.0000, 1.0000, -1.0000 + 3.2888*10^-10 I, -1.0000 - 3.2888*10^-10 I, 1.0000 + 1.6608*10^-9 I, 1.0000 - 1.6608*10^-9 I, 1.0000 + 3.6260*10^-10 I, 1.0000 - 3.6260*10^-10 I, -1.0000, -1.0000}

{-1.000000000, 1.000000000 + 1.941032805*10^-11 I, 1.000000000 - 1.941032805*10^-11 I, 1.000000000, -1.000000000 + 6.539486873*10^-20 I, -1.000000000 - 6.539486873*10^-20 I, -1.000000000, 1.000000000 + 8.442441251*10^-20 I, 1.000000000 - 8.442441251*10^-20 I, -1.000000000}

{-1.000000000000000000000000000000000000000000000000000000000110954162\ 5915143645283700322398250961229805, \ 1.00000000000000000000000000000000000000000000000000000000000000000000\ 00000000000000000000000000000000, \ -1.0000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000 + 2.976853631677875691628349618766633756396719842482978185187919770933\ 308169805544158257069331975894975*10^-58 I, \ -1.0000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000 - 2.976853631677875691628349618766633756396719842482978185187919770933\ 308169805544158257069331975894975*10^-58 I, 1.0000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000 + 3.558020626382386620531022727918765817478289860513062693986158315078\ 213370359938021098572273950853424*10^-58 I, 1.0000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000 - 3.558020626382386620531022727918765817478289860513062693986158315078\ 213370359938021098572273950853424*10^-58 I, 1.0000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000 + 3.328624877745430935851100967194752883689414382046741134022348191945\ 432724468803967697080267882337017*10^-58 I, 1.0000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000 - 3.328624877745430935851100967194752883689414382046741134022348191945\ 432724468803967697080267882337017*10^-58 I, \ -1.0000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000, \ -0.9999999999999999999999999999999999999999999999999999999998890458374\ 0848563547162996776017490387701952}

Of course, we can use Chop or N at the end to get rid of the tails, but that does not solve the problem of getting these results out of Eigenvalues command.

One dirty way to solve this is to introduce an auxiliary symbolic variable to the process to avoid possible inner cancellations, and then take that variable to zero. For the same examples above, if we now run

Eigenvalues[h[1] + ConstantArray[\[Epsilon], {10, 10}]] /. \[Epsilon] -> 0
Eigenvalues[h[1.5] + ConstantArray[\[Epsilon], {10, 10}]] /. \[Epsilon] -> 0
Eigenvalues[h[1.10] + ConstantArray[\[Epsilon], {10, 10}]] /. \[Epsilon] -> 0
Eigenvalues[h[1.100] + ConstantArray[\[Epsilon], {10, 10}]] /. \[Epsilon] -> 0


one gets

{-1, -1, -1, -1, -1, 1, 1, 1, 1, 1}

{-1.00, -1.00, -1.00, -1.00, -1.00, 1.00, 1.00, 1.00, 1.00, 1.00}

{-1., -1., -1., -1., -1., 1., 1., 1., 1., 1.}

{-1.000000000000000000, -1.000000000000000000, -1.000000000000000000, \ -1.000000000000000000, -1.000000000000000000, 1.000000000000000000, \ 1.000000000000000000, 1.000000000000000000, 1.000000000000000000, \ 1.000000000000000000}

which is way better than the previous ones.

I am sure this method would fail under certain scenarios, e.g. adding a symbolic number may prevent finding all roots of characteristic polynomial. Still, this may be a quick work-around.

Another workaround is to just construct the characteristic polynomial and find the zeros yourself:

Solve[CharacteristicPolynomial[h[a], x] == 0, x]
Solve[CharacteristicPolynomial[h[1.], x] == 0, x]


{{x -> -1}, {x -> -1}, {x -> -1}, {x -> -1}, {x -> -1}, {x -> 1}, {x -> 1}, {x -> 1}, {x -> 1}, {x -> 1}}

{{x -> -1.}, {x -> -1.}, {x -> -1.}, {x -> -1.}, {x -> -1.}, {x -> 1.}, {x -> 1.}, {x -> 1.}, {x -> 1.}, {x -> 1.}}