Inconsistent meaning of success in MINRES and LSMR for systems without a solution #629

For a linear system A x = b which has no solution (singular A, with b outside its range), MinimumResidual and LSMR in psydac.linalg.solvers disagree on info['success']:

import numpy as np
from psydac.ddm.cart import DomainDecomposition, CartDecomposition
from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix
from psydac.linalg.solvers import inverse

n, p = 6, 1
ddm  = DomainDecomposition([n - p], [False])
cart = CartDecomposition(ddm, [n], [np.array([0])], [np.array([n - 1])], [p], [1])
V    = StencilVectorSpace(cart)

A = StencilMatrix(V, V)        # A = diag(1, 1, 1, 1, 1, 0)
for i in range(n - 1):
    A[i, 0] = 1.0
b = V.zeros()
b[:] = 1.0                     # b is not in the range of A

for solver in ['MINRES', 'LSMR']:
    S = inverse(A, solver, tol=1e-10, maxiter=100)
    x = S @ b
    r = b - A @ x
    print(solver, S.get_info()['success'], S.get_info()['res_norm'], np.sqrt(r.inner(r).real))

MINRES False 0.9165... 1.0408...
LSMR   True  1.0       1.0

• LSMR reports success because it found a least-squares solution, even though A x = b has no solution. This follows SciPy, whose lsmr stops with istop = 2 ("A least-squares solution was found").
• MinimumResidual reports failure because its estimate of the residual is not below tol. However, this estimate is not reliable: here it is 0.92 while the true residual is 1.04, and for a zero operator (A = 0) it collapses to 0 at the first iteration, so that MINRES reports success with x = 0 (with the fixes of PR Fix several bugs in iterative linear solvers #604; before, this case raises a NameError). Note that SciPy's minres also returns info = 0 (success) for a zero operator.

A user who solves A x = b cannot tell from success whether the system was actually solved.

Possible solutions

1. Keep the SciPy semantics and document them: success means that the stopping criteria were met, which for MINRES and LSMR may be a least-squares solution
2. Base success on the true residual at the end (‖b - A x‖ < tol ‖b‖), at the cost of one extra matrix-vector product per solve
3. Keep success as it is, and add the reason for stopping to info (e.g. istop and its message, as in SciPy), so that users can tell a solution from a least-squares solution