Inconsistent meaning of niter and maxiter in iterative solvers #630

Problem

The value niter returned by get_info() does not mean the same thing in every solver of psydac.linalg.solvers. The meaning of maxiter differs as well:

Hence:

• With the same maxiter, CG performs one iteration fewer than the other solvers.
• The values of niter cannot be compared between solvers.
• In GMRES, niter changes meaning depending on whether the solver stopped because of tol or because of maxiter.

Reproducer

The matrix is symmetric positive definite and has only two distinct eigenvalues, so every solver reaches the exact solution in 2 iterations:

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 = 8, 1
domain_decomposition = DomainDecomposition([n - p], [False])
cart = CartDecomposition(domain_decomposition, [n], [np.array([0])], [np.array([n - 1])], [p], [1])
V = StencilVectorSpace(cart)

# Block diagonal matrix with blocks [[2, 1], [1, 2]]: symmetric positive definite
# with eigenvalues 1 and 3 only, hence the exact solution is reached in 2 iterations
A = StencilMatrix(V, V)
A[:, 0] = 2
A[0:n:2, 1] = 1
A[1:n:2, -1] = 1
A.remove_spurious_entries()

b = V.zeros()
b[0:n] = np.arange(1, n + 1)

for solver in ['CG', 'BiCG', 'BiCGSTAB', 'MINRES', 'LSMR', 'GMRES']:
    for maxiter in [2, 100]:
        A_inv = inverse(A, solver, tol=1e-10, maxiter=maxiter)
        x = A_inv @ b
        info = A_inv.get_info()
        print(f"{solver:8s} maxiter={maxiter:3d}: niter={info['niter']}, success={info['success']}")

Output:

CG       maxiter=  2: niter=2, success=False
CG       maxiter=100: niter=3, success=True
BiCG     maxiter=  2: niter=2, success=True
BiCG     maxiter=100: niter=2, success=True
BiCGSTAB maxiter=  2: niter=2, success=True
BiCGSTAB maxiter=100: niter=2, success=True
MINRES   maxiter=  2: niter=2, success=True
MINRES   maxiter=100: niter=2, success=True
LSMR     maxiter=  2: niter=2, success=True
LSMR     maxiter=100: niter=2, success=True
GMRES    maxiter=  2: niter=2, success=True
GMRES    maxiter=100: niter=3, success=True

• CG needs maxiter=3 to perform its 2 iterations, and reports niter=3.
• GMRES performs 2 iterations in both cases, but reports niter=2 or niter=3. (With maxiter=2, its success=True is also wrong: the returned x has a true residual of about 0.94, because of the off-by-one bug listed in Several bugs in iterative solvers #603.)

History

• 2018: CG was added in e9a0c01 with niter = iterations performed, the same meaning as in the other solvers today.
• 2021: FEEC with Dirichlet BCs. #89 (a0c432b) started printing the initial residual as row 1 of the verbose table, and shifted the loops of cg and pcg to start at 2. As a result, niter became the number of matrix-vector products (one for A x0, plus one per iteration), and maxiter allowed one iteration fewer.
• 2023: New Linear Operators #261 and Make linear algebra layer compatible with complex numbers #273 kept this convention in the class-based solvers, and GMRES (Add GMRES solver #301) followed it only partly.

Proposed solution

Use the meaning of the other four solvers everywhere, as in SciPy:

• niter is the number of iterations performed (0 if x0 already satisfies tol).
• maxiter is the maximum number of iterations.
• If the number of matrix-vector products is useful, it could be reported separately in get_info(), since it depends on the solver: 1 + niter for CG, MINRES and GMRES, 1 + 2 niter for BiCG and BiCGSTAB, 2 + 2 niter for LSMR.

This is intentionally a change of behaviour, to be listed under "Changed" in the CHANGELOG:

• CG: the same maxiter allows one more iteration, and niter is 1 lower.
• GMRES: the iterations performed are unchanged, but niter is 1 lower when GMRES converges before maxiter.

The verbose tables can keep printing the initial residual, labelled as iteration 0.

Notes

This issue was left out of scope in #603.