The Picard-Lagrange Framework for Higher-Order Langevin Monte Carlo

TL;DR

This work introduces a Picard–Lagrange discretization for general $K$-th order Langevin dynamics to sample from smooth, strongly log-concave targets using only a gradient oracle. By coupling higher-order dynamics with a Lagrange polynomial drift surrogate and solving the induced implicit updates via Picard iterations, the authors establish non-asymptotic convergence in Wasserstein distance and derive gradient-query complexities that improve as the order $K$ grows. The main result shows that for $K\ge 3$, the required gradient evaluations to attain $W_2\le\varepsilon$ scale as $\tilde{O}\big(d^{\frac{K-1}{2K-3}}\varepsilon^{-\frac{2}{2K-3}}\big)$, with the convergence rate improving with $K$ and approaching $d^{1/2}$-like behavior for large $K$. The framework remains within the standard first-order gradient oracle model, offering a principled, scalable path to higher-order Langevin Monte Carlo methods and sharpening the theoretical understanding of discretization benefits in high-dimensional sampling.

Abstract

Sampling from log-concave distributions is a central problem in statistics and machine learning. Prior work establishes theoretical guarantees for Langevin Monte Carlo algorithm based on overdamped and underdamped Langevin dynamics and, more recently, some third-order variants. In this paper, we introduce a new sampling algorithm built on a general $K$th-order Langevin dynamics, extending beyond second- and third-order methods. To discretize the $K$th-order dynamics, we approximate the drift induced by the potential via Lagrange interpolation and refine the node values at the interpolation points using Picard-iteration corrections, yielding a flexible scheme that fully utilizes the acceleration of higher-order Langevin dynamics. For targets with smooth, strongly log-concave densities, we prove dimension-dependent convergence in Wasserstein distance: the sampler achieves $\varepsilon$-accuracy within $\widetilde O(d^{\frac{K-1}{2K-3}}\varepsilon^{-\frac{2}{2K-3}})$ gradient evaluations for $K \ge 3$. To our best knowledge, this is the first sampling algorithm achieving such query complexity. The rate improves with the order $K$ increases, yielding better rates than existing first to third-order approaches.

Figures (2)

• Figure 1: Schematic representation of Picard iterations. The procedure begins with constant values, which are used to form a Lagrange interpolant. Applying the update operator $\widehat{\mathcal{T}}$ at the interpolation nodes yields the next set of values, which then feed into the following Picard iteration.
• Figure 2: An illustration for one–step error decomposition: Starting from the current iterate $\widehat{X}(nh)$, (a) runs the continuous process (equivalently, the fixed–point flow of $\mathcal{T}$), yielding $\widehat{X}(nh) \to X^\ast_{\widehat{X}(nh)}$; (b) runs the fixed–point flow of the discretized operator $\widehat{\mathcal{T}}$, yielding $\widehat{X}(nh) \to \widehat{X}^\ast_{\widehat{X}(nh)}$; (c) is the algorithmic update, $\widehat{X}(nh) \to \widehat{X}((n+1)h)$. The three labeled distances match the terms in \ref{['eq:one-step-ED']}: (I) contraction of continuous process, (II) interpolation error $\|X^\ast_{\widehat{X}(nh)}(h)-\widehat{X}^\ast_{\widehat{X}(nh)}(h)\|_{S}^{2}$, (III) Picard convergence error $\|\widehat{X}^\ast_{\widehat{X}(nh)}(h)-\widehat{X}((n+1)h)\|_{S}^{2}$. Hence $E_{n+1}\le\text{(I)}+\text{(II)}+\text{(III)}$.

Theorems & Definitions (77)

• Definition 1: Tensor spectral norm
• Theorem 1
• Remark 5.1
• Proposition 1: One–step error bound
• Proposition 2
• Proposition 3: Interpolation error within one step
• Lemma 4
• Remark 6.1
• Remark 6.2
• Proposition 5: Picard iteration error