High-Order Langevin Monte Carlo Methods

• High-Order Langevin Monte Carlo is a family of sampling methods that use higher-order discretizations, refined continuous dynamics, or geometric proposals to target distributions more accurately.
• These methods overcome the limitations of Euler–Maruyama by reducing discretization bias, improving convergence rates, and enabling efficient sampling in high-dimensional settings.
• They incorporate techniques such as tamed Itô–Taylor schemes, stochastic Runge–Kutta methods, and lifted dynamics to achieve provable improvements in mixing times and invariant measure accuracy.

High-Order Langevin Monte Carlo denotes a family of sampling methods that improve on first-order Langevin Monte Carlo by using higher-order numerical discretizations, higher-order continuous-time dynamics, or higher-order proposal constructions while targeting a distribution of the form $\pi(x)\propto e^{-U(x)}$. In the recent literature, the label covers strong-order $1.5$ overdamped schemes, third- and $K$-th order lifted Langevin diffusions, randomized high-order integrators for kinetic Langevin and unadjusted Hamiltonian Monte Carlo, microcanonical projected dynamics, and Metropolis-adjusted proposals designed through higher-order cancellation in the acceptance-ratio expansion (Li et al., 2019, Sabanis et al., 2018, Mou et al., 2019, Mahajan et al., 21 Oct 2025, Dang et al., 24 Aug 2025, Durmus et al., 2015).

1. Scope and meanings of “high-order”

The canonical starting point is the overdamped Langevin diffusion

$dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$

whose invariant law is $\pi(x)\propto e^{-U(x)}$. For its discretizations, the literature distinguishes strong and weak convergence. A scheme has strong order $p$ if $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$ on a fixed horizon, and weak order $p$ if $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$ for smooth test functions $\varphi$ (Yang et al., 8 May 2026).

Across the literature, “high-order” is used in several non-equivalent senses. One usage refers to higher-order discretization of the overdamped SDE, as in Itô–Taylor, stochastic Runge–Kutta, or splitting methods (Sabanis et al., 2018, Li et al., 2019, Yang et al., 8 May 2026). A second usage refers to higher-order continuous dynamics: third-order and more general $1.5$0-th order Langevin systems introduce auxiliary variables so that the physical position becomes smoother in time and is more amenable to accurate integration (Mou et al., 2019, Mahajan et al., 21 Oct 2025, Dang et al., 24 Aug 2025). A third usage is geometric: microcanonical Langevin Monte Carlo employs a second-order minimum-norm integrator and projected noise on a constant-energy manifold, so “high-order” refers to the order of the integrator rather than to higher derivatives (Bayer et al., 2023). A fourth usage appears in Metropolis-adjusted proposals such as fMALA, where the improvement comes from cancelling low-order terms in the expansion of the Metropolis–Hastings log-ratio; in that sense the proposal is “high-order” even though the underlying discretization still has weak order $1.5$1 (Durmus et al., 2015).

This multiplicity of meanings is not merely terminological. It reflects distinct algorithmic strategies for overcoming the limitations of Euler–Maruyama: discretization bias, dimension-dependent step-size restrictions, instability under superlinear drift, and poor asymptotic scaling in Wasserstein distance or effective sample size.

2. Overdamped high-order discretizations

A first major branch of High-Order Langevin Monte Carlo keeps the overdamped diffusion as the base process and replaces Euler–Maruyama by a more accurate scheme. The 2018 Higher Order Langevin Monte Carlo algorithm constructs a tamed strong-order $1.5$2 Itô–Taylor discretization of the overdamped Langevin SDE. In the superlinear case it uses tamed versions of $1.5$3, $1.5$4, and $1.5$5, and proves convergence of order $1.5$6 in Wasserstein-2 and order $1.5$7 in weighted total variation; in the globally Lipschitz strongly convex setting the Wasserstein-2 bias becomes order $1.5$8 (Sabanis et al., 2018). The same paper emphasizes that the $1.5$9 term is essential: without it, the rate drops to order $K$0 in Wasserstein-2.

A second branch uses stochastic Runge–Kutta constructions. “Stochastic Runge-Kutta Accelerates Langevin Monte Carlo and Beyond” studies an overdamped Langevin sampler, SRK-LD, obtained from a mean-square order $K$1 stochastic Runge–Kutta integrator. For strongly convex potentials that are smooth up to a certain order, it yields $K$2 iterations in $K$3-Wasserstein distance, improving the $K$4-dependence over Euler-based overdamped methods while using only a gradient oracle in the algorithm (Li et al., 2019). “Accelerating Langevin Monte Carlo via Efficient Stochastic Runge--Kutta Methods beyond Log-Concavity” develops a Hessian-free strong-order $K$5 method, RKLMC-2G, that uses two gradient evaluations per step, achieves the uniform-in-time bound

$K$6

and therefore mixing time $K$7 under dissipativity, gradient/Hessian/third-derivative Lipschitzness, and a Log-Sobolev inequality (Yang et al., 8 May 2026).

A third line of work emphasizes invariant-measure accuracy rather than finite-horizon strong order. “Ergodicity and error estimate of laws for a random splitting Langevin Monte Carlo” shows that randomizing the order of drift and diffusion substeps in overdamped Langevin cancels the first-order bias in the invariant measure. The resulting random splitting Langevin Monte Carlo has $K$8 invariant-measure bias in Wasserstein-1, obtained through sharp local relative-entropy estimates and geometric ergodicity via reflection coupling (Li et al., 9 Oct 2025).

More recent work extends these ideas to non-convex and superlinear settings. “Non-asymptotic estimates for accelerated high order Langevin Monte Carlo algorithms” proposes aHOLA and aHOLLA, both based on the strong-order $K$9 Kloeden–Platen scheme. Under a local Hölder condition of exponent $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$0, convexity at infinity, and superlinear growth controlled by polynomial taming, aHOLA attains $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$1-error of order $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$2 and $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$3-error of order $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$4; aHOLLA obtains the same exponents under global continuity and dissipativity conditions in the linear-growth regime (Neufeld et al., 2024).

These overdamped methods separate into Hessian-free and derivative-rich regimes. HOLA and aHOLA explicitly use Hessian and third-derivative information, whereas SRK-LD and RKLMC-2G are analyzed using higher derivatives but evaluate only gradients in the algorithm (Sabanis et al., 2018, Yang et al., 8 May 2026).

3. Lifted third-order and $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$5-th order Langevin dynamics

A second major branch changes the continuous-time dynamics rather than only the discretization. “High-Order Langevin Diffusion Yields an Accelerated MCMC Algorithm” introduces a third-order Langevin diffusion on $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$6: $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$7 with invariant density proportional to $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$8 (Mou et al., 2019). The construction makes the position variable smoother in time because the Brownian forcing acts two integrations away from $dX_t = -\nabla U(X_t)\,dt + \sqrt{2}\,dW_t,$9. For ridge-separable generalized linear model potentials it yields mixing time $\pi(x)\propto e^{-U(x)}$0; for general strongly convex potentials with $\pi(x)\propto e^{-U(x)}$1-th order smoothness it yields

$\pi(x)\propto e^{-U(x)}$2

up to condition-number factors and logarithms (Mou et al., 2019).

Two later works generalize this lifting principle. “The Picard-Lagrange Framework for Higher-Order Langevin Monte Carlo” defines $\pi(x)\propto e^{-U(x)}$3-th order Langevin dynamics on $\pi(x)\propto e^{-U(x)}$4, with Hamiltonian $\pi(x)\propto e^{-U(x)}$5, noise injected only in the last block, and a $\pi(x)\propto e^{-U(x)}$6-diffusion structure ensuring that the $\pi(x)\propto e^{-U(x)}$7-marginal is the target distribution (Mahajan et al., 21 Oct 2025). Its discretization uses Lagrange interpolation in time together with Picard iterations to solve the resulting fixed-point equations, while requiring only gradient evaluations. Under strong log-concavity and bounded higher-order derivatives, the resulting query complexity is

$\pi(x)\propto e^{-U(x)}$8

with improving $\pi(x)\propto e^{-U(x)}$9-dependence as $p$0 increases (Mahajan et al., 21 Oct 2025).

“High-Order Langevin Monte Carlo Algorithms” develops a related but distinct family of $p$1-th order Langevin dynamics with $p$2 auxiliary variables, discretized by splitting and accurate integration. For $p$3, its mixing time satisfies

$p$4

under strong convexity, smoothness, and a derivative-growth condition on high-order derivatives (Dang et al., 24 Aug 2025). For $p$5, this becomes $p$6, and the paper presents explicit Gaussian one-step updates for fourth-order Bayesian linear and logistic regression models (Dang et al., 24 Aug 2025).

The common theme is that higher-order lifted dynamics trade extra state variables for smoother temporal behavior of the position variable. This suggests, and the cited results substantiate, that part of the acceleration comes from making polynomial or splitting approximations of the drift substantially more accurate than they are for the overdamped diffusion.

4. Geometric, kinetic, and microcanonical variants

High-order ideas also enter kinetic and geometric samplers. “Randomized Runge-Kutta-Nyström Methods for Unadjusted Hamiltonian and Kinetic Langevin Monte Carlo” constructs $p$7- and $p$8-order $p$9-accurate randomized Runge–Kutta–Nyström integrators for Hamiltonian flows and inserts them into unadjusted HMC and unadjusted kinetic Langevin algorithms (Bou-Rabee et al., 2023). The $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$0-order method uses two force evaluations per step, the $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$1-order method uses three, and both are analyzed under gradient and Hessian Lipschitz assumptions. In numerical experiments, the resulting unadjusted samplers reduce bias per computational cost relative to Verlet and earlier randomized integrators, with especially clear gains for kinetic Langevin and well-behaved high-dimensional targets (Bou-Rabee et al., 2023).

A different geometric route is microcanonical sampling. “Field-Level Inference with Microcanonical Langevin Monte Carlo” studies the stochastic dynamics

$$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$2

where $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$3 is constrained to approximately unit norm and both drift and noise are projected orthogonally to $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$4 (Bayer et al., 2023). The method uses a second-order minimum-norm integrator, has no Metropolis–Hastings correction, and relies on controlling energy fluctuations per dimension. In cosmological field-level inference at dimension $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$5, it achieved effective sample size per gradient evaluation improvements of $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$6 for field modes and $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$7 for cosmological parameters relative to HMC at $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$8 resolution (Bayer et al., 2023). Here the “high-order” aspect is explicitly geometric: second-order integration and benign energy-error growth, not higher derivatives.

These geometric schemes broaden the scope of High-Order Langevin Monte Carlo beyond overdamped discretization accuracy. They show that high-order integration can be exploited in non-reversible and constrained phase-space dynamics, sometimes with empirical gains that are more strongly dimension-dependent than the asymptotic Wasserstein theory available for overdamped schemes.

5. Rates, scaling laws, and representative results

The literature exhibits several distinct complexity regimes. Some methods improve the dependence on the error tolerance $$\sup_{0\le n\le N}(\mathbb{E}|X_{t_n}-Y_n|^2)^{1/2}=O(h^p)$$9, some improve the dependence on dimension $p$0, and some primarily reduce invariant-measure bias at a fixed step size. The following summary collects representative non-asymptotic or empirical results exactly as stated in the cited works.

Method family	Representative result	Source
fMALA	complexity $p$1 versus $p$2 for standard MALA; asymptotical optimal acceptance $p$3	(Durmus et al., 2015)
SRK-LD	$p$4 iterations for strongly convex overdamped Langevin	(Li et al., 2019)
RKLMC-2G	$p$5 with two gradient evaluations per step under LSI and dissipativity	(Yang et al., 8 May 2026)
Third-order Langevin diffusion	$p$6 in the ridge-separable case	(Mou et al., 2019)
Picard–Lagrange $p$7-th order LMC	$p$8	(Mahajan et al., 21 Oct 2025)
$p$9-th order LMC algorithms	$$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$0, $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$1	(Dang et al., 24 Aug 2025)
Random splitting LMC	invariant-measure bias $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$2 in Wasserstein-1	(Li et al., 9 Oct 2025)
MCLMC	improvement factor $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$3 for modes and $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$4 for cosmological parameters at $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$5	(Bayer et al., 2023)

These rates are not directly comparable without attention to oracle model, regularity assumptions, and whether “cost” means gradient evaluations, full transitions, or asymptotic bias at fixed step size. fMALA is Metropolis-adjusted and exploits higher derivatives in the proposal; SRK-LD and RKLMC-2G are overdamped, unadjusted, and gradient-only in implementation; third-order and $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$6-th order lifted schemes change the continuous dynamics; microcanonical and randomized Nyström methods are non-reversible phase-space samplers (Durmus et al., 2015, Li et al., 2019, Yang et al., 8 May 2026, Mou et al., 2019, Mahajan et al., 21 Oct 2025, Bou-Rabee et al., 2023).

A recurring pattern is a trade-off between dimension and accuracy exponents. The Picard–Lagrange framework improves the $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$7-dependence as $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$8 increases but its dimension exponent tends to $$|\mathbb{E}\varphi(X_T)-\mathbb{E}\varphi(Y_N)|=O(h^p)$$9 (Mahajan et al., 21 Oct 2025). The $\varphi$0-th order Langevin algorithms similarly improve both exponents as $\varphi$1 grows, but under strong derivative-growth assumptions and with increasingly elaborate per-step Gaussian updates (Dang et al., 24 Aug 2025). By contrast, random splitting and microcanonical methods emphasize invariant-measure bias reduction or empirical effective sample size rather than explicit asymptotic minimax scaling (Li et al., 9 Oct 2025, Bayer et al., 2023).

6. Assumptions, misconceptions, and open directions

The strongest common limitation is regularity. High-order overdamped schemes typically assume Lipschitz gradient plus bounded or Lipschitz Hessian and third derivatives; the newer non-log-concave SRK analysis still requires dissipativity, gradient/Hessian/third-derivative Lipschitzness, and a Log-Sobolev inequality (Yang et al., 8 May 2026). aHOLA extends to superlinear drifts, but only through intricate polynomial taming and high moment assumptions such as $\varphi$2 (Neufeld et al., 2024). Lifted $\varphi$3-th order methods demand bounded high-order derivatives of $\varphi$4 up to order $\varphi$5, while the $\varphi$6-th order LMC algorithms impose a derivative-growth condition strong enough that polynomial approximation errors remain subordinate to splitting errors (Mahajan et al., 21 Oct 2025, Dang et al., 24 Aug 2025).

A second limitation concerns exactness. Many of these methods are unadjusted. HOLA, SRK-LD, RKLMC-2G, random splitting LMC, aHOLA, and the lifted higher-order algorithms all trade Metropolis correction for higher-order bias control (Sabanis et al., 2018, Li et al., 2019, Yang et al., 8 May 2026, Li et al., 9 Oct 2025, Neufeld et al., 2024). MCLMC goes further by removing accept–reject entirely and replacing exact canonical invariance with approximate microcanonical invariance regulated through energy fluctuations per dimension (Bayer et al., 2023). This suggests strong practical gains in some regimes, but it also means that finite-step bias must be analyzed or diagnosed rather than automatically eliminated.

Several recurring misconceptions are addressed explicitly in the literature. High-order need not mean higher derivatives: MCLMC uses only first-order derivatives, and RKLMC-2G is explicitly Hessian-free (Bayer et al., 2023, Yang et al., 8 May 2026). Conversely, higher derivatives need not imply higher weak order in the sampling sense: fMALA is designed through higher-order cancellation of the Metropolis ratio, but as an SDE integrator fULA still has weak order $\varphi$7 (Durmus et al., 2015). Likewise, a method may be high-order for invariant-measure bias yet not for pathwise strong approximation, as in random splitting LMC (Li et al., 9 Oct 2025).

Open problems recur across papers. The most prominent are extending high-order theory beyond strong log-concavity or LSI, reducing dependence on higher derivatives in lifted schemes, designing adaptive step-size and noise-control rules without Metropolis correction, and obtaining sharper dimension dependence for high-order methods in genuinely high-dimensional regimes (Yang et al., 8 May 2026, Mahajan et al., 21 Oct 2025, Neufeld et al., 2024). Another plausible implication is that multilevel Monte Carlo frameworks for invariant measures could provide an additional acceleration layer once a high-order integrator admits a suitable contractive coupling; the MLMC paper explicitly frames its theory as a backend for more advanced discretizations, including high-order ones (Giles et al., 2016).

Taken together, the field no longer treats Langevin Monte Carlo as synonymous with Euler-discretized overdamped diffusion. It has become a broader program in numerical stochastic analysis and MCMC design, spanning higher-order SDE integrators, lifted non-reversible dynamics, geometric structure preservation, and carefully engineered proposal mechanisms.