Stochastic Systems

Motivation

Many equilibrium problems in operations research have the form of a stochastic Nash game: each of $N$ players minimizes an expected cost that depends on the strategies of the others, and the expectation can only be sampled. The paper's introduction lists applications in generation-capacity and power markets and in communication networks (Abada, de Maere d'Aertrycke and Smeers, 2017; Başar, 2007). A natural way to compute an equilibrium is a best-response scheme: each player repeatedly re-optimizes against the current strategies of its rivals. In a large network the players cannot synchronize. They update at different times and see their rivals' strategies with delays.

Lei, Shanbhag, Pang and Sen (arXiv:1704.04578v2; Mathematics of Operations Research, 2020) analyze inexact proximal best-response schemes in three regimes: synchronous, randomized and asynchronous. This mission formalizes the asynchronous one (§5 and Appendix C). It adapts the partially asynchronous iterations of Bertsekas and Tsitsiklis (1989) to stochastic games with inexact, sampled best responses, and it gives an explicit bound on the number of stochastic gradient steps each player needs.

Setting

Player $i$ chooses $x_i$ in a nonempty compact convex set $X_i \subseteq \mathbb R^{n_i}$. A profile is $x = (x_1,\dots,x_N) \in X = \prod_i X_i$. Player $i$'s cost is $f_i(x_i, x_{-i}) = \mathbb E[\psi_i(x_i, x_{-i};\xi)]$, which is convex and twice continuously differentiable in a neighbourhood of $X$. A sampling oracle returns $\nabla_{x_i}\psi_i(x;\xi)$, whose second moment is at most $M_i^2$ (Assumption 1). A Nash equilibrium $x^*$ is a profile in which each $x^*_i$ minimizes $f_i(\cdot, x^*_{-i})$ over $X_i$.

For $\mu > 0$, the proximal best response of player $i$ to a profile $y \in X$ is

$$\hat x_i(y) = \operatorname*{argmin}_{x_i \in X_i}\Big[f_i(x_i, y_{-i}) + \tfrac{\mu}{2}\|x_i - y_i\|^2\Big].$$

The curvature constants $\zeta_{i,\min} = \inf_{x\in X}\lambda_{\min}(\nabla^2_{x_i}f_i(x))$ and $\zeta_{ij,\max} = \sup_{x\in X}\|\nabla^2_{x_ix_j}f_i(x)\|$ define the matrix $\Gamma$ with $\gamma_{ii} = \mu/(\mu+\zeta_{i,\min})$ and $\gamma_{ij} = \zeta_{ij,\max}/(\mu+\zeta_{i,\min})$. Assumption 5 (strict diagonal dominance) asks $\zeta_{i,\min} > \sum_{j\ne i}\zeta_{ij,\max}$. It makes $a_\infty = \|\Gamma\|_\infty$, the maximum absolute row sum, smaller than one.

The asynchronous scheme (Algorithm 3) runs on a deterministic schedule. At time $k$ the players in $I_k$ update. Player $i \in I_k$ sees the outdated profile $y^i_k = (x_{1,k-\tau_{i1}(k)},\dots,x_{N,k-\tau_{iN}(k)})$ and computes $x_{i,k+1} \in X_i$ with

$$\mathbb E\big[\|x_{i,k+1} - \hat x_i(y^i_k)\|^2 \mid \mathcal F_k\big] \le \alpha_{i,k}^2 .$$

The other players keep their strategies. Assumption 4 requires every player to update at least once in every $B_1$ consecutive times, and every delay to be at most $B_2$. The accuracy is $\alpha_{i,k} = \eta^{\beta_{i,k}}$, where $\beta_{i,k}$ counts player $i$'s updates so far. The inexact response is computed by $j_{i,k} = \lceil Q_i/\eta^{2(k+1)}\rceil$ projected stochastic gradient steps (43), with step size $1/(\mu(t+1))$ and $Q_i = 2M_i^2/\mu^2 + 2D_{X_i}^2$. Write $n_0 = \lceil B_2/B_1\rceil$, $\rho = (\max\{a_\infty,\eta\})^{1/(n_0+1)}$ and $c = \rho^{1/B_1}$.

Formalization targets

Goal: Theorem 3, (45)

Let $q \in (c,1)$, $D \ge 1/(e\ln(q/c))$, $\epsilon > 0$, and

$$\hat\epsilon = \frac{\epsilon}{C+D}\,\rho^{\frac{B_1-1}{B_1}} < 1,$$

where $\mathbb E\|x_{i,0}-x^*_i\|^2 \le C^2$. After $K = \lceil \ln(1/\hat\epsilon)/\ln(1/q)\rceil$ major iterations, $\max_i \mathbb E\|x_{i,K}-x^*_i\| \le \epsilon$, and player $i$ has taken at most

$$\ell^{(1)}_i(\eta) = \frac{Q_i}{\eta^4\ln(1/\eta^2)}\Big(\frac1{\hat\epsilon}\Big)^{\frac{\ln(1/\eta^2)}{\ln(1/q)}} + \Big\lceil\frac{\ln(1/\hat\epsilon)}{\ln(1/q)}\Big\rceil$$

projected gradient steps.

Milestones

In the order the proof uses them:

• the $\Gamma$-contraction (5) of the proximal best response;
• the fixed-point identity $x^*_i = \hat x_i(x^*)$;
• $a_\infty < 1$ under Assumption 5;
• its expected $\infty$-norm form (40), $\mathbb E\|\hat x_i(y)-x^*_i\| \le a_\infty\max_j\mathbb E\|y_j-x^*_j\|$;
• the one-step recursion (C.2) and the delay step (C.3) of Appendix C;
• Lemma 2, $zc^z \le Dq^z$;
• Lemma 7, the rate

$$\max_i\mathbb E\|x_{i,k}-x^*_i\| \le (C+k)\rho^{\lfloor k/B_1\rfloor},\qquad \max_i\mathbb E\|x_{i,k}-x^*_i\| \le \rho^{-\frac{B_1-1}{B_1}}(C+D)q^k;$$

• Lemma 8, the $Q_i/(t+1)$ mean-square error of the inner stochastic approximation loop;
• the summation bound (30).

Significance

Theorem 3 makes the cost of asynchrony explicit. Bounded delays and a bounded update window enter only through the rate per window, $\max\{a_\infty,\eta\}^{1/(n_0+1)}$. The overall effort stays polynomial in $1/\epsilon$, and the paper's Table 1 summarizes the exponent as $2B_1(1+\lceil B_2/B_1\rceil)+\delta$ for a suitable choice of $q$. The result connects classical partially asynchronous fixed-point theory with the sample-complexity analysis of stochastic approximation. Lemma 7 holds for any inexact solver that achieves the accuracy (38), not only for the stochastic gradient loop, so it applies to other subproblem solvers as well.

The paper proves these results by hand. As far as is known, no part of them has been machine-checked: there is no formal treatment of stochastic Nash games, of proximal best-response maps, or of asynchronous iterations with delays in Mathlib or on this platform. A formalization would check the induction over windows in Appendix C, whose index bookkeeping has several printed slips. It would also make precise which properties of the delays the argument needs.

Difficulty

The obvious argument does not carry over from the synchronous case. There, one step of the scheme contracts the whole error vector by $a$ in one norm. Here a step updates only some players, and each of them uses information up to $B_2$ steps old. No single step contracts anything. The proof has to show that every window of $B_1$ steps contracts the worst-case expected error by a factor $\rho$, with the delays absorbed into the root $1/(n_0+1)$. This requires a nested induction with case distinctions on the window position (Appendix C). The analysis also mixes three layers: a deterministic contraction of the exact best response; conditional expectations, Jensen's inequality and adaptedness for the inexact responses; and an $O(1/t)$ stochastic-approximation bound for the inner loop, which in turn needs strong convexity, the first-order optimality conditions and nonexpansiveness of the projection.

Formalization scope

Players are indexed by Fin N. Player $i$'s space is EuclideanSpace ℝ (Fin (n i)), and costs are functions of the whole profile: $f_i(z, y_{-i})$ is f i (Function.update y i z). The projection $\Pi_{X_i}$ is any map satisfying the published predicate SpectralProjGrad.Shared.IsProjOnto. The proximal best response is a map constrained by an argmin predicate on $X$. $\zeta_{i,\min}$ and $\zeta_{ij,\max}$ are Rayleigh-quotient infima and suprema of second Fréchet derivatives, and $\|\Gamma\|_\infty$ is the maximum absolute row sum. The probability space carries a filtration $(\mathcal F_k)$ to which the iterates are adapted. The inner-loop $\sigma$-algebras $\sigma\{\mathcal F_k,\xi^{[t-1]}_{i,k}\}$ are parameters satisfying the inclusions the proof uses. Conditional expectations are Mathlib's condExp.

Statement repairs and implicit hypotheses, each recorded in the item it affects:

• Deterministic delays. Assumption 4(c) allows random delays, but step (C.5) of the proof bounds an expectation at a random past time by the maximum over past times. That step is valid only when the delays do not depend on the iterates. The delays here are deterministic functions $\tau_{ij}(k) \le k$ with $\tau_{ii}(k) = 0$.
• Oracle second moment. Lemma 8 prints $\psi_i$ for $\nabla_{x_i}\psi_i$, and its hypothesis $\mathbb E[\|\nabla_{x_i}\psi_i\|^2\mid\cdot] = \|\nabla_{x_i}f_i\|^2$ forces zero noise. The bound $\le M_i^2$ that the proof uses replaces it.
• Regularity. Joint $C^2$ regularity of $f_i$ is assumed, which the mixed Hessian blocks in $\Gamma$ need, and $X_i$ is nonempty.
• Constants. Theorem 3 assumes $C \ge 0$ and $\hat\epsilon < 1$. Its "$D \ge /\ln((q/c)^e)$" is read as $D \ge 1/(e\ln(q/c))$, as in Lemma 7.
• Step count. Theorem 3 counts the steps over player $i$'s update times; the proof's bound is for the larger sum over all times.

The special case $B_1 = 1$, $B_2 = 0$ of Theorem 3, Corollary 3 (cyclic updates) and Theorem 4 are not stated.

The goal cannot be met trivially. The schedule is a parameter, not "every player at every step". The sampled gradients may be genuinely random, since the second-moment hypothesis is an inequality. $\hat x$ is tied to $f_i$ by its argmin property, and $x^*$ is a Nash equilibrium. Each expectation in a conclusion is of a bounded measurable quantity, so it cannot vanish by non-integrability.

The development needs the following:

• the contraction (5), which the paper only cites from Facchinei and Pang (2009, §12.6.1) and which uses the mean-value theorem along segments in $X$;
• strong-convexity optimality conditions on convex sets;
• nonexpansiveness of Euclidean projections;
• conditional Jensen and tower arguments;
• an integer-index induction over windows.

The projection and strong-convexity lemmas, and (5) itself, are reusable for the synchronous and randomized missions of this series. Proofs of individual milestones, or of general Mathlib-level facts such as the first-order optimality condition over a convex set, are welcome.

Selected references

• J. Lei, U. V. Shanbhag, J.-S. Pang, S. Sen, On Synchronous, Asynchronous, and Randomized Best-Response Schemes for Stochastic Nash Games, arXiv:1704.04578v2, 2018; Mathematics of Operations Research, 2020, https://doi.org/10.1287/moor.2018.0986. https://arxiv.org/abs/1704.04578v2
• D. P. Bertsekas, J. N. Tsitsiklis, Parallel and Distributed Computation: Numerical Methods, Prentice Hall, 1989 (reference [9] of the paper).
• F. Facchinei, J.-S. Pang, Nash equilibria: the variational approach, in Convex Optimization in Signal Processing and Communications, Cambridge University Press, 2009 (reference [19] of the paper).
• I. Abada, G. de Maere d'Aertrycke, Y. Smeers, On the multiplicity of solutions in generation capacity investment models with incomplete markets, Mathematical Programming 165(1):5–69, 2017 (reference [1] of the paper).

Motivation

Many stochastic differential equations used in physics, biology and Bayesian sampling are ergodic: their law converges to an invariant measure $\pi$, and quantities of interest are averages $\pi(\varphi)=\int\varphi\,\mathrm d\pi=\lim_{t\to\infty}\mathbb E[\varphi(X_t)]$. Estimating such averages by simulation requires a numerical method whose accuracy does not deteriorate as the simulated time grows. Two obstacles meet here. First, the drifts of interest (the FENE polymer model, Langevin dynamics with polynomial potentials) are not globally Lipschitz, and for such drifts the explicit Euler–Maruyama method with a uniform timestep has moments that blow up as the step tends to zero (Hutzenthaler, Jentzen, Kloeden 2011). Second, standard strong error bounds grow exponentially with the time horizon, which makes them useless for long-time averages.

W. Fang and M. B. Giles (Ann. Appl. Probab. 30 (2020) 526–560) answer both with an adaptive timestep: the step is chosen from the current state, small where the drift is large. On a finite interval they prove strong convergence of order $\tfrac12$; for a class of contractive ergodic SDEs they prove that the moments of the numerical solution and its strong error are bounded uniformly in time. The second result is the subject of this mission. It underlies the paper's multilevel Monte Carlo estimator for invariant measures.

Earlier work on the same problem used tamed schemes, which modify the drift instead of the step: Hutzenthaler, Jentzen, Kloeden 2012 on finite intervals and Sabanis 2016 with a time-uniform rate. An adaptive Euler–Maruyama scheme of a different form, with error control, was shown to converge and to be stable by Lamba, Mattingly, Stuart 2007.

Setting

Let $(\Omega,\mathcal F,\mathbb P)$ carry a filtration $(\mathcal F_t)_{t\ge0}$ and a $d$-dimensional standard Brownian motion $W$ with respect to it. On $\mathbb R^m$ (Euclidean norm $\|\cdot\|$, inner product $\langle\cdot,\cdot\rangle$) consider

$$\mathrm dX_t=f(X_t)\,\mathrm dt+g(X_t)\,\mathrm dW_t,\qquad X_0=x_0,$$

with drift $f:\mathbb R^m\to\mathbb R^m$ and diffusion $g:\mathbb R^m\to\mathbb R^{m\times d}$; matrices carry the Frobenius norm $\|A\|=(\sum_{i,j}A_{ij}^2)^{1/2}$.

Given a timestep function $h:\mathbb R^m\to(0,\infty)$, the adaptive Euler–Maruyama scheme is

$$t_0=0,\ \widehat X_0=x_0,\qquad t_{n+1}=t_n+h(\widehat X_{t_n}),\qquad \widehat X_{t_{n+1}}=\widehat X_{t_n}+f(\widehat X_{t_n})h(\widehat X_{t_n})+g(\widehat X_{t_n})(W_{t_{n+1}}-W_{t_n}).$$

For $t\ge0$ let $\underline t$ be the last grid time not after $t$. The piecewise-constant interpolant is $\overline X_t=\widehat X_{\underline t}$ and the continuous interpolant is $\widehat X_t=\widehat X_{\underline t}+f(\widehat X_{\underline t})(t-\underline t)+g(\widehat X_{\underline t})(W_t-W_{\underline t})$.

The hypotheses are:

• Assumption 7 (dissipative). $f,g$ are locally Lipschitz, and for constants $\alpha,\beta>0$: $\langle x,f(x)\rangle\le-\alpha\|x\|^2+\beta$ and $\|g(x)\|^2\le\beta$.
• Assumption 8 (timestep). $h$ is continuous with values in $(0,h_{\max}]$, and $\langle x,f(x)\rangle+\tfrac12h(x)\|f(x)\|^2\le-\alpha\|x\|^2+\beta$.
• Assumption 9 (contractive). For some $p^*>2$ and $\lambda,\eta>0$: $\langle x-y,f(x)-f(y)\rangle+\tfrac{p^*-1}2\|g(x)-g(y)\|^2\le-\lambda\|x-y\|^2$, $\|g(x)-g(y)\|^2\le\eta\|x-y\|^2$, and $\|f(x)-f(y)\|\le(\gamma(\|x\|^q+\|y\|^q)+\mu)\|x-y\|$ with $\gamma,\mu,q>0$.
• Assumption 3 (refinement). For $\delta\in(0,1]$ and a fixed $T>0$, the refined timestep $h^\delta$ satisfies $\delta\min(T,h(x))\le h^\delta(x)\le\min(\delta T,h(x))$.

Formalization targets

Goal: Theorem 6 (p. 534)

Under Assumption 9, with $g$ bounded, $h$ satisfying Assumption 8 and $h^\delta$ satisfying Assumption 3, for every $p\in(0,p^*]$ there is $C_p$ such that, with $\widehat X$ computed with timestep $h^\delta$,

$$\mathbb E\big[\|\widehat X_t-X_t\|^p\big]\le C_p\,\delta^{p/2}\qquad\text{for all }\delta\in(0,1],\ t\ge0.$$

The constant is fixed before $\delta$ and $t$; uniformity in $t$ is the content of the theorem.

Milestones, in the order the proof uses them

1. Lemma 3 (p. 532): under Assumption 7, $\sup_{t\ge0}\mathbb E\|X_t\|^p<\infty$ for every $p>0$.
2. (45) (p. 554): the weighted one-step inequality bounding $e^{2\alpha t_{n+1}}\|\widehat X_{t_{n+1}}\|^2$ by $e^{2\alpha t_n}\|\widehat X_{t_n}\|^2$ plus explicit drift, noise and martingale terms.
3. (47) (p. 554): its partial-step analogue from $\underline t$ to $t$.
4. §6.3, Step 2 (p. 556): $\mathbb E[\sup_{0\le s\le t}e^{\alpha ps}\|\widehat X_s\|^p]\le2C^4_p\|x_0\|^p+2C^5_pe^{\alpha pt}$ for $p\ge4$, with $C^4_p,C^5_p$ independent of $x_0$ and uniform over the coefficients $f,g$ (given the constants of Assumption 7) and over timestep functions.
5. Theorem 5 (p. 533): under Assumptions 7 and 8, $\mathbb E\|\widehat X_t\|^p<C_p$ and $\mathbb E\|\overline X_t\|^p<C_p$ for all $t\ge0$, with $C_p$ independent of $f$, $g$ (given the constants of Assumption 7) and of $h$.
6. (52) (p. 558): for $p\in(2,p^*]$, with $e_t=\widehat X_t-X_t$ and an explicit function $\widehat L$,

$$\mathbb E\big[e^{\lambda pt/2}\|e_t\|^p\big]\le2\Big(\tfrac{8(p-2)}{p\lambda}\Big)^{p/2-1}\int_0^t\mathbb E\big[\widehat L(\widehat X_s,\overline X_s)^{p/2}\|\widehat X_s-\overline X_s\|^p\big]e^{\lambda ps/2}\,\mathrm ds.$$

The remaining step from (52) to Theorem 6 is a time-uniform bound $\mathbb E\|\widehat X_s-\overline X_s\|^{2p}\le C\delta^p$, obtained as in the finite-time analysis of §6.2.

Significance

Theorem 6 says that the adaptive scheme tracks the true solution with error of order $\delta^{1/2}$ forever, not merely on a fixed window. Combined with ergodicity, it bounds the bias of long-time averages of the scheme uniformly in the averaging horizon, which is what the paper's adaptive multilevel Monte Carlo method for invariant measures (§3) needs. Theorem 5 is of independent use: it shows that the adaptive scheme preserves the time-uniform moment bounds of a dissipative SDE, a property the uniform-step Euler scheme lacks for superlinear drift.

All results here are proved in the source; none is open. To our knowledge none has a machine-checked proof. A complete formalization would give the first verified strong convergence result for an adaptive-step SDE scheme, and the infinite-horizon energy estimates (45)–(47) and the error inequality (52) would be the first verified time-uniform estimates of this kind.

Difficulty

The timestep is random and depends on the state, so the grid times $t_n$ are stopping times rather than fixed numbers, and the number of steps up to time $t$ is random. Standard discrete Grönwall arguments, which iterate over a deterministic grid, do not apply directly; neither does the usual argument that the scheme is an Itô process with adapted, piecewise-constant coefficients, which must be rebuilt for random grids.

The obvious route to a time-uniform bound — a finite-time bound with constant $C_T$ followed by $T\to\infty$ — fails because $C_T$ grows exponentially. Uniformity has to come from an exponential weight ($e^{2\alpha t}$ in (45)–(47), $e^{\lambda pt/2}$ in (52)) that converts dissipation or contractivity into a bound that survives division by the weight, with every constant independent of $t$ and of the timestep function.

Formalization scope

The state space is EuclideanSpace ℝ (Fin m); diffusion values are EuclideanSpace ℝ (Fin m × Fin d), whose norm is the Frobenius norm. Times are in ℝ≥0, grid indices start at $0$, and $x_0$ is deterministic. The grid is defined by recursion; the index $n_t$ is the least $n$ with $t<t_{n+1}$ and equals the junk value $0$ on the event that the grid never passes $t$ (a null event under the hypotheses). Moments are lower Lebesgue integrals ∫⁻ … ‖·‖ₑ ^ p in $[0,\infty]$ with real exponents, never Bochner integrals. A solution on $[0,\infty)$ is a process that solves the SDE on every horizon $[0,T]$, in the sense of the published definition SabanisEuler.Shared.IsSolution. Existential constants are placed before every quantifier they must not depend on: before $t$ in Lemma 3, before $f$, $g$, $h$ and $t$ in Theorem 5, before $f$, $g$, $x_0$, $h$ and $t$ in Step 2, before $\delta$ and $t$ in Theorem 6. A constant chosen after $t$, or a moment that is a Bochner integral of a possibly non-integrable function, would make these statements trivial; both are ruled out. Pages are the journal's printed pages (PDF page + 525).

Two deviations from the printed text are deliberate. Theorem 6 additionally assumes $\|g(x)\|^2\le\beta_g$: its proof invokes Theorem 5, whose Assumption 7 contains this bound, which Assumption 9 does not imply, and the introduction (p. 527) assumes "a bounded and non-degenerate diffusion coefficient". Each $h^\delta$ is assumed measurable. Separately, the printed (47) has $\phi(\widehat X_{\underline t})=\widehat X_{\underline t}+h(\widehat X_{\underline t})f(\widehat X_{\underline t})$ in its last term; the formalized inequality uses $\widehat X_{\underline t}+(t-\underline t)f(\widehat X_{\underline t})$, as the expansion and the source's (48) require. The undefined word "nondegenerate" in Assumption 7 is not encoded.

A full development needs the Itô formula for $e^{ct}\|Y_t\|^p$ applied to Itô processes with random-grid coefficients, the Burkholder–Davis–Gundy inequality, and moment bounds for Brownian increments over stopping-time intervals. These are reusable well beyond this mission. Contributions of any of them, or of proofs of the deterministic inequalities (45) and (47), are welcome.

Selected references

• W. Fang, M. B. Giles, Adaptive Euler–Maruyama method for SDEs with nonglobally Lipschitz drift, Ann. Appl. Probab. 30(2) (2020) 526–560. https://doi.org/10.1214/19-AAP1507
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong and weak divergence in finite time of Euler's method for SDEs with non-globally Lipschitz coefficients, Proc. R. Soc. A 467 (2011) 1563–1576. https://doi.org/10.1098/rspa.2010.0348
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Ann. Appl. Probab. 22(4) (2012) 1611–1641. https://doi.org/10.1214/11-AAP803
• S. Sabanis, Euler approximations with varying coefficients: the case of superlinearly growing diffusion coefficients, Ann. Appl. Probab. 26(4) (2016) 2083–2105. https://doi.org/10.1214/15-AAP1140
• H. Lamba, J. C. Mattingly, A. M. Stuart, An adaptive Euler–Maruyama scheme for SDEs: convergence and stability, IMA J. Numer. Anal. 27(3) (2007) 479–506. https://doi.org/10.1093/imanum/drl032
• D. J. Higham, X. Mao, A. M. Stuart, Strong convergence of Euler-type methods for nonlinear stochastic differential equations, SIAM J. Numer. Anal. 40(3) (2002) 1041–1063. https://doi.org/10.1137/S0036142901389530

Why strong rates for non-monotone SDEs matter

Stochastic differential equations (SDEs) $dX_t=\mu(X_t)\,dt+\sigma(X_t)\,dW_t$ model physical, biological and financial systems whose drift $\mu$ often grows superlinearly: the stochastic Lorenz equation, the stochastic van der Pol and Duffing–van der Pol oscillators, Langevin dynamics. Solutions are rarely explicit, so they are simulated by time-stepping schemes. For multilevel Monte Carlo methods (Giles 2008) the cost of a simulation is governed by the strong convergence rate, the speed at which $\|X_t-Y_t\|_{L^r(\Omega)}$ tends to $0$ as the step size shrinks.

For superlinearly growing coefficients the classical Euler–Maruyama method diverges in $L^p$ (Hutzenthaler, Jentzen, Kloeden 2011). Tamed schemes repair this, and rate $\tfrac12$ was proved for them under a global monotonicity condition: $\langle x-y,\mu(x)-\mu(y)\rangle+\frac{p-1}{2}\|\sigma(x)-\sigma(y)\|^2\le c\|x-y\|^2$ for all $x,y$ (Hutzenthaler, Jentzen, Kloeden 2012; Sabanis 2016). Most of the examples above violate that condition. Before this paper no strong rate was known for any time-discrete scheme for a multidimensional SDE with non-globally monotone coefficients.

Timeline. 2011: divergence of explicit Euler for superlinear drifts. 2012: tamed Euler, rate $\tfrac12$ under one-sided Lipschitz drift and globally Lipschitz diffusion. 2012–2013: strong convergence without rates for non-globally monotone SDEs (stopped and tamed schemes). 2014: this paper (arXiv:1401.0295), rate $\tfrac12$ for the stopped-tamed Euler–Maruyama scheme under the conditions below; published in Ann. Probab. 48 (2020).

Setting

Fix $d,m\in\mathbb N=\{1,2,\dots\}$ and $T\in(0,\infty)$. On a probability space with a normal filtration $(\mathcal F_t)$, $W$ is a standard $m$-dimensional $(\mathcal F_t)$-Brownian motion. Matrices $A\in\mathbb R^{d\times m}$ carry the Hilbert–Schmidt norm $\|A\|_{HS}=(\sum_{ij}A_{ij}^2)^{1/2}$. Given coefficients $\mu:\mathbb R^d\to\mathbb R^d$ and $\sigma:\mathbb R^d\to\mathbb R^{d\times m}$, $X$ is an adapted continuous process with

$$X_t=X_0+\int_0^t\mu(X_s)\,ds+\int_0^t\sigma(X_s)\,dW_s,\qquad t\in[0,T].$$

The taming map is $\psi(v)=v/(1+\|v\|^2)$. The stopped-tamed Euler–Maruyama scheme with $N$ steps is $Z^N_0=X_0$ and

$$Z^N_{n+1}=Z^N_n+\mathbb 1_{\{\|Z^N_n\|<\exp(|\ln(T/N)|^{1/2})\}}\;\psi\Big(\mu(Z^N_n)\tfrac TN+\sigma(Z^N_n)\big(W_{(n+1)T/N}-W_{nT/N}\big)\Big).$$

The whole increment, drift and noise together, is tamed, and the scheme freezes once it leaves a ball whose radius grows slowly as $N\to\infty$.

The conditions use a Lyapunov-type function $U_0\ge1$ and a function $U_1\ge0$. The generator is $(\mathcal G_{\mu,\sigma}\varphi)(x)=\varphi'(x)\mu(x)+\frac12\sum_j\varphi''(x)(\sigma_j(x),\sigma_j(x))$, where $\sigma_j(x)$ is the $j$-th column of $\sigma(x)$. For partitions $\theta=(0=t_0<\dots<t_n=T)$ with mesh $|\theta|=\max_k(t_{k+1}-t_k)$, the same recursion run in continuous time on $[t_k,t_{k+1}]$ defines interpolations $Y^\theta$.

Formalization targets

Goal: Theorem 1.3

Let $c,r\in(0,\infty)$, $q_0,q_1\in(0,\infty]$, $\alpha\ge0$ and $p,q\in[2,\infty)$ with $\frac1p+\frac1{q_0}+\frac1{q_1}=\frac1r$. Let $\mu,\sigma,U_1$ be $C^1$ with polynomially growing derivatives. Let $U_0\in C^3$ satisfy $U_0\ge1$, $\sum_{i\le3}\|U_0^{(i)}\|\le cU_0^{1-1/q}$, $\|x\|^{1/c}\le c(1+U_0(x))$ and $\mathbb E[e^{U_0(X_0)}]<\infty$, and assume for $x\neq y$

$$\frac{\langle x-y,\mu(x)-\mu(y)\rangle+\frac{(p-1)(1+1/c)}{2}\|\sigma(x)-\sigma(y)\|_{HS}^2}{\|x-y\|^2}\le c+\frac{U_0(x)+U_0(y)}{2q_0Te^{\alpha T}}+\frac{U_1(x)+U_1(y)}{2q_1e^{\alpha T}},$$ $$(\mathcal G_{\mu,\sigma}U_0)(x)+\tfrac12\|\sigma(x)^*\nabla U_0(x)\|^2+U_1(x)\le\alpha U_0(x)+c.$$

Then there is $C\in[0,\infty)$ such that for all $N\in\mathbb N$,

$$\sup_{n\in\{0,\dots,N\}}\|X_{nT/N}-Z^N_n\|_{L^r(\Omega)}\le C\,N^{-1/2}.$$

The constant $C$ is not specified. Only the rate is asserted, and that is the stable content of the result.

Milestones

In the proof's order:

• Proposition 2.9: a pathwise Itô inequality for $\|X-Y\|^p$ between a solution and an arbitrary Itô process.
• Theorem 2.10: the perturbation estimate, which bounds $\|X_\tau-Y_\tau\|_{L^r}$ by an exponential moment of the monotonicity ratio times the local errors.
• Corollary 2.12 (= Theorem 1.2): the same estimate with $L^p$ local errors.
• Lemma 3.1: derivative bounds for $\psi$.
• Lemma 3.2: the error estimate on one partition for the scheme stopped at its first grid exit from a set $\mathcal O$.
• Proposition 3.3: $\sup_{t\le T}\|X_t-Y^\theta_t\|_{L^r}\le C|\theta|^{1/2}$ uniformly over all partitions $\theta$. Theorem 1.3 is its uniform-grid case.

Significance

Theorem 1.3 gives strong rate $\tfrac12$, and with it the complexity bounds of multilevel Monte Carlo, for SDEs with non-globally monotone coefficients. The paper's applications include the stochastic Lorenz equation with bounded noise, the stochastic van der Pol and Duffing–van der Pol oscillators, a model from experimental psychology and overdamped Langevin dynamics. The perturbation estimate, Theorem 2.10, is a separate tool. It compares a solution with any Itô process, so it also gives local Lipschitz dependence on the initial value and, in the paper's §3.2, rates for Galerkin approximations of SPDEs.

All results listed here are proved in the paper, and none has been machine-checked so far. This mission produces Lean statements of the goal and of every numbered step of its proof. The development needs Itô calculus for $\|x-y\|^p$, a localization argument, exponential-moment bounds via Lyapunov functions, and calculus for the taming map. Several of these pieces are reusable well beyond this paper.

Difficulty

The obvious argument is Gronwall's lemma applied to $\mathbb E\|X_t-Y_t\|^p$. It needs the global bound $\langle x-y,\mu(x)-\mu(y)\rangle\le c\|x-y\|^2$, which fails for the Lorenz and van der Pol drifts: there the one-sided Lipschitz "constant" grows with $|x|+|y|$. The difficulty is to control a random, state-dependent Gronwall exponent $\int_0^t[\text{ratio}(X_s,Y_s)]^+ds$. That requires exponential integrability of $U_0(X_t)$ and $U_1$ along both the solution and the scheme. For the scheme it holds only because of the taming and the stopping, and only up to an error of the right order.

Formalization scope

• Representation. States are EuclideanSpace ℝ (Fin d). Diffusion values are the published Diffusion d m, whose norm is the Hilbert–Schmidt norm the paper uses throughout. Time is $\mathbb R_{\ge0}$.
• Probability model. Solutions and Itô processes are the published SabanisEuler.Shared.IsSolution / IsItoProcess. The normal filtration is right-continuity plus completion, and $W$ is the published IsWienerMartingale, a Brownian motion on $[0,\infty)$. For $H=\mathbb R^d$, $U=\mathbb R^m$ the paper's cylindrical Wiener process is exactly such a Brownian motion. All processes are hypotheses; existence is never claimed.
• Finite-dimensional specialization. The paper states §2 for separable Hilbert spaces. Propositions 2.9, Theorem 2.10 and Corollary 2.12 are formalized for $H=\mathbb R^d$, $U=\mathbb R^m$, the only case the proof of Theorem 1.3 uses. Proposition 2.9 is further restricted to $\varepsilon\in(0,\infty)$. Proposition 2.5, Theorem 1.4 and the SPDE results of §3.2 are not part of the mission.
• Extended arithmetic. Every norm, moment and exponential moment lies in $[0,\infty]$, with $0\cdot\infty=0$, $0/0=0$, $1/0=\infty$ and $\exp(\infty)=\infty$. The monotonicity ratio is computed in extended reals, so $\varepsilon=\infty$ and $q_0=\infty$ are allowed. Predictable processes are predictable with respect to $(\mathcal F_t)$, and stopping times are taken with respect to the completed filtration.
• Ruling out trivial versions. The exponential-moment hypothesis is a Lebesgue integral in $[0,\infty]$, never a Bochner integral that would default to $0$. Stochastic integrals are pinned down by the published Itô-integral relation, never left as free witnesses. The classes $\mathcal C^1_{\mathcal P}$ (locally Lipschitz with polynomial constant) and $\mathcal C^3_{\mathcal D}$ of Proposition 3.3 are not replaced by $C^1$/$C^3$. A sorry-free check confirms that the goal's hypotheses are jointly satisfiable.
• Contributions welcome. Itô's formula for $C^2$ functions of two Itô processes, Burkholder–Davis–Gundy-type moment bounds for the scheme increments, and exponential-moment estimates from Lyapunov conditions (the paper relies on Cox, Hutzenthaler, Jentzen 2013). Each is reusable beyond this mission.

Selected references

• M. Hutzenthaler, A. Jentzen, On a perturbation theory and on strong convergence rates for stochastic ordinary and partial differential equations with non-globally monotone coefficients, arXiv:1401.0295v1, 2014; Ann. Probab. 48(1), 2020. https://arxiv.org/abs/1401.0295
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Ann. Appl. Probab. 22(4), 2012. https://doi.org/10.1214/11-AAP803
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong and weak divergence in finite time of Euler's method for SDEs with non-globally Lipschitz coefficients, Proc. R. Soc. A 467, 2011. https://doi.org/10.1098/rspa.2010.0348
• S. Sabanis, Euler approximations with varying coefficients: the case of superlinearly growing diffusion coefficients, Ann. Appl. Probab. 26(4), 2016. https://doi.org/10.1214/15-AAP1140
• S. Cox, M. Hutzenthaler, A. Jentzen, Local Lipschitz continuity in the initial value and strong completeness for nonlinear SDEs, arXiv:1309.5595, 2013. https://arxiv.org/abs/1309.5595
• M. B. Giles, Multilevel Monte Carlo path simulation, Oper. Res. 56(3), 2008. https://doi.org/10.1287/opre.1070.0496

Motivation

Many equilibrium problems in operations research, such as networked Cournot competition, power markets and communication networks, are stochastic Nash games: each player minimizes an expected cost $f_i(x_i,x_{-i})=\mathbb E[\psi_i(x_i,x_{-i};\xi)]$ that depends on the rivals' strategies and on a random vector $\xi$ whose law is accessible only through samples. Best-response schemes, in which every player in turn solves its own optimization problem given the others' latest strategies, are the natural distributed method for such games. In the stochastic setting no player can compute an exact best response, because evaluating $f_i$ requires the expectation itself. Lei, Shanbhag, Pang and Sen (arXiv:1704.04578v2; Mathematics of Operations Research, 2020) analyze inexact proximal best-response schemes in which each best response is approximated by a finite number of stochastic gradient steps, and they ask how many sampled gradients each player needs to reach an approximate equilibrium.

The deterministic theory of proximal best responses, including the contraction matrix $\Gamma$ used here, is due to Facchinei and Pang (Nash equilibria: the variational approach, 2009, §12.6; reference [19] of the paper). The present mission covers the paper's first scheme, the synchronous one, in which all players update simultaneously.

Setting

There are $N$ players. Player $i$ chooses $x_i$ from a closed, compact, convex, nonempty set $X_i\subseteq\mathbb R^{n_i}$; a profile is $x=(x_1,\dots,x_N)\in X=\prod_iX_i$, and $(z,y_{-i})$ is the profile $y$ with its $i$-th block replaced by $z$. A Nash equilibrium is a profile $x^*\in X$ such that every $x^*_i$ minimizes $f_i(\cdot,x^*_{-i})$ over $X_i$.

Assumption 1 asks that each $f_i$ be convex in $x_i$ and twice continuously differentiable near $X$, that the sampled gradients $\nabla_{x_i}\psi_i(x;\xi)$ be unbiased for $\nabla_{x_i}f_i(x)$, and that $\mathbb E\|\nabla_{x_i}\psi_i(x;\xi)\|^2\le M_i^2$. For $\mu>0$ the proximal best response is

$$\hat x_i(y)=\operatorname*{argmin}_{x_i\in X_i}\Big[f_i(x_i,y_{-i})+\tfrac\mu2\|x_i-y_i\|^2\Big].$$

The $N\times N$ matrix $\Gamma$ has entries $\gamma_{ii}=\mu/(\mu+\zeta_{i,\min})$ and $\gamma_{ij}=\zeta_{ij,\max}/(\mu+\zeta_{i,\min})$, where $\zeta_{i,\min}$ is the least eigenvalue of $\nabla^2_{x_i}f_i$ over $X$ and $\zeta_{ij,\max}$ the largest spectral norm of the mixed block $\nabla^2_{x_ix_j}f_i$ over $X$. Assumption 2 is $a:=\|\Gamma\|<1$ (spectral norm).

Algorithm 1 starts from a deterministic $x_0\in X$. At major iteration $k$ every player computes $x_{i,k+1}\in X_i$ with $\mathbb E[\|x_{i,k+1}-\hat x_i(x_k)\|^2\mid\mathcal F_k]\le\alpha_{i,k}^2$, by running $j_{i,k}$ projected stochastic gradient steps

$$z_{i,t+1}=\Pi_{X_i}\Big[z_{i,t}-\tfrac{1}{\mu(t+1)}\big(\nabla_{x_i}\psi_i(z_{i,t},x_{-i,k};\xi^t_{i,k})+\mu(z_{i,t}-x_{i,k})\big)\Big],\qquad z_{i,1}=x_{i,k},$$

and setting $x_{i,k+1}=z_{i,j_{i,k}}$. A random profile $x$ is an $\epsilon$-NE$_2$ if $\mathbb E\big[(\sum_i\|x_i-x^*_i\|^2)^{1/2}\big]\le\epsilon$.

Formalization targets

Goal: Theorem 1(a)

With $\alpha_{i,k}=\eta^{k+1}$, $\eta\in(0,1)$, $\|x_{i,0}-x^*_i\|\le C$, $c=\max\{a,\eta\}$, $q\in(c,1)$, $D=1/\ln((q/c)^e)$, $Q_i=2M_i^2/\mu^2+2D_{X_i}^2$ and $j_{i,k}=\lceil Q_i/\eta^{2(k+1)}\rceil$, the iterate $x_K$ after $K=\lceil\ln(\sqrt N(C+D)/\epsilon)/\ln(1/q)\rceil$ major iterations is an $\epsilon$-NE$_2$, and player $i$ uses at most

$$\ell_i(\eta)=\frac{Q_i}{\eta^4\ln(1/\eta^2)}\left(\frac{\sqrt N(C+D)}{\epsilon}\right)^{\frac{\ln(1/\eta^2)}{\ln(1/q)}}+\left\lceil\frac{\ln(\sqrt N(C+D)/\epsilon)}{\ln(1/q)}\right\rceil$$

projected stochastic gradient steps. The goal fixes the paper's explicit expression, not its $O(\cdot)$ summary.

Milestones

1. (5): $\|\hat x_i(y')-\hat x_i(y)\|\le\sum_j\gamma_{ij}\|y'_j-y_j\|$ for $y,y'\in X$.
2. (9): the same in the Euclidean norm of block distances, with factor $a$.
3. A Nash equilibrium is a fixed point: $\hat x(x^*)=x^*$.
4. Lemma 2: $zc^z\le Dq^z$ for $z\ge0$, $0<c<q<1$, $D\ge1/(e\ln(q/c))$.
5. (20): $u_{k+1}\le au_k+\sqrt N\max_i\alpha_{i,k}$, where $u_k$ is the mean Euclidean norm of block errors.
6. Proposition 4: $u_k\le\sqrt N(C+D)q^k$.
7. Lemma 3: $\mathbb E[\|z_{i,t}-\hat x_i(x_k)\|^2\mid\mathcal F_k]\le Q_i/(t+1)$.
8. (30): $\sum_{k=1}^K\beta^k\le\int_1^{K+1}\beta^xdx\le\beta^{K+1}/\ln\beta$ for $\beta>1$.

Significance

Theorem 1(a) shows that the synchronous inexact scheme keeps the linear rate of the exact best-response iteration, with factor arbitrarily close to $\max\{\|\Gamma\|,\eta\}$, and turns it into a sample complexity explicit in $\epsilon$, in the number of players and in the game's constants. Choosing $\eta=\|\Gamma\|$ gives a bound of order $(\sqrt N/\epsilon)^{2+\delta}$ for any $\delta>0$ (Theorem 1(b)), close to the $O(1/\epsilon^2)$ of a single stochastic convex program; the dependence on $N$ is mild, which matters for games with many players. The same building blocks (the contraction (5), Lemma 2, Lemma 3 and (30)) are reused by the randomized and asynchronous schemes of the same paper.

The result is proved in the paper; it has not been machine-checked. A formalization adds a verified proof of the complexity bound together with fully explicit hypotheses: several conditions that the paper leaves implicit or states imprecisely (listed under Formalization scope) are made exact, and the Lean statement records which version of each is used.

Difficulty

The arithmetic of Theorem 1(a) from Proposition 4, Lemma 3 and (30) is routine. The weight lies elsewhere. The contraction (5) is quoted from Facchinei–Pang and not proved in the paper; it needs the variational optimality conditions of the two proximal problems and a mean-value argument along a segment of $X$ for the partial gradient, which mixes the Hessian blocks $\nabla^2_{x_i}f_i$ and $\nabla^2_{x_ix_j}f_i$. Lemma 3 is a conditional-expectation argument about a stochastic approximation recursion: the iterate $z_{i,t}$ must be shown measurable for the nested $\sigma$-fields generated by the samples, and the conditional unbiasedness and moment hypotheses must be combined with the tower property and the projection's nonexpansiveness. Proposition 4 then needs conditional Jensen to pass from (8) to the Euclidean norm of block errors.

Formalization scope

Players are Fin N, player $i$'s space is EuclideanSpace ℝ (Fin (n i)), profiles are dependent functions, and $(z,y_{-i})$ is Function.update y i z. The Euclidean norm of the vector of block distances is written out (blockDist), since Mathlib's norm on profiles is the sup norm; $\|\Gamma\|$ is the operator norm of Matrix.toEuclideanCLM Γ. $\zeta_{i,\min}$ and $\zeta_{ij,\max}$ are the infimum and supremum of Rayleigh-type quotients of the second derivative over $X$ and unit vectors. The proximal BR is any map minimizing the proximal objective on $X$ (unique there), and $\Pi_{X_i}$ is any map satisfying the published SpectralProjGrad.Shared.IsProjOnto. The $\sigma$-fields $\mathcal F_k$ and $\sigma\{\mathcal F_k,\xi^{[t-1]}_{i,k}\}$ are generated by the samples. Every expectation carries an integrability hypothesis.

Statement repairs, each recorded in the item's Formalization Note:

• the initial point is deterministic, as Algorithm 1 states, so $\mathbb E\|x_{i,0}-x^*_i\|\le C$ reads $\|x_{i,0}-x^*_i\|\le C$ (the proof's $u_0\le\sqrt NC$ needs this);
• Assumption 1(b) is twice continuous differentiability in the whole profile, which (4) needs for the mixed blocks, and convexity on an open set containing $X_i$;
• in Lemma 3, $\psi_i$ is read as $\nabla_{x_i}\psi_i$, and the second hypothesis is the bound $\mathbb E[\|\nabla_{x_i}\psi_i\|^2\mid\cdot]\le M_i^2$ that the proof uses (the printed equality would force zero sampling noise);
• $0<\epsilon<\sqrt N(C+D)$, without which the ceiling in (27) can be non-positive;
• $D=1/(e\ln(q/c))$, the reading of $\ln((q/c)^e)$ confirmed by Lemma 2's proof.

The goal is not trivialized by deterministic gradients: the sampled gradients are random, constrained only by conditional unbiasedness and a conditional second-moment bound, the proximal BR is tied to $f_i$ by its argmin property, $x^*$ must be a Nash equilibrium, and the norms are Euclidean. A deterministic instance ($N=1$, $f_1(x)=\|x\|^2$ on a ball, exact gradients) satisfies every hypothesis jointly.

Needed infrastructure: optimality conditions for constrained strongly convex minimization, the mean-value inequality for partial gradients of $C^2$ functions on convex sets, conditional expectation of vector-valued functions with the tower property and conditional Jensen, and the nonexpansiveness of Euclidean projection. The contraction (5) and Lemma 3 are reusable beyond this mission. Proofs of any milestone are welcome.

Selected references

• J. Lei, U. V. Shanbhag, J.-S. Pang, S. Sen, On Synchronous, Asynchronous, and Randomized Best-Response Schemes for Stochastic Nash Games, arXiv:1704.04578v2, 2018; published in Mathematics of Operations Research, 2020. https://arxiv.org/abs/1704.04578v2
• F. Facchinei, J.-S. Pang, Nash equilibria: the variational approach, in Convex Optimization in Signal Processing and Communications, Cambridge University Press, 2009 (reference [19] of arXiv:1704.04578v2).
• B. T. Polyak, Introduction to Optimization, Optimization Software, 1987 (reference [38] of the paper).

Why adaptive timesteps

The explicit Euler–Maruyama scheme is the standard way to simulate a stochastic differential equation

$$dX_t=f(X_t)\,dt+g(X_t)\,dW_t,\qquad X_0=x_0,$$

on $\mathbb R^m$ driven by a $d$-dimensional Brownian motion. When $f$ and $g$ are globally Lipschitz it converges strongly with order $\tfrac12$ in the step size (Kloeden and Platen 1992). Many models of interest in physics, chemistry and statistics have a drift that grows faster than linearly but points inwards, such as $f(x)=-x^3$. For these the SDE is well posed and its solution has finite moments of all orders, but the uniform-step scheme does not: for $dX_t=-X_t^3dt+dW_t$, Hutzenthaler, Jentzen and Kloeden proved that $\mathbb E\|X_T\|^p\to\infty$ as the step size $h\to0$, for every $T>0$ and $p\ge2$ (Proc. R. Soc. A 2011).

Several remedies have been proposed: implicit schemes (Higham, Mao and Stuart 2002; Mao and Szpruch 2013), the tamed Euler scheme (Hutzenthaler, Jentzen and Kloeden 2012), truncated schemes (Mao 2015). Fang and Giles (Ann. Appl. Probab. 2020) keep the explicit Euler–Maruyama update and change only the step: the timestep is a function of the current state, small where the drift is large. This mission formalizes their finite-time results: stability of the adaptive scheme (Theorem 1) and strong convergence of order $\tfrac12$ (Theorem 3).

The adaptive scheme

Let $\|x\|$ and $\langle x,y\rangle$ be the Euclidean norm and inner product on $\mathbb R^m$, and $\|A\|=(\sum_{i,j}A_{ij}^2)^{1/2}$ the Frobenius norm of an $m\times d$ matrix. Given a timestep function $h:\mathbb R^m\to(0,\infty)$, the scheme (5) starts from $t_0=0$, $\widehat X_0=x_0$ and sets $h_n=h(\widehat X_{t_n})$,

$$t_{n+1}=t_n+h_n,\qquad \widehat X_{t_{n+1}}=\widehat X_{t_n}+f(\widehat X_{t_n})h_n+g(\widehat X_{t_n})(W_{t_{n+1}}-W_{t_n}).$$

The grid is random: it depends on the path. For a time $t$, $\underline t$ is the last grid time not after $t$; the piecewise constant interpolant is $\overline X_t=\widehat X_{\underline t}$ and the continuous interpolant is

$$\widehat X_t=\widehat X_{\underline t}+f(\widehat X_{\underline t})(t-\underline t)+g(\widehat X_{\underline t})(W_t-W_{\underline t}).$$

The hypotheses (Assumptions 1–4, pp. 528–530):

1. Assumption 1: $f,g$ locally Lipschitz, $\langle x,f(x)\rangle\le\alpha_1\|x\|^2+\beta_1$ and $\|g(x)\|^2\le\alpha_1\|x\|^2+\beta_1$.
2. Assumption 2: $h$ continuous and positive with $\langle x,f(x)\rangle+\tfrac12h(x)\|f(x)\|^2\le\alpha\|x\|^2+\beta$ for some $\alpha,\beta>0$.
3. Assumption 3: a family $h^\delta$, $0<\delta\le1$, with $\delta\min(T,h(x))\le h^\delta(x)\le\min(\delta T,h(x))$; $\delta\to0$ plays the role of $h\to0$.
4. Assumption 4: $\langle x-y,f(x)-f(y)\rangle\le\tfrac12\alpha\|x-y\|^2$, $\|g(x)-g(y)\|^2\le\tfrac12\alpha\|x-y\|^2$, and $\|f(x)-f(y)\|\le(\gamma(\|x\|^q+\|y\|^q)+\mu)\|x-y\|$.

For $f(x)=-x^3$ in one dimension, $h(x)=\min(1,1/x^2)$ satisfies Assumption 2.

Formalization targets

Goal: Theorem 3 (strong convergence order, p. 531)

Under Assumption 4, with $h^\delta$ satisfying Assumption 3 for an $h$ satisfying Assumption 2, for every $p>0$ there is $C_{p,T}$ such that for all $\delta\in(0,1]$

$$\mathbb E\Big[\sup_{0\le t\le T}\|\widehat X_t-X_t\|^p\Big]\le C_{p,T}\,\delta^{p/2}.$$

Milestones

• Lemma 1 (p. 529): under Assumption 1, $\mathbb E[\sup_{0\le t\le T}\|X_t\|^p]<\infty$ for all $p>0$.
• (35), (37) (p. 548): one-step and partial-step energy inequalities for the projected K-scheme $\widehat X^K_{t_{n+1}}=P_K(\widehat X^K_{t_n}+f(\widehat X^K_{t_n})h_n+g(\widehat X^K_{t_n})\Delta W_n)$, $P_K(Y)=\min(1,K/\|Y\|)Y$.
• §6.1 Step 3 (p. 550): $\mathbb E[\sup_{0\le t\le T}\|\widehat X^K_t\|^p]\le C_{p,T}$ uniformly in $K>\|x_0\|$.
• Theorem 1 (p. 529): $T$ is almost surely attained and $\mathbb E[\sup_{0\le t\le T}\|\widehat X_t\|^p]<C_{p,T}$, with $C_{p,T}$ independent of $h$.
• §6.2, after (41) (p. 553): $\mathbb E\|\widehat X_s-\overline X_s\|^{2p}\le C^3_{p,T}\delta^p$.
• (40) (p. 552): the integral inequality for $e_t=\widehat X_t-X_t$ to which Grönwall's inequality is applied.

Significance

Theorem 1 shows that an explicit scheme, at the cost of one drift and one diffusion evaluation per step, is stable for one-sided Lipschitz drifts of polynomial growth, the class where the uniform-step scheme fails. Theorem 3 shows that refining the timestep function keeps the classical order $\tfrac12$. Combined with the paper's bound on the expected number of steps (Lemma 2) this gives order $\tfrac12$ in terms of computational cost, and the results underpin the paper's adaptive multilevel Monte Carlo estimators. The uniform-in-$h$ form of Theorem 1 is what lets Theorem 3 apply the stability bound to the whole family $h^\delta$.

The paper's proofs are complete on paper; none of these results has a machine-checked proof. Formalizing them needs a theory of Euler schemes on a random, state-dependent grid, which is not in Mathlib and is not covered by uniform-grid formalizations of tamed or implicit schemes. The definitions here (the scheme on random times, the projected K-scheme) are reusable for the infinite-time companion mission and for other adaptive schemes.

Difficulty

The obvious argument for the uniform-step scheme bounds moments step by step, using deterministic step sizes and independent Gaussian increments. Here the step $h_n$ depends on $\widehat X_{t_n}$, so the grid times are stopping times and the increments $W_{t_{n+1}}-W_{t_n}$ have random length. Sums over steps become Itô integrals only after the scheme is written as the solution of $d\widehat X_t=f(\overline X_t)dt+g(\overline X_t)dW_t$, and moment bounds for those integrals need the Burkholder–Davis–Gundy inequality on a random grid. A second difficulty is that $T$ need not be reached: if $h(\widehat X_{t_n})$ decays fast enough, $t_n$ could converge below $T$. Excluding this needs a positive lower bound on $h$ over bounded sets together with moment bounds, which the paper obtains by first analysing a version of the scheme projected onto a ball.

Formalization scope

Lean conventions: $\mathbb R^m$ is the published SDEState m (EuclideanSpace ℝ (Fin m)); matrices are the published Diffusion m d, whose norm is the Frobenius norm; $W$ is a standard $d$-dimensional $(\mathcal F_t)$-Brownian motion in the sense of the published IsWienerMartingale; the SDE solution is any process satisfying the published IsSolution on $[0,T]$ (existence and uniqueness are not claimed). Times are in $[0,\infty)$, indices start at $0$, and $x_0$ is deterministic. The scheme is defined pathwise; the index of $\underline t$ is a junk value $0$ when the grid never passes $t$, an event of probability zero (Theorem 1). $h$ and $h^\delta$ are real-valued and their positivity is a hypothesis; each $h^\delta$ is assumed measurable. Moments are $[0,\infty]$-valued integrals of extended norms, never Bochner integrals. Every existential constant is chosen before the parameter it must not depend on: before $\delta$ in Theorem 3 and its milestones, before $h$ (and $K$) in Theorem 1 and Step 3. Results the paper proves for $p\ge4$ inside its proofs are stated for $p\ge4$. Pages are the journal's printed pages (PDF page + 525).

A formalization that places the constant after $\delta$, uses a deterministic grid, or allows $h=0$ would be trivially true or a different theorem; these are ruled out by the statements.

Contributions welcome: proofs of the milestones, and of general facts they need (Burkholder–Davis–Gundy, moments of Brownian increments over stopping-time intervals, Grönwall's inequality for $[0,\infty]$-valued functions).

Selected references

• W. Fang, M. B. Giles, Adaptive Euler–Maruyama method for SDEs with nonglobally Lipschitz drift, Ann. Appl. Probab. 30(2), 526–560, 2020. https://doi.org/10.1214/19-AAP1507 (preprint arXiv:1609.08101)
• D. J. Higham, X. Mao, A. M. Stuart, Strong convergence of Euler-type methods for nonlinear stochastic differential equations, SIAM J. Numer. Anal. 40, 1041–1063, 2002. https://doi.org/10.1137/S0036142901389530
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients, Proc. R. Soc. A 467, 1563–1576, 2011. https://doi.org/10.1098/rspa.2010.0348
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Ann. Appl. Probab. 22, 1611–1641, 2012. https://doi.org/10.1214/11-AAP803
• X. Mao, L. Szpruch, Strong convergence and stability of implicit numerical methods for stochastic differential equations with non-globally Lipschitz continuous coefficients, J. Comput. Appl. Math. 238, 14–28, 2013. https://doi.org/10.1016/j.cam.2012.08.015
• X. Mao, The truncated Euler–Maruyama method for stochastic differential equations, J. Comput. Appl. Math. 290, 370–384, 2015. https://doi.org/10.1016/j.cam.2015.06.002
• P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, 1992. https://doi.org/10.1007/978-3-662-12616-5