Multi-Step Richardson Procedure

• Multi-step Richardson procedure is a numerical method that uses recursive extrapolation from multi-fidelity approximations to cancel leading error terms and enhance convergence.
• It is applied in systems such as ODE solvers, stochastic approximation, and quadrature by combining structured asymptotic evaluations on refined grids.
• Its implementations, including Richardson–Romberg extrapolation and residual correction, achieve higher order accuracy while preserving stability properties.

“Multi-step Richardson procedure” denotes a family of extrapolative or residual-correction constructions that use several approximations at different fidelities, discretization levels, or bias scales in order to cancel leading error terms and accelerate convergence. In some literatures the procedure is a repeated Richardson–Romberg tableau based on asymptotic expansions in $h$, $n^{-\alpha}$, or $\epsilon$; in others it is a long-step polynomial iteration, a residual-driven correction cycle, or a probabilistic multi-fidelity estimator. The common structure is the same: an observable admits a structured expansion, and several evaluations are combined so that lower-order terms vanish, leaving a higher-order residual (Fekete et al., 2023, Lemaire et al., 2014, Dick et al., 2017, Genans et al., 19 May 2026, Oates et al., 2024).

1. Core mathematical template

The basic prerequisite is an asymptotic error model. In one widely used form, a numerical approximation satisfies

$Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$

Richardson cancellation then combines values at refined steps so that one or more of the terms $c_1 h^p,c_2 h^{p+1},\dots$ disappear. For linear multistep ODE solvers, the first repeated extrapolants are written explicitly as

$$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$

$$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$

$$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$

with resulting orders $p+1$, $p+2$, and $n^{-\alpha}$0 (Fekete et al., 2023).

A more abstract formulation uses Vandermonde conditions. In stochastic debiasing with missing covariates, the $n^{-\alpha}$1-th order estimator is

$n^{-\alpha}$2

where $n^{-\alpha}$3 and the coefficients satisfy

$n^{-\alpha}$4

These constraints preserve the target term and annihilate successive homogeneous bias orders (Genans et al., 19 May 2026). The same linear-algebraic principle appears in multistep Richardson–Romberg for stochastic approximation, where block Vandermonde weights eliminate powers $n^{-\alpha}$5 from the implicit discretization bias (Frikha et al., 2014).

The concept also extends beyond scalar discretization error. For vector sequences $n^{-\alpha}$6 with asymptotic model

$n^{-\alpha}$7

a vectorized generalized Richardson process solves

$n^{-\alpha}$8

and sets

$n^{-\alpha}$9

In that setting the first $\epsilon$0 asymptotic components are eliminated exactly, and the remaining error admits its own asymptotic expansion (Sidi, 2016).

2. Global and repeated Richardson extrapolation for ODE solvers

In numerical ODEs, a particularly explicit multi-step Richardson procedure is the global Richardson extrapolation of linear multistep methods (LMMs). The underlying initial-value problem is

$\epsilon$1

and the base $\epsilon$2-step LMM has form

$\epsilon$3

The global/passive construction computes the same LMM independently on nested uniform grids with steps $\epsilon$4, $\epsilon$5, $\epsilon$6, and so on, then combines the whole-grid solutions only after they are available. For one level, the extrapolated approximation is

$\epsilon$7

and, under the asymptotic global error expansion

$\epsilon$8

the order increases from $\epsilon$9 to at least $Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$0 (Fekete et al., 2022).

The repeated version applies Richardson extrapolation $Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$1 times across the dyadic grid chain $Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$2. Under assumptions $Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$3-$Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$4, the resulting sequence has order $Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$5 (Fekete et al., 2023). This is a genuine multi-step Richardson procedure in the classical sense: each additional extrapolation stage cancels one more power of the global error expansion. The construction is global rather than local; the paper explicitly distinguishes this from local/active Richardson extrapolation, which would feed extrapolated values back into the multistep history and is left as future work (Fekete et al., 2022).

A major theoretical feature is the preservation of stability structure. For one-step GRE,

$Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$6

and for repeated GRE,

$Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$7

If the base stability region is convex, equality follows (Fekete et al., 2022, Fekete et al., 2023). For BDF methods this yields the same $Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$8-stability angle as the underlying method, while increasing the order. One concrete consequence is that $Y(h)=Y^\star + c_1 h^p + c_2 h^{p+1}+\cdots.$9-GRE is a third-order $c_1 h^p,c_2 h^{p+1},\dots$0-stable method, and $c_1 h^p,c_2 h^{p+1},\dots$1-2GRE is a fourth-order $c_1 h^p,c_2 h^{p+1},\dots$2-stable method; there is no contradiction with the Dahlquist barrier because the extrapolated method is not itself a standard fixed-coefficient LMM (Fekete et al., 2022, Fekete et al., 2023).

The practical attraction is that existing LMM codes can be used without modification. For one extrapolation level the cost is roughly three times that of one LMM run; for $c_1 h^p,c_2 h^{p+1},\dots$3 repeated levels it is approximately $c_1 h^p,c_2 h^{p+1},\dots$4 times the base cost (Fekete et al., 2022, Fekete et al., 2023). This suggests a characteristic Richardson trade-off: extra solves in exchange for higher order and, in favorable cases, preserved stability geometry.

3. Bias cancellation in stochastic approximation and learning

A second major branch uses multi-step Richardson procedures for bias reduction rather than deterministic truncation error. In stochastic approximation, the target is the zero $c_1 h^p,c_2 h^{p+1},\dots$5 of

$c_1 h^p,c_2 h^{p+1},\dots$6

but only an approximation $c_1 h^p,c_2 h^{p+1},\dots$7 is simulable, leading to $c_1 h^p,c_2 h^{p+1},\dots$8 and its zero $c_1 h^p,c_2 h^{p+1},\dots$9. Under the expansion

$$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$0

the multistep Richardson–Romberg estimator combines $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$1 stochastic approximation runs at levels $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$2: $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$3 The weights solve a block Vandermonde system and produce residual bias of order $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$4 instead of $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$5 (Frikha et al., 2014). In that setting the optimized cost exponent improves from $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$6 for crude SA to $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$7 for the Richardson–Romberg version (Frikha et al., 2014).

In learning with missing covariates, the parameter of interest is not a root of a weak discretization but a population gradient. The paper proves that the imputation-induced bias has leading order $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$8 in the missingness vector $$r_n^{[1]}(h)= \frac{2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)}{2^p-1},$$9, and proposes to deliberately add missingness, generating a further-thinned observation at scale $$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$0. The one-step corrected gradient

$$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$1

cancels the first-order bias, reducing it from $$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$2 to $$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$3 (Genans et al., 19 May 2026). Under independent missing indicators, the bias is an exact multilinear polynomial in $$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$4, and the higher-order estimator

$$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$5

has bias $$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$6; if $$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$7, the population gradient bias is canceled exactly (Genans et al., 19 May 2026). This is a genuine multi-step Richardson construction, but it depends critically on the polynomial structure induced by independence of missing indicators.

A related but narrower case is constant-stepsize linear stochastic approximation under Markovian noise. There the averaged iterate has asymptotic bias

$$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$8

and the paper applies the two-level Richardson–Romberg correction

$$r_n^{[2]}(h)= \frac{2^{2p+1}\,y_{4n}\!\left(\frac h4\right)-3\cdot 2^p\,y_{2n}\!\left(\frac h2\right)+y_n(h)} {(2^p-1)(2^{p+1}-1)},$$9

to eliminate the linear term (Levin et al., 7 Aug 2025). That work explicitly analyzes only the first extrapolation step rather than a full multi-step tableau, although its $$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$0 bias decomposition suggests a higher-order continuation (Levin et al., 7 Aug 2025). A common misconception is therefore inaccurate: not every paper invoking Richardson–Romberg in stochastic approximation actually develops a repeated multistep hierarchy.

4. Quadrature, integral equations, and discretization-parameter extrapolation

Multi-step Richardson procedures also appear when the extrapolated variable is neither a timestep nor a Monte Carlo bias scale, but another discretization or regularization parameter. In quasi-Monte Carlo integration over weighted Sobolev spaces of dominating mixed smoothness $$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$1, extrapolated polynomial lattice rules are built from a chain of classical polynomial lattice rules with sizes

$$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$2

The recursive Richardson step is

$$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$3

and after $$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$4 stages the first $$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$5 powers in the error expansion are canceled (Dick et al., 2017). The resulting extrapolated polynomial lattice rule achieves worst-case error

$$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$6

equivalently almost $$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$7, while retaining the classical polynomial-lattice structure needed for the fast QMC matrix-vector method (Dick et al., 2017). The total construction cost is $$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$8, improving on the $$r_n^{[3]}(h)= \frac{2^{3p+3}\,y_{8n}\!\left(\frac h8\right)-7\cdot 2^{2p+1}\,y_{4n}\!\left(\frac h4\right)+7\cdot 2^p\,y_{2n}\!\left(\frac h2\right)-y_n(h)} {(2^p-1)(2^{p+1}-1)(2^{p+2}-1)},$$9 cost quoted for interlaced polynomial lattice rules (Dick et al., 2017).

For the method of regularized stokeslets, the extrapolated parameter is the regularization radius $p+1$0. The paper uses the three-value sequence

$p+1$1

and the matrix formula

$p+1$2

which cancels the $p+1$3 and $p+1$4 regularization terms and leaves

$p+1$5

(Gallagher et al., 2021). The paper interprets this as repeated-in-spirit Richardson extrapolation in $p+1$6, with concrete gains: for the unit sphere, the minimum relative error improves from $p+1$7 for the raw Nyström method to $p+1$8 for the extrapolated version, and $p+1$9 error is achieved in $p+2$0 seconds walltime; for the prolate spheroid, the small-$p+2$1 plateau at $p+2$2 drops from $p+2$3 to $p+2$4 (Gallagher et al., 2021).

A simpler pedagogical instance arises when a double integral is transformed into a second-order IVP,

$p+2$5

then solved by Euler’s method. Since the raw approximation has expansion

$p+2$6

successively refined meshes $p+2$7 can be combined into $p+2$8, canceling one more power of $p+2$9 at each stage (Prentice, 2023). In the reported example, the final value at $n^{-\alpha}$00 is within $n^{-\alpha}$01 of the exact integral (Prentice, 2023). This case makes explicit that a multi-step Richardson procedure may be attached to a very low-order base scheme and still produce near-machine-precision results, provided the expansion is regular enough.

5. Residual-correction, long-step, and iterative-linear-algebra interpretations

In linear algebra, the phrase “multi-step Richardson procedure” becomes more heterogeneous. Some papers use it for genuine long-step polynomial methods, others for repeated one-step Richardson updates, and others only by analogy. For computing $n^{-\alpha}$02, the matrix exponential can be reinterpreted as the ODE

$n^{-\alpha}$03

Given an approximation $n^{-\alpha}$04, the residual is

$n^{-\alpha}$05

and the exact error solves

$n^{-\alpha}$06

The Richardson update is

$n^{-\alpha}$07

where $n^{-\alpha}$08 is an approximate solution of the residual-driven correction equation (Botchev, 2011). Repeated residual corrections yield a multi-stage process, and in the Krylov version each stage becomes a restart mechanism based on the separable residual form $n^{-\alpha}$09 (Botchev, 2011). This is multi-step in the sense of repeated correction stages, not in the sense of a fixed multistep recurrence.

By contrast, “Richardson($n^{-\alpha}$10)” for large linear systems is a genuine long-step formulation. Starting from the classical iteration

$n^{-\alpha}$11

the paper bundles $n^{-\alpha}$12 consecutive Richardson steps into one outer step with weight schedule $n^{-\alpha}$13: $n^{-\alpha}$14 This makes the outer iteration a degree-$n^{-\alpha}$15 polynomial method (Wang et al., 2024). Momentum and preconditioning are then added, producing MOM-Richardson($n^{-\alpha}$16), NAG-Richardson($n^{-\alpha}$17), and NAGex-Richardson($n^{-\alpha}$18) (Wang et al., 2024). Numerically, increasing $n^{-\alpha}$19 gives strong acceleration: at $n^{-\alpha}$20 in the anisotropic diffusion test, Richardson(1)-NS requires $n^{-\alpha}$21 iterations, Richardson(3)-NS $n^{-\alpha}$22, Richardson(7)-NS $n^{-\alpha}$23, and Richardson(15)-NS $n^{-\alpha}$24 (Wang et al., 2024).

A different residual-refinement interpretation appears in fixed-point inverse solvers. There the base method is standard Richardson iteration on the normal equations, but the full solver performs outer residual updates

$n^{-\alpha}$25

where each correction $n^{-\alpha}$26 is itself computed by an inner Richardson solve (Zhu et al., 2021). The one-stage asymptotic error floor is $n^{-\alpha}$27, and after $n^{-\alpha}$28 residual-correction stages the paper proves

$n^{-\alpha}$29

(Zhu et al., 2021). This is again a multi-step Richardson procedure in the residual-correction sense, designed to break the fixed-point precision barrier rather than to alter the spectral polynomial explicitly.

Not every Richardson paper in linear algebra belongs to this category. For ill-conditioned least-squares systems, one paper explicitly states that it does not introduce a genuine multi-step Richardson recurrence and instead studies repeated one-step preconditioned updates

$n^{-\alpha}$30

with residual-based stopping (Stotsky, 2023). That distinction is important: repeated application of one-step Richardson is not automatically the same object as a multi-step Richardson procedure with a designed hierarchy of bias cancellations or polynomial filters.

6. Probabilistic, geometric, and terminological generalizations

Recent work generalizes Richardson logic in two opposite directions: toward highly abstract extrapolation theory and toward conceptually distinct uses of the Richardson name. In probabilistic Richardson Extrapolation, the numerical output at fidelity $n^{-\alpha}$31 is modeled by a Gaussian process with covariance

$n^{-\alpha}$32

The continuum estimate is the posterior mean at $n^{-\alpha}$33,

$n^{-\alpha}$34

which is a weighted linear combination of all available fidelities (Oates et al., 2024). In dimension $n^{-\alpha}$35, if $n^{-\alpha}$36 reproduces polynomials and $n^{-\alpha}$37, the posterior mean recovers polynomial extrapolation to $n^{-\alpha}$38, hence classical Richardson (Oates et al., 2024). The framework also converts fidelity selection into an optimization problem over

$n^{-\alpha}$39

something absent from classical fixed-grid tableaux (Oates et al., 2024). This suggests a broad reinterpretation: a multi-step Richardson procedure can be viewed as a particular optimal linear estimator over a fidelity-dependent function class.

At the abstract end, the vectorized generalized Richardson process already mentioned shows that extrapolation need not be tied to powers of a scalar mesh width. The transformed sequence $n^{-\alpha}$40 cancels the first $n^{-\alpha}$41 model terms in a general asymptotic scale $n^{-\alpha}$42, and the resulting error has its own full asymptotic expansion (Sidi, 2016). This places many domain-specific multistep Richardson procedures inside a wider asymptotic-elimination framework.

The term can also refer to a conceptually distinct construction outside numerical extrapolation. In stochastic geometric mechanics, a “Richardson triple” is the composition

$n^{-\alpha}$43

with fast, intermediate, and slow flow maps. The procedure is sequential: first homogenize the fast/intermediate composition $n^{-\alpha}$44 into a stochastic flow, then incorporate the slow component by moving into the large-scale non-inertial frame via $n^{-\alpha}$45 (Holm, 2017). This is a two-stage Richardson reduction of a three-scale flow, not a Richardson extrapolation. Its presence in the literature illustrates that “Multi-step Richardson Procedure” is not a single universal term and must be interpreted from context.

Across these variants, several constraints recur. A structured expansion is essential; without the clean leading-order form, cancellation can weaken or fail. Startup accuracy matters in multistep ODE extrapolation; nonconvex stability regions can shrink the extrapolated stability set; dependence among missing indicators blocks the multilinear polynomial bias structure needed for higher-order Richardson-SGD; and higher-order combinations often increase computational cost or variance (Fekete et al., 2022, Fekete et al., 2023, Genans et al., 19 May 2026). A plausible implication is that the most robust use of a multi-step Richardson procedure is not “apply more levels whenever possible,” but “apply as many levels as the available asymptotic structure, coupling mechanism, and cost model can genuinely support.”