Joint work with Sonja Cox and Roy Schieven (University of Amsterdam). Paper: arXiv:2605.31088.
Two copies of the same diffusion, started apart. How can you get them to coalesce? For kinetic Langevin dynamics — the SDE behind the diffusive-to-ballistic speedup in modern sampling — the surprising answer is that the optimal coupling, with respect to the natural joint filtration, is provably non-Markovian: it strictly beats every Markovian strategy in TV. We call this the non-Markovian advantage.
This work was announced at the FIM workshop on Scalable MCMC Sampling at ETH Zurich, a terrific event organized by Yuansi Chen, Francesco Pedrotti, and Peter Whalley. The program drew researchers from probability, statistics, and machine learning, and the conversations across those communities made the week particularly rewarding. I also got a chance to walk around Lake Zurich and hike up Uetliberg.
1. The setting: sampling via kinetic Langevin
Fix a target distribution μ_target(dx) ∝ exp(-U(x)) dx on R^d, with U in C^1. Two natural diffusions leave μ_target invariant, possibly after marginalization.
Overdamped Langevin.
dX_t = -∇U(X_t) dt + √2 dW_t.
Kinetic Langevin. Lift to R^(2d) by introducing a velocity coordinate,
dX_t = V_t dt,
dV_t = -∇U(X_t) dt – γ V_t dt + √(2γ) dW_t,
which preserves the Boltzmann-Gibbs measure μ = μ_target ⊗ N(0, I_d).
Two structural facts are worth pausing on. First, noise enters only through the velocity coordinate; it reaches the position indirectly, via the kinematic constraint dX = V dt. This makes kinetic Langevin a hypoelliptic SDE: the noise covariance is degenerate, yet the law of (X_t, V_t) admits a smooth density in all 2d coordinates for t > 0. Second, kinetic Langevin can be viewed as a stochastic relaxation of Newton’s equations of motion in a potential U, augmented by friction γ V_t dt and noise √(2γ) dW_t in fluctuation-dissipation balance.
2. Why hypoellipticity is worth the trouble
The degeneracy is exactly what makes the dynamics powerful. Under a Poincaré inequality with constant m:
- Overdamped Langevin (Bakry, Gentil & Ledoux, 2014, Thm. 4.2.5) relaxes on a diffusive timescale T_relax ≍ m^(-1).
- Critically-tuned kinetic Langevin (Cao, Lu & Wang, 2023; Eberle & Lörler, 2024) relaxes on a ballistic timescale T_relax ≍ m^(-1/2).
A square-root speedup in ill-conditioned problems: the probabilistic analog of Nesterov acceleration in optimization. This is the motivation for studying kinetic Langevin and, more practically, its splitting discretizations.
3. From SDE to MCMC: the OBABO splitting
For sampling, we discretize. The kinetic Langevin generator admits a clean decomposition into three solvable flows:
- O_h, the Ornstein–Uhlenbeck step, exact in distribution;
- θ^(A)_h, the kinematic flow (X, V) ↦ (X + hV, V);
- θ^(B)_h, the velocity kick by -h∇U(X).
The OBABO splitting (Bussi & Parrinello, 2007; Leimkuhler & Matthews, 2013) assembles them symmetrically,
Φ^(OBABO)_h = O_(h/2) ∘ θ^(B)_(h/2) ∘ θ^(A)_h ∘ θ^(B)_(h/2) ∘ O_(h/2).
The composition is second-order weakly accurate. Each step is a deterministic map of two Gaussian increments ξ_k = (ξ_k^(1), ξ_k^(2)) ~ N(0, I_(2d)),
Z^h_(k+1) = Ψ^h_(Z^h_k)(ξ_(k+1)), Z^h_k = (X^h_k, V^h_k).
The central question is when the corresponding Markov chain (Z_h^k) realizes the acceleration. Writing κ = L/m for the condition number of the target (m·I ⪯ ∇²U ⪯ L·I), and t_mix(ν, ε) = inf{n : d_TV(ν P^n, μ) ≤ ε} for the TV mixing time from initial law ν, the paper is motivated by:
Question. For which initial laws ν, if any, does OBABO attain accelerated mixing t_mix(ν, ε) ≍ √κ · log(1/ε)?
4. Couplings and total variation
Two pieces of standard probabilistic machinery underlie our approach.
Coupling characterization of TV. The identity
d_TV(ν_1, ν_2) = inf P(X ≠ X̃),
where the infimum ranges over couplings (X, X̃) with X ~ ν_1 and X̃ ~ ν_2, reduces upper bounds on TV to the construction of a coupling and a bound on the disagreement probability.
Wasserstein-to-TV regularization. If, for every pair of point masses,
d_TV(δ_z P^n, δ_z̃ P^n) ≤ C_n · |z – z̃|,
then by joint convexity of TV along a W_1-optimal coupling,
d_TV(ν P^n, ν̃ P^n) ≤ C_n · W_1(ν, ν̃).
Combined with Wasserstein contraction in the chain, this transfers Wasserstein control to TV control, which is what one ultimately wants for mixing.
The problem now reads: build a coupling of two OBABO chains started at z and z̃ that brings their laws together in TV. Hypoellipticity makes this hard. The very degeneracy that yields the acceleration frustrates attempts to force two copies to meet.
5. A brief history of hypoelliptic couplings
The hypoelliptic coupling program has a long arc.
- Stochastic oscillator (Ben Arous, Cranston & Kendall, 1995). The first hypoelliptic coupling: a co-adapted switching between synchronous and antithetic Brownian drivers makes two copies of the stochastic oscillator coincide in finite time almost surely.
- Kolmogorov diffusion (W_t, ∫_0^t W_s ds) (Banerjee & Kendall, 2016). The first impossibility theorem: no Markovian coupling matches the asymptotic TV decay rate.
- Wasserstein contraction for kinetic Langevin (Eberle, Guillin & Zimmer, 2019). A hybrid coupling — synchronous on a contractive hyperplane in phase space, reflection transverse to it — gives quantitative Wasserstein convergence without requiring global convexity of U. Discrete-time analogues followed in work of Cheng et al. (2018, 2020), Dalalyan & Riou-Durand (2020), Leimkuhler-Paulin-Whalley (2024), and Schuh & Whalley (2025).
- TV mixing via Wasserstein-to-TV (Roberts–Rosenthal, 2002; Madras–Sezer, 2010; Monmarché, 2021; Gouraud et al., 2025; Chak–Monmarché, 2025). Hypoelliptic smoothing transfers Wasserstein bounds to TV bounds. Chak & Monmarché’s recent work constructs an explicit coalescence map Ψ^n_(z, z̃) from an ansatz, with the resulting Wasserstein-to-TV bound closed by Lemma 15 of B.-R. & Eberle (2023) – our starting point below.
6. Preliminaries
6.1 A general lemma: TV bound between a reference and a perturbed Gaussian
Closing a Wasserstein-to-TV bound for OBABO comes down to controlling the total variation distance d_TV between two laws on noise space — the reference Gaussian Law(ξ) and its pushforward Law(Ψⁿ_(z, z̃)(ξ)) under the coalescence map. The following sharp lemma of B.-R. & Eberle (2023) provides the required bound.
Lemma (B.-R. & Eberle, 2023; Lem. 15). Let ξ ~ N(0, I_d) and Φ : ℝ^d → ℝ^d be a C^1 diffeomorphism. Then
d_TV( Law(ξ), Law(Φ(ξ)) ) ≤ √(2 KL),
where
KL = E[ (1/2) |Φ(ξ) − ξ|^2 + tr(DΦ(ξ) − I) − log|det DΦ(ξ)| ].
Proof sketch. Pinsker’s inequality gives d_TV ≤ √(KL/2). A change of variables on the pushforward density yields
KL = E[ (1/2) |Φ(ξ) − ξ|^2 + (Φ(ξ) − ξ) · ξ − log|det DΦ(ξ)| ].
Gaussian integration by parts (Stein’s identity) handles the cross term:
E[ (Φ(ξ) − ξ) · ξ ] = E tr( DΦ(ξ) − I ),
since for any C^1 vector field F : ℝ^d → ℝ^d, E[ F(ξ) · ξ ] = E[ div F(ξ) ]. Substituting recovers the displayed KL. □
Remark (OBABO). For the coalescence map Ψⁿ_(z, z̃) considered next, the i-th coordinate ξ̃_i depends only on ξ_1, …, ξ_(i-1), so the Jacobian DΨⁿ_(z, z̃) is block lower-triangular with identity diagonal. Hence
tr( DΨⁿ_(z, z̃) − I ) ≡ 0, det DΨⁿ_(z, z̃) ≡ 1,
and the lemma collapses to a pure second-moment bound,
d_TV( Law(ξ), Law(Ψⁿ_(z, z̃)(ξ)) ) ≤ E[ |Ψⁿ_(z, z̃)(ξ) − ξ|^2 ]^(1/2).
Choosing the interior gap trajectory y_1, …, y_(n-1) to minimize this second moment is a classical minimum-energy LQ control problem, presented later in §8 — and the resulting bound is one of our main theorems in §7.1.
6.2 The coalescence map and its Malliavin derivative
The three results that follow share a common object — the coalescence map — and a common framing in terms of differentiation on noise space. We fix both here.
Write Ψⁿ_z : ℝ^(2dn) → ℝ^(2d) for the chain map of OBABO: given a noise sequence ξ = (ξ_1, …, ξ_n), the value Ψⁿ_z(ξ) is the n-step state started at z and driven by ξ. The coalescence map Ψⁿ_(z, z̃) : ℝ^(2dn) → ℝ^(2dn) is the corresponding noise transport: given ξ driving the chain started at z, the image ξ̃ := Ψⁿ_(z, z̃)(ξ) is the noise that, applied to the chain started at z̃, makes the two chains meet at time n,
Ψⁿ_(z̃)( Ψⁿ_(z, z̃)(ξ) ) = Ψⁿ_z(ξ).
The figure below records the construction via the gap trajectory y_k := z̃_k − z_k. The endpoints are fixed by y_0 = z̃ − z and y_n = 0; the interior values y_1, …, y_(n-1) are free, and chosen later by LQ optimization.
The natural framework for the analysis is Malliavin calculus. The chain map Ψⁿ_z is a smooth functional of the Gaussian noise on ℝ^(2dn), and the Jacobian DΨⁿ_(z, z̃) — its Malliavin derivative with respect to ξ — encodes how the coalescing noise depends on the driving noise. Two structural facts about OBABO make this Jacobian especially well-behaved. First, the i-th component ξ̃_i depends only on ξ_1, …, ξ_(i-1), so DΨⁿ_(z, z̃) is block lower-triangular. Second, the diagonal blocks are the identity. Consequently det DΨⁿ_(z, z̃) ≡ 1 and tr(DΨⁿ_(z, z̃) − I) ≡ 0, and the KL bound of Lemma 15 in B.-R. & Eberle (2023) collapses to a second-moment cost on the gap trajectory. Optimizing that cost gives the explicit non-Markovian coupling.
7. Three results
7.1 A quantitative TV bound for OBABO
Theorem (B.-R.–Cox–Schieven, 2026; Thm. 3.2). For γ, h > 0, n a positive integer, and U in C^2(R^d) with ∇U being L-Lipschitz, for all z, z̃ in R^(2d),
d_TV(π_n(δ_z), π_n(δ_z̃))
≤ γ^(-1/2) · [ 5·(hn)^(-3/2) + (12 + 5γ)·(hn)^(-1/2)
+ (1 + hn/(1 + γhn)) · L · (γh)^(1/2) · (1 – exp(-γh))^(-1/2) · (hn)^(1/2) ]
· |z̃ – z|.
A few features of the bound are worth noting:
- It holds for all h > 0 and all positive integers n — no hn ≤ 1 restriction and no h ≤ h_0 smallness condition.
- It only requires ∇U Lipschitz; the Hessian-Lipschitz assumption of Chak & Monmarché (2025) is not needed.
- It yields a Wasserstein-to-TV regularization for arbitrary initial laws ν, ν̃ via the convexity argument above.
- It is realized by an explicit non-Markovian coupling.
7.2 An impossibility theorem
The natural place to test optimality is a quadratic potential. Take U(x) = α|x|^2, α ≥ 0, in the overdamped regime γ^2 > 4α. The drift matrix has eigenvalues
λ_± = -γ/2 ± (1/2)·√(γ^2 – 4α), λ_- < λ_+ ≤ 0.
If the initial gap Δz = (Δx, Δv) lies in the λ_- eigenspace — that is, λ_- · Δx = Δv — then
d_TV(Law(Z_t), Law(Z̃_t)) ≲ exp(λ_- t) · |Δz|.
The question is whether a Markovian coupling can match this rate. The answer is no.
Theorem (B.-R.–Cox–Schieven, 2026; continuous time Thm. 4.3). Under γ^2 > 4α ≥ 0: if λ_- · Δx = Δv and Δx ≠ 0, then d_TV ≤ C · exp(λ_- t) · |Δz| (upper). For every Markovian coupling μ and every Δz,
μ(Z_t ≠ Z̃_t) ≥ c_μ · min( t^(-1/2), exp(λ_+ t) ) for t ≥ t_μ.
A discrete-time version yields, for every h > 0 and every Markovian μ_h,
μ_h(Z^h_k ≠ Z̃^h_k) ≥ c · min( c_(μ_h) · (hk+1)^(-1/2), c_(μ_h) · exp(λ_+ hk), h^(-1) · exp(λ_- hk) ).
Since λ_- < λ_+, the Markovian lower bound is strictly slower than the upper bound for non-Markovian couplings. The result extends Banerjee-Kendall (2016) from the Kolmogorov diffusion to kinetic Langevin, ruling out an entire class of strategies.
7.3 An exact meeting probability for iterated one-shot
The canonical Markovian candidate is the iterated one-shot coupling: at each step, maximize the meeting probability via a reflection coupling. It satisfies the now-equals-forever property (Z^h_s = Z̃^h_s implies Z^h_t = Z̃^h_t for t ≥ s, almost surely) and is asymptotically optimal for overdamped Euler-Maruyama (Durmus & Moulines, 2019). For kinetic Langevin, the impossibility result says it must be suboptimal — but it is natural to ask by how much.
Theorem (B.-R.–Cox–Schieven, 2026; Thm. 5.1). Let (Z_k, Z̃_k) be the iterated one-shot coupling of the linear chain Z_(k+1) = A_(k+1) Z_k + B_(k+1) ξ_(k+1) with A_k, B_k non-singular. Then
P(Z_n ≠ Z̃_n) = 2 · Φ( 1 / (2 · Θ_n^(1/2)) ) – 1,
where Θ_n = Σ_(k=1)^n [ 1 / |B_k^(-1) Π_k Δz|^2 ] and Π_k = A_k · A_(k-1) · ⋯ · A_1.
This sharpens Durmus–Moulines (2019, Thm. 19) from inequality to equality in the linear-drift case. It is used as a lower bound on P(Z_n ≠ Z̃_n), not a TV upper bound. For free kinetic Langevin (α = 0) with initial gap Δz = (Δx, -γ · Δx), the probability of not meeting degrades like h^(-1) · exp(λ_- hk) as the step size shrinks — exactly saturating the corresponding term in the impossibility lower bound.
8. The non-Markovian construction
The coupling that realizes the upper bound is non-Markovian.
- Sample ξ ~ N(0, I_(2dn)) — the noise driving the first chain over the entire horizon [0, n].
- Compute the proposal ξ̃* := Ψ^n_(z, z̃)(ξ) — the noise the second chain would need to coalesce with the first at time n.
- Maximally couple ξ and ξ̃*. On the acceptance event, set ξ̃ = ξ̃*: the chains meet at time n. On rejection, draw ξ̃ independently: the chains evolve independently.
By construction,
P(Z^h_n ≠ Z̃^h_n) = d_TV( Law(ξ), Law(Ψ^n_(z, z̃)(ξ)) ).
The coupling is non-Markovian because ξ̃ depends on the entire ξ at once. Intermediate chains generally disagree; the chains meet only at the terminal horizon, and generically never before.
Choosing Ψ optimally. The trajectory Ψ^n_(z, z̃) is chosen to minimize the TV cost. In the force-free case (∇U ≡ 0), Lemma 15 of B.-R. & Eberle (2023) reduces the TV bound to a controlled L^2 cost on the gap trajectory y_k := z̃_k – z_k,
d_TV ≤ (1/2) · ( Σ_(k=1)^n |E_k|^2 )^(1/2), E_(k+1) = -L_h^(-1) · (y_(k+1) – A_h · y_k),
with boundary conditions y_0 = z̃ – z and y_n = 0. This is a classical minimum-energy linear–quadratic control problem, solvable in closed form via the controllability Gramian Σ_(h,n):
Σ_(k=1)^n |E_k|^2 = | Σ_(h,n)^(-1/2) · A_h^n · Δz |^2.
For a general potential, the same trajectory is used and ∇U is handled as a Lipschitz perturbation — the design principle of Eberle–Guillin–Zimmer.
One technical point makes the OBABO bound clean. The Jacobian DΨ^n_(z, z̃) is block lower-triangular with identity diagonal, because the i-th coalescence noise depends only on ξ_1, …, ξ_(i-1). Hence det DΨ^n_(z, z̃) ≡ 1 and tr( DΨ^n_(z, z̃) – I ) ≡ 0, so the KL bound from Lemma 15 reduces to the second-moment term
d_TV ≤ E[ |Ψ^n_(z, z̃)(ξ) – ξ|^2 ]^(1/2).
9. The non-Markovian advantage
Taken together, the three theorems give a complete picture. The natural strategy for coupling kinetic Langevin in total variation — the iterated reflection coupling that works so cleanly for overdamped dynamics — provably cannot capture the sharp asymptotic rate. The coupling that does is global: it solves a classical minimum-energy control problem on the gap trajectory, and forces the chains to meet only at the terminal horizon. The very hypoellipticity that delivers the ballistic speedup is what makes the Markov property an obstruction.
Establishing the corresponding accelerated TV mixing bound t_mix(ν, ε) ≍ κ^(1/2) log(1/ε) reduces, via the Wasserstein-to-TV regularization of §6.1, to proving W_1 contraction for OBABO at rate κ^(-1/2). This remains open.
Related directions include KL and Rényi divergence bounds (B.-R., Mitra & Wibisono, 2026), other splittings (Schuh & Whalley, 2025), and Metropolis-adjusted variants (B.-R. & Oberdörster, 2024, EJP).
Acknowledgements
Many thanks to Pierre Monmarché, Andreas Eberle, and Stefan Oberdörster for fruitful discussions.
References
- Bakry, D., Gentil, I., & Ledoux, M. (2014). Analysis and Geometry of Markov Diffusion Operators. Grundlehren der mathematischen Wissenschaften, vol. 348. Springer. DOI
- Banerjee, S., & Kendall, W. S. (2016). Coupling the Kolmogorov diffusion: maximality and efficiency considerations. Advances in Applied Probability, 48(A), 15–35. DOI
- Ben Arous, G., Cranston, M., & Kendall, W. S. (1995). Coupling constructions for hypoelliptic diffusions: two examples. In Stochastic Analysis (Ithaca, NY, 1993), Proc. Sympos. Pure Math., vol. 57, 193–212. American Mathematical Society. DOI
- Bou-Rabee, N., Cox, S., & Schieven, R. (2026). On couplings for kinetic Langevin diffusions. arXiv preprint. arXiv:2605.31088
- Bou-Rabee, N., & Eberle, A. (2023). Mixing time guarantees for unadjusted Hamiltonian Monte Carlo. Bernoulli, 29(1), 75–104. DOI
- Bou-Rabee, N., Mitra, S., & Wibisono, A. (2026). Tail-sensitive KL and Rényi convergence of unadjusted Hamiltonian Monte Carlo via one-shot couplings. arXiv preprint. arXiv:2601.09019
- Bou-Rabee, N., & Oberdörster, S. (2024). Mixing of Metropolis-adjusted Markov chains via couplings: the high acceptance regime. Electronic Journal of Probability, 29, Paper No. 89. DOI
- Bussi, G., & Parrinello, M. (2007). Accurate sampling using Langevin dynamics. Physical Review E, 75(5), 056707. DOI
- Cao, Y., Lu, J., & Wang, L. (2023). On explicit L²-convergence rate estimate for underdamped Langevin dynamics. Archive for Rational Mechanics and Analysis, 247, Paper No. 90. DOI
- Chak, M., & Monmarché, P. (2025). Reflection coupling for unadjusted generalized Hamiltonian Monte Carlo in the nonconvex stochastic gradient case. IMA Journal of Numerical Analysis, draf045. DOI
- Cheng, X., Chatterji, N. S., Bartlett, P. L., & Jordan, M. I. (2018). Underdamped Langevin MCMC: A non-asymptotic analysis. In Conference on Learning Theory (COLT), 300–323. PMLR.
- Cheng, X., Chatterji, N. S., Abbasi-Yadkori, Y., Bartlett, P. L., & Jordan, M. I. (2020). Sharp convergence rates for Langevin dynamics in the nonconvex setting. arXiv preprint. arXiv:1805.01648
- Dalalyan, A. S., & Riou-Durand, L. (2020). On sampling from a log-concave density using kinetic Langevin diffusions. Bernoulli, 26(3), 1956–1988. DOI
- Durmus, A., & Moulines, É. (2019). High-dimensional Bayesian inference via the unadjusted Langevin algorithm. Bernoulli, 25(4A), 2854–2882. DOI
- Eberle, A., Guillin, A., & Zimmer, R. (2019). Couplings and quantitative contraction rates for Langevin dynamics. Annals of Probability, 47(4), 1982–2010. DOI
- Eberle, A., & Lörler, F. (2024). Non-reversible lifts of reversible diffusion processes and relaxation times. Probability Theory and Related Fields. DOI
- Gouraud, N., Le Bris, P., Majka, A., & Monmarché, P. (2025). HMC and underdamped Langevin united in the unadjusted convex smooth case. SIAM/ASA Journal on Uncertainty Quantification, 13(1), 278–303. DOI
- Leimkuhler, B., & Matthews, C. (2013). Rational construction of stochastic numerical methods for molecular sampling. Applied Mathematics Research Express. AMRX, 2013(1), 34–56. DOI
- Leimkuhler, B. J., Paulin, D., & Whalley, P. A. (2024). Contraction and convergence rates for discretized kinetic Langevin dynamics. SIAM Journal on Numerical Analysis, 62(3), 1226–1258. DOI
- Madras, N., & Sezer, D. (2010). Quantitative bounds for Markov chain convergence: Wasserstein and total variation distances. Bernoulli, 16(3), 882–908. JSTOR
- Monmarché, P. (2021). High-dimensional MCMC with a standard splitting scheme for the underdamped Langevin diffusion. Electronic Journal of Statistics, 15(2), 4117–4166. DOI
- Roberts, G. O., & Rosenthal, J. S. (2002). One-shot coupling for certain stochastic recursive sequences. Stochastic Processes and their Applications, 99(2), 195–208. DOI
- Schuh, K., & Whalley, P. A. (2025). Convergence of kinetic Langevin samplers for non-convex potentials. arXiv preprint. arXiv:2405.09992