Numerical Analysis

Motivation

Sampling from a Gibbs distribution $\pi_U\propto e^{-U}\,d\lambda_d$ on $\mathbb R^d$ is a basic task in Bayesian statistics, statistical physics and machine learning. A standard approach simulates the Langevin diffusion $dX_t=-\nabla U(X_t)\,dt+\sqrt2\,dB_t$, whose invariant distribution is $\pi_U$, through its Euler scheme, and averages a test function $f$ along one long trajectory. The discretization introduces a bias of the order of the step size, so a single-step-size average with accuracy $\varepsilon$ costs more than $\varepsilon^{-2}$ steps. Dalalyan (arXiv:1412.7392) and Durmus–Moulines (arXiv:1605.01559) quantified this for strongly log-concave targets; for occupation measures of the Euler scheme the best previous bound was of order $d\varepsilon^{-4}$ (Durmus–Moulines, Bernoulli 2019, denoted [DM19] in the paper).

Egéa and Panloup (arXiv:2109.07753, Math. Oper. Res., 2025) combine Giles's multilevel Monte Carlo idea (Oper. Res. 2008) with pathwise averages. Each level is one long trajectory, not many independent samples. This mission formalizes their result for a uniformly strongly contractive, twice differentiable drift. In that case the multilevel occupation measure reaches accuracy $\varepsilon$ at a cost of order $d\varepsilon^{-2}$, the order of an unbiased Monte Carlo method.

Setting

Let $b:\mathbb R^d\to\mathbb R^d$ be $L$-Lipschitz and satisfy Assumption $(\mathbf C_\alpha)$: $\langle b(y)-b(x),y-x\rangle\le-\alpha|x-y|^2$ for all $x,y$, with $\alpha\in(0,L\wedge1]$. Let $\sigma>0$. The additive diffusion is $X^x_t=x+\int_0^tb(X^x_s)\,ds+\sigma B_t$, driven by a standard $d$-dimensional Brownian motion $B$. Let $\pi$ be its invariant distribution. Its Euler scheme with step $\gamma$ is $\bar X^{\gamma,x}_0=x$, $\bar X^{\gamma,x}_{(n+1)\gamma}=\bar X^{\gamma,x}_{n\gamma}+\gamma b(\bar X^{\gamma,x}_{n\gamma})+\sigma(B_{(n+1)\gamma}-B_{n\gamma})$, with invariant distribution $\pi^\gamma$. Write $\underline s_\eta=\eta\lfloor s/\eta\rfloor$ and $[f]_1$ for the Lipschitz constant of $f$.

Given a first step $\gamma_0$, steps $\gamma_r=\gamma_02^{-r}$, horizons $T_0>\dots>T_R$, a warm start $\tau$, and independent Brownian motions $B^{(0)},\dots,B^{(R)}$, the multilevel occupation measure (7) is

$$\mathcal Y(f)=\frac1{T_0-\tau}\int_\tau^{T_0}f(\bar X^{\gamma_0}_{\underline s_{\gamma_0}})ds+\sum_{r=1}^R\frac1{T_r-\tau}\int_\tau^{T_r}\big(f(\bar X^{\gamma_r}_{\underline s_{\gamma_{r-1}}})-f(\bar X^{\gamma_{r-1}}_{\underline s_{\gamma_{r-1}}})\big)ds,$$

where both schemes of level $r$ start at $x_0$ and are driven by $B^{(r)}$. Its cost is $\mathcal C(\mathcal Y)=T_0/\gamma_0+\sum_{r=1}^R(T_r/\gamma_r+T_r/\gamma_{r-1})$.

When $b$ is $\mathcal C^2$ with bounded derivatives, set $\|\nabla b\|_{2,\infty}=\sup_x\|\nabla b(x)\|_F$ (Frobenius norm of the Jacobian) and $\|\Delta b\|^2_{2,\infty}=\sup_x\sum_i|\Delta b_i(x)|^2$. The constant (19) is

$$\Upsilon_2^2=\max\Big(1,\ \frac{|b(x_0)|^2}{\alpha^2}+\frac{\sigma\sqrt\alpha}{L^3}\|\nabla b\|_{2,\infty}|b(x_0)|+\frac{\sigma^4}{L^4}\|\Delta b\|^2_{2,\infty}+\frac{\sigma^2\alpha\sqrt d}{L^3}\|\nabla b\|_{2,\infty}+\frac{\sigma^2d}{\alpha}\Big).$$

Formalization targets

Goal: Theorem 2.3

Take $\gamma_0\in(0,\frac\alpha{2L^2}\wedge\frac12]$, $R_\varepsilon=\lceil\log_2(\Upsilon_2/\varepsilon)\rceil$, $T_r=\Upsilon_2^2\log(\gamma_0^{-1})\alpha^{-1}\varepsilon^{-2}2^{-3r/2}$ and $\tau=\alpha^{-1}|\log\varepsilon|+(2\alpha)^{-1}\log\Upsilon_2$. There is a universal constant $\mathfrak C_1$ such that for every $\varepsilon\in(0,1)$ with $\sqrt2\log(\Upsilon_2\varepsilon^{-2})\le\log(\gamma_0^{-1})(\Upsilon_2/\varepsilon)^{1/2}$,

$$\|\mathcal Y(f)-\pi(f)\|_2\le\mathfrak C_1[f]_1\varepsilon,\qquad\mathcal C(\mathcal Y)\le\frac{\frac12+\sqrt2}{\sqrt2-1}\,\frac{\Upsilon_2^2\log(\gamma_0^{-1})}{\gamma_0\alpha}\,\varepsilon^{-2}.$$

The constant $\mathfrak C_1$ is the same for all $d$, $b$, $\sigma$, $f$, $x_0$, $\gamma_0$ and $\varepsilon$.

Companion: (22)

Under $\alpha/L^2\le1$, $\sigma^2\alpha^{-1}d\ge1$ and dimension-free bounds on $|b(x_0)|$, $\|\nabla b\|_{2,\infty}$ and $\|\Delta b\|_{2,\infty}$, the choice $\gamma_0=\alpha/(2L^2)$, $\tilde\Upsilon_2^2=\sigma^2\alpha^{-1}d$ gives $\mathcal C(\mathcal Y)\le\frac{2\sqrt2+1}{\sqrt2-1}\frac{\sigma^2L^2}{\alpha^3}\log\frac{2L^2}{\alpha}\,d\varepsilon^{-2}$.

Milestones

The milestones follow the paper's proof in order:

• Lemma 5.1 (i), (iii): stability of the scheme and the assumptions $(\mathbf H_1)$, $(\mathbf H_4)$.
• Lemma 5.2 (ii): a one-step error of order $\gamma^3$.
• Proposition 5.1 (ii), (iii): a uniform-in-time strong error $\beta_2(x)\gamma^2$ and $(\mathbf H_2)$ with $\mathfrak a=2$.
• Proposition 5.2: $\mathcal W_1(\pi,\pi^\gamma)\le\sqrt{\beta_2(x^\star)}\gamma$.
• Proposition 2.2: $(\mathbf H_2)$–$(\mathbf H_4)$ with constants controlled by $\Upsilon_2$.
• Remark 4.1: the paper's general complexity theorem, with constants known up to universal factors.

Significance

The result. Theorem 2.3 gives an explicit algorithm whose cost to reach $L^2$ accuracy $\varepsilon$ on $\pi(f)$ is of order $\varepsilon^{-2}$, uniformly over Lipschitz $f$. Under (22) the cost is $O(d\varepsilon^{-2})$. Apart from constants depending on $\alpha$ and $L$, this is the complexity of a Monte Carlo method without bias, an improvement over the $d\varepsilon^{-4}$ known for single-level occupation measures. Corollary 2.2 of the paper applies it to Gibbs distributions with a $\mathcal C^3$ potential.

Formalizing it. The theorem is proved in the paper. A machine-checked version needs a formal multilevel estimator over infinite horizons, uniform-in-time strong error bounds for the Euler scheme, and an application of Itô's formula to $b(y+\sigma B_t)$, none of which exists in Lean. Formalizing also checks the paper's constants. For instance, the printed exponential rate in the $\mathcal W_2$ contraction behind $(\mathbf H_1)$ is $\alpha$, while its proof gives $\alpha/2$, and the milestone states $\alpha/2$. No part of this result has been machine-checked.

Difficulty

Two steps resist a direct approach. First, the order-one strong error (Lemma 5.2(ii), Proposition 5.1(ii)). Bounding the one-step error by the Lipschitz constant alone gives order $\gamma^2$ per step, hence order $\gamma^{1/2}$ globally; this is Lemma 5.2(i) and the $\mathfrak a=1$ regime, with an extra $\log^3\varepsilon$ in the cost. The gain to $\gamma^3$ needs the second-order structure of $b$ along the Brownian path, through a martingale term and the Laplacian of $b$. Second, the error must be controlled uniformly over infinite horizons. Every estimate has to combine the contraction $(\mathbf C_\alpha)$ with the step-size restriction $\gamma\le\alpha/(2L^2)$; a Gronwall argument on a finite horizon gives constants that grow with time.

For the goal itself, the cost bound (21) is arithmetic once the parameters are fixed. The content is the error bound with a constant that is universal, not merely independent of $\varepsilon$.

Formalization scope

• $\mathbb R^d$ is the published EthierKurtz.SDEState d, and Brownian motions are the published EthierKurtz.IsStandardBrownian.
• The diffusion is a pathwise solution of $X_t=x+\int_0^tb(X_s)ds+\sigma W_t$ with continuous paths. No stochastic integral appears in any statement.
• $\pi$ is a probability measure invariant for the diffusion's transition laws. $\pi^\gamma$ is any probability invariant for the Euler step $x\mapsto x+\gamma b(x)+\sigma\sqrt\gamma Z$, with $Z$ standard Gaussian.
• $\|\cdot\|_2$ is eLpNorm · 2 P, valued in $[0,\infty]$. The schemes are evaluated at grid times, except in $(\mathbf H_4)$ and (47), which use the continuous-time scheme.
• $\mathcal C^2$ with bounded partial derivatives is ContDiff ℝ 2 b plus a uniform bound on the first and second Fréchet derivatives.
• $\|\nabla b\|_{2,\infty}$ and $\|\Delta b\|^2_{2,\infty}$ follow pp. 31–32, which correct the misplaced squares of (18).
• "$\lesssim_{uc}$" is an existential constant quantified before every datum.
• The $\varepsilon_0$ of Theorem 2.3 is read as membership in its defining set.
• The cost formula is used for every $T_r$.
• Added hypotheses: $\sigma>0$ wherever the page leaves it implicit; $\mathcal C^2$ regularity in Lemma 5.2(ii) and Propositions 5.1–5.2; $\gamma<1$ in Proposition 5.2; $d\ge1$ in Proposition 5.1(iii); $\alpha\le1$ in Proposition 2.2.

The goal does not assume $(\mathbf H_1)$–$(\mathbf H_4)$, Proposition 2.2 or any variance or bias estimate. It assumes only $(\mathbf C_\alpha)$, Lipschitz and $\mathcal C^2$ regularity of $b$, $\sigma>0$, invariance of $\pi$ and the parameter choices. Those results appear only as milestones.

Contributions that would be reusable elsewhere: uniform-in-time moment and strong-error bounds for Euler schemes of contractive SDEs, Itô's formula for a $\mathcal C^2$ function of Brownian motion, and the general multilevel complexity theorem (Remark 4.1). Proofs of the milestones are welcome independently.

Selected references

• M. Egéa, F. Panloup, Multilevel-Langevin pathwise average for Gibbs approximation, Math. Oper. Res., 2025; arXiv:2109.07753v2. https://arxiv.org/abs/2109.07753
• M. B. Giles, Multilevel Monte Carlo path simulation, Oper. Res. 56(3), 2008. https://doi.org/10.1287/opre.1070.0496
• A. Durmus, É. Moulines, High-dimensional Bayesian inference via the unadjusted Langevin algorithm, Bernoulli 25(4A), 2019. https://doi.org/10.3150/18-BEJ1073
• A. S. Dalalyan, Theoretical guarantees for approximate sampling from smooth and log-concave densities, J. R. Stat. Soc. B 79(3), 2017. https://arxiv.org/abs/1412.7392
• A. Durmus, É. Moulines, Nonasymptotic convergence analysis for the unadjusted Langevin algorithm, Ann. Appl. Probab. 27(3), 2017. https://arxiv.org/abs/1507.05021

Motivation

Sampling from a Gibbs distribution $\pi_U \propto e^{-U}\,d\lambda_d$ on $\mathbb R^d$, or computing integrals $\pi(f)$ against it, is a basic task of Bayesian statistics, statistical physics and stochastic optimization. A standard route is to simulate a diffusion whose invariant law is the target, for instance the overdamped Langevin equation $dX_t = -\nabla U(X_t)\,dt + \sqrt2\,dB_t$, discretize it with the Euler scheme, and average a function along one long trajectory. Discretization introduces a bias of order of the step size, and the long-time average introduces a variance that decays like the inverse of the horizon, so reaching root-mean-square accuracy $\varepsilon$ with a single Euler scheme costs more than $\varepsilon^{-2}$ steps.

Egéa and Panloup (arXiv:2109.07753, Mathematics of Operations Research, 2025) adapt Giles's multilevel Monte Carlo idea (Giles 2008) to this infinite-horizon setting: they combine a coarse pathwise average with corrections between coupled Euler schemes of successively halved steps, run over decreasing horizons. Their general complexity theorem (Theorem 2.1) is stated under four abstract assumptions (H₁)–(H₄). This mission formalizes the first concrete instance the paper gives: an additive diffusion with a strongly contractive drift, where the assumptions are verified with explicit dependence on the dimension and the cost is of order $\varepsilon^{-2}\log^3(\varepsilon^{-1})$.

Setting

Let $|\cdot|$ and $\langle\cdot,\cdot\rangle$ be the Euclidean norm and inner product on $\mathbb R^d$. Let $b:\mathbb R^d\to\mathbb R^d$ be $L$-Lipschitz and satisfy the contractivity assumption $(\mathbf C_\alpha)$:

$$\langle b(y)-b(x),y-x\rangle\le-\alpha|x-y|^2\qquad(x,y\in\mathbb R^d),$$

with $0<\alpha\le L$. For $b=-\nabla U$ this is uniform strong convexity of $U$. Let $\sigma>0$ and let $W$ be a standard $d$-dimensional Brownian motion. The additive diffusion started at $x$ is $X^x_t = x+\int_0^t b(X^x_s)\,ds+\sigma W_t$; a probability $\pi$ is invariant for it if $\int\mathbb P(X^x_t\in A)\,\pi(dx)=\pi(A)$ for all $t$ and Borel $A$.

The Euler scheme with step $\gamma>0$ is $\bar X^{\gamma,x}_0=x$, $\bar X^{\gamma,x}_{(n+1)\gamma}=\bar X^{\gamma,x}_{n\gamma}+\gamma b(\bar X^{\gamma,x}_{n\gamma})+\sigma(W_{(n+1)\gamma}-W_{n\gamma})$, extended to continuous time by $\bar X_t=x+\int_0^t b(\bar X_{\underline s})\,ds+\sigma W_t$, $\underline s=\gamma\lfloor s/\gamma\rfloor$. A probability $\pi^\gamma$ invariant for one step of the scheme is its invariant distribution.

The multilevel occupation measure with start $x_0$, coarse step $\gamma_0$, levels $r=0,\dots,R$, steps $\gamma_r=\gamma_02^{-r}$, horizons $T_r$ and warm-start time $\tau$ is

$$\mathcal Y(f)=\frac1{T_0-\tau}\int_\tau^{T_0}f(\bar X^{\gamma_0}_{\underline s_{\gamma_0}})\,ds+\sum_{r=1}^R\frac1{T_r-\tau}\int_\tau^{T_r}\big(f(\bar X^{\gamma_r}_{\underline s_{\gamma_{r-1}}})-f(\bar X^{\gamma_{r-1}}_{\underline s_{\gamma_{r-1}}})\big)\,ds,$$

where level $r$ uses its own Brownian motion $B^{(r)}$, independent across levels, and both schemes of a level share it. Its cost is $\mathcal C(\mathcal Y)=T_0/\gamma_0+\sum_{r=1}^R(T_r/\gamma_r+T_r/\gamma_{r-1})$, the number of simulated Euler steps. Write $[f]_1$ for the Lipschitz constant of $f$ and

$$\Upsilon_1^2=1\vee\big(\alpha^{-2}|b(x_0)|^2+\alpha^{-1}\sigma^2d\big).$$

Formalization targets

Goal: Theorem 2.2

There is a universal constant $\mathfrak C_1$ such that, for $\alpha\in(0,L\wedge1]$, $\gamma_0\in(0,\frac\alpha{2L^2}\wedge\frac12]$, the parameters

$$r_0=\Upsilon_1,\quad R_\varepsilon=\lceil2\log_2(r_0\varepsilon^{-1})\rceil,\quad T_r=\frac{\Upsilon_1^2\log(\gamma_0^{-1})}{\alpha}\varepsilon^{-2}R_\varepsilon^22^{-r},\quad \tau=\frac2\alpha|\log\varepsilon|+\frac1\alpha\log\Upsilon_1,$$

and every $\varepsilon\in(0,1)$ with $2\log(\Upsilon_1\varepsilon^{-2})\le\log(\gamma_0^{-1})R_\varepsilon^2$,

$$\|\mathcal Y(f)-\pi(f)\|_2\le\mathfrak C_1[f]_1\varepsilon,\qquad\mathcal C(\mathcal Y)\le\frac{5\log(\gamma_0^{-1})}{2\alpha\gamma_0}\Upsilon_1^2\varepsilon^{-2}\lceil\log_2(\Upsilon_1^2\varepsilon^{-2})\rceil^3.$$

The constant $\mathfrak C_1$ does not depend on $d$, $b$, $\sigma$, $f$, $x_0$ or $\varepsilon$; all dependence on the problem sits in the explicit cost.

Milestones

1. Lemma 5.1 (i)–(iii): uniform second-moment bound (47), existence and uniqueness of $\pi^\gamma$, the $L^2$ contraction $\mathbb E|\bar X^{\gamma,x}_t-\bar X^{\gamma,y}_t|^2\le|x-y|^2e^{-\alpha t}$, the bound $\mathcal W_2(\delta_x,\pi^\gamma)\le c_1(x)$, and from these (H₁) and (H₄).
2. Lemma 5.2 (i): the one-step error between diffusion and Euler step.
3. Proposition 5.1 (i) and (iii) with $\mathfrak a=1$: $\|X^x_{n\gamma}-\bar X^{\gamma,x}_{n\gamma}\|_2^2\le\beta_1(x)\gamma$ uniformly in $n$, hence (H₂).
4. Proposition 5.2 with $\mathfrak a=1$, $\delta=\frac12$: $\mathcal W_1(\pi,\pi^\gamma)\le\sqrt{\beta_1(x^\star)\gamma}$, hence (H₃).
5. Proposition 2.1: all four assumptions with $\max(\alpha c_2^2/L^2,\alpha c_3^2/L^2,c_4^2)\lesssim\Upsilon_1^2$.
6. Remark 4.1: Theorem 2.1 (ii) for the additive scheme when the constants are known up to universal factors.
7. Theorem 2.2 (15): the refinement with $\tilde\Upsilon_1^2=\sigma^2\alpha^{-1}d$ and $\gamma_0=\alpha/(2L^2)$.

Significance

The goal gives a complexity of order $\varepsilon^{-2}\log^3(\varepsilon^{-1})$ for approximating $\pi(f)$ by a pathwise average of Euler schemes, within a logarithmic factor of the $\varepsilon^{-2}$ of an unbiased Monte Carlo method, with constants explicit in $d$, $\alpha$, $L$, $\sigma$. When $|b(x_0)|\lesssim\sqrt d$, the cost is of order $d\varepsilon^{-2}(\log^3d+\log^3\varepsilon)$. The intermediate results (uniform-in-time strong error of the Euler scheme, Wasserstein bias of its invariant law, step-independent contraction) are standard tools of Langevin Monte Carlo and are reusable on their own.

The results are proved in the paper; none of them is machine-checked. A formal proof would certify both the abstract multilevel argument as instantiated here and the explicit constants. The formal statements also pin three printed slips: the ceiling in (14) (cube of the ceiling, not ceiling of the cube), the first line of (12) (the function $c_1$ of Lemma 5.1 (ii)), and the factor in (15).

Difficulty

The obvious argument controls bias and variance separately for each level. The bias of a pathwise average involves the long-time behaviour of the scheme, and the variance of a correcting level involves the $L^2$ distance between two coupled schemes over an unbounded horizon: both require estimates uniform in time, which fail for a merely Lipschitz drift (the strong error of the Euler scheme grows exponentially on finite horizons). Contractivity is what makes the strong error $O(\sqrt\gamma)$ uniformly in $n$, and the delicate point is to obtain this with constants that do not degrade with $d$ or with the step. The complexity bound then requires balancing $R_\varepsilon$, the horizons and the warm start so that the long-time bias, the discretization bias and the variances all sit at level $\varepsilon$, with a universal constant.

Formalization scope

$\mathbb R^d$ is EuclideanSpace ℝ (Fin d) (published EthierKurtz.SDEState); Brownian motions are published EthierKurtz.IsStandardBrownian, one per level, mutually independent as path-valued variables. The diffusion is a pathwise solution of the integral equation with continuous paths; no Itô integral is needed because the noise is additive. Invariance of $\pi$ and $\pi^\gamma$ is stated through transition probabilities. $L^2$ norms are eLpNorm … 2 in $[0,\infty]$, so an infinite moment is never read as zero. $\mathcal W_1$ is used in its Kantorovich–Rubinstein form, $\mathcal W_2(\delta_x,\cdot)$ as an explicit integral, and the clause $\mathcal W_2(\mathcal L(\bar X_{\underline t}),\pi^\gamma)\le\mathcal W_2(\delta_x,\pi^\gamma)e^{-\alpha\underline t}$ as the existence of a coupling with that quadratic cost (Mathlib has no Wasserstein distance; optimal couplings exist, so the forms agree). (H₁)–(H₃) are stated at grid times, (H₄) over all $t\ge0$. "$\lesssim_{uc}$" is an existential constant quantified before every datum.

Standing assumptions carried as hypotheses: $b$ $L$-Lipschitz, $(\mathbf C_\alpha)$, $0<\alpha\le L$ (and $\alpha\le1$ where the paper requires it), $\sigma>0$. Added and disclosed: the §5.1 hypotheses in Lemma 5.2, $\gamma<1$ in Proposition 5.2 (needed by Proposition 5.1), integrability of $f$ under $\pi$, existence of $\pi^\gamma$ and $\tilde c_4>0$ in Remark 4.1. The goal's $\varepsilon_0$ (a maximum of a set that need not be an interval) is replaced by the inequality defining it.

The goal assumes only the contractive setting, the diffusion's invariance of $\pi$ and the parameter choices; it does not assume (H₁)–(H₄), a bias or a variance bound, which appear only as milestones, so it cannot be discharged by restating its hypotheses.

A complete development needs: second-moment and Grönwall estimates for the Euler scheme and the diffusion, existence and uniqueness of invariant laws of Feller chains under a Lyapunov condition, independence and variance decompositions for the multilevel estimator, and the abstract Theorem 2.1. The Gaussian-increment and Wasserstein tools are reusable well beyond this mission. Proofs of any milestone, and a Mathlib-style $\mathcal W_2$ for general laws, are welcome.

Selected references

• M. Egéa, F. Panloup, Multilevel-Langevin pathwise average for Gibbs approximation, Mathematics of Operations Research, 2025; preprint arXiv:2109.07753v2, 2024. https://arxiv.org/abs/2109.07753
• M. B. Giles, Multilevel Monte Carlo path simulation, Operations Research 56(3), 2008. https://doi.org/10.1287/opre.1070.0496
• A. S. Dalalyan, Theoretical guarantees for approximate sampling from smooth and log-concave densities, J. R. Stat. Soc. B 79(3), 2017. https://doi.org/10.1111/rssb.12183
• A. Durmus, É. Moulines, Nonasymptotic convergence analysis for the unadjusted Langevin algorithm, Ann. Appl. Probab. 27(3), 2017. https://doi.org/10.1214/16-AAP1238

Motivation

Sampling from a Gibbs distribution $\pi(dx)\propto e^{-U(x)}dx$ on $\mathbb R^d$, or computing integrals $\pi(f)$ against it, is a basic task in Bayesian statistics and machine learning. A standard approach simulates an ergodic diffusion whose invariant distribution is $\pi$, such as the overdamped Langevin diffusion $dX_t=-\nabla U(X_t)dt+\sqrt2\,dB_t$, discretizes it with an Euler scheme of step $\gamma$, and averages $f$ along one long trajectory. Two errors compete: a bias of order $\gamma^\delta$, because the scheme's invariant distribution $\pi^\gamma$ differs from $\pi$, and a variance of order $1/T$ for a time horizon $T$. With a single step, balancing them requires of the order of $\varepsilon^{-3}$ Euler iterations for an $\varepsilon$-approximation when the bias is of order $\gamma$ (Egéa–Panloup, p. 3).

Multilevel Monte Carlo (Giles 2008, doi:10.1287/opre.1070.0496) removes this extra factor in finite-horizon problems by combining coarse and fine discretizations whose differences have small variance. Pagès and Panloup (arXiv:1609.00412) and Giles, Majka, Szpruch, Vollmer and Zygalakis (arXiv:1805.10974) adapted the idea to invariant distributions. Egéa and Panloup (arXiv:2109.07753, Math. Oper. Res. 2025) build a multilevel pathwise average, a combination of occupation measures of Euler schemes, and prove an $L^2$-error of order $\varepsilon$ at cost $\varepsilon^{-2}$, or $\varepsilon^{-2}|\log\varepsilon|^3$ in a borderline case, under four abstract assumptions. This mission formalizes that general result (their Theorem 2.1 (ii)) and the bias and variance estimates behind it.

Setting

Let $b:\mathbb R^d\to\mathbb R^d$ and $\sigma:\mathbb R^d\to\mathbb M_{d,d}$ be Lipschitz. For a step $\gamma>0$, a starting point $x$ and a $d$-dimensional standard Brownian motion $W$, the Euler scheme is $\bar X^{\gamma,x}_0=x$ and

$$\bar X^{\gamma,x}_{(n+1)\gamma}=\bar X^{\gamma,x}_{n\gamma}+\gamma b(\bar X^{\gamma,x}_{n\gamma})+\sigma(\bar X^{\gamma,x}_{n\gamma})(W_{(n+1)\gamma}-W_{n\gamma}).$$

Write $\underline{s}_\eta=\eta\lfloor s/\eta\rfloor$ and $[f]_1$ for the Lipschitz constant of $f:\mathbb R^d\to\mathbb R$.

Fix $x_0$, a number of levels $R\ge1$, steps $\gamma_r=\gamma_02^{-r}$, horizons $T_0,\dots,T_R$, a warm start $\tau$, and independent Brownian motions $B^{(0)},\dots,B^{(R)}$. The multilevel occupation measure is

$$\mathcal Y(f)=\frac1{T_0-\tau}\int_\tau^{T_0}f(\bar X^{\gamma_0,x_0,B^{(0)}}_{\underline{s}_{\gamma_0}})ds+\sum_{r=1}^R\frac1{T_r-\tau}\int_\tau^{T_r}\Big(f(\bar X^{\gamma_r,x_0,B^{(r)}}_{\underline{s}_{\gamma_{r-1}}})-f(\bar X^{\gamma_{r-1},x_0,B^{(r)}}_{\underline{s}_{\gamma_{r-1}}})\Big)ds .$$

Inside level $r\ge1$ both schemes use the same Brownian motion and are read at the coarse times. The cost is $\mathcal C(\mathcal Y)=T_0/\gamma_0+\sum_{r=1}^R(T_r/\gamma_r+T_r/\gamma_{r-1})$, the number of Euler iterations.

The assumptions involve $\alpha>0$, $\mathfrak a\in[1,2]$, $\eta_0>0$, the invariant distributions $\pi$ of the diffusion and $\pi^\gamma$ of the scheme, and $\|Z\|_2=\mathbb E[|Z|^2]^{1/2}$. For every $\gamma\in(0,\eta_0]$:

• $(\mathbf H_1)$ $|\mathbb E f(\bar X^{\gamma,x}_{n\gamma})-\pi^\gamma(f)|\le c_1(x)[f]_1e^{-\alpha n\gamma}$ (exponential convergence to equilibrium);
• $(\mathbf H_2)$ $\|\bar X^{\gamma,x_0}_{n\gamma}-\bar X^{\gamma/2,x_0}_{n\gamma}\|_2\le c_2\gamma^{\mathfrak a/2}$ ($L^2$-confluence);
• $(\mathbf H_3)$ $|\pi(f)-\pi^\gamma(f)|\le c_3[f]_1\gamma^\delta$, with $\delta\in[1/2,1]$ if $\mathfrak a=1$ and $\delta\in((1+\mathfrak a)/4,1]$ if $\mathfrak a>1$;
• $(\mathbf H_4)$ $\|\bar X^{\gamma,x_0}_{t}-x_0\|_2+\|c_1(\bar X^{\gamma,x_0}_{t})\|_2\le c_4$.

Formalization targets

Goal: Theorem 2.1 (ii)

With $\alpha\le1$, $\eta_0\le\frac12$, $\gamma_0\le\eta_0$, $r_0\ge1\vee c_3\gamma_0^\delta$, $\mathfrak t\ge\mathfrak t_0=\frac{\mathfrak d_{\mathfrak a}}\alpha\max(c_2^2\gamma_0^{\mathfrak a}\log\gamma_0^{-1},c_4^2)$, $\varepsilon\in(0,1)$, $R_\varepsilon=\lceil\delta^{-1}\log_2(r_0/\varepsilon)\rceil$, $T_r=\mathfrak t\varepsilon^{-2}2^{-\frac{1+\mathfrak a}2r}$ ($\mathfrak a>1$) or $\mathfrak t\varepsilon^{-2}R_\varepsilon^22^{-r}$ ($\mathfrak a=1$), and $\tau\in[\tau_1|\log\varepsilon|+\tau_2,\frac12T_{R_\varepsilon}]$:

$$\|\mathcal Y(f)-\pi(f)\|_2\le\mathfrak C_1[f]_1\varepsilon,\qquad \mathcal C(\mathcal Y)\le\mathfrak c_{\mathfrak a}\gamma_0^{-1}\mathfrak t\,\varepsilon^{-2}\ (\mathfrak a>1),\quad\le\mathfrak c_{\mathfrak a}\gamma_0^{-1}\mathfrak t\,\varepsilon^{-2}R_\varepsilon^3\ (\mathfrak a=1),$$

with $\mathfrak C_1$ a universal constant. The cost constant $\mathfrak c_{\mathfrak a}$ is explicit; the error constant is left existential, which is how the paper states it.

Companion: Theorem 2.1 (i)

With the same steps and horizons but $r_0\ge1$, $\mathfrak t>0$ arbitrary and $\tau\in[\tau_1\log(\varepsilon^{-1})\wedge\frac12T_{R_\varepsilon},\frac12T_{R_\varepsilon}]$ for a fixed $\tau_1>\frac{1+\mathfrak a-2\delta}{2\alpha\delta}$, the same two bounds hold with an error constant $\mathfrak C_1$ that depends on the parameters but not on $f$, $\varepsilon$ or $\tau$. It is a separate theorem item, not a milestone.

Milestones

Lemma 3.1 (long-time bias of one average), Proposition 3.2 (bias of $\mathcal Y$), the variance decomposition (32), Lemma 3.2 (variance of one correcting level), Proposition 3.3 (variance of $\mathcal Y$), the cost bound (28), and Proposition 3.1, the non-asymptotic squared-error bound for given steps and horizons.

Significance

Theorem 2.1 separates the analysis of the multilevel estimator from the analysis of any particular diffusion. Once a diffusion satisfies $(\mathbf H_1)$–$(\mathbf H_4)$, the complexity follows with explicit parameters. The paper then checks the assumptions for Langevin diffusions with strongly convex potentials and obtains dimension-explicit costs $\mathcal O(d\varepsilon^{-2}\log^3(d\varepsilon^{-2}))$ and $\mathcal O(d\varepsilon^{-2})$. Those are the companion missions II and III of this series; they rest on the same estimator and cost.

The paper's proofs are on paper only. No machine-checked proof of a multilevel ergodic estimator is known to us. This mission produces a statement-level formalization of the estimator, its cost, the four assumptions and the bias–variance bounds. A complete development would add discrete-time Markov-property arguments for Euler schemes driven by a Brownian motion, covariance bounds for ergodic averages, and the elementary but delicate parameter calculus of §4.

Difficulty

The bias step is close to routine once $(\mathbf H_1)$ and $(\mathbf H_3)$ are given. The variance of a correcting level is not. Bounding it by the second moment alone gives $(T-\tau)^{-1}\cdot T\gamma^{\mathfrak a}$, which does not decay with $T$. The proof must also use mixing: covariances between times $s$ and $u$ far apart are controlled by $(\mathbf H_1)$ through the Markov property, which needs the conditional law of each scheme given its past and a moment bound on $c_1$ along the path $(\mathbf H_4)$. The two regimes are matched at a separation of order $|\log\gamma|/\alpha$, and this is where the factor $\log(1/\gamma)$ comes from. At the level of the theorem, the parameters (8) and the warm start $\tau$ must make three error terms simultaneously $\lesssim\varepsilon$ with a constant independent of $d$, $f$, the diffusion and every assumption constant. Getting the uniformity right is what the explicit thresholds $\mathfrak t_0$, $\tau_1$, $\tau_2$ are for.

Formalization scope

• $\mathbb R^d$ is the published EthierKurtz.SDEState d, and a Brownian motion is the published EthierKurtz.IsStandardBrownian. The levels' Brownian motions form one family B : ℕ → ℝ≥0 → Ω → ℝ^d, mutually independent as path-valued variables. $\sigma(x)$ is a continuous linear map $\mathbb R^d\to\mathbb R^d$.
• The Euler scheme is evaluated at grid times. $(\mathbf H_1)$–$(\mathbf H_4)$ are stated at grid times $n\gamma$, for every level's Brownian motion. For $(\mathbf H_4)$, whose supremum the paper takes over all real $t$, this is a weaker assumption, so the theorems are stronger.
• $\|\cdot\|_2$ is eLpNorm · 2 P and variances are evariance, both valued in $[0,\infty]$, so an infinite moment is never read as $0$. Expectations that are concluded come with an integrability conjunct.
• $\pi$ and $\pi^\gamma$ are arbitrary probability measures satisfying $(\mathbf H_1)$ and $(\mathbf H_3)$, with $f$ integrable. The SDE and the invariance of $\pi$, $\pi^\gamma$ are not formalized: the results never use invariance, so the statements generalize the paper's.
• $[f]_1$ is the least Lipschitz constant. The cost is the formula of p. 6 for all horizons. $(\mathfrak a-1)^{-2}\wedge R^2$, $\mathfrak d_{\mathfrak a}$ and $\mathfrak c_{\mathfrak a}$ are defined case by case, so $\mathfrak a=1$ never divides by zero.
• Universal constants ($\mathfrak C_1$, $\mathfrak c_u$) are quantified before the dimension, the probability space and every datum.
• Added hypotheses: $\gamma\le1/2$ and $\alpha\le1$ in Lemma 3.2, and $\eta_0\le1/2$ in Proposition 3.3. These are the standing ranges of Theorem 2.1; without them the printed bounds fail for steps $\ge1$. The restriction $\varepsilon<\varepsilon_0$ of the goal is replaced by quantifying over every $\tau$ in the admissible interval, which is the same claim.
• A formalization in which the bias or variance bounds, Proposition 3.1 or any covariance estimate appear as hypotheses of the goal is excluded. The goal assumes only $(\mathbf H_1)$–$(\mathbf H_4)$ and the parameter choices (8), (39)–(40).

Contributions are welcome on any milestone, on Markov-property and moment lemmas for Euler schemes driven by IsStandardBrownian, and on the Theorem 2.1 (i) companion with non-explicit constants.

Selected references

• M. Egéa, F. Panloup, Multilevel-Langevin pathwise average for Gibbs approximation, arXiv:2109.07753v2, 2024; Math. Oper. Res., 2025. https://arxiv.org/abs/2109.07753
• G. Pagès, F. Panloup, Weighted multilevel Langevin simulation of invariant measures, Ann. Appl. Probab. 28(6), 2018. https://arxiv.org/abs/1609.00412
• M. B. Giles, M. B. Majka, L. Szpruch, S. Vollmer, K. Zygalakis, Multi-level Monte Carlo methods for the approximation of invariant measures of stochastic differential equations, Stat. Comput. 30, 2020. https://arxiv.org/abs/1805.10974
• M. B. Giles, Multilevel Monte Carlo path simulation, Oper. Res. 56(3), 2008. https://doi.org/10.1287/opre.1070.0496

Motivation

Many-server queues are the basic model of call centres and hospital wards, and the Halfin–Whitt regime (the number of servers $n$ grows while the traffic intensity approaches one at rate $n^{-1/2}$) is where such systems operate efficiently and with short delays. In this regime the centred and scaled queue process of a $G/Ph/n+GI$ queue (general arrivals, phase-type service, $n$ servers, generally distributed patience) converges to a $d$-dimensional diffusion, $d$ being the number of service phases (Dai, He and Tezcan, Ann. Appl. Probab. 2010). For $GI/Ph/n+M$ queues the stationary distribution of the queue converges to the invariant measure $\mu$ of that diffusion (Dai, Dieker and Gao, Queueing Systems 2014). So $\mu$ is the object one would like to compute.

With exponential service the diffusion is one-dimensional and $\mu$ is explicit. With phase-type service it is multidimensional with a piecewise-linear drift; Dieker and Gao (Ann. Appl. Probab. 2013) proved that it is exponentially ergodic with a unique $\mu$, but no closed form is known, and the characterization of multidimensional piecewise-linear diffusions was left open in Browne and Whitt (1995). Jin, Pang, Xu and Xu (arXiv:2209.07361) approximate $\mu$ by the invariant measure of the Euler–Maruyama discretization and give an explicit error rate in the Wasserstein-1 distance.

Setting

Fix $d\ge1$, a patience rate $\alpha>0$, a slack $\beta\in\mathbb R$, a probability vector $p\in\mathbb R^d$, service rates $v_1,\dots,v_d>0$ and a sub-stochastic phase-transition matrix $P$ with $P_{ii}=0$ and $I-P$ invertible. Let $R=(I-P')\,\mathrm{diag}(v)$, $e=(1,\dots,1)'$, and assume the phase-type distribution has mean one, $e'R^{-1}p=1$. The drift is the piecewise-linear map

$$g(x)=-\beta p-Rx+(R-\alpha I)p\,(e'x)^+,\qquad y^+=\max\{0,y\}.$$

Let $\sigma$ be a $d\times d$ matrix with $\xi'\sigma\sigma'\xi\ge c\,\xi'\xi$ for some $c>0$ (in the queueing model $\sigma\sigma'$ has an explicit formula in terms of the arrival variability $c_a^2$, $p$, $P$, $v$; only the ellipticity is used). The limiting diffusion is

$$X_t=x+\int_0^t g(X_s)\,ds+\sigma B_t,$$

with $B$ a standard $d$-dimensional Brownian motion; $P_t(x,\cdot)$ is the law of $X_t$ started at $x$, and $\mu$ is a probability measure with $\int P_t(x,\cdot)\,\mu(dx)=\mu$ for all $t\ge0$.

The Euler–Maruyama chain with step $\eta>0$ is

$$\tilde X^\eta_{k+1}=\tilde X^\eta_k+g(\tilde X^\eta_k)\,\eta+\sqrt\eta\,\sigma\xi_{k+1},\qquad k\in\mathbb N_0,$$

with independent standard Gaussian $\xi_k$; its one-step kernel is $\tilde P_\eta(x,\cdot)=\mathcal N(x+\eta g(x),\eta\sigma\sigma')$, and $\tilde\mu_\eta$ denotes an invariant probability measure of it. The Wasserstein-1 distance is taken in its dual form

$$d_W(\mu_1,\mu_2)=\sup_{h\in\mathrm{Lip}(1)}\Big|\int h\,d\mu_1-\int h\,d\mu_2\Big|,$$

where $\mathrm{Lip}(1)$ is the set of $h$ with $|h(x)-h(y)|\le|x-y|$.

Formalization targets

Goal: Theorem 1.2(i)

There is $\eta_0>0$ such that for every $\varsigma\in(0,1)$ there is $C_\varsigma>0$, depending on $\varsigma$ but not on $\eta$, such that for every $\eta\in(0,\eta_0]$ the chain has a unique invariant probability measure $\tilde\mu_\eta$ and

$$d_W(\mu,\tilde\mu_\eta)\le C_\varsigma\,\eta^{\frac{1-\varsigma}{2}}.$$

The constants are left unspecified, so the goal survives any improvement of them.

Milestones

In the order the proof uses them:

1. (A.4)–(A.7): a $C^2$ Lyapunov function $V$ with $\hat c_1|y|^2\le V(y)\le\hat C|y|^2+\hat C'$, $|\nabla V(y)|\le C(1+|y|)$ and $\mathcal AV\le-c_1V+\breve c_1$, where $\mathcal A$ is the generator.
2. Lemma A.1: moment bounds $\mathbb EV^\ell(X^x_t)\le e^{-c_1t}V^\ell(x)+\breve c_\ell(1-e^{-c_1t})/c_1$, and exponential ergodicity of the diffusion in $d_W$ and in weighted total variation.
3. Lemma A.2: $\tilde P_\eta f$ is Lipschitz for bounded measurable $f$ (strong Feller), and the chain is irreducible.
4. Lemma A.3: the one-step error $\mathbb E|X^x_\eta-\tilde X^{\eta,x}_1|^{2\ell}\le\tilde C_\ell(1+V^\ell(x))\eta^{3\ell}$.
5. Lemma A.4: the drift condition $\tilde P_\eta V^\ell\le(e^{-c_1\eta}+\tilde C_\ell\eta^{3/2})V^\ell+\check K_\ell$ for the chain.
6. Proposition 2.1: existence, uniqueness and exponential ergodicity of $\tilde\mu_\eta$, and $\tilde\mu_\eta(|\cdot|^\ell)\le C$ uniformly in $\eta$.
7. Lemma 2.2: for the solution $f$ of the Stein equation $\mathcal Af=h-\mu(h)$ with $h\in\mathrm{Lip}_0(1)$, the polynomial bounds on $f,\nabla f,\nabla^2f$ and the local $(1-\varsigma)$-Hölder bound on $\nabla^2 f$.
8. (2.7)–(2.8): the second-order Taylor remainder in Stein's argument is $O(\eta^{(3-\varsigma)/2})$.

Significance

The theorem gives a computable approximation of the stationary distribution of the diffusion limit, with an explicit rate. Combined with the convergence of the scaled queue's stationary distribution to $\mu$ at rate $n^{(\varsigma-1)/2}$ for $M/Ph/n+M$ queues (Braverman and Dai), it yields an approximation of the steady state of the queue itself, by simulating the chain for $O(n\log n)$ steps (Remark 2.3 of the paper). Proposition 2.1 also shows that the explicit Euler scheme does not blow up for this non-gradient, non-differentiable drift, which is not automatic: the Euler scheme of an ergodic SDE may blow up.

The result is proved in the paper; nothing here is open mathematically, though Lemma 2.2 is imported from Gurvich's regularity theory rather than proved there. To our knowledge none of it is formalized. A formalization provides the first machine-checked Wasserstein bound between the invariant measures of a multidimensional diffusion and its discretization, and on the way the Lyapunov and Harris-type ergodicity theory of Euler–Maruyama chains.

Difficulty

Standard error analyses of Euler schemes for invariant measures assume a gradient drift, strong dissipativity, or several bounded derivatives; none holds here. The drift is not differentiable on the hyperplane $\{e'x=0\}$ and is dissipative only in the Lyapunov sense of Dieker and Gao, so contraction-based couplings do not apply. The Stein argument instead needs second derivatives of the Stein solution $f$ that are Hölder continuous with polynomially growing constants; this is where the loss $\eta^{-\varsigma/2}$ comes from, and it rests on elliptic regularity for an operator with only Lipschitz coefficients. The ergodicity of the chain must also hold with constants uniform in $\eta$, which the classical Lyapunov criteria do not give directly.

Formalization scope

• The state space is EuclideanSpace ℝ (Fin d), so $|\cdot|$ is the Euclidean norm; matrices act by the usual product, and $\|A\|_{\mathrm{op}}=\sup_{|u|=1}|Au|$.
• $\sigma$ is a general uniformly elliptic matrix, not the paper's explicit covariance; the general case contains the paper's.
• Brownian motion is the published EthierKurtz.IsStandardBrownian (time ℝ≥0). The diffusion is a pathwise integral equation with continuous paths, so no Itô integral is needed. The transition function $P_t$ is characterized by: on every probability space, every solution from $x$ has law $P_t(x,\cdot)$ at time $t$. $\mu$ is a hypothesis (an invariant probability measure of $P_t$), never assumed smooth or with moments.
• The chain's kernel is the image of stdGaussian under $z\mapsto x+\eta g(x)+\sqrt\eta\,\sigma z$; $k\in\mathbb N_0$ starts at $0$.
• $d_W$, $\|\cdot\|_{\mathrm{TV}}$ and $\|\cdot\|_{\mathrm{TV},V}$ take values in $[0,\infty]$ and equal $+\infty$ unless the needed first moment (or $\int(1+V)$) is finite; so a bound on $d_W$ cannot hold because an integral of a non-integrable test function defaults to $0$. All moment bounds are lower Lebesgue integrals.
• "For any small enough $\varsigma$" is read as every $\varsigma\in(0,1)$; the step threshold $\eta_0$ is chosen before $\varsigma$, and $C_\varsigma$ before $\eta$.
• The Stein solution is $f=-\int_0^\infty(P_th-\mu(h))\,dt$ (the paper's (3.2)); Lemma 2.2 concludes that it is $C^2$ and solves $\mathcal Af=h-\mu(h)$, so its derivatives are never junk.
• Lemma A.2's gradient bound is stated in Lipschitz form, since $\tilde P_\eta f$ need not be differentiable where $g$ is not.

Trivializing encodings are excluded: the existence and uniqueness of $\tilde\mu_\eta$ is a conclusion, not a hypothesis, and the distance is infinite rather than zero when moments fail.

A complete development needs Brownian motion with Gaussian increments, a pathwise SDE solution theory for Lipschitz drift, Itô's formula (or Dynkin's formula) for $C^2$ functions of additive-noise diffusions, Harris-type ergodic theorems for Feller chains and diffusions, and elliptic regularity for the Stein equation. The Harris theorems and Wasserstein machinery are reusable well beyond this mission. Contributions to any milestone, and to general infrastructure for Euler–Maruyama chains, are welcome.

Selected references

• X. H. Jin, G. Pang, L. Xu, X. Xu, An approximation to the invariant measure of the limiting diffusion of G/Ph/n+GI queues in the Halfin–Whitt regime and related asymptotics, arXiv:2209.07361v1, 2022. https://arxiv.org/abs/2209.07361
• A. B. Dieker, X. Gao, Positive recurrence of piecewise Ornstein–Uhlenbeck processes and common quadratic Lyapunov functions, Annals of Applied Probability 23(4): 1291–1317, 2013.
• J. G. Dai, S. He, T. Tezcan, Many-server diffusion limits for G/Ph/n+GI queues, Annals of Applied Probability 20(5): 1854–1890, 2010.
• J. G. Dai, A. B. Dieker, X. Gao, Validity of heavy-traffic steady-state approximations in many-server queues with abandonment, Queueing Systems 78(1): 1–29, 2014.
• S. Browne, W. Whitt, Piecewise-linear diffusion processes, in J. Dshalalow (ed.), Advances in Queueing: Theory, Methods, and Open Problems, 463–480, 1995.
• I. Gurvich, Diffusion models and steady-state approximations for exponentially ergodic Markovian queues, Annals of Applied Probability 24(6): 2527–2559, 2014.
• A. Braverman, J. G. Dai, Stein's method for steady-state diffusion approximations of M/Ph/n+M systems, Annals of Applied Probability 27(1): 550–581, 2017.
• M. Hairer, J. C. Mattingly, Spectral gaps in Wasserstein distances and the 2D stochastic Navier–Stokes equations, Annals of Probability 36(6): 2050–2091, 2008.

Bibliographic data as in the reference list of arXiv:2209.07361v1, pp. 35–37.