LAPACK exception thrown from `eigen(::SymTridiagonal)`

[kbarros]

I'm getting a LAPACK error in the eigen-decomposition of a symmetric tridiagonal matrix,

using LinearAlgebra

mat = SymTridiagonal([1.9081958235744878e-17, 2.4936649967166602e-17, 0.07337398524939442, -0.037808419232506454, -0.019817057861955256, 0.02058130717549595, 0.0954455291925313, -0.14462978163628976, -0.03603900926534216, 0.10574604437162191, 0.0032290029813105, -0.060091988785585734, 0.15705171399550666, -0.14949686924737385, -0.07781372420855004, 0.03553277340671564, -0.016676626898563265, 0.041446049523777104, 0.10027425296079334, -0.10113668506693083, -0.1489473358740452, 0.1254312751044567, 0.34898062637027333, -0.33073386964350204, 0.2590576535805837, -0.25350583802618526, 0.43822601856881893, -0.4307646088223485, 0.779796594155107, -0.8224956491946511, 0.9089301791662714, -0.8619090557448831, 0.2318453071607881, -0.23320912386320142, 1.3140530211788515, -1.3139256900796155, -0.08390903751033216, 0.08392205661902721, 0.22240033452483454, -0.22246103504490297, 0.029954811490418787, -0.027812830389636778, 0.5263757772376259, -0.5365156423835895, 0.39737424822653156, -0.37936748640222284, 0.9333689384817099, -0.9703849367869426, 1.4044213669711443, -1.3377550567558987, 0.13172978466921414, -0.13519350037871353, 0.4999843222780838, -0.5567430191155477, 1.581258524837276, -1.5341564150622489, 0.8586746358146186, -0.8739294540224967, -0.17581567809216825, 0.15338653603265703, 0.780052785514151, -0.7342698870670235, -0.3173535579981542, 1.3364975241689332, -2.007075402840063, 1.0579643498230522, 0.31139215377118973, -0.39386762936936015, 1.8984071979199968, -1.920544791889951, 0.49358795540233863, -0.4613352087464951, 2.455277289025621, -2.4875372089056578, 0.9570609479497718, -0.9570638796447994, -1.1779359208299722, 1.163074811976796, -1.9301179244054274, 1.9450298141233324, -1.1378646951262403, 0.6350150906475753, 0.08802515616268328, 0.44163280306776015, -0.32773792909920113, 0.32952024534152274, 1.4740728008471815, -1.5038699938845244, 0.6847048634935279, -0.6874625497241057, 1.8751566802154134, -1.8547286700155103, 1.2640098482403956, -1.2174383863708917], [4.0208830900619335, 2.571546145080502e-15, 3.6768867672571894, 2.1275700268385056, 2.6591525597707077, 1.5630051292577454, 3.169596676339297, 1.4476618215762946, 3.1868304175740785, 1.424991114261858, 3.4173704890401693, 1.081122203108007, 3.8388891688820155, 0.7188503200238101, 3.628959037595643, 1.2527365479664625, 3.655577511598652, 0.5795107784623137, 3.9154918603711404, 0.7305602983115349, 4.204984775263665, 0.4667273034707952, 4.086056779436833, 0.7605098745710153, 3.765995021766351, 0.684123305522116, 4.180477278972855, 0.22109649225044853, 4.1892495574511095, 0.5187290599147804, 4.267105890594963, 0.38956213309842774, 4.395771702058166, 0.012968327575103515, 4.370023869671557, 9.368149211179955e-5, 1.7478905204265782, 0.012088046413693718, 2.5400351007261417, 0.02355280591484052, 2.966532136162213, 0.1157714633509318, 3.4610631638791225, 0.21831668966320836, 3.6404645575284627, 0.2911241355881879, 4.444637765978322, 0.46127544435962, 3.9217404428871023, 0.4673958058890463, 3.6185829343499414, 0.6185854819553929, 4.063721954940776, 0.4953177806242663, 4.017151353166678, 0.4035551805802303, 4.318876887238373, 0.17260927157837835, 4.41629817034136, 0.251496875950196, 4.017922104317157, 0.47982184361329844, 3.1370471883575424, 2.4334652498110247, 1.9981994742982454, 0.8095271305424269, 4.40837268034087, 0.18738580231835025, 4.147189662985785, 0.02664936509207341, 2.521497611392047, 0.4752292595225756, 3.795983764313679, 0.004636199088478572, 2.8022242947838416, 0.0004261020805263317, 3.2721199558188827, 0.16661903081014354, 3.1196340292912215, 0.03715487568175004, 2.6328978561492336, 1.5366796979500048, 2.1884275437040994, 0.7509912208995596, 4.056762760327159, 0.3470101947911033, 4.267342895676449, 0.4303372808790385, 4.051346888184243, 0.14621115143865943, 4.2045273832895145, 0.2540956144073191, 4.064260328749705])
eigen(mat) # LAPACKException(22) from stegr!

There is no problem for the same matrix in a dense representation,

mat2 = collect(mat)
values, vectors = eigen(mat2) # Works fine
extrema(values) # ≈ ± 4.6
norm(vectors) # ≈ 9.7

The problem seems independent of the BLAS backend.

BLAS.get_config()
#=
LinearAlgebra.BLAS.LBTConfig
Libraries: 
└ [ILP64] libopenblas64_.dylib
=#

using AppleAccelerate
BLAS.get_config() # LinearAlgebra.BLAS.LBTConfig, Libraries: [ILP64] libopenblas64_.dylib
#=
LinearAlgebra.BLAS.LBTConfig
Libraries: 
├ [ILP64] libopenblas64_.dylib├ [ LP64] Accelerate└ [ILP64] Accelerate
=#

eigen(mat) # Same error, LAPACKException(22)

The error persists between Julia versions (I tested 1.10 and 1.12).

[stevengj]

I can reproduce (on both macos/ARM and linux/x86_64).

It's possible that this is a bug in LAPACK itself. Computing all of the eigenvalues and taking their differences by diff(eigvals(Symmetric(Matrix(mat)))), it seems that this matrix is quite special: lots of eigenvalues come in degenerate pairs (to machine precision), so maybe one of the algorithms is hitting a bad corner case:

However, the fact that eigen(Symmetric(Matrix(mat)) works for the dense matrix is telling. The first step of the dense eigensolver is to reduce it to tridiagonal form via dsytrd, but in this case that should be the identity operation, and that seems to be true here as well:

julia> hessenberg(Symmetric(Matrix(mat))).H ≈ mat
true

I suspect that what is going on is that the generic eigensolver is doing some kind of balancing preprocessing that is not happening with the SymTridiagonal eigensolver. Update: nope, it's weirder than that, see below.

[stevengj]

Note that eigvals(mat) (which calls dstev) succeeds, but not eigen(mat). eigen(mat) calls dstegr, which is supposed to be just a compatibility wrapper around dstemr.

As I mentioned, eigen(Symmetric(Matrix(mat)) succeeds. It calls dsyevr, which looks like it simply calls dsytrd followed by dstemr, with no other preprocessing or balancing (it does scaling if maximum(abs, matrix) is very tiny or very huge, but that's not the case here). This is with the default RobustRepresentations() algorithm.

• Note that hessenberg(Symmetric(Matrix(mat))) also calls dsytrd, as noted above, and that this is the identity operation for the SymTridiagonal matrix mat here.

With eigen(Symmetric(Matrix(mat)); alg=LinearAlgebra.QRIteration()), which also succeeds, it calls dsytrd followed by dorgtr and dsteqr.

With eigen(Symmetric(Matrix(mat)); alg=LinearAlgebra.DivideAndConquer()), which also succeeds, it calls dsyevd, which calls dsytrd followed by dstedc.

In particular, note that eigen(Symmetric(Matrix(mat)) calls dstemr, which is supposed to be what dstegr is calling as well these days. But I wonder if our dstegr is somehow using an older algorithm? Should we try calling dstemr directly?

It would also be nice to provide a choice of algorithms for the SymTriangular eigensolver, like we do for eigen on dense symmetric matrices.

[stevengj]

Aha, wait a minute, dsyevr has an extra step: if dstemr fails, it calls dstebz and dstein.

I'm guessing that this is the problem: our SymTridiagonal eigensolver simply calls dstemr (via dstegr), but doesn't call the fallback algorithm on errors.

[stevengj]

I've confirmed that implementing the fallback algorithm works for this example. I'll have a PR shortly.

[stevengj]

As an interim workaround, you could directly copy and use the definition of syevr_tri_eigen from #1493.

[kbarros]

Thanks a lot for tracking this down!