Accuracy issues with Arpack in Julia for eigenvalues of smallest magnitude

Question

Following the documentation of Julia's Arpack package (Cf. https://julialinearalgebra.github.io/Arpack.jl/stable/eigs/) I have computed some largest and smallest magnitude eigenvalues of sparse matrices encoded in the required CSC format and noticed major discrepancies concerning the smallest magnitude eigenvalues $\theta$ giving me residuals $\mid \mid Av - \theta v \mid \mid / \mid \mid v \mid \mid$ in the order of $10^7$ though all of those eigenvalues were claimed to converge.

I compared the results to a computation with Python's scipy.sparse.linalg package (Cf. https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html) which gives appropriate results, i.e. residuals in the order of $10^{-5}$.

The sparse matrices arise from Macaulay matrices in Groebner basis theory so are kind of complicated to explain but I was able to reproduce the odd behavior with randomly generated sparse matrices of the same size and similar entries.

This Julia code

using Arpack, LinearAlgebra, SparseArrays

rf(n) = rand(-2^15+1:2^15, n)
S = sprand(229375, 16384, 0.0005, rf)
S = convert(SparseMatrixCSC{Float64,Int64},S)
A = transpose(S)*S
l = eigs(A, nev=4, which=:LM)
s = eigs(A, nev=4, which=:SM)

println("Number of converged eigenvalues of largest magnitude: ", l[3])
println("Number of converged eigenvalues of smallest magnitude: ", s[3])

# Tests

println( "Epsilon according to documentation: ", eps(real(eltype(A)))/2 )

# Largest Eigenvalue

x = 1
theta = l[1][x]
v = l[2][:,x]

println( "Residual of largest eigenvalue: ", norm(A * v - v * theta) )
println( "Convergence criterion for largest eigenvalue: ", (residual <= epsilon * max(cbrt(epsilon^2),abs(theta))) )
println( "Upper bound in convergence criterion for largest eigenvalue: ", epsilon * max(cbrt(epsilon^2),abs(theta)))

# Smallest Eigenvalue

y = 1
iota = s[1][x]
w = s[2][:,x]

println( "Residual of smallest eigenvalue: ", norm(A * w - w * iota) )
println( "Convergence criterion test for smallest eigenvalue: ", (residual <= epsilon * max(cbrt(epsilon^2),abs(iota))) )
println( "Upper bound in convergence criterion for smallest eigenvalue: ", epsilon * max(cbrt(epsilon^2),abs(iota)))

gives the following output

Number of converged eigenvalues of largest magnitude: 4
Number of converged eigenvalues of smallest magnitude: 4
Epsilon according to documentation: 1.1102230246251565e-16
Residual of largest eigenvalue: 0.00018192639656053935
Convergence criterion for largest eigenvalue: false
Upper bound in convergence criterion for largest eigenvalue: 8.018644716606565e-6
Residual of smallest eigenvalue: 0.00011757473085532065
Convergence criterion test for smallest eigenvalue: false
Upper bound in convergence criterion for smallest eigenvalue: 2.056712217644155e-6

The eigenvalues seem fine but the convergence criterion described in the above mentioned documentation of Julia's Arpack package seems to be different from the actual implementation.

I am using Julia 1.5.3 with

Installed Arpack_jll ─── v3.5.0+3
Installed OpenBLAS_jll ─ v0.3.9+5
Installed Arpack ─────── v0.5.1

In my specific use case mentioned above the computed smallest magnitude eigenvalues are far of anything useful. Also tuning the parameter tol didn't help. Can someone explain this strange behavior or help me out with assessing the right convergence criterion and tuning the computation towards useful residuals?

Update

The following Julia code provides an example for the above mentioned BAD residuals.

using Arpack, LinearAlgebra, SparseArrays

rr = [1,2,3,5,1,2,3,4,5,6,1,2,3,4,5,7,2,3,4,6,7,8,1,2,3,5,6,7,2,4,5,6,7,8,3,4,5,6,7,4,6,8]
cc = [1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,7,8,8,8]
dd = [4.02018713e8,1.42895766e8,1.91997182e8,-5.94015895e8,1.42895766e8,9.88669461e8,9.1763724e7,-3.932278e6,-2.3465949e8,-1.198600777e9,1.91997182e8,9.1763724e7,3.156467582e9,-6.25317095e8,-1.03130053e9,-2.317164234e9,-3.932278e6,-6.25317095e8,4.565270498e9,-4.93246362e8,6.25317095e8,-1.691847139e9,-5.94015895e8,-2.3465949e8,-1.03130053e9,2.883067709e9,-7.60572644e8,-1.257751284e9,-1.198600777e9,-4.93246362e8,-7.60572644e8,4.2947594e9,7.60572644e8,-4.9717864e8,-2.317164234e9,6.25317095e8,-1.257751284e9,7.60572644e8,3.574915518e9,-1.691847139e9,-4.9717864e8,5.815904688e9]

A = sparse(rr,cc,dd)

l = Arpack.eigs(A, nev=6, which=:LM)
s = Arpack.eigs(A, nev=6, which=:SM)

println("Number of converged largest magnitude Ritz values: ", l[3])
println("Number of converged smallest magnitude Ritz values: ", s[3])

epsilon = eps(real(eltype(A)))/2 # LAPACK's epsilon
tol = epsilon # tol default

println("Machine epsilon and tolerance default: ", epsilon)

l_res = [norm( A * l[2][:,i] - l[1][i] * l[2][:,i]) / norm( l[2][:,i] ) for i in 1:size(l[1])[1] ]
println("Residues of largest magnitude Ritz values: ", l_res)

s_res = [norm( A * s[2][:,i] - s[1][i] * s[2][:,i]) / norm( s[2][:,i] ) for i in 1:size(s[1])[1] ]
println("Residues of smallest magnitude Ritz values: ", s_res)

l_con = [(l_res[i] <= tol * max(cbrt(epsilon^2),abs(l[1][i]))) for i in 1:size(l[1])[1]]
println("Convergence criteria of largest magnitude Ritz values: ", l_con)

s_con = [(s_res[i] <= tol * max(cbrt(epsilon^2),abs(s[1][i]))) for i in 1:size(s[1])[1]]
println("Convergence criteria of smallest magnitude Ritz values: ", s_con)

The output is:

Number of converged largest magnitude Ritz values: 6
Number of converged smallest magnitude Ritz values: 6
Machine epsilon and tolerance default: 1.1102230246251565e-16
Residues of largest magnitude Ritz values: [2.5183513490019574e-6, 2.0857933315848784e-6, 2.3053979403750557e-6, 1.0323827311807143e-6, 4.182955279236495e-7, 1.3601799608960461e-6]
Residues of smallest magnitude Ritz values: [1.4473263034070623e-6, 3.3449964080033436e9, 2.50443508742856e9, 4.42611459230314e9, 1.1814710823584437e9, 6.968580232683431e8]
Convergence criteria of largest magnitude Ritz values: Bool[0, 0, 0, 0, 0, 0]
Convergence criteria of smallest magnitude Ritz values: Bool[0, 0, 0, 0, 0, 0]

If you analyze the (claimed) smallest magnitude Ritz value you will notice that the corresponding Ritz vector is an estimate of the zero vector. Also all other smallest magnitude Ritz values claimed to have converged are bogus. One can use the arrays above specifying the indices and data to instantiate a sparse matrix with scipy and compute trustworthy smallest magnitude Ritz values.

Seems pretty buggy. Or is there another explanation?

Answer 1

OK, so with new version in edit: it's not your fault, it's https://github.com/JuliaLinearAlgebra/Arpack.jl/issues/87. You can either call Arpack manually yourself, or even better use a pure-julia solver (KrylovKit.jl and IterativeSolvers.jl are good choices)

Answer 2

Looks like everything is working relatively OK? Your matrix is of order 1e10, so residuals of 1e-4 are actually close to machine precision. The convergence criterion is indeed violated, but not by much; not sure what's going on there, are you sure Arpack really guarantees it or is it a best effort kind of thing? It surprises me that you get different results in python because both call into the same underlying library (although may do so a bit differently). Also eigs is more directly comparable to numpy's eigs; if you want eigsh, you should do eigs(Symmetric(A)) in julia.

Answer 3

eigsh( ... sigma=0 ) looks like a bug in python-scipy arpack too. Here's a comparison of arpack with numpy eigh (Lapack syevd, gold standard), run on your matrix / 1e6:

sigma=0 --
evals numpy eigh:      [-8.23e-13 370 515 2.82e3 3.64e3 5.21e3 5.9e3 7.22e3]
evals arpack:     [-146 -5.68e-14 96.1 844 3.43e3 5.06e3 5.18e3] ❓

eigcheck AV - VΛ numpy: max 2.5e-12 at λ 3.64e3  [1.1e-12 1e-12 2.4e-12 1.9e-12 2.5e-12 1.4e-12 1.4e-12 1.8e-12]
eigcheck AV - VΛ arpack: max 1.7e3 at λ -146  [1.7e3 7.1e-13 9.3e2 1.3e3 3.1e2 7.5e2 1.9e2]

--------------------------------------------------------------------------------
arpack sigma=-1 agrees pretty well with numpy --
evals numpy eigh: [-8.23e-13 370 515 2.82e3 3.64e3 5.21e3 5.9e3 7.22e3]
evals arpack:     [7.73e-14 370 515 2.82e3 3.64e3 5.21e3 5.9e3]  # neig = n - 1

eigcheck AV - VΛ numpy: max 2.5e-12 at λ 3.64e3  [1.1e-12 1e-12 2.4e-12 1.9e-12 2.5e-12 1.4e-12 1.4e-12 1.8e-12]
eigcheck AV - VΛ arpack: max 4.2e-10 at λ 2.82e3  [6.4e-13 4.6e-11 7.5e-11 4.2e-10 4.8e-11 4.3e-11 3.8e-10]

Can anyone try julia arpack with sigma=-1 ?

Notes on Arpack

• Use eigsh( ... v0=np.ones ), not the default guess eigenvector v0 = random. Better yet, get v0 from a previous run or a coarser grid.

• for eigenvalues near a given sigma, scipy arpack calls different inner solvers of Ax = b for A dense / sparse / LinearOperator, by default LU / sparse LU / GMRES. But the complex back-and-forth arpack $\leftrightarrow$ solver is hidden, completely murky. (One can make a LinOp with callbacks to track this, or to use e.g. LGMRES instead of GMRES.)

• Shift clouds of near-0 eigenvalues to the right half-plane if you can: sigma=-1 does shift-invert on A + I, better conditioned, then shifts back. In this case, sigma=-1 works. (Also shift $L L^T$: if $L$ is numerically singular, $L L^T$ can well have tiny negative eigenvalues, tough for solvers.)

I stumbled across this in pybinding code:

   if sigma == 0:
        # eigsh can cause problems when sigma is exactly zero
        sigma = np.finfo(model.hamiltonian.dtype).eps