Fix GMRES restart residual reconstruction

[po-nuvai]

Summary

Fixes the restarted-GMRES residual reconstruction in pyfr/integrators/implicit/krylov/gmres.py. Closes #585.

The restart was reconstructing the residual as parallel to the last Arnoldi vector only, which has the correct magnitude but the wrong direction. In exact arithmetic the residual after a GMRES cycle is

r_m = β_m · V_{m+1} · (G₀ᵀ G₁ᵀ … G_{m-1}ᵀ e_{m+1})

i.e. a linear combination of all the Arnoldi vectors v_0 … v_m, with coefficients cascading through every Givens rotation — not just v_m. Dropping the earlier components caused restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

The fix

Cascade the Givens rotations back through the Krylov basis to reconstruct the full residual direction. It reuses the already-stored cs/sn rotations and Krylov vectors, so it costs no extra matvec:

# reconstruct the residual direction *before* the solution update, because
# the preconditioned update reuses v[0] as scratch for M⁻¹(V y) and would
# otherwise clobber it; stash the result in v[j+1], which survives.
if will_restart:
    self._add(1 / h_jp1_j, v[j + 1])          # normalise the last Arnoldi vector

    z = np.zeros(j + 2)
    z[j + 1] = 1.0
    for i in range(j, -1, -1):
        c, s = self._cs[i], self._sn[i]
        z[i], z[i + 1] = c*z[i] - s*z[i + 1], s*z[i] + c*z[i + 1]
    z *= np.copysign(1.0, self._beta[j + 1])

    self._addv([z[j + 1], *z[:j + 1].tolist()], [v[j + 1], *v[:j + 1]])

# ... solution update ...

if not will_restart:
    break

rnorm = abs(self._beta[j + 1])
self._add(0, v[0], 1, v[j + 1])               # install the reconstructed residual

The one subtlety worth calling out: the reconstruction has to happen before the solution update, not after. The preconditioned solution update reuses v[0] as a scratch register for M⁻¹(V y), which would destroy the first Arnoldi vector that the cascade needs as a source.

Validation

Replicating the algorithm in NumPy (MGS/CGS Arnoldi + incremental Givens) and comparing against scipy.sparse.linalg.gmres on a well-conditioned SPD matrix (cond ≈ 8), targeting rtol = 1e-10, with and without a right preconditioner:

restart m	old	fix (no precond)	fix (preconditioned)	scipy
2	1.3e-01	9.2e-11	5.0e-11	9.2e-11
4	3.3e-02	6.5e-11	3.6e-11	6.5e-11
8	1.0e-03	8.8e-11	5.2e-11	8.8e-11
15	3.9e-06	1.0e-10	9.8e-11	1.0e-10

Notes

• Masked by default — solver-gmres restart defaults to 0 (single cycle), so the buggy path only triggers when restart < linear-max-iter.
• Independent of the Arnoldi variant (cgs vs mgs) and of preconditioning.

[po-nuvai]

Reviewed this change. The math is right — the old restart kept only the v[j+1] component of the residual (correct magnitude, wrong direction), and cascading the Givens rotations back through the basis is the proper reconstruction of r_m = β_m V_{m+1} (Gᵀ e_{m+2}).

One thing I caught while going back over it, and have now fixed in the branch: the reconstruction has to happen before the solution update, not after. With a preconditioner active the update reuses v[0] as scratch for M⁻¹(V y), which destroys the first Arnoldi vector that the cascade reads as a source. So the residual direction is now computed ahead of the update, stashed in v[j+1] (which survives), and installed into v[0] afterwards.

Non-blocking nits, if you want them:

• the cascade tuple has redundant outer parens (z[i], z[i + 1] = (...)).
• np.zeros(j + 2) could be dtype=self._beta.dtype for explicitness; it defaults to float64 which matches _cs/_sn/_beta anyway.

Numerically I replicated the algorithm and it matches scipy's restarted GMRES to ~1e-10 with and without a right preconditioner (the old code sat at 1e-1 … 1e-6 and stagnated/breakdown'd). Default config is unaffected since restart defaults to 0 → single cycle.

[FreddieWitherden]

Overall looks good. Can probably clean up the z assignments a bit z[i:i+1] = and see if it is easier to enumerate(zip(...)) over the c and s terms.

[po-nuvai]

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:

for k, (c, s) in enumerate(zip(self._cs[j::-1], self._sn[j::-1])):
    i = j - k
    z[i:i + 2] = (c*z[i] - s*z[i + 1], s*z[i] + c*z[i + 1])

Reads a bit cleaner, and it's still the same backward cascade (applies G_j^T … G_0^T in order). I re-checked the equivalence numerically — identical to machine precision against the previous form.

[FreddieWitherden]

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:

for k, (c, s) in enumerate(zip(self._cs[j::-1], self._sn[j::-1])):
    i = j - k
    z[i:i + 2] = (c*z[i] - s*z[i + 1], s*z[i] + c*z[i + 1])

Reads a bit cleaner, and it's still the same backward cascade (applies G_j^T … G_0^T in order). I re-checked the equivalence numerically — identical to machine precision against the previous form.

I think it ends up a bit uglier in practice due to the extra line breaks. Maybe just revert and just tweak the assignment.

[FreddieWitherden]

Aim for single line comments.

[po-nuvai]

Agreed on both — reverted the enumerate(zip(...)) (the line breaks made it worse) and just kept the slice-assignment tweak:

for i in range(j, -1, -1):
    c, s = self._cs[i], self._sn[i]
    z[i:i + 2] = (c*z[i] - s*z[i + 1], s*z[i] + c*z[i + 1])

Also condensed that comment block down to a single line. :)

[FreddieWitherden]

Can drop the ()

[FreddieWitherden]

Do we need copysign here or will np.sign work?

[FreddieWitherden]

Given we not this in one if later on, we may want to remove the not and rename the variable accordingly.

[FreddieWitherden]

done maybe?