eigs with singular matrix

[alainchau]

I get NaNs when I try to use eigs on a singular matrix with which=:SM. 🙁
Is this behavior expected? which=:SR works for my use-case but I thought this was strange..
(I'm on Julia v1.3.0 and Arpack v0.4.0.)

julia> eigs(zeros(3,3), nev=2, which=:LM)
([0.0, 0.0], [0.5941386191532766 0.48617004898294336; -0.7258278379742127 0.6673894993800032; 0.3466601373857654 0.5641187282737836], 2, 1, 3, [0.0, 0.0, 0.0])

julia> eigs(zeros(3,3), nev=2, which=:SM)
([NaN, NaN], [NaN NaN; NaN NaN; NaN NaN], 2, 1, 3, [0.0, 0.0, 0.0])

julia> eigs(zeros(3,3), nev=2, which=:SR)
([0.0, 0.0], [0.45524873566374424 0.6593503833676924; -0.3016328607164829 -0.5180634739681107; -0.837714871547923 0.5448553100530177], 2, 1, 3, [0.0, 0.0, 0.0])

[ViralBShah]

The behaviour is slightly changed now - but by and large we just send this all to ARPACK and return whatever it gives back.

[RalphAS]

It's not that simple.

For which = :SM, the Julia wrapper silently converts the problem into shift-and-invert with zero shift; i.e. it does linear solves. If the small eigenvalues are not tiny and are in the interior of the spectrum, this is reasonable (albeit paternalistic), since the true "SM" approach may converge slowly or not at all in this case. On the other hand, for small values at the edge of the spectrum or in an isolated, loose cluster, the Arpack library's true "SM" approach often works without the need for a factorization, but it is not currently available through this package.

If the small eigenvalues are tiny, the current bait-and-switch is a trap. You are lucky to get NaN since it is more likely to converge to an incorrect result (inverse eigenvalues of an unstable factorization). Of course, similar problems can affect any use of shift-and-invert, but here it happens with no clue in the documentation.

In general, shift-and-invert is powerful for interior eigenvalues, but users should be advised to watch out for near-singularity. The current handling of which = :SM makes it an important subcase.

[antoine-levitt]

That's very annoying. Arpack's methods are sufficiently complex that we shouldn't add another layer of complexity... I was confused by #4, I'm even more confused now. Can't we just remove all of this complexity and just call Arpack with the provided arguments?

[antoine-levitt]

Ping @jarlebring and @stevengj : shouldn't we just remove every fancy feature of this library and just call ARPACK directly without explicit support for shift-invert, leaving it to the user to do it?

[ViralBShah]

What do you mean by calling ARPACK directly? I mean you can just do that by directly using the aupd and eupd wrappers, and then doing the full reverse communication on your own.

[antoine-levitt]

What I mean is that we should do minimal preprocessing and just convert julia's matrices/functions to something arpack understands. This particular issue is because in the :SM case (and possibly in other cases) this package tries to be clever and does a transformation (shift-invert) although the user didn't ask for it: that's what I meant by "fancy" above. I've always been thoroughly confused by both ARPACK and ARPACK.jl's handling of shift-invert (eg I never know whether the which argument refers to the eigenvalues of A or of the shift-inverted version), to the point that I don't trust it and always do it manually (ie for finding the SM eigenvalues I'd do B = factorize(A); arpack(x -> B\x, :LM)).

Btw this came up at https://scicomp.stackexchange.com/questions/36767/accuracy-issues-with-arpack-in-julia-for-eigenvalues-of-smallest-magnitude/36782

[ViralBShah]

I see what you mean. Provide a straightforward path from what ARPACK.jl's API to ARPACK itself - without doing these fancy things. That seems like a good idea. That way the user can depend on the ARPACK documentation.

[jarlebring]

I agree that a minimalistic interface sounds appealing, in particular considering that there are other Krylov-solver packages for Julia nowadays. A word of caution: There is also some "fancy" stuff in ARPACK itself which also can generate issues like above. I think we will get more difficult to interpret errors, which are directly related to the ARPACK-code. In packages like MATLAB there was high-level code to lessen the impact of these issues. Not a pretty sight. It was wise of MATLAB to replace eigs with their own standalone implementation of Krylov-Schur. (But their way of doing it is still not great because of legacy limitations for eigs, e.g. visible in restrictions in target specification and warmstarting.)

[antoine-levitt]

Hm. Should we just soft-deprecate (ie put a big redirect first thing in the README) ARPACK in favor of KrylovKit, which is pure-julia and has a clearer API?

[jarlebring]

It would be okay with me. If so, I would suggest we provide links to both KrylovKit.jl and ArnoldiMethod.jl. They both have advantages / disadvantages.

[ViralBShah]

I'm happy to add folks who can help with this to this repo.

[stevengj]

I'm generally in favor of having the caller do shift-and-invert explicitly if that is what is desired.

You are lucky to get NaN since it is more likely to converge to an incorrect result (inverse eigenvalues of an unstable factorization).

My understanding is that there is nothing generally wrong with solving a nearly singular system arising in shift-and-invert eigensolvers.

Although solving nearly-singular systems typically raises fears of an inaccurate solution because of ill-conditioning, in the specific case of eigensolvers it isn't a problem (at least in theory): the "error" term from ill-conditioning is in the direction of the desired eigenvector(s), and since we don't care about overall scale factors in eigensolvers you still get an accurate final result. There is a homework problem about this in Trefethen & Bau.

Of course, the implementation may get confused if it tries to solve a system that is exactly singular (hence NaNs), but in this case ideally the implementation would be a bit smarter and realize that it has an "exact" eigenvector from the nullspace in that situation.

[antoine-levitt]

There is a homework problem about this in Trefethen & Bau.

I remember that one! I had Trefethen in class and asked him exactly this question, and he pointed me to the exercice.

I agree, ideally, the factorization API would be somewhat richer and allow the factorization to return more information (eg a vector in the nullspace), but that seems too much to ask of a factorization algorithm. Also how would you get all the eigenvectors in something like diag(0,0,eps,1)? Surfing on the edge of instability like that is something best left to advanced users and manual intervention.

[ViralBShah]

Yup definitely remember learning about the error being in the right direction from my advisor John Gilbert, in class. :-)

[RalphAS]

@stevengj I think your argument applies to a nearly singular eigenpair, but if the user also wants some other (small but not tiny) eigenvalues, the factorization will be corrupt so those results will be unreliable. How else do you explain the problem in the stackexchange link?