Add check_finite and check_cond in NumpyMatrixOperator and use SciPy solvers elsewhere

[pmli]

Closes #1052.

Since #1603, NumpyMatrixOperator caches the LU factorization and checks the condition number. To enable faster apply_inverse, especially for small matrices, this PR adds check_finite and check_cond default in NumpyMatrixOperator.apply_inverse.

Regarding the choice between NumPy and SciPy linalg routines, it seems that SciPy is used predominantly (see the two greps below, which was enabled by forbidding import scipy and from numpy.linalg import). Therefore, this PR replaces NumPy's solve, lstsq, inv, qr, eigvals by SciPy's (in src/).

---

NumpyMatrixOperator uses lu_factor and lu_solve.

Places where numpy.linalg.solve is used:

$ grep -rn 'np.linalg.solve(' src/
src/pymortests/operators/operators.py:503:    u = np.linalg.solve(O.T, v.ravel())
src/pymor/reductors/spectral_factor.py:113:        L = np.linalg.solve(M.T, C - B.T @ X @ E)
src/pymordemos/phlti.py:114:    Q = np.linalg.solve(J - R, A)

Places where scipy.linalg.solve is used:

$ grep -rn 'spla.solve(' src/
src/pymortests/algorithms/riccati.py:65:        quadratic = quadratic @ spla.solve(R, quadratic.T)
src/pymortests/algorithms/to_matrix.py:116:    assert_type_and_allclose(spla.solve(S, A.T.dot(R)), Aadj, 'dense')
src/pymortests/algorithms/to_matrix.py:122:    assert_type_and_allclose(spla.solve(S, A.T.dot(R)), Aadj, 'dense')
src/pymortests/algorithms/to_matrix.py:128:    assert_type_and_allclose(spla.solve(S, A.T.dot(R)), Aadj, 'dense')
src/pymortests/algorithms/to_matrix.py:134:    assert_type_and_allclose(spla.solve(S, A.T.dot(R)), Aadj, 'sparse')
src/pymortests/algorithms/to_matrix.py:229:    assert_type_and_allclose(L @ spla.solve(C, R.T), LR, 'dense')
src/pymortests/operators/block.py:92:    assert np.allclose(spla.solve(C, v), w)
src/pymortests/operators/block.py:112:    assert np.allclose(spla.solve(C.T, v), w)
src/pymortests/operators/low_rank_op.py:40:                       L.to_numpy().T @ spla.solve(C, R.to_numpy() @ U.to_numpy().T))
src/pymortests/operators/low_rank_op.py:53:                       R.to_numpy().T @ spla.solve(C.T, L.to_numpy() @ V.to_numpy().T))
src/pymortests/operators/low_rank_op.py:63:    assert np.allclose(U.to_numpy().T, spla.solve(mat, V.to_numpy().T))
src/pymortests/operators/low_rank_op.py:68:    mat = A.matrix + L.to_numpy().T @ spla.solve(C, R.to_numpy())
src/pymortests/operators/low_rank_op.py:69:    assert np.allclose(U.to_numpy().T, spla.solve(mat, V.to_numpy().T))
src/pymortests/operators/low_rank_op.py:79:    assert np.allclose(V.to_numpy().T, spla.solve(mat.T, U.to_numpy().T))
src/pymortests/operators/low_rank_op.py:84:    mat = A.matrix + L.to_numpy().T @ spla.solve(C, R.to_numpy())
src/pymortests/operators/low_rank_op.py:85:    assert np.allclose(V.to_numpy().T, spla.solve(mat.T, U.to_numpy().T))
src/pymor/reductors/sobt.py:382:            W1TV1invW1TV2 = spla.solve(self.W1.inner(self.V1), self.W1.inner(self.V2))
src/pymor/reductors/h2.py:789:        b = spla.solve(EX, B)
src/pymor/algorithms/lrradi.py:148:            ImBNKL = spla.solve(ImBN, B.inner(L))
src/pymor/algorithms/basic.py:95:        coeffs = spla.solve(gramian, rhs, assume_a='pos', overwrite_a=True, overwrite_b=True).T
src/pymor/algorithms/to_matrix.py:126:                res = spla.solve(source_product, res)
src/pymor/algorithms/to_matrix.py:183:            res = op.left.to_numpy().T @ spla.solve(op.core, op.right.to_numpy())
src/pymor/operators/constructions.py:499:            V = spla.solve(self.core, V)
src/pymor/operators/constructions.py:508:            U = spla.solve(self.core.T.conj(), U)
src/pymor/operators/constructions.py:572:        U.axpy(-beta, AinvL.lincomb(spla.solve(mat, RhAinvV).T))
src/pymor/operators/constructions.py:590:        V.axpy(-beta, AinvhR.lincomb(spla.solve(mat, LhAinvhU).T))
src/pymor/bindings/pymess.py:499:        RinvC = spla.solve(R, C) if R is not None else C
src/pymor/bindings/pymess.py:502:        RinvBT = spla.solve(R, B.T) if R is not None else B.T
src/pymor/bindings/slycot.py:239:                G = B @ spla.solve(R, B.T)
src/pymor/bindings/slycot.py:246:                G = C.T @ spla.solve(R, C)
src/pymor/models/iosys.py:3121:    b = spla.solve(EX, B)
$ grep -rn 'from scipy.linalg import.*solve' src/
src/pymor/algorithms/ei.py:17:from scipy.linalg import solve
src/pymor/analyticalproblems/functions.py:8:from scipy.linalg import solve, solve_triangular
src/pymor/operators/numpy.py:24:from scipy.linalg import lu_factor, lu_solve
src/pymor/operators/ei.py:8:from scipy.linalg import solve, solve_triangular
src/pymor/bindings/scipy.py:7:from scipy.linalg import solve, solve_continuous_are, solve_continuous_lyapunov, solve_discrete_lyapunov

[sdrave]

I disagree with the new import conventions imposed by the PR. Could you move this part into a separate PR?

[pmli]

I disagree with the new import conventions imposed by the PR. Could you move this part into a separate PR?

Two of them (banning from numpy.linalg and import scipy) make it easier to find which solve (and lstsq) method is used where. Is it ok to at least keep those?

[sdrave]

I'm fine with the changes. However, we should add a bit more documentation:

• Add a comment to the docstring of NumpyMatrixOperator discussing that disabling the condition number check might be favorable for performance in the ROM case.
• Document somewhere why using from scipy.linalg import foo is forbidden.
• Better explain in the first PR comment what has been changed.
• Choose better title for PR.