Fix GMRES restart residual reconstruction by po-nuvai · Pull Request #586 · PyFR/PyFR

po-nuvai · 2026-08-14T06:52:03Z

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.ifwill_restart:
self._add(1/h_jp1_j, v[j+1]) # normalise the last Arnoldi vectorz=np.zeros(j+2)
z[j+1] =1.0foriinrange(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 ...ifnotwill_restart:
breakrnorm=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 · 2026-08-14T07:14:36Z

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 · 2026-08-14T10:20:37Z

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 · 2026-08-14T14:42:15Z

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

fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:11Z

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:
fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:23Z

+ # preconditioner is active. The residual is
+ # r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2}), i.e. a combination of all
+ # of the Arnoldi vectors v_0..v_{j+1}; stash the result in v[j+1],
+ # which survives the solution update.


Aim for single line comments.

The restarted-GMRES restart was reconstructing the residual as being parallel to the last Arnoldi vector only, dropping the components in the earlier Krylov directions. The correct residual is a combination of all the Arnoldi vectors, obtained by cascading the Givens rotations back through the basis (r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2})). The old form had the correct magnitude but wrong direction, causing restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

po-nuvai · 2026-08-14T15:26:13Z

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

foriinrange(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 · 2026-08-14T18:30:19Z

+ z[j + 1] = 1.0
+ 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])


Can drop the ()

FreddieWitherden · 2026-08-14T18:31:55Z

+ 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])
+ z *= np.copysign(1.0, self._beta[j + 1])


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

FreddieWitherden · 2026-08-16T00:41:05Z

 y = np.linalg.solve(self._H[:j + 1, :j + 1], self._beta[:j + 1])

+ # Determine if a restart is needed after this cycle
+ will_restart = not (err < rtol or h_jp1_j < self._breakdown_tol


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

done maybe?

FreddieWitherden · 2026-08-28T11:33:07Z

Are you okay if I get this upstream w/corrections?

po-nuvai · 2026-08-14T06:52:03Z

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.ifwill_restart:
self._add(1/h_jp1_j, v[j+1]) # normalise the last Arnoldi vectorz=np.zeros(j+2)
z[j+1] =1.0foriinrange(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 ...ifnotwill_restart:
breakrnorm=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 · 2026-08-14T07:14:36Z

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 · 2026-08-14T10:20:37Z

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 · 2026-08-14T14:42:15Z

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

fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:11Z

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:
fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:23Z

+ # preconditioner is active. The residual is
+ # r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2}), i.e. a combination of all
+ # of the Arnoldi vectors v_0..v_{j+1}; stash the result in v[j+1],
+ # which survives the solution update.


Aim for single line comments.

The restarted-GMRES restart was reconstructing the residual as being parallel to the last Arnoldi vector only, dropping the components in the earlier Krylov directions. The correct residual is a combination of all the Arnoldi vectors, obtained by cascading the Givens rotations back through the basis (r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2})). The old form had the correct magnitude but wrong direction, causing restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

po-nuvai · 2026-08-14T15:26:13Z

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

foriinrange(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 · 2026-08-14T18:30:19Z

+ z[j + 1] = 1.0
+ 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])


Can drop the ()

FreddieWitherden · 2026-08-14T18:31:55Z

+ 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])
+ z *= np.copysign(1.0, self._beta[j + 1])


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

FreddieWitherden · 2026-08-16T00:41:05Z

 y = np.linalg.solve(self._H[:j + 1, :j + 1], self._beta[:j + 1])

+ # Determine if a restart is needed after this cycle
+ will_restart = not (err < rtol or h_jp1_j < self._breakdown_tol


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

done maybe?

FreddieWitherden · 2026-08-28T11:33:07Z

Are you okay if I get this upstream w/corrections?

po-nuvai · 2026-08-14T06:52:03Z

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.ifwill_restart:
self._add(1/h_jp1_j, v[j+1]) # normalise the last Arnoldi vectorz=np.zeros(j+2)
z[j+1] =1.0foriinrange(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 ...ifnotwill_restart:
breakrnorm=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 · 2026-08-14T07:14:36Z

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 · 2026-08-14T10:20:37Z

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 · 2026-08-14T14:42:15Z

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

fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:11Z

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:
fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:23Z

+ # preconditioner is active. The residual is
+ # r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2}), i.e. a combination of all
+ # of the Arnoldi vectors v_0..v_{j+1}; stash the result in v[j+1],
+ # which survives the solution update.


Aim for single line comments.

The restarted-GMRES restart was reconstructing the residual as being parallel to the last Arnoldi vector only, dropping the components in the earlier Krylov directions. The correct residual is a combination of all the Arnoldi vectors, obtained by cascading the Givens rotations back through the basis (r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2})). The old form had the correct magnitude but wrong direction, causing restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

po-nuvai · 2026-08-14T15:26:13Z

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

foriinrange(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 · 2026-08-14T18:30:19Z

+ z[j + 1] = 1.0
+ 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])


Can drop the ()

FreddieWitherden · 2026-08-14T18:31:55Z

+ 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])
+ z *= np.copysign(1.0, self._beta[j + 1])


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

FreddieWitherden · 2026-08-16T00:41:05Z

 y = np.linalg.solve(self._H[:j + 1, :j + 1], self._beta[:j + 1])

+ # Determine if a restart is needed after this cycle
+ will_restart = not (err < rtol or h_jp1_j < self._breakdown_tol


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

done maybe?

FreddieWitherden · 2026-08-28T11:33:07Z

Are you okay if I get this upstream w/corrections?

po-nuvai · 2026-08-14T06:52:03Z

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.ifwill_restart:
self._add(1/h_jp1_j, v[j+1]) # normalise the last Arnoldi vectorz=np.zeros(j+2)
z[j+1] =1.0foriinrange(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 ...ifnotwill_restart:
breakrnorm=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 · 2026-08-14T07:14:36Z

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 · 2026-08-14T10:20:37Z

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 · 2026-08-14T14:42:15Z

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

fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:11Z

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:
fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:23Z

+ # preconditioner is active. The residual is
+ # r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2}), i.e. a combination of all
+ # of the Arnoldi vectors v_0..v_{j+1}; stash the result in v[j+1],
+ # which survives the solution update.


Aim for single line comments.

The restarted-GMRES restart was reconstructing the residual as being parallel to the last Arnoldi vector only, dropping the components in the earlier Krylov directions. The correct residual is a combination of all the Arnoldi vectors, obtained by cascading the Givens rotations back through the basis (r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2})). The old form had the correct magnitude but wrong direction, causing restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

po-nuvai · 2026-08-14T15:26:13Z

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

foriinrange(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 · 2026-08-14T18:30:19Z

+ z[j + 1] = 1.0
+ 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])


Can drop the ()

FreddieWitherden · 2026-08-14T18:31:55Z

+ 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])
+ z *= np.copysign(1.0, self._beta[j + 1])


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

FreddieWitherden · 2026-08-16T00:41:05Z

 y = np.linalg.solve(self._H[:j + 1, :j + 1], self._beta[:j + 1])

+ # Determine if a restart is needed after this cycle
+ will_restart = not (err < rtol or h_jp1_j < self._breakdown_tol


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

done maybe?

FreddieWitherden · 2026-08-28T11:33:07Z

Are you okay if I get this upstream w/corrections?

po-nuvai · 2026-08-14T06:52:03Z

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.ifwill_restart:
self._add(1/h_jp1_j, v[j+1]) # normalise the last Arnoldi vectorz=np.zeros(j+2)
z[j+1] =1.0foriinrange(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 ...ifnotwill_restart:
breakrnorm=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 · 2026-08-14T07:14:36Z

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 · 2026-08-14T10:20:37Z

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 · 2026-08-14T14:42:15Z

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

fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:11Z

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:
fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:23Z

+ # preconditioner is active. The residual is
+ # r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2}), i.e. a combination of all
+ # of the Arnoldi vectors v_0..v_{j+1}; stash the result in v[j+1],
+ # which survives the solution update.


Aim for single line comments.

The restarted-GMRES restart was reconstructing the residual as being parallel to the last Arnoldi vector only, dropping the components in the earlier Krylov directions. The correct residual is a combination of all the Arnoldi vectors, obtained by cascading the Givens rotations back through the basis (r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2})). The old form had the correct magnitude but wrong direction, causing restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

po-nuvai · 2026-08-14T15:26:13Z

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

foriinrange(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 · 2026-08-14T18:30:19Z

+ z[j + 1] = 1.0
+ 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])


Can drop the ()

FreddieWitherden · 2026-08-14T18:31:55Z

+ 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])
+ z *= np.copysign(1.0, self._beta[j + 1])


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

FreddieWitherden · 2026-08-16T00:41:05Z

 y = np.linalg.solve(self._H[:j + 1, :j + 1], self._beta[:j + 1])

+ # Determine if a restart is needed after this cycle
+ will_restart = not (err < rtol or h_jp1_j < self._breakdown_tol


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

done maybe?

FreddieWitherden · 2026-08-28T11:33:07Z

Are you okay if I get this upstream w/corrections?

po-nuvai · 2026-08-14T06:52:03Z

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.ifwill_restart:
self._add(1/h_jp1_j, v[j+1]) # normalise the last Arnoldi vectorz=np.zeros(j+2)
z[j+1] =1.0foriinrange(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 ...ifnotwill_restart:
breakrnorm=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 · 2026-08-14T07:14:36Z

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 · 2026-08-14T10:20:37Z

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 · 2026-08-14T14:42:15Z

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

fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:11Z

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:
fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:23Z

+ # preconditioner is active. The residual is
+ # r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2}), i.e. a combination of all
+ # of the Arnoldi vectors v_0..v_{j+1}; stash the result in v[j+1],
+ # which survives the solution update.


Aim for single line comments.

The restarted-GMRES restart was reconstructing the residual as being parallel to the last Arnoldi vector only, dropping the components in the earlier Krylov directions. The correct residual is a combination of all the Arnoldi vectors, obtained by cascading the Givens rotations back through the basis (r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2})). The old form had the correct magnitude but wrong direction, causing restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

po-nuvai · 2026-08-14T15:26:13Z

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

foriinrange(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 · 2026-08-14T18:30:19Z

+ z[j + 1] = 1.0
+ 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])


Can drop the ()

FreddieWitherden · 2026-08-14T18:31:55Z

+ 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])
+ z *= np.copysign(1.0, self._beta[j + 1])


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

FreddieWitherden · 2026-08-16T00:41:05Z

 y = np.linalg.solve(self._H[:j + 1, :j + 1], self._beta[:j + 1])

+ # Determine if a restart is needed after this cycle
+ will_restart = not (err < rtol or h_jp1_j < self._breakdown_tol


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

done maybe?

FreddieWitherden · 2026-08-28T11:33:07Z

Are you okay if I get this upstream w/corrections?

po-nuvai · 2026-08-14T06:52:03Z

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.ifwill_restart:
self._add(1/h_jp1_j, v[j+1]) # normalise the last Arnoldi vectorz=np.zeros(j+2)
z[j+1] =1.0foriinrange(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 ...ifnotwill_restart:
breakrnorm=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 · 2026-08-14T07:14:36Z

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 · 2026-08-14T10:20:37Z

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 · 2026-08-14T14:42:15Z

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

fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:11Z

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:
fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:23Z

+ # preconditioner is active. The residual is
+ # r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2}), i.e. a combination of all
+ # of the Arnoldi vectors v_0..v_{j+1}; stash the result in v[j+1],
+ # which survives the solution update.


Aim for single line comments.

The restarted-GMRES restart was reconstructing the residual as being parallel to the last Arnoldi vector only, dropping the components in the earlier Krylov directions. The correct residual is a combination of all the Arnoldi vectors, obtained by cascading the Givens rotations back through the basis (r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2})). The old form had the correct magnitude but wrong direction, causing restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

po-nuvai · 2026-08-14T15:26:13Z

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

foriinrange(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 · 2026-08-14T18:30:19Z

+ z[j + 1] = 1.0
+ 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])


Can drop the ()

FreddieWitherden · 2026-08-14T18:31:55Z

+ 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])
+ z *= np.copysign(1.0, self._beta[j + 1])


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

FreddieWitherden · 2026-08-16T00:41:05Z

 y = np.linalg.solve(self._H[:j + 1, :j + 1], self._beta[:j + 1])

+ # Determine if a restart is needed after this cycle
+ will_restart = not (err < rtol or h_jp1_j < self._breakdown_tol


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

done maybe?

FreddieWitherden · 2026-08-28T11:33:07Z

Are you okay if I get this upstream w/corrections?

po-nuvai · 2026-08-14T06:52:03Z

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.ifwill_restart:
self._add(1/h_jp1_j, v[j+1]) # normalise the last Arnoldi vectorz=np.zeros(j+2)
z[j+1] =1.0foriinrange(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 ...ifnotwill_restart:
breakrnorm=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 · 2026-08-14T07:14:36Z

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 · 2026-08-14T10:20:37Z

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 · 2026-08-14T14:42:15Z

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

fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:11Z

Done — switched the cascade to slice assignment and enumerate(zip(...)) over the c/s pairs:
fork, (c, s) inenumerate(zip(self._cs[j::-1], self._sn[j::-1])):
i=j-kz[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 · 2026-08-14T15:12:23Z

+ # preconditioner is active. The residual is
+ # r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2}), i.e. a combination of all
+ # of the Arnoldi vectors v_0..v_{j+1}; stash the result in v[j+1],
+ # which survives the solution update.


Aim for single line comments.

The restarted-GMRES restart was reconstructing the residual as being parallel to the last Arnoldi vector only, dropping the components in the earlier Krylov directions. The correct residual is a combination of all the Arnoldi vectors, obtained by cascading the Givens rotations back through the basis (r_m = beta[j+1]*V_{j+2}*(G^T e_{j+2})). The old form had the correct magnitude but wrong direction, causing restarted GMRES (restart < linear-max-iter) to stagnate or break down instead of converging.

po-nuvai · 2026-08-14T15:26:13Z

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

foriinrange(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 · 2026-08-14T18:30:19Z

+ z[j + 1] = 1.0
+ 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])


Can drop the ()

FreddieWitherden · 2026-08-14T18:31:55Z

+ 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])
+ z *= np.copysign(1.0, self._beta[j + 1])


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

FreddieWitherden · 2026-08-16T00:41:05Z

 y = np.linalg.solve(self._H[:j + 1, :j + 1], self._beta[:j + 1])

+ # Determine if a restart is needed after this cycle
+ will_restart = not (err < rtol or h_jp1_j < self._breakdown_tol


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

done maybe?

FreddieWitherden · 2026-08-28T11:33:07Z

Are you okay if I get this upstream w/corrections?

po-nuvai force-pushed the fix/gmres-restart branch 2 times, most recently from 7894c26 to fd9a7f9Compare August 14, 2026 07:12

po-nuvai force-pushed the fix/gmres-restart branch 2 times, most recently from b0eb448 to 80463c5Compare August 14, 2026 14:41