Higher-Order SDE Discretization Methods

• Higher-order SDE discretization methods are numerical schemes using stochastic Taylor expansions to achieve convergence orders beyond the baseline 0.5, with both strong and weak accuracy improvements.
• They encompass approaches like Milstein-type, stochastic Runge–Kutta, and derivative-free exponential methods (e.g., ERKM1.5) that balance precision with computational cost in high-dimensional and nonlinear systems.
• Applications span SPDEs, mathematical finance, and machine learning, where adaptive multi-scale time-stepping and energy-stable formulations enhance practical performance and robustness.

Higher-order SDE-discretization methods comprise a class of numerical schemes designed to approximate solutions of stochastic differential equations (SDEs) and stochastic partial differential equations (SPDEs) with strong or weak convergence rates that exceed the canonical order $1/2$ achieved by the Euler–Maruyama method. Motivated by applications in stochastic analysis, mathematical finance, SPDEs, machine learning, and statistical physics, such methods aim to balance increased stochastic Taylor expansion accuracy with computational tractability, especially in systems with potentially high-dimensional state spaces and complex nonlinearities.

1. Methodological Principles of Higher-Order SDE Discretization

The construction of higher-order SDE discretization methods is fundamentally rooted in stochastic Taylor expansions, which represent the exact solution to an SDE as a series in powers of the timestep $h$, involving iterated Itô or Stratonovich integrals. To achieve a desired global convergence order $p$, the discretization must replicate the expansion of the exact solution up to $O(h^{p+1})$ in the strong sense (mean-square convergence) or $O(h^{p})$ in the weak sense.

For strong order 1.5 ($p=3/2$), as exemplified by the exponential stochastic Runge–Kutta (ERKM1.5) method for SPDEs of Nemytskii (multiplicative noise) type, the algorithm is derived by matching the series expansion up to $h^{2.5}$ and imposing algebraic conditions on the Butcher tableau entries for both drift and diffusion contributions. This matches not only the increments but selected covariances and mixed stochastic terms, ensuring that local mean-square errors are $O(h^{2.5})$ and that global strong convergence $O(h^{1.5-\epsilon})$ is achieved for arbitrarily small $\epsilon > 0$ (Hallern et al., 2024).

2. Algorithmic Frameworks and Main Schemes

Higher-order SDE discretization schemes include:

• Milstein-type methods: Add a corrective term involving the Lie derivative of the diffusion, attaining strong order one under commutativity. For systems with non-commutative noise, achieving higher order requires including further iterated integrals.
• Stochastic Runge–Kutta (SRK) and exponential SRK: These perform multiple (explicit or implicit) drift and diffusion stage evaluations per step, with algebraically determined stages and random weights.
• Derivative-free exponential SRK (ERKM1.5): For SPDEs with multiplicative noise, six-stage updates combine exponential propagators, drift/diffusion evaluations, and finite-difference approximations of derivatives, utilizing Gaussian random weights simulating required iterated-integral statistics. The ERKM1.5 method is characterized by (Hallern et al., 2024):
  • Six-stage derivative-free exponential update.
  • Drift/diffusion increments weighted by simulated standard normals via explicit covariance forms.
  • Full strong mean-square convergence up to order $1.5-\epsilon$ (for solutions with regularity index $\gamma$ arbitrarily close to $1.5$).
• High-order splitting and path-based methods: Replace the stochastic driving signal with tailored piecewise linear interpolants that precisely reproduce key iterated integrals in law or expectation for commutative noise. Foster–dos Reis–Strange establish that if-iterated integrals up to order $p$ are matched appropriately, strong order $p-\tfrac12$ is obtained. For commutative SDEs, this gives strong order 1.5 with only ODE integration on small pieces (Foster et al., 2022).
• Variational and energy-stable methods for hyperbolic SPDEs: In the stochastic wave equation setting with multiplicative noise, higher-order methods based on mid-point and BDF-style time integration, together with a discrete Itô correction, achieve strong convergence order $3/2$ for the principal component (Feng et al., 2022).
• Adaptive multi-scale time-stepping: Adaptive methods deploy locally refined step sizes based on the proximity to domain boundaries or other error indicators. For exit time problems in Itô diffusions, combining Milstein or strong Itô–Taylor ($\gamma=1.5$) integrators with a multi-layer adaptive step hierarchy reduces strong exit-time error to $O(h^{\gamma-\xi})$ with log cost overhead, exemplifying the synergy of higher-order strong schemes with adaptivity (Hoel et al., 2022).

3. Order Conditions, Regularity Requirements, and Limitations

The realization of strong convergence orders beyond one demands stringent smoothness conditions:

• Coefficient regularity: To attain order $1.5$, coefficients must generally be twice continuously Fréchet differentiable, with globally bounded first and second derivatives. For example, the ERKM1.5 construction requires $F\in C_b^2$ and $B\in C_b^2$ (Hallern et al., 2024). In some settings (e.g., weakly regular SPDEs or SDEs with non-globally monotone coefficients), specialized tamed, stopped, or finite-difference approximations are introduced to control moments and maintain integrability (Dai et al., 2024).
• Iterated stochastic integrals: For noncommutative or multidimensional noise, access to simulated or approximated multiple Wiener–Itô integrals is required, increasing algorithmic complexity.
• Solution regularity: To propagate high regularity from initial data to solutions, the problem must ensure sufficient spatial and pathwise smoothness, often forcing $H^\gamma$-regularity up to order $1.5$ for SPDEs, or high moment bounds for SDEs.

Barriers: For certain SDEs—notably models with non-Lipschitz or discontinuous drift (e.g., financial SDEs such as the log-Heston model)—there exist sharp theoretical barriers. Even the best possible equidistant-grid method is limited by a maximum strong order $\min\{\nu,1/2\}$, with $\nu$ the Feller index of the underlying diffusion (Mickel et al., 2022). Thus, higher-order discretizations offer no additional asymptotic gain unless the underlying regularity supports it.

4. Computational Cost and Practical Implementation

Higher-order SDE schemes necessarily increase per-step computational effort, but recent designs focus on computationally efficient, derivative-free, or adaptive variants:

Method	Strong Order	Function Calls per Step	Derivatives Required	Notable Features
Euler–Maruyama	0.5	$2d$	None	Baseline, lowest cost
Derivative-free Milstein	1.0	$\sim 4d$	None	Efficient, commutative noise
Exponential Wagner–Platen	1.5	$\sim 2(d+d^2+d^3)$	1st and 2nd	High cost, high order
ERKM1.5	up to 1.5	$11d$	None	Efficient, high-order, derivative-free

For multi-component (large $d$) systems, ERKM1.5 dramatically improves over classical derivative-based schemes by eliminating the cubic complexity in $d$ while retaining order $1.5$ strong convergence (Hallern et al., 2024).

Adaptive methods retain the optimal order even near domain boundaries for exit problems, at $O(h^{-1}\log(h^{-1}))$ expected cost (Hoel et al., 2022). In splitting methods, each high-order step typically requires two or three Gaussian random draws and $O(m)$ ODE-flow evaluations (Foster et al., 2022).

5. Applications and Comparative Performance

SPDEs with Multiplicative Noise

The ERKM1.5 scheme for infinite-dimensional SPDEs with Nemytskii-type multiplicative noise combines exponential semigroup propagators with derivative-free high-order stochastic integration. Numerical benchmarks on the stochastic heat equation confirm that ERKM1.5 matches the 1.5-order convergence slope (strictly, $1.5-\epsilon$) of the best Wagner–Platen schemes but at vastly reduced computational cost and without high-order derivatives (Hallern et al., 2024).

Higher-Order Methods for Stochastic Wave Equations

For stochastic wave equations with general drift and multiplicative noise, tailored implicit BDF-type time stepping with a discrete Itô correction achieves the strong order $3/2$ under noise structures depending only on the position, with variational stability analysis ensuring robust error control (Feng et al., 2022).

High-Dimensional and Stiff SDEs

Stability-optimized explicit and semi-implicit Runge–Kutta methods (SOSRA/SOSRI, SKenCarp) with embedded stiffness detection allow for large time steps and switch to L-stable integrators in pathwise stiff regimes. These methods guarantee strong order $1.5$ and deliver substantial CPU speed-ups for large and stiff biological models (Rackauckas et al., 2018).

Diffusion Models in Machine Learning

In generative diffusion models for high-dimensional data, higher-order SDE discretizations (e.g., 1.5-order Chang–Platen) are proven to maintain their theoretical $O(h^{1.5}+\varepsilon)$ advantage in 2-Wasserstein distance over Euler–Maruyama, provided the score function estimators are sufficiently smooth. Numerical results confirm the asymptotic advantage in latent-space and certain data regimes, motivating further adoption in computational learning (Pfarr et al., 26 Jan 2026).

6. Extensions, Limitations, and Outlook

Higher-order schemes are limited by:

• Coefficient regularity and commutativity: Noncommutative noise and irregular coefficients (e.g., in finance SDEs or jump-diffusions with discontinuous drift) impose theoretical order barriers, with $3/4$ provably sharp in the jump-adapted quasi-Milstein context (Przybyłowicz et al., 2022, Mickel et al., 2022).
• High-dimensionality: Complexity of required stochastic iterated integrals may offset gains for very large systems in the absence of efficient simulation techniques.
• SPDE discretization: Combined high-order space–time methods are an active area of research, with techniques such as high-order path-based splittings, high-regularity finite elements, and stabilized interpolation for nonlinear drifts advancing the state of the art (Li et al., 2023, Hallern et al., 2024).

Emerging directions include adaptive multi-stage, high-order Taylor and splitting schemes for broad classes of SPDEs, hybrid derivative-free formulations for non-globally monotone systems (Dai et al., 2024), and the systematic integration of high-order strong schemes into multilevel Monte Carlo and stochastic optimal control.

---

References:

• "An Exponential Stochastic Runge-Kutta Type Method of Order up to 1.5 for SPDEs of Nemytskii-type" (Hallern et al., 2024)
• "High order splitting methods for SDEs satisfying a commutativity condition" (Foster et al., 2022)
• "Higher order time discretization for the stochastic semilinear wave equation with multiplicative noise" (Feng et al., 2022)
• "Stability-Optimized High Order Methods and Stiffness Detection for Pathwise Stiff Stochastic Differential Equations" (Rackauckas et al., 2018)
• "Analyzing the Error of Generative Diffusion Models: From Euler-Maruyama to Higher-Order Schemes" (Pfarr et al., 26 Jan 2026)
• "Higher order numerical methods for SDEs without globally monotone coefficients" (Dai et al., 2024)
• "A higher order approximation method for jump-diffusion SDEs with discontinuous drift coefficient" (Przybyłowicz et al., 2022)
• "Second order discretization of Backward SDEs" (Crisan et al., 2010)
• "Higher-order adaptive methods for exit times of Itô diffusions" (Hoel et al., 2022)
• "The order barrier for the $L^1$-approximation of the log-Heston SDE at a single point" (Mickel et al., 2022)
• "Higher order time discretization method for a class of semilinear stochastic partial differential equations with multiplicative noise" (Li et al., 2023)