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

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

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


It looks 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

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... Commented Sep 1, 2018 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]; Commented Sep 1, 2018 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? Commented Sep 1, 2018 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. Commented Sep 1, 2018 at 11:43
• With Julia I don't get large errors (with MATLAB I do). Commented Sep 1, 2018 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.

• The matrix is defective, so Eigenvalues[] will surely have trouble. Commented May 11, 2020 at 1:39

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