# Accuracy issues with Arpack in Julia for eigenvalues of smallest magnitude

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)
println("Number of converged eigenvalues of smallest magnitude: ", s)

# Tests

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

# Largest Eigenvalue

x = 1
theta = l[x]
v = l[:,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[x]
w = s[:,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)
println("Number of converged smallest magnitude Ritz values: ", s)

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[:,i] - l[i] * l[:,i]) / norm( l[:,i] ) for i in 1:size(l) ]
println("Residues of largest magnitude Ritz values: ", l_res)

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

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

s_con = [(s_res[i] <= tol * max(cbrt(epsilon^2),abs(s[i]))) for i in 1:size(s)]
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?

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.