remove automatic shift-invert when which==:SM

[tomhaber]

PR against v0.4.0

Can probably do something with the explicittransform parameter.

resolves #87

[ViralBShah]

@stevengj @antoine-levitt @andreasnoack Can you review?

[ViralBShah]

@tomhaber I messed up the merge conflict. Is it possible for you to revert my commit and resolve the merge conflict? It really should be rebased to master. Sorry about this.

[ViralBShah]

What should we do with explicittransform?

[antoine-levitt]

Thanks, will take a look. Also pinging @jarlebring

[tomhaber]

I don't really understand the meaning of the explicittransform. Especially since sigma does a similar thing.
If you want shift-invert, why not have Arpack deal with that?

[ViralBShah]

I'm not sure about explicittransform myself. Since we are removing all the automatic stuff, it seems like we should remove this as well. Perhaps others more familiar should chime in.

[tomhaber]

+1 on removing it. Seems to make the interface even more complicated:

• largest/smallest after transform or before?
• what happens combined with sigma?

[antoine-levitt]

This PR looks good to me. It doesn't solve all the confusion, but it's a start. explicittransform was introduced by @jarlebring to deal with quirks in generalized eigenproblems if I remember well, let's see what he says.

[tomhaber]

I've fixed the ordering of the eigenvalues/vectors. In the SM case this was broken

Not sure if you like the which[1] == 'L' bit, but I can fold that into the conditions as well

Might need some additional testing for the SM case

[ViralBShah]

Should we get this merged for now?

[ViralBShah]

@tomhaber Can I give you admin on this repo so that you don't have to wait for one of us?

[jarlebring]

The reason for explicittransform is documented here JuliaLang/julia#24668 and #4. In simple terms we needed more detailed control of what mode was running. I'm totally in favor of reducing the interface, which would speak in favor of removing this kwarg. On the other hand, a major update with only plain wrappers of the fortran code would seem more natural than removing/fixing features one by one.

[ViralBShah]

Should we merge this as is now, and then work with explicittransform and other wrappers in a separate PR?

The main challenge here is that this package does not have a maintainer at the moment - and it would be great if one of the folks who uses it routinely can own the maintenance.

I try to keep things moving the best I can for now (and so do some of the others).

[antoine-levitt]

Since this is technically breaking maybe it would make more sense to break everything just once? I lack the motivation to do it because I use KrylovKit now...

[ViralBShah]

Yes, breaking everything once is a good idea. Could be in this PR, or a series of PRs and then tag a new minor release. Perhaps could even be 1.0. At that point, the package's capabilities will also become minimal and stable.

[jarlebring]

The main challange with 1.0 clean rewrite that forwards operations to ARPACK as I see it is that the ARPACK fortran code does not have a clean interface. Forwarding the "reverse communication" directly to julia user would be easy (and it is actually essentially available in the functions dsaupd_, dseupd_, dnaupd_, dneupd_) but too complicated for users.

More precisely: The Arpack fortran code does not own the main iteration loop, but only does a CPU-critical part of one step of the loop. The user is expected to write the for loop, which is in the julia code eupd_wrapper, and that's one point where it gets complicated in terms of special cases (symmetry, real vs complex, target, standard EV vs generalized EV).

A much cleaner interface would be easier to construct if we restrict the cases, which would cripple the functionality, e.g., if we would remove generalized EV mode and target can only be :LM or :LR. This would "cover" GEPs and other targets since the user can do shift and invert by hand. However it would be slightly slower / less memory efficient than using the fortran code.

[antoine-levitt]

I think what we'd want from the julia pov is something like KrylovKit's geneigsolve (https://jutho.github.io/KrylovKit.jl/latest/man/eig/#Generalized-eigenvalue-problems) : eigs(opA, opB, howmany, which) where which is the two-letter things that (can) correspond to extremal points : LM, LR, SR, LI, SI, and where opA and opB can be functions. Shift-invert has to be done by the user explicitly. Do you agree?

[jarlebring]

Okay. That's a different idea, but that may be better. I see your suggestion as removing the :SM-mode, and special case code. (I thought we wanted to reveal a more raw interface to the julia user.)

For the record, your suggestion would also "cripple" Arpack.f-functionality, since the mode :SM does something different in Arpack.f. Calling Arpack on a "by-hand shift-and-invert" is not quite the same as currently using the :SM-mode with a shift. I think it effectively means a slightly different Ritz value extraction. As far as I remember, it is due to the fact that symmetry information for GEP is lost: If A and B are symmetric, (A-\sigma B)\B is in general not symmetric, but Arpack.f can still exploit it in :SM-mode.

For me this is okay, since I think it would make the interface cleaner and allow us to remove a lot of code, which is very hard to maintain, and the above case can be handled separately if deemed necessary.

You also suggest to have function handles instead of matrix-objects. I see it as something separate. It would require real-ness, symmetry to be inserted as kwargs, which goes against the idea of reducing the interface. I usually use LinearMaps.jl for instead of function handles. More juliaesque.

[antoine-levitt]

I think KrylovKit doesn't support SM since it is unlikely to ever converge? But if Arpack does support it then sure, no reason we should forbid it.

KrylovKit's interface seems simple, clean and complete to me. I don't actually understand Arpack's fortran interface (I looked at it a long time ago, was thoroughly confused and never had the courage to come back to it...), is what I suggest different from what the Fortran Arpack does? Do you think we should have an interface closer to Fortran Arpack?

You also suggest to have function handles instead of matrix-objects. I see it as something separate. It would require real-ness, symmetry to be inserted as kwargs, which goes against the idea of reducing the interface. I usually use LinearMaps.jl for instead of function handles. More juliaesque.

OK I agree: essentially we should not type A and B and just allow the user to pass whatever object they want that supports A*x and a few methods like that (so it can be used with Base matrices plus LinearMaps). KrylovKit taking functions is very useful but we don't need it, I shouldn't have mentioned it.

[jarlebring]

KrylovKit's interface seems simple, clean and complete to me

It maybe complete, but considering the information handed in, it cannot exploit the situation I explained above. The symmetry information about for GEP is lost if you do shift-and-invert by hand, so I think Arpack.f can be faster for symmetric GEPs.

I think KrylovKit doesn't support SM since it is unlikely to ever converge? But if Arpack does support it then sure, no reason we should forbid it.

It's more complicated than that. Quite unintuative at first sight, the target parameter does not only influence how the Ritz value extraction method used. If :SM would only take the smallest Ritz values in the restart it would indeed rarely work. I think, when running in :SM mode in Arpack.f, you are expected to also/instead provide shift and invert operations. If we want support for :SM, we would need to handle all the cases that we currently handle and I don't see how the rewrite would simplify anything.

[antoine-levitt]

It maybe complete, but considering the information handed in, it cannot exploit the situation I explained above. The symmetry information about for GEP is lost if you do shift-and-invert by hand, so I think Arpack.f can be faster for symmetric GEPs.

I think I don't understand the kind of problems you have in mind. If you can cholesky B, then you can just convert your problem to standard form while keeping the symmetry and so you don't need anything fancy from Arpack, just do shift-invert by hand. If you can't cholesky B, how do you do shift-invert while keeping the symmetry?

[jarlebring]

If you can cholesky B, then you can just convert your problem to standard form while keeping the symmetry and so you don't need anything fancy from Arpack, just do shift-invert by hand. If you can't cholesky B, how do you do shift-invert while keeping the symmetry?

I agree that a factorization can be computed (by hand) to carry out matrix vector products. There is more to Arpack (and IRAM in general) than the matrix vector product operations. I already mentioned one point. I think the Ritz value extraction depends on the original matrices.

[antoine-levitt]

OK, I see, sorry I'm slow. Is there reason to think the Ritz values of the original pencil are better than those of the transformed one? I did play with this kind of things at one point but couldn't see any difference. Otherwise I'd be tempted to say we just ignore this (messy) part of Arpack and just expose the "standard" mode.

[jarlebring]

I don't know a specific situation / paper where this setup is a huge advantage, but I think it is along the same lines as the difference between Ritz values and harmonic Ritz values (which are not supported in Arpack). Removing the messy modes and do the standard modes properly 👍. Once we figure out how to simplify the standard mode, it might be clearer how we want to implement the remaining parts, but possibly as separate functions (e.g. eigs_M_innerproduct referring to the case in §3.2.1 in
the Arpack manual)

[antoine-levitt]

That seems like a good plan! So to sum up: simply expose eigs(A, B, howmany, which) and forget about shift-invert completely.

[ViralBShah]

Won't that break people's codes? Maybe we ask them to clamp on an older version then?

What's the plan on this PR?

[antoine-levitt]

I think the plan is to break everything and request people to do shift invert manually with a much simple code and API. It should be OK if we tag a new major version, right? If we do that we should do it only once, so probably not merge this PR as is but have a bigger one where we do all of it. Or we merge this into a separate branch and work from there?

[jarlebring]

This level of the arpack call is very tricky. @tomhaber why you prefer your version of dmap combined with rev? Over just dmap including a sign in previous version?

[tomhaber]

For me, it simplifies the logic.
If I remember correctly, this version allows shift-invert to also work with non magnitude modes and fixes wrong ordering for SM and no shift-invert.

[jarlebring]

I understand the removal of if-block in Arpack.jl. That part looks fine. The change of ordering logic in eupd_wrapper in libarpack.jl seems better than current, because rev is more natural than sorting by adding minus signs, but there is a substantial difference if we consider all variations of which and iparam[7]. In contrast to previous version of the code, the composition with spectral transform is always applied if iparam[7]==3:

    if iparam[7] == 3 # shift-and-invert
       dmap = dmap ∘ (x -> 1 / (x - sigma))

Can you provide more details on the reasoning on this?

[ViralBShah]

How do folks feel about removing the automatic shift-invert for now? I suppose if we go down this path, we at least want other PRs ready for the full API simplification before we merge this one - otherwise we could end up with a partial API update and it will be difficult to make other bugfix releases.

[tomhaber]

I was working on some unit tests for the ordering logic (as asked by @jarlebring). It was pretty much done, but it totally got away from me.

[tomhaber]

I've added the unit test for the ordering modes.
The problems with the ordering are:

• SM and no shift-invert gives the wrong order of eigenvalues (I'm expecting the first one to be the smallest in magnitude)
• SM and shift-invert gives the wrong eigenvalues (I'm expecting the same as LM and no shift-invert)
• ARPACK seems to sort by the magnitude of the imaginary part for LI and SI: therefore I'm also using abs(imag(x))

Excerpt from the ARPACK sorting code

c          'LI' -> sort XIMAG into increasing order of magnitude.
c          'SI' -> sort XIMAG into decreasing order of magnitude.

I hope this makes sense and is correct

[ViralBShah]

Wondering what the next steps are here to get it merged.

[tomhaber]

I did a quick check on registered packages depending on Arpack and there are some that use explicit_transform and some that might be affected by a change in ordering from which=SM

explicit_transform

• FinEtoolsDeforLinear
• FinEtoolsFlexBeams

which=SM:

• FinEtoolsAcoustics
• Superfluids
• ExactDiagonalization
• FinEtoolsDeforLinear
• LightGraphs
• UMAP # My reason for wanting to fix this, so with this fix it's actually working correctly
• FrankWolfe
• ShellModel
• Mozi
• Laplacians

I also compared the results for the SM and no-shift-invert case with KrylovKit and it agrees with my test and fix

[ViralBShah]

That's great! We should bump the version number in Project.toml to 0.6.0 in this PR here and also make a discourse announcement.

Perhaps we should also update the project docs to mention these changes and how people should be able to reinstate some of this functionality on their own using this package.

[ViralBShah]

Shall we go ahead and merge?

[tomhaber]

I have nothing to add

[ViralBShah]

@tomhaber @jarlebring I have invited you both to this org and will make sure you have write permissions here to make it easier to contribute.

[ViralBShah]

Thank you for your contributions!

[matbesancon]

@ViralBShah should there be a new release after this?

[ViralBShah]

Yes, and it is a breaking change. I was trying to fix a few other things but was unable to fix. It would be giod to make a new release and we should announce it as well.

This could potentially be breaking for some packages.

We should discuss with @KristofferC how to do this in a sane way to minimize breakage.

[PetrKryslUCSD]

Related: #140

[PetrKryslUCSD]

I think :SM was broken somewhere along the way here: please see #147

[CodeLenz]

Hi.

I want to share a recent experience. Last week I was in an important Engineering conference. In the last day we had a "Computational and Numerical Methods in Practice – Scientific Dissemination" where I was showing some Julia code to the participants. To my surprise, my code to evaluate buckling response was broken and I did not have a clue about it (I had just shown the package manager and some packages were updated). Latter on I found this thread and it happens I use :SM and matrix B is not spd.

I really understand the burden of maintaining the code and the decisions you made, but such modification may have an important impact to the general user. Also, It would be advised to present some official workaround in the documentation, since not every user knows all the linear algebra theory needed to fix it (specially for large problems, where computing all the eigenvalues is not viable). Finally, we must remember that many users come from matlab/scilab/... where an easy interface for Arpack is provided (as was the Julia interface prior to this modification).

Thank you,
Eduardo.