---
title: Multi-Step Richardson Procedure
url: https://www.emergentmind.com/topics/multi-step-richardson-procedure
type: topic
---

# Multi-Step Richardson Procedure

“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 [2307.01345] [1401.1177] [1707.03989] [2605.19641] [2401.07562].

## 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 \(p+3\) [2307.01345].

A more abstract formulation uses Vandermonde conditions. In stochastic debiasing with missing covariates, the \(k\)-th order estimator is
\[
\hat g^{[k]}(w):=\sum_{\ell=0}^{k}\alpha_\ell\,\hat g^{(C_\ell p)}(w),
\]
where \(1=C_0<C_1<\cdots<C_k\) and the coefficients satisfy
\[
\sum_{\ell=0}^k\alpha_\ell=1,\qquad \sum_{\ell=0}^k\alpha_\ell C_\ell^m=0,\quad m=1,\dots,k.
\]
These constraints preserve the target term and annihilate successive homogeneous bias orders [2605.19641]. The same linear-algebraic principle appears in multistep Richardson–Romberg for stochastic approximation, where block Vandermonde weights eliminate powers \(n^{-\alpha},\dots,n^{-\alpha(R-1)}\) from the implicit discretization bias [1409.4748].

The concept also extends beyond scalar discretization error. For vector sequences \(\{\mathbf{x}_m\}\) with asymptotic model
\[
\mathbf{x}_m\sim \mathbf{s} + \sum_{i=1}^\infty \alpha_i \mathbf{g}_i(m),
\]
a vectorized generalized Richardson process solves
\[
\sum_{i=1}^{k}\langle \mathbf{y},\Delta \mathbf{g}_i(m)\rangle \widetilde{\alpha}_i=\langle \mathbf{y},\Delta \mathbf{x}_m\rangle,\qquad n\le m\le n+k-1,
\]
and sets
\[
\mathbf{s}_{n,k}=\mathbf{x}_n-\sum_{i=1}^k \widetilde{\alpha}_i\mathbf{g}_i(n).
\]
In that setting the first \(k\) asymptotic components are eliminated exactly, and the remaining error admits its own asymptotic expansion [1605.02630].

## 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
\[
y'(t)=f(t,y(t)), \qquad y(t_0)=y_0,
\]
and the base \(k\)-step LMM has form
\[
\sum_{j=0}^k \alpha_j y_{n+j}=\sum_{j=0}^k h \beta_j f_{n+j}.
\]
The global/passive construction computes the same LMM independently on nested uniform grids with steps \(h\), \(h/2\), \(h/4\), and so on, then combines the whole-grid solutions only after they are available. For one level, the extrapolated approximation is
\[
r_n(h):=\frac{2^p}{2^p-1}\cdot y_{2n}\left(\frac{h}{2}\right)-\frac{1}{2^p-1}\cdot y_n(h),
\]
and, under the asymptotic global error expansion
\[
y_n(h)-y(t_n)=C_{k,p} h^p \mathbf e(t_n)+{\cal O}(h^{p+1}),
\]
the order increases from \(p\) to at least \(p+1\) [2206.10220].

The repeated version applies Richardson extrapolation \(\ell\) times across the dyadic grid chain \(h,h/2,\dots,h/2^\ell\). Under assumptions \(\boxed{a1}\)-\(\boxed{a4}\), the resulting sequence has order \(p+\ell\) [2307.01345]. 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 [2206.10220].

A major theoretical feature is the preservation of stability structure. For one-step GRE,
\[
{\mathcal S}_\text{LMM}\cap 2{\mathcal S}_\text{LMM}\subseteq {\mathcal S}_\text{GRE}\subseteq {\mathcal S}_\text{LMM},
\]
and for repeated GRE,
\[
\bigcap_{j=0}^{\ell} \bigl(2^j\mathcal S_{\mathrm{LMM}}\bigr) \subseteq \mathcal S^{[\ell]}_{\mathrm{RGRE}} \subseteq \mathcal S_{\mathrm{LMM}}.
\]
If the base stability region is convex, equality follows [2206.10220] [2307.01345]. For BDF methods this yields the same \(A(\alpha)\)-stability angle as the underlying method, while increasing the order. One concrete consequence is that \(\mathrm{BDF}2\)-GRE is a third-order \(A\)-stable method, and \(\mathrm{BDF}2\)-2GRE is a fourth-order \(A\)-stable method; there is no contradiction with the Dahlquist barrier because the extrapolated method is not itself a standard fixed-coefficient LMM [2206.10220] [2307.01345].

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 \(\ell\) repeated levels it is approximately \(2^{\ell+1}-1\) times the base cost [2206.10220] [2307.01345]. 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 \(\theta^\*\) of
\[
h(\theta)=\mathbb E[H(\theta,U)],
\]
but only an approximation \(U^n\) is simulable, leading to \(h^n(\theta)=\mathbb E[H(\theta,U^n)]\) and its zero \(\theta^{*,n}\). Under the expansion
\[
\theta^{*,n}-\theta^\* = \frac{C_1}{n^\alpha}+\cdots+\frac{C_p}{n^{\alpha p}}+o(n^{-\alpha p}),
\]
the multistep Richardson–Romberg estimator combines \(R\) stochastic approximation runs at levels \(n,2n,\dots,Rn\):
\[
\Theta_M^{(R)}(n):=\sum_{r=1}^R \mathbf w_r \theta_M^{rn}.
\]
The weights solve a block Vandermonde system and produce residual bias of order \(n^{-\alpha R}\) instead of \(n^{-\alpha}\) [1409.4748]. In that setting the optimized cost exponent improves from \(\frac{2}{\beta}+\frac{1}{\alpha}\) for crude SA to \(\frac{2}{\beta}+\frac{1}{\alpha R}\) for the Richardson–Romberg version [1409.4748].

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 \(O(\|p\|)\) in the missingness vector \(p\), and proposes to deliberately add missingness, generating a further-thinned observation at scale \(Cp\). The one-step corrected gradient
\[
\hat g_C^{\mathrm R}(w):=\frac{C\,\hat g^{(p)}(w)-\hat g^{(Cp)}(w)}{C-1}
\]
cancels the first-order bias, reducing it from \(O(\|p\|)\) to \(O(\|p\|^2)\) [2605.19641]. Under independent missing indicators, the bias is an exact multilinear polynomial in \(p\), and the higher-order estimator
\[
\hat g^{[k]}(w):=\sum_{\ell=0}^{k}\alpha_\ell\,\hat g^{(C_\ell p)}(w)
\]
has bias \(O(\|p\|^{k+1})\); if \(k=d_{\mathrm{miss}}\), the population gradient bias is canceled exactly [2605.19641]. 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
\[
\Pi_\alpha(\theta_0)-\theta^\*=\alpha\Delta+O(\alpha^{3/2}),
\]
and the paper applies the two-level Richardson–Romberg correction
\[
\bar\theta_n^{\alpha,\mathrm{RR}}=2\bar\theta_n^\alpha-\bar\theta_n^{2\alpha}
\]
to eliminate the linear term [2508.05570]. That work explicitly analyzes only the first extrapolation step rather than a full multi-step tableau, although its \(J_n^{(1,\alpha)},J_n^{(2,\alpha)},\dots\) bias decomposition suggests a higher-order continuation [2508.05570]. 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 \(\alpha\ge 2\), extrapolated polynomial lattice rules are built from a chain of classical polynomial lattice rules with sizes
\[
b^{m-\alpha+1},\, b^{m-\alpha+2},\,\ldots,\, b^m.
\]
The recursive Richardson step is
\[
I_n^{(\tau+1)}=\frac{b^\tau I_n^{(\tau)}-I_{n-1}^{(\tau)}}{b^\tau-1},\qquad 1\le \tau<\alpha,
\]
and after \(\alpha-1\) stages the first \(\alpha-1\) powers in the error expansion are canceled [1707.03989]. The resulting extrapolated polynomial lattice rule achieves worst-case error
\[
O(N^{-1/\lambda}) \qquad \text{for any } \frac1\alpha<\lambda\le 1,
\]
equivalently almost \(O(N^{-\alpha+\varepsilon})\), while retaining the classical polynomial-lattice structure needed for the fast QMC matrix-vector method [1707.03989]. The total construction cost is \(O((s+\alpha)N\log N)\), improving on the \(O(s\alpha N\log N)\) cost quoted for interlaced polynomial lattice rules [1707.03989].

For the method of regularized stokeslets, the extrapolated parameter is the regularization radius \(\epsilon\). The paper uses the three-value sequence
\[
(\epsilon_1,\epsilon_2,\epsilon_3)=\left(\epsilon,\sqrt{2}\epsilon,2\epsilon\right)
\]
and the matrix formula
\[
\widetilde M(\epsilon_1,\epsilon_2,\epsilon_3;h):=
\begin{pmatrix}1&0&0\end{pmatrix}
B^{-1}
\begin{pmatrix}
M(\epsilon_1,h)\\
M(\epsilon_2,h)\\
M(\epsilon_3,h)
\end{pmatrix},
\]
which cancels the \(O(\epsilon)\) and \(O(\epsilon^2)\) regularization terms and leaves
\[
E_d(h;\epsilon_0)+O(\epsilon_1^3+\epsilon_2^3+\epsilon_3^3)
\]
[2101.09286]. The paper interprets this as repeated-in-spirit Richardson extrapolation in \(\epsilon\), with concrete gains: for the unit sphere, the minimum relative error improves from \(0.6\%\) for the raw Nyström method to \(0.05\%\) for the extrapolated version, and \(0.1\%\) error is achieved in \(250\) seconds walltime; for the prolate spheroid, the small-\(h\) plateau at \(\epsilon=0.2\) drops from \(8.7\%\) to \(0.059\%\) [2101.09286].

A simpler pedagogical instance arises when a double integral is transformed into a second-order IVP,
\[
C''(x)=g(x),
\]
then solved by Euler’s method. Since the raw approximation has expansion
\[
N(h)=M+K_1(x)h+K_2(x)h^2+K_3(x)h^3+\cdots,
\]
successively refined meshes \(h,h/2,h/4,\dots\) can be combined into \(M_2,M_3,M_4,M_5\), canceling one more power of \(h\) at each stage [2305.07777]. In the reported example, the final value at \(x=5\) is within \(8.3\times 10^{-13}\) of the exact integral [2305.07777]. 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 \(e^{-tA}v\), the matrix exponential can be reinterpreted as the ODE
\[
y'(t)=-Ay(t),\qquad y(0)=v.
\]
Given an approximation \(y_k(t)\), the residual is
\[
r_k(t)=-Ay_k(t)-y_k'(t),
\]
and the exact error solves
\[
\epsilon_k'(t)=-A\epsilon_k(t)+r_k(t),\qquad \epsilon_k(0)=0.
\]
The Richardson update is
\[
y_{k+1}(t)=y_k(t)+\tilde\epsilon_k(t),
\]
where \(\tilde\epsilon_k\) is an approximate solution of the residual-driven correction equation [1112.5670]. 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 \(r_k(t)=\psi_k(t)w_k\) [1112.5670]. This is multi-step in the sense of repeated correction stages, not in the sense of a fixed multistep recurrence.

By contrast, “Richardson(\(m\))” for large linear systems is a genuine long-step formulation. Starting from the classical iteration
\[
\bm u_k = \bm u_{k-1} + \omega(\bm f-\bm A\bm u_{k-1}),
\]
the paper bundles \(m\) consecutive Richardson steps into one outer step with weight schedule \(\bm\omega=(\omega_1,\dots,\omega_m)\):
\[
\bm T_m(\bm A,\bm \omega) = (\bm I-\omega_m\bm A)\cdots(\bm I-\omega_1\bm A).
\]
This makes the outer iteration a degree-\(m\) polynomial method [2412.08076]. Momentum and preconditioning are then added, producing MOM-Richardson(\(m\)), NAG-Richardson(\(m\)), and NAGex-Richardson(\(m\)) [2412.08076]. Numerically, increasing \(m\) gives strong acceleration: at \(\epsilon=10^{-6}\) in the anisotropic diffusion test, Richardson(1)-NS requires \(27257\) iterations, Richardson(3)-NS \(3315\), Richardson(7)-NS \(685\), and Richardson(15)-NS \(184\) [2412.08076].

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
\[
r^{(l)}=y-Ax^{(l)},\qquad x^{(l+1)}=x^{(l)}+\delta x^{(l+1)},
\]
where each correction \(\delta x^{(l+1)}\) is itself computed by an inner Richardson solve [2105.02106]. The one-stage asymptotic error floor is \(\theta\), and after \(M\) residual-correction stages the paper proves
\[
\theta^{(M)}\le \theta^M
\]
[2105.02106]. 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
\[
V_i = V_{i-1} - G_i(A_kV_{i-1}-b_k),
\]
with residual-based stopping [2302.08868]. 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 \(\mathbf{x}\) is modeled by a Gaussian process with covariance
\[
k(\mathbf{x},\mathbf{x}')=\sigma^2\left[k_0^2+b(\mathbf{x})b(\mathbf{x}')k_e(\mathbf{x},\mathbf{x}')\right].
\]
The continuum estimate is the posterior mean at \(\mathbf{0}\),
\[
m_n[f](\mathbf{0})=\frac{\mathbf{1}^\top \mathbf{K}_b^{-1}f(X_n)}{\mathbf{1}^\top \mathbf{K}_b^{-1}\mathbf{1}},
\]
which is a weighted linear combination of all available fidelities [2401.07562]. In dimension \(d=1\), if \(k_e\) reproduces polynomials and \(k_0^2\to\infty\), the posterior mean recovers polynomial extrapolation to \(0\), hence classical Richardson [2401.07562]. The framework also converts fidelity selection into an optimization problem over
\[
\argmax_{X\subset \mathcal D}\ \mathbf{1}^\top \mathbf{K}_b^{-1}\mathbf{1}
\quad \text{s.t.}\quad \sum_{\mathbf{x}\in X} c(\mathbf{x})\le C,
\]
something absent from classical fixed-grid tableaux [2401.07562]. 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 \(\mathbf{s}_{n,k}\) cancels the first \(k\) model terms in a general asymptotic scale \(\{\mathbf{g}_i(m)\}\), and the resulting error has its own full asymptotic expansion [1605.02630]. 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
\[
g_{t/\varepsilon,t,\varepsilon t}=h_{t/\varepsilon}k_t l_{\varepsilon t},
\]
with fast, intermediate, and slow flow maps. The procedure is sequential: first homogenize the fast/intermediate composition \(h_{t/\varepsilon}k_t\) into a stochastic flow, then incorporate the slow component by moving into the large-scale non-inertial frame via \((\mathrm{d}g)l_{\varepsilon t}\) [1708.04183]. 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 [2206.10220] [2307.01345] [2605.19641]. 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.”

Source: https://www.emergentmind.com/topics/multi-step-richardson-procedure