PaperlayerView exact revision on arXiv

Source-linked AI summary

Mean Square Error Analysis of Stochastic Runge-Kutta Integrators

Xuda Ye

arXiv:2609.11528v4math.NAcs.NA

TL;DR

Higher-order Langevin integrators lacked finite-step mean square error bounds for time averages under conditions that do not require global convexity. This paper analyzes two stochastic Runge–Kutta schemes using a local-error decomposition and discrete Poisson equation, obtaining second-order MSE and Wasserstein-1 guarantees.

  • Problem

    For overdamped Langevin dynamics, prior work had not established a second-order finite-N MSE bound under conditions on the potential alone without global convexity.

    Source: §1 · p. 2, §1.2 · p. 6

  • Method

    The analysis splits each integrator’s local error into a mean-zero term and a smaller remainder, then uses the discrete Poisson equation to bound time-averaged error.

    Source: §4.1 · p. 13

  • Results

    The MSE is O(1/(Nh) + h^4), and the chain’s law approaches the exact solution’s law in Wasserstein-1 at second order up to a logarithm, uniformly in the number of steps.

    Source: §4.2 · pp. 14–15, §6 · p. 22

  • Takeaways & Limitations

    The expected second-order MSE rates hold at every finite N without global convexity, for integrators that evaluate ∇U alone.

    Source: §6 · p. 22

    The potential is assumed convex only outside a bounded region, and the dependence of the bounds on dimension is not tracked.

    Source: §1.1 · p. 2, §4.2 · p. 15

Abstract

from arXiv · show

We analyze the mean square error of stochastic Runge-Kutta integrators for overdamped Langevin dynamics whose potential is strongly convex outside a bounded region. A decomposition splits the local error into a mean-zero term and a smaller remainder, and the discrete Poisson equation turns their moments into a bound on the error of a time average. We carry this out for two stochastic Runge-Kutta integrators proposed by Yang & Wang (2026), both of strong order $\frac32$ and weak order $2$, which evaluate the gradient of the potential and no higher derivative. We show that their laws approach the law of the exact solution at second order in Wasserstein-1 distance, up to a logarithm, uniformly in the number of steps. For a test function with bounded derivatives up to third order, we prove that the mean square error over $N$ steps with step size $h$ is $\mathcal O ( \frac{1}{Nh} + h^4 )$, which is the optimal order in the discretization. Experiments measure the strong and weak orders and the sampling bias on a nonconvex potential in $\mathbb R^2$, and compare the integrators on a diffusion model of CIFAR-10.

1. Introduction

The introduction frames a gap in finite-step mean square error guarantees for higher-order Langevin integrators and presents an analysis of two stochastic Runge–Kutta schemes. The paper uses a local-error decomposition and discrete Poisson equation to obtain second-order bounds without global convexity.

1. Introduction

  • A non-asymptotic mean square error bound for a single trajectory was lacking for higher-order integrators, despite existing invariant-measure and Wasserstein analyses.The gap is whether these methods provide a finite-step bound on the error of a time average.Source: §1 · p. 2

1.1. Contributions

  • The analysis covers SRK-I and SRK-II, two stochastic Runge–Kutta schemes, and establishes a second-order MSE bound at every finite step count without global convexity or higher potential derivatives.The schemes evaluate ∇U alone; the potential need only be convex outside a bounded region.Source: §1.1 · p. 2
  • A local-error split into a mean-zero term and a smaller remainder, combined with the discrete Poisson equation, transfers one-step moment estimates to time-averaged error bounds.The Poisson solution requires derivative bounds up to third order, uniform in step size.Source: §1.1 · pp. 2–3

1.2. Related work

  • For second-order MSE bounds, weak-error and coupling methods bound test-function expectations at each step but cannot provide the variance of a time average along one trajectory; contractivity or the Poisson equation is needed.Source: §1.1 · p. 2

1.1. Contributions

  • The discrete Poisson equation extends to potentials convex only outside a bounded region; exponential time decay of the Kolmogorov solution and its derivatives supplies the needed estimates without requiring scheme contraction.The decay is established on expectations using reflection coupling and a derivative formula.Source: §1.1 · p. 3

1.2. Related work

  • For overdamped Langevin dynamics, prior results had not combined assumptions on U alone, finite-N MSE bounds, and second-order accuracy; earlier works achieved at most two of these three properties.The paper’s stated gap is specific to the overdamped dynamics.Source: §1.2 · p. 6

2. Preliminaries

The preliminaries set assumptions for the Langevin dynamics, introduce two gradient-only stochastic Runge–Kutta schemes, and establish regularity estimates for the discrete Poisson solution.

2.1. Notation and assumptions

  • Assumption 1 permits nonconvexity inside |x|≤R while requiring convexity outside; its ∇5U bound serves the integrator analysis.Source: §2.1 · p. 8

2.2. Stochastic Runge–Kutta integrators

  • SRK-I uses two gradient evaluations per step and SRK-II uses three; both have strong order 3/2 and weak order 2.Source: §2.2 · p. 8

2.3. Elliptic estimates of discrete Poisson equation

  • The discrete Poisson solution sums the Kolmogorov solution over steps; its exponential decay ensures convergence and yields derivative bounds through third order, uniform in step size.Source: §2.3 · p. 9

3. Local error decomposition

The section splits one-step error into a mean-zero term and a smaller remainder, then uses this structure and a related product estimate to establish finite-time convergence orders.

3.2. Local error moment estimates

  • The local error is split into a mean-zero term and a remainder; the product error receives the same treatment, with its estimate essential to second-order weak analysis.The numerical and exact steps share the same Brownian increment.Source: §3 · p. 10, §3.2 · p. 11

3.1. Heuristic Taylor expansion

  • Expanding beyond a direct bound sharpens the time-averaged error: the product term has conditional mean of order h^3 after its order-h^(5/2) mean-zero part is split off.Source: §3.1 · p. 10

3.3. Finite-time error estimates

  • The integrators have strong order 3/2 and finite-horizon weak order 2; the weak error compares the laws and does not depend on coupling.Source: §3.3 · p. 12
  • The finite-horizon weak error is second order because one step reproduces the exact step to third order, and summing the local bound gives weak order 2.Source: §3.3 · p. 12

4. Mean square error analysis

The analysis combines one-step consistency and chain stability through a time-averaged error formula. For SRK-I and SRK-II, this yields a finite-step mean square error bound and a uniform-in-step Wasserstein convergence result.

4. Mean square error analysis

  • The analysis combines one-step consistency, chain stability, and the time-averaged error formula; stability uses uniform-in-time chain moment bounds.The constants depend on m, M, R, and d, and the stated step-size restriction is tied to those moment bounds.Source: §4 · p. 12, §4.2 · p. 15

4.1. Time-averaged error formula

  • The time-averaged error formula uses the discrete Poisson solution to express the trajectory average through a boundary term, discretization terms, and martingale differences.The identity holds pathwise, and the Poisson-solution differences telescope when summed over steps.Source: §4.1 · p. 13
  • The local decomposition makes the martingale component mean zero and order h, while its remainder is order h^2 and accumulates linearly in N.The martingale-difference sums have additive variance; the remainder terms accumulate linearly.Source: §4.1 · p. 14

4.2. Main results

  • MSE(N,h) is bounded by a multiple of 1/(Nh) + h^4 for SRK-I and SRK-II under Assumptions 1 and 2.The 1/(Nh) term is sampling variance, while h^4 is the squared O(h^2) discretization bias.Source: §4.2 · pp. 14–15, §4.1 · p. 14
  • At fixed N, optimizing the bound gives h of order N^-1/5 and MSE O(N^-4/5); the result requires bounded derivatives of f through order three.The passage contrasts this with weak order 1, whose optimized step size and MSE are N^-1/3 and O(N^-2/3), respectively.Source: §4.2 · p. 15
  • The laws of the numerical chain approach the exact-solution law at second order up to a logarithm, uniformly in the number of steps.The finite-horizon weak-error bound has a constant containing e^CT; decay of propagated-test-function derivatives removes that factor for the uniform result.Source: §4.2 · p. 15

5. Numerical experiments

The experiments compare gradient-only integrators on a nonconvex planar Langevin problem and a CIFAR-10 diffusion model. SRK-I and SRK-II perform best in finite-time tests, while Leimkuhler–Matthews has the smallest long-time sampling error.

5. Numerical experiments

  • The planar benchmark uses a potential with three Gaussian wells that is nonconvex between wells but convex outside a bounded region; weak errors use 20 Lipschitz test functions.Finite-time comparisons use a fine-step SRK-LD reference on the same Brownian path, while the benchmarks also measure long-time time-average bias and CIFAR-10 image quality.Source: §5 · pp. 15–16, §5.1 · p. 16

5.1. Finite-time strong and weak error

  • SRK-I and SRK-II show measured strong orders 1.56 and 1.55 and weak orders 1.82 and 1.91; at equal cost, SRK-I is more accurate in strong error.The stated bound has slope 1.80, between the measured weak orders and order 2.00.Source: Figure 1 · p. 16, §5.1 · p. 17

5.2. Uniform-in-time sampling error

  • For SRK-I and SRK-II, the fitted long-time sampling-bias orders are 1.96 and 2.01; the coarse curves settle early at the step-size bias, while the finest settle only late.Leimkuhler–Matthews has the smallest sampling error at every cost, while its two finest curves approach the single-trajectory sampling floor.Source: §5.2 · p. 17, Figure 2 · p. 18
  • At equal step size, stochastic Heun's sampling-bias level is about 1.3 times SRK-I's, while random splitting's is several times SRK-II's.The second-order schemes differ in their long-time bias by constants; Leimkuhler–Matthews is about one fifth of SRK-I at equal step size.Source: §5.2 · pp. 17, 19

5.3. Sampling a diffusion model

  • At 60 diffusion steps, SRK-I's image error is 0.021 versus 0.172 for random splitting; stochastic Heun and SRK-II are similar, and random splitting performs worst.The ranking follows the strong error, which governs the quality of a single generated image.Source: §5.3 · p. 19, Figure 4 · p. 21

5.4. Comparison of integrators

  • SRK-I and SRK-II lead in finite-time tests and diffusion-image quality, while Leimkuhler–Matthews has the smallest long-time sampling error, an advantage confined to long-time averages.At equal cost, Leimkuhler–Matthews has an order-of-magnitude lower sampling error than both SRK schemes; SRK-I generates the closest reference images.Source: §5.4 · pp. 20–21, Figure 3 · p. 20, Figure 4 · p. 21

6. Conclusion

The paper develops a local-error route to mean-square bounds for time averages and applies it to SRK-I and SRK-II. It establishes second-order results under convexity only outside a bounded region and reports experiments consistent with those rates.

6. Conclusion

  • The local-error decomposition separates a mean-zero term from a remainder, and the time-averaged error formula converts their moments into a mean-square bound using derivative bounds from Theorem 1.The results for SRK-I and SRK-II include an MSE bound of 1/(Nh) + h^4, finite-horizon weak order 2, and a uniform-in-step Wasserstein-1 bound of order h^2 up to a logarithm.Source: §6 · p. 22
  • These rates hold without global convexity, at finite N, and for integrators evaluating only ∇U; the elliptic estimates bound Poisson-solution derivatives through third order when the potential is convex only outside a bounded set.The experiments confirm the Wasserstein-1 and sampling-bias second orders and rank four diffusion-model schemes consistently with their strong errors.Source: §6 · p. 22
  • Higher weak order would require more derivatives than Theorem 1 supplies; stochastic gradients add first-order bias proportional to their variance, and extending the analysis to underdamped dynamics requires hypoelliptic estimates.A mean-square bound for Leimkuhler–Matthews is another open direction: its second-order bias is long-time, whereas its finite-horizon weak order is one.Source: §6 · p. 22

Data availability

The paper provides code and symbolic coefficient verification, then proves derivative estimates for the discrete Poisson equation under convexity outside a bounded region. The appendix combines reflection coupling with short-time smoothing to control derivatives up to third order.

Data availability

  • The source code is available in the paper's repository, which also contains symbolic verification of the coefficients in Lemmas 1 and 2 from Brownian-motion covariance.The diffusion model is trained on the public CIFAR-10 dataset.Source: Data availability · p. 23

A. Proof of Theorem 1

  • Theorem 1's proof obtains exponential decay of the Kolmogorov solution and its derivatives through third order, then sums those estimates over time steps to bound the discrete Poisson solution.The proof is organized around reflection coupling, short-time derivative bounds, and summation over t = kh.Source: A. Proof of Theorem 1 · p. 25

A.1. Gradient decay and pointwise bound on u

  • Reflection coupling yields gradient decay despite nonconvexity inside a bounded region, but the rate and prefactor are weaker than in the globally convex case.The first variation can grow while trajectories spend time in the nonconvex region, so the proof does not retain the global-convexity rate and unit prefactor.Source: A. Proof of Theorem 1 · p. 25, §A.1 · p. 26

A.2. Short-time derivative bounds

  • Short-time smoothing converts a Lipschitz bound into bounds on the second and third derivatives of the semigroup, which combine with gradient decay to bound the corresponding derivatives of the Poisson solution.The estimates apply to functions with bounded derivatives through third order, and the stated lemma's constant depends on M and dimension.Source: §A.2 · pp. 27, 30

A.3. Conclusion of proof

The proof combines semigroup smoothing and geometric decay to control derivatives of the discrete Poisson solution, while short-time singularities retain a third-derivative requirement. The appendices set up and estimate local-error and product decompositions for the integrators.

A.3. Conclusion of proof

  • For t≥1, semigroup splitting applies unit-time smoothing after gradient decay, yielding the derivative bounds needed to prove (2.8).The discrete Poisson solution is bounded by summing the Kolmogorov-solution bounds over t=kh.Source: §A.3 · pp. 32–33

A.4. Short-time smoothing of semigroup

  • Short-time smoothing uses only a gradient bound, but its third-derivative singularity is nonintegrable and leaves the bounded-third-derivative assumption in place.The second-derivative singularity is integrable, whereas the third-derivative sum grows like log(1/h).Source: §A.4 · p. 34

B. Proofs for Section 3

  • The appendix derives Section 3’s decompositions and finite-time orders for SRK-I, with SRK-II differing through stage displacements and scalars c1 and c3.Remainders are recorded for both integrators.Source: B. Proofs for Section 3 · p. 34

B.1. Notation

  • The notation treats each one-step quantity as measurable from the current state and that step’s Brownian increments, which are independent of the past.This conditional setup applies in particular to the local error and its mean-zero and remainder terms.Source: §B.1 · pp. 34–35

B.2. Explicit remainders

  • The explicit remainders include a fourth-order Taylor remainder of ∇U involving ∇⁵U, alongside remainders for the exact-step gradient and displacement.The product decomposition also uses the local error multiplied by the numerical-versus-exact displacement sum.Source: §B.2 · pp. 36–37

B.3. Proof of Lemma 1: local error decomposition

The local-error proof expands the gradient around the initial state and separates Brownian-driven mean-zero terms from remainders. It establishes this decomposition for both SRK integrators.

B.3. Proof of Lemma 1: local error decomposition

  • For SRK-I, collecting the Gaussian terms yields the mean-zero local-error component, while the remaining terms form the remainder.The mean-zero component’s first, third, and fourth terms are odd in the Brownian path; its second is a centered quadratic functional.Source: §B.3 · pp. 40–41
  • For SRK-II, the stage terms differ from SRK-I, but collecting them gives the same mean-zero-plus-remainder decomposition.The two integrators’ mean-zero terms use different scalar coefficients: c1=3/2 and c3=5/2 for SRK-II, versus c1=3/4 and c3=9/4 for SRK-I.Source: §B.3 · pp. 40–41
  • The resulting moment bounds give the mean-zero term order h²(1+|Z0|) and control the remainder for both integrators.The proof uses the same remainder-moment argument for SRK-II and SRK-I.Source: §B.3 · pp. 41–42

B.4. Proof of Lemma 2: product decomposition

The product-error decomposition centers a leading product term by subtracting its conditional mean. The surrounding proofs use these local decompositions to establish strong- and weak-error bounds.

B.4. Proof of Lemma 2: product decomposition

  • The leading product term is not conditionally centered, so the decomposition subtracts its conditional mean to define M(∆Q0) and preserve the identity.The other product terms have order at most h³.Source: §B.4 · p. 43
  • The product decomposition’s moment estimates bound both its mean-zero component and its remainder for h≤1.The deterministic terms in the remainder are bounded by Ch³(1+|Z0|).Source: §B.4 · p. 43

B.6. Proof of Corollary 2: weak error analysis

  • The weak-error proof bounds each of the N telescoping summands and obtains the stated finite-time estimate when Nh≤T.Each summand is bounded by Ce^(CT)Lh³(1+|Z0|)⁶.Source: §B.6 · pp. 47–48

C. Proofs for Section 4

  • The proof overview identifies Theorem 2 as the mean-square-error bound and Theorem 3 as the uniform Wasserstein bound, combining stability with consistency estimates.The mean-square analysis bounds the three terms in the time-averaged error formula; the Wasserstein analysis sums the local weak error.Source: C. Proofs for Section 4 · p. 48

C.1. Proof of Lemma 3

A Lyapunov argument establishes uniform moment and displacement bounds for the integrator chain. It relies on dissipativity outside a bounded region and a sufficiently small step size.

C.1. Proof of Lemma 3

  • The Lyapunov proof uses that each integrator’s stage drift averages gradients at nearby stage points with nonnegative weights summing to one.The argument uses the stages’ distance bound from the step-start state.Source: §C.1 · p. 48
  • Dissipativity outside a bounded region supplies the drift condition used in the Lyapunov estimate.The proof applies the dissipativity bound to pairs involving the origin and controls the bounded region separately.Source: §C.1 · p. 48
  • For a sufficiently small step size, iteration of the Lyapunov estimate gives moments of Zk uniformly in k, along with displacement bounds.The step-size threshold depends on m, M, R, d, and the moment exponent q.Source: §C.1 · p. 50

C.2. Proof of Theorem 2

The proof uses the time-averaged error formula and the discrete Poisson equation to turn one-step error estimates into a bound, treating mean-zero terms and remainders separately.

C.2. Proof of Theorem 2

  • The time-averaged error is split into a boundary term, discretization terms, and martingale differences, which are estimated in turn.The proof uses a Taylor expansion of the propagated test function to express the discretization terms and their remainder.Source: §C.2 · pp. 50–51
  • The Taylor-remainder estimates use bounds on derivatives of the discrete Poisson solution through third order and moment bounds for the one-step displacements.The proof applies moment estimates to remove conditioning and then combines the resulting bounds in the time-averaged error estimate.Source: §C.2 · pp. 51–52
  • The mean-zero terms have zero conditional mean and are orthogonal in L2, while the remainders are bounded in L2 and accumulate linearly in the number of steps.The proof also bounds the martingale-difference contribution using conditional second moments.Source: §C.2 · pp. 51–53
  • The final estimate treats Nh ≥ 1 and Nh < 1 separately, absorbing smaller terms into the stated bound in either case.For Nh < 1, (Nh)^-1 exceeds 1 and the proof uses a bound that holds for every N.Source: §C.2 · p. 53

C.3. Proof of Theorem 3

The proof propagates test functions with the exact semigroup and sums local weak errors; derivative decay controls that sum uniformly in the number of steps.

C.3. Proof of Theorem 3

  • Decay of the derivatives of the exactly propagated test function makes the sum of local weak errors converge uniformly in N.The proof uses the local weak-error estimate together with derivative bounds for the propagated function.Source: §C.3 · pp. 53–54
  • To obtain the Wasserstein-1 bound, the proof tests against smooth functions with gradient norm at most one and extends to 1-Lipschitz functions by mollification.Finite first moments of both compared laws justify the Wasserstein-1 comparison.Source: §C.3 · p. 54
  • An averaging argument bounds the propagated test-function discrepancy using invariance of π, its finite first moment, and exponential decay in Nh.The triangle inequality for W1 then yields the claimed comparison of the laws.Source: §C.3 · p. 54

D. Integrators of Table 1 and SRK-LD

This section describes the integrators and the conditions behind the convergence comparisons. It also explains how randomized splitting and the reference scheme are used.

D. Integrators of Table 1 and SRK-LD

  • The schemes' convergence entries rely on different hypotheses, and the table reports established rather than nominal orders.The conditions for this work include convexity outside a bounded region; the table's weak entries are the least orders guaranteed by the cited finite-horizon analyses, with stated exceptions.Source: D. Integrators of Table 1 and SRK-LD · p. 55
  • Leimkuhler–Matthews reuses each Gaussian across successive steps, so its chain is not Markov, although its accumulated noise variance agrees to leading order.Its per-step noise variance is h per coordinate rather than 2h, while over n steps the variances are 2nh − h and 2nh, respectively.Source: D. Integrators of Table 1 and SRK-LD · pp. 55–56
  • Either fixed ordering of the drift and diffusion substeps is a first-order Lie–Trotter splitting, while randomizing between the two orderings cancels the leading commutator and raises the order.Source: D. Integrators of Table 1 and SRK-LD · p. 56
  • Random splitting LMC has strong order 1; solving its drift flow with a solver of order at least two leaves the table's order unchanged.The experiments solve the drift flow using the third-order Runge–Kutta method RK3.Source: D. Integrators of Table 1 and SRK-LD · p. 56
  • SRK-I and SRK-II use both the Brownian increment and its time integral over each step, with that pair simulated exactly.The schemes evaluate the gradient and no higher derivatives.Source: D. Integrators of Table 1 and SRK-LD · pp. 55, 57
  • SRK-LD is a three-stage scheme with three gradient evaluations per step, used to compute references rather than included in Table 1.Its cited analysis assumes global strong convexity, and its strong order 3/2 suffices for references with fine steps.Source: D. Integrators of Table 1 and SRK-LD · p. 57
Loading 2609.11528v4…