Problem in Beyn Contour Rank Revealing

[wrs28]

The way that contour_beyn deals with more or fewer eigenvalues inside the contour than asked for leads to inaccurate answers.

As written, it seems that if you ask for k eigenvalues and the contour actually encloses neigs<k, then all k eigenvalues of the matrix B are returned, and are not marked as fictitious. Furthermore, it seems that the warning saying a rank drop is detected triggers in this scenario, and it encourages a smaller k. However, asking for k>neigs is actually not a problem so long as the extraneous eigenvalues/vectors are either not returned or marked as such.

Next, if k<neigs, the algorithm is in principle inaccurate, but there is to my knowledge no way to diagnose this from the results. If k=neigs than the algorithm should be accurate, but again there is no way to distinguish this from the k<neigs case just by looking at the answer.

One way to ensure the algorithm is accurate is to reveal the rank in the SVD by looking for negligible singular values, which indicates that k>neigs which is the limit where the algorithm works. This is exactly what the current warning discourages.

I think an easy fix would be to look for the small singular values (which is what the warning is currently triggered by), and to throw them away, or perhaps later mark them with NaNs, and to reverse the sense of the warning. In other words, the new warning would be triggered if the rank is not revealed and suggest the user increase k.

[jarlebring]

Thanks. Indeed, the recommendation / warning seems off. I agree. Returning incorrect results is bad, and should be avoidable in this case. In order to reliably investigate this issue, it would be nice to first implement the correct eigenvector extraction (from the SVD as described in Beyn's paper) rather than our current "naive" eigenvector extraction 4f32fed. This will anyway be necessary if we continue with #154 and remove compute_eigvec_from_eigval.

[jarlebring]

Wish-list / todo-list:

1. Introduce new kwarg: Neigs=no'f wanted eigvals and set default number of columns to k=Neigs+1.
2. Let p be the number of (essentially) non-zero singvals.
3. If Neigs>p, we did not find enough eigvals. Generate warning: Insufficient eigvals found (and suggest to try to increase the domain)
4. Do errmeasure check on all eigpairs. Let Neigs_found be the nof okay eigpairs.
5. No rank drop detected <=> p = k. (Note: By construction p<= k): if p=k and Neigs_found<Neigs, generate warning and suggest to increase k.

I think 4 should be made optional through a kwarg, since I fear it might get computationally demanding in an HPC-setting (depending on problem and which errmeasure we use).