Source-linked AI summary
Convergence of Numerical Time-Averaging and Stationary Measures via Poisson Equations
Jonathan C. Mattingly, Andrew M. Stuart, M. V. Tretyakov
TL;DR
The paper addresses how to approximate the long-time behavior and stationary measure of an SDE using numerical trajectories. It uses an associated Poisson equation to obtain time-averaging error estimates and transfer finite-time weak accuracy to stationary-measure convergence. The approach applies broadly to explicit and implicit schemes, including hypoelliptic cases, within a torus and smooth-test-function setting.
Problem
The paper asks how numerical time averages and stationary measures can be shown to approximate those of an SDE over long times.
Method
The analysis uses an associated Poisson equation for the SDE to control numerical time averages and stationary measures.
Results
The numerical method’s stationary measures are close to the SDE’s stationary measure, with convergence matching the finite-time weak convergence rate.
Takeaways & Limitations
The Poisson-equation approach provides a simple, universal analysis applicable to explicit and implicit schemes, simple random variables, and hypoelliptic SDEs.
Takeaways & Limitations
The main results are presented for a compact torus and smooth test functions; extending them to R^d requires controlled solutions of the relevant Poisson equation.
Abstract
from arXiv · showhide
Numerical approximation of the long time behavior of a stochastic differential equation (SDE) is considered. Error estimates for time-averaging estimators are obtained and then used to show that the stationary behavior of the numerical method converges to that of the SDE. The error analysis is based on using an associated Poisson equation for the underlying SDE. The main advantage of this approach is its simplicity and universality. It works equally well for a range of explicit and implicit schemes including those with simple simulation of random variables, and for hypoelliptic SDEs. To simplify the exposition, we consider only the case where the state space of the SDE is a torus and we study only smooth test functions. However we anticipate that the approach can be applied more widely. An analogy between our approach and Stein's method is indicated. Some practical implications of the results are discussed.
1. Introduction.
The paper develops a simple Poisson-equation analysis for transferring finite-time weak-approximation accuracy to long-time averages and stationary measures. Its scope includes explicit and implicit schemes and hypoelliptic SDEs, while the exposition is restricted to a torus and smooth test functions.
- Contribution: The paper estimates time-averaging errors and uses them to show numerical stationary measures approach the SDE’s stationary measure.The convergence rate matches the numerical integrator’s weak convergence rate on finite time intervals.
- Contribution: The analysis uses an associated Poisson equation and does not require the numerical Markov chain to be uniquely ergodic.Every stationary measure of the numerical method is shown to be close to the unique stationary measure of the SDE.
- Generality: The approach applies to explicit and implicit schemes, simple random variables, and hypoelliptic SDEs.The authors emphasize simplicity, universality, and reliance on classical PDE results.
- Scope: The exposition restricts the state space to a compact torus and considers relatively smooth test functions.Extending the results to R^d requires control of time spent outside the phase-space center.
- Organization: The paper organizes its analysis around the SDE setting, numerical methods, the auxiliary Poisson equation, main error results, and statistical implications.The later discussion classifies numerical time-averaging errors and addresses statistical error.
2. SDE Setting.
The SDE is modeled as a Markov process on a torus with Lipschitz drift and diffusion, whose randomness must spread sufficiently to yield smoothing. Under elliptic or hypoelliptic assumptions, the process has a unique stationary measure with a density.
- SDE dynamics: The SDE has drift f and independent random kicks in the directions given by the columns of g.The generator and associated Markov semigroup describe the resulting dynamics.
- Hypoellipticity: The noise must spread across directions sufficiently to produce smoothing of probability densities, with hypoellipticity covering cases beyond uniform ellipticity.Lie brackets of vector fields characterize how additional directions are generated.
- Applications: The torus setting also represents periodic-boundary applications, including noisy gradient systems used to sample Gibbs distributions.The authors also mention possible applications on compact smooth manifolds.
- Assumptions: The standing assumptions include either positive-definite diffusion or smooth hypoelliptic coefficients whose generated directions span R^d.In the hypoelliptic case, the spanning condition is combined with existence of a unique stationary measure.
- Stationarity: Under either elliptic or hypoelliptic assumptions, the SDE has a unique stationary measure μ with a density on the torus.This stationary measure is the reference distribution for the numerical analysis.
- Purpose: The paper studies numerical approximations whose time averages remain close to the SDE’s corresponding ergodic limits.The numerical methods are introduced to approximate the dynamics on the torus.
3. Numerical approximations.
The paper considers generalized Euler–Maruyama-type methods on the torus, characterized by drift, diffusion, and random variables satisfying local consistency conditions. These conditions yield first-order weak convergence on finite time intervals and motivate comparison with the SDE dynamics.
- General framework: The numerical methods are defined using a time increment Δ, coefficient functions F and G, and i.i.d. random variables η_n.The random variables satisfy moment and adaptation requirements.
- Moment conditions: The moment order r is left unspecified because common schemes use random variables with bounded moments of every order.The required finite moment order can still depend on the proof.
- Consistency: The local-error assumption is imposed so that the numerical dynamics approximate the SDE over one time step.Under this assumption, X_n approximates X(t_n) with t_n=nΔ.
- Accuracy: First-order weak convergence on finite time intervals follows from the stated local consistency conditions.This provides the finite-time accuracy used in the paper’s long-time analysis.
- Generator comparison: The associated operator captures the leading-order part of the numerical method’s generator, although it is not itself the Markov process generator.Its closeness to the SDE generator motivates expecting close distributions.
- Motivation: The numerical method’s leading generator part is close to that of the SDE, supporting closeness of their distributions.This is presented as a reasonable expectation rather than the main convergence theorem.
- Examples: The framework includes explicit Euler–Maruyama and an implicit split-step method.Both methods are presented as examples satisfying the framework’s assumptions.
- Higher-order methods: The local consistency characterization also provides a starting point for analyzing higher-order numerical methods.The paper later considers more general methods in its higher-order section.
4. Poisson equation.
The paper uses a Poisson equation for the SDE to establish laws of large numbers and quantitative convergence estimates for time averages. Under the stated assumptions, the equation has a unique regular solution, enabling control of stochastic and numerical long-time averages.
- Poisson-equation approach: Solving the relevant Poisson equation provides a general route to proving limits of time averages.The solution controls the fluctuation of the semigroup from the stationary average.
- Poisson-equation approach: Under the stated assumptions, the Poisson equation has a unique solution with enhanced regularity for smooth right-hand sides.Elliptic assumptions yield two additional derivatives, while hypoelliptic assumptions yield a positive regularity gain.
- Scope of the method: Extending the analysis from the torus to general ergodic SDEs on R^d requires well-controlled solutions of the corresponding Poisson equation.The paper identifies this regularity and existence theory as a principal technical issue.
- Continuous-time illustration: The proof uses the Poisson equation with Itô’s formula to decompose time averages into vanishing boundary terms and martingale terms.Boundedness of the solution and its derivatives controls the resulting terms.
- Continuous-time illustration: The resulting strong law holds from every initial condition, rather than only almost every initial condition under the stationary measure.The argument also supplies quantitative convergence-rate estimates.
5. Error analysis for the numerical time-average.
The section derives error estimates for discrete time-averaging estimators using sample-path analysis, including mean, L2, and almost-sure results. These estimates quantify numerical and statistical errors and support convergence analysis for first- and higher-order schemes.
- Estimator and setup: The analysis studies discrete time averages over a simulation horizon T = N∆ and compares them with stationary averages.The estimator is introduced for the stationary average, with constants independent of ∆ and T and linear dependence on the test-function norm.
- Mean convergence: Theorem 5.1 establishes a mean convergence result for the numerical time-average under elliptic or hypoelliptic regularity assumptions.The required test-function space is W 2,∞ in the elliptic setting and W 4,∞ in the hypoelliptic setting.
- Proof structure: The proofs decompose the time-average error into local numerical-error terms and martingale contributions controlled through moment estimates, inequalities, and the Borel-Cantelli lemma.The argument also uses bounded derivatives of the Poisson-equation solution and boundedness of the scheme coefficients.
- L2 and almost-sure convergence: Theorem 5.2 extends the analysis to L2 convergence, while Theorem 5.3 gives an almost-sure result for sufficiently small ∆ and sufficiently large T.The almost-sure bound uses an a.s. bounded random variable depending on ε and the particular test function.
- Stationarity assumptions: The results do not require the numerical Markov chain to be ergodic, so any stationary measure can be analyzed through the time-average estimates.If the long-time estimator limit exists independently of the initial condition, the almost-sure result yields the corresponding stationary-average conclusion.
- Error interpretation: The error analysis separates numerical and statistical effects, with the statistical contribution scaling as T = 1/N^1/2∆^1/2 in the stated discussion.The section notes that large-scale simulations may be dominated by statistical error when numerical error is relatively small.
6. Error analysis for the numerical stationary measures.
The paper proves that stationary measures of the numerical method are close to the SDE’s stationary measure, without requiring unique ergodicity of the numerical method. The analysis uses a metric based on smooth test functions and connects the argument to Stein’s method.
- Stationary-measure error: The metric ρ is defined through test functions in a smoothness class H, with different Sobolev regularity for elliptic and hypoelliptic settings.The paper uses H = W 2p,∞ in the elliptic setting and H = W 2(p+1),∞ in the hypoelliptic setting.
- Stationary-measure error: Any stationary measure of the numerical method is close to the SDE’s unique stationary measure in the metric ρ.The result applies without assuming that the numerical method is uniquely ergodic.
- Proof strategy: The stationary-measure result follows by combining finite-time approximation estimates with a uniform time-average bound and then letting the averaging horizon tend to infinity.The proof invokes the relevant finite-time error theorem and uses uniformity over test functions in H.
- Proof strategy: The Poisson-equation approach does not require unique ergodicity of the numerical Markov chain, so every stationary measure satisfies the same closeness conclusion.This is presented as a central distinction from arguments that first establish uniqueness for the numerical method.
- Relationship to Stein’s method: The analysis parallels Stein’s method: a small generator discrepancy and a solvable Poisson equation yield a bound on the distance between stationary measures.Taking a supremum over a determining class of test functions produces a metric estimate.
7. Variance of the Empirical Time Average.
The variance of the numerical time-average estimator contains discretization and finite-time contributions, while mixing assumptions control its statistical component. The paper also separates practical error sources for time-averaging and ensemble-averaging estimators.
- Error decomposition: The practical time-averaging error has three components: numerical integration error, finite-time bias from distance to stationarity, and statistical error.The paper associates these with bounds involving K∆ or K∆p, K/T, and sampling fluctuations.
- Variance estimate: Var(ˆφN) = O(∆2 + 1/T), combining a time-step contribution with a contribution that decreases as the integration time T grows.A mixing-type condition is needed to obtain this variance estimate.
- Mixing assumptions: A relaxed mixing condition for the numerical chain is sufficient for the variance proposition, and faster decorrelation does not improve the resulting estimate.The condition is weaker than exponential mixing of the continuous-time process.
- Practical estimation: The paper recommends block-based variance estimation for long trajectories, treating sufficiently long blocks as approximately uncorrelated.Confidence intervals can use c = 2 for probability 0.95 or c = 3 for probability 0.997.
- Ensemble averaging: For ensemble averaging, the total error separates into approximation, numerical integration, and Monte Carlo terms controlled respectively by t, ∆, and the number of independent trajectories L.The numerical integration term is C∆p, where p is the weak order of the method.