Occasional `LAPACKException(22)` in `eigsolve`

[yakovbraver]

Consider diagonalisation of a real symmetric matrix in the following MWE:

import BandedMatrices as bm
import KrylovKit as kk
import LinearAlgebra as la

function test_krylov(a, b, N, n_j)
    h = bm.BandedMatrix(bm.Zeros(2n_j + 1, 2n_j + 1), (2N, 2N))
    h[bm.band(0)] .= [(2j/N)^2 + (a + b)/2 for j = -n_j:n_j]
    h[bm.band(-2N)] .= h[bm.band(2N)] .= a / 4
    h[bm.band(-N)] .= h[bm.band(N)] .= b / 4
    kk.eigsolve(h, 14*2N, :SR; krylovdim=120) # sometimes triggers `ERROR: LinearAlgebra.LAPACKException(22)`
    # H = Array(h); la.eigvals(H) # always succeds
end

test_krylov(-2000, 3, 4, 112)

Approximately once in 20 runs, eigsolve tirggers LinearAlgebra.LAPACKException(22):

ERROR: LinearAlgebra.LAPACKException(22)
Stacktrace:
[1] chklapackerror(ret::Int64)
   @ LinearAlgebra.LAPACK \julia\stdlib\v1.7\LinearAlgebra\src\lapack.jl:43
 [2] stegr!(dv::Vector{Float64}, ev::Vector{Float64}, Z::SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
   @ KrylovKit \packages\KrylovKit\kWdb6\src\dense\linalg.jl:482
 [3] tridiageigh!
   @ \packages\KrylovKit\kWdb6\src\dense\linalg.jl:115 [inlined]
 [4] eigsolve(A::BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, x₀::Vector{Float64}, howmany::Int64, which::Symbol, alg::KrylovKit.Lanczos{KrylovKit.ModifiedGramSchmidt2, Float64})
   @ KrylovKit \packages\KrylovKit\kWdb6\src\eigsolve\lanczos.jl:49
 [5] #eigsolve#38
   @ \packages\KrylovKit\kWdb6\src\eigsolve\eigsolve.jl:202 [inlined]
 [6] #eigsolve#36
   @ \packages\KrylovKit\kWdb6\src\eigsolve\eigsolve.jl:176 [inlined]
 [7] test_krylov(a::Int64, b::Int64, N::Int64, n_j::Int64)
   @ Main \KrylovTest.jl:12
 [8] top-level scope
   @ \KrylovTest.jl:16

When eigsolve does succeed, it always converges to the same result, which coincides with the one returned by LinearAlgebra.eigvals. Perhaps the problem is caused by the groups (of size 2N) of very closely spaced eigenvalues of the constructed matrix. The problem shows up as well for e.g., krylovdim=40 (with howmany=30), provided we keep the matrix of the same size as in the example above (where n_j=112). Tested on Julia 1.7.2 and KrylovKit 0.5.4.

[Jutho]

Thanks for this detailed and reproducible bug report. Nonetheless, this is a difficult one to debug, because of the LAPACK error, which is outside of my own control. Even digging into the LAPACK documentation did not make me much wiser as to what this error code exactly indicates. Nonetheless, I do indeed assume that it has to do with the clustered structure of the eigenvalues as you point out.

Perhaps an a not strongly related but nonetheless important note is that, theoretically speaking, Krylov methods cannot find degenerate eigenvalues. In principle, the starting vector (which in your case is chosen as a random vector by KrylovKit) has a particular component in each eigenspace, irrespective of its dimensionality, and Krylov methods than find ways to extract this individual components, which are thus eigenvectors of the matrix. However, this would be one eigenvector in every eigenspace, irrespective of that eigenspace being higher-dimensional. However, it's only because numerical precision introduces small errors in the orthonormalization process, that new, independent directions in higher-dimensional eigenspaces can be generated. But it is somewhat fragile to rely on this behaviour.

But I assume that in your case these eigenvalues are not exactly degenerate? Nonetheless, this behaviour is indicative of the kind of issues that you might expect to encounter.

I will do some further investigation, but I am afraid that I will not have an easy fix.

[yakovbraver]

Thank you for the detailed answer. The matrix in the above example does possess degenerate eigenvalues. We can make the MWE even more minimal by considering a matrix of only three non-zero bands:

function test_krylov2(;b, N, n_j)
    h = bm.BandedMatrix(bm.Zeros(2n_j + 1, 2n_j + 1), (N, N))
    h[bm.band(0)] .= [j^2 for j = -n_j:n_j]
    h[bm.band(-N)] .= h[bm.band(N)] .= b
    kk.eigsolve(h, 80, :SR; krylovdim=90)
end

test_krylov2(b=1000, N=3, n_j=100)

For N > 2, the constructed matrix will possess pairs of identical eigenvalues, triggering the exception. Note that the larger the value of b, the higher the likelihood of the exception. I think it would be nice if eigsolve could detect degeneracies and inform the user, returning the best-so-far eigenvectors.

Additionally, I have noticed that if one constructs a matrix of integers (e.g. using bm.Zeros{Int} in the above code), then the exception is never thrown and the degenerate eigenvalues are always calculated correctly. However, the eigenvalues that are not degenerate vary considerably between runs, and the relative error sometimes reaches 0.1 unit.

Anyway, thanks for the package!

[gaspardbb]

+1, stumbled on this issue too and took me a while to figure it out has it is randomly reproducible. I also had some degenerate eigenvalues. However, in my case I am surprised that the degenerate eigenvalues at the end of the spectrum makes estimating the single, top eigenvalue (with :LR) difficult. In any case, I can provide some data if that helps!

[draftman9]

@yakovbraver @Jutho @gaspardbb

I also encountered similar error in ExponentialUtilities.jl. May these info can help. See https://github.com/SciML/ExponentialUtilities.jl/issues/171