Pith. sign in

REVIEW 3 major objections 3 minor 26 references

Forward and backward error bounds for a mixed precision preconditioned conjugate gradient algorithm

T0 review · 3 major / 3 minor · reviewed 2026-08-04 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read The preconditioned conjugate gradient algorithm can reach backward error O(u) and forward error O(u)√κ(A) without assuming the recursive residual norm drops below machine precision.

desk verdict Genuinely new result with a real conditional gap: the proof removes the residual-smallness assumption, but the abstract's O(u) claim overstates Theorem 2's dependence on an unverified k* satisfying condition (76). read the letter →

arxiv 2510.11379 v3 pith:AD6DB4QL submitted 2025-10-13 math.NA cs.NA

Thomas Bake , Erin Carson , Yuxin Ma This is my paper · ORCID
classification math.NAcs.NA MSC 65F1065F0865G5065Y20
keywords conjugategradientpreconditioningmixedprecisionbackwarderrorforwardroundinganalysisfiniteKrylovsubspacemethods
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper proves that the preconditioned conjugate gradient (PCG) algorithm, when run in finite precision, can reach relative backward error on the order of the machine roundoff unit u and relative forward error on the order of u times the square root of the condition number of the matrix A — and it does so without the usual assumption that the recursively updated residual norm falls far below u. Previous analyses of CG-type methods took that small-recursive-residual condition as a starting point and then bounded the residual gap; this work instead uses the fact that CG minimizes a quadratic function, showing that the difference in function values at two consecutive iterates eventually becomes tiny. From that, the authors derive a quadratic inequality whose solution bounds the residual, and hence the backward and forward errors. They also show that applying preconditioners in reduced precision does not degrade these accuracy levels under stated conditions, and they introduce a split-preconditioned variant that avoids the accuracy loss that occurs when the left preconditioner is applied in low precision in the classical split algorithm.

What carries the argument

The load-bearing object is the quadratic objective f(y)=½yᵀAy−yᵀb of the linear system, together with the specific iteration index k* at which the difference f(x̄_{k*)})−f(x̄_{k*+1}) first becomes as small as O(n²(k*)²u²). Lemma 4 expresses this difference as −½α̂_k(1+δα̂_k)‖r̂_k‖²_{M^{-1}} plus a perturbation; Lemma 5 proves such a k* exists because the monotone sequence f(x̄_k) converges; Lemma 6 combines these into a quadratic inequality in ‖r̂_k‖_{M^{-1}} and solves it. The standard residual-gap bound (Theorem 1) is then used only after the recursive residual is already known to be small, not as an assumption.

What would settle it

Run PCG (with M=I for simplicity) on an SPD matrix whose eigenvalues are clustered so that the CG polynomial decreases slowly, e.g., n=100, condition number 10^8, with eigenvalues packed near the extremes; track the actual quadratic difference f(x̄_{k+1})−f(x̄_k) in floating point. Find the smallest k where this difference falls below O(n²k²u²)‖A‖²‖x‖²; if that k also satisfies n k u κ(A) > 1/2 (so condition (76) fails), then check whether the relative backward error at that iterate is indeed larger than the bound O(n k² u) max(‖x̂_j‖/‖x‖). If the backward error remains O(u) anyway, the theore

Watch

Extended reading notes

Core claim

The central result is Theorem 2: under the assumption that O(nu)κ(A)+2ε_pre ≤ 1/2 and that the specific iteration k* whose existence is guaranteed by Lemma 5 also satisfies condition (76), there is an index i ≤ k* for which the relative true-residual (backward) error satisfies ‖b−Ax̂_i‖/(‖A‖‖x‖) ≤ O(n(k*)²u) κ(M^{-1})^{1/2} max_{j≤k*+1}(‖x̂_j‖/‖x‖,1), the same bound holds for the recursively updated residual, and the A-norm forward error satisfies ‖x̂_i−x‖_A/(‖A‖^{1/2}‖x‖) ≤ O(n(k*)²u) κ(M^{-1})^{1/2} κ(A)^{1/2} max(...). These are the first bounds in this literature that reach the levels O(u) and O(u)κ(A)^{1/2} without assuming the recursively updated residual norm has already dropped below

Load-bearing premise

The proof needs the first iteration at which the quadratic function values stop changing to occur before the accumulated rounding error terms O(n k* u) grow past 1/2; Lemma 5 shows such an iteration exists but does not bound how early it happens.

Editorial extensions

If this is right

  • For any SPD system satisfying the stated conditions, PCG is provably backward stable in the O(u) sense and forward stable in the O(u)κ(A)^{1/2} sense; no separate assumption about the recursive residual norm is needed.
  • Preconditioners can be stored and applied in lower precision than the working precision without changing the asymptotic error bounds, as long as the precision-dependent parameters ε_s, ε_q, ε_z keep condition (99) satisfied; Corollary 1 gives explicit thresholds for Cholesky-factor preconditioners.
  • The split-preconditioned variant of Algorithm 2 does not suffer the accuracy degradation seen in the classical split PCG when the left preconditioner is applied in low precision, because the left-preconditioner solve is not embedded in the residual recurrence.
  • The bounds degrade polynomially with n and k*, so the theorem predicts that very slow convergence can eventually destroy the O(u) guarantee; the authors observe that in practice no failure is seen when the initial guess has a normal size.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The monotone-convergence argument suggests a practical stopping rule: monitor f(x̄_{k+1})−f(x̄_k) and stop when it drops below a threshold proportional to n²u²‖A‖²‖M^{-1}‖‖x‖²; this would directly trigger the regime where Theorem 2's bounds hold, and is implementable without forming the true residual.
  • Because condition (76) bounds k* in terms of n, u, κ(A), and κ(M), the analysis implies a trade-off between allowed iterations and the working precision: for a fixed precision, there is an effective maximum number of PCG iterations beyond which the theorem no longer guarantees accuracy. A testable extension is to construct an SPD system where the optimal CG polynomial converges so slowly that the
  • The same optimization-based proof strategy might extend to other iterative solvers with a variational formulation, such as symmetric indefinite MINRES, where analogous residual-gap assumptions are currently used.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 3 minor

Summary. The paper develops a finite-precision analysis of a unified preconditioned conjugate gradient (PCG) framework (Algorithm 2) covering left, right, and split preconditioning. The analysis derives local rounding-error bounds (Lemmas 1–3), a relation between the quadratic objective decrease and the M^{-1}-norm of the recursively updated residual (Lemma 4), an existence statement for an iteration with arbitrarily small objective decrease (Lemma 5), and a quadratic-inequality argument yielding a bound on the residual at that iteration (Lemma 6). The main result, Theorem 2, claims that under O(nu)κ(A)+2ε_pre≤1/2 and condition (76), there is i≤k* such that the true residual and A-norm error satisfy O(n(k*)^2u)κ(M^{-1})^{1/2}-type bounds. The paper also derives precision conditions for Cholesky-factor preconditioners (Lemma 7, Corollary 1) and reports numerical experiments, including comparisons with Saad's split PCG variant.

Significance. If Theorem 2 were unconditional, it would be a significant advance: it would give backward and forward error bounds for PCG at levels O(u) and O(u)κ(A)^{1/2} without the customary assumption that the recursively updated residual norm decreases below machine precision, and it would provide a theoretical basis for applying preconditioners in low precision. The proof is carefully structured, the individual lemmas are clearly stated, no fitted parameters enter the derivation, and the numerical experiments are supported by a publicly available Julia package. However, the central result is conditional on the existence of an iteration k* satisfying (76), and the paper does not prove that such an iteration exists under the stated assumptions. Since the experiments themselves show convergence of bfloat16 cases in which the analogous assumption (99) fails for every k, the scope of the theorem is substantially narrower than the abstract suggests. The contribution is valuable as a framework and as a conditional analysis, but the advertised unconditional claim is not established.

major comments (3)
  1. [§4.3, Lemmas 5–6 and Theorem 2] Theorem 2's O(n(k*)²u) bounds require the k* from Lemma 5 to satisfy condition (76), but Lemma 5 only proves the existence of a k* with f(bar x_k*)−f(bar x_{k*+1}) ≤ ε(k*)² αhat_k* together with (66); it gives no quantitative control on k*. The hypotheses O(nu)κ(A)+2ε_pre≤1/2 do not imply that any such k* satisfies (76). For M=I, (76) requires k* ≲ 1/(n u κ(A)), while the first sufficiently small CG objective decrease in exact arithmetic can occur at k* of order sqrt(κ(A)) log(1/u); for κ(A) ≫ 1/(nu) these regimes are disjoint. Thus the paper's central 'O(u)' claim is only conditional on an unverified existential bound. Please either prove that such a k* exists or state (76) as an explicit assumption in Theorem 2 and revise the abstract/conclusions accordingly. The current wording 'there exists k* ... if k* satisfies (76)' does not authorize the abstract's unconditional reading.
  2. [§5, footnote 4] The advertised improvement of Algorithm 2 over Saad's split PCG (Section 3 and Figure 2) is attributed to a claimed O(u+u_L) limitation of Saad's variant, but footnote 4 says the derivation is omitted. This claim is load-bearing for the paper's comparative conclusion. Without a proof, the comparison is only an empirical observation. Please include the derivation or explicitly state that the O(u+u_L) bound is a heuristic supported only by the experiments.
  3. [§5, Figures 1/4 and Corollary 1] The bfloat16 experiments converge even though the corresponding assumption (99) fails for every k≥1. The text mentions this, but it is not just a harmless observation: it shows that the theorem's hypotheses are not necessary and, more importantly, that the converged bfloat16 runs are not covered by Theorem 2. The manuscript should clearly state that these numerical results are outside the theorem's scope and explain what the experiments are intended to demonstrate—e.g., the sufficiency, but not necessity, of the conditions, or the practical looseness of the bounds.
minor comments (3)
  1. [Abstract and Theorem 2] The abstract's O(u) and O(u)κ(A)^{1/2} wording suppresses the factors n(k*)², κ(M^{-1})^{1/2}, and max_{j≤k*+1}(||xhat_j||/||x||,1). Please state the full scaling in the abstract or add a sentence making clear that 'O(u)' here depends polynomially on n and k*.
  2. [Keywords] The keyword 'conjugated gradient algorithm' should be 'conjugate gradient algorithm'.
  3. [Figure 3 caption] The caption says 'first row' and 'second row', but the figure shows three separate heatmaps (a), (b), (c). Please correct the caption to match the actual layout.

Circularity Check

0 steps flagged · score 1.0 of 10

No circular derivation: the main theorem is a conditional forward/backward error bound built from external rounding-error and residual-gap results; the unproved iteration condition (76) is a hypothesis, not a fitted or self-referential input.

full rationale

The proof chain does not reduce to its own inputs. Theorem 2's residual and error bounds are derived from Lemmas 1-6, which in turn use standard rounding-error estimates (Higham), Greenbaum's residual-gap bound [9], and a monotone-convergence existence argument for an iteration k*. No parameter is fitted to data and then reported as a prediction; the O(.) constants are left unspecified in the usual rounding-analysis sense. The self-citations [3,5] appear only in the literature survey and are not used as load-bearing justification for Theorem 2; the central external inputs, especially Theorem 1/Greenbaum [9] and the error analysis in [13], are independent. The main weakness is that condition (76) is assumed of the particular k* whose existence Lemma 5 asserts, and the paper does not prove that such a k* exists for all SPD A, M satisfying the opening assumptions; the authors themselves note in Remark 1 and in the bfloat16 experiments that the analogous condition (99) can fail while the method still converges. This is a gap between the theorem's hypotheses and its advertised scope, and a correctness-risk concern, but it is not circularity: the theorem is conditional, and the unproved existence claim is not obtained by defining k* in terms of the desired bound. Likewise, the comparisons with Saad's split PCG are algorithmic and experimental, not a renaming of a prior result. Accordingly, no circular step of any of the enumerated kinds is present.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

The proof is largely self-contained and introduces no fitted parameters. The main external inputs are the standard rounding error model and the explicit, self-referential condition (76).

assumptions (4)
  • standard math Standard floating-point model: each arithmetic operation satisfies fl(a op b) = (a op b)(1+delta) with |delta| <= u, and preconditioner applications satisfy (14)-(16) with error parameters eps_s, eps_q, eps_z.
    Section 4.1 equations (8)-(18); the entire proof is built on the Higham-style rounding error model.
  • domain assumption M_L^{-1} A M_R^{-1} and M^{-1} = M_R^{-1} M_L^{-1} are symmetric positive definite, so the preconditioned problem is SPD.
    Algorithm 2 input and Section 3; needed for CG optimality and for the M-inner products used throughout.
  • standard math The function values f(bar x_k) are bounded below by f(x) for every computed iterate bar x_k, so a monotone sequence of f-values converges.
    Lemma 5, step 2 relies on f(x) <= f(bar x_k) <= f(x0) for every k after equation (7).
  • ad hoc to paper Condition (76) holds for the k* whose existence is asserted in Lemma 5.
    Theorem 2 and Lemma 6: without (76), the quadratic inequality for ||rhat_k*||_{M^{-1}} has no positive root and the bounds (82)-(84) do not follow; the paper does not prove that the k* from Lemma 5 satisfies it.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Forward and backward error bounds for a mixed precision preconditioned conjugate gradient algorithm." pith.science (2026). https://pith.science/paper/AD6DB4QL

@misc{pith2026251011379,
  author       = {Pith},
  title        = {Pith review of: Forward and backward error bounds for a mixed precision preconditioned conjugate gradient algorithm},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/AD6DB4QL}},
  note         = {Machine review of arXiv:2510.11379}
}
abstract

The preconditioned conjugate gradient (PCG) algorithm is one of the most popular algorithms for solving large-scale linear systems $Ax = b$, where $A$ is a symmetric positive definite matrix. Rather than computing residuals directly, it updates the residual vectors recursively. Current analyses of the conjugate gradient (CG) algorithm in finite precision typically assume that the norm of the recursively updated residual goes orders of magnitude below the machine precision, focusing mainly on bounding the residual gap thereafter. This work introduces a framework for the PCG algorithm and provides rigorous proofs that the relative backward and forward errors of the computed results of PCG can reach the levels $O(u)$ and $O(u)\kappa(A)^{1/2}$, respectively, after a sufficient number of iterations without relying on an assumption concerning the norm of the recursively updated residual, where $u$ represents the unit roundoff and $\kappa(A)$ is the condition number of $A$. Our PCG framework further shows that applying preconditioners in low precision does not compromise the accuracy of the final results, provided that reasonable conditions are satisfied. Moreover, this framework introduces a new split PCG variant that improves upon the classical split PCG algorithm when the left preconditioner is applied in low precision. Our theoretical results are illustrated through a set of numerical experiments.

Figures

Figures reproduced from arXiv: 2510.11379 by the authors.

Figure 1
Figure 1. Convergence curves for left PCG using the arithmetics listed in Table 1 for the appli [PITH_FULL_IMAGE:figures/full_fig_p029_1.png] view at source ↗
Figure 2
Figure 2. Comparison of convergence curves for split PCG variants using combinations of double [PITH_FULL_IMAGE:figures/full_fig_p030_2.png] view at source ↗
Figure 3
Figure 3. Achieved backward and forward errors by split PCG for (first row) and corresponding [PITH_FULL_IMAGE:figures/full_fig_p032_3.png] view at source ↗
Figures from the paper (1 more)
Figure 4
Figure 4. Figure 4: Convergence curves for split PCG using arithmetics listed in Table 1 excluding fp16. [PITH_FULL_IMAGE:figures/full_fig_p033_4.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

26 extracted references · 10 canonical work pages

  1. [1]

    Abdelfattah, H

    A. Abdelfattah, H. Anzt, E. G. Boman, E. Carson, T. Cojean, J. Dongarra, A. Fox, M. Gates, N. J. Higham, X. S. Li, J. Loe, P. Luszczek, S. Pranesh, S. Rajamanickam, T. Ribizel, B. F. Smith, K. Swirydowicz, S. Thomas, S. Tomov, Y. M. Tsai, and U. M. Yang. A survey of numerical linear algebra methods utilizing mixed-precision arithmetic.Int. J. High Perform...

    work page doi:10.1177/10943420211003313 2021

  2. [2]

    Bj¨ orck, T

    ˚A. Bj¨ orck, T. Elfving, and Z. Strakos. Stability of conjugate gradient and Lanczos methods for linear least squares problems.SIAM J. Matrix Anal. Appl., 19(3), 1998.doi:10.1137/ S089547989631202X

    1998

  3. [3]

    Carson and J

    E. Carson and J. Demmel. A residual replacement strategy for improving the maximum attainable accuracy ofs-step Krylov subspace methods.SIAM J. Matrix Anal. Appl., 35(1), 2014.doi:10.1137/120893057

    doi:10.1137/120893057 2014

  4. [4]

    Carson, J

    E. Carson, J. Liesen, and Z. Strakoˇ s. Towards understanding CG and GMRES through examples.Linear Algebra Appl., 692:241–291, 2024.doi:10.1016/j.laa.2024.04.003

    doi:10.1016/j.laa.2024.04.003 2024

  5. [5]

    E. C. Carson, M. Rozloˇ zn ´ ık, Z. Strakoˇ s, P. Tich´ y, and M. T ˚ uma. The numerical stabil- ity analysis of pipelined conjugate gradient methods: Historical context and methodology. SIAM J. Sci. Comput., (5), 2018.doi:10.1137/16M1103361

    doi:10.1137/16m1103361 2018

  6. [6]

    M. P. Connolly, N. J. Higham, and T. Mary. Stochastic rounding and its probabilistic backward error analysis.SIAM J. Sci. Comput., 43(1):A566–A585, 2021.doi:10.1137/ 20M1334796

    2021

  7. [7]

    El Arar, D

    E.-M. El Arar, D. Sohier, P. de Oliveira Castro, and E. Petit. Stochastic rounding variance and probabilistic bounds: A new approach.SIAM J. Sci. Comput., 45(5):C255–C275, 2023. doi:10.1137/22M1510819

    doi:10.1137/22m1510819 2023

  8. [8]

    Greenbaum

    A. Greenbaum. Behavior of slightly perturbed Lanczos and conjugate-gradient recurrences. Linear Algebra Appl., 113:7–63, 1989.doi:10.1016/0024-3795(89)90285-1

    work page doi:10.1016/0024-3795(89)90285-1 1989

Show all 26 references
  1. [9]

    Greenbaum

    A. Greenbaum. Estimating the attainable accuracy of recursively computed residual meth- ods.SIAM J. Matrix Anal. Appl., 18(3):535–551, 1997.doi:10.1137/S0895479895284944. 33

    1997 doi

  2. [10]

    Y. Guo, E. de Sturler, and T. Warburton. An adaptive mixed precision and dynamically scaled preconditioned conjugate gradient algorithm, 2025. URL:https://arxiv.org/abs/ 2505.04155

    2025 arXiv

  3. [11]

    M. H. Gutknecht and Z. Strakos. Accuracy of two three-term and three two-term recurrences for Krylov space solvers.SIAM J. Matrix Anal. Appl., 22(1):213–229, 2000.doi:10.1137/ S0895479897331862

  4. [12]

    M. R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. J. Res. Natl. Bur. Stand., 49(6):409–436, 1952.doi:10.6028/jres.049.044

    1952 doi

  5. [13]

    N. J. Higham.Accuracy and Stability of Numerical Algorithms. SIAM, Philadelphia, PA, USA, 2nd edition, 2002.doi:10.1137/1.9780898718027

    2002 doi

  6. [14]

    N. J. Higham and T. Mary. Mixed precision algorithms in numerical linear algebra.Acta Numer., 31:347–414, 2022.doi:10.1017/S0962492922000022

    2022 doi

  7. [15]

    I. C. F. Ipsen and H. Zhou. Probabilistic error analysis for inner products.SIAM J. Matrix Anal. Appl., 41(4):1726–1741, 2020.doi:10.1137/19M1270434

    2020 doi

  8. [16]

    Jir´ anek and M

    P. Jir´ anek and M. Rozloˇ zn ´ ık. Maximum attainable accuracy of inexact saddle point solvers. SIAM J. Matrix Anal. Appl., 29(4):1297–1321, 2008.doi:10.1137/060659727

    2008 doi

  9. [17]

    Liesen and Z

    J. Liesen and Z. Strakos.Krylov subspace methods: principles and analysis. Oxford Univer- sity Press, Oxford, UK, 2013.doi:10.1093/acprof:oso/9780199655410.003.0002

    2013

  10. [18]

    Meurant and Z

    G. Meurant and Z. Strakoˇ s. The Lanczos and conjugate gradient algorithms in finite preci- sion arithmetic.Acta Numer., 15:471–542, 2006.doi:10.1017/S096249290626001X

    2006 doi

  11. [19]

    C. C. Paige, M. Rozlozn ´ ık, and Z. Strakos. Modified Gram–Schmidt (MGS), least squares, and backward stability of MGS-GMRES.SIAM J. Matrix Anal. Appl., 28(1):264–284, 2006. doi:10.1137/050630416

    2006 doi

  12. [20]

    C. C. Paige and Z. Strakos. Residual and backward error bounds in minimum residual Krylov subspace methods.SIAM J. Sci. Comput., 23(6):1898–1923, 2002.doi:10.1137/ S1064827500381239

    1923

  13. [21]

    Saad.Iterative methods for sparse linear systems

    Y. Saad.Iterative methods for sparse linear systems. SIAM, Philadelphia, PA, USA, 2nd edition, 2003.doi:10.1137/1.9780898718003

    2003 doi

  14. [22]

    Scott and M

    J. Scott and M. T ˚ uma. Avoiding breakdown in incomplete factorizations in low precision arithmetic.ACM Trans. Math. Software, 50(2):1–25, 2024.doi:10.1145/3651155

    2024 doi

  15. [23]

    G. L. G. Sleijpen, H. A. Van der Vorst, and J. Modersitzki. Differences in the effects of rounding errors in Krylov solvers for symmetric indefinite linear systems.SIAM J. Matrix Anal. Appl., 22(3), 2001.doi:10.1137/S0895479897323087

    2001 doi

  16. [24]

    Strakoˇ s and P

    Z. Strakoˇ s and P. Tich` y. On error estimation in the conjugate gradient method and why it works in finite precision computations.Electron. Trans. Numer. Anal., 13:56–80, 2002

    2002

  17. [25]

    Strakoˇ s and P

    Z. Strakoˇ s and P. Tich` y. Error estimation in preconditioned conjugate gradients.BIT, 45:789–817, 2005.doi:10.1007/s10543-005-0032-1

    2005 doi

  18. [26]

    Vieubl´ e.Mixed precision iterative refinement for the solution of large sparse linear sys- tems

    B. Vieubl´ e.Mixed precision iterative refinement for the solution of large sparse linear sys- tems. Theses, INP Toulouse, Nov. 2022. URL:https://hal.science/tel-03975935. 34

    2022

Pith tools

Reviewed August 4, 2026 · model on record in the stance chip above.