Beyn rankdrop checking

[jarlebring]

This aims to fix #162.

• Added (optional) sanity checking with errmeasure and tol
• Added neigs=wanted eigs
• Added user recommendations in case of insufficient eigvals found

Currently, the rankdrop is only used to decide on what the user recommendation should be in case insufficient eigvals are found.

[wrs28]

Hmm, still seems to not quite capture the right behavior. For example, the following problem: (A+1k)v= λv, where A = diag(1,2,3,4,...) should have eigenvalues (-1,-2,-3,-4,...).

But the following triggers a warning:

using NonlinearEigenproblems, LinearAlgebra, SparseArrays
N=1000
As = [sparse(1.0I,N,N), spdiagm(0=>1.0*(1:N))]
fns = [identity, one]
nep = SPMF_NEP(As,fns)
contour_beyn(nep;neigs=3, σ=-3, radius=.5)

It gives the correct answer (-3 eval) but it should not trigger the warning.
Also with the same inputs as above the call contour_beyn(nep;neigs=3, σ=-3.5, radius=.1) should not trigger a warning, it should just return an empty eval/evec, as there are no evals between -3.6 and -3.4. But, as written, it does trigger a warning, since all the singular values are small.

I think there needs to be both a measure of a drop in rank (as measured by S[:]/S[1]) as well as an absolute rank tolerance (S[1] > abs_rank_tol to keep).

I like the neigs and k keywords, and the evec computation!

[wrs28]

I've added an absolute rank tolerance on my branch, and I see now how this logic works, so I slightly revise my earlier comment.

I think as written (having added an absolute tolerance), if neigs≤p it returns all p without triggering a warning, which seems less than ideal.
If neigs>p and p==k a warning is triggered.
If neigs>p and p!==k, which is equivalent to p<k, it also triggers a warning (insufficient number of eigenvalues).

So, once an absolution rank tolerance is implemented (which handles the case of no enclosed eigenvalues), I think the only scenario that really needs addressing is the case neigs<p, as there is not an unambiguous way to return neigs evals/evecs.

I guess we should think about how this should operate. Should it always throw a warning if neigs≂̸p? I think that is reasonable, as long as the warning is descriptive. How many evals/evecs to return if neigs≂̸p? The case neigs≤p seems easy, but the case neigs<p seems like it should either return all p with a warning or perhaps throw an error. However, in this last case, if it returns all p, one must question what the real meaning of the neigs keyword is, so perhaps an error is the most appropriate.

[wrs28]

Ah, I can think of an easy solution. If we treat it like eigs, it returns the neigs evals/evecs that are closest to σ. Therefore so long as neigs≤p<k no warnings are thrown, and it returns the first neigs evals/evecs sorted by their distance from σ. I think the rest of the logic is fine.

How does that sound?

[jarlebring]

Thanks. Indeed there is still some improvements possible. I have changed the logic (and added sorting + only taking max neigs eigpairs) in the commit 182ecb5

With the nep in your example:

1. (λ,v)=contour_beyn(nep;neigs=3, σ=-3.1, radius=0.5); -> p=1. One eigval returned + warning: "We found less eigvals than requested."
2. (λ,v)=contour_beyn(nep;neigs=3, σ=-3.1, radius=1.5); -> p=3. Three eigvals returned. No warning.
3. (λ,v)=contour_beyn(nep;neigs=2, σ=-3.5, radius=0.6,k=2) -> p=2. Two eigvals returned. Warning: "No rankdrop detected"
4. (λ,v)=contour_beyn(nep;neigs=2, σ=-3.1, radius=1.5); -> p=3. Two eigvals returned. Warning: "No rankdrop detected"
5. (λ,v)=contour_beyn(nep;neigs=1, σ=-3.4, radius=0.8, k=3) -> p=2. One eigval (λ=-3) returned. No warning.
6. (λ,v)=contour_beyn(nep;neigs=typemax(Int), σ=-3.4, radius=0.8,k=3) -> p=2. Two eigvals returned. No warning

It should be consistent with what you write. Did I forget something?

Point 1 corresponds to your first example where you wanted no warning. I think a warning that insufficient eigvals were computed is okay. The warning can be suppressed with sanity_check=false.

The sorting would have to be replaced by some form of inside contour checking once we add support for domains other than disks.

[wrs28]

I think this is great, thanks!

[jarlebring]

👍