RFC: selection of shift-and-invert mode in eigs

[acroy]

This changes the selection of shift-and-invert mode in eigs as discussed in #4474:

• if sigma==0 and which="LM" (default) the eigenvalues with largest magnitude are calculated
• if sigma==0 and which="SM" inverse iteration is used to find eigenvalues with smallest magnitude
• if sigma!=0 shift-and-invert is used to find eigenvalues closest to sigma

The docs are updated and I have added a test for the case with which="SM".
(cc: @ViralBShah, @stevengj)

[ViralBShah]

How about naming this variable arpack_which_lm?

[ViralBShah]

This looks good. The build failure seems to be not relevant.

[acroy]

One thing that might be problematic is that eigs(A, sigma=0.) will find the largest eigenvalues (since by default which="LM") although one could argue that it should give the smallest ones? One solution could be to change the default of sigma to Inf (does one actually need which at all in that case?).

[stevengj]

@acroy, that's a good point. Explicitly specifying sigma=0 should give the smallest eigenvalues.

[stevengj]

But it would still be nice to have a which argument if we eventually support finding the smallest/largest-algebraic (sa/la) eigenvalues for Hermitian matrices.

[acroy]

@stevengj : I see. But for SA/LA we don't need sigma, right? Maybe we could use sigma=Nothing() (or some symbol) as default and keep which as an argument. The behavior would be as follows:

• if isa(sigma, Nothing) and which="LM" (default) the eigenvalues with largest magnitude are calculated
• if isa(sigma, Nothing) and which="SM" inverse iteration is used to find eigenvalues with smallest magnitude
• if !isa(sigma, Nothing) shift-and-invert is used to find eigenvalues closest to sigma. The case sigma==0 could be handled separately to avoid the shifting.

[acroy]

I have updated my commit according to the discussion above. I am not sure about using nothing here though?

[ViralBShah]

I am not excited about nothing too. How come the build is still failing?

[acroy]

I am not excited about nothing too.

@ViralBShah : How about using sigma=() as a default?

How come the build is still failing?

This is strange. Last time only the gcc-build was failing in the tests (like in #5839), but now clang isn't even building.

[acroy]

I just figured out that eigs already supports which="LA" and "SA" for symmetric matrices! It's just not documented (#4476 can probably be closed). The question is if we want to have the MATLAB behavior, where specifying sigma always implies which="LM" (and where which="SM" is the same as sigma=0)? In that case "LA" and "SA" would not use shift-and-invert.

[acroy]

New attempt: Now sigma=() by default and in this case which is simply forwarded to ARPACK, otherwise shift-and-invert (or for sigma=0 inverse iteration) is used together with which="LM" to find the EV closest in magnitude to sigma. In contrast to MATLAB, which="SM" doesn't imply inverse iteration anymore (more similar to SciPy).

[stevengj]

I find sigma=nothing much more intuitive than sigma=(). (Julians instinctively recoil from the type-instability of using different types for the sigma parameter, but the performance impact should be negligible here.)

[acroy]

@stevengj @ViralBShah : I could live with both sigma=nothing or sigma=(). IMO the more important question is if which="SM" should imply inverse iteration or not.

[pao]

To help guide the sigma bikeshed, @StefanKarpinski recently commented that he found the use of () as an "empty" or "nil" signifier in Julia code to be odd in a comment on one of Jeff's recent commits.

[StefanKarpinski]

And importantly, Jeff agreed. Using () to mean nothing is not Julian. In a few places, I've used the empty string as a default or to mean nothing in contexts where a string values is expected. I felt only mildly guilty about it.

[JeffBezanson]

That's right; () is something, not nothing. For example it is the size of a 0-d array.

[acroy]

Ok, I will change it in the next commit. I just had the impression that nothing is not widely (at all?) used to indicate an undefined (=empty?) keyword.

[StefanKarpinski]

It's better to use a value of the same type, or redesign the interface so that there's no need for a nothing value, but in a pinch, you can use nothing this way.

[StefanKarpinski]

Part of the lack of clarity here is what the type of sigma is. The function signature doesn't indicate. Is it a float, an integer?

[acroy]

It may be float or complex. I think the problem is that one wants sigma to do two things: indicate the use of shift and invert and provide the value of the shift. Maybe it would be better to have two functions eigs(A; args...) and eigs(A, sigma; args...)?

[jiahao]

sigma now carries the double meaning that when it is nothing, it specifies ordinary (forward) iteration mode, and when it is a BlasFloat (real or complex), it specifies the level shift to be used in inverse iteration mode. sigma==0 is a special case where no shift is used but inverse iteration is still specified, which is semantically distinct from sigma==nothing where no shift is used since a different algorithm which requires no shift is used.

[jiahao]

I read the ARPACK manual and I don't see where exactly shift-and-inverse mode implies which="LM", so I took that out.

[jiahao]

Hm, I meant to update this PR but didn't do this right. #6053