Accelerated Mixing Time of Randomized Hamiltonian Monte Carlo Siddharth Mitra∗
Vishwak Srinivasan†
Xiuyuan Wang‡
Andre Wibisono§
July 15, 2026
arXiv:2607.12902v1 [stat.ML] 14 Jul 2026
Abstract We show the Randomized Hamiltonian Monte Carlo (RHMC) algorithm has accelerated mixing time guarantees for sampling from log-concave probability distributions. RHMC proceeds by repeatedly simulating the continuous-time Hamiltonian dynamics for some random integration times, and resetting the velocity to be an independent Gaussian random variable between each simulation. We show that when the target distribution is log-concave and satisfies an αTalagrand inequality (for example, if the target distribution is α-strongly log-concave), if we use a random integration time from either the triangular or the exponential distribution with mean Θ(α−1/2 ), then RHMC converges exponentially fast in KL divergence, and the total integration time to reach error ε in KL divergence scales as O(α−1/2 log(ε−1 )). We also show that when the target distribution is log-concave, if we use a sequence of random integration times from the triangular distribution with exponentially increasing means, then the total integration time to reach error ε in KL divergence scales as O(ε−1/2 ). Our analysis relies on a bound on the average KL divergence along Hamiltonian dynamics, which is inspired by an analogous result on accelerated optimization methods based on Hamiltonian dynamics.
∗
Department of Computer Science, Yale University. Email: [email protected]. Department of Electrical Engineering and Computer Science, MIT. Email: [email protected]. ‡ Department of Computer Science, Yale University. Email: [email protected]. § Department of Computer Science, Yale University. Email: [email protected]. This work was supported by NSF awards CCF-2403391 and CAREER CCF-2443097. †
1
Contents 1 Introduction 1.1 Randomized Hamiltonian Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Our contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Background 2.1 Probability distributions and statistical divergences . . . . . . . . . . . . . . . . . . . 3 Main results 3.1 Convergence result under Talagrand inequality and semi log-concavity . . . . . . . . 3.2 Convergence result under log-concavity . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Proof technique: Average KL divergence along Hamiltonian dynamics . . . . . . . . 4 Proofs of the main results 4.1 Exponential distribution as a mixture of triangular distribution . . . . . . . . . . . . 4.2 Proof of Theorem 1 and Corollary 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Proof of Theorem 2 and Corollary 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Discussion A Additional related work B Connections to Hamiltonian Flow for optimization C Additional preliminaries C.1 General notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 Further discussion of probability distributions and statistical distances . . . . . . . . C.3 Convergence of probability distributions . . . . . . . . . . . . . . . . . . . . . . . . . D Details on properties of the Hamiltonian flow D.1 Deterministic properties of the Hamiltonian flow . . . . . . . . . . . . . . . . . . . . D.2 Distributional properties of the Hamiltonian flow . . . . . . . . . . . . . . . . . . . . D.3 Moments along Hamiltonian flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E Properties of the Hamiltonian flow under regular initialization E.1 Propagation of warmness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.2 Implications of regular initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . F Bound on average KL divergence under regular initialization F.1 Formula for the acceleration field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . F.2 Formula for the upper second time derivative of Wasserstein distance . . . . . . . . . F.3 A bound on the second time derivative of the Wasserstein distance . . . . . . . . . . F.4 Bound on average KL divergence under regularity assumption . . . . . . . . . . . . . G Approximation argument for regularity of initial distribution G.1 Approximation of the initial distribution . . . . . . . . . . . . . . . . . . . . . . . . . G.2 Finite-time Wasserstein stability of the Hamiltonian flow . . . . . . . . . . . . . . . . G.3 Lower semicontinuity of the integral functional . . . . . . . . . . . . . . . . . . . . . G.4 Convergence of the initial functional . . . . . . . . . . . . . . . . . . . . . . . . . . . G.5 Proof of the bound on average KL divergence along Hamiltonian flow . . . . . . . .
2
3 5 7 7 7 10 10 11 12 13 13 14 17 20 22 24 27 27 28 30 30 31 34 38 41 42 42 50 50 51 53 60 61 61 67 68 69 69
1
Introduction
Drawing samples from a probability distribution is an essential algorithmic task in many areas of science and engineering. A popular class of techniques for performing sampling is by running an ergodic Markov chain whose stationary distribution is the target probability distribution; such techniques are referred to as Markov Chain Monte Carlo (MCMC) [BGJM11]. The efficiency of an MCMC algorithm is dictated by its mixing time — this quantifies how quickly the Markov chain converges to the stationary distribution. In practice, many sampling algorithms are based on timediscretizations of continuous-time dynamics, including the overdamped or underdamped Langevin dynamics, which are stochastic processes, and the Hamiltonian dynamics, which is deterministic. In this paper, we focus on the Randomized Hamiltonian Monte Carlo (RHMC) algorithm, which is based on simulating the Hamiltonian dynamics for times that are randomized, and show that it achieves accelerated mixing time guarantees. To contextualize our results for RHMC, we briefly review relevant guarantees for the overdamped and underdamped Langevin dynamics, which have been studied in more detail. We denote the target distribution that we aim to draw samples from by ν X . The overdamped Langevin dynamics (LD) [RT96] is a basic stochastic process for sampling ν X by taking into account the gradient of the log density of ν X . Langevin dynamics has a natural optimization interpretation as the gradient flow dynamics (in the space of probability measures equipped with the Wasserstein metric) to minimize the relative entropy or Kullback-Leibler (KL) divergence with respect to the target distribution [JKO98]. Partly motivated by this optimization perspective, there have been many works that result in a rich understanding of the mixing time of LD and its algorithmic implementations under various structural assumptions on the target distribution. The mixing times discussed here depend on (a) the measure of discrepancy, and (b) the error level ε to the target. We highlight two key results that are particularly relevant to our work: • when ν X satisfies an α-log-Sobolev inequality (see Section 2.1.3 for a definition), the KL-mixing time for error ε of LD scales as α−1 log(ε−1 ) [OV00, Proof of Lemma 3], and • when ν X is log-concave (see Section 2.1.2 for a definition), the KL-mixing time for error ε scales as ε−1 [OV01, Corollary 2]. Notably, these rates match the time complexities of the Euclidean gradient flow for minimizing a function f that satisfies α-gradient domination (implied by α-strong convexity) and convexity, respectively, to within ε of the optimal value. There are also guarantees for time discretizations of LD, which we do not discuss further since in this paper we focus on the continuous-time complexity. Motivated by the algorithmic applications, as well as taking inspiration from the theory of acceleration in convex optimization, there have been several efforts made towards identifying other dynamics that mix faster. To be more specific, we seek analogues of the following results from convex optimization for the accelerated gradient flow [SBC16, WRJ21] for the time complexities to within ε of than the optimal value: • for a function that is α-strongly convex, the time complexity scales as α−1/2 log(ε−1 ), and • for a function that is convex, the time complexity scales as ε−1/2 . A candidate dynamics is the underdamped Langevin dynamics (ULD) which is a stochastic process defined on the phase space of both position and velocity. In ULD, the evolution of the 3
position is deterministically governed by the velocity, and the evolution of the velocity is encoded as a stochastic process involving the gradient of the log density of ν X . The ULD is a natural dynamics to consider primarily owing to the structural similarities to the accelerated gradient flow referenced earlier, and many works have studied whether ULD has faster convergence guarantees compared to LD (see for instance, [CCBJ18, MCC+ 21]). Indeed ULD has been shown to have an accelerated convergence guarantee; in particular, [CLW23] show ULD has a χ2 -mixing time for error ε that scales as α−1/2 log(ε−1 ) when ν X is log-concave and satisfies an α-Poincaré inequality (implied when ν X satisfies α-log Sobolev inequality). [Lu26] recently shows that ULD also has KL-mixing time for error ε that scales as α−1/2 log(ε−1 ) when ν X is both log-concave and satisfies an α-log-Sobolev inequality. This strictly improves on the KL-mixing time scaling of LD (while additionally assuming log-concavity), and also matches the accelerated convergence rate we expect from optimization. When ν X is only log-concave, [ACZ26] show that ULD has a KL-mixing time for error ε that scales as ε−1 , which is the same scaling achieved by LD; the accelerated rate of ε−1/2 seems unknown for ULD. The Hamiltonian dynamics is a deterministic process on the phase space of position and velocity that conserves the Hamiltonian or energy function. The Hamiltonian Monte Carlo (HMC) algorithm [DKPR87] is a piecewise-deterministic Markov chain that proceeds iteratively by (a) simulating the Hamiltonian dynamics for some integration time, and (b) resetting the velocity to be an independent Gaussian random variable. HMC and its variants are some of the most widely used sampling algorithms in practice, and underlie probabilistic programming systems such as Stan [CGH+ 17] and PyMC [PHF10]. Despite its practical importance, the theoretical guarantees of HMC and its variants are less developed than LD or ULD. Many existing results focus on the short integration time, where the time to simulate the Hamiltonian flow between velocity resetting scales inversely with the smoothness of the target distribution. With short integration time, HMC has unaccelerated convergence guarantees similar to LD and with additional dependence on the smoothness. [CV22] show that when ν X is strongly log-concave (see Section 2.1.2 for a definition) and L-log smooth, the W2 -mixing time of HMC (i.e., iterations) for error ε scales as Lα−1 log(ε−1 ) where each iteration simulates the Hamiltonian dynamics for integration time L−1/2 (the “short” integration time) [CV22, Theorem 1.3] and thus implies a total integration time that scales as L1/2 α−1 log(ε−1 ). This result is also shown to be tight via a matching lower bound. Later work by [Mon24] improve this result by generalizing to target distributions that satisfy a α-log-Sobolev inequality and are L-log-smooth, and providing a KL-mixing time of the same order. Hence, additional ideas are required to obtain faster convergence rates for HMC; in fact, not only are deterministic integration times unable to obtain accelerated mixing time guarantees [CV22, Theorem 1.4], deterministic and long integration times that scale inversely with the curvature lower bound as α−1/2 can fail to make any progress at all, as this integration time coincides with the natural oscillation frequencies of the dynamics, e.g., even for Gaussian target distributions. This motivates the Randomized Hamiltonian Monte Carlo (RHMC) [BRSS17], which randomizes the integration times between velocity refreshments. RHMC with suitably randomized integration times are conjectured to achieve accelerated convergence guarantees, see e.g., [Jia23]. Toward this conjecture, [LW22] show that when the target distribution is log-concave and satisfies α-Poincaré inequality, RHMC reaches error ε in chi-square divergence in an expected total integration time that scales as α−1/2 log(ε−1 ). A recent work by [MW26], extending the technique of [Lu26, LL26], show that when the target distribution is log-concave and satisfies an α-log-Sobolev inequality, RHMC reaches error ε in KL divergence in an expected total integration time α−1/2 log(ε−1 ), which
4
matches the desired accelerated rate from optimization. In this work, we establish accelerated mixing-time guarantees for RHMC in KL divergence in two settings: (a) when the target distribution is semi-log-concave and satisfies Talagrand inequality, we obtain a result which matches [MW26]; and (b) when the target distribution is log-concave, we obtain a new result. The accelerated mixing time guarantees for RHMC in all of these works also overcome the dependence on the smoothness parameter which is present in works studying deterministic integration times. We provide additional discussion on related work in Section A. In addition to the RHMC algorithm presented in [BRSS17] and studied in [LW22, MW26], which considers the integration time to be drawn from an exponential distribution, we also consider RHMC with integration time drawn from a triangular distribution; this has the benefit of being compactly supported, and therefore guarantees on the total integration time hold deterministically instead of in expectation. We describe these RHMC algorithms in Section 1.1 and describe our results in more detail in Section 1.2.
1.1
Randomized Hamiltonian Monte Carlo
We introduce the Randomized Hamiltonian Monte Carlo (RHMC) algorithm studied in this paper; see Algorithm 1 below. Our goal is to sample from a target distribution ν X with full support on Rd . We assume the target distribution has density function ν X ∝ e−f for some potential function f : Rd → R. The RHMC algorithm is based on the Hamiltonian dynamics. Consider the phase space R2d = Rd × Rd , which consists of the joint position and velocity variables. We define the Hamiltonian function H : R2d → R by, for all (x, y) ∈ R2d , 1 (1) H(x, y) := f (x) + ∥y∥2 . 2 Given the Hamiltonian function Eq. (1), the Hamiltonian dynamics or Hamiltonian flow is the evolution of the joint variables (Xt , Yt ) ∈ R2d following the system of ordinary differential equations: Ẋt = ∇y H(Xt , Yt ) = Yt Ẏt = −∇x H(Xt , Yt ) = −∇f (Xt )
(HF)
starting from any initial configuration (X0 , Y0 ) ∈ R2d at time t = 0. The RHMC algorithm proceeds by evolving along the dynamics (HF) combined with periodic velocity refreshment from the standard Gaussian distribution N (0, Id ). We describe the RHMC algorithm formally in Algorithm 1. Algorithm 1 Randomized Hamiltonian Monte Carlo (RHMC) Require: Number of iterations K ∈ N ; Integration time distributions D1 , D2 , . . . , DK ; Initial distribution ρX 0 . Draw initial iterate x0 ∼ ρX 0 . for k = 1, . . . , K do (k) (k) Set X0 = xk−1 , and draw Y0 ∼ N (0, Id ) independently. Draw integration time τk ∼ Dk . (k) (k) (k) (k) Run Hamiltonian flow (HF) from (X0 , Y0 ) for time τk to obtain (Xτk , Yτk ). (k) Set next iterate xk = Xτk which is a random variable xk ∼ ρX k . end for return xK which is a random variable xK ∼ ρX K. 5
Randomized Hamiltonian Monte Carlo (RHMC), presented in Algorithm 1, corresponds to a suite of algorithms, with different choices of integration time distributions D1 , . . . , DK corresponding to different algorithms. When Dk = δT for some T > 0, or equivalently, τk = T is deterministic for all k ∈ {1, . . . , K}, the algorithm is referred to as Hamiltonian Monte Carlo (HMC); see Section A for a discussion of works studying HMC. We study the RHMC algorithm where the integration time distributions D1 , . . . , DK are not degenerate, as studied in [Mac89, CLS07, Nea11, BRSS17]. We study Algorithm 1 with the following choices of integration time distributions: [Set-Tri]
In the setting [Set-Tri], we study Algorithm 1 with Dk := TriT with T > 0 for all k ∈ {1, . . . , K} . Here TriT is the triangular distribution with parameter T > 0, which is supported on [0, T ] with density function at t ∈ [0, T ] given by: TriT (t) :=
2 (T − t) . T2
(2)
The choice of triangular integration time distributions in [Set-Tri] is motivated by studying the analogue of Algorithm 1 for optimization; see Section B for a discussion. We also note that the choice of triangular integration time is similar to the adaptive stopping time of Hamiltonian dynamics-based No-U-Turn Sampler (NUTS) which is widely used in practice [HG14]. [Set-Exp]
In the setting [Set-Exp], we study Algorithm 1 with Dk := Exp1/T with T > 0 for all k ∈ {1, . . . , K} . Here Exp1/T is the exponential distribution supported on [0, ∞) with parameter T1 > 0, with density function given by, for all t ∈ [0, ∞), Exp1/T (t) :=
t 1 exp − . T T
(3)
The mean of Exp1/T is T , i.e., Eτ ∼Exp1/T [τ ] = T . We will use the fact that Exp1/T can be represented as a mixture of triangular distributions; see Lemma 4. The RHMC algorithm with exponentially distributed integration times was proposed in [BRSS17]. f [Set-Tri]
f we study Algorithm 1 with In the setting [Set-Tri], f T for some Tk > 0 for all k ∈ {1, . . . , K} . Dk := Tri k f T refers to the endpoint-biased triangular distribution given by Here for any T > 0, Tri f T := 1 δT + 1 TriT . Tri 2 2 Like [Set-Tri], the choice of the integration time distribution in this setting is motivated by studying the analogue of Algorithm 1 for convex optimization; see Section B for more details.
6
1.2
Our contributions
In this work, we prove convergence guarantees for the idealized RHMC algorithm in Algorithm 1 under several assumptions on the target distribution. Here idealized means we assume we can simulate the Hamiltonian flow (HF) exactly. We measure the complexity of RHMC by the total amount of integration time required to simulate the Hamiltonian flow in order to reach a prescribed error in KL divergence. Our first result concerns target distributions ν X that satisfy α-Talagrand inequality (see Section 2.1.3) and are M -semi-log-concave (see Section 2.1.2) for some 0 ≤ M < α. In this setting, we prove a convergence guarantee in KL divergence for RHMC with integration times drawn from either the triangular or exponential distribution with mean Θ((α − M )−1/2 ). As we show in TheX orem 1, to obtain KL(ρX K ∥ ν ) ≤ ε, it suffices to simulate Hamiltonian flow for a total integration time X KL(ρX 1 0 ∥ν ) . log O √ ε α−M When using [Set-Tri], this upper bound is deterministic. When using [Set-Exp], this is an upper bound on the expected total integration time. In particular, when M = 0 (i.e., when ν X is log-concave), this gives a total integration time O(α−1/2 log(ε−1 )), which improves over the O(α−1 log(ε−1 )) time of the overdamped Langevin dynamics in the same setting, and matches the result of [MW26]. Our second result concerns the log-concave case, without assuming Talagrand inequality. In this f T = 1 δT + 1 TriT , with integrationsetting, we use an endpoint-biased triangular distribution Tri k k 2 k 2 time parameter Tk increasing with iteration k. We show in Theorem 2 a convergence guarantee in KL divergence, which implies (see Corollary 2) that to obtain error ε in KL divergence, it suffices to simulate Hamiltonian flow for total integration time s 1 2 X X X X KL(ρ0 ∥ ν ) + 3 W2 (ρ0 , ν ) . O ε This improves the O(ε−1 ) time of the overdamped Langevin dynamics in this setting, and matches what we expect from the theory of accelerated convex optimization. Our analysis is inspired by translating the proofs from the Hamiltonian dynamics-based optimization method of [WSF+ 26] to the sampling setting, in particular to the space of probability distributions with the Wasserstein geometry; see Section B for further discussion. In particular, although our end result in the Talagrand case recovers the recent result of [MW26], our proof technique is different, and has a clear optimization analogue.
2
Background
We briefly recall relevant definitions we use in this paper. We provide more details in Section C.
2.1
Probability distributions and statistical divergences
Let P(Rd ) denote the space of all probability distributions on Rd , and P2 (Rd ) denote the space of probability distributions on Rd with finite second moment, so Eρ [∥X∥2 ] < ∞ for all ρ ∈ P2 (Rd ). Let P2,ac (Rd ) denote the subspace of P2 (Rd ) consisting of probability distributions which are absolutely 7
continuous with respect to the Lebesgue measure dx on Rd . We identify a probability distribution ρ ∈ P2,ac (Rd ) with its probability density function (or Radon-Nikodym derivative) which we also R denote by ρ : Rd → [0, ∞), so ρ(x) ≥ 0 for all x ∈ Rd and Rd ρ(x) dx = 1. For ρ ∈ P(Rd ), let supp(ρ) denote the support of ρ, which is the smallest closed subset A ⊆ Rd with ρ(A) = 1. Let P2,ac,fs (Rd ) denote the subspace of P2,ac (Rd ) consisting of probability distributions with full support and strictly positive density function, so for all ρ ∈ P2,ac,fs (Rd ), supp(ρ) = Rd and the density function satisfies ρ(x) > 0 for all x ∈ Rd . For ρ, ν ∈ P(Rd ), we write ρ ≪ ν to denote ρ is absolutely continuous with respect to ν, which means if ν(A) = 0 for some A ⊆ Rd , then ρ(A) = 0. In particular, if ρ ≪ ν, then supp(ρ) ⊆ supp(ν). For ρ, ν ∈ P(Rd ), let Π(ρ, ν) denote the set of couplings of ρ and ν, i.e., joint probability distributions on R2d with marginal distributions ρ and ν. We denote by γ := N (0, Id ) the standard Gaussian distribution on Rd . 2.1.1
Wasserstein distance and KL divergence
We recall the following definitions of statistical distances and divergences. For ρ, ν ∈ P2 (Rd ), the Wasserstein–2 distance between ρ and ν is W2 (ρ, ν) :=
inf ω∈Π(ρ,ν)
h i 21 2 , E(x,y)∼ω ∥x − y∥
where the infimum is over all couplings of ρ and ν. For ρ, ν ∈ P2,ac,fs (Rd ) with ρ ≪ ν, the Kullback-Leibler (KL) divergence of ρ with respect to ν is Z h ρi ρ(x) KL(ρ ∥ ν) := Eρ log = ρ(x) log dx . ν ν(x) Rd For ρ, ν ∈ P2,ac,fs (Rd ) with ρ ≪ ν and where ρ and ν have differentiable density functions, the relative Fisher information (FI) of ρ with respect to ν is ρ 2 FI(ρ ∥ ν) = Eρ ∇ log . ν We recall the Wasserstein–2 distance, KL divergence, and relative Fisher information between any two distributions are non-negative, and they are equal to 0 if and only if the two distributions are the same. They are related via functional inequalities such as Talagrand or log-Sobolev inequality, see Section 2.1.3. 2.1.2
Convexity and Smoothness
Let f : Rd → R be a differentiable function. Recall the Bregman divergence Df : Rd × Rd → R is Df (x, y) = f (x) − f (y) − ⟨∇f (y), x − y⟩ for all x, y ∈ Rd . We say that f : Rd → R is convex if Df (x, y) ≥ 0 for all x, y ∈ Rd . We say f is α-strongly convex for some α > 0 if the function x 7→ f (x) − α2 ∥x∥2 is convex, or equivalently, Df (x, y) ≥ α2 ∥x − y∥2 for all x, y ∈ Rd . We say f is M -semi-convex for some M ≥ 0 if the function M 2 2 d x 7→ f (x) + M 2 ∥x∥ is convex, or equivalently, Df (x, y) ≥ − 2 ∥x − y∥ for all x, y ∈ R . We note α-strong convexity implies convexity (the case α = 0), and convexity implies M -semi-convexity (the 8
case M = 0). We remark that the class of semi-convex functions is quite broad, and refer readers to [DP19] for a thorough review. We say f is L-smooth for some L ∈ (0, ∞) if |Df (x, y)| ≤ L2 ∥x − y∥2 for all x, y ∈ Rd . We say that f is smooth if f is L-smooth for some L ∈ (0, ∞). When f is twice-differentiable, the definitions above have equivalent characterization in terms of the Hessian matrix ∇2 f , which we summarize in Table 1. Property Convexity α-strong convexity M -semi-convexity L-smoothness
Characterization in terms of Df Df (x, y) ≥ 0 Df (x, y) ≥ α2 ∥x − y∥2 2 Df (x, y) ≥ − M 2 ∥x − y∥ L 2 |Df (x, y)| ≤ 2 ∥x − y∥
Characterization in terms of ∇2 f ∇2 f (x) ⪰ 0 ∇2 f (x) ⪰ α · Id ∇2 f (x) ⪰ −M · Id L · Id ⪰ ∇2 f (x) ⪰ −L · Id
Table 1: Characterizations of convexity, semi-convexity, strong convexity, and smoothness. For a probability distribution ν ∝ e−f ∈ P2,ac,fs (Rd ), we say that ν is log-concave if f : Rd → R is a convex function. We say ν ∝ e−f is α-strongly log-concave (α-SLC) if f is an α-strongly convex function. We say ν ∝ e−f is M -semi-log-concave if f is an M -semi-convex function. We note that α-strong log-concavity implies log-concavity (the case α = 0), and log-concavity implies M -semilog-concavity (the case M = 0). We say ν ∝ e−f is L-log-smooth if f is an L-smooth function. We say ν ∝ e−f is log-smooth if f is a smooth function, i.e., f is L-smooth for some L ∈ (0, ∞). Throughout the paper, we assume that the target distribution ν X ∝ exp(−f ) is log-smooth and f is twice continuously differentiable; which from the early discussion in Section 2.1.2 means that there exists an ∞ > L > 0 such that L · Id ⪰ ∇2 f (x) ⪰ −L · Id for all x ∈ Rd . We note the log-smoothness of f is only used for regularity reasons to ensure the validity of some steps in the proofs (see Section 3.3 for more details); in particular, the quantity L does not appear in the main results presented in Section 3. 2.1.3
Functional inequalities
We state the key functional inequalities for distributions that we use in this work. We say ν ∈ P2,ac,fs (Rd ) satisfies Talagrand inequality with constant α > 0 if for any ρ ∈ P2,ac,fs (Rd ), we have α 2 W (ρ, ν) ≤ KL(ρ ∥ ν) . 2 2 We say ν ∈ P2,ac,fs (Rd ) satisfies an log-Sobolev inequality with constant α > 0 (α-LSI) if for any ρ ∈ P2,ac,fs (Rd ), we have 1 FI(ρ ∥ ν) . KL(ρ ∥ ν) ≤ 2α Consider ρ 7→ F(ρ) = KL(ρ ∥ ν). The Talagrand and the log-Sobolev inequalities stated above can be viewed as quadratic growth and gradient domination conditions on F respectively, in the space of probability measures equipped with the 2-Wasserstein metric; see [OV00, Wib18, VW23] for an extended discussion. In particular, we recall the following relationships, which we recall from [OV00]: • If ν is α-SLC, then ν satisfies α-LSI. • If ν satisfies α-LSI, then ν satisfies α-Talagrand inequality. 9
• If ν satisfies α-Talagrand and o is M -semi-log-concave for α ≥ M ≥ 0, then ν satisfies n inequality α M 2 β-LSI with β = max 4 1 − α , −M . Our first main result in Theorem 1 holds when the target distribution ν X satisfies α-Talagrand inequality and is M -semi-log-concave for some α > M ≥ 0, so by the preceding discussion, in fact ν X also satisfies an log-Sobolev inequality. This assumption (semi-log-concavity and Talagrand/logSobolev inequality) is consistent with the prior results in the accelerated convergence rates of the underdamped Langevin dynamics in [Lu26, LL26] and RHMC in [MW26].
3
Main results
We present the convergence guarantees for RHMC (Algorithm 1) in the settings [Set-Tri], [Set-Exp], f as defined in Section 1.1, as well as the corollaries on their iteration complexities. and [Set-Tri],
3.1
Convergence result under Talagrand inequality and semi log-concavity
We analyze Algorithm 1 in settings [Set-Tri] and [Set-Exp], assuming the target distribution ν X satisfies α-Talagrand inequality (Section 2.1.3) and is M -semi-log-concave (Section 2.1.2) for some 0 ≤ M < α < ∞. In particular, this includes the log-concave setting when M = 0. We provide the proof of Theorem 1 in Section 4.2.1. Theorem 1. Suppose ν X ∈ P2,ac,fs is a log-smooth distribution that satisfies α-Talagrand inequald ity and is M -semi-log-concave for 0 ≤ M < α < ∞. Let the initial distribution ρX 0 ∈ P2,ac,fs (R ) X X be such that KL(ρ0 ∥ ν ) < ∞. For K ∈ N and T > 0, define the sequence {Dk }k∈[K] in a settingspecific manner as: (a) for [Set-Tri], set Dk = TriT , and (b) for [Set-Exp], set Dk = Exp1/T for all k ∈ {1, . . . , K}. Then, in both settings, the distribution ρX K of xK output by Algorithm 1 satisfies 1 2 X X K X X 1+ . KL ρK ∥ ν ≤ C · KL ρ0 ∥ ν , where C := 3 − M/α αT 2 From the above theorem, we can immediately infer that the coefficient C is strictly less than 1 2 . This implies that for a sufficiently large integration time, the sequence provided T > Tmin := √α−M of distributions {ρX k } of random variables {xk } produced within Algorithm 1 converge to the target distribution ν X exponentially quickly. Using this result, we obtain the following corollary on the P continuous-time complexity of RHMC, which is the total integration time K k=1 τk of Algorithm 1 X ∥ ν X ) ≤ ε. We provide to output a sample XK ∼ ρX that has ε error in KL divergence, i.e., KL(ρ K K the proof of Corollary 1 in Section 4.2.2. Corollary 1. Consider the same assumptions on ν X and ρX 0 as in Theorem 1. For any ε > 0, X X X to output a sample XK ∼ ρK with guarantee KL(ρK ∥ ν ) ≤ ε, it suffices to run Algorithm 1 with 2 either [Set-Tri] or [Set-Exp] with T = √α−M , for the number of iterations
X ∥ νX ) KL(ρ 0 . log K= α−M ε log 1 + 1
5α−M
K X
"K # X 2K 2K The total integration time is τk ≤ √ in [Set-Tri], and E τk = √ in [Set-Exp]. α − M α − M k=1 k=1 10
In [Set-Exp] above, the expectation is over the randomness of the integration times τk . Although the bounds on the total integration time are the same in both [Set-Tri] and [Set-Exp], we note the difference between them. In [Set-Tri], when the integration time τk is drawn from the triangular 2K . In [Set-Exp], when τk is drawn distribution, the total integration time is always bounded by √α−M 2K , from the exponential distribution, only the expected total integration time is bounded by √α−M but the total integration time itself is a random variable which can be arbitrarily large. We also remark on the implication of the result above for log-concave sampling, i.e., the case M = 0. When the target distribution ν X is log-concave and satisfies α-Talagrand inequality, Corollary 1 states that the total integration time of Algorithm 1 scales as Õ(α−1/2 log(ε−1 )). This is the accelerated rate in continuous time, compared to the Õ(α−1 log(ε−1 )) complexity of the continuoustime Langevin dynamics in the same setting, and also matches what we expect from the theory of acceleration in convex optimization. Thus, our result above shows that RHMC indeed achieves an accelerated mixing time guarantee for log-concave sampling under Talagrand inequality. This is also consistent with the results of [Lu26, LL26] for underdamped Langevin dynamics and of [MW26] for RHMC.
3.2
Convergence result under log-concavity
f assuming log-concavity and log-smoothness of the target We now analyze Algorithm 1 in [Set-Tri] distribution. We prove Theorem 2 in Section 4.3.1. d Theorem 2. Assume ν X ∝ e−f ∈ P2,ac,fs (Rd ) is log-concave and log-smooth. Let ρX 0 ∈ P2,ac,fs (R ) k/2 9 X f f for all with KL(ρX 0 ∥ ν ) < ∞. For K ∈ N and in [Set-Tri], set Dk = TriTk where Tk = 8 X f k ∈ {1, . . . , K}. Then the distribution ρK of xk output by Algorithm 1 in [Set-Tri] satisfies: K 1 8 K 2 X X 1 2 X X 8 X X X X KL ρK ∥ ν + W2 ρK , ν ≤ KL ρ0 ∥ ν + W2 ρ0 , ν . (4) 3 9 9 3
The result above shows that the KL divergence to the target distribution converges exponentially fast in the number of iterations K, with a choice of integration time which increases exponentially in each iteration. We can extract the continuous-time complexity of RHMC to obtain a sample with a desired error in KL divergence. We provide the proof of Corollary 2 in Section 4.3.2. Corollary 2. Consider the same assumptions on ν X and ρX 0 as in Theorem 2. For any ε > 0, X X X to output a sample XK ∼ ρK with guarantee KL(ρK ∥ ν ) ≤ ε, it suffices to run Algorithm 1 with f with Tk = 9 k/2 for k ∈ {1, . . . , K}, for the number of iterations [Set-Tri] 8 & ' X + 1 W2 ρX , ν X KL ρX 1 0 ∥ν 0 3 2 K= log . log(9/8) ε The total integration time satisfies K X k=1
s τk ≤ 19
X + 1 W2 ρX , ν X KL ρX 0 ∥ν 0 3 2 . ε
Corollary 2 states that in the log-concave case, RHMC reaches error ε in KL divergence in a total integration time which scales as O(ε−1/2 ). We note this improves on the O(ε−1 ) complexity of the overdamped Langevin dynamics in the same log-concave setting, and also matches the improved complexity we expect from accelerated convex optimization. Thus, our result above shows that RHMC also achieves an accelerated mixing time guarantee in the log-concave case. 11
3.3
Proof technique: Average KL divergence along Hamiltonian dynamics
The proofs of the main algorithmic results crucially rely on a bound on the average KL divergence of the X-marginal along the trajectory of the deterministic Hamiltonian flow when the target distribution ν X is semi-log-concave; we present this bound in Lemma 3 in Section 3.3.2. We first review general properties of the KL divergence along the Hamiltonian dynamics (HF) in Section 3.3.1. 3.3.1
Distributional properties of the Hamiltonian dynamics
Recall from Eq. (1) the Hamiltonian function H : R2d → R is defined as H(x, y) := f (x) + 12 ∥y∥2 , which induces the joint probability distribution ν XY ∝ e−H on R2d . Note that ν XY = ν X ⊗ γ where γ = N (0, Id ) is the standard Gaussian distribution. Since we assume f is L-smooth for some L ∈ (0, ∞), from any initial condition (X0 , Y0 ) ∈ R2d , by Picard-Lindelöf theorem from standard ODE theory, the Hamiltonian flow (HF) is well defined and admits a unique solution (Xt , Yt ) ∈ R2d for all t ∈ R. A key property of the Hamiltonian flow (HF) is that it conserves the Hamiltonian function: H(Xt , Yt ) = H(X0 , Y0 )
for all t ∈ R ,
see Lemma 5 in Section D.1.1. Another key property is that the Hamiltonian flow conserves volume (Lebesgue measure) on the phase space R2d ; see Lemma 6 in Section D.1.2. Suppose we run the deterministic Hamiltonian flow (HF) from a random variable (X0 , Y0 ) ∼ ρXY 0 drawn from some initial distribution ρXY . Then at each time t ∈ R, we obtain another random 0 XY ∝ e−H is a stationary distribution along the variable (Xt , Yt ) ∼ ρXY t . We observe that ν XY XY XY Hamiltonian flow: If ρ0 = ν , then ρt = ν XY for all t ∈ R. Furthermore, the Hamiltonian flow conserves the KL divergence to ν XY . We provide the proof of Lemma 1 in Section D.2.2. Lemma 1. Let ρXY ∈ P2,ac,fs (R2d ) with KL(ρXY ∥ ν XY ) < ∞. For t ∈ R, let (Xt , Yt ) ∼ ρXY be t 0 0 the solution to the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρXY . Then we have 0 KL(ρXY ∥ ν XY ) = KL(ρXY ∥ ν XY ) . t 0
(5)
From Lemma 1, we see that ν XY is a stationary distribution of the Hamiltonian flow (HF). However, ν XY is not a unique stationary distribution. In fact, for any function ϕ : R → R with R −ϕ(H(x,y)) dx dy < ∞, the probability distribution ν̃ XY ∝ exp(−ϕ(H)) is also stationary along R2d e Hamiltonian flow (HF), and the KL divergence to this probability distribution is also conserved: KL(ρXY ∥ ν̃ XY ) = KL(ρXY ∥ ν̃ XY ) for all t ∈ R. t 0 What distinguishes the joint distribution ν XY = ν X ⊗ γ ∝ e−H in HMC is when we initialize X d the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρXY with ρXY = ρX 0 0 0 ⊗ γ for some ρ0 ∈ P2,ac,fs (R ). In this case, since the Y -marginal is initially chosen correctly from γ = N (0, Id ), which is the same as the Y -marginal of ν XY , this identifies the scale of the Hamiltonian function and isolates ν XY as a distinguished stationary distribution. In particular, if at some time t ∈ R we drop the Y -marginal from (Xt , Yt ) and only return the X-marginal Xt ∼ ρX t , then we obtain the following descent X property in KL divergence to the target ν . We provide the proof of Lemma 2 in Section D.2.3. X d X X Lemma 2. Let ρXY = ρX 0 0 ⊗ γ for some ρ0 ∈ P2,ac,fs (R ) with KL(ρ0 ∥ ν ) < ∞. For t ∈ R, let XY X (Xt , Yt ) ∼ ρt be the solution to the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρXY 0 , and let Xt ∼ ρt denote the X-marginal. Then we have X X X KL(ρX t ∥ ν ) ≤ KL(ρ0 ∥ ν ) .
12
X In Lemma 9 in Section D.2.4, we show the map t 7→ KL(ρX t ∥ ν ) is continuous. In Lemma 11 in Section D.2.5, we show that the second moment remains finite along Hamiltonian flow.
3.3.2
Bounding the average KL divergence
We now present the key lemma bounding the average of the KL divergence of the X-marginal along the trajectory of the Hamiltonian flow when the target distribution ν X is semi-log-concave. Lemma 3. Assume ν X is log-smooth and M -semi-log-concave for some 0 ≤ M < ∞. Assume d X X XY = ρX ⊗ γ. For t ≥ 0, let (X , Y ) ∼ ρXY ρX t t t 0 ∈ P2,ac,fs (R ) satisfies KL(ρ0 ∥ ν ) < ∞, and let ρ0 0 X be the solution to the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρXY 0 , and let Xt ∼ ρt denote the X-marginal. Then for all 0 ≤ T < ∞, the following holds: Z T Z M T 1 2 X X X X X (T − t) KL(ρt ∥ ν ) dt − (T − t) W22 (ρX W (ρ , ν ) + 3 t , ν ) dt 2 2 T 2 0 0 1 X 2 X X ≤ W22 (ρX 0 , ν ) + T KL(ρ0 ∥ ν ) . (6) 2 We provide the proof of Lemma 3 in Section G.5. The proof proceeds via two steps, that we sketch here. First, we assume the initial distribution ρX 0 satisfies a warmness and smoothness regularity condition (see Assumption 1). Under this regularity condition, we show the following differential inequality along the Hamiltonian flow: 1 d2 2 X X M 2 X X X X X (7) W2 (ρt , ν ) ≤ 2KL(ρX W (ρ , ν ) , 0 ∥ ν ) − 3KL(ρt ∥ ν ) + 2 2 dt 2 2 t see Lemma 22 in Section F.4 for the precise statement and proof (where we replace the second time derivative of Wasserstein distance by the limit of the second-order finite difference). Integrating this inequality in time twice results in the claimed inequality (6), under the regularity assumption on ρX 0 . To remove the regularity assumption, we show that we can approximate any initial distribution ρX 0 by a sequence of regular distributions, for which the claimed inequality holds, and we can take limits to obtain the conclusion of Lemma 3; see Section G for details. We remark that the derivation of Eq. (7) and its use to prove Lemma 3 are motivated by our earlier development in Hamiltonian-based accelerated optimization algorithms in [WSF+ 26], which follows exactly the same proof structure; see Section B for more details and comparison. The approximation argument to remove the regularity assumption follows a similar strategy used in [Lu26].
4
Proofs of the main results
4.1
Exponential distribution as a mixture of triangular distribution
We show that the exponential distribution (with density function given by Eq. (3)) can be viewed as a mixture of the triangular distribution (with density function given by Eq. (2)), with mixture proportion given by the Gamma distribution. This fact allows us to analyze Algorithm 1 with [Set-Exp] using the results of Algorithm 1 with [Set-Tri], as can be seen in Section 4.2. Let Γ(θ, λ) denote the Gamma distribution supported on [0, ∞) with shape parameter θ > 0 and rate parameter λ > 0, with density function at s ∈ [0, ∞) given by: λθ sθ−1 −λs e , Γ(θ) 13
where the normalization constant is the Gamma function Γ(θ) =
R∞ 0
λθ sθ−1 e−λs ds =
R ∞ θ−1 −t e dt. 0 t
Lemma 4. Let T ∈ (0, ∞). Let (S, τ ) ∈ [0, ∞) × [0, ∞) be a joint random variable drawn from the following process: S ∼ Γ(3, 1/T ) τ | S ∼ TriS . Then marginally, τ ∼ Exp1/T . Proof. Let pS denote the density of S (from the Gamma distribution), pτ |S denote the conditional density of τ given S (from the triangular distribution), and pτ denote the marginal density of τ . Note for θ = 3, the Gamma function is Γ(3) = 2! = 2. Then we can calculate the marginal density of τ at t ∈ [0, ∞) as: Z ∞ Z ∞ 2 s −s/T 2(s − t) e · 1[0,s] (t) ds pτ (t) = pS (s) pτ |S (t | s) ds = 2T 3 s2 0 0 Z ∞ 1 = 3 e−s/T (s − t) ds T t Z ∞ 1 −t/T 2 = 3e ·T e−u u du T 0 1 = e−t/T , T where in the above we have used change of variable u = (s−t)/T , du = ds/T , and the last equality follows by integration by parts. Thus, we conclude τ has marginal distribution Exp1/T .
4.2 4.2.1
Proof of Theorem 1 and Corollary 1 Proof of Theorem 1
Proof of Theorem 1. We prove a contraction of the KL divergence in each iteration of the algorithm. d X X We will show inductively that for each k ∈ {1, . . . , K}, we have ρX k ∈ P2,ac,fs (R ) and KL(ρk ∥ν ) < ∞, so we can apply Lemma 3 in each iteration. (k) (k) In iteration k ∈ {1, . . . , K}, we start the Hamiltonian flow from X0 , Y0 ∼ ρX k−1 ⊗ γ. Let (k) (k) (k) (k) = Ψt X0 , Y0 denote the solution of the Hamiltonian flow at time t ≥ 0, and let Xt , Yt (k)
X X ρX t,k denote the law of Xt . Note that ρ0,k = ρk−1 . By Lemma 3, for all T ∈ [0, ∞) we have: Z T Z 1 2 X M T X X X X W2 ρT,k , ν + 3 (T − t) KL ρt,k ∥ ν dt − (T − t) W22 ρX dt t,k , ν 2 2 0 0 1 X X + T 2 KL ρX . ≤ W22 ρX 0,k , ν 0,k ∥ ν 2 Assuming ν X satisfies α-Talagrand inequality, we can bound the W22 -term on the left-hand side by: Z Z M T M T 2 X X X − (T − t) W2 ρt,k , ν dt ≥ − (T − t) KL ρX dt . t,k ∥ ν 2 0 α 0
Therefore, the following inequality holds for all T ∈ [0, ∞): Z T 1 2 X M 1 X X X X W2 ρT,k , ν + 3 − (T −t) KL ρX dt ≤ W22 ρX +T 2 KL ρX . (8) t,k ∥ ν 0,k , ν 0,k ∥ ν 2 α 2 0 We now analyze the settings [Set-Tri] and [Set-Exp] separately. 14
1. Setting [Set-Tri]: In [Set-Tri], we draw the integration time τk from the triangular distribution TriT supported on [0, T ]. Therefore, the output density ρX k is equal to Z T X 2(T − t) X ρt,k (x) dt . ρX = E ρ = τk ∼TriT k τk ,k T2 0 By the convexity of KL divergence and Jensen’s inequality, Z T Z T 2(T − t) 2(T − t) X X X X X X ≤ ρ (x) dt ν KL ρ ∥ ν dt . KL ρk ∥ ν = KL t,k t,k T2 T2 0 0
(9)
Plugging this into Eq. (8) and simplifying, we get: 1 2 X M T2 M T2 X X X X 3− KL ρk ∥ ν ≤ W2 ρT,k , ν + 3 − KL ρX k ∥ν α 2 2 α 2 Z T (9) 1 M X X (T − t) KL ρX dt ≤ W22 ρX + 3− t,k ∥ ν T,k , ν 2 α 0 (8) 1 X X + T 2 KL ρX ≤ W22 ρX 0,k ∥ ν 0,k , ν 2 1 X 2 ≤ , + T KL ρX 0,k ∥ ν α X where the last inequality follows from applying α-Talagrand inequality. Recalling ρX 0,k = ρk−1 and T2 dividing both sides above by 3 − M α 2 > 0, we conclude that 2 1 X X X KL ρk ∥ ν ≤ 1+ KL ρX . (10) k−1 ∥ ν 2 3 − M/α αT
Iterating the bound above shows the claim in Theorem 1 for the setting [Set-Tri]. X d X X We now check that ρk ∈ P2,ac,fs (R) with KL ρk ∥ ν < ∞. Assume inductively that ρX k−1 = d ) for all X d X X ∈ P (R < ∞. By Lemma 11, we know ρX ρ0,k ∈ P2,ac,fs (R ) with KL ρk−1 ∥ ν t,k i h2,ac,fs i h 2 d X X t ∈ [0, T ]. Hence, ρk = Eτk ∼TriT ρτk ,k ∈ P2,ac,fs (R ), since EρX [∥X∥ ] = Eτk ∼TriT EρX [∥X∥2 ] ≤ k
τk ,k
has full support and positive density. By Lemma 2, we know max0≤t≤T EρX [∥X∥2 ] < ∞, and ρX t,k k X ≤ KL ρX ∥ ν X < ∞. Hence, by Eq. (9), KL ρX ∥ ν X ≤ KL ρX ∥ ν X < ∞. ∥ ν KL ρX k 0,k 0,k t,k 2. Setting [Set-Exp]: In [Set-Exp], we draw the integration time τk from the exponential distribution Exp1/T with rate parameter 1/T . Recall from Lemma 4 that we can write the density of the exponential distribution as a mixture of the triangular and Gamma distributions: For each t ∈ [0, ∞), Z ∞ 2 1 −t/T s −s/T 2(s − t) e = e 1[0,s] (t) ds . T 2T 3 s2 0 Z s 2(s − t) X X For s > 0, define ρ̄s,k := ρt,k dt. Then the density of ρX k in this setting becomes 2 s 0 Z ∞ X 1 −t/T X ρk = Eτk ∼Exp1/T ρt,k = ρX e dt t,k T 0 Z ∞Z ∞ s2 −s/T 2(s − t) = ρX e 1[0,s] (t) ds dt t,k 2T 3 s2 0 0 15
Z ∞ Z ∞ = 0
Z ∞
0
2(s − t) ρX 1[0,s] (t) dt t,k 2 s
s2 −s/T e ds 2T 3
s2
e−s/T ds 3 2T 0 = ES∼Γ(3,1/T ) ρ̄X S,k =
ρ̄X s,k
where in the third line above we have exchanged the order of integration. Next, by the convexity of KL divergence and Jensen’s inequality, Z ∞ Z ∞ s2 −s/T s2 −s/T X X X X X X KL ρ̄ ∥ ν KL ρk ∥ ν = KL ρ̄s,k e ds ν ≤ e ds . (11) s,k 2T 3 2T 3 0 0 Note that ρ̄X s,k is the result of applying one iteration of the algorithm with the setting [Set-Tri], so the bound Eq. (10) holds for ρ̄X s,k (with T in Eq. (10) replaced by s). Then continuing from the above, we obtain: Z s2 −s/T (11) ∞ X X X ≤ KL ρ̄X e ds KL ρk ∥ ν s,k ∥ ν 2T 3 0 Z ∞ (10) s2 −s/T 1 2 X 1 + 2 KL ρX e ds ≤ k−1 ∥ ν 3 − M/α αs 2T 3 0 Z ∞ 1 2 1 1 −s/T X = KL ρX + s2 e ds k−1 ∥ ν 2 3 − M/α 2T α T 0 1 2 X X ∥ ν 1+ KL ρ = k−1 3 − M/α 2αT 2 2 1 X ≤ . 1+ KL ρX k−1 ∥ ν 2 3 − M/α αT In the last equality above, we have used the fact that the exponential distribution Exp1/T with density T1 e−s/T has second moment equal to 2T 2 . Iterating the bound above shows the claim in Theorem 1 for the setting [Set-Exp]. d X X We now check that ρX < ∞. Assume inductively that k ∈ P2,ac,fs (R ) with KL ρk ∥ ν d ) with KL ρX ∥ ν X X ∈ P (R < ∞. As in the setting [Set-Tri], we know = ρ ρX k−1 k−1 0,k 2,ac,fs X X X < ∞ for any s ∈ (0, ∞). By (11), we also have ρ̄X ≤ KL ρX 0,k ∥ ν s,k satisfies KL ρ̄s,k ∥ ν X ≤ KL ρX ∥ ν X < ∞. Next, recall from Lemma 11 the second moment of ρX grows KL ρX t k ∥ν 0,k at most quadratically hin t,iso it is integrable against any exponentially decreasing function. ThereX X fore, ρX k = Eτk ∼Exp1/T ρt,k has finite second moment, and ρk has full support and positive density, d so ρX k ∈ P2,ac,fs (R ).
Conclusion: In both settings gives the claimed con [Set-Tri] and [Set-Exp], iterating the bound 2 X . Plugging in the choice T ≥ T √ vergence rate on KL ρX ∥ ν = gives the contraction min K α−M constant as 2 5α − M α−M α − M −1 C≤ 1+ = = 1+ < 1, (12) 3 − M/α 4α 6α − 2M 5α − M where the last inequality holds since we assume M < α. 16
4.2.2
Proof of Corollary 1
2 Proof of Corollary 1. Take T = Tmin = √α−M , and plug this in to the formula of the contraction ratio C from Eq. (12) into Theorem 1 to obtain:
KL
X ρX K ∥ν
α − M −K X · KL ρX ≤ 1+ 0 ∥ν 5α − M X X 1 KL(ρ0 ∥ ν ) α−M X · log = exp − · log 1 + · KL ρX 0 ∥ν α−M ε 5α − M log 1 + 5α−M X KL(ρX X 0 ∥ν ) · KL ρX ≤ exp − log 0 ∥ν ε = ε.
In [Set-Tri], since τk ∼ TriTmin which is supported on [0, Tmin ], we have τk ≤ Tmin . Therefore, the total integration time satisfies K X
τk ≤ K · Tmin = √
k=1
2K . α−M
In [Set-Exp], since τk ∼ Exp1/Tmin , we have E[τk ] = Tmin . Therefore, the expected total integration time is "K # X 2K E τk = K · Tmin = √ . α−M k=1
4.3
Proof of Theorem 2 and Corollary 2
4.3.1
Proof of Theorem 2
Proof of Theorem 2. We prove the following more general result under log-concavity and logf if we set Dk = Tri f T for any smoothness of the target distribution: For Algorithm 1 under [Set-Tri], k X f Tk > 0 for all k ∈ {1, . . . , K} (with T0 = 1), then the output ρ of Algorithm 1 under [Set-Tri] K
satisfies KL ρX K ∥ν
X
)# "K ( 2 Y 1 2 X X T 1 8 k−1 2 X X X X + KL ρ0 ∥ ν + W2 ρ0 , ν ≤ max , . 2 W2 ρK , ν 9 Tk2 3 3TK k=1
(13) k/2 Setting Tk = 98 yields the claimed bound in Theorem 2. To show Eq. (13), we will prove a contraction in each iteration. We show inductively below d X X that ρX k ∈ P2,ac,fs (R ) and KL(ρk ∥ ν ) < ∞ for all k ∈ {1, . . . , K}, so we can apply Lemma 3 in each iteration. (k) (k) In iteration k ∈ {1, . . . , K}, we start the Hamiltonian flow from X0 , Y0 ∼ ρX k−1 ⊗ γ. Let (k) (k) (k) (k) Xt , Yt = Ψt X0 , Y0 denote the solution of the Hamiltonian flow at time t ≥ 0, and let
17
(k)
X X X ρX t,k denote the law of Xt . Note that ρ0,k = ρk−1 . By Lemma 3 (with M = 0 since we assume ν is log-concave), for all T ∈ [0, ∞) we have:
1 2 X W2 ρT,k , ν X + 3 2
Z T
1 X X X (T − t) KL ρX dt ≤ W22 ρX + T 2 KL ρX . t,k ∥ ν 0,k , ν 0,k ∥ ν 2 0
(14)
f T = 1 TriT + 1 δT , the output distribution ρX at Since we draw the integration time from Tri k k k 2 2 k iteration k is 1 X 1 X ρX k = ρ̄Tk ,k + ρTk ,k 2 2 where Z Tk X 2(Tk − t) X X ρt,k dt . ρ̄Tk ,k := Eτk ∼TriTk ρτk ,k = Tk2 0 We now derive several inequalities that we will combine to obtain the result. First inequality: By the convexity of KL divergence and Jensen’s inequality, and by applying Eq. (14) at T = Tk , we obtain Z Tk 1 2 X 3Tk2 3Tk2 1 2(Tk − t) X X X X X X = + W22 ρX KL ρ̄Tk ,k ∥ ν + W2 ρTk ,k , ν KL ρt,k dt ν Tk ,k , ν 2 2 2 2 2 Tk 0 Z Tk 1 X X ≤3 (Tk − t) KL ρX dt + W22 ρX t,k ∥ ν Tk ,k , ν 2 0 (14) 1 2 X X X + W2 ρ0,k , ν ≤ Tk2 KL ρX . 0,k ∥ ν 2 X Multiplying both sides by 2/(3Tk2 ) and recalling ρX 0,k = ρk−1 yield:
2 1 1 X X X X + KL ρ̄X + ≤ KL ρX . W22 ρX W22 ρX Tk ,k ∥ ν k−1 ∥ ν Tk ,k , ν k−1 , ν 2 2 3 3Tk 3Tk
(15)
Second inequality:
By the convexity of squared Wasserstein distance, we can bound: 1 X 1 2 X W2 ρ̄Tk ,k , ν X = W22 Eτk ∼TriTk ρX τk ,k , ν 2 2 1 2 X X ≤ Eτk ∼TriTk W ρ ,ν 2 2 τk ,k Z τk 1 2 X X X X ≤ Eτk ∼TriTk W ρ , ν +3 (τk − t) KL ρt,k ∥ ν dt 2 2 τk ,k 0 (14) 1 2 X X 2 X X W ρ , ν + τk KL ρ0,k ∥ ν ≤ Eτk ∼TriTk 2 2 0,k Tk2 1 X X , ν + KL ρX . = W22 ρX 0,k 0,k ∥ ν 2 6
In the above, we have applied the bound from (14) at each T = τk ∈ [0, T ], and we have used the X fact that Eτ ∼TriT [τ 2 ] = T 2 /6. Multiplying both sides by 2/(3Tk2 ) and recalling ρX 0,k = ρk−1 yield: 1 1 1 X X X W22 ρ̄X ≤ W22 ρX + KL ρX . Tk ,k , ν k−1 , ν k−1 ∥ ν 2 2 9 3Tk 3Tk 18
(16)
Third inequality: we have:
By the descent property in KL divergence along one step of HMC (Lemma 2),
X X X KL ρX ≤ KL ρX = KL ρX . Tk ,k ∥ ν 0,k ∥ ν k−1 ∥ ν Combining:
(17)
Summing Eq. (15), Eq. (16), and Eq. (17) gives:
1 X X X X W22 ρ̄X + W22 ρX KL ρ̄X + KL ρX + Tk ,k , ν Tk ,k , ν Tk ,k ∥ ν Tk ,k ∥ ν 2 3Tk 16 2 X ≤ KL(ρX W 2 (ρX , ν X ) , k−1 ∥ ν ) + 9 3Tk2 2 k−1
(18)
2 1 where in the above, 16 9 = 1 + 3 + 9. 1 X 1 X Since ρX k = 2 ρ̄Tk ,k + 2 ρTk ,k , by the convexity of KL divergence and the convexity of the squared Wasserstein distance, we can bound:
1 X 2 X X KL ρX + W ρ , ν 2 k k ∥ν 3Tk2 1 X 1 1 X 1 X X X 2 1 X + ρ ν + + ρ ,ν = KL ρ̄ W ρ̄ 2 Tk ,k 2 Tk ,k 3Tk2 2 2 Tk ,k 2 Tk ,k 1 1 X X X X X X 2 X 2 X ≤ KL ρ̄Tk ,k ∥ ν + KL ρTk ,k ∥ ν + W2 ρ̄Tk ,k , ν + W2 ρTk ,k , ν 2 3Tk2 (18) 8 1 X X ≤ KL ρX + W22 ρX . k−1 , ν k−1 ∥ ν 2 9 3Tk 2 2 8 8 Tk−1 1 1 8 Tk−1 Finally, bounding 9 ≤ max 9 , T 2 and 3T 2 ≤ 3T 2 · max 9 , T 2 yields: k
KL
X ρX k ∥ν
1 X W22 ρX ≤ max + k ,ν 2 3Tk
k
(
k−1
2 8 Tk−1 , 2 9 Tk
k
) KL
X ρX k−1 ∥ ν
1 X + W22 ρX k−1 , ν 2 3Tk−1
! .
Telescoping this bound over k ∈ {1, . . . , K} and recalling we define T0 = 1 prove Eq. (13). d ) with KL ρX ∥ ν X < ∞. Assume inductively that ρX We now check that ρX ∈ P (R 2,ac,fs k k−1 = k X d X X X d ρ0,k ∈ P2,ac,fs (R ) with KL ρk−1 ∥ ν < ∞. By Lemma 11 we know that ρt,k ∈ P2,ac,fs (R ) for all h i 1 1 X X d t ∈ [0, Tk ], and hence ρk = 2 Eτk ∼TriTk ρX τk ,k + 2 ρTk ,k ∈ P2,ac,fs (R ), since h i 1 1 EρX [∥X∥2 ] = Eτk ∼TriTk EρX [∥X∥2 ] + EρX [∥X∥2 ] ≤ max EρX [∥X∥2 ] < ∞ , k τk ,k t,k 0≤t≤Tk 2 2 Tk ,k X X < ∞ for all and ρX has full support and positive density. By Lemma 2, we know KL ρ ν k t,k t > 0, and hence by the convexity of KL divergence, we also have 1 1 X X X X X X KL ρX ≤ Eτk ∼TriTk KL(ρX τk ,k ∥ ν ) + KL(ρTk ,k ∥ ν ) ≤ max KL(ρt,k ∥ ν ) < ∞ . k ∥ν 0≤t≤Tk 2 2 This completes the proof.
19
4.3.2
Proof of Corollary 2
Proof of Corollary 2. Since Tk =
9 k/2 , by Theorem 2, after K iterations, 8
K 1 2 X X 8 X X KL(ρ0 ∥ ν ) + W2 (ρ0 , ν ) . 9 3 & ' X + 1 W2 ρX , ν X KL ρX 1 2 0 0 ∥ν 3 X Plugging in the choice K = gives KL(ρX · log K ∥ ν ) ≤ ε. ε log( 98 ) It remains to calculate the total integration time. For this choice of K, we first calculate that !! K/2 X + 1 W2 ρX , ν X KL ρX 9 1 1 9 0 0 ∥ν 3 2 · ≤ exp log +1 · log 8 2 8 ε log( 98 ) r s 1 2 X X X KL(ρX 9 0 ∥ ν ) + 3 W2 ρ0 , ν = · . 8 ε X KL(ρX K ∥ν )≤
f T = 1 TriT + 1 δT , and hence τk ≤ Tk = In iteration k, the integration time is drawn from Tri k k 2 2 k Therefore, ! K/2 K K k/2 X X 9 9 3 √ · τk ≤ = −1 8 8 3 − 8 k=1 k=1 r s 1 2 X X X X KL ρ0 ∥ ν + 3 W2 ρ0 , ν 3 9 √ · · − 1 ≤ 8 ε 3− 8 r s X + 1 W2 ρX , ν X √ KL ρX 9 0 ∥ν 2 0 3 ≤ 3(3 + 8) · · 8 ε s X + 1 W2 ρX , ν X KL ρX 0 0 ∥ν 3 2 ≤ 19 . ε
5
9 k/2 . 8
Discussion
In this work, we prove accelerated mixing time guarantees for the idealized Randomized Hamiltonian Monte Carlo (RHMC) algorithm in two settings: when the target distribution is semi-logconcave and satisfies Talagrand inequality, and when the target distribution is log-concave. The resulting continuous-time complexity guarantees of RHMC improve on the guarantees for the overdamped Langevin dynamics under the same settings, and match the rates that we expect from the theory of accelerated convex optimization. The key technical ingredient underlying our result is Lemma 3, which shows that the average KL divergence of the position marginal along Hamiltonian flow is decreasing with an explicit contraction factor. To apply this result, it is important that the integration time in RHMC be randomized, with a sufficiently long expected value. Our analysis is motivated by our prior work on Hamiltonian dynamics-based optimization [WSF+ 26].
20
Our results complement the recent results on the accelerated mixing time guarantees for underdamped Langevin dynamics [Lu26, LL26] and RHMC [MW26], which proceed via different analysis techniques. Our results in this paper are for the idealized RHMC, where we assume we can simulate the Hamiltonian flow exactly. An important future direction is to study how to extend these guarantees to discrete-time implementations of RHMC, where Hamiltonian flow is implemented using a numerical integrator, and whether we can obtain discrete-time iteration complexity guarantees that match what we can obtain from accelerated convex optimization. It would also be interesting to investigate whether similar guarantees can be established for the No-U-Turn Sampler (NUTS) [HG14], an adaptive variant of HMC that is widely used in practice, and whose theory is still being developed [BRO24, Obe25, BRCKL26, GKG+ 26].
21
A
Additional related work
Guarantees for Hamiltonian dynamics-based sampling algorithms such as RHMC (Algorithm 1) are an active area of research. As described in Section 1.1, when the integration time in RHMC is taken to be deterministic, the algorithm is usually referred to as HMC. Guarantees for HMC correspond to Rows 7-9 in Table 2 and are all unaccelerated; as shown in [CV22, Theorem 1.4], this is unavoidable whenever the integration time is deterministic. Table 2 also includes a summary of prior works studying the overdamped Langevin dynamics (LD) and the underdamped Langevin dynamics (ULD), and in particular, includes all of the works discussed in Section 1. These results presented in Section 1 are part of a broader question in sampling which is to obtain a diffusive-to-ballistic speedup: when the target distribution is α-strongly log-concave, the goal is to improve the total simulation time required to obtain ε-accurate samples from the diffusive scale or unaccelerated rate of O(α−1 log(ε−1 )), achieved by the Langevin dynamics (Table 2, Row 1), to the ballistic scale or accelerated rate of O(α−1/2 log(ε−1 )). One line of work is based on hypocoercivity and space-time Poincaré inequalities. The space-time Poincaré approach was developed in [AAMN24]. Building on this framework, [CLW23] obtain accelerated convergence rates for underdamped Langevin dynamics in χ2 divergence (Table 2, Row 3). This area of research has been particularly active recently, with multiple works adapting these techniques to obtain guarantees in relative entropy. For instance, [Lu26] generalize this line of results and obtain accelerated guarantees for underdamped Langevin dynamics in KL divergence (Table 2, Row 5), while [LL26] further extend it to Rényi divergence (Table 2, Row 6). We remark that moving from χ2 divergence to KL divergence is not purely of theoretical interest, but is also often desirable in high-dimensional d settings due to the milder dependence on initialization: for typical initial distribution ρX 0 ∈ P(R ) X X d X and reference distribution ν ∈ P(R ), the initial KL divergence KL(ρ0 ∥ ν ) scales linearly in the X dimension d, whereas χ2 (ρX 0 ∥ ν ) can scale exponentially in d. A further perspective on acceleration is provided by the non-reversible lifting framework of [EL26]. They formalize the underdamped Langevin dynamics and RHMC as second-order non-reversible lifts of overdamped Langevin dynamics, and established accelerated convergence of RHMC in χ2 divergence when the target distribution satisfies a Poincaré inequality (Table 2, Row 11). There are other attempts at translating the accelerated gradient flow dynamics from optimization to the space of probability distributions for sampling [WL22, CLTW25], resulting in mean-field dynamics which have accelerated convergence guarantees, but may be more challenging to implement algorithmically than ULD. When the target distribution is Gaussian, improved guarantees for HMC can be obtained by carefully leveraging the structural properties of Gaussian distributions. [WW23] construct a deterministic, time-varying integration-time schedule from the roots of a Chebyshev polynomial constructed from the covariance matrix of the Gaussian distribution; they show that using this integration time schedule, HMC reaches ε error in Wasserstein-2 distance using total integration time O α−1/2 log(ε−1 ) . [Jia23] show the same continuous-time complexity can be achieved using RHMC for Gaussian distributions either with exponentially distributed integration times or with partial velocity refreshment. [AGS24] analyze a Metropolis-adjusted HMC implementation with long ran e L1/2 α−1/2 d1/4 log(ε−1 ) in domized integration times and obtain a gradient-query complexity of O total variation distance.
22
1
Reference
Algorithm/Dynamics
Divergence
Total Time
Accelerated
[OV00]
LD
KL
O(α−1 log(ε−1 ))
No No
2
[ACZ26]
ULD
KL, Rq
O(L1/2 α−1 log(ε−1 ))
3
[CLW23]
ULD
χ2
O(α−1/2 log(ε−1 ))
Yes
4
[FLL26]
ULD
χ2
O(α−1/2 log(ε−1 ))
Yes Yes
5
[Lu26]
ULD
KL
O(α−1/2 log(ε−1 ))
6
[LL26]
ULD
Rq
O(α−1/2 log(ε−1 ))
Yes
7
[MS21]
HMC; T = O(α1/2 L−1 )
W2
O(L α−3/2 log(ε−1 ))
No
8
[CV22]
HMC; T = O(L−1/2 )
W2
O(L1/2 α−1 log(ε−1 ))
No
9
[Mon24]
HMC; T = O(L−1/2 )
KL
O(L1/2 α−1 log(ε−1 ))
No
O(α−1/2 log(ε−1 ))
Yes
10
[LW22]
RHMC; [Set-Exp]
χ2
11
[EL26]
RHMC; [Set-Exp]
χ2
O(α−1/2 log(ε−1 ))
Yes
12
[MW26]
RHMC; [Set-Exp]
KL
O(α−1/2 log(ε−1 ))
Yes Yes
13
Theorem 1
RHMC; [Set-Tri], [Set-Exp]
KL
O(α−1/2 log(ε−1 ))
14
[OV01]
LD
KL
O(ε−1 )
No
15
[ACZ26]
ULD
KL, Rq
O(L1/2 ε−1 )
No
16
Theorem 2
f RHMC; [Set-Tri]
KL
O(ε−1/2 )
Yes
Table 2: Summary of total time complexity of prior works to output a sample within error ε > 0 (Column 4) in various distance or divergences (Column 3). In Column 3, W2 and KL denote the Wasserstein–2 distance and KL divergence respectively (Section 2.1.1); and χ2 and Rq denote the chi-squared divergence and Rényi–q divergence respectively. Rows 1-13 correspond to works where α-strong log-concavity (Section 2.1.2) or isoperimetry such as α-LSI (Section 2.1.3) is assumed, and Rows 14-16 correspond to guarantees under log-concavity. The dependence on the L-smoothness is presented as and when it arises. The algorithms LD and ULD correspond to the overdamped and underdamped Langevin dynamics respectively, as discussed in Section 1; algorithms RHMC and HMC are defined in Section 1.1. As discussed in Section 1, [CV22, Theorem 1.4] implies that obtaining accelerated or ballistic guarantees for HMC for general α-strongly convex and L-smooth target distributions is not feasible due to the deterministic integration time. To address this concern, [CLS07, BRSS17] consider randomized integration times, i.e., RHMC (Algorithm 1) with non-degenerate integration time distributions. [LW22] provides rigorous mixing time guarantees for RHMC and they obtain accelerated rates of convergence for RHMC in χ2 divergence under a space-time Poincaré inequality (Table 2, Row 10). Very recently, [MW26] establish diffusive-to-ballistic acceleration in KL divergence for continuous-time RHMC with exponentially distributed integration times (Table 2, Row 12). Their proof adapts the approach of [Lu26] (Table 2, Row 5) by considering the RHMC semigroup, introducing a similar Lyapunov functional, and studying the evolution of this functional under the semigroup of RHMC.
23
A related approach to avoiding a poorly chosen deterministic integration time is the No-U-Turn Sampler (NUTS) [HG14], which locally adapts the length of a leapfrog trajectory using a geometric U-turn criterion and then selects the next state from the resulting orbit. For the standard Gaussian target, [BRO24] establish the first quantitative mixing-time guarantee for NUTS, and there are follow-up works including [Obe25, BRCKL26, GKG+ 26]. Another line of work, including [TPL19, Wan24, FW25, WSF+ 26], studies Hamiltonian dynamics as an algorithmic primitive for optimization. The present work is specifically inspired by our recent contribution [WSF+ 26] to this line of research. In this optimization setting, Hamiltonian flow conserves energy in the phase space, and therefore its trajectories can oscillate and may not converge. The key observation in [WSF+ 26] is that a suitable time average of the trajectory of Hamiltonian flow satisfies an accelerated convergence guarantee. The same principle drives the current work. Hamiltonian flow preserves the KL divergence to the joint stationary distribution in the phase space, but the time-averaged position marginal along the Hamiltonian trajectory admits a descent guarantee in KL divergence. We provide a more detailed discussion in Section B. The preceding discussion has all been for continuous-time sampling methods. There is a rich body of work focusing on discrete-time algorithmic implementations of these ideal continuous-time methods. We briefly mention that discrete-time implementations of HMC largely focus on short integration times of T = O(L−1/2 ), as this corresponds to the regime when discretization schemes such as velocity Verlet are stable [BRSS18]. Discretizing the Hamiltonian dynamics introduces a bias in the Markov chain, which, when not adjusted using a Metropolis-Hastings filter, leads to unadjusted Hamiltonian Monte Carlo; some works studying this are [BREZ20, BRE23, BRS23, CDMS23, GBMM25, BRMW26]. Adding the Metropolis-Hastings filter restores the target distribution as the stationary distribution and some works studying adjusted Hamiltonian Monte Carlo include [BREZ20, CDWY20, CGJ26]. It is interesting to study if the accelerated mixing time guarantees we present here for the continuous-time RHMC can be extended to discrete-time algorithms. We conclude by mentioning that a partial progress toward discrete-time acceleration was made in [ACZ26] for an algorithm based on ULD, albeit still for a low-accuracy guarantee.
B
Connections to Hamiltonian Flow for optimization
In this section, we present a more detailed discussion about how the main results of this work relate to the recent results of [WSF+ 26] who focus on developing a Hamiltonian flow based algorithm for accelerated convex optimization. We preface the discussion by highlighting a classical connection between optimization and sampling that we alluded to in Section 1. In optimization, given a function f : Rd → R, the Euclidean gradient flow is a continuous-time dynamics that converges to a stationary point of f . Consequently, when f is convex, running gradient flow leads to a minimizer of f . Analogously, one can treat the algorithmic task of sampling from ν X as finding a distribution that minimizes a discrepancy to ν X . When the space of probability measures is endowed with the Wasserstein metric, and the discrepancy is chosen to be the KL divergence, the gradient flow of ρ 7→ KL(ρ∥ν X ) precisely coincides with the overdamped Langevin dynamics (LD), and structural assumptions on ν X such as log-concavity and isoperimetry have natural geometric interpretations that allow efficient optimization of KL divergence. This connection was originally discovered in the seminal work of [JKO98], and has recently been popular for many applications, see e.g. [Wib18] for further discussion. In contrast, the present work is using Hamiltonian dynamics, which has a conservation property, rather than 24
gradient flow which is dissipative. We begin the discussion with a key conceptual similarity between using HF for sampling and HF for optimization. Recall Lemma 2 states that the law of the position Xt obtained by running X X X X HF from (X0 , Y0 ) ∼ ρXY = ρX 0 0 ⊗ γ satisfies KL(ρt ∥ ν ) ≤ KL(ρ0 ∥ ν ). Intriguingly, a similar implication of HF can be shown for optimization with a different initialization for Y0 ; more precisely, when initializing from Y0 = 0. While Lemma 2 results from the volume preserving structure of the Hamiltonian flow, the property for optimization results from the conservation of the Hamiltonian H(x, y) by HF, which implies that for any time t: 1 1 f (Xt ) − f (x⋆ ) + ∥Yt ∥2 = f (X0 ) − f (x⋆ ) + ∥Y0 ∥2 . 2 2 Using the non-negativity of the squared norms and the initialization Y0 = 0 leads to the property: f (Xt ) − f (x⋆ ) ≤ f (X0 ) − f (x⋆ ). Both of these implications can be viewed as establishing an unconditional descent; for sampling, this is descent in the KL divergence of the X-marginal, while for optimization, this is descent in the optimality gap of the X-iterate. The implication for optimization was originally discovered in [TPL19], and leads to an iterative algorithm for optimization involving: (1) simulating HF for a certain amount of integration time, and (2) resetting the velocity to 0. Note the sole difference to HMC is in how the velocity is reset: in HMC we reset the velocity to be a fresh Gaussian. Later work by [Wan24] develop this insight of [TPL19] and show that for specific choices of integration times, the scheme for optimization described above can yield accelerated convergence for minimizing convex quadratic functions. The result of [Wan24] was based on a previous result by [WW23], where similar integration times led to accelerated convergence of HMC for sampling from multivariate Gaussian distributions. However, it remained unclear if algorithms based on HF can lead to accelerated algorithms for general differentiable convex objectives. More recent work by [FW25] demonstrated that by choosing randomized integration times from an exponential distribution — borrowing inspiration from RHMC [BRSS17] — one can minimize differentiable convex objectives in an accelerated manner, albeit in expectation over the randomness of the integration times. This leads to the work by [WSF+ 26] who develop new insights about HF, and propose an accelerated algorithm for differentiable convex objectives by averaging trajectories as shown in the algorithm stated below. Algorithm 2 Hamiltonian flow for optimization with averaging (HFA-opt) [WSF+ 26] Require: Initialization: X0 ∈ Rd ; number of iterations: K ∈ N; parameter: λ ≥ 0; integration time sequence: {T1 , . . . , TK }. for k = 1, . . . , K do (k) (k) Set X0 = xk−1 , and set Y0 = 0. (k) (k) (k) (k) Solve HF from (X0 , Y0 ) = (Xk−1 , 0) for time Tk to obtain trajectory ((Xt , Yt ))t∈[0,Tk ] . Z Tk 2 (k) Compute X avg (Xk−1 ; Tk ) = 2 (Tk − t)Xt dt . Tk 0 λ 1 (k) X avg (Xk−1 ; Tk ) + XTk . Set next iterate Xk = λ+1 λ+1 end for return xK . We note that while the outline of the above algorithm is similar to RHMC (Algorithm 1), a key 25
difference is that RHMC (Algorithm 1) chooses the integration time randomly according to TriT or f T , whereas in HFA-opt (Algorithm 2), the trajectory is averaged according to TriT (when λ = 0) Tri f T (when λ = 1). However, the random selection of integration times in RHMC indeed leads or Tri to averaging the distributions in the space of measures. For example, in [Set-Tri] we have: ρX k = Eτk ∼TriT
X ρτk ,k =
Z T 0
2(T − t) X ρt,k (x) dt , T2
see Section 4.2.1 and Section 4.3.1 for details. As a result, Algorithm 1 can be interpreted as the sampling analogue of Algorithm 2. The connection to [WSF+ 26] does not end at the algorithmic level, but also appears in the proof techniques used to obtain guarantees for Algorithm 1 from this work and Algorithm 2 from [WSF+ 26]. As discussed in Section 3.3, Lemma 3 is the key lemma that paves the road to deriving Theorems 1 and 2. The statement of Lemma 3 results from integrating the differential inequality stated in Eq. (7). A similar differential inequality is derived in [WSF+ 26, Lemma 2], which states that when f is a convex function, for any X0 ∈ Rd , the trajectory (Xτ , Yτ )τ ∈[0,t] obtained by solving HF with initial conditions (X0 , 0) satisfies 1 d2 ∥Xt − x⋆ ∥2 ≤ 2(f (X0 ) − f (x⋆ )) − 3(f (Xt ) − f (x⋆ )) . 2 dt2 Integrating the above inequality leads to the optimization version of the statement of Lemma 3 stated as [WSF+ 26, Remark 1]. In summary, both [WSF+ 26] and this work derive an inequality of the following form 1 dist2 (qT , q ⋆ ) + 3 2
Z T
1 (T − t)F (qt ) dt ≤ dist2 (q0 , q ⋆ ) + T 2 F(q0 ) , 2 0
(19)
for convex f and log-concave ν X , respectively. The quantities F, dist(·, ·), (qt )t∈[0,T ] , q ⋆ in each result are highlighted below. Domain
F
dist
(qt )t∈[0,T ]
q⋆
This work (Lemma 3)
P2,ac,fs (Rd )
ρ 7→ KL(ρ ∥ ν X )
W2
(ρX t )t∈[0,T ]
νX
[WSF+ 26] (Remark 1)
Rd
x 7→ f (x) − f (x⋆ )
Euclidean
(Xt )t∈[0,T ]
x⋆
Both [WSF+ 26] and this work rely on (19) to derive convergence guarantee for Algorithm 2 and Algorithm 1, respectively. In [WSF+ 26] the authors leveraged the convexity of F when f is convex, and in this work we leveraged the intrinsic convexity of KL divergence (see Eq. (22)). Lower bounding the left-side of Eq. (19) using convexity and dividing both sides by 23 T 2 gives: 1 dist2 (qT , q ⋆ ) + F 3T 2
Z T 0
2(T − t) 1 2 qt dt ≤ dist2 (q0 , q ⋆ ) + F(q0 ) . 2 2 T 3T 3 Z T
(20)
2(T − t) qt dt corresponds to taking weighted T2 0 time-average over trajectory of Hamiltonian flow in Algorithm 2 and randomizing the integration As previously discussed, the weighted average
26
time in Algorithm 1, respectively. Both [WSF+ 26] and this work involve averages q avg (q0 ; T ) and q mix (q0 ; T ) in the notation of the table above defined as Z T 2(T − t) q avg (q0 ; T ) := Et∼TriT [qt ] = qt dt , T2 0 Z 1 T 2(T − t) 1 mix qt dt . q (q0 ; T ) := Et∼Tri f T [qt ] = qT + 2 2 0 T2 In both settings, Eq. (20) yields convergence guarantees that can be succinctly summarized in the following theorems, in the notation of the table above. We first summarize Theorem 1 (for setting [Set-Tri]) and [WSF+ 26, Theorem 1]. Theorem. Let K ∈ N and T > 0. If F is convex and satisfies the quadratic-growth condition: α ∀ q ∈ dom(F) , F(q) ≥ · dist(q, q ⋆ )2 2 then F(q
(K)
)≤
2 2 + 3 3αT 2
K F(q0 )
where q (k) = q avg (q (k−1) ; T ) for k ≥ 1 and q (0) = q0 . We remark that the quadratic growth condition above is referred to by the same name for F ← f (·) − f (x⋆ ) and dist being the Euclidean distance. On the other hand, when F ← KL(· ∥ ν X ) and dist is the Wasserstein distance, this is equivalent to ν X satisfying the α-Talagrand inequality. Similarly, we can summarize the results of Theorem 2 and [WSF+ 26, Theorem 2] as follows. 3 T with T0 = 1. If F is Theorem. Let K ∈ N and {Tk }k≥0 be a sequence such that Tk ≥ 2√ 2 k−1 convex, then K 8 1 (K) ⋆ 2 F(q ) ≤ F(q0 ) + dist(q0 , q ) 9 3
where q (k) = q mix (q (k−1) ; Tk ) for k ≥ 1 and q (0) = q0 . In summary, the analysis and results we present in this paper for RHMC follow by translating the continuous-time analysis and results from [WSF+ 26] to the space of probability measures equipped with the Wasserstein distance. The work [WSF+ 26] was also able to show accelerated rates in discrete time for optimization, and it would be interesting to study how to translate their analysis to a time-discretization of RHMC, which we leave for future work.
C
Additional preliminaries
C.1
General notations
We work on the Euclidean state space Rd of dimension d ∈ N, or the phase space R2d . Let [d] := {1, . . . , d}. For vectors u, v ∈ Rd with u = (u1 , . . . , ud )⊤ and v = (v1 , . . . , vd )⊤ , we denote Pd ⊤v = their ℓ2 -inner product by ⟨u, v⟩ = u i=1 ui vi . We denote the ℓ2 -norm of u by ∥u∥ = qP p d 2 ⟨u, u⟩ = i=1 ui . Let Id ∈ Rd×d denote the identity matrix. We say a matrix A ∈ Rd×d is positive semi-definite, denoted by A ⪰ 0, if A is symmetric and u⊤ Au ≥ 0 for all u ∈ Rd . For a symmetric matrix 27
A ∈ Rd×d with entries A = (Aij )di,j=1 and eigenvalues λ1 , . . . , λd ∈ R, the trace of A is Tr(A) = Pd Pd Qd d×d , the i=1 Aii = i=1 λi , and the determinant of A is det(A) = i=1 λi . For a matrix A ∈ R operator norm of A is ∥A∥op = sup∥u∥=1 ∥Au∥. In particular, if A is positive semi-definite with eigenvalues λ1 , . . . , λd ≥ 0, then ∥A∥op = max{λi : i ∈ [d]}. For A, B ∈ Rd×d , we write A ⪰ B to denote A − B ⪰ 0, and denote their Frobenius inner product by ⟨A, B⟩F := Tr(A⊤ B). For a twice-differentiable function f : Rd → R, we denote the gradient and Hessian map as ∇f : Rd → Rd and ∇2 f : Rd → Rd×d respectively. The Laplacian of f is ∆f (x) = Tr(∇2 f (x)) ∈ R. For a differentiable vector field v : Rd → Rd , the Jacobian of v at x is denoted by ∇v(x) ∈ Rd×d , and the divergence is defined as (∇ · v)(x) = Tr(∇v(x)) ∈ R. For a map (x, y) 7→ f (x, y), we use ∇x f (x, y) and ∇y f (x, y) to denote the partial derivative with respect to x and y respectively while keeping the other fixed. For a time-dependent vector field vt : Rd → Rd , let ∂t vt (x) denote the time derivative vector at a fixed x ∈ Rd : (∂t vt (x))i = ∂(vt∂t(x))i . We also write ∂t vt (x) = v̇t (x). For r ∈ N, we use C r to denote functions with continuous derivatives up to order r, and C ∞ to denote smooth (infinitely-differentiable) functions. The subscript c means compact support; for example, Cc∞ is the class of smooth and compactly supported functions. The subscript b means boundedness, for example, Cb1 is the class of continuously differentiable and uniformly bounded functions.
C.2
Further discussion of probability distributions and statistical distances
We review additional definitions and facts related to probability distributions that are not introduced in the main text. Let P2 (Rd ) denote the space of probability distributions on Rd with finite second moment, so Eρ [∥X∥2 ] < ∞ for all ρ ∈ P2 (Rd ). For ρ ∈ P2 (Rd ), let Covρ (X) = Eρ [(X − µ)(X − µ)⊤ ] ∈ Rd×d denote its covariance matrix, where µ = Eρ [X] ∈ Rd is its mean vector. The variance of ρ is Varρ (X) = Tr(Covρ (X)) = Eρ [∥X − µ∥2 ] ∈ R, and note Varρ (X) ≤ Eρ [∥X∥2 ] < ∞. For a measurable map S : Rd → Rd and a probability distribution ρ, the pushforward distribution S# ρ is defined by (S# ρ)(A) = ρ(S −1 (A)) for any measurable set A ⊆ Rd . If S is a diffeomorphism, i.e., both S and S −1 are continuously differentiable and ρ has a density, then by the change-of-variable formula, the density of the pushforward distribution S# ρ is given by (S# ρ)(x) = ρ(S −1 (x)) det ∇S −1 (x) . Entropy. For ρ ∈ P2,ac (Rd ), the (differential) entropy of ρ is defined as the following, with the convention that 0 log 0 := 0: Z Ent(ρ) := −Eρ [log ρ] = − ρ(x) log ρ(x) dx . Rd
We recall the property of the Gaussian distribution as the maximum entropy distribution for a given covariance matrix [CT06, Theorem 8.6.5]. For any ρ ∈ P2,ac (Rd ), we have 2πeVarρ (X) d 1 d Ent(ρ) ≤ Ent (N (0, Covρ (X))) = log(2πe) + log det Covρ (X) ≤ log . 2 2 2 d In particular, if ρ ∈ P2,ac (Rd ), then Ent(ρ) < ∞.
28
KL divergence. For ρ, ν ∈ P2,ac,fs (Rd ) with ρ ≪ ν, recall the Kullback–Leibler (KL) divergence or relative entropy of ρ with respect to ν is defined as Z h ρ(x) ρi ρ(x) log = dx , KL(ρ ∥ ν) := Eρ log ν ν(x) d R and we define KL(ρ ∥ ν) := ∞ if ρ ̸≪ ν. We recall the following properties of KL divergence that we use in this work. First, we recall the chain rule for KL divergence, see [PW24, Theorem 2.15] for a review. Given joint probability distributions ρXY , ν XY ∈ P2,ac,fs (R2d ), let ρX and ν X denote their X-marginals, and ρY |X=x and ν Y |X=x denote their conditional distributions of Y given X = x, so we can factorize ρXY (x, y) = ρX (x) · ρY |X=x (y) and ν XY (x, y) = ν X (x) · ν Y |X=x (y). Then we have the following decomposition (chain rule): Z XY XY X X KL ρY |X=x ∥ ν Y |X=x ρX (x) dx . (21) KL ρ ∥ν = KL ρ ∥ ν + Rd
Second, we recall the KL divergence is jointly convex in both argument, see [PW24, Section 5.1] for a review. In particular, it implies the following. Given ρXY ∈ P2,ac,fs (R2d ) with a facR torization ρXY (x, y) = ρX (x) · ρY |X=x (y) as above, with Y -marginal ρY (y) = Rd ρXY (x, y) dx = R Y |X=x (y)ρX (x) dx = E Y |X=x (y)], and for any ν ∈ P d 2,ac,fs (R ), by Jensen’s inequality x∼ρX [ρ Rd ρ we have: h i (22) KL ρY ∥ ν ≤ Ex∼ρX KL ρY |X=x ∥ ν . Third, we recall the KL divergence is lower semicontinuous in weak convergence [PW24, Theorem 4.9]. We will review the definition of weak convergence in Section C.3. In particular, let (Pn )n∈N and (Qn )n∈N be sequences of probability measures such that KL(Pn ∥ Qn ) < ∞ for every n. If Pn and Qn converge weakly to probability measures P and Q, respectively, then KL(P ∥ Q) ≤ lim inf KL(Pn ∥ Qn ) . n→∞
(23)
Relative Fisher information. For probability distributions ρ, ν ∈ P2,ac,fs (Rd ) with ρ ≪ ν and differentiable density functions, we recall the relative Fisher information of ρ with respect to ν is Z ρ 2 ρ(x) 2 FI(ρ ∥ ν) := Eρ ∇ log dx . = ρ(x) ∇ log ν ν(x) Rd Wasserstein-2 distance. Recall the Wasserstein-2 distance between probability distributions ρ, ν ∈ P2 (Rd ) is defined by: 1/2 W2 (ρ, ν) = inf E ∥X − Y ∥2 , ω∈Π(ρ,ν)
where the infimum is taken over all couplings between ρ and ν, i.e., joint distributions of (X, Y ) ∼ ω with the correct marginal distributions X ∼ ρ and Y ∼ ν. If ρ ∈ P2,ac (Rd ), then Brenier’s theorem guarantees the existence of a unique optimal transport R map T : Rd → Rd pushing ρ forward to ν, i.e., T# ρ = ν, so that W22 (ρ, ν) = Rd ∥x−T (x)∥2 ρ(x) dx; moreover, for ρ-a.e. x, T (x) = ∇φ(x) for some convex function φ : Rd → R. If ρ, ν ∈ P2,ac (Rd ), then they satisfy the change-of-variable formula (Monge–Ampere equation) for ρ-a.e. x [Vil09, Example 11.2] ρ(x) = ν(T (x)) det(∇T (x)) . 29
Total variation distance. P2,ac,fs (Rd ) is defined by:
The total variation distance between probability distributions ρ, ν ∈
1 TV(ρ, ν) = sup |ρ(A) − ν(A)| = 2 A⊆Rd
C.3
Z |ρ(x) − ν(x)| dx . Rd
Convergence of probability distributions
Throughout this subsection, let (ρn )n∈N ⊂ P2,ac (Rd ) and ρ ∈ P2,ac (Rd ). We say that ρn converges weakly to ρ, if for every bounded continuous function φ : Rd → R, Z Z φ(x)ρ(x) dx . (24) φ(x)ρn (x) dx = lim n→∞ Rd
Rd
We say that ρn converges in total variation to ρ if lim TV(ρn , ρ) = 0 .
n→∞
(25)
We say that ρn converges in Wasserstein-2 distance to ρ if lim W2 (ρn , ρ) = 0 .
n→∞
(26)
We recall that ρn converges in total variation to ρ implies ρn convergence weakly to ρ. We recall from [Vil09, Theorem 6.9] that ρn converges in Wasserstein-2 distance to ρ if and only if ρn converges to ρ weakly and in second moment: Z Z 2 ∥x∥ ρn (x) dx → ∥x∥2 ρ(x) dx . Rd
D
Rd
Details on properties of the Hamiltonian flow
In this appendix we review the properties of the Hamiltonian flow. Recall from Eq. (1) the Hamiltonian of interest: 1 H(x, y) = f (x) + ∥y∥2 , 2 for (x, y) ∈ Rd × Rd . We refer to R2d = Rd × Rd as the phase space, and write z = (x, y), z(0) = (X0 , Y0 ), and z(t) = (Xt , Yt ) interchangeably. We introduce the skew-symmetric symplectic matrix 0 I Ω := . −I 0 Recall from (HF) that the Hamiltonian dynamics can be written equivalently as Ẋt = Yt , Ẏt = −∇f (Xt ), or, in phase-space form, ż(t) = Ω∇H(z(t)). It is convenient to denote solutions of (HF) by the corresponding flow map, and we use this notation throughout the appendix. 30
Definition 1. Whenever (HF) admits a continuously differentiable solution for all t ∈ R and (X0 , Y0 ) ∈ R2d , we denote the corresponding Hamiltonian flow map by Ψ : R × R2d → R2d ,
(t, X0 , Y0 ) 7→ Ψ(t, X0 , Y0 ) ≡ Ψt (X0 , Y0 ) .
That is, Ψt (X0 , Y0 ) = (Xt , Yt ), where (Xt , Yt ) is the solution of (HF) at time t with initial condition (X0 , Y0 ). We also define the marginal variables ΨX t (X0 , Y0 ) = ΠX (Ψt (X0 , Y0 )) ,
ΨYt (X0 , Y0 ) = ΠY (Ψt (X0 , Y0 )) ,
where ΠX : R2d → Rd and ΠY : R2d → Rd denote projections onto the first and second variables, respectively: ΠX (x, y) = x and ΠY (x, y) = y. Proposition 1. Let f be a L-smooth function. Then, the vector field b(x, y) := (y, −∇f (x)) is (1 + L)-Lipschitz. Additionally, by the Picard–Lindelöf theorem, there exists a unique solution at any time t for HF from any initial condition. Proof. We show that b(x, y) is (1 + L)-Lipschitz. For any z = (x, y) and z̃ = (x̃, ỹ) we calculate: ∥b(z) − b(z̃)∥ = ∥(y − ỹ, −∇f (x) + ∇f (x̃))∥ ≤ ∥y − ỹ∥ + L∥x − x̃∥ ≤ (L + 1)∥z − z̃∥ .
(27)
As a consequence of the above proposition, the Hamiltonian flow satisfies the flow property Ψ0 = Id ,
Ψt+s = Ψt ◦ Ψs ,
Ψ−1 t = Ψ−t .
Roadmap. We briefly describe the organization of this appendix. In Section D.1, we review the properties of the trajectory of the Hamiltonian flow (Xt , Yt ) in the phase space. In Section D.2, XY evolve according we randomize the initial condition by taking (X0 , Y0 ) ∼ ρXY 0 , let (Xt , Yt ) ∼ ρt to the Hamiltonian flow, and review the resulting distributional properties of ρXY t . Finally, in along Section D.3, we derive the differential equations governing the evolution of moments of ρXY t the Hamiltonian flow.
D.1
Deterministic properties of the Hamiltonian flow
D.1.1
Conservation of the Hamiltonian function
Lemma 5. Let (Xt , Yt )t≥0 evolve along (HF) with initial condition (X0 , Y0 ). For any t ∈ R and (X0 , Y0 ) ∈ R2d , H(Xt , Yt ) = H(X0 , Y0 ) . Proof. Let z(t) = Ψt (z0 ). By the chain rule, d H(z(t)) = ∇H(z(t))⊤ ż(t) = ∇H(z(t))⊤ Ω∇H(z(t)) = 0 , dt where the last equality follows from the anti-symmetry of Ω.
31
D.1.2
Conservation of volume
Recall Ψt : R2d → R2d is the Hamiltonian flow map that sends the initial state (X0 , Y0 ) ∈ R2d at time 0 to the solution (Xt , Yt ) ∈ R2d of the Hamiltonian flow (HF) at time t ∈ R. Lemma 6. For any t ∈ R and (x, y) ∈ R2d , we have det (∇Ψt (x, y)) = 1 . Consequently, for any measurable set A ⊆ R2d , vol(Ψt (A)) = vol(A). This is standard property of Hamiltonian flow, but we include a proof for the readers’ convenience. Proof. Fix z = (x, y) ∈ R2d , and write J(t) := ∇Ψt (z) with J(0) = I2d . The Hamiltonian flow (HF) dynamics is Ψ̇t (z) = Ω∇H(Ψt (z)) . Taking gradient of this relation with respect to z and using the chain rule gives ˙ = ∇ Ψ̇t (z) = Ω∇2 H(Ψt (z))J(t) . J(t) Define M (t) := J(t)⊤ Ω J(t) . By the product rule, ˙ ⊤ ΩJ(t) + J(t)⊤ ΩJ(t) ˙ Ṁ (t) = J(t) h ⊤ i = J(t)⊤ Ω∇2 H(Ψt (z)) Ω + Ω Ω∇2 H(Ψt (z)) J(t) h i = J(t)⊤ −∇2 H(Ψt (z)) Ω2 + Ω2 ∇2 H(Ψt (z)) J(t) h i = J(t)⊤ ∇2 H(Ψt (z)) − ∇2 H(Ψt (z)) J(t) =0 where in the computation above we have used the fact that ∇2 H(Ψt (z)) is symmetric, Ω⊤ = −Ω, and Ω2 = −I2d . Therefore Ṁ (t) = 0. Since M (0) = I⊤ 2d Ω I2d = Ω, we get J(t)⊤ ΩJ(t) = Ω
∀ t ∈ R.
Taking determinants on both sides gives det(J(t))2 · det(Ω) = det(J(t)⊤ ) det(Ω) det(J(t)) = det(Ω) , 2 and since det(Ω) ̸= 0, this implies det J(t) = 1. Since t 7→ det J(t) is continuous and det J(0) = 1, we conclude det J(t) = 1 for all t ∈ R. Finally, the change-of-variables formula yields Z Z vol(Ψt (A)) = |det(∇Ψt (z))| dz = 1 dz = vol(A) . A
A
32
D.1.3
Bi-Lipschitzness of Hamiltonian flow map
If f is smooth, then the Hamiltonian flow map is a bi-Lipschitz map. Lemma 7. If f is L-smooth, then for all t ∈ R, z = (x, y) ∈ R2d and z̃ = (x̃, ỹ) ∈ R2d , e−(1+L)|t| ∥z − z̃∥ ≤ ∥Ψt (z) − Ψt (z̃)∥ ≤ e(1+L)|t| ∥z − z̃∥ . Furthermore, ∥∇Ψt (z)∥op ≤ e(1+L)|t| . Proof. We show the estimate for t ≥ 0; the case t ≤ 0 follows by applying the same argument to the backward flow. Let b(x, y) := Ω∇H(x, y) = (y, −∇f (x)). Since f is L-smooth, ∇f is L-Lipschitz, and recall from Eq. (27) that b is therefore (1 + L)-Lipschitz. Fix z0 , z̃0 ∈ R2d . By the integral form of the Hamiltonian flow: Z t Ψt (z0 ) − Ψt (z̃0 ) = z0 − z̃0 + (b(Ψs (z0 )) − b(Ψs (z̃0 ))) ds . 0
Therefore, Z t ∥b(Ψs (z0 )) − b(Ψs (z̃0 ))∥ ds
∥Ψt (z0 ) − Ψt (z̃0 )∥ ≤ ∥z0 − z̃0 ∥ + 0
Z t ≤ ∥z0 − z̃0 ∥ +
(1 + L)∥Ψs (z0 ) − Ψs (z̃0 )∥ ds . 0
By Grönwall’s inequality, ∥Ψt (z0 ) − Ψt (z̃0 )∥ ≤ e(1+L)t ∥z0 − z̃0 ∥ ,
t ≥ 0.
On the other hand, note that the backward flow s 7→ Ψ−s is generated by the vector field −b, which is also (1 + L)-Lipschitz. Hence, the same Grönwall argument gives, for all u, v ∈ R2d and t ≥ 0, ∥Ψ−t (u) − Ψ−t (v)∥ ≤ e(1+L)t ∥u − v∥ . Applying this inequality with u = Ψt (z0 ) and v = Ψt (e z0 ), and using Ψ−t = (Ψt )−1 , yields ∥z0 − ze0 ∥ = ∥Ψ−t (Ψt (z0 )) − Ψ−t (Ψt (z̃0 ))∥ ≤ e(1+L)t ∥Ψt (z0 ) − Ψt (z̃0 )∥ . Therefore, ∥Ψt (z0 ) − Ψt (e z0 )∥ ≥ e−(1+L)t ∥z0 − ze0 ∥ . The flow map z → Ψt (z) is continuously differentiable with respect to z under L-smoothness of f . Therefore for any z ∈ R2d and v ∈ R2d , ∥∇Ψt (z)v∥ = lim
ε→0
Ψt (z + εv) − Ψt (z) ≤ e(L+1)t ∥v∥ . ε
Taking supremum over ∥v∥ ≤ 1 yields ∥∇Ψt (z)∥op ≤ e(L+1)t .
33
D.2
Distributional properties of the Hamiltonian flow
In this section, suppose we run the Hamiltonian flow (HF) from an initial joint random variable (X0 , Y0 ) ∼ ρXY for some ρXY ∈ P2,ac,fs (R2d ), to obtain (Xt , Yt ) ∼ ρXY for all t ∈ R. Recall t 0 0 2d 2d 2d Ψt : R → R is the flow map that sends (X0 , Y0 ) ∈ R to the solution (Xt , Yt ) ∈ R2d of the Hamiltonian flow (HF) at time t ∈ R, so we have ρXY = (Ψt )# ρXY t 0 . D.2.1
Conservation of entropy
We first note the Hamiltonian flow conserves entropy. This is the distributional analog of the fact that the Hamiltonian flow conserves volume (Lemma 6). Lemma 8. For all t ∈ R, Ent ρXY = Ent ρXY . t 0 Proof. Since (Ψt )−1 = Ψ−t , by the change-of-variable formula for ρXY = (Ψt )# ρXY t 0 , we have for 2d all (x, y) ∈ R , XY ρXY (Ψ−t (x, y)) · |det ∇Ψ−t (x, y)| = ρXY t (x, y) = ρ0 0 (Ψ−t (x, y)) ,
note that |det ∇Ψ−t (x, y)| = 1 follows from Lemma 6. Therefore, XY − log ρXY t (x, y) = − log ρ0 (Ψ−t (x, y)) . XY Now let (Xt , Yt ) ∼ ρXY t , so Ψ−t (Xt , Yt ) = (X0 , Y0 ) ∼ ρ0 . Taking expectation on both sides gives
Ent ρXY = E(Xt ,Yt )∼ρXY − log ρXY t t (Xt , Yt ) t XY = E(Xt ,Yt )∼ρXY − log ρ (Ψ (X , Y )) −t t t 0 t = E(X0 ,Y0 )∼ρXY − log ρXY (X0 , Y0 ) 0 0 = Ent ρXY . 0
D.2.2
Conservation of KL divergence to the joint distribution
Since the Hamiltonian flow conserves both the Hamiltonian function and entropy, it also conserves the KL divergence to the joint distribution ν XY ∝ exp(−H). R Proof of Lemma 1. Let ZH := R2d exp(−H(x, y)) dx dy, so ν XY (x, y) = exp(−H(x, y))/ZH . Lemma 8 gives Ent ρXY = Ent ρXY . Furthermore, by Lemma 5, H(Xt , Yt ) = H(X0 , Y0 ). So for (X0 , Y0 ) ∼ t 0 XY XY ρ0 , (Xt , Yt ) ∼ ρt , we also have EρXY [H] = EρXY [H]. Therefore, t 0
KL ρXY ∥ ν XY = −Ent ρXY + EρXY [H] + log ZH t t t = −Ent ρXY + EρXY [H] + log ZH 0 0 XY XY = KL ρ0 ∥ ν .
34
D.2.3
Descent property of KL divergence in the X-marginal
Proof of Lemma 2. For t ∈ R and the solution (Xt , Yt ) ∼ ρXY along the Hamiltonian flow, let ρX t t Y |X=x denote the marginal distribution of Xt , and let ρt denote the conditional distribution of Yt given Xt = x. By Lemma 1, we know that KL ρXY ∥ ν XY = KL ρXY ∥ ν XY . t 0 XY = ν X ⊗ γ, we have Since ρXY = ρX 0 0 ⊗ γ and ν X KL ρXY ∥ ν XY = KL ρX . 0 0 ∥ν
Furthermore, by the chain rule for KL divergence, we have h i Y |X=x X KL ρXY ∥ ν XY = KL ρX + Ex∼ρX KL ρ ∥ γ t t ∥ν t t X X = KL ρt ∥ ν + KL ρXY ∥ ρX t t ⊗γ X ≥ KL ρX , t ∥ν
(28)
where the last inequality follows by dropping the KL divergence of the conditional distributions. Combining the three calculations above, we obtain the desired result: X X KL ρX = KL ρXY ∥ ν XY = KL ρXY ∥ ν XY ≥ KL ρX . 0 ∥ν 0 t t ∥ν
D.2.4
Continuity of KL divergence
XY ) < ∞. Lemma 9. Assume ν X ∝ exp(−f ) is log-smooth. Let ρXY ∈ P2,ac,fs (R2d ) with KL(ρXY 0 0 ∥ν be the solution to the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρXY For t ∈ R, let (Xt , Yt ) ∼ ρXY t 0 , X ) is continuous. X ∥ ν and let Xt ∼ ρt denote the X-marginal. Then t 7→ KL(ρX t
Proof. Since ν X is log-smooth, there exists 0 < L < ∞ such that f is L-smooth. Consequently, for any (X0 , Y0 ) ∈ R2d , the solution (Xt , Yt ) ∈ R2d of HF starting from X0 , Y0 is uniquely defined for all t ∈ R and the map t 7→ (Xt , Yt ) is continuous with respect to t. Equivalently, we know that the flow map Ψt (x, y) is continuous in t for all (x, y) ∈ R2d . Therefore, for all sequences such that {tn } → t and for any continuous and bounded test function φ : R2d → R, we have Z Z XY XY lim φ(x, y)ρtn (x, y) dx dy = lim φ(ΨXY tn (x, y))ρ0 (x, y) dx dy n→∞ R2d n→∞ R2d Z XY = φ(ΨXY t (x, y))ρ0 (x, y) dx dy R2d Z = φ(x, y)ρXY t (x, y) dx dy . XY (see Section C.3), or equivalently we may Therefore, for all tn → t, ρXY tn converges weakly to ρt XY say that ρt is weakly continuous in t. Taking marginal gives that t 7→ ρX t is also weakly continuous in t. The remainder of this proof aims to show that X X X lim KL(ρX tn ∥ ν ) = KL(ρt ∥ ν )
n→∞
35
X for any arbitrary sequence {tn } that converges to t. Since ρX tn converges weakly to ρt , by the lower semicontinuity of KL divergence (see Eq. (23)), X X X KL(ρX . (29) t ∥ ν ) ≤ lim inf KL ρtn ∥ ν n→∞
XY X X We also know ρXY tn converges weakly to ρt , and ρtn ⊗ γ converges weakly to ρt ⊗ γ, so by XY X XY the lower semicontinuity of KL divergence, KL ρt ∥ ρt ⊗ γ ≤ lim inf n→∞ KL ρtn ∥ ρX tn ⊗ γ . Equivalently, XY X XY X −KL ρXY ∥ ρX (30) t t ⊗ γ ≥ − lim inf KL ρtn ∥ ρtn ⊗ γ = lim sup −KL ρtn ∥ ρtn ⊗ γ . n→∞
n→∞
From Eq. (28), we have X KL ρX = KL ρXY ∥ ν XY − KL ρXY ∥ ρX t ∥ν 0 t t ⊗γ (30)
X ≥ KL ρXY ∥ ν XY + lim sup −KL ρXY 0 tn ∥ ρtn ⊗ γ n→∞ XY XY X = lim sup KL ρ0 ∥ ν − KL ρXY tn ∥ ρtn ⊗ γ n→∞ X = lim sup KL ρX . tn ∥ ν n→∞
where the last step follows from applying the Eq. (28) identity at time tn . By combining this with X X X X Eq. (29), we conclude limn→∞ KL ρtn ∥ ν = KL ρt ∥ ν . Since (tn )n∈N is arbitrary, this shows X is continuous. that t 7→ KL ρX ∥ ν t D.2.5
Finiteness of second moment
We first show that the normalizability of ν X ∝ e−f , together with the L-smoothness of f , implies that f is bounded below. Lemma 10. Assume ν X ∝ exp(−f ) is L-log-smooth for some L ∈ (0, ∞), and assume Zf := R d Rd exp(−f (u)) du < ∞. Then, for all x ∈ R , d 2π f (x) ≥ log − log Zf . 2 L In particular, there exists x⋆ ∈ Rd such that f (x⋆ ) = inf x∈Rd f (x) > −∞. Proof. By the L-smoothness of f , for all x, u ∈ Rd we have: f (u) ≤ f (x) + ⟨∇f (x), u − x⟩ +
L ∥u − x∥2 . 2
Therefore, for any x ∈ Rd , Z Zf = exp(−f (u)) du d ZR L 2 ≥ exp −f (x) − ⟨∇f (x), u − x⟩ − ∥u − x∥ du 2 Rd ! Z ∥∇f (x)∥2 L ∇f (x) 2 = exp −f (x) + exp − · u−x+ du 2L 2 L Rd 36
d/2 2π ∥∇f (x)∥2 · = exp −f (x) + , 2L L where the last step holds by evaluating the Gaussian integral. Taking logarithm on both sides gives: ∥∇f (x)∥2 d 2π − log Zf + f (x) ≥ log > −∞ . 2 L 2L Therefore, f is uniformly lower bounded. By rearranging the inequality above, we get: 1 ν (x) = exp(−f (x)) ≤ Zf X
L 2π
d/2
∥∇f (x)∥2 . exp − 2L
Therefore, ∥∇ν X (x)∥ = ν X (x)∥∇f (x)∥ ≤ r2 because supr∈R r exp(− 2L ) =
q
L e.
L 2π
d/2
d/2 r ∥∇f (x)∥2 L L exp − ∥∇f (x)∥ ≤ , 2L 2π e
Hence, ν X (x) is globally Lipschitz, and the finiteness of its
second moment implies that lim∥x∥→∞ ν X (x) = 0. Since ν X (x) = Z1f e−f (x) , this implies lim f (x) = − log(Zf ) − log(ν X (x)) = ∞ .
∥x∥→∞
Then there exists R ∈ (0, ∞) such that f (x) > f (0) whenever ∥x∥ > R. Since f is continuous, it attains its minimum on the compact set BR (0) at some x⋆ ∈ Rd . Because 0 ∈ BR (0), this also implies f (x⋆ ) ≤ f (0) < f (x) for every ∥x∥ > R. This shows x⋆ is a global minimizer of f . Next, we show the solution of the Hamiltonian flow remains in P2,ac,fs (Rd ), and the second moment grows at most quadratically in time. Lemma 11. Assume ν X ∝ exp(−f ) is L-log-smooth for some L ∈ (0, ∞). For any finite t ∈ R, XY ∈ let (Xt , Yt ) ∼ ρXY be the solution to the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρXY t 0 . If ρ0 2d XY 2d X d P2,ac,fs (R ), then ρt ∈ P2,ac,fs (R ) and ρt ∈ P2,ac,fs (R ). In particular, h i 2 2 2 E(Xt ,Yt )∼ρXY ∥(X , Y )∥ ≤ (2t + 1)(2L + 2)E + (4t2 + 2)L∥x⋆ ∥2 , XY ∥(X0 , Y0 )∥ t t (X ,Y )∼ρ 0 0 t 0
where x⋆ ∈ Rd satisfies f (x⋆ ) = inf x∈Rd f (x). Proof. Since the Hamiltonian flow map Ψt is a diffeomorphism and ρXY ∈ P2,ac,fs (R2d ), we have 0 ρXY = (Ψt )# ρXY is also absolutely continuous with respect to Lebesgue measure on R2d and has t 0 full support with positive density. Therefore, its X-marginal ρX t is also absolutely continuous with 2d respect to Lebesgue measure on R and has full support with positive density. We next verify the finiteness of the second moment for the case t ≥ 0; the other case follows similarly. For any (X0 , Y0 ) ∈ R2d , let (Xs , Ys ) be the solution of the HF at time s with initial condition (X0 , Y0 ). Recall Lemma 10 guarantees f (x∗ ) = inf x∈Rd f (x) > −∞. By the conservation of the Hamiltonian function, the following holds for any s ∈ R: ∥Ys ∥2 = 2 (H(Xs , Ys ) − f (Xs )) = 2 (H(X0 , Y0 ) − f (Xs )) ≤ 2 (H(X0 , Y0 ) − f (x⋆ )) .
37
For any 0 ≤ t < ∞, by HF, we know that Xt − X0 = 2
Ys ds 0
0 Ys ds. This gives
2
Z t
∥Xt ∥ = X0 +
Rt
≤ 2∥X0 ∥2 + 2t2 max ∥Ys ∥2 ≤ 2∥X0 ∥2 + 4t2 (H(X0 , Y0 ) − f (x⋆ )) . s∈[0,t]
Combining the two estimates and use L-smoothness of f gives: ∥(Xt , Yt )∥2 ≤ 2∥X0 ∥2 + (4t2 + 2) (H(X0 , Y0 ) − f (x⋆ )) L 1 ≤ 2∥X0 ∥2 + (4t2 + 2) ∥X0 − x⋆ ∥2 + ∥Y0 ∥2 2 2 ≤ (4t2 L + 2L + 2)∥X0 ∥2 + (4t2 + 2)L∥x⋆ ∥2 + (2t2 + 1)∥Y0 ∥2 ≤ (2t2 + 1)(2L + 2)∥(X0 , Y0 )∥2 + (4t2 + 2)L∥x⋆ ∥2 . Taking expectation on both sides: h i 2 2 2 + (4t2 + 2)L∥x⋆ ∥2 . E(Xt ,Yt )∼ρXY ∥(X , Y )∥ ≤ (2t + 1)(2L + 2)E XY ∥(X0 , Y0 )∥ t t (X ,Y )∼ρ 0 0 t 0
h
i 2 ∥(X , Y )∥ < ∞, and therefore Since ρXY has finite second moment, E 0 0 0 (X0 ,Y0 )∼ρXY 0
h
i 2 E(Xt ,Yt )∼ρXY ∥(X , Y )∥ < ∞. t t t
D.3
Moments along Hamiltonian flow
In this section, suppose we run the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρXY for some ρXY ∈ 0 0 2d XY XY XY P2,ac,fs (R ), to obtain (Xt , Yt ) ∼ ρt for all t ∈ R, so ρt = (Ψt )# ρ0 . Let Xt ∼ ρX denote the t X-marginal at time t. We introduce some notations for the moments along the Hamiltonian flow. At each t ∈ R, we define the conditional mean of the velocity variable Yt given Xt to be ut : Rd → Rd given by Z 1 ut (x) := E[Yt | Xt = x] = X yρXY (x, y) dy . (31) ρt (x) Rd t We define the conditional second moment of Yt given Xt to be Mt : Rd → Rd×d given by Z 1 ⊤ Mt (x) := E[Yt Yt | Xt = x] = X yy ⊤ ρXY (32) t (x, y) dy . ρt (x) Rd We define the conditional covariance of Yt given Xt to be Σt : Rd → Rd×d given by Σt (x) := Cov(Yt | Xt = x) = Mt (x) − ut (x)ut (x)⊤ .
(33)
Whenever the derivatives below are justified, we define the acceleration field at : Rd → Rd by at (x) := ∂t ut (x) + ∇ut (x) ut (x) .
(34)
We refer to this as the “acceleration field” for the following reason. If Φt (x) solves the ordinary differential equation Φ̇t (x) = ut (Φt (x)), then Φ̈t (x) = ∇ut (Φt (x))ut (Φt (x)) + ∂t ut (Φt (x)) = at (Φt (x)). 38
d Since ρXY ∈ P2,ac,fs (R2d ), we have (Lemma 11) that ρXY ∈ P2,ac,fs (R2d ) and ρX t t ∈ P2,ac,fs (R ) 0 for all t ∈ R. We define the optimal transport map Rt : Rd → Rd that pushes forward ρX t to ν X = (Rt )# ρX and satisfies t h i 2 X W22 (ρX , ν ) = E ∥X − R (X )∥ . X t t t t Xt ∼ρt
We also define the displacement map vt : Rd → Rd by vt (x) := x − Rt (x) .
(35)
We note the definitions of the optimal transport map Rt and the displacement map vt depend on the initial distribution ρXY of the Hamiltonian flow (HF). 0 D.3.1
Continuity equation
We have the following continuity equation for the quantities defined above. Lemma 12. Assume ν X is log-smooth. Along the Hamiltonian flow (HF) (Xt , Yt ) ∼ ρXY from t XY XY 2d (X0 , Y0 ) ∼ ρ0 with ρ0 ∈ P2,ac,fs (R ), we have the following continuity equations (in the sense of distribution) for all t ≥ 0: X ∂t ρX t + ∇ · (ρt ut ) = 0 , X X ∂t (ρX t ut ) + ∇ · ρt Mt + ρt ∇f = 0 .
Remark 1. We recall the notion of an identity to hold in the sense of distributions from [AGS05, Section 8.1]. Fix T > 0. We write Cc∞ (Rd × (0, T )) for the class of all compactly supported smooth functions ϕ : Rd × (0, T ) → R. In general, for a function F that is integrable on every compact subset of Rd × (0, T ), we say the equation F = 0 holds in the sense of distributions on Rd × (0, T ) if Z F(x, t)ϕ(x, t) dx dt = 0 , Rd ×(0,T )
for all ϕ ∈ Cc∞ (Rd × (0, T )). Note that the compact-support property of ϕ makes integration by parts convenient, because no boundary terms appear: ⟨∂t F, ϕ⟩ := −⟨F, ∂t ϕ⟩ ,
⟨∂xi F, ϕ⟩ := −⟨F, ∂xi ϕ⟩ .
We extend these definitions to vector-valued functions componentwise, with scalar multiplication replaced by the Euclidean inner product. XY X Y Proof. Fix T < ∞. Draw z0 = (x0 , y0 ) ∼ ρXY 0 , and write Ψt (z0 ) = (Ψt (z0 ), Ψt (z0 )). We prove Y the first equation. Let ϕ ∈ Cc∞ (Rd × (0, T )). Since Ψ̇X t (z0 ) = Ψt (z0 ), by chain rule:
d X X Y ϕ(ΨX t (z0 ), t) = ∂t ϕ(Ψt (z0 ), t) + ⟨∇x ϕ(Ψt (z0 ), t), Ψt (z0 )⟩ . dt Since ϕ is compactly supported on a compact subset of Rd ×(0, T ), we know that 0 = ϕ(ΨX 0 (z0 ), 0) = X ϕ(ΨT (z0 ), T ). Therefore: X 0 = ϕ(ΨX T (z0 ), T ) − ϕ(Ψ0 (z0 ), 0) =
Z T 0
d ϕ(ΨX t (z0 ), t) dt dt 39
Z T =
X Y ∂t ϕ(ΨX t (z0 ), t) + ⟨∇x ϕ(Ψt (z0 ), t), Ψt (z0 )⟩ dt .
0
Since z0 = (x0 , y0 ) ∼ ρXY 0 , taking expectation on both sides gives: Z Z T X X Y ∂t ϕ(Ψt (z0 ), t) + ⟨∇x ϕ(Ψt (z0 ), t), Ψt (z0 )⟩ dt ρXY 0= 0 (z0 ) dz0 . R2d
0
Since ϕ is compactly supported on a compact subset of Rd ×(0, T ) and is continuously differentiable, there exists a constant Cϕ dependent only on ϕ such that Y X Y ∂t ϕ(ΨX t (z0 ), t) + ⟨∇x ϕ(Ψt (z0 ), t), Ψt (z0 )⟩ ≤ Cϕ (1 + ∥Ψt (z0 )∥) .
Therefore, the integrand is bounded by Cϕ (1 + ∥ΨYt (z0 )∥), and by Lemma 11, Z ∥ΨYt (z0 )∥2 ρXY sup 0 (z0 ) dz0 < ∞ . t∈[0,T ]
R2d
2d Hence the integrand is integrable with respect to ρXY 0 (z0 ) dz0 dt on R × [0, T ]. Fubini’s theorem allows us to switch the order of integrations: Z Z T X X Y 0= ∂t ϕ(Ψt (z0 ), t) + ⟨∇x ϕ(Ψt (z0 ), t), Ψt (z0 )⟩ dt ρXY 0 (z0 ) dz0 R2d
0
Z TZ
= 0 (a)
R2d
Z TZ
=
0 R2d Z TZ
= 0
Rd
XY X Y ∂t ϕ(ΨX t (z0 ), t) + ⟨∇x ϕ(Ψt (z0 ), t), Ψt (z0 )⟩ ρ0 (z0 ) dz0 dt
[∂t ϕ(x, t) + ⟨∇x ϕ(x, t), y⟩] ρXY t (x, y) dx dy dt
[∂t ϕ(x, t) + ⟨∇x ϕ(x, t), ut (x)⟩] ρX t (x) dx dt ,
XY ∞ d where step (a) follows because ρXY = (ΨXY t t )# ρ0 . Since this holds for every ϕ ∈ Cc (R × (0, T )), we obtain X ∂t ρX t + ∇x · (ρt ut ) = 0
in the sense of distributions. Y We now prove the second equation. Let ψ ∈ Cc∞ (Rd ×(0, T ); Rd ), we consider ⟨ψ(ΨX t (z0 ), t), Ψt (z0 )⟩. X Y Y X Since Ψ̇t (z0 ) = Ψt (z0 ) and Ψ̇t (z0 ) = −∇f (Ψt (z0 )), by chain rule: d Y X Y X Y Y ψ(ΨX t (z0 ), t), Ψt (z0 ) = ∂t ψ(Ψt (z0 ), t), Ψt (z0 ) + ∇x ψ(Ψt (z0 ), t)Ψt (z0 ), Ψt (z0 ) dt X − ψ(ΨX t (z0 ), t), ∇f (Ψt (z0 )) . Since ψ is supported on a compact subset of Rd × (0, T ), we know that 0 = ψ(ΨX 0 (z0 ), 0) = ψ(ΨX (z ), T ). Therefore: 0 T Y X Y 0 = ψ(ΨX T (z0 ), T ), ΨT (z0 ) − ψ(Ψ0 (z0 ), 0), Ψ0 (z0 ) Z T d Y = ψ(ΨX t (z0 ), t), Ψt (z0 ) dt 0 dt Z Th i Y X Y Y X X = ∂t ψ(ΨX (z ), t), Ψ (z ) + ∇ ψ(Ψ (z ), t)Ψ (z ), Ψ (z ) − ψ(Ψ (z ), t), ∇f (Ψ (z )) dt . 0 0 x 0 0 0 0 0 t t t t t t t 0
40
Since z0 = (x0 , y0 ) ∼ ρXY 0 , taking expectation on both sides gives: Z Z Th Y X Y Y 0= ∂t ψ(ΨX t (z0 ), t), Ψt (z0 ) + ∇x ψ(Ψt (z0 ), t)Ψt (z0 ), Ψt (z0 ) R2d 0 i X − ψ(ΨX (z ), t), ∇f (Ψ (z )) dt ρXY 0 0 t t 0 (z0 ) dz0 . Since ψ is supported on a compact subset of Rd × (0, T ) and is continuously differentiable, there exists a constant Cψ dependent only on ψ such that Y X Y Y X X ∂t ψ(ΨX t (z0 ), t), Ψt (z0 ) + ∇x ψ(Ψt (z0 ), t)Ψt (z0 ), Ψt (z0 ) − ψ(Ψt (z0 ), t), ∇f (Ψt (z0 )) Y Y 2 ≤ Cψ ∥∇f (ΨX . t (z0 ))∥ + ∥Ψt (z0 )∥ + ∥Ψt (z0 )∥ X By L-smoothness of f or equivalently L-Lipschitzness of ∇f , ∥∇f (ΨX t (z0 ))∥ ≤ L∥Ψt (z0 )∥ + Y Y 2 ∥∇f (0)∥. Therefore, the integrand is bounded by Cψ 1 + ∥∇f (0)∥ + ∥ΨX t (z0 )∥ + ∥Ψt (z0 )∥ + ∥Ψt (z0 )∥ . By Lemma 11, Z
sup t∈[0,T ]
R2d
2 Y 2 XY (∥ΨX t (z0 )∥ + ∥Ψt (z0 )∥ )ρ0 (z0 )dz0 < ∞ .
2d × [0, T ]. Fubini’s theorem Hence the integrand is integrable with respect to ρXY 0 (dz0 ) dt on R allows switching the order of integration: Z Z Th Y X Y Y 0= ∂t ψ(ΨX t (z0 ), t), Ψt (z0 ) + ∇x ψ(Ψt (z0 ), t)Ψt (z0 ), Ψt (z0 ) 2d R 0 i X dt ρXY − ψ(ΨX (z ), t), ∇f (Ψ (z )) 0 0 0 (z0 )dz0 t t Z TZ h Y X Y Y = ∂t ψ(ΨX t (z0 ), t), Ψt (z0 ) + ∇x ψ(Ψt (z0 ), t)Ψt (z0 ), Ψt (z0 ) 0 R2d i X ρXY − ψ(ΨX 0 (z0 )dz0 dt t (z0 ), t), ∇f (Ψt (z0 )) Z TZ h i (a) = ⟨∂t ψ(x, t), y⟩ + ⟨∇x ψ(x, t)y, y⟩ − ⟨ψ(x, t), ∇f (x)⟩ ρXY t (x, y) dx dy dt 0 R2d Z TZ h i (b) = ⟨∂t ψ(x, t), ut (x)⟩ + ⟨∇x ψ(x, t), Mt (x)⟩F − ⟨ψ(x, t), ∇f (x)⟩ ρX t (x) dx dt . 0
Rd
XY and use In the above, step (a) follows from ρXY = (Ψt )# ρXY t 0 , and in step (b), we disintegrate ρt the fact that ⟨∇x ψ(x, t)y, y⟩ = ∇x ψ(x, t), yy ⊤ F and Mt (x) = E[Yt Yt⊤ | Xt = x]. Since this holds for every ψ ∈ Cc∞ (Rd × (0, T ); Rd ), we obtain X X ∂t (ρX t ut ) + ∇x · (ρt Mt ) + ρt ∇f = 0
in the sense of distributions.
E
Properties of the Hamiltonian flow under regular initialization
We derive properties of the Hamiltonian flow under the assumption that the initial distribution is “regular”, which means it is warm and smooth relative to the target distribution; see Assumption 1. We remark that the design of Assumption 1 is inspired by [Lu26, Definition 2.5]. Throughout, we assume ν X ∈ P2,ac,fs (Rd ) is L-log-smooth. 41
d d Assumption 1. Let ρX 0 ∈ P2,ac,fs (R ) be the initial X-distribution, and define q0 : R → R by X q0 (x) := ρX 0 (x)/ν (x). We assume q0 is continuously differentiable with bounded first derivative, and there exist constants 0 < ζ ≤ ξ < ∞ such that for all x ∈ Rd ,
ζ ≤ q0 (x) ≤ ξ .
E.1
Propagation of warmness
Recall one step of HMC starts from (X0 , Y0 ) ∼ ρXY = ρX 0 0 ⊗γ, evolves via the Hamiltonian flow (HF) XY for some time t ≥ 0 to reach (Xt , Yt ) ∼ ρt , and returns Xt ∼ ρX t as the next iterate. Since the X , then ρX = ν X ), the warmness target distribution ν X is conserved under this operation (if ρX = ν t 0 of the X-iterate is propagated. XY be Lemma 13. Assume ρX 0 satisfies Assumption 1 for some 0 < ζ ≤ ξ < ∞. Let (Xt , Yt ) ∼ ρt XY X X the solution of the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρ0 = ρ0 ⊗ γ, and let Xt ∼ ρt be the X-iterate. Define ρXY (x, y) ρX (x) , gt (x, y) := tXY . qt (x) := tX ν (x) ν (x, y)
Then for all t ≥ 0 and (x, y) ∈ R2d , we have ζ ≤ qt (x) ≤ ξ,
ζ ≤ gt (x, y) ≤ ξ .
Proof. Recall Ψt is the Hamiltonian flow map at time t, so ρXY = (Ψt )# ρXY t 0 . Then by the change2d of-variable formula, for all (x, y) ∈ R : XY XY ρXY (x, y) = | det ∇Ψ−t (x, y)| · ρXY t (x, y) = (Ψt )# ρ0 0 (Ψ−t (x, y)) = ρ0 (Ψ−t (x, y)) , since | det ∇Ψ−t (x, y)| = 1 along Hamiltonian flow (Lemma 6). On the other hand, since ν XY is stationary along the Hamiltonian flow, we have ν XY (x, y) = ν XY (Ψ−t (x, y)) . XY = ν X ⊗ γ, we have by Assumption 1, Since ρXY = ρX 0 0 ⊗ γ and ν
gt (x, y) =
X ρX ρXY ρXY 0 (Ψ−t (x, y)) t (x, y) 0 (Ψ−t (x, y)) = = ∈ [ζ, ξ] . ν XY (x, y) ν XY (Ψ−t (x, y)) ν X (ΨX −t (x, y))
Integrating over y ∈ Rd gives Z Z Z ρX ρXY ρXY t (x) t (x, y) t (x, y) qt (x) = X = dy = γ(y) dy = gt (x, y)γ(y) dy ∈ [ζ, ξ] . XY (x, y) ν (x) ν X (x) Rd Rd ν Rd
E.2
Implications of regular initialization
From a regular initialization (Assumption 1), we derive the following regularity properties of the iterates along the Hamiltonian flow (HF). We recall the definitions of the conditional mean ut (31), conditional second moment Mt (32), conditional covariance Σt (33), acceleration field at (34), and displacement field vt (35). 42
Lemma 14. Assume ν X is log-smooth, and ρX 0 satisfies Assumption 1 for some 0 < ζ ≤ ξ < ∞. XY Let (Xt , Yt ) ∼ ρt be the solution of the Hamiltonian flow (HF) from (X0 , Y0 ) ∼ ρXY = ρX 0 0 ⊗ γ, X and let Xt ∼ ρt be the X-iterate. Then for each T ∈ (0, ∞), the following regularity conditions hold for all t ∈ (0, T ): (R1) ρX t (x), ut (x), and Mt (x) are continuously differentiable with respect to t and x. (R2) The following tail integrability condition holds: Z lim sup 1 + ∥vt (x)∥2 + ∥ut (x)∥2 + ∥at (x)∥2 + ∥Σt (x)∥2F ρX t (x) dx = 0 . R→∞ t∈[0,T ] ∥x∥≥R
(36)
(R3) The continuity equation X ∂t ρX t + ∇ · (ut ρt ) = 0
admits a classical characteristic representation on [0, T ]. Namely, there exists a flow map Φt : Rd → Rd such that for all t ∈ [0, T ], d Φt (x) = ut (Φt (x)) , dt
Φ0 (x) = x ,
X (Φt )# ρX 0 = ρt .
Moreover, for all t ∈ (0, T ), Z lim
h→0 Rd
Z lim
h→0 Rd
2 Φt+h (x) − Φt (x) − ut (Φt (x)) ρX 0 (x) dx = 0 , h
2 Φt+h (x) + Φt−h (x) − 2Φt (x) − a (Φ (x)) ρX t t 0 (x) dx = 0 . h2
Proof. We prove property R1 in Section E.2.1. We prove property R2 in Section E.2.2. We prove property R3 in Section E.2.3. E.2.1
Proof of R1 ρX
Proof of Lemma 14 (property R1). By Assumption 1, q0 = ν 0X is continuously differentiable with ∥∇q0 ∥∞ := supx∈Rd ∥∇q0 (x)∥ < ∞ and ζ ≤ q0 (x) ≤ ξ for all x ∈ Rd . Recall from Lemma 13 that qt (x) =
ρX t (x) ∈ [ζ, ξ] , ν X (x)
gt (x, y) =
ρXY t (x, y) ∈ [ζ, ξ] , ν XY (x, y)
∀(x, y) ∈ R2d , t ∈ [0, T ] .
We define the notations: Z jt (x) :=
y gt (x, y)γ(y) dy , ZR
d
Kt (x) :=
yy ⊤ gt (x, y)γ(y) dy .
Rd
Using these, we may rewrite ut (x) = E[Yt | Xt = x] and Mt (x) = E[Yt Yt⊤ | Xt = x] as: 1 ut (x) = X ρt (x)
Z ν X (x) jt (x) XY yρt (x, y) dy = X y gt (x, y) γ(y) dy = , qt (x) ρt (x) Rd Rd
Z
43
(37a)
1 Mt (x) = X ρt (x)
Z ν X (x) Kt (x) ⊤ XY yy ρt (x, y) dy = X yy ⊤ gt (x, y) γ(y) dy = . qt (x) ρt (x) Rd Rd
Z
(37b)
X 1 d Since ν X ∈ C 1 (Rd ), in order to show that ρX t (x) = qt (x)ν (x) ∈ C (R × [0, ∞)), it suffices 1 d to show that qt (x) ∈ C (R × [0, ∞)). For ut and Mt , since 0 < ζ ≤ qt (x), it suffices to show jt (x), Kt (x) ∈ C 1 ([0, ∞) × Rd ]). We first consider the differentiability of gt (x, y), since qt (x), jt (x) and Kt (x) are defined using gt (x, y). Recall from the proof of Lemma 13 that gt (x, y) = q0 (ΨX −t (x, y)). By chain rule,
⊤ ∇x gt (x, y) = ∇x ΨX ∇q0 (ΨX −t (x, y) −t (x, y)) , X X Y ∂t gt (x, y) = ⟨∇q0 (ΨX −t (x, y)), ∂t Ψ−t (x, y)⟩ = −⟨∇q0 (Ψ−t (x, y)), Ψ−t (x, y)⟩ .
Thus, gt (x, y) is differentiable with respect to t and x, by continuous differentiability of ∇ΨX t (x, y) and q0 (x) with respect to t and x. In particular, since ∥∇q0 ∥∞ = supx∈Rd ∥∇q0 (x)∥ < ∞, X ∥∇x gt (x, y)∥ ≤ ∇x ΨX −t (x, y) op · ∥∇q0 (Ψ−t (x, y))∥ ≤ CT ∥∇q0 ∥∞ ,
|∂t gt (x, y)| ≤ ∥ΨY−t (x, y)∥ · ∥∇q0 (ΨX −t (x, y))∥ ≤ CT (1 + ∥x∥ + ∥y∥)∥∇q0 ∥∞ . Here, CT is the constant given by Lemma 7, which satisfies sup ∥∇Ψt (x, y)∥op ≤ CT ,
∥Ψt (x, y)∥ ≤ CT (1 + ∥x∥ + ∥y∥) .
|t|≤T
We now show that qt (x), jt (x) and Kt (x) are continuously differentiable with respect to t and x. Chain rule gives: Z Z ∇x qt (x) = ∇x gt (x, y) γ(y) dy , ∂t qt (x) = ∂t gt (x, y) γ(y) dy , d d ZR ZR ∇x jt (x) = y ∇x gt (x, y) γ(y) dy , ∂t jt (x) = y ∂t gt (x, y) γ(y) dy , d d R R Z Z ∇x Kt (x) = yy ⊤ ∇x gt (x, y) γ(y) dy , ∂t Kt (x) = yy ⊤ ∂t gt (x, y) γ(y) dy . Rd
Rd
Note that differentiation under the integration sign is valid, since ∥∇x gt (x, y)∥ and ∥∂t gt (x, y)∥ are upper bounded by polynomials of x and y that are integrable against the Gaussian distribution γ(y). Moreover, since ∇x gt (x, y) and ∂t gt (x, y) are continuous in (t, x) for each fixed y, dominated convergence theorem shows qt (x), jt (x) and Kt (x) are continuously differentiable in t and x. Remark 2. By Lemma 13, gt (x, y) ≤ ξ. Therefore, there exists a constant C0 dependent only on ξ and dimension d, such that |qt (x)| + ∥jt (x)∥ + ∥Kt (x)∥F ≤ C0 ,
∀t ∈ [0, T ], x ∈ Rd .
(38)
The bounds are uniform in t and x because qt , jt , Kt are the zeroth, first, and second Gaussian moments weighted by the uniformly bounded factor gt (x, ·), and the standard Gaussian distribution has finite moments of all orders: Z Z qt (x) = gt (x, y)γ(y) dy ≤ ξ γ(y) dy < ∞ , Rd
Rd
44
Z
Z
∥y∥γ(y) dy < ∞ , gt (x, y)∥y∥γ(y) dy ≤ ξ ∥jt (x)∥ ≤ d d R R Z Z ∥y∥2 γ(y) dy < ∞ . gt (x, y)∥y∥2 γ(y) dy ≤ ξ ∥Kt (x)∥F ≤ Rd
Rd
Since ζ ≤ qt (x) by Lemma 13, ∥ut (x)∥ and ∥Mt (x)∥F are also upper bounded: ∥ut (x)∥ + ∥Mt (x)∥F ≤
C0 , ζ
∀t ∈ [0, T ], x ∈ Rd .
The proof follows from the representation of ut and Mt in Eq. (37): ut (x) = E.2.2
jt (x) , qt (x)
Mt (x) =
Kt (x) . qt (x)
Proof of R2
Proof of Lemma 14 (property R2). We verify the limit in (36) is 0 term by term. Throughout the proof, CT denotes a finite constant that may change from line to line and may depend on T , ζ, ξ, and q0 , but not on R, t, or x. X By Lemma 13, we have ρX t (x) ≤ ξν (x) for all t ∈ [0, T ]. Therefore, Z Z ρX (x) dx ≤ ξ lim ν X (x) dx = 0 . lim sup t
First term of (36).
R→∞ ∥x∥≥R
R→∞ t∈[0,T ] ∥x∥≥R
Second term of (36). Recall vt (x) = x − Rt (x) is the displacement field, where Rt is the optimal X transport map that satisfies (Rt )# ρX t = ν . We can bound: Z Z Z 2 X 2 X ∥vt (x)∥ ρt (x) dx ≤ 2 ∥x∥ ρt (x) dx + 2 ∥Rt (x)∥2 ρX (39) t (x) dx . ∥x∥≥R
∥x∥≥R
∥x∥≥R
The first term in (39) is controlled by warmness: Z Z lim sup ∥x∥2 ρX (x) dx ≤ ξ lim t R→∞ t∈[0,T ] ∥x∥≥R
R→∞ ∥x∥≥R
∥x∥2 ν X (x) dx = 0 .
X For the second term in (39), let x ∼ ρX t , and let x̃ := Rt (x) ∼ ν . Given any B > 0, the following inequality holds for any x: ∥Rt (x)∥2 1{∥x∥≥R} (x) = ∥Rt (x)∥2 1{∥x∥≥R} (x) 1{∥Rt (x)∥<B} (x) + 1{∥Rt (x)∥≥B} (x)
≤ B 2 1{∥x∥≥R} (x) + ∥Rt (x)∥2 1{∥Rt (x)∥≥B} (x) . Taking expectation on both sides with respect to x ∼ ρX t gives: Z Z Z 2 X 2 X ∥Rt (x)∥ ρt (x) dx ≤ B ρt (x) dx + ∥Rt (x)∥2 ρX t (x) dx ∥x∥≥R ∥x∥≥R ∥Rt (x)∥≥B Z Z ≤ B2 ρX (x) dx + ∥x̃∥2 ν X (x) dx t ∥x∥≥R ∥x̃∥≥B Z Z ≤ B2ξ ν X (x) dx + ∥x̃∥2 ν X (x) , ∥x∥≥R
45
∥x̃∥≥B
where the last step follows from warmness. For fixed B, the first term vanishes upon taking R → ∞. Taking B → ∞, the second term vanishes since ν X ∈ P2 (Rd ). Therefore, Z lim sup ∥Rt (x)∥2 ρX t (x) dx = 0 . R→∞ t∈[0,T ] ∥x∥≥R
Third term of (36).
Recall from Remark 2 that ∥ut (x)∥ is uniformly upper bounded:
C0 , ζ here, C0 is a constant dependent only on ξ and dimension d. Therefore, Z Z C02 lim sup ρX (x) dx lim sup ∥ut (x)∥2 ρX (x) dx ≤ t R→∞ t∈[0,T ] ∥x∥≥R ζ 2 R→∞ t∈[0,T ] ∥x∥≥R t Z C02 ≤ 2 · ξ lim ν X (x) dx R→∞ ∥x∥≥R ζ ∥ut (x)∥ ≤
= 0. Fourth term of (36). Recall from (34) that at (x) = ∂t ut (x) + ∇ut (x)ut (x). Since ut = qjtt , we can calculate using chain rule and Cauchy-Schwartz inequality: ∥∇jt (x)∥ ∥jt (x)∥ · ∥∇qt (x)∥ + , qt (x) qt (x)2 ∥∂t jt (x)∥ ∥jt (x)∥ · |∂t qt (x)| ∥∂t ut (x)∥ ≤ + . qt (x) qt (x)2 R R Recall qt (x) = Rd gt (x, y) γ(y) dy and jt (x) = Rd y gt (x, y) γ(y) dy. Recall also from Remark 2 that for all t ∈ [0, T ], ∥∇x gt (x, y)∥ ≤ CT , |∂t gt (x, y)| ≤ CT (1 + ∥x∥ + ∥y∥), and ∥jt (x)∥ ≤ CT . Thus, differentiating under the Gaussian integral as justified in Section E.2.1, Z ∥∇x qt (x)∥ + ∥∇x jt (x)∥ ≤ CT (1 + ∥y∥)γ(y) dy ≤ CT , d ZR |∂t qt (x)| + ∥∂t jt (x)∥ ≤ CT (1 + ∥y∥)(1 + ∥x∥ + ∥y∥)γ(y) dy ≤ CT (1 + ∥x∥) . ∥∇ut (x)∥ ≤
Rd
Finally, since qt ≥ ζ, and jt is bounded, the quotient rule gives ∥∇jt (x)∥ ∥jt (x)∥ · ∥∇qt (x)∥ + ≤ CT , qt (x) qt (x)2 ∥∂t jt (x)∥ ∥jt (x)∥ · |∂t qt (x)| ∥∂t ut (x)∥ ≤ + ≤ CT (1 + ∥x∥) . qt (x) qt (x)2 ∥∇ut (x)∥ ≤
Together with the boundedness of ut , this implies ∥at (x)∥ ≤ ∥∂t ut (x)∥ + ∥∇ut (x)∥ · ∥ut (x)∥ ≤ CT (1 + ∥x∥) . Therefore, Z lim sup
R→∞ t∈[0,T ] ∥x∥≥R
∥at (x)∥2 ρX t (x) dx ≤ CT
Z lim sup
1 + ∥x∥2 ρX t (x) dx
Z
1 + ∥x∥2 ν X (x) dx = 0 ,
R→∞ t∈[0,T ] ∥x∥≥R
≤ CT ξ lim
R→∞ ∥x∥≥R
where the last equality holds because ν X ∈ P2 (Rd ). 46
Fifth term of (36).
Recall that Σt (x) = Mt (x) − ut (x)ut (x)⊤ . Cauchy-Schwartz gives: ∥Σt (x)∥F ≤ ∥Mt (x)∥F + ∥ut (x)ut (x)⊤ ∥F .
Recall from Remark 2 that ∥ut (x)∥ and ∥Mt (x)∥F are uniformly upper bounded: ∥ut (x)∥ ≤
C0 , ζ
∥Mt (x)∥F ≤
C0 . ζ
here, C0 is a constant dependent only on ξ and dimension d. Then C0 C0 ⊤ 1+ . ∥Σt (x)∥F ≤ ∥Mt (x)∥F + ∥ut (x)ut (x) ∥F ≤ ζ ζ Therefore, Z lim sup
R→∞ t∈[0,T ] ∥x∥≥R
2 Z C0 2 1+ ρX lim sup t (x) dx R→∞ t∈[0,T ] ∥x∥≥R ζ 2 Z C0 C0 2 1+ ξ lim ≤ ν X (x) dx = 0 . R→∞ ∥x∥≥R ζ ζ
∥Σt (x)∥2F ρX t (x) dx ≤
C0 ζ
Combining the five estimates above gives the claim in R2: Z 1 + ∥vt (x)∥2 + ∥ut (x)∥2 + ∥at (x)∥2 + ∥Σt (x)∥2F ρX lim sup t (x) dx = 0 . R→∞ t∈[0,T ] ∥x∥≥R
E.2.3
Proof of R3
Proof of Lemma 14 (property R3). By R1, as shown in Section E.2.2, we have that (t, x) 7→ ut (x) is continuously differentiable in t and x. Moreover by Remark 2, ut (x) is uniformly bounded on (t, x) ∈ [0, T ] × Rd . Hence, the ordinary differential equation d Φt (x) = ut (Φt (x)) , dt
Φ0 (x) = x ,
admits a solution Φt (x) for all t ∈ [0, T ]. By Lemma 12, ρX t solves the continuity equation X ∂t ρX t + ∇ · (ut ρt ) = 0 .
Since the velocity field ut is bounded and continuously differentiable, the classical method of charX acteristics applies (see e.g. [AGS05, Theorem 8.1.8]), yielding (Φt )# ρX 0 = ρt . We now prove the differentiability results. Fix t ∈ (0, T ) and take h such that t − h ∈ [0, T ] and t + h ∈ [0, T ]. Assume without loss of generality that h > 0. For the first-order differentiability, we calculate the following: Z Rd
Z 2 1 h ut+s (Φt+s (x)) ds − ut (Φt (x)) ρX 0 (x) dx h d R 0 Z Z h (2) 1 ≤ ∥ut+s (Φt+s (x)) − ut (Φt (x))∥2 ρX 0 (x) ds dx h Rd 0
2 Φt+h (x) − Φt (x) (1) − ut (Φt (x)) ρX 0 (x) dx = h
Z
47
(3) 1
=
h
Z hZ 0
Rd
∥ut+s (Φt+s (x)) − ut (Φt (x))∥2 ρX 0 (x) dx ds ,
where (1) follows from fundamental theorem of calculus, (2) follows from Jensen’s inequality, (3) follows from switching the order of integration, which is valid by the boundedness of ut (x). Therefore, R it suffices to show lims→0 Rd ∥ut+s (Φt+s (x)) − ut (Φt (x))∥2 ρX 0 (x) dx = 0. We can calculate: Z Z t+s 2 (1) 2 X lim ∥ut+s (Φt+s (x)) − ut (Φt (x))∥ ρ0 (x) dx = lim ar (Φr (x))dr ρX 0 (x) dx s→0 Rd s→0 Rd t Z Z t+s (2) ≤ lim ∥ar (Φr (x))∥2 ρX 0 (x) dr dx s→0 Rd t Z t+s Z (3) = lim ∥ar (x)∥2 ρX r (x) dx dr , s→0 t d R Z
where (1) follows from fundamental theorem of calculus, (2) follows from Jensen’s inequality, and (3) follows from switching the order of integration and applying a change of variable X ∼ ρX 0 , so Φr (X) ∼ ρX . By the derivations in Section E.2.2, ∥a (x)∥ is uniformly bounded by C (1 + ∥x∥) t T r d X for (t, x) ∈ [0, T ] × R . By the derivations in Section E.2.1, the second moment of ρt is uniformly R t+s R 2 X bounded for t ∈ [0, T ]. Thus, lims→0 t Rd ∥ar (x)∥ ρr (x) dx dr = 0. Therefore, the first-order differentiability holds. We now prove the second-order differentiability. We first calculate: 2 Φt+h (x) − 2Φt (x) + Φt−h (x) − a (Φ (x)) ρX t t 0 (x) dx 2 h d R Z Z h 2 1 (1) = [ut+s (Φt+s (x)) − ut−s (Φt−s (x))] ds − at (Φt (x)) ρX 0 (x) dx 2 Rd h 0 Z Z h 2 1 (2) = (h − s) [a (Φ (x)) + a (Φ (x))] ds − a (Φ (x)) ρX t+s t+s t−s t−s t t 0 (x) dx 2 h d R 0 Z h Z 2 1 (3) (h − s) [at+s (Φt+s (x)) − 2at (Φt (x)) + at−s (Φt−s (x))] ds ρX = 0 (x) dx 2 h d 0 R Z Z h 2 (4) 1 ≤2 (h − s) [at−s (Φt−s (x)) − at (Φt (x))] ds ρX 0 (x) dx 2 Rd h 0 Z Z h 2 1 +2 (h − s) [a (Φ (x)) − a (Φ (x))] ds ρX t+s t+s t t 0 (x) dx 2 h d R 0 Z h Z (5) 1 2(h − s) ≤ ∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx ds 2 0 h2 d R Z Z 1 h 2(h − s) + ∥at−s (Φt−s (x)) − at (Φt (x))∥2 ρX 0 (x) dx ds , 2 0 h2 d R
Z
where in (1) we used the fact that Φ̇t (x) = ut (Φt (x)); in (2) we used the fact that u̇t (Φt (x)) = at (Φt (x)); in (3) we rearranged the integrand; in (4) we used Cauchy-Schwartz inequality and in (5) we used Jensen’s inequality. Since 0 < s ≤ h, to show Z lim h↓0
Rd
2 Φt+h (x) − 2Φt (x) + Φt−h (x) − a (Φ (x)) ρX t t 0 (x) dx = 0 , h2
48
it suffices to show that Z lim
s→0 Rd
∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx = 0 .
We prove this by decomposing the domain of integration into a compact set and its tail integral. We X will use the following tail estimate, which follows by R2 and the transport identity (Φr )# ρX 0 = ρr : Z lim sup ∥ar (Φr (x))∥2 ρX (40) 0 (x) dx = 0 . R→∞ r∈[0,T ] {∥Φr (x)∥≥R}
Fix R > 0. We split Z
∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx
RdZ
= {∥Φt (x)∥≤R}
∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx +
Z {∥Φt (x)∥>R}
∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx .
We first treat the first term (compact part) above. Since ut (x) is uniformly bounded on (t, x) ∈ [0, T ] × Rd (see the calculation for the third term in Section E.2.2), we have Z t+s ∥Φt+s (x) − Φt (x)∥ ≤ ∥ur (Φr (x))∥ dr ≤ CT |s| . t
Therefore, Φt (x) is continuous with respect to t. Recall the definition of acceleration field: at (x) = ∂t ut (x) + ∇ut (x) ut (x). We have shown in R1 that ut (x) is continuously differentiable with respect to t and x, and hence at (x) is continuous with respect to t and x. Since at (x) is continuous in t and x and Φt (x) is continuous in t, their composition at (Φt (x)) is continuous in t, for any x. Therefore, Z lim ∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx = 0 . s→0 {∥Φt (x)∥≤R}
We now treat the second term (tail part) above. We have Z ∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx {∥Φt (x)∥>R} Z Z 2 X ≤2 ∥at (Φt (x))∥ ρ0 (x) dx + 2 ∥at+s (Φt+s (x))∥2 ρX 0 (x) dx . {∥Φt (x)∥>R}
{∥Φt (x)∥>R}
For the first term above, we use the trivial bound Z Z ∥at (Φt (x))∥2 ρX (x) dx ≤ sup 0 {∥Φt (x)∥>R}
r∈[0,T ] {∥Φr (x)∥>R/2}
∥ar (Φr (x))∥2 ρX 0 (x) dx
which we will control below. For the second term, using the estimate ∥Φt+s (x) − Φt (x)∥ ≤ CT |s|, we have that ∥Φt (x)∥ > R implies ∥Φt+s (x)∥ > R − CT |s|. Thus, for |s| sufficiently small such that R − CT |s| ≥ R/2, Z Z 2 X ∥at+s (Φt+s (x))∥ ρ0 (x) dx ≤ ∥at+s (Φt+s (x))∥2 ρX 0 (x) dx {∥Φt (x)∥>R} {∥Φt+s (x)∥>R/2} Z ≤ sup ∥ar (Φr (x))∥2 ρX 0 (x) dx . r∈[0,T ] {∥Φr (x)∥>R/2}
49
Therefore, Z {∥Φt (x)∥>R}
∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx ≤ 4 sup
Z
r∈[0,T ] {∥Φr (x)∥>R/2}
Combining the two estimates above gives Z Z 2 X lim ∥at+s (Φt+s (x)) − at (Φt (x))∥ ρ0 (x) dx ≤ 4 sup s→0 Rd
r∈[0,T ] {∥Φr (x)∥>R/2}
∥ar (Φr (x))∥2 ρX 0 (x) dx .
∥ar (Φr (x))∥2 ρX 0 (x) dx .
Recall R > 0 is arbitrary, and it does not appear on the left-hand side above. Sending R → ∞ and using the uniform tail estimate (40) yields Z ∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX lim 0 (x) dx = 0 . s→0 Rd
Therefore, 2 Φt+h (x) + Φt−h (x) − 2Φt (x) − a (Φ (x)) ρX t t 0 (x) dx h→0 Rd h2 Z h Z h − |s| ≤ lim ∥at+s (Φt+s (x)) − at (Φt (x))∥2 ρX 0 (x) dx ds h→0 −h h2 d R
Z
lim
= 0, where the last step holds because the triangular kernel h−|s| 1{|s|≤h} has total mass 1, and the inner h2 integral converges to 0 as s → 0.
F
Bound on average KL divergence under regular initialization
In this section, we prove in Lemma 22 the bound on the average KL divergence along the Hamiltonian flow claimed in Lemma 3, under the warmness and smoothness assumption (Assumption 1). We first provide some preliminary results that we will use in the calculation. In Lemma 15 in Section F.1, we provide a formula for the acceleration field at . In Lemma 16 in Section F.2, we provide a formula for the second time derivative of W22 . In Lemma 18 in Section F.3.1, we provide a bound to the second time derivative of W22 in terms of the difference between KL divergence at the initial and final times along the Hamiltonian flow. This provides a second-order differential inequality relating Wasserstein distance and KL divergence, that we can integrate twice to obtain the claimed bound on the average KL divergence along Hamiltonian flow; see Lemma 22 in Section F.4. We recall the definitions of the conditional mean ut (31), conditional second moment Mt (32), conditional covariance Σt (33), acceleration field at (34), and displacement field vt (35).
F.1
Formula for the acceleration field
Lemma 15. Assume Assumption 1, in addition to the set up in Lemma 3. Then the following equation holds for all x ∈ Rd : at (x) = −∇f (x) −
1 ρX t (x) 50
∇ · (ρX t Σt )(x) .
Proof. By R1, the acceleration field at (x) is well-defined. We first show the equation holds in the sense of distributions, using the identities shown in Lemma 12. First, we rewrite the continuity equation using Mt (x) = ut (x)ut (x)⊤ + Σt (x): X ⊤ X X 0 = ∂t (ρX t ut ) + ∇ · (ρt ut ut ) + ∇ · (ρt Σt ) + ρt ∇f (1)
X X X X X = ρX t ∂t ut + ut ∂t ρt + ρt ∇ut ut + ut ∇ · (ρt ut ) + ∇ · (ρt Σt ) + ρt ∇f X X X X = ρX t ∂t ut + ∇ut ut + ut ∂t ρt + ∇ · (ρt ut ) + ∇ · (ρt Σt ) + ρt ∇f (2) X X = ρt ∂t ut + ∇ut ut + ∇ · (ρX t Σt ) + ρt ∇f (3)
X X = ρX t at + ∇ · (ρt Σt ) + ρt ∇f,
where (1) holds because ut (x) is continuously differentiable by R1, (2) holds by the continuity X equation ∂t ρX t + ∇ · (ρt ut ) = 0, and (3) holds by the definition of at = ∂t ut + ∇ut ut . Therefore, we have the distributional identity X X 0 = ρX t at + ∇ · (ρt Σt ) + ρt ∇f,
(41)
which means the identity above holds when we integrate both sides against test functions in space and time. We claim that Eq. (41) in fact holds pointwise for all (t, x). Indeed, by R1, all terms appearing in Eq. (41) are continuous in (t, x). If the expression in Eq. (41) were nonzero at some point (t0 , x0 ), then one of its components would be nonzero there. By continuity, that component would have a fixed sign on a small space-time neighborhood of (t0 , x0 ). Therefore, integrating (41) against a nonnegative smooth bump function supported in this neighborhood would give a strictly nonzero integral, contradicting the distributional identity. Thus for every t ∈ (0, T ) and x ∈ Rd , X X 0 = ρX t (x)at (x) + ∇ · (ρt Σt )(x) + ρt (x)∇f (x). X for all x ∈ Rd . Since ν X ∈ P d X X By Lemma 13, ρX 2,ac,fs (R ), ρt (x) > ζν (x) > 0 holds for t (x) ≥ ζν d X all x ∈ R . Dividing by ρt (x) from both sides gives:
1 ∇ · (ρX t Σt )(x) + ∇f (x), ρX (x) t 1 at (x) = −∇f (x) − X ∇ · (ρX t Σt )(x). ρt (x) 0 = at (x) +
F.2
Formula for the upper second time derivative of Wasserstein distance
The following estimate is in the same spirit as [CCT18, Theorem 1, 2], where it is formulated for a.e. t ∈ (0, T ). Our assumptions include a pointwise-in-time L2 differentiability condition in R3, and therefore the statement holds for each t ∈ (0, T ). Lemma 16. Assume Assumption 1, in addition to the set up in Lemma 3. Then the map t 7→ 1 2 X X 2 W2 (ρt , ν ) is differentiable for every t ∈ (0, T ), with one-sided derivatives at the endpoints, and Z d 1 2 X X W2 (ρt , ν ) = ⟨vt (x), ut (x)⟩ρX t (x) dx . dt 2 d R 51
Moreover, for every t ∈ (0, T ), Z 1 2 X 1 2 X X X 2 X X 2 W2 (ρt+h , ν ) + 2 W2 (ρt−h , ν ) − W2 (ρt , ν ) ≤ lim sup h2 h↓0 Rd
∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx .
d Proof. Note ν X ∈ P2,ac,fs (Rd ) by assumption, and by Lemma 11, we also have ρX t ∈ P2,ac,fs (R ). X Therefore, W22 (ρX t , ν ) < ∞ for all t ≥ 0. X X First derivative: We first show t 7→ W22 (ρX t , ν ) is differentiable. By R3, ρt satisfies the X continuity equation ∂t ρX t + ∇ · (ρt ut ) = 0. Moreover, recall from Remark 2 that ∥ut (x)∥ ≤ C0 for d all (t, x) ∈ [0, T ] × R , and C0 is a constant dependent only on ξ and dimension d. Therefore by X [Vil09, Theorem 23.9, Proof Step 2], t 7→ 21 W22 (ρX t , ν ) is differentiable for every t ∈ (0, T ), with one-sided derivatives at the endpoints. In particular, the derivative is equal to the following: Z d 1 2 X X ⟨vt (x), ut (x)⟩ ρX W (ρ , ν ) = t (x) dx. dt 2 2 t d R
Second derivative: We now show a bound for the upper second time derivative. Recall Φt X X X pushes ρX 0 forward to ρt for all t ∈ (0, T ), so ρt = (Φt )# ρ0 . Fix t ∈ (0, T ), and let h > 0 be small enough such that 0 < t − h < t + h < T . We construct (suboptimal) couplings between ρX t+h and X ), as follows: and ν ν X ), and between ρX t−h + := (Φt+h , Rt ◦ Φt )# ρX γt,h 0 , − := (Φt−h , Rt ◦ Φt )# ρX γt,h 0 . X X X X These are valid couplings, since (Φt±h )# ρX 0 = ρt±h , and (Rt ◦ Φt )# ρ0 = (Rt )# ρt = ν . Then by the definition of the Wasserstein distance: Z Z 2 X X 2 ± W2 (ρt±h , ν ) ≤ ∥x − y∥ γt,h (x, y) dx dy = ∥Φt±h (x) − Rt (Φt (x))∥2 ρX 0 (x) dx . R2d
Rd
X X X Since (Φt )# ρX 0 = ρt and Rt is the optimal transport map from ρt to ν , we know that: Z Z 2 X X 2 X W2 (ρt , ν ) = ∥y − Rt (y)∥ ρt (y) dy = ∥Φt (x) − Rt (Φt (x))∥2 ρX 0 (x) dx . Rd
Rd
Then we can compute, recalling that vt (x) = x − Rt (x): X 2 X X 2 X X W22 (ρX t+h , ν ) + W2 (ρt−h , ν ) − 2W2 (ρt , ν ) Z 2 2 2 ≤ ∥Φt+h (x) − Rt (Φt (x))∥ + ∥Φt−h (x) − Rt (Φt (x))∥ − 2∥Φt (x) − Rt (Φt (x))∥ ρX 0 (x) dx d ZR = ∥Φt+h (x) − Φt (x) + vt (Φt (x))∥2 + ∥Φt−h (x) − Φt (x) + vt (Φt (x))∥2 − 2∥vt (Φt (x))∥2 ρX 0 (x) dx d ZR (x) dx . = ∥Φt+h (x) − Φt (x)∥2 + ∥Φt−h (x) − Φt (x)∥2 + 2⟨Φt+h (x) + Φt−h (x) − 2Φt (x), vt (Φt (x))⟩ ρX {z } | {z } 0 Rd | :=Ah (x)
:=Bh (x)
By R3, we have Z Z Z 1 2 X lim 2 Ah (x)ρX (x) dx = 2 ∥u (Φ (x))∥ ρ (x) dx = 2 ∥ut (x)∥2 ρX t t t (x) dx . 0 0 h→0 h d d d R R R 52
X ∈ P d 2 X X On the other hand, since ρX 2,ac (R ), W2 (ρt , ν ) < ∞, and therefore, vt (Φt ) is squaret ,ν integrable under ρX 0 . By R1, at (x) is well-defined. Combining this with R3, we obtain Z Z Z 1 X X ⟨at (x), vt (x)⟩ρX ⟨at (Φt (x)), vt (Φt (x))⟩ρ0 (x) dx = 2 Bh (x)ρ0 (x) dx = 2 lim t (x) dx . h→0 h2 Rd Rd Rd
Combining the two steps above gives: lim sup h→0
Z X 2 X X 2 X X W22 (ρX 1 t+h , ν ) + W2 (ρt−h , ν ) − 2W2 (ρt , ν ) ≤ lim (Ah (x) + Bh (x)) ρX 0 (x) dx h→0 h2 Rd h2 Z ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX =2 t (x) dx . Rd
Dividing by 2 on both sides gives the desired inequality.
F.3
A bound on the second time derivative of the Wasserstein distance
In this section, we provide a bound on the second time derivative of the Wasserstein X distance, i.e., R 2 the right-hand side of the identity in Lemma 16: Rd ∥ut (x)∥ + ⟨vt (x), at (x)⟩ ρt (x) dx. Similar bounds were previously derived in [Lu26, Lemma 4.1], under slightly different regularity assumptions. To give an intuition for how to bound this quantity, we first consider the following heuristic derivation. First, using the representation of at ρX t from Lemma 15: X X ρX t (x)at (x) = −ρt (x)∇f (x) − ∇ · (ρt Σt )(x) .
Therefore, Z Rd
⟨vt (x), ρX t (x)at (x)⟩ dx = −
Z Rd
X ⟨vt (x), ρX t (x)∇f (x) + ∇ · (ρt Σt )(x)⟩ dx .
To simplify the second term above, we want to apply integration by parts to obtain: Z Z Z X X − ⟨vt (x), ∇·(ρt Σt )(x)⟩ dx = ⟨∇vt (x), Σt (x)⟩F ρt (x) dx = ⟨Id −∇Rt (x), Σt (x)⟩F ρX t (x) dx . Rd
Rd
Rd
The proof then proceeds by bounding the right-hand side above However, the heuristic computation above is not rigorous due to several issues. First, vt is not compactly supported, so an integration by parts on Rd may produce boundary terms at infinity. Second, by Brenier’s theorem, Rt is characterized ρX t -a.e. as Rt = ∇φt for a convex function φt . In general, the Hessian of a convex function φ is not a classical map, but may contain singularities (e.g., φ(x) = |x| in one dimension has second derivative which is singular at x = 0). Thus, ∇Rt is not an ordinary matrix-valued function defined everywhere, so the formal identity ∇vt = Id − ∇Rt does not hold everywhere; see e.g. [Vil03, Remark 4.11] and [Vil09, pp. 273–274] for further discussions. To make the proof rigorous, we follow the treatment of [Lu26] using the notion of the distributional Hessian of convex functions. We recall the standard characterization from [EG15, Theorem 6.8]. Lemma 17. [EG15, Theorem 6.8] Let φ : Rd → R be convex. Then there exists a signed Radon measure µ = (µij )di,j=1 satisfying µij = µji such that for every smooth compactly supported test function ψ : Rd → R, Z Z ψ(x) µij ( dx) .
φ(x) ∂ij ψ(x) dx = Rd
Rd
53
The matrix-valued Radon measure [D2 φ] := (µij )di,j=1 is called the distributional Hessian of φ, and P it is positive semidefinite in sense that for every ξ ∈ Rd , the scalar Radon measure di,j=1 ξi ξj µij is nonnegative. Finally, the gradient of φ is locally of bounded variation, and therefore its distributional derivative is the distributional Hessian [D2 φ]. In our setting, we denote the distributional Hessian of φt , or equivalently the distributional derivative of the Brenier map Rt = ∇φt , by d∇Rt := [D2 φ]. By the Lebesgue decomposition theorem, each component of this matrix-valued measure decomposes into an absolutely continuous part and a singular part with respect to Lebesgue measure. Thus, the matrix-valued measure itself decomposes as d∇Rt (x) = Gt (x) dx + d∇s Rt (x) = Gt (x) dx + Nt (x)σt ( dx) ,
(42)
where the absolutely continuous part admits an integrable density Gt : Rd → Rd×d , and the singular part is d∇s Rt (x) = Nt (x)σt ( dx), where Nt (x) is a matrix-valued function and σt is singular with respect to the Lebesgue measure. Since the distributional Hessian of a convex function is positive semidefinite as a matrix-valued measure, both its absolutely continuous and singular parts are positive semidefinite. Therefore, Gt (x) ⪰ 0 for Lebesgue-a.e. x, and Nt (x) ⪰ 0 for σt -a.e. x. In particular, Gt satisfies the following Monge-Ampere equation for ρX t -a.e. x: X ρX t (x) = ν (Rt (x)) det Gt (x) ,
(43)
see [Vil09, Theorem 11.1, Example 11.2] for a characterization. To handle the aforementioned issue of non-vanishing boundary terms, we use a standard cutoff argument. Let η : R → [0, ∞) be a smooth, non-increasing, compactly supported function satisfying: 0 ≤ η ≤ 1,
η(s) = 1 for 0 ≤ s ≤ 1,
η(s) = 0 for s ≥ 4 .
See e.g. [Eva10, Appendix C.5] for the construction of a cutoff satisfying the required properties above. For R > 0, define ∥x∥2 := . (44) χR (x) η R2 Then we have: 0 ≤ χR ≤ 1,
χR (x) = 1 for ∥x∥ ≤ R,
χR (x) = 0 for ∥x∥ ≥ 2R .
Moreover, by chain rule, ∇χR (x) = η ′
∥x∥2 R2
2x . R2
Since η ′ is compactly supported on [1, 4], ∇χR is supported on the annulus AR := {x ∈ Rd : R ≤ ∥x∥ ≤ 2R} . d Hence, for some constant C > 0 independent of R, we have ∥∇χR (x)∥ ≤ C R for all x ∈ R .
54
F.3.1
Bound on the second derivative of the Wasserstein distance
With the preparation above, we now provide a bound on the second derivative of the Wasserstein distance. Lemma 18. Assume Assumption 1 and the set up in Lemma 3. In particular, assume ν X ∝ exp(−f ) is M -semi-log-concave. Then for any T ∈ (0, ∞), the following holds for all t ∈ (0, T ): Z M 2 X X X X X X W (ρ , ν ) . ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx ≤ 2KL(ρ0 ∥ ν ) − 3KL(ρt ∥ ν ) + 2 2 t Rd Proof. Let χR be the cutoff function defined in (44). We perform the following calculation: Z ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx Rd Z (1) χR (x) ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX = lim t (x) dx R→∞ Rd Z Z X (2) 2 X = lim χR (x) ∥ut (x)∥ − ⟨vt (x), ∇f (x)⟩ ρt (x) dx − vt (x), ∇ · (χR ρt Σt )(x) dx R→∞
(3)
Rd
Rd
Z
χR (x) ∥ut (x)∥2 + ⟨Id − Gt (x), Σt (x)⟩F − ⟨vt (x), ∇f (x)⟩ ρX t (x) dx R→∞ Rd Z (4) ρX M Y |X=x t (x) 2 ≤ lim χR (x) 2KL ρt ∥ γ − log X + ∥vt (x)∥ ρX t (x) dx R→∞ Rd ν (x) 2 Z M ρX (5) Y |X=x t (x) 2 = 2KL ρt + ∥vt (x)∥ ρX ∥ γ − log X t (x) dx ν (x) 2 Rd M 2 X X (6) X X X = 2KL(ρX W (ρ , ν ) . 0 ∥ ν ) − 3KL(ρt ∥ ν ) + 2 2 t ≤ lim
In the derivation above, (1) and (2) follow from Lemma 19, proved in Section F.3.2 below; (3) follows from Lemma 20, proved in Section F.3.3 below; (4) follows from Lemma 21, proved in Section F.3.4 below; (5) is justified below; and (6) follows from the chain rule for KL divergence, the conservation of joint KL divergence (Lemma 1), and the initialization ρXY = ρX 0 0 ⊗ γ: Z Y |X=x XY X KL ρt ∥ γ ρX ∥ ν XY − KL ρX t (x) dx = KL ρt t ∥ν Rd X = KL ρXY ∥ ν XY − KL ρX 0 t ∥ν X X = KL ρX − KL ρX . 0 ∥ν t ∥ν We now justify the convergence in (5). First, since 0 ≤ χR ≤ 1 and limR→∞ χR (x) = 1 for all x ∈ Rd , by the monotone convergence theorem, Z Z Y |X=x Y |X=x lim χR (x) 2KL ρt ∥ γ ρX (x) dx = 2KL ρ ∥ γ ρX t t (x) dx t R→∞ Rd
Rd
X X X = 2KL(ρX 0 ∥ ν ) − 2KL(ρt ∥ ν ) .
For the marginal term, note that X X Z Z ρX ρt (x) ρt (x) t (x) X ρt (x) dx = χR (x) X log ν X (x) dx . χR (x) log X (x) X (x) ν ν (x) ν d d R R 55
ρX
ρX
X −1 for s > 0, the positive and negative parts of t log t Since KL(ρX t ∥ ν ) < ∞, and s log s ≥ −e νX νX X are integrable under ν . Applying the dominated convergence theorem to the positive and negative parts yields Z Z ρX (x) X ρX (x) X X χR (x) log tX lim − log tX ρt (x) dx = − ρt (x) dx = −KL(ρX t ∥ ν ). R→∞ ν (x) ν (x) d d R R
Finally, recall that vt (x) = x − Rt (x) is the displacement map for the optimal coupling between ρX t R 2 X X 2 and ν X , and Rd ∥vt (x)∥2 ρX t (x) dx = W2 (ρt , ν ) < ∞. The convergence of the ∥vt (x)∥ term then also follows from monotone convergence theorem, since χR (x) is non-decreasing with respect to R for any fixed x. F.3.2
Helper result 1: Cutoff approximation
Lemma 19. Assume the setting of Lemma 18. We have: Z ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx Rd Z (1) = lim χR (x) ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx R→∞ Rd Z Z X (2) 2 = lim χR (x) ∥ut (x)∥ − ⟨vt (x), ∇f (x)⟩ ρt (x) dx − R→∞
Rd
Rd
vt (x), ∇ · (χR ρX t Σt )(x)
dx .
Proof. We first verify the convergence of the cutoff integral: Z (1 − χR (x)) ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx Rd
(1)
Z
≤ ∥x∥≥R (2)
Z
≤ ∥x∥≥R
∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx
1 1 2 2 ∥ut ∥ + ∥vt (x)∥ + ∥at (x)∥ ρX t (x) dx , 2 2 2
where (1) follows because 1 − χR (x) = 0 for ∥x∥ ≤ R and 0 ≤ 1 − χR (x) ≤ 1 for ∥x∥ ≥ R, according to the definition of χR ; and (2) follows from Young’s inequality. According to the property R2, the right-side of (2) above converges to 0 as R → ∞. Therefore, Z lim (1 − χR (x)) ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx = 0 . R→∞
Rd
This proves the first step (1): Z Z X 2 ∥ut (x)∥ + ⟨vt (x), at (x)⟩ ρt (x) dx = lim
R→∞ Rd
Rd
χR (x) ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx .
Next, using the representation of at from Lemma 15, we have: Z χR (x) ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx RdZ Z X 2 = χR (x) ∥ut (x)∥ − ⟨vt (x), ∇f (x)⟩ ρt (x) dx − χR (x) vt (x), ∇ · (ρX t Σt )(x) dx . Rd
Rd
56
We handle the last term above. By chain rule, X X χR ∇ · (ρX t Σt ) = ∇ · (χR ρt Σt ) − ρt Σt ∇χR .
Therefore, we have: Z χR (x) vt (x), ∇ · (ρX − t Σt )(x) dx d RZ Z X vt (x), ∇ · (χR ρt Σt )(x) dx + =−
Rd
Rd
vt (x), ρX t (x)Σt (x)∇χR (x) dx .
We claim the second term above is the boundary error that vanishes as R → ∞. Indeed, since ∇χR is supported on the annulus AR = {x ∈ Rd : ∥x∥ ∈ [R, 2R]} and ∥∇χR (x)∥ ≤ C/R, Cauchy– Schwarz inequality gives Z vt (x), ρX t (x)Σt (x)∇χR (x) dx d R Z C ≤ ∥vt (x)∥ · ∥Σt (x)∥op ρX t (x) dx R AR 1/2 Z 1/2 Z C 2 X 2 X ≤ ∥vt (x)∥ ρt (x) dx · ∥Σt (x)∥op ρt (x) dx R AR AR !1/2 !1/2 Z Z C 2 X 2 X ≤ ∥vt (x)∥ ρt (x) dx · ∥Σt (x)∥op ρt (x) dx . R ∥x∥≥R ∥x∥≥R By property R2, the right-hand side tends to 0 as R → ∞. Combining the arguments above, we conclude that step (2) is also valid, and thus: Z ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx Rd Z (1) = lim χR (x) ∥ut (x)∥2 + ⟨vt (x), at (x)⟩ ρX t (x) dx R→∞ Rd Z Z X (2) 2 X = lim χR (x) ∥ut (x)∥ − ⟨vt (x), ∇f (x)⟩ ρt (x) dx − vt (x), ∇ · (χR ρt Σt )(x) dx R→∞
Rd
Rd
as desired. F.3.3
Helper result 2: Integration by parts
We recall the measure decomposition in (42). Lemma 20. Assume the setting of Lemma 18. For every t ∈ (0, T ), we have: Z Z X 2 X lim χR (x) ∥ut (x)∥ − ⟨vt (x), ∇f (x)⟩ ρt (x) dx − vt (x), ∇ · (χR ρt Σt )(x) dx R→∞ Rd Rd Z ≤ lim χR (x) ∥ut (x)∥2 + ⟨Id − Gt (x), Σt (x)⟩F − ⟨vt (x), ∇f (x)⟩ ρX t (x) dx . R→∞ Rd
Proof. Since vt (x) = x − Rt (x) = x − ∇φt (x), by Lemma 17 it is locally of bounded variation, with distributional derivative characterized by the following: d∇vt (x) = (Id − Gt (x)) dx − Nt (x)σt ( dx) . 57
Since χR (x) ρX t (x) Σt (x) is compactly supported and is continuously differentiable, we may use the integration by parts formula for functions of bounded variation locally [EG15, Theorem 5.1]: Z Z X − vt (x), ∇ · (χR ρt Σt )(x) dx = d∇vt (x), χR (x) ρX t (x) Σt (x) F . Rd
Rd
Using the measure decomposition as discussed in (42), this further equals to: Z d∇vt (x), χR (x) ρX t (x) Σt (x) F d R Z Z χR (x)ρX = χR (x)ρX (x) ⟨Σ (x), I − G (x)⟩ dx − t t d t (x) ⟨Σt (x), Nt (x)⟩F σt ( dx) . t F Rd
Rd
Since Σt (x) is a conditional covariance matrix, it satisfies Σt (x) ⪰ 0. Since we have also shown that Nt (x) ⪰ 0, we have ⟨Σt (x), Nt (x)⟩F ≥ 0, and hence the singular term above is nonpositive. Therefore, Z Z X − vt (x), ∇ · (χR ρt Σt )(x) dx ≤ χR (x)ρX t (x) ⟨Σt (x), Id − Gt (x)⟩F dx . Rd
Rd
Since this holds for every R > 0, we can set R → ∞ to conclude: Z Z X 2 X lim χR (x) ∥ut (x)∥ − ⟨vt (x), ∇f (x)⟩ ρt (x) dx − vt (x), ∇ · (χR ρt Σt )(x) dx R→∞ Rd Rd Z ≤ lim χR (x) ∥ut (x)∥2 + ⟨Id − Gt (x), Σt (x)⟩F − ⟨vt (x), ∇f (x)⟩ ρX t (x) dx , R→∞ Rd
as desired. F.3.4
Helper result 3: Pointwise bound
We recall the measure decomposition in (42). Lemma 21. Assume the setting of Lemma 18. For all t ∈ (0, T ) and for ρX t -a.e. x, we have: ρX (x) M Y |X=x ∥ut (x)∥2 + ⟨Id − Gt (x), Σt (x)⟩F − ⟨∇f (x), vt (x)⟩ ≤ 2KL ρt ∥ γ − log tX + ∥vt (x)∥2 . ν (x) 2 Proof. We upper bound each term above. All statements in this proof hold for ρX t -a.e. x. First term.
We will show the following: Y |X=x ∥ut (x)∥2 ≤ 2KL ρt ∥ γ − Tr(Σt (x)) + log det Σt (x) + d .
Indeed, we have: KL
Y |X=x ρt ∥γ
Z
Y |X=x
ρt
(y) = log γ(y) Rd Z Y |X=x = −Ent ρt +
!
Rd
58
Y |X=x
ρt
(y) dy
d 1 log(2π) + ∥y∥2 2 2
Y |X=x
ρt
(y) dy
1 d Y |X=x = −Ent ρt + Tr (Mt (x)) + log(2π) 2 2 (1) 1 ≥ (Tr (Mt (x)) − log det Σt (x) − d) 2 1 ∥ut (x)∥2 + Tr (Σt (x)) − log det Σt (x) − d . = 2 In the above, (1) follows from the fact that Gaussian distribution maximizes entropy given fixed Y |X=x covariance, so Ent ρt ≤ Ent (N (0, Σt (x))) = d2 log(2πe) + 21 log det Σt (x). Rearranging gives the desired inequality. Second term.
We will show the following:
⟨Id − Gt (x), Σt (x)⟩F ≤ Tr (Σt (x)) − log det Σt (x) − log
ρX t (x) − f (Rt (x)) + f (x) − d . ν X (x)
Note the identity ⟨Id , Σt (x)⟩F = Tr (Σt (x)). Next, we consider −⟨Gt (x), Σt (x)⟩F . By (43), we have ρX X have positive t (x) det Gt (x) = ν X (R . Since Gt ⪰ 0 by the convexity of φt , and since ρX t and ν t (x)) density, we have det Gt (x) > 0, and hence Gt (x) ≻ 0. Next, we note Σt (x) ≻ 0 by Lemma 11. Y |X=x ρXY (x,y) Indeed, for all y ∈ Rd , we have ρt (y) = tρX (x) > 0. Then for any v ∈ Rd \ {0}, we have: t
⊤
v Σt (x)v = E
⊤
⊤
v Yt − v ut (x)
2
Z Xt = x =
v ⊤ y − v ⊤ ut (x)
2
Rd
Y |X=x
ρt
(y) dy > 0 ,
where the last inequality holds since the integrand is strictly positive outside the hyperplane {y ∈ Rd : v ⊤ y = v ⊤ ut (x)}, which has Lebesgue measure zero. Therefore, (1)
⟨Σt (x), Gt (x)⟩F ≥ log det Σt (x) + log det Gt (x) + d (2)
= log det Σt (x) + log
ρX t (x) + f (Rt (x)) − f (x) + d , ν X (x)
where (1) follows from the inequality ⟨A, B⟩F ≥ log det A + log det B + d which holds for A, B ≻ 0 [Bha97, Proposition II.3.20]; and (2) follows from the Monge–Ampere identity for Brenier maps and the definition ν X (x) ∝ e−f (x) . Rearranging gives the desired inequality. Term three.
We will show the following: −⟨∇f (x), vt (x)⟩ ≤ f (Rt (x)) − f (x) +
M ∥vt (x)∥2 . 2
Indeed, by M -semi-convexity of f , the following holds for any y, z ∈ Rd : Df (y, z) = f (y) − f (z) − ⟨∇f (z), y − z⟩ ≥ −
M ∥y − z∥2 . 2
Taking y = Rt (x) and z = x gives: f (Rt (x)) − f (x) − ⟨∇f (x), Rt (x) − x⟩ ≥ −
M ∥Rt (x) − x∥2 . 2
Rearranging and noting that vt (x) = x − Rt (x) gives the desired inequality. Summing the three terms above gives the result. 59
F.4
Bound on average KL divergence under regularity assumption
Lemma 22. Assume the setting of Lemma 3, and assume further that ρX 0 satisfies Assumption 1. Then for all 0 ≤ T < ∞, the following holds: 1
lim sup 2 h→0
1 2 X X X 2 X X W22 (ρX t+h , ν ) + 2 W2 (ρt−h , ν ) − W2 (ρt , ν ) h2 M 2 X X X X X ≤ 2KL(ρX W (ρ , ν ) . 0 ∥ ν ) − 3KL(ρt ∥ ν ) + 2 2 t
Therefore, 1 2 X X W (ρ , ν ) + 3 2 2 T
Z T
M X (T − t) KL(ρX t ∥ ν ) dt − 2
0
Z T
X (T − t) W22 (ρX t , ν ) dt
0
1 X 2 X X ≤ W22 (ρX 0 , ν ) + T KL(ρ0 ∥ ν ) . 2 Proof. Fix 0 < T < ∞. For t ∈ (0, T ), combining the results of Lemma 16 and Lemma 18 gives the claimed differential inequality: 1
lim sup 2 h→0
1 2 X X 2 X X X W22 (ρX M 2 X X t+h , ν ) − W2 (ρt , ν ) + 2 W2 (ρt−h , ν ) X X X ≤ 2KL(ρX W (ρ , ν ) . 0 ∥ν )−3KL(ρt ∥ν )+ 2 h 2 2 t
We wish to integrate this differential inequality twice in time. Since the left-hand side above is not a true second derivative, we proceed via concavity. We recall a standard characterization of concavity from [NP18, Theorem 1.4.7]: If a real-valued function g on an open interval (a, b) is continuous and satisfies lim inf h→0
g(t + h) − 2g(t) + g(t − h) ≤0 h2
for all t ∈ (a, b) ,
then g is concave on (a, b). Concretely, we define F : [0, T ] → R by: 1 X F (t) := W22 (ρX t , ν ) − G(t) 2 where G : [0, T ] → R is defined by: Z t M 2 X X X X X X G(t) := (t − s) 2KL(ρ0 ∥ ν ) − 3KL(ρs ∥ ν ) + W (ρ , ν ) ds . 2 2 s 0 X By Lemma 16, the map t 7→ 12 W22 (ρX t , ν ) is continuous on [0, T ] and is differentiable with one-sided X derivatives at the endpoints. By Lemma 11, W22 (ρX t , ν ) is uniformly bounded on the compact X X interval t ∈ [0, T ]. By Lemma 9, t 7→ KL(ρt ∥ ν ) is finite and continuous on [0, T ], and thus integrable. Hence, F is continuous on [0, T ]. Note G(t) is twice-continuously differentiable on t ∈ (0, T ), with
G(t + h) − 2G(t) + G(t − h) M 2 X X X X X W (ρ , ν ) . = G′′ (t) = 2KL(ρX 0 ∥ ν ) − 3KL(ρt ∥ ν ) + 2 h→0 h 2 2 t lim
X Combining the above, we have for F (t) = 21 W22 (ρX t , ν ) − G(t):
lim sup h→0
F (t + h) − 2F (t) + F (t − h) ≤ 0. h2 60
By the characterization from [NP18, Theorem 1.4.7], this shows F is concave on (0, T ). Since F is continuous on [0, T ], it extends as a concave function on the closed interval [0, T ]. Next, we will show the right derivative of F at 0 is 0, which will imply the desired inequality. By the one-sided version of the first-order derivative formula in Lemma 16, Z 1 2 X X d ⟨v0 (x), u0 (x)⟩ ρX W2 (ρt , ν ) = 0 (x) dx = 0 , dt t=0+ 2 d R where the last equality holds since ρXY = ρX 0 0 ⊗ γ, so u0 (x) = E[Y0 | X0 = x] = 0. On the other X X X hand, by the continuity of t 7→ KL(ρt ∥ ν ) and t 7→ W22 (ρX t , ν ), Z G(h) − G(0) M 2 X X 1 h X X X X lim (h − s) 2KL(ρ0 ∥ ν ) − 3KL(ρs ∥ ν ) + = lim W (ρ , ν ) ds = 0 . h→0 h→0 h 0 h 2 2 s Thus, F+′ (0) = 0. Since F is concave on [0, T ], its secant slopes are nonincreasing, and therefore, F (T ) − F (0) ≤ F+′ (0) = 0 , T which shows that F (T ) ≤ F (0). Expanding the definition of F gives the desired inequality.
G
Approximation argument for regularity of initial distribution
In this section, we provide an approximation argument to remove the regularity Assumption 1 on the initial distribution. Our treatment is inspired by the approach in [Lu26, Section 7].
G.1
Approximation of the initial distribution
In Lemma 23, we show how to approximate the initial distribution by a sequence of distributions which satisfy Assumption 1, with convergence in Wasserstein distance and KL divergence. d X X Lemma 23. Assume ν X ∈ P2,ac,fs (Rd ) is log-smooth, and ρX 0 ∈ P2,ac,fs (R ) satisfies KL(ρ0 ∥ ν ) < X X ∞. Then there exists a sequence {ρ0,n }n∈N such that each ρ0,n satisfies Assumption 1 for some X 0 < ζn < ξn < ∞, and the sequence ρX 0,n converges to ρ0 in the following sense: X lim W22 (ρX 0,n , ρ0 ) = 0 ,
n→∞
X X X lim KL(ρX 0,n ∥ ν ) = KL(ρ0 ∥ ν ) .
n→∞
X X Proof. For brevity, in this proof we omit the time index 0; that is, write ρX := ρX 0 and ρn := ρ0,n . We follow a standard regularization procedure of truncating the relative density, mollifying, and adding a positive floor.
Step 1: Truncation. For n ≥ 1, define
X
Define q := νρX . We first truncate q(x) in the argument and in the value. q̄n (x) := min {q(x), n} 1{∥x∥≤n} , Z mn := q̄n (x) ν X (x) dx , Rd
q̂n (x) :=
q̄n (x) , mn 61
X ρ̂X n (x) := q̂n (x) ν (x) .
R R X X Note that ρ̂X n is a probability density function: Rd ρ̂n (x) dx = Rd q̂n (x) ν (x) dx = 1. By cond struction, 0 ≤ q̄n (x) ≤ q(x), and limn→∞ q̄n (x) = q(x) for all x ∈ R . Since n 7→ qn (x) is non-decreasing for each x ∈ Rd , by the monotone convergence theorem we have limn→∞ mn = 1. Furthermore, by the construction of q̄n , we can bound: 1 + ∥x∥2 |q̄n (x) − q(x)| ≤ 1 + ∥x∥2 (q̄n (x) + q(x)) ρX (x) . ≤ 2 1 + ∥x∥2 q(x) = 2 1 + ∥x∥2 X ν (x) Since ρX ∈ P2,ac,fs (Rd ), the right-hand side above is integrable with respect to ν X (x). Thus, by the dominated convergence theorem, Z lim 1 + ∥x∥2 |q̄n (x) − q(x)| ν X (x) dx = 0 . (45) n→∞ Rd
X We now show the convergence of ρ̂X n to ρ . We can bound: Z X 1 + ∥x∥2 ρ̂X n (x) − ρ (x) dx RdZ (1) = 1 + ∥x∥2 |q̂n (x) − q(x)| ν X (x) dx d ZR q̄n (x) (2) = 1 + ∥x∥2 − q(x) ν X (x) dx m d n R Z Z (3) q̄ (x) n − q̄n (x) ν X (x) dx + 1 + ∥x∥2 |q̄n (x) − q(x)| ν X (x) dx ≤ 1 + ∥x∥2 mn d Rd Z R Z 1 = −1 1 + ∥x∥2 |q̄n (x) − q(x)| ν X (x) dx 1 + ∥x∥2 q̄n (x) ν X (x) dx + mn Rd Rd Z Z (4) 1 ≤ 1 + ∥x∥2 q(x) ν X (x) dx + 1 + ∥x∥2 |q̄n (x) − q(x)| ν X (x) dx , −1 mn Rd Rd
where (1) and (2) follow from definitions, (3) follows from triangle inequality, and (4) follows from the bound q̄n (x) ≤ q(x). R R Since mn → 1 and Rd 1 + ∥x∥2 q(x) ν X (x) dx = Rd 1 + ∥x∥2 ρX (x) dx < ∞, the first term above converges to 0 as n → ∞. By (45), the second term above also converges to 0 as n → ∞. R X Therefore, limn→∞ Rd 1 + ∥x∥2 ρ̂X n (x) − ρ (x) dx = 0. This implies: Z 1 X lim TV(ρ̂n , ρ) = lim ρ̂X n (x) − ρ (x) dx = 0 , n→∞ n→∞ 2 Rd Z Z lim ∥x∥2 ρ̂X (x) dx = ∥x∥2 ρX (x) dx . n n→∞ Rd
Rd
Since convergence in W2 distance is equivalent to weak convergence (which is implied by convergence in total variation) and convergence of second moment (see Section C.3 for a review), we conclude that X lim W2 ρ̂X = 0. n ,ρ n→∞
62
X → KL ρX ∥ ν X . Let ψ(z) := z log z, with ψ(0) := 0. Then We now show that KL ρ̂X n ∥ν Z X X X KL(ρ̂X ∥ ν ) = KL q̂ ν ∥ ν = ψ(q̂n (x)) ν X (x) dx n n d ZR q̄n (x) q̄n (x) log ν X (x) dx = mn Rd mn Z 1 = ψ(q̄n (x)) ν X (x) dx − log mn . mn Rd We bound the two terms separately. Define ψ(z)+ = max{ψ(z), 0} and ψ(z)− = min{ψ(z), 0}, so ψ(z) = ψ(z)+ + ψ(z)− . We apply the dominated convergence theorem to both parts separately. For ψ+ , note that z log z is increasing whenever z log z > 0. Moreover, 0 ≤ q̄n (x) ≤ q(x) and limn→∞ q̄n (x) = q(x). Therefore, limn→∞ ψ(q̄n (x))+ = ψ(q(x))+ . Since KL(ρX ∥ ν X ) < ∞, ψ(q)+ is integrable under ν X . Therefore, by the dominated convergence theorem: Z Z ψ(q(x))+ ν X (x) dx . lim ψ(q̄n (x))+ ν X (x) dx = n→∞ Rd
Rd
For ψ− , note that −e−1 ≤ ψ(q(x))− ≤ 0. By the dominated convergence theorem: Z Z lim ψ(q̄n (x))− ν X (x) dx = ψ(q(x))− ν X (x) dx . n→∞ Rd
Rd
R R Combining both parts gives limn→∞ Rd ψ(q̄n (x)) ν X (x) dx = Rd ψ(q(x)) ν X (x) dx. Finally, since limn→∞ mn = 1, we obtain Z 1 X X lim KL ρ̂X ∥ ν = lim ψ(q̄ (x)) ν (x) dx − log m n n n n→∞ n→∞ mn Rd Z = ψ(q(x)) ν X (x) dx − 0 Rd = KL ρX ∥ ν X . This shows we can truncate the density and maintain convergence in W2 and KL divergence. Step 2: Mollification. We now smooth the density to guarantee differentiability. By the construction in Step 1, we know that supp(ρ̂n ) = Bn := {x ∈ Rd : ∥x∥ ≤ n}. Let η : Rd → R denote the following standard mollifier: 1 C exp if ∥x∥ < 1 , ∥x∥2 −1 η(x) := 0 if ∥x∥ ≥ 1 , R where C ∈ (0, ∞) is a constant such that Rd η(x) dx = 1. By construction, η is smooth (infinitely differentiable) and compactly supported, with η(x) ≥ 0 for all x ∈ Rd , and supp(η) ⊆ B1 . For ε > 0, define ηε : Rd → R by x ηε (x) := ε−d η , ε X and define ρ̂X n,ε as the convolution of ηε and ρ̂n : X ρ̂X n,ε := ηε ∗ ρ̂n .
63
Then by construction, we have the following properties (see [Eva10, Theorem 7, Appendix C.5]): ∞ d 1 X the probability density function ρ̂X n,ε (x) ∈ Cc (R ) , supp(ρ̂n,ε ) ⊆ Bn+ε , and satisfies Z X (46) lim ρ̂X n,ε (x) − ρ̂n (x) dx = 0 . ε→0 Rd
We define the corresponding relative density qn,ε with respect to ν X by qn,ε (x) :=
ρ̂X n,ε (x) . ν X (x)
Then qn,ε is supported on Bn+ε , and since ν X ∈ P2,ac,fs (Rd ) and its density ν X (x) ∈ C 1 (Rd ), qn,ε ∈ Cc1 (Rd ). Since qn,ε is compactly supported, there exists ξn,ε ∈ (0, ∞) such that for all x ∈ Rd : 0 ≤ qn,ε (x) ≤ ξn,ε . We will show the convergence in W2 and KL divergence. For W2 , we use the property (46) X above. For 0 < ε < 1, since both ρ̂X n,ε and ρ̂n are supported in Bn+1 , we have: Z Z X 2 X lim 1 + ∥x∥2 ρ̂X (x) − ρ̂ (x) dx ≤ lim 1 + (n + 1) ρ̂X n,ε n n,ε (x) − ρ̂n (x) dx = 0 . ε→0 Rd
ε→0
Rd
As in Step 1, this implies convergence in total variation distance and in second moment: X lim TV ρ̂X n,ε , ρ̂n = 0 , ε→0 Z Z lim ∥x∥2 ρ̂X (x) dx = ∥x∥2 ρ̂X n,ε n (x) dx . ε→0 Rd
Rd
Therefore, this implies convergence in Wasserstein distance: X lim W2 (ρ̂X n,ε , ρ̂n ) = 0 .
ε→0
We now show convergence in the KL divergence. We can split: Z Z X X X X X KL ρ̂n,ε ∥ ν = ρ̂n,ε (x) log ρ̂n,ε (x) dx − ρ̂X n,ε (x) log ν (x) dx . Rd
(47)
Rd
We will show both terms converge. For the first term in (47), recall ψ(z) = z log z. We can bound: Z Z Z Z X X X ψ(ρ̂n,ε (x)) dx = ψ (ηε ∗ ρ̂n )(x) dx ≤ ηε ∗ ψ(ρ̂n ) (x) dx = ψ(ρ̂X n (x)) dx , Rd
Rd
Rd
Rd
where the inequality follows from applying Jensen’s inequality pointwise to the integrand since ψ is convex and ηε is a probability distribution, and the last equality follows from expanding the convolution and again using the fact that ηε is a probability distribution. Therefore, expanding the definition of ψ and sending ε → 0, we obtain: Z Z X X X lim sup ρ̂n,ε (x) log ρ̂n,ε (x) dx ≤ ρ̂X n (x) log ρ̂n (x) dx . ε→0
Rd
Rd
1
This notation means that ρ̂X n,ε (x) is smooth (infinitely differentiable and compactly supported), see Section C.1 for a review of the notations.
64
We now show the reverse inequality. Set K := Bn+1 , and let |K| = Vol(K). Then for 0 < ε < 1, X d both ρ̂X n,ε and ρ̂n are supported in K. Define a probability distribution µK on R with density 1 X X µK (x) := |K| 1K (x). Since ρ̂n,ε converges weakly to ρ̂n as ε → 0, the lower semicontinuity of KL divergence gives Z X X ρ̂X n (x) log ρ̂n (x) dx + log |K| = KL ρ̂n ∥ µK Rd ≤ lim inf KL ρ̂X n,ε ∥ µK ε→0 Z X X = lim inf ρ̂n,ε (x) log ρ̂n,ε (x) dx + log |K| . ε→0
Rd
Canceling log |K| on both sides yields: Z Z X X ρ̂n (x) log ρ̂n (x) dx ≤ lim inf ε→0
Rd
Rd
Therefore, we conclude that Z Z X X lim ρ̂n,ε (x) log ρ̂n,ε (x) dx = ε→0 Rd
Rd
X ρ̂X n,ε (x) log ρ̂n,ε (x) dx .
X ρ̂X n (x) log ρ̂n (x) dx .
X For the second term in (47), since ρ̂X n,ε is supported in Bn+1 and log ν is bounded on this compact X X set, the weak convergence of ρ̂n,ε to ρ̂n implies Z Z X X X lim ρ̂n,ε (x) log ν (x) dx = ρ̂X n (x) log ν (x) dx . ε→0 Rd
Rd
Combining both terms, we conclude that X X lim KL ρ̂X = KL ρ̂X . n,ε ∥ ν n ∥ν
ε→0
The convergence above holds for each fixed n ≥ 1, as ε → 0. Now for each n ≥ 1, we may choose εn ∈ (0, 1) sufficiently small such that 1 X W2 ρ̂X , n,εn , ρ̂n ≤ n
and
1 X X KL ρ̂X − KL ρ̂X ≤ . n,εn ∥ ν n ∥ν n
Then define X ρ̃X n (x) := ρ̂n,εn (x) ,
and
q̃n (x) :=
(48)
ρ̂X n,εn (x) . ν X (x)
∞ d X 1 d By the preceding argument, we know that ρ̂X n,ε (x) ∈ Cc (R ) and ν (x) ∈ C (R ). Therefore, 1 d q̃n ∈ Cc (R ). Since q̃n is continuous and has compact support, it is uniformly bounded, and 0 ≤ q̃n (x) ≤ ξn < ∞ for all x ∈ Rd , for some ξn ∈ (0, ∞). Furthermore, by (48), ρ̃X n satisfies:
X lim W2 ρ̃X = 0, n ,ρ
n→∞
and
X lim KL ρ̃X = KL ρX ∥ ν X . n ∥ν
n→∞
65
Step 3: Adding a positive floor. Finally, we add the lower bound to the relative density to make the distribution bounded below. For ζ ∈ (0, 1), define qn,ζ (x) := (1 − ζ)q̃n (x) + ζ , X ρX n,ζ (x) := qn,ζ (x) ν (x) .
R X Note that ρX n,ζ is still a probability density function: Rd ρn,ζ (x) dx = 1, and qn,ζ satisfies 0 < ζ ≤ qn,ζ (x) ≤ (1 − ζ)ξn + ζ < ∞, so qn,ζ ∈ Cb1 (Rd ). We next show the convergence in W2 and KL divergence as ζ → 0, for each fixed n. Note ρX n,ζ X + ζν X . Then the joint convexity of W2 gives: is a mixture distribution: ρX = (1 − ζ)ρ̃ n 2 n,ζ X 2 X X X W22 ρX n,ζ , ρ̃n = W2 (1 − ζ)ρ̃n + ζν , ρ̃n X 2 X X ≤ (1 − ζ)W22 ρ̃X n , ρ̃n + ζW2 ρ̃n , ν X = ζW22 ρ̃X . n ,ν X = 0. 2 X X X Since ρ̃X n has finite second moment, W2 (ρ̃n , ν ) < ∞. Therefore, limζ→0 W2 ρn,ζ , ρ̃n For the convergence in KL divergence, recall ψ(z) = z log z. By the convexity of ψ and since ψ(1) = 0, we have ψ (qn,ζ (x)) = ψ ((1 − ζ)q̃n (x) + ζ) ≤ (1 − ζ)ψ(q̃n (x)) + ζψ(1) = (1 − ζ)ψ(q̃n (x)). Integrating against ν X gives Z Z X X X X = ψ (qn,ζ (x)) ν (x) dx ≤ (1 − ζ) ψ (q̃n (x)) ν X (x) dx = (1 − ζ) KL ρ̃X . KL ρn,ζ ∥ ν n ∥ν Rd
Rd
Thus, X X lim sup KL ρX ≤ KL ρ̃X . n ∥ν n,ζ ∥ ν ζ→0
X On the other hand, since we have already shown that ρX n,ζ converges to ρ̃n in W2 , and hence weakly, by the lower semicontinuity of KL divergence, we have
X X . KL ρ̃X ≤ lim inf KL ρX n ∥ν n,ζ ∥ ν ζ→0
Combining the previous two inequalities shows that for each n ≥ 1, X X lim KL ρX = KL ρ̃X . n ∥ν n,ζ ∥ ν
ζ→0
Now for each fixed n ≥ 1, we may choose ζn ∈ (0, 1/n) sufficiently small such that 1 X W2 ρX , n,ζ , ρ̃n ≤ n
and
1 X X KL ρX − KL ρ̃X ≤ . n ∥ν n,ζ ∥ ν n
Finally, we define qn := qn,ζn = (1 − ζn )q̃n + ζn ,
and
X ρX n := qn ν .
Since qn,ζn ∈ Cb1 (Rd ), qn ∈ Cb1 (Rd ) as well, and 0 < ζn ≤ qn (x) ≤ (1 − ζn )ξn + ζn < ∞ for all x ∈ Rd . Furthermore, by the triangle inequality, X X X X X lim W2 ρX = lim W2 ρX ≤ lim W2 ρX n ,ρ n,ζn , ρ n,ζn , ρ̃n + W2 ρ̃n , ρ
n→∞
n→∞
n→∞
66
≤ lim
n→∞
1 X + W2 ρ̃X n ,ρ n
= 0.
Similarly, lim KL ρX ∥ ν X − KL ρX ∥ ν X n n→∞ X − KL ρX ∥ ν X = lim KL ρX n,ζn ∥ ν n→∞ ≤ lim KL ρX ∥ ν X − KL ρ̃X ∥ ν X + KL ρ̃X ∥ ν X − KL ρX ∥ ν X n n n,ζ n n→∞ 1 X X X X ≤ lim + KL ρ̃n ∥ ν − KL ρ ∥ ν n→∞ n = 0. Thus, we have constructed an approximating sequence of continuously differentiable distributions X X ρX n with bounded relative density qn (x) = ρn (x)/ν (x) satisfying the desired convergence in W2 distance and KL divergence.
G.2
Finite-time Wasserstein stability of the Hamiltonian flow
We recall that Ψt is the solution of the Hamiltonian flow at time t and ΨX t is its X-marginal, i.e., if we flow via the Hamiltonian flow from (X0 , Y0 ) to reach (Xt , Yt ), then Ψt (X0 , Y0 ) = (Xt , Yt ) and ΨX t (X0 , Y0 ) = Xt . Lemma 24. Assume ν X is log-smooth and M -semi-log-concave for some 0 ≤ M < ∞. Given d X X X ρX 0 ∈ P2,ac,fs (R ) with KL(ρ0 ∥ ν ) < ∞, let {ρ0,n }n∈N be the approximating sequence of regular distributions constructed in Lemma 23. For T ∈ (0, ∞) and t ∈ [0, T ], define X X ρX t,n := (Ψt )# (ρ0,n ⊗ γ) ,
and
X X ρX t := (Ψt )# (ρ0 ⊗ γ) .
Then X lim sup W22 (ρX t,n , ρt ) = 0 .
n→∞ t∈[0,T ]
X 2 = W2 (ρX , ρX ). Proof. Let (X0,n , X0 ) be an optimal coupling of ρX 0,n and ρ0 , so that E ∥X0,n − X0 ∥ 2 0,n 0 Let Y0 ∼ γ be independent of (X0,n , X0 ). We use the same velocity Y0 for both initial positions. Define (Xt,n , Yt,n ) := Ψt (X0,n , Y0 ) , and (Xt , Yt ) := Ψt (X0 , Y0 ) . X Then by definition, Xt,n ∼ ρX t,n and Xt ∼ ρt . By assumption, ν X is L-log-smooth for some L ∈ (0, ∞). Then by Lemma 7, for every t ∈ [0, T ],
∥Xt,n − Xt ∥ ≤ ∥Ψt (X0,n , Y0 ) − Ψt (X0 , Y0 )∥ ≤ e(1+L)t ∥(X0,n , Y0 ) − (X0 , Y0 )∥ = e(1+L)t ∥X0,n − X0 ∥ . Therefore, for every t ∈ [0, T ], X X W22 ρX ≤ E ∥Xt,n − Xt ∥2 ≤ e2(1+L)T E ∥X0,n − X0 ∥2 = e2(1+L)T W22 ρX . t,n , ρt 0,n , ρ0 X ≤ e2(1+L)T W2 ρX , ρX . Taking supremum over t ∈ [0, T ], we get supt∈[0,T ] W22 ρX t,n , ρt 2 0,n 0 Therefore, by Lemma 23, X X lim sup W22 ρX ≤ e2(1+L)T · lim W22 ρX t,n , ρt 0,n , ρ0 = 0 . n→∞ t∈[0,T ]
n→∞
67
G.3
Lower semicontinuity of the integral functional
Lemma 25. Assume the same set-up and definitions as in Lemma 24. For T ∈ (0, ∞) and n ≥ 1, define: Z Z T 1 2 X X M T X X X (T − t) W22 ρX dt , FT := W2 ρT , ν + 3 (T − t) KL ρt ∥ ν dt − t ,ν 2 2 0 0 Z Z T 1 M T X X X X (T − t) W22 ρX dt . FT,n := W22 ρX , ν + 3 (T − t) KL ρ ∥ ν dt − t,n , ν T,n t,n 2 2 0 0 Then FT < ∞, FT,n < ∞, and FT ≤ lim inf n→∞ FT,n . X Proof. We first check the functionals are well-defined. By Lemma 9, the map t 7→ KL(ρX t ∥ ν ) is X X X continuous on [0, T ], and KL(ρX t ∥ ν ) ≤ KL(ρ0 ∥ ν ) < ∞ for all t ∈ [0, T ]. Therefore, the integral RT X X 0 (T − t) KL(ρt ∥ ν ) dt is finite. On the other hand, the sequence of approximating distribution X X ≤ ξ , therefore KL(ρX ∥ ν X ) ≤ KL(ρX ∥ ν X ) < ∞ ρX n t,n 0,n we constructed satisfies 0 < ζn ≤ ρ0,n /ν 0,n for every n, and FT and FT,n are finite. X X 2 X Next, we prove the lower semicontinuity. By Lemma 24, limn→∞ W22 (ρX T,n , ν ) = W2 (ρT , ν ). X For each t ∈ [0, T ], by Lemma 24 we know ρX t,n converges to ρt in W2 distance, and hence also X X . weakly. By the lower semicontinuity of KL divergence, KL(ρt ∥ ν X ) ≤ lim inf n→∞ KL ρX t,n ∥ ν Therefore, Z T Z T X (T − t)KL(ρX ∥ ν ) dt ≤ lim inf (T − t) KL ρX ∥ ν X dt t t,n n→∞ 0 0 Z T X dt , ≤ lim inf (T − t) KL ρX t,n ∥ ν n→∞
0
where the second inequality follows from Fatou’s Lemma. Using the bound |a2 − b2 | = |a − b| · |a + b| ≤ 2|a| · |a − b| + (a − b)2 , and using triangle inequality X X X X X to bound W2 (ρX t,n , ν ) − W2 (ρt , ν ) ≤ W2 (ρt,n , ρt ), we have: X 2 X X lim sup W22 (ρX t,n , ν ) − W2 (ρt , ν )
n→∞ t∈[0,T ]
X X X X X ≤ 2 lim sup W2 (ρX t , ν ) · W2 (ρt,n , ν ) − W2 (ρt , ν ) + lim sup
n→∞ t∈[0,T ]
n→∞ t∈[0,T ]
X X X W2 (ρX t,n , ν ) − W2 (ρt , ν )
2
X X X 2 X X ≤ 2 lim sup W2 (ρX t , ν ) · W2 (ρt,n , ρt ) + lim sup W2 (ρt,n , ρt ) n→∞ t∈[0,T ]
n→∞ t∈[0,T ]
= 0, where the last step again follows from Lemma 24. Therefore, Z Z M T M T 2 X X X lim (T − t) W2 (ρt,n , ν ) dt = (T − t) W22 (ρX t , ν ) dt . n→∞ 2 2 0 0 Combining the arguments above, we obtain the desired claim: Z T Z M T 1 2 X X X X X FT = W2 (ρT , ν ) + 3 (T − t) KL(ρt ∥ ν ) dt − (T − t) W22 (ρX t , ν ) dt 2 2 0 0 Z T Z 1 M T X X X X ≤ lim W22 (ρX , ν ) + 3 lim inf (T − t) KL ρ ∥ ν dt − lim (T − t) W22 (ρX T,n t,n t,n , ν ) dt n→∞ 2 n→∞ n→∞ 2 0 0 68
≤ lim inf n→∞
1 2 X W (ρ , ν X ) + 3 2 2 T,n
Z T (T − t) KL
X ρX t,n ∥ ν
0
M dt − 2
Z T
X (T − t) W22 (ρX t,n , ν ) dt
0
= lim inf FT,n . n→∞
G.4
Convergence of the initial functional
Lemma 26. Assume the same set-up and definitions as in Lemma 24. For T ∈ (0, ∞), define 1 X 2 X X GT := W22 (ρX 0 , ν ) + T KL(ρ0 ∥ ν ) , 2 1 X 2 X X GT,n := W22 (ρX . 0,n , ν ) + T KL ρ0,n ∥ ν 2 Then GT < ∞, GT,n < ∞, and limn→∞ GT,n = GT . X X X have finite second Proof. Clearly GT < ∞ since we assume KL(ρX 0 ∥ ν ) < ∞ and both ρ0 and ν X X moments. Similarly, GT,n < ∞ since KL(ρ0,n ∥ ν ) < ∞ by the construction in Lemma 23, and ρX 0,n has a finite second moment. Next, by Lemma 23 and the triangle inequality, X X X X X lim W2 (ρX 0,n , ν ) − W2 (ρ0 , ν ) ≤ lim W2 (ρ0,n , ρ0 ) = 0 .
n→∞
n→∞
X 2 X X X X = Therefore, limn→∞ W22 (ρX 0,n , ν ) = W2 (ρ0 , ν ). Again by Lemma 23, limn→∞ KL ρ0,n ∥ ν X KL(ρX 0 ∥ ν ). Therefore, 1 1 2 X X 2 X X X 2 X X W2 (ρ0,n , ν ) + T KL ρ0,n ∥ ν = W22 (ρX lim GT,n = lim 0 , ν ) + T KL(ρ0 ∥ ν ) = GT . n→∞ n→∞ 2 2
G.5
Proof of the bound on average KL divergence along Hamiltonian flow
With the preparation above, we finally prove Lemma 3. Proof of Lemma 3. We recall the definitions of FT and FT,n from Lemma 25, and the definitions of GT and GT,n from Lemma 26. By Lemma 25, FT ≤ lim inf n→∞ FT,n . By applying Lemma 22 to the approximating sequence of distributions (which satisfy the regularity Assumption 1), we have FT,n ≤ GT,n , so lim inf n→∞ FT,n ≤ limn→∞ GT,n . Furthermore, by Lemma 26, limn→∞ GT,n = GT . Combining the three inequalities above, we obtain FT ≤ lim inf FT,n ≤ lim GT,n = GT . n→∞
n→∞
Expanding the definitions of FT and GT completes the proof.
69
References [AAMN24] Dallas Albritton, Scott Armstrong, Jean-Christophe Mourrat, and Matthew Novack. Variational methods for the kinetic Fokker–Planck equation. Analysis & PDE, 17(6):1953–2010, 2024. [ACZ26] Jason M Altschuler, Sinho Chewi, and Matthew S Zhang. Shifted composition IV: Toward ballistic acceleration for log-concave sampling. In Proceedings of the 58th Annual ACM Symposium on Theory of Computing, pages 1739–1750, 2026. [AGS05] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows: in metric spaces and in the space of probability measures. Springer, 2005. [AGS24] Simon Apers, Sander Gribling, and Dániel Szilágyi. Hamiltonian monte carlo for efficient gaussian sampling: long and random steps. Journal of Machine Learning Research, 25(348):1–30, 2024. [BGJM11] Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Handbook of Markov chain Monte Carlo. CRC press, 2011. [Bha97] Rajendra Bhatia. Matrix analysis, volume 169 of Graduate Texts in Mathematics. Springer-Verlag, New York, 1997. [BRCKL26] Nawaf Bou-Rabee, Bob Carpenter, Tore Selland Kleppe, and Sifan Liu. The withinorbit adaptive leapfrog No-U-turn sampler. Journal of Machine Learning Research, 27(113):1–64, 2026. [BRE23] Nawaf Bou-Rabee and Andreas Eberle. Mixing time guarantees for unadjusted Hamiltonian Monte Carlo. Bernoulli, 29(1):75–104, 2023. [BREZ20] Nawaf Bou-Rabee, Andreas Eberle, and Raphael Zimmer. Coupling and convergence for Hamiltonian Monte Carlo. The Annals of applied probability, 30(3):1209–1250, 2020. [BRMW26] Nawaf Bou-Rabee, Siddharth Mitra, and Andre Wibisono. Tail-Sensitive KL and Rényi Convergence of Unadjusted Hamiltonian Monte Carlo via One-Shot Couplings. arXiv preprint arXiv:2601.09019, 2026. [BRO24] Nawaf Bou-Rabee and Stefan Oberdörster. Mixing of the No-U-Turn Sampler and the Geometry of Gaussian Concentration. arXiv preprint arXiv:2410.06978, 2024. [BRS23] Nawaf Bou-Rabee and Katharina Schuh. Convergence of unadjusted Hamiltonian Monte Carlo for mean-field models. Electronic Journal of Probability, 28:1–40, 2023. [BRSS17] Nawaf Bou-Rabee and Jesús Marı́a Sanz-Serna. Randomized Hamiltonian Monte Carlo. The Annals of Applied Probability, pages 2159–2194, 2017. [BRSS18] Nawaf Bou-Rabee and Jesús Maria Sanz-Serna. Geometric integrators and the Hamiltonian Monte Carlo method. Acta Numerica, 27:113–206, 2018.
70
[CCBJ18] Xiang Cheng, Niladri S Chatterji, Peter L Bartlett, and Michael I Jordan. Underdamped Langevin MCMC: A non-asymptotic analysis. In Conference on Learning Theory, pages 300–323. PMLR, 2018. [CCT18] José A. Carrillo, Young-Pil Choi, and Oliver Tse. Convergence to equilibrium in Wasserstein distance for damped Euler equations with interaction forces. Communications in Mathematical Physics, 365(1):329–361, October 2018. [CDMS23] Evan Camrud, Alain Durmus, Pierre Monmarché, and Gabriel Stoltz. Second order quantitative bounds for unadjusted generalized Hamiltonian Monte Carlo. arXiv preprint arXiv:2306.09513, 2023. [CDWY20] Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu. Fast mixing of Metropolized Hamiltonian Monte Carlo: Benefits of multi-step gradients. Journal of Machine Learning Research, 21(92):1–72, 2020. [CGH+ 17] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of Statistical Software, 76:1–32, 2017. [CGJ26] Yuansi Chen, Khashayar Gatmiry, and Minhui Jiang. When does metropolized Hamiltonian Monte Carlo provably outperform Metropolis-adjusted Langevin algorithm? arXiv preprint arXiv:2304.04724 v3, 2026. [CLS07] Eric Cancès, Frédéric Legoll, and Gabriel Stoltz. Theoretical and numerical comparison of some sampling methods for molecular dynamics. ESAIM: Mathematical Modelling and Numerical Analysis, 41(2):351–389, 2007. [CLTW25] Shi Chen, Qin Li, Oliver Tse, and Stephen J. Wright. Accelerating optimization over the space of probability measures. Journal of Machine Learning Research, 26(31):1–40, 2025. [CLW23] Yu Cao, Jianfeng Lu, and Lihan Wang. On explicit L2 -convergence rate estimate for underdamped langevin dynamics. Archive for Rational Mechanics and Analysis, 247(5), August 2023. [CT06] Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. WileyInterscience, second edition, 2006. [CV22] Zongchen Chen and Santosh S. Vempala. Optimal convergence rate of Hamiltonian Monte Carlo for strongly logconcave distributions. Theory of Computing, 18(9):1–18, 2022. [DKPR87] Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics Letters B, 195(2):216–222, 1987. [DP19] Dmitriy Drusvyatskiy and Courtney Paquette. Efficiency of minimizing compositions of convex functions and smooth maps. Mathematical Programming, 178(1–2):503–558, 2019. 71
[EG15] Lawrence C. Evans and Ronald F. Gariepy. Measure Theory and Fine Properties of Functions. Chapman and Hall/CRC, revised edition, 2015. [EL26] Andreas Eberle and Francis Lörler. Non-reversible lifts of reversible diffusion processes and relaxation times. Probability Theory and Related Fields, 194(1):173–203, 2026. [Eva10] Lawrence C. Evans. Partial Differential Equations. American Mathematical Society, 2 edition, 2010. [FLL26] Zexi Fan, Bowen Li, and Jianfeng Lu. Sharp hypocoercive convergence estimates for underdamped Langevin dynamics via the modified L2 method. arXiv preprint arXiv:2604.10068, 2026. [FW25] Qiang Fu and Andre Wibisono. Hamiltonian descent algorithms for optimization: Accelerated rates via randomized integration time. In Advances in Neural Information Processing Systems, volume 38, pages 93628–93685. Curran Associates, Inc., 2025. [GBMM25] Nicolaı̈ Gouraud, Pierre Le Bris, Adrien Majka, and Pierre Monmarché. HMC and underdamped Langevin united in the unadjusted convex smooth case. SIAM/ASA Journal on Uncertainty Quantification, 13(1):278–303, 2025. [GKG+ 26] Samuel Gruffaz, Kyurae Kim, Fares Guehtar, Hadrien Duval-Decaix, and Pacôme Trautmann. A Theoretical Comparison of No-U-Turn Sampler Variants: Necessary and Sufficient Convergence Conditions and Mixing Time Analysis under Gaussian Targets. arXiv preprint arXiv:2603.18640, 2026. [HG14] Matthew D Hoffman and Andrew Gelman. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res., 15(1):1593–1623, 2014. [Jia23] Qijia Jiang. On the dissipation of ideal hamiltonian monte carlo sampler. Stat, 12(1):e629, 2023. [JKO98] Richard Jordan, David Kinderlehrer, and Felix Otto. The variational formulation of the Fokker–Planck equation. SIAM J. Math. Anal., 29(1):1–17, 1998. [LL26] Bowen Li and Jianfeng Lu. Space-time log-Sobolev inequality and hypocoercive hypercontractivity for underdamped Langevin dynamics. arXiv preprint arXiv:2605.25083, 2026. [Lu26] Jianfeng Lu. A sharp hypocoercive entropy decay estimate for underdamped Langevin dynamics. arXiv preprint arXiv:2605.01933, 2026. [LW22] Jianfeng Lu and Lihan Wang. On explicit L2 -convergence rate estimate for piecewise deterministic Markov processes in MCMC algorithms. The Annals of Applied Probability, 32(2), April 2022. [Mac89] Paul B. Mackenzie. An Improved Hybrid Monte Carlo Method. Physics Letters B, 226:369–371, 1989.
72
[MCC+ 21] Yi-An Ma, Niladri S. Chatterji, Xiang Cheng, Nicolas Flammarion, Peter L. Bartlett, and Michael I. Jordan. Is there an analog of Nesterov acceleration for gradient-based MCMC? Bernoulli, 27(3):1942 – 1992, 2021. [Mon24] Pierre Monmarché. An entropic approach for Hamiltonian Monte Carlo: the idealized case. The Annals of Applied Probability, 34(2):2243–2293, 2024. [MS21] Oren Mangoubi and Aaron Smith. Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions: Continuous dynamics. The Annals of Applied Probability, 31(5):2019–2045, 2021. [MW26] Pierre Monmarché and Lihan Wang. On the entropic convergence for piecewise deterministic samplers: speedup and obstruction. arXiv preprint arXiv:2606.26086, 2026. [Nea11] Radford M. Neal. MCMC using hamiltonian dynamics. In Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng, editors, Handbook of Markov Chain Monte Carlo, chapter 5. Chapman and Hall/CRC, 2011. [NP18] Constantin P. Niculescu and Lars-Erik Persson. Convex Functions and Their Applications: A Contemporary Approach. CMS Books in Mathematics. Springer, Cham, second edition, 2018. [Obe25] Stefan Oberdörster. On Accelerated Mixing of the No-U-turn Sampler. arXiv preprint arXiv:2507.13259, 2025. [OV00] Felix Otto and Cédric Villani. Generalization of an Inequality by Talagrand and Links with the Logarithmic Sobolev Inequality. Journal of Functional Analysis, 173(2):361– 400, 2000. [OV01] Felix Otto and Cédric Villani. Comment on: “Hypercontractivity of Hamilton–Jacobi equations”’, by S. Bobkov, I. Gentil and M. Ledoux. Journal de Mathématiques Pures et Appliquées, 80(7):697–700, 2001. [PHF10] Anand Patil, David Huard, and Christopher J Fonnesbeck. PyMC: Bayesian stochastic modelling in Python. Journal of Statistical Software, 35:1–81, 2010. [PW24] Yury Polyanskiy and Yihong Wu. Information Theory: From Coding to Learning. Cambridge University Press, 2024. [RT96] Gareth O. Roberts and Richard L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 1996. [SBC16] Weijie Su, Stephen Boyd, and Emmanuel J. Candès. A Differential Equation for Modeling Nesterov’s Accelerated Gradient Method: Theory and Insights. Journal of Machine Learning Research, 17(153):1–43, 2016. [TPL19] Andrew R Teel, Jorge I Poveda, and Justin Le. First-order optimization algorithms with resets and Hamiltonian flows. In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 5838–5843. IEEE, 2019.
73
[Vil03] Cédric Villani. Topics in Optimal Transportation, volume 58. American Mathematical Soc., 2003. [Vil09] Cédric Villani. Optimal Transport: Old and New. Springer, 2009. [VW23] Santosh S. Vempala and Andre Wibisono. Rapid convergence of the Unadjusted Langevin Algorithm: Isoperimetry Suffices. In Geometric Aspects of Functional Analysis, pages 381–438. Springer, 2023. [Wan24] Jun-Kun Wang. Frictionless Hamiltonian Descent and Coordinate Hamiltonian Descent for Strongly Convex Quadratic Problems and Beyond. arXiv preprint arXiv:2402.13988, 2024. [Wib18] Andre Wibisono. Sampling as optimization in the space of measures: The Langevin dynamics as a composite optimization problem. In Conference on Learning Theory, pages 2093–3027. PMLR, 2018. [WL22] Yifei Wang and Wuchen Li. Accelerated information gradient flow. Journal of Scientific Computing, 90:11, 2022. [WRJ21] Ashia C Wilson, Ben Recht, and Michael I Jordan. A Lyapunov analysis of accelerated methods in optimization. Journal of Machine Learning Research, 22(113):1–34, 2021. [WSF+ 26] Xiuyuan Wang, Vishwak Srinivasan, Qiang Fu, Siddharth Mitra, Ashia Wilson, and Andre Wibisono. Accelerated convex optimization via Hamiltonian dynamics with deterministic integration time. In Proceedings of Thirty Ninth Conference on Learning Theory, volume 336 of Proceedings of Machine Learning Research, pages 6692–6742. PMLR, 2026. [WW23] Jun-Kun Wang and Andre Wibisono. Accelerating Hamiltonian Monte Carlo via Chebyshev Integration Time. In Proceedings of the Eleventh International Conference on Learning Representations, 2023.
74