Samdp likes to fail with InversionError when theta approaches pole

[sdrave]

The issue pops up in CI periodically, like here.

The following code snippet reproduces the issue for me on Linux with PyPI packages:

import numpy as np
from pymortests.algorithms.samdp import conv_diff_1d_fd
from pymor.algorithms.samdp import samdp
from pymor.operators.numpy import NumpyMatrixOperator

def test_samdp(n, m, k, wanted, which, rng):
    A = conv_diff_1d_fd(n, 1, 1)
    Eop = None

    B = rng.standard_normal((n, m))
    C = rng.standard_normal((k, n))

    Aop = NumpyMatrixOperator(A)
    Bva = Aop.source.from_numpy(B.T)
    Cva = Aop.source.from_numpy(C)

    return samdp(Aop, Eop, Bva, Cva, wanted, which=which)

rng = np.random.default_rng(595)
test_samdp(50, 1, 3, 15, 'NR', rng)

@lbalicki, any ideas?

[sdrave]

At the moment, CI failures only appear in the coda tests on macos. I have xfailed the samdp tests and the mt tests using samdp in #2360. The respective code has to be removed when this issue has been fixed.

[lbalicki]

The reason this happens is because an exact eigenvalue is passed to _twosided_rqi. I think a first step to fix this would be to do a residual check before trying to invert the singular matrix theta * E - A for the first time. This will only fix the case where both, the eigenvalue and the eigenvector are exact. I will have to think about if there is a way to detect the case where the eigenvalue is exact but the eigenvector is not.

[lbalicki]

So it looks like there already is a residual check before starting _twosided_rqi. Passing exact eigenvalues to these types of algorithms is a tricky issue and I don't think other implementations have good ways of dealing with it either. The code below for example also fails with a "factor is exactly singular" error for me and the reason is exactly the same as in _twosided_rqi:

from scipy.sparse.linalg import eigs
from scipy.sparse import identity
eigs(identity(10),k=1,sigma=1)

I can only think of some "hacky" ways to fix this. E.g., always slightly perturb any initial guesses passed to _twosided_rqi or somehow catch the inversion error and add the perturbation to the initial eigenvalue guess then. Any preferences or other suggestions?