Do deterministic midpoint methods outperform randomized ones for Langevin sampling?

Deterministic quadrature versus randomized midpoint schemes

The paper compares two families of discretizations of Langevin dynamics for sampling from a target π ∝ e^{-V} on R^d, where the Hessian satisfies 0 ≺ αI ⪯ ∇²V(x) ⪯ βI and κ := β/α. No bounds on third or higher derivatives of V are assumed. Accuracy is measured by the number of gradient evaluations needed to produce a distribution μ with √α W₂(μ, π) ≤ ε, where W₂ is the 2-Wasserstein distance.

On one side sit deterministic quadrature rules: Heun for overdamped dynamics and the underdamped exponential midpoint method. On the other sit randomized midpoint (RLMC and RULMC) and Poisson midpoint (PLMC and PULMC) schemes. The paper's question is whether the randomization that earlier analyses credited with improving dimension and accuracy exponents actually helps once sharper tools are used.

Heun reaches d^{1/3} ε^{-2/3} gradient queries

For the deterministic Heun method, the paper proves that if 0 < βh ≤ h₀ then the kernel has a unique invariant distribution π̂_h^{Heun} and α W₂²(π̂_h^{Heun}, π) ≤ C κ d (βh)³. Iterating from the minimizer x⋆, for 0 < ε ≤ √d/2 the number of gradient queries sufficient for the accuracy criterion is at most C κ^{4/3} d^{1/3} ε^{-2/3} log(d/ε²).

The proof does not rely on the standard coupling argument, which would only give κ^{-1} W₂²(π̂_h^{Heun}, π) ≲ κ d h. Instead it separates the conditional mean of the one-step error from its centered fluctuation. The systematic part is a trapezoidal quadrature defect for the semigroup, and because it is a deterministic function of the current position, fresh noise does not damp it. The authors cancel it in advance by shifting the comparison point using a corrector χ that solves the discrete Poisson equation (I − P_h)χ = e_h, giving a reference distribution ν := (I + χ)_#π with W₂(ν, π) ≤ C √d h³. At this shifted reference the conditional mean of the error is only O(√d h⁵), while the centered fluctuation, though of size O(√d h³), costs only O(√d h⁵) in W₂ because a centered perturbation affects a Gaussian convolution to second order.

Exponential midpoint cuts the dimension exponent to d^{1/4}

For the underdamped exponential midpoint method with γ = 2√β, the paper proves that for 0 < √β h ≤ c κ^{-1/4} the kernel has a unique invariant distribution π_h^{ExpMid} in P₂(R^{2d}) and √α W₂((π_h^{ExpMid})_x, π) ≤ C √(κd) (√β h)². Starting from X̂₀ = x⋆ and P̂₀ ∼ N(0, I_d), the iterate satisfies √α W₂(Law(X̂_n), π) ≤ C √d e^{-cn√β h/κ} + C √(κd) (√β h)². Consequently, for 0 < ε ≤ √d/2, the gradient-query count satisfies N_{∇V} ≤ C κ^{5/4} d^{1/4} ε^{-1/2} log(d/ε²).

The mechanism mirrors the overdamped case. The exact solution expresses the endpoint through two gradient integrals, and the method approximates each by h times its integrand at t = h/2 using a predictor X^{pred}_{h/2}. The predictor differs from the exact midpoint by O(√d h³) in L², and the gradient enters with a coefficient of order h, so replacing the predictor by the exact midpoint changes the output by only O(√d h⁴). Contraction is measured in a twisted hypocoercive metric where friction acts on momentum and a cross term transfers dissipation to position; the one-step error at stationarity carries the full order h³, and a fixed-point argument divides it by the contraction gap to bound the invariant distribution.

Lower bounds match the d and ε powers at bounded κ

At bounded condition number, the paper proves a lower bound that matches the d and ε powers of both deterministic methods. For Heun, the lower bound contains κ d^{1/3} ε^{-2/3} for d^{-1} ≤ ε ≤ ε₀, so the exponents of d and ε agree with the upper bound while a factor κ^{1/3} remains. The authors attribute this remaining gap to accumulating the distance between consecutive iterates by the triangle inequality; improving that accumulation would require further cancellation across steps, which the present proof does not establish.

The same kind of condition-number discrepancy appears in the underdamped case, where the authors state it cannot easily be rectified in the current proof. Readers comparing the two deterministic methods should therefore treat the κ exponents as upper-bound artifacts rather than as settled separations.

Randomized midpoint methods face a d^{1/3} ε^{-1/3} barrier

For the randomized midpoint methods and Poisson midpoint with at least two grid points, in both overdamped and underdamped variants, a simple Gaussian calculation yields a lower bound of d^{1/3} ε^{-1/3} to obtain an ε-close sample, even when starting from a benign initialization. The calculation uses a scale-Gaussian decomposition: the conditional state covariance is represented as an independent Gaussian multiplied by a random scalar or matrix, and a Wasserstein lower bound for Gaussian scale perturbations converts the variance of that scale into a distance from the standard Gaussian.

For the overdamped randomized midpoint method, the leading centered covariance defect is h²(1 − 2Θ) for one uniform evaluation time, and the invariant distribution satisfies W₂(π̂_h^{RLMC}, N(0, I_d)) ≥ c min{d h³, √(d h³)}. For the underdamped Poisson midpoint with k ≥ 2, the variance coefficient ϑ_k := (k−1)²(2k−1)/(64k³) lies in [3/512, 1/32), and the same machinery produces the corresponding lower bound. The contrast with the deterministic upper bounds is what the paper calls surprising: in high dimensions, deterministic discretizations can outperform their random counterparts.

Stationarity cancellation and Malliavin smoothing drive the gain

The improvement over previous upper bounds based on standard couplings comes from two analytic ingredients. The first is cancellation at stationarity: the proofs exploit the fact that the target is invariant, so the stationary path X_t has distribution π for every t, and the systematic error can be isolated as a quadrature defect for the semigroup rather than accumulated step by step. The second is smoothing using techniques from Malliavin calculus, which the appendix uses to verify the bounds; the convolution estimate shows that a centered perturbation affects a Gaussian convolution to second order.

These tools are what allow the analysis to recover two further powers of h beyond the elementary displacement estimate. The paper also notes that the recently proposed Poisson midpoint method obtained accuracy rates ε^{-1/3} exceeding the previously conjectured exponent, and that LMC was shown to match the classical randomized overdamped guarantee; both results revealed cancellations missed by conventional bounds and reopened the comparison between deterministic and randomized quadrature.

Where the deterministic advantage actually applies

The comparison is stated for a specific setting: a target π ∝ e^{-V} with 0 ≺ αI ⪯ ∇²V ⪯ βI, no bound on third or higher derivatives, accuracy measured as √α W₂ ≤ ε, and query counts counted as gradient evaluations. The deterministic upper bounds hold under step-size restrictions — 0 < βh ≤ h₀ for Heun and 0 < √β h ≤ c κ^{-1/4} for exponential midpoint — and the randomized lower bound applies to randomized midpoint and Poisson midpoint with at least two grid points, in both overdamped and underdamped variants, even from a benign initialization.

The practical reading is a dimensional one. When d is large and the target is strongly log-concave with bounded condition number, the deterministic methods carry the better d and ε exponents: Heun at d^{1/3} ε^{-2/3} and exponential midpoint at d^{1/4} ε^{-1/2}, against the randomized barrier of d^{1/3} ε^{-1/3}. The condition-number exponents are where the deterministic upper bounds are not yet tight, so the choice between Heun and exponential midpoint should be read through the d and ε powers first, and the κ factors second.

Sources used in this answer