Implement eigenvalue estimation for eigs

[cshaa]

The current eigenvalue algorithm supports shifting however as of now it is not used. The gist of shifting is this:

1. The smaller (the absolute value of) the eigenvalue, the faster the convergence.
2. Let A be a matrix with eigenvalues λⱼ and k a number. The eigenvalues of the matrix (A − kE) are λⱼ−k.
3. If we set k close to an eigenvalue, the smallest eigenvalue of (A − kE) will be really small.
4. This can be exploited to speed up the convergence.

As I said, the shifting is already implemented, the only thing that needs to be done is to set k close to an eigenvalue. The relevant code is on this line.

[josdejong]

So this would improve performance, that's always nice!

[gwhitney]

Looking at the literature, one simple choice (that used to be commonly used, but apparently has been superseded in more recent even better algorithms) is just to take k to be the bottom right corner of the current approximation to the desired triangular matrix. I just tried this, and it passes all tests and vastly improves the precision which one can use in computing the eigenvalues of the "bad" matrix in #3036 from 1e-4 to 1e-15.

[Digression: This "tiny" precision 1e-15, which is smaller than the "standard" precision config.epsilon, does not produce any more accurate eigenvalues than config.epsilon does for this bad case; in fact config.epsilon produces the closest to the true eigenvalues, which are all exactly 2 -- at config.epsilon precision, they are off by around 5e-6. And shrinking both config.epsilon and the precision used doesn't improve the results either; in fact, the computed eigenvalues seem to be converging on fixed values slightly farther away from the "true" ones as config.epsilon and the precision parameter shrink below 1e-12. Perhaps we have simply reached the limit of what the combination of IEEE double arithmetic and this algorithm can accomplish with this numerically tricky matrix. Perhaps unsurprisingly, though, if one simply "guesses" the correct shift k=2 for this matrix, it not only converges right away to essentially the correct values, it also realizes that the matrix is defective, which it never does with the simple choice of the bottom right corner for k, at any precision level. Didn't try with bignumbers though; it might stand a chance there.]

Anyhow, all this suggests that I should add the simple one-line change from k=0 to k = (bottom right corner of current approximation) to the current PR #3037, even though it's actually orthogonal to the issue that PR addresses, and then close this issue because it seems to be a big improvement to the convergence of the algorithm, for almost no time investment on our part, and looking into superior methods would presumably require significant immersion in the literature on this, time which might be better spent just using software that others have optimized for eigenvalue/eigenvector problems. Or I could do a separate subsequent tiny PR with this small change once #3037 is merged but still before v12, if you prefer. Let me know.

[josdejong]

That is awesome news @gwhitney 😎

If it's such a small change in the code, feel free to directly apply it in PR #3037.

Is there still a reason to specify a custom precision for eigs with this improvement, or can we deprecate that?

Would this in itself be a breaking change though? (Since you're discussing it in the v12 breaking changes topic?)

[gwhitney]

OK will do that after I respond to all the code review in #3037. Yes we still need precision, convergence in eigenvalue problems is always touchy. My understanding is that with the best modern methods (which we have definitely not yet implemented, this simple shift is now very old-school) it's possible to get all 4x4 matrices to converge in double precision in about 40 iterations. But for bigger matrices things could still go weirdly, I believe. So some explicit control is worthwhile.

This is not a breaking change, just an improvement of results, but it only works on top of #3037, which is breaking.

[gwhitney]

Oh and on precision: it's also totally necessary with bignumbers, and I think always will be.

[josdejong]

Thanks!

OK all is clear to me 👍