JuliaLinearAlgebra
Arpack 0.5.4 does not converge for the smallest eigenvalue (:SM)
d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, tol=tol, maxiter=maxiter, explicittransform=:none, check = 1)
This call converges for 0.5.3 just fine. No eigenvalues converge for 0.5.4.
These are the relevant changes in libarpack.jl (0.5.3 on the left, master on the right):
The rest of the changes has to do with
Ptr -> Ref.
These are the changes in Arpack.jl (0.5.3 on the left, master on the right):
Reverting the first change in Arpack.jl fixed the failures for :SM eigenvalue problems in my tests.
@stevengj @antoine-levitt @andreasnoack @tomhaber I believe this is caused by the "shift" PR (https://github.com/JuliaLinearAlgebra/Arpack.jl/pull/120)? Could you review please?
So should we revert that PR then- it seems to have broken things for a few people.
The old arpack did very confusing things, and in particular inverted the matrix implicitly in a way that was hard to figure out if I remember correctly. It was decided in that PR to be explicit, so do something like F = factorize(A) ; eigs(x->F\x, :LM) to get back the old behavior. We should maybe document it more but it seems to me the right thing to do. The rationale is that eigs is somewhat of an advanced method for large scale systems, where you want to have control over what you factorize.
I think we should have a very good reason for breaking semantics that a call to Arpack eigs shares with Matlab. That currently does not work in 0.5.4.
I wanted the interface to be more explicit and mostly to "do as it was told". For my problematic case, I asked it to find the SM but it was actually computing the LM of the inverse which didn't work and gave unexpected wrong answers. Scipy has a similar explicit interface and did return the correct eigenvals.
I personally think whatever choice is made, it will not work for some specific problems and therefore might be better to be explicit and have the user decide how to solve the problem. An alternative would be to add an implicit/explicit options to disable automatic changing of the problem.
The current version (0.5.4) is in disagreement with its own documentation. The generalized eigenvalue problem for the smallest eigenvalues is supposed to work, and it doesn't.
The documentation should definitely clarify that restarted Arnoldi and similar methods are best for eigenvalues near the edges of the convex hull of the spectrum in the complex plane, and if you request eigenvalues in the interior of the spectrum (which can be the case for :SM) it may not converge (the standard solution being to replace the operator with the inverted operator, or more generally shift-and-invert). It can also converge slowly (or not at all) if you have a lot of other eigenvalues clustered near the desired eigenvalues, and again the solution is shift-and-invert.
@stevengj This is indubitably true. However, my problem was definitely one where the smallest eigenvalues were at one end of the spectrum. 0.5.3 converges just fine, 0.5.4 fails.
I am experiencing the same issue with some code breaking as a result of 0.5.4. While the issue can be fixed by explicitly rewriting the problem using LinearMaps and doing some post-processing, the easiest solution was to just downgrade to 0.5.3 which is not ideal.
It is not really my field of expertise, so I do not know what the best/recommended solution here is. Should one instead switch to geneigsolve from KrylovKit.jl?
Perhaps try https://github.com/PetrKryslUCSD/VibrationGEPHelpers.jl (Arpack is pinned at 0.5.3). This would allow you to evaluate the various packages for eigenvalue problems.