Delocalization of bias in unadjusted Hamiltonian Monte Carlo and underdamped Langevin

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Delocalization of bias in unadjusted Hamiltonian Monte Carlo and underdamped Langevin

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

DELOCALIZATION OF BIAS IN UNADJUSTED HAMILTONIAN MONTE CARLO AND UNDERDAMPED LANGEVIN

arXiv:2607.15208v1 [stat.CO] 16 Jul 2026

YIFAN CHEN† , XIAOOU CHENG‡ , JONATHAN NILES-WEED‡ , AND JONATHAN WEARE‡ Abstract. Unadjusted samplers such as unadjusted Hamiltonian Monte Carlo and underdamped Langevin are well-known to be biased. Metropolis–Hastings adjustment has been conventionally incorporated into Hamiltonian Monte Carlo to eliminate the bias. However, this adjustment can significantly increase the iteration complexity due to the small step size required for reasonable Metropolis acceptance rates. In this work, we extend the delocalization of bias phenomenon, previously established for the overdamped Langevin algorithm, to these two unadjusted algorithms. We show that to control √ the W2 bias of any K-dimensional marginal of a high-dimensional distribution, O( K) integration steps suffice up to log d terms, assuming either weak or sparse interactions among variables. The discrete-time integrators here introduce technical difficulties beyond those of the overdamped setting, which we address through a broadly applicable matrix-polynomial framework that characterizes their propagators. Our result for the underdamped Langevin algorithm is valid for all large friction parameters, implying that the Leimkuhler-Matthews integrator for the overdamped Langevin dynamics also exhibits delocalization of bias.

1. Introduction Hamiltonian (or Hybrid) Monte Carlo (HMC) and underdamped Langevin (UL)1 are widely used to sample from high dimensional probability distributions. For a target distribution π ∝ exp(−V ), where V : Rd → R denotes the potential energy, both sampling dynamics evolve in terms of a position variable q, on which the potential V is defined, and an auxiliary momentum variable p. HMC relies on the Hamiltonian ordinary differential equation dqt = pt , dt (1.1) dpt = −∇V (qt ) . dt corresponding to the Hamiltonian H(q, p) = V (q)+|p|2ℓ2 /2, where |·|ℓ2 is the ℓ2 norm of a vector. The closely related continuous-time UL is described by the stochastic differential equation dqt = pt dt , p (1.2) dpt = −∇V (qt )dt − γpt dt + 2γdBt . †

Department of Mathematics, University of California, Los Angeles, CA 90095, USA Courant Institute, New York University, NY 10012, USA E-mail addresses: [email protected], [email protected], [email protected], [email protected]. 1 Note that the unadjusted Langevin algorithm (ULA) is instead a discretization of the overdamped Langevin dynamics. ‡

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

where γ > 0 is the friction parameter and Bt is a standard Brownian motion. Both continuous-time dynamics admit π ⊗ N (0, I) as their invariant measure on phase space (q, p), allowing samplers built upon them to target π via the q-marginal. Applying these continuous-time dynamics numerically requires time discretization. For HMC, one must also refresh the momentum regularly to ensure sampling ergodicity. The most commonly used discretization of (1.1) is the leap-frog scheme (also called velocity Verlet) [31], under which the (k + 1)-th step of HMC is (1) Sample p0 ∼ N (0, I) and set q0 = Xk . (2) Inner loop (leap-frog integration): for j = 0 to m − 1, h2 ∇V (qj ) , 2 h h pj+1 = pj − ∇V (qj ) − ∇V (qj+1 ) . 2 2 qj+1 = qj + hpj −

(1.3)

Set Xk+1 = qm . Here, h is the step size, and we also refer to one HMC iteration as one outer loop. For (1.2), a commonly used discretization is the BAOAB scheme [31], which is written as h2 hp h 1 − η 2 ξk , qk+1 = qk + (1 + η)pk − (1 + η)∇V (qk ) + 2 4 2 (1.4) p h h pk+1 = η(pk − ∇V (qk )) + 1 − η 2 ξk − ∇V (qk+1 ) . 2 2 The scheme (1.4) is constructed through an operator splitting. We denote the solution maps for the three components of (1.2) over a step size h as follows: Bh (q, p) = (q, p − h∇V (q)) , (1.5)

Ah (q, p) = (q + hp, p) , p Oh (q, p; ξ) = (q, ηp + 1 − η 2 ξ) .

Here, ξ ∼ N (0, I) is an independent Gaussian vector, and η = exp(−γh). Splitting methods for underdamped Langevin dynamics are defined with different strings of A, B, and O. By standard convention, operators employed symmetrically in the string are applied with a half step h/2, while those that appear only once are applied with a full step h. Specifically, the BAOAB scheme (1.4) is derived from the composition Bh/2 ◦ Ah/2 ◦ Oh (·; ξ) ◦ Ah/2 ◦ Bh/2 , which we abbreviate as BAOAB. We use sans-serif A, B, O and their compositions for the solution maps, while the same letters in roman type identify the corresponding scheme. For the frictionless case γ = 0, we have η = 1, so the Oh step becomes the identity, and the BAOAB scheme (1.4) reduces to the leap-frog scheme (1.3). 1.1. Motivation. The unadjusted HMC and UL algorithms are biased. Because of the finite step size h, when the marginal ρk of qk converges, its limit πh differs from the target π. Both HMC and UL can be modified with a Metropolis–Hastings (MH) accept/reject step to eliminate bias. In fact, HMC is most frequently applied with MH adjustment [42, 13, 47]. However, adjusted HMC and UL are well known to degrade as d increases and the step size h must be decreased to maintain a non-zero acceptance probability.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

In [17], we have identified a “delocalization of bias” phenomenon for the unadjusted overdamped Langevin algorithm: To control the bias of low-dimensional marginals, a step size scaling with the low dimension suffices, with only a logarithmic dependence on the full dimension d. The analysis in [17] relies on a novel W2,ℓ∞ metric, which captures the accuracy of low-dimensional marginals. For two probability measures µ and ν on Rd , their W2,ℓ∞ distance is defined as Z Å ã1/2 2 (1.6) W2,ℓ∞ (µ, ν) = min |x − y|ℓ∞ γ(dx, dy) , γ∈Π(µ,ν)

where | · |ℓ∞ is the ℓ∞ norm of a vector and Π(µ, ν) represents the set of measures on the joint space Rd × Rd that have marginals µ and ν. We note that, for any choice of K

(K) and ν (K) of µ and ν, satisfy coordinates, the corresponding marginal √ distributions µ P (jt ) − (K) (K) the upper bound W2 (µ , ν ) ≤ KW2,ℓ∞ (µ, ν), since K|x − y|2ℓ∞ ≥ K t=1 |x y (jt ) |2 for any 1 ≤ jt ≤ d. In particular, W2,ℓ∞ (µ, ν) upper bounds the distance between any one-dimensional marginals. Our analysis for unadjusted overdamped Langevin in [17] proceeds via a coupling argument under the W2,ℓ∞ distance. In this paper, we extend that analysis to unadjusted HMC and UL, which are both much more commonly used for complicated, high-dimensional sampling applications than overdamped Langevin. This extension introduces significant technical challenges, as the dynamics now evolve on the full phase space (q, p) with an auxiliary momentum variable. The coupling framework in [17] requires analyzing the discrete dynamics over multiple steps to bound the accumulation of discretization error. However, the propagator matrices on the (q, p) space exhibit intricate position-momentum coupling, making their multi-step composition difficult to control. To this end, we introduce a novel matrix polynomial representation of the propagators. Such a representation allows for precise characterization of the multi-step dynamics and enables sharp error bounds under the ℓ∞ norm.

1.2. Literature review. Hamiltonian and underdamped Langevin dynamics have long been foundational to the simulation of physical systems. HMC was originally introduced in [23] for particle physics, and was later adapted and popularized for modern statistics and machine learning by [42]. Today, it serves as one of the standard samplers in probabilistic programming packages like Stan [13] and PyMC [47]. In the realm of molecular dynamics (MD), unadjusted HMC is closely related to the Andersen thermostat [4]. Similarly, UL is commonly employed as a robust thermostat [27, 3, 31]. The BAOAB splitting scheme was introduced in [32] and variants have since been implemented widely in MD software packages [25, 1, 29, 14, 44, 43, 50]. To maintain stability in integration of Hamilton’s ODE (1.1) or the underdamped Langevin SDE (1.2), the step size h must be chosen smaller than the shortest vibrational period among all processes in the system (e.g. about a femtosecond in MD). For adjusted HMC, by contrast, [6] establishes the scaling h ∼ d−1/4 as the requirement for a nonvanishing acceptance probability as d → ∞, under a product measure with bounded fourth-order derivatives of the potential. The same scaling was later shown for an adjusted UL scheme in [45], and also appears in mixing time guarantees for adjusted HMC under warm starts [18]. Recent work [51] shows that unadjusted HMC generates a warm start in O(d1/4 ) steps up to logarithmic factors.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Unadjusted algorithms avoid the rejection bottleneck but incur a discretization bias. A rich body of literature has emerged focusing on the non-asymptotic analysis of unadjusted HMC and UL algorithms, establishing bias bounds as the number of iterations goes to infinity. These analyses rely on various distance metrics and target measure assumptions. The most fundamental setting, detailed in Assumption 1.1, requires the potential to be a C 2 function that is strongly convex with Lipschitz gradient. However, because commonly employed integrators possess weak second-order accuracy [33, 8], fully demonstrating their accuracy often requires bounding the higher-order derivatives of the potential [39]. Under these frameworks, the theoretical guarantees show power laws h ∼ d−c for some order c > 0 depending on the specific assumptions. For unadjusted HMC, [36] requires h ∼ d−1/2 to bound the Wasserstein distance, while [38] achieves an h ∼ d−1/4 scaling for a high-probability bound under a third-order regularity condition. Furthermore, bias bounds arising from convergence guarantees have been derived in total variation distance (h ∼ d−3/4 ) [7], Kullback-Leibler (KL) divergence (h ∼ d−3/4 ) [9], and Rényi divergence (h ∼ d−3/2 ) [9]. For unadjusted UL, [34] establishes the Wasserstein bias for the BAOAB scheme, requiring h ∼ d−1/2 under Assumption 1.1, with improved dimensional dependence achieved under stricter smoothness assumptions. Wasserstein bias for different integrators of UL was studied in [20, 22, 39]. [39] also establishes a total variation bias bound through a Wassersteintotal variation regularization property. Under KL divergence, [35] similarly requires h ∼ d−1/2 . These analyses have also been generalized to non-convex settings using techniques such as reflection coupling [26, 48, 16] and entropic approaches [41, 12]. Recent literature has increasingly unified the theoretical analysis of HMC and UL algorithms. For instance, under the umbrella framework of “generalized HMC” (gHMC), [28] determines the optimal scaling of the friction parameter and the HMC integration time to attain accelerated convergence regarding the condition number for Gaussian targets, while [12] develops KL convergence guarantees for gHMC. Moving beyond traditional fixed-step integrators, [49] introduces the Randomized Midpoint Method (RMM) for UL dynamics. They establish an improved dimension dependence h ∼ d−1/3 in the Wasserstein distance. Inspired by this mechanism, [8] achieves the same d−1/3 scaling for unadjusted HMC by employing a stratified Monte Carlo integrator. Crucially, the aforementioned guarantees quantify the error across all variables of the joint system. In contrast, our work builds upon [17, 30] to focus exclusively on controlling the error of low-dimensional marginals, which is a more practical notion of error in high-dimensional applications where the majority of coordinates act as nuisance variables. We previously introduced the term “delocalization of bias” in [17] to describe the phenomenon where the W2 bias of low-dimensional marginals remains nearly independent of the full system dimension. Our analysis there focused on the unadjusted (overdamped) Langevin algorithm, employing a coupling framework to establish Wasserstein bias bounds for the marginals. [30] investigates this phenomenon under KL divergence using a hierarchical analysis of marginal relative entropies. In [21], a related question is studied for graphical models, where dimension-independent bounds are established, through Stein’s method, on the marginal error between a target distribution and its localized approximation. Prior to our work in [17], [11] and [24]

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

observed that large step sizes remain viable in high dimensions if the observable of interest has a Lipschitz constant that decays with dimension, thereby canceling out the dimensional scaling of the full system’s Wasserstein bias (e.g., an averaged observable P f (x) = d1 di=1 ϕ(x(i) ) with an ℓ2 -Lipschitz ϕ). In such regimes, a constant step size is sufficient to achieve a bounded error regardless of dimension. Besides these analysis efforts, there are ample numerical evidence demonstrating the viability of dimension-free step size choices [15, 46]. 1.3. Main results. In [17], we have shown for overdamped Langevin that the delocalization of bias holds for product measures, for Gaussian measures, and for distributions with sparse interactions between variables. We also provided a negative example showing that this effect does not hold universally. In this paper, we establish the delocalization effect for the unadjusted HMC and UL algorithms under analogous assumptions, namely weak or sparse interactions, or Gaussian targets. Our analysis is conducted under the following strongly log-concave and log-smooth assumption. Assumption 1.1. Let π ∝ exp(−V ) with V ∈ C 2 (Rd ). Assume V is α-strongly convex and β-smooth such that αI ⪯ ∇2 V (x) ⪯ βI for any x ∈ Rd , where 0 < α ≤ β < ∞. Throughout, unadjusted HMC uses the leap-frog integrator (1.3) and unadjusted UL uses the BAOAB scheme (1.4). As a warm-up, we present the delocalization of bias result for Gaussian distributions. Since the BAOAB scheme is unbiased for any Gaussian distribution [33], only unadjusted HMC is discussed in this example. We employ the explicit formula for the biased Gaussian distribution and obtain a W2,ℓ∞ bias bound directly. The proof appears in Appendix A. Example 1.2 (Gaussian distributions). Consider π = N (µ, Σ) where µ ∈ Rd and √ −1 αI ⪯ Σ ⪯ βI. Then, for h < 1/ β, unadjusted HMC, when it converges, has the following bias: » » (1.7) W2,ℓ∞ (πh , π) = O(h2 β log(2d)) = O(h log(2d)) . As in the case of overdamped Langevin dynamics [17], the proof of Example 1.2 offers no insight into which target distributions beyond the Gaussian family evince the delocalization of bias effect for unadjusted HMC or UL. Our main results describe two different sets of structural assumptions on V under which this effect appears. For general distributions satisfying Assumption 1.1, our first results address the “weak interactions” setting of [30], made precise below. Assumption 1.3 (weak interactions). The Hessian ∇2 V admits the decomposition ∇2 V = ∇2 V (D) + ∇2 V (O) where ∇2 V (D) = diag(∇2 V ) is the diagonal part and ∇2 V (O) is the off-diagonal part. We assume there exists a sufficiently small constant C (O) > 0, independent of α, β, h, and d, such that the the off-diagonal part satisfies |∇2 V (O) |ℓ∞ ≤ C (O) α uniformly over Rd . Assumption 1.1 ensures that every diagonal entry satisfies [∇2 V (D) ]ii ≥ α. The interactions defined in Assumption 1.3 are therefore weak, as the total intensity of

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

interactions between any given variable and all other variables, measured by off-diagonal row sums, is small compared to the independent components on the diagonal. We now state our main result under weak interactions. The proofs for the two algorithms are given in Appendix F. Theorem 1.4 (weak interactions). Let Assumptions 1.1 and 1.3 hold. For either un√ √ adjusted HMC with mh = 1/ 20β and C (O) = 1/50, or unadjusted UL with h ≤ 21−η β and C (O) = 1/20, the resulting invariant distribution πh satisfies Å » ã β (1.8) W2,ℓ∞ (πh , π) = O h log(2d) . α Essentially, the distribution with weak interactions constitutes a perturbation of the product measure. The delocalization effect holds for product measures, since the dynamics are separable for each variable, leading to a dimension-free bias in each coordinate. This effect persists as long as the interactions between variables are sufficiently weak. The second setting of interest is the case of “sparse interactions” studied in [17], where each variable interacts directly with only a few others. Because the interactions are local, the influence of one variable spreads only gradually through the interaction structure over multiple steps, keeping the marginal bias small. This interaction structure is described by a graphical model, with each variable a node, as illustrated in Figure 1. Let G be an undirected graph with d nodes, labeled by 1 ≤ i ≤ d, and write i ∼ j when nodes i and j are connected by an edge. The neighborhood N(i) is the set of nodes connected to i. Here the neighborhood relationship is symmetric, and without loss of generality, we take i ∼ i for all 1 ≤ i ≤ d. We quantify the locality of interactions through a growth condition on the neighborhood sizes. Define recursively Nk (i) = {1 ≤ j ≤ d : ∃ l ∈ Nk−1 (i), such that j ∼ l}, for k ≥ 2, with N1 (i) := N(i), and write sk := max1≤i≤d |N2k (i)| for the relevant neighborhood cardinality where sk ≥ 1 by definition. We make the following assumption on the potential. Assumption 1.5. The potential has the form (1.9)

V (x) =

d X

Vi (X i ) ,

i=1

where X i = {x(j) : j ∈ N(i)} and each Vi depends only on the variables in X i . Moreover, the sparsity parameter sk satisfies the polynomial growth bound sk ≤ C(k + 1)n for constants C, n > 0. The main result under sparse interactions is as follows. The proofs are given in Appendix F. Theorem 1.6 (sparse interactions). Let Assumptions 1.1 and 1.5 hold. For either √ √ , the resulting unadjusted HMC with mh = 1/ 20β, or unadjusted UL with h ≤ 21−η β invariant distribution πh satisfies Å Å ãã n +1 » 2 β log(2d) . (1.10) W2,ℓ∞ (πh , π) = h log(2d) O α

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Some ith variable x(i) 1st layer: N1 (i) 2nd layer: N2 (i) P Potential V (x) = di=1 Vi (x) and Vi only depends on N1 (i) Sparsity parameter sk = maxi N2k (i) .

This example: sk = O(k2 ).

Figure 1. Illustration of a potential V (x) with sparse interactions. Adapted from [17]. Compared to ourpresults for the p unadjusted Langevin algorithm in [17], we have an improvement from h log(2d) to h log(2d). The√new scaling aligns with [34], where the Wasserstein-2 bias of the BAOAB scheme is h d. It does not show a second order accuracy due to lack of higher order regularity conditions on V [39, 8]. We also note that the step-size conditions for both methods align with√existing results. For HMC, the chosen length of the integration time scales with 1/ β for both the continuous and discrete dynamics [19, 11], while we choose the specific constant to optimize the contraction rate in our analysis. For UL, we allow for the same step-size condition as in [34]. Specifically, the condition indicates that h ≤ 2√1 β . Besides, the relation √ η =√ exp(−γh) √ implies the condition γ > 2 β. In comparison, [40] establishes that γ > β − α is the sharp condition to guarantee Wasserstein contraction for general distributions, whereas for Gaussian targets, contraction holds for any γ. Therefore, in the highly ill-conditioned regime (α ≪ β), our friction requirement matches the scaling of this tight condition up to a factor of 2. This scaling also aligns with continuous-time dynamics, where the optimal Wasserstein-2 convergence √rate of O(α/γ) is attained at √ γ = β + α, a value that is asymptotically equivalent to β when α ≪ β. Notably, the step-size condition accommodates the entire high friction regime, and ensures that the bias bounds remain valid in the overdamped limit γ → ∞, where the BAOAB scheme reduces to the Leimkuhler-Matthews (LM) integrator, a method for the overdamped dynamics that achieves weak second-order accuracy [32]. Our results thus also imply the delocalization of bias effect for the LM scheme. The key technical challenge of extending our analysis in [17] is that the propagator matrices of HMC and UL are more involved. The coupling argument, a common device for obtaining Wasserstein bounds, relies on the propagator of the dynamics to induce contraction and to control the accumulation of discretization error. The lack of one-step contraction in the W2,ℓ∞ metric necessitated the use of a multi-step coupling argument in [17], where the sparsity of the potential was leveraged to control the accumulated errors over these steps. In this paper, we employ a novel matrix polynomial framework that applies to both HMC and UL, to characterize the multi-step propagator matrices

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

on the joint (q, p) space. The framework takes a global-in-time view, giving the full multi-step propagator at once as a polynomial of the Hessians along the trajectory, with coefficients we can characterize and bound. On the joint (q, p) space, position and momentum differences mix as the trajectory evolves, so the position-difference propagator is not a simple product of one-step maps but a structured polynomial pairing the two. The matrix-polynomial framework makes this structure explicit. The polynomial’s coefficients control both the contraction and the error accumulation in the ℓ∞ norm. Under sparse interactions, in particular, the ℓ∞ norm of the multi-step propagator is magnified over its ℓ2 norm by only a factor governed by the potential’s sparsity. In the easier weak interaction regime, the propagators are nearly diagonal, so contraction follows from a more elementary, near-single-step argument. For HMC this contraction is still read off from the multi-step polynomial, since one outer loop comprises several integration steps, whereas for UL it follows directly from the splitting structure of the scheme without invoking the polynomial. 1.4. Organization of this paper. In Section 2, we present our analysis framework based on the coupling argument and introduce the notation for the dynamics. Section 3 introduces a matrix polynomial framework to characterize the propagator matrices, and Section 4 establishes bounds for these propagators. Discretization error estimates are provided in Section 5. We conclude the paper in Section 6. All proofs are deferred to the appendices. 1.5. Notation. For a random vector X, we define |X|2,ℓ∞ = (E[|X|2ℓ∞ ])1/2 . We write A = O(B) or A ≲ B to denote that there exists a constant C independent of α, β, h, d such that A ≤ CB. We use | · |ℓ∞ and | · |ℓ2 to represent the ℓ∞ and ℓ2 norms for vectors and matrices. When applied to matrices, they stand for the corresponding operator norms. We use ⪯ for the Loewner order such that if M ⪯ N where M, N are symmetric matrices, then the matrix N −M is positive semi-definite. We use ⌊x⌋ and ⌈x⌉ to denote the greatest integer less than or equal to x, and the least integer greater than or equal to x, respectively. 2. Sketch of Techniques In this section, we outline our framework of the coupling argument and identify the propagator matrices as the key components for analysis. 2.1. Wasserstein bias. Our proof utilizes a coupling argument to establish both the Wasserstein convergence rate and the accumulated discretization error across iterations. These results then yield the final Wasserstein bias bounds. The following proposition outlines the Wasserstein convergence for HMC when the potential describes weak interactions. Proposition 2.1 (HMC, weak interactions). Suppose Assumptions 1.1 and 1.3 hold with C (O) = 1/50. Let √ ρk be the marginal distribution of Xk , the k-th step in the HMC chain. For mh = 1/ 20β, we have Å ã » α (2.1) W2,ℓ∞ (ρk+1 , π) ≤ 1 − W2,ℓ∞ (ρk , π) + O(h log(2d)) . 400β

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

By taking k → ∞, the iterative inequality (2.1) establishes the final bias bound in Theorem 1.4. The inequality consists of two components. The first characterizes the Wasserstein contraction rate of the HMC chain, while the second bounds the discretization error relative to an ideal HMC process that preserves π exactly. The detailed proof is in Appendix F. For the case of sparse interactions, we instead rely on a multi-step recursive inequality. Proposition 2.2 (HMC, √ sparse interactions). Suppose Assumptions 1.1 and 1.5 hold, √ log d and set ri := ⌈ih βe+ log(5/3) ⌉, a neighborhood threshold chosen to control the ℓ∞ norm √ of the propagators as in Section 4.3. When mh = 1/ 20β, we have (2.2) ! ! ãN √ Å NX m−1 p 2» α √ dW2,ℓ∞ (ρk , π)+O sri βh log(2d) . W2,ℓ∞ (ρk+N , π) ≤ 1 − 200β i=0

See Appendix F for the detailed proof. When deriving a bias bound, to preserve contraction in the first term, we take N ∼ (β/α) log d, which necessitates a multi-step coupling. The discretization error remains controlled due to the sparsity of the Hessians, √ which manifests as the sri terms in the bound. By taking k → ∞, (2.2) results in the final bias bound in Theorem 1.6. We now lay out the coupling framework underlying Propositions 2.1 and 2.2. The framework is common to HMC and UL. Once a contraction is established, over one step or several, the total bias is a sum of single-step discretization errors, each propagated forward under the discrete-time scheme. The object governing this propagation, which is the propagator matrix acting on the difference between two paths, is therefore central to the argument. Working with the discrete-time scheme is itself a deliberate choice, as the existing ℓ2 analyses in continuous time [19, 37] exploit the inner-product structure, but the ℓ∞ norm is less compatible with continuous-time differentiation. For HMC, the propagator matrices contract in the ℓ∞ norm in one outer loop under weak interactions, while the ℓ∞ contraction emerges over multiple outer loops under sparse interactions. Moreover, the same class of matrices controls how the discretization error propagates across inner steps and subsequent outer loops. The UL analysis follows a similar template, but its iterative inequalities are more involved due to a weighted norm and the absence of momentum refreshment. We therefore present the framework for both algorithms but give the iterative inequalities only for HMC. 2.2. HMC. We compare two chains driven by the same momentum refreshments ξk , . . . , ξk+N −1 . The exact chain Y has an inner loop following the continuous Hamiltonian dynamics (1.1), while the approximate chain X has an inner loop running the leap-frog scheme (1.3). We assume Y0 ∼ π is at stationarity. For each chain, one outer loop maps the current position to the next by drawing a fresh momentum ξ, integrating the inner dynamics for time mh (over m leap-frog steps in the approximate case), and reading off the resulting position. We write Φmh for one exact outer loop and Φm ξ h,ξ for its leap-frog counterpart, the subscript ξ recording the Gaussian refreshment used in that loop. After N outer loops, the approximate position is the N -fold composition along the shared

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

refreshment sequence, m m Xk+N = Φm h,ξk+N −1 Φh,ξk+N −2 · · · Φh,ξk (Xk ) ,

(2.3)

and the exact position Yk+N is given analogously by composing the Φmh maps. We ξ omit the composition symbol ◦ throughout for brevity. One-step bound. Coupling the two chains over a single outer loop, we have mh |Xk+1 − Yk+1 |2,ℓ∞ = |Φm h,ξk (Xk ) − Φξk (Yk )|2,ℓ∞ m mh m ≤ |Φm h,ξk (Xk ) − Φh,ξk (Yk )|2,ℓ∞ + |Φh,ξk (Yk ) − Φξk (Yk )|2,ℓ∞ , | {z } | {z }

(2.4)

(a)

(b)

where (a) measures how the leap-frog loop changes the distance between two paths and ∞ norm, (b) is the discretization error in one outer loop. If Φm h,ξ is a contraction in the ℓ then it suffices to bound the part (b). Weakly interacting potentials guarantee this contraction. 2 ∞ norm. Multi-step bound. In general, Φm h,ξ contracts in the ℓ norm but not in the ℓ However, an ℓ∞ contraction can be induced by the ℓ2 contraction through coupling over multiple steps, which motivated the multi-step coupling argument for the sparse interactions in [17]. The challenge persists here, where the sparsity of the potential controls the multi-step error accumulation. Iterating over N outer loops yields the same split into a contraction term (a) and a discretization term (b), mh mh m |Xk+N − Yk+N |2,ℓ∞ = |Φm h,ξk+N −1 · · · Φh,ξk (Xk ) − Φξk+N −1 · · · Φξk (Yk )|2,ℓ∞

(2.5)

m m m ≤ |Φm h,ξk+N −1 · · · Φh,ξk (Xk ) − Φh,ξk+N −1 · · · Φh,ξk (Yk )|2,ℓ∞ | {z } (a)

m mh mh + |Φm h,ξk+N −1 · · · Φh,ξk (Yk ) − Φξk+N −1 · · · Φξk (Yk )|2,ℓ∞ ,

(b)

}

where now the map in (a) becomes an ℓ∞ contraction once N ∼ (β/α) log d. Discretization error. It remains to quantify the discretization term (b), for both a single outer loop and for N of them. We begin with (b) of (2.4). To expose its structure m we need the finer building blocks of an outer loop. Let Umh hmc and Uhmc,h denote the exact and leap-frog solution maps on the full phase space R2d for time mh, and Π1 the mh m m projection onto position, so that Φmh ξ (q) = Π1 Uhmc (q, ξ) and Φh,ξ (q) = Π1 Uhmc,h (q, ξ). mh Telescoping the difference Π1 Um hmc,h − Π1 Uhmc over the m integration steps and using that Uih hmc (Yk , ξk ) ∼ π ⊗ N (0, I) is at stationarity (since Y0 ∼ π and the inner loop is exact), term (b) of (2.4) is bounded by (2.6) (b) ≤

m−1 X

1 h ∗ ∗ Π1 Um−i−1 hmc,h (Uhmc,h − Uhmc )(q , p ) 2,ℓ∞ ,

(q ∗ , p∗ ) ∼ π ⊗ N (0, I) .

i=0

Each term represents the error incurred in a single integration step, propagated over the subsequent m − i − 1 leap-frog steps. The discretization error over N outer loops, term (b) of (2.5), telescopes in the same way at the level of outer loops. It is a sum over loops l = 0, . . . , N − 1, of the oneloop error incurred at loop l, propagated through the remaining N − l − 1 outer loops.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Combined with a single-loop decomposition similar to (2.6), the total error is a sum of single integration-step errors, each propagated forward through the remaining leap-frog steps within its own loop and through all subsequent outer loops. Term (b) in (2.5) can be therefore bounded as (2.7)

(b) ≤

N −1 m−1 X X

m−i−1 1 h ∗ ∗ m |Φm h,ξk+N −1 · · · Φh,ξk+l+1 Π1 Uhmc,h (Uhmc,h − Uhmc )(q , p )|2,ℓ∞ ,

l=0 i=0

where (q ∗ , p∗ ) ∼ π ⊗ N (0, I). Here i indexes the inner integration steps and l the outer loops. Characterizing the propagator matrices that govern how a local error evolves into the position coordinate is therefore the core of the discretization-error analysis. 2.3. UL. The framework for UL follows that for HMC, with two modifications. First, the dynamics evolve without momentum refreshment, so we do not project onto the position between steps. Second, since contraction in the ℓ2 norm requires a weighted metric [39, 34], we adopt a weighted ℓ∞ norm throughout. 2.3.1. Notation and weighted norm. We use solution-map notation analogous to HMC. The one-step BAOAB map with step size h and O-step noise ξ is Uξbaoab,h , and we write Uξabao,h , Uξbao,h , and Uab,h for the ABAO, BAO, and AB sub-schemes, defined by the operator compositions (2.8)

Uξbaoab,h := Bh/2 Ah/2 Oh (·; ξ)Ah/2 Bh/2 , Uξbao,h := Oh (·; ξ)Ah/2 Bh/2 ,

Uξabao,h := Oh (·; ξ)Ah/2 Bh Ah/2 ,

Uab,h := Bh/2 Ah/2 .

We use half steps for each A and B in BAO and AB rather than the conventional full steps, since we invoke them as the components Uξbaoab,h = Uab,h ◦ Uξbao,h of the BAOAB scheme. In place of the continuous-time UL dynamics (1.2), we measure the discretization error against the HOH splitting Uξhoh,h := Hh/2 O(·; ξ) Hh/2 , where Hh := Uhmc,h is the Hamiltonian solution map. Following [34], this keeps the error estimate accurate in the high-friction regime, and HOH preserves the target π ⊗ N (0, I). We consider the weighted ℓ∞ norm ï ò ï ò I √ bI q := W , W = (2.9) |(q, p)|ℓ∞ , w p ℓ∞ 0 a − b2 I h . The weight matrix W is taken from [34], where the same where a = β1 and b = 1−η weight makes BAOAB contract in ℓ2 . Equivalence with the Euclidean norm requires b2 ≤ a/4 [39, 34], under which 2 2 2 2 2 1 3 (2.10) 2 |q|ℓ2 + a|p|ℓ2 ≤ |(q, p)|ℓ2w ≤ 2 |q|ℓ2 + a|p|ℓ2 . √ for UL. For probabilistic This condition is met under the step-size requirement h ≤ 21−η β » 2 ] and its unweighted := bounds we use the expected norms |(q, p)|2,ℓ∞ E[ |(q, p)| ℓ∞ w w counterpart | · |2,ℓ∞ . Starting from (Xk , Pk ), the BAOAB chain evolves over N steps to

(2.11)

k+N −1 k (Xk+N , Pk+N ) = Ubaoab,h · · · Uξbaoab,h (Xk , Pk ) ,

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

and the HOH chain evolves analogously from (Yk , Pk′ ) by composing the Uξhoh,h maps, where ξk , . . . , ξk+N −1 are the O-step Gaussian noises, shared by both chains to match the coupling below.

2.3.2. Analysis framework via coupling. Given a fixed k, we couple the two chains and initialize (Y0 , P0′ ) from the HOH stationary distribution π ⊗ N (0, I). Multi-step bound. As in HMC, the position gap splits into a contraction term (a) and a discretization term (b),

(2.12)

|Xk − Yk |2,ℓ∞ Ä ξ ä ξk−1 k−1 0 0 ≤ Π1 Ubaoab,h · · · Uξbaoab,h (X0 , P0 ) − Ubaoab,h · · · Uξbaoab,h (Y0 , P0′ ) 2,ℓ∞ {z } | (a)

Ä ξ ä ξk−1 k−1 0 0 . + Π1 Ubaoab,h · · · Uξbaoab,h (Y0 , P0′ ) − Uhoh,h · · · Uξhoh,h (Y0 , P0′ ) 2,ℓ∞ {z } | (b)

In both the weak and sparse interaction cases, the contraction is established over multiple steps, though for different reasons. For weak interactions, we exploit the one-step weighted contraction of the ABAO operator, writing (BAOAB)k = BAO (ABAO)k−1 AB, so that the contraction of the inner ABAO operators drives the bound while BAO and AB contribute only constant prefactors. These prefactors are what preclude a strict one-step bound and require accumulating contraction over several steps. This decomposition follows the ℓ2 argument for the BAOAB scheme in [34]. For sparse interactions, ABAO is no longer a weighted-ℓ∞ contraction, so we instead analyze the full BAOAB map and upgrade the ℓ2 contraction of [34] to an ℓ∞ contraction over multiple steps. Since we only need the asymptotic bias in the position variable, one option there is to track the propagator matrices governing error propagation into the position coordinates, whose marginal is non-Markovian across intermediate steps. In both cases, we establish the multi-step contraction from the initial state through to the final step k rather than step by step. In both regimes, (2.12) is taken in the unweighted ℓ∞ norm of the position coordinate to match the Wasserstein bias structure. The weighted contraction enters by inserting W and W −1 at the intermediate steps and invoking the equivalence in (2.10). In Section 2.1, we omit an explicit iterative Wasserstein bias inequality for UL analogous to Propositions 2.1 and 2.2, because the insertion of the weight matrices and the end-to-end bias bound complicate the statement. Discretization error. The discretization term (b) is handled through the same propagator-matrix mechanism as HMC, following the asymptotic-bias analysis of BAOAB in [34] with the weighted ℓ2 norm replaced by the ℓ∞ norm. To bound the error accurately across regimes of γ, [34] analyzes the process over blocks of ˜l = ⌈ 2√1βh ⌉ steps. Organizing

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

(b) in (2.12) block by block gives ⌊k/l̃⌋−1

(b) ≤ (2.13)

k−(i+1)l̃

Π1 Ubaoab,h (Ul̃baoab,h − Ul̃hoh,h )(q ∗ , p∗ )

i=0

k−⌊k/l̃⌋l̃ k−⌊k/l̃⌋l̃ (q ∗ , p∗ ) + Π1 Ubaoab,h − Uhoh,h

2,ℓ∞

where (q ∗ , p∗ ) ∼ π ⊗ N (0, I), and Ulbaoab,h , Ulhoh,h denote the l-fold compositions of the BAOAB and HOH solution maps, with the coupled Gaussian noises from (2.12) omitted for brevity. Each term is the error over at most ˜l steps, propagated through the subsequent BAOAB steps. Characterizing these propagator matrices for the evolution of the difference between two paths is therefore central to the accumulated discretizationerror analysis. 3. A Matrix Polynomial Framework for Propagator Analysis As shown in the previous section, the propagator matrices are fundamental to characterizing both the contraction rates and the discretization error under the ℓ∞ norm. In the absence of continuous-time differentiation, we track the evolution of path differences by introducing a framework that expresses these propagators as matrix polynomials. The underlying polynomial coefficients then facilitate the derivation of the required ℓ∞ bounds. We illustrate this structure with a Gaussian distribution. 3.1. Propagators: Gaussian examples. The propagators of the leap-frog scheme are closely related to Chebyshev polynomials, as illustrated by the following proposition. Proposition 3.1. For π = N (0, Σ), after k leap-frog steps, the difference between the position components xk = Π1 Ukhmc,h (x0 , p0 ) and yk = Π1 Ukhmc,h (y0 , p′0 ) satisfies h2 −1 h2 Σ )(x0 − y0 ) + hUk−1 (I − Σ−1 )(p0 − p′0 ) , 2 2 where Tk is the k-th Chebyshev polynomial of the first kind, and Uk is the k-th Chebyshev polynomial of the second kind. We adopt the convention U−1 = 0, which is consistent with the standard three-term recurrence relation Uk+1 (x) = 2xUk (x)−Uk−1 (x) for k ≥ 0 when initializing the sequence at U0 (x) = 1. (3.1)

xk − yk = Tk (I −

Proof sketch. Writing down the iteration (1.3) in two consecutive steps and eliminating the momentum difference yields, for ∆qk := xk − yk , the three-term recurrence (3.2)

∆qk+1 = (2I − h2 Σ−1 )∆qk − ∆qk−1 ,

k ≥ 1.

Deriving the corresponding recurrence for the propagator matrices reveals the Chebyshev polynomials. □ This is a special case of Proposition 3.5 for which we provide the complete proof. It is natural for ∆qk to satisfy a three-term recurrence, as the leap-frog scheme is a discretization of a second-order ordinary differential equation in the position variable. In fact, the position-only iteration (3.2) is the Störmer form of the leap-frog scheme [31]. The relation to the Chebyshev polynomials is also well-motivated. For instance, in a one-dimensional system with the quadratic potential V (q) = 21 q 2 , the exact solution

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

is q(t) = q(0) cos(t) + p(0) sin(t). The Chebyshev polynomials provide the discrete analogue to these trigonometric solutions via the identities Tk (cos θ) = cos(kθ) and Uk−1 (cos θ) = sin(kθ)/ sin θ. With momentum refreshment between outer loops, multiplying their matrix polynomials gives the following HMC propagator across multiple loops. Proposition 3.2. For π = N (0, Σ), let xN m+k and yN m+k denote the position components after initial k leap-frog steps, followed by N full HMC outer loops (each conm k sisting of m steps). Specifically, let xN m+k := Φm h,ξN −1 · · · Φh,ξ0 Π1 Uhmc,h (x0 , p0 ) and m k ′ yN m+k := Φm h,ξN −1 · · · Φh,ξ0 Π1 Uhmc,h (y0 , p0 ). The difference ∆qN m+k := xN m+k − yN m+k satisfies ãN Å h2 −1 ∆qN m+k = Tm (I − Σ ) 2 (3.3) ã Å h2 −1 h2 −1 ′ · Tk (I − Σ )(x0 − y0 ) + hUk−1 (I − Σ )(p0 − p0 ) . 2 2 Here, xN m+k and yN m+k are coupled by utilizing identical momentum refreshments ξ0 , . . . , ξN −1 at each outer loop. The UL dynamics are closely related to the Hamiltonian dynamics but crucially involve a friction term. The friction induces a modification to the polynomial structure. Specifically, the resulting three-term recurrence for the BAOAB scheme incorporates damped coefficients, which recovers the leap-frog coefficients by taking γ = 0. Proposition 3.3. For π = N (0, Σ), let xk and yk denote the position components ξk−1 0 after k BAOAB steps. Specifically, let xk := Π1 Ubaoab,h · · · Uξbaoab,h (x0 , p0 ) and yk := ξ

k−1 0 Π1 Ubaoab,h · · · Uξbaoab,h (y0 , p′0 ). The difference ∆qk := xk − yk then satisfies

h2 −1 Σ )∆qk − η∆qk−1 , k ≥ 1. 2 3.2. Propagators: general cases. Inspired by the Gaussian case, we establish that the propagator matrices can be represented as matrix polynomials, whose coefficients can be accurately characterized. However, the Hessian varies along the trajectory, and the resulting matrices do not commute. The polynomial must therefore be generalized to one in several ordered, non-commuting arguments, each monomial being a product taken in a fixed order. We view the general propagator as a “path-dependent” polynomial, where the recurrence remains identical to the Gaussian case but the arguments are updated in order at each step. We now define these polynomials. (3.4)

∆qk+1 = (1 + η)(I −

Definition 3.4 (Multivariate Damped Chebyshev Polynomials). For a sequence of matrices Mi ∈ Rd×d , we define the multivariate damped Chebyshev polynomial of the first kind, T̃kη (M1 , . . . , Mk ), via the three-term recurrence (3.5)

η η T̃k+1 (M1 , . . . , Mk+1 ) = (1 + η) Mk+1 T̃kη (M1 , . . . , Mk ) − η T̃k−1 (M1 , . . . , Mk−1 )

1−η for k ≥ 1, with the initial conditions T̃0η ≡ I and T̃1η (M1 ) = 1+η 2 M1 + 2 I. Note that the new argument Mk+1 enters on the left. The multivariate damped Chebyshev

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

polynomial of the second kind, Ũkη (M1 , . . . , Mk ), is defined by the same recurrence (3.6)

η η Ũk+1 (M1 , . . . , Mk+1 ) = (1 + η)Mk+1 Ũkη (M1 , . . . , Mk ) − η Ũk−1 (M1 , . . . , Mk−1 )

η for k ≥ 0, with different initial conditions Ũ−1 ≡ 0, Ũ0η ≡ 1+η 2 I. For notation brevity, we also adopt the sequence notation T̃kη ({Mi }ki=1 ) and Ũkη ({Mi }ki=1 ).

We note that when η = 1, which corresponds to the frictionless Hamiltonian regime γ = 0, setting all Mi := xI with a scalar x gives T̃kη (xI, . . . , xI) = Tk (x)I and Ũkη (xI, . . . , xI) = Uk (x)I, recovering the classical Chebyshev polynomials Tk and Uk . We also note that, even for scalar arguments, our polynomials differ from the conventional multivariate Chebyshev polynomials on Rd , which are built from tensor products of one-dimensional polynomials rather than a recurrence in ordered arguments. In the following proposition, we show that the multivariate damped Chebyshev polynomials govern the general propagator matrices in the leap-frog and BAOAB schemes. In particular, it includes the Gaussian propositions stated above. We present its proof in Appendix B. Proposition 3.5. The propagators of the leap-frog and BAOAB schemes are characη . Throughout, terized by the multivariate damped Chebyshev polynomials T̃kη and Ũk−1 let (x0 , p0 ) and (y0 , p′0 ) be two initial states evolved under synchronous coupling, sharing the same momentum refreshments in HMC and the same O-step Gaussian noise in BAOAB, and let xi , yi denote the resulting coupled positions after i steps along the trajectory specified in each item below. We set Z 1 2 h Hi = ∇2 V τ xi + (1 − τ )yi dτ . (3.7) Mi := I − 2 Hi , 0

(1) (One scheme, k steps.) For the leap-frog iterations in HMC (η = 1) and the BAOAB iterations in UL (η = e−γh ), write Ui for i steps of the scheme, namely ξi−1 ··· the leap-frog map Uihmc,h for HMC and the BAOAB composition Ubaoab,h 0 Uξbaoab,h for UL, so that xi = Π1 Ui (x0 , p0 ) and yi = Π1 Ui (y0 , p′0 ). The position difference after k steps satisfies

(3.8)

η (M1 , . . . , Mk−1 )(p0 − p′0 ) . xk − yk = T̃kη (M0 , . . . , Mk−1 )(x0 − y0 ) + h Ũk−1

(2) (HMC, initial leap-frog run plus N outer loops.) After k leap-frog steps followed by N outer loops of m leap-frog steps each, the position difference between m k m m xN m+k := Φm h,ξN −1 · · · Φh,ξ0 Π1 Uhmc,h (x0 , p0 ) and yN m+k := Φh,ξN −1 · · · Φh,ξ0 Π1 Ukhmc,h (y0 , p′0 ) satisfies (3.9)   ←− Y 1 xN m+k − yN m+k =  T̃m (Mj , Mj+1 , . . . , Mj+m−1 ) −1 j∈{k+im}N i=0

i 1 · T̃k1 (M0 , . . . , Mk−1 )(x0 − y0 ) + h Ũk−1 (M1 , . . . , Mk−1 )(p0 − p′0 ) , where xj , yj are the coupled positions after j total steps (the first k within Qthe initial leap-frog run, the remainder within subsequent outer loops), and ←− orders the factors so that those with smaller index j stand to the right.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

4. Bounds for Propagators In this section, we present ℓ2 and ℓ∞ bounds for propagators, instantiating the technical framework sketched in Section 2. While the ℓ2 bounds are obtained by adapting existing techniques to our setting, the ℓ∞ bounds represent our primary analytical contribution. These estimates hinge on the structural properties of the multivariate damped Chebyshev polynomials, which allow us to bridge from ℓ2 stability to ℓ∞ guarantees under the assumption of weak or sparse interactions. 4.1. ℓ2 bounds. The relevant ℓ2 bounds for HMC are as follows. Proposition 4.1 (ℓ2 bounds for HMC propagators). Suppose Assumption 1.1 hold. If R 2 Hi = ∇ V dνi for a probability measure νi on Rd , we have the following ℓ2 bounds for √ the leap-frog propagators T̃k and Ũk−1 provided kh ≤ 1/ 20β: (4.1)

|T̃k1 ({I −

h2 α Hi }k−1 (kh)2 , i=0 )|ℓ2 ≤ 1 − 2 10

1 ({I − |hŨk−1

h2 Hi }k−1 i=1 )|ℓ2 ≤ 2

2 . 5β

The proof is provided in Appendix C. Combined with Proposition 3.5, Proposition 4.1 states that the leap-frog scheme induces a contraction in the position components, provided the trajectories are initialized with identical momentum and the total integration √ time remains within a specific threshold scaled as 1/ β. The momentum difference may √ expand, though this growth remains well-controlled by a magnitude scaled as 1/ β. These scalings are consistent with the Gaussian setting under the continuous-time dynamics within the time range. We adapt the proof of an existing discrete-time contraction result [11, Lemma 19] to obtain the bound for T̃k1 . Applying a similar framework for the momentum components, we establish the Ũk1 bound. For UL, contraction is induced in the weighted norms of Section 2, so we build the weight matrix W into the propagator representation. For two coupled trajectories with initial states z0 = (x0 , p0 ) and z0′ = (y0 , p′0 ), viewed as column vectors in R2d , a one-step update U (any of BAOAB, ABAO, BAO, AB) gives z1 = U(z0 ) and z1′ = U(z0′ ). The resulting state difference is z1 − z1′ = M (H0 )(z0 − z0′ ), where M (H0 ) is the one-step propagator of the unweighted difference, parametrized by an integrated Hessian H0 . Inserting I = W −1 W gives (4.2) W (z1 − z1′ ) = W M (H0 )W −1 W (z0 − z0′ ) . The weighted propagator M w (H0 ) := W M (H0 )W −1 acts on the weighted difference W (z0 − z0′ ), so a weighted ℓ2 contraction follows from bounding |M w (H0 )|ℓ2 , and discretization errors are measured in the weighted norm |z0 − z0′ |ℓ2w = |W (z0 − z0′ )|ℓ2 . Over k steps the propagator is the composition Mk := M (Hk−1 ) · · · M (H0 ), with weighted counterpart Mkw := W Mk W −1 , and a k-step weighted contraction follows from bounding the full-space operator norm |Mkw |ℓ2 . Our goal, however, is the bias of the unweighted position marginal, recovered by the projection Π1 = I 0 ∈ Rd×2d . After one step, we have (4.3) Π1 (z1 − z1′ ) = I 0 W −1 M w (H0 ) W (z0 − z0′ ) . The factor I 0 M (H0 ) captures the dependence of x1 − y1 on the initial position and η momentum differences. The same structure holds over k steps, with T̃kη and hŨk−1 as

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

the position-row blocks of Mk . We bound both the full-space propagator and these weighted blocks for the BAOAB scheme as follows. Proposition 4.2 (ℓ2 bounds for UL propagators). Suppose Hi satisfies the same as2 √ , and letting c(h) = αh , sumptions as in Proposition 4.1. For a step size h ≤ 21−η 4(1−η) β the full-space weighted BAOAB propagator satisfies (4.4)

k−1

w |Mbaoab,k |ℓ2 ≤ 7(1 − c(h)) 2 ,

and its weighted position-row blocks satisfy (4.5)

|T̃kη ({I −

√ k−1 h2 2 , Hi }k−1 i=0 )|ℓ2 ≤ 7 2 (1 − c(h)) 2

(4.6) √ k−1 1 h2 h2 η −bT̃kη ({I − Hi }k−1 ) + hŨk−1 ({I − Hi }k−1 ) ≤ 7 2 (1 − c(h)) 2 , i=0 i=1 2 2 a − b2 ℓ2 The proof is detailed in Appendix D. Proposition 4.2 is an adaptation of the existing contraction result [34, Theorem 5.1]. The specific combinations of matrices are of interest because we have ó î η ó î η η η η b √ 1 + T̃ h Ũ (4.7) W −1 = T̃k − √a−b T̃k hŨk−1 2 k k−1 . a−b2 √

4.2. ℓ∞ bounds under weak interactions. The weak interactions defined in Assumption 1.3 indicate that ∇2 V is a small perturbation of a diagonal matrix. The resulting propagators therefore have ℓ∞ norms close to their ℓ2 norms, since the two norms coincide for a diagonal matrix, and a contraction is preserved in ℓ∞ whenever the interactions are weak. For HMC, since one outer loop is m leap-frog steps, the relevant object is the multi-step propagator T̃k1 . Proposition 4.3 (ℓ∞ bounds under weak interactions, HMC). Suppose Assumptions 1.1 and 1.3 hold with the weak interaction parameter C (O)√= 1/50. Let the integrated Hessians Hi be defined as in Proposition 4.1. For kh ≤ 1/ 20β, we have (4.8) h2 α khα h2 3 1 √ , |hŨk−1 |T̃k1 ({I − Hi }k−1 (kh)2 + ({I − Hi }k−1 i=0 )|ℓ∞ ≤ 1 − i=1 )|ℓ∞ ≤ √ . 2 10 2 100 β 2 β 2

1 α In particular, for kh = √20β , we have the contraction |T̃k1 ({I − h2 Hi }k−1 i=0 )|ℓ∞ ≤ 1− 400β .

The proof can be found in Appendix C. For UL, by contrast, the contraction is obtained at the level of a single step rather than through a multi-step polynomial. We decompose the full BAOAB operator into BAO, ABAO, and AB, and establish a one-step contraction for the ABAO operator. w As in the general construction (4.2), Mabao (H0 ) is the weighted one-step propagator ξ ξ w of the Uabao,h map, satisfying W (Uabao,h (z0 ) − Uξabao,h (z0′ )) = Mabao (H0 ) W (z0 − z0′ ) ′ ′ for states z0 = (x0 , p0 ), z0 = (y0 , p0 ) coupled through the same O-step noise. Here R1 H0 = 0 ∇2 V τ (x0 + h2 p0 ) + (1 − τ )(y0 + h2 p′0 ) dτ , and   1 + η 2 4 η (1−η) h2 2 3 h √ H0  I − 2 H0 − 1−η h H0 w 1/β−h2 /(1−η)2 (4.9) Mabao (H0 ) =  q  . h2 1 −ηh β1 − (1−η) η I + (− 12 + 1−η )h2 H0 2 H0

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

We work with this full phase-space propagator rather than its position projection because the ABAO operator drives contraction in the joint space over several iterations w and M w are defined analobefore the position marginal is extracted. The maps Mbao ab gously by their one-step propagation of the weighted state differences, though the state combination over which H0 integrates may differ between schemes. Proposition 4.4 (ℓ∞ bounds under weak interactions, UL). Suppose Assumptions 1.1 and 1.3 hold with C (O) = 1/20. Let the integrated Hessian H0 be defined as in 2 √ and recalling c(h) = αh Proposition 4.1. For h ≤ 21−η 4(1−η) from Proposition 4.2, we have β 1 w |Mabao (H0 )|ℓ∞ ≤ 1 − c(h) , 2

(4.10)

w |Mbao (H0 )|ℓ∞ ≤ 3 ,

w |Mab (H0 )|ℓ∞ ≤ 2 .

These bounds are the ℓ∞ analogues of the ℓ2 ingredients used to prove BAOAB contraction in [34]. See Appendix D for its proof. 4.3. ℓ∞ bounds under sparse interactions. We now turn to the propagator bounds under sparse interactions. Recall from Assumption 1.5 that each term Vi depends only on the coordinates in the neighborhood N(i). A single Hessian ∇2 V therefore has a nonzero (i, j) entry only when some term Vl involves both i and j, that is, i, j ∈ N(l), so i and j are within two edges in the interaction graph. More generally, Lemma C.2 in Appendix C shows that a product of r Hessians, the basic building block of the matrix polynomials, has its (i, j) entry supported on j ∈ N2r (i), with at most sr nonzero entries per row. This locality is what yields the ℓ∞ bounds. To bound a row of a matrix polynomial, we fix a degree threshold r and split the entries by neighborhood. For j ∈ N2r (i), there are at most sr such indices, and a Cauchy–Schwarz estimate over them bounds their √ / N2r (i), the entry only arises from prodcontribution by sr times the ℓ2 norm. For j ∈ ucts of degree above r, which carry high powers of the step size and decay geometrically. Bounding sr together with this decay converts the ℓ2 propagator estimates into ℓ∞ ones and controls error accumulation. The following proposition is proved in Appendix C. Proposition 4.5 (ℓ∞ bounds under sparse interactions, HMC). Suppose Assumptions 1.1 and 1.5 hold. Let Hi be defined as in Proposition 4.1. Take √ the integrated Hessians √ √ log d ri = ⌈ih βe + log(5/3) ⌉. For mh ≤ 1/ 20β, integers N ≥ 0 and k < m, we have (4.11)

←− Y

1 T̃m ({I −

−1 j∈{k+im}N i=0

(4.12)

←− Y

−1 j∈{k+im}N i=0

1 T̃m ({I −

h2 h2 √ Hj+l }m−1 ) T̃k1 ({I − Hi }k−1 i=0 ) ∞ ≤ 4 srN m+k , l=0 2 2 ℓ

h2 h2 4 √ 1 Hj+l }m−1 ) hŨk−1 ({I − Hi }k−1 srN m+k . i=1 ) ∞ ≤ √ l=0 2 2 ℓ β

For the BAOAB scheme, the weighted matrices relevant to the position projection have the following ℓ∞ bounds. The proof is provided in Appendix D. Proposition 4.6 (ℓ∞ bounds under sparse interactions, UL). Suppose Assumptions 1.1 and 1.5 hold. Let the integrated Hessians Hi be defined as in Proposition 4.1, and

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

√ 2β 1−η αh2 2 √ let ri = ⌈ ie 1−η h + log( d)⌉. Recall that c(h) = 4(1−η) . For h ≤ 2 β , we have (4.13)

|T̃kη ({I −

√ h2 k−1 √ Hi }k−1 c(h)) , i=0 )|ℓ∞ ≤ (7 2 + 1) srk exp(− 2 2

h2 1 h2 η −bT̃kη ({I − Hi }k−1 Hi }k−1 ) + h Ũ ({I − i=0 i=1 ) k−1 2 2 a − b2 ℓ∞ √ k−1 2 √ c(h)) . ≤(7 2 + √ ) srk exp(− 2 3 √

(4.14)

A contraction factor decaying exponentially in N could also be included in (4.11) of Proposition 4.5, as in Proposition 4.6. We omit it, as our analysis does not require it. 5. Discretization Error and Sampling Bias Besides the propagator bounds of Section 4, the other critical component of the coupling framework in Section 2 is the discretization error. In this section we bound it for HMC and UL, then combine it with the propagator bounds to obtain the sampling bias. 5.1. Discretization error in HMC. To analyze the position and momentum errors separately, we split each solution map into its position and momentum components, writing Uthmc = (Qt , Pt ) for the continuous Hamiltonian dynamics (1.1), and Um hmc,h = m m (Qh , Ph ) for m leap-frog steps (1.3). The single-step discretization error then splits as (5.1) (U1hmc,h − Uhhmc )(q ∗ , p∗ ) = (Q1h − Qh )(q ∗ , p∗ ), (P1h − Ph )(q ∗ , p∗ ) , which is the common building block of (2.6) and (2.7). The following proposition bounds each component. We provide the proof to a generalized version in Appendix E. Proposition 5.1. SupposeR for any probability measures ν1 , ν2 on Rd , it holds that R 2 R 2 2 | ∇ V dν1 |ℓ∞ < C3 and |( ∇ V dν1 )( √∇ V dν2 )|ℓ∞ < C4 for constants C3 , C4 > 0 that may depend on β. Then for h ≤ 1/ 20β and (q ∗ , p∗ ) ∼ π ⊗ N (0, I), we have » 1 (5.2) |Q1h (q ∗ , p∗ ) − Qh (q ∗ , p∗ )|2,ℓ∞ ≤ √ h3 C3 log(2d) , 2 » » 1 (5.3) |P1h (q ∗ , p∗ ) − Ph (q ∗ , p∗ )|2,ℓ∞ ≤ 2h2 C3 log(2d) + √ h4 C4 log(2d) . 2 5 5.2. Discretization error in UL. For the BAOAB integrator we adapt the framework of [34], comparing the BAOAB chain with an HOH chain for an accurate discretization error estimate even in the large-γ limit. Our analysis differs mainly in measuring the error in the ℓ∞ norm rather than ℓ2 . The key component of (2.13) is (Ulbaoab,h − Ulhoh,h )(q ∗ , p∗ ). Run l steps of each scheme from the same initial state (q ∗ , p∗ ) under the same O-step noise ξ0 , . . . , ξl−1 : (5.4)

i−1 (xi , pi ) = Uhoh,h (xi−1 , pi−1 ) ,

i−1 (yi , p′i ) = Ubaoab,h (yi−1 , p′i−1 ) ,

with (x0 , p0 ) = (y0 , p′0 ) = (q ∗ , p∗ ). Writing ∆iq = xi −yi and ∆ip = pi −p′i for the accumulated discrepancies, the difference above is exactly (∆lq , ∆lp ). The following proposition bounds it in the weighted norm. See Appendix E for the proof.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Proposition 5.2. Assume the shared initial state is (q ∗ , p∗ ) ∼ π ⊗ N (0, I) and supd pose R 2Assumption 1.1 holds. R If 2for any measures ν1 , ν2 on R , we have Rprobability 2 ∇ V dν2 |ℓ∞ ≤ C6 , and the step size sat| ∇ V dν1 |ℓ∞ ≤ C5 and | ∇ V dν1 1−η isfies h ≤ 2√β , then the accumulated discretization error across l steps can be bounded by ï p ò 3 C5 1 C6 l l ∞ + β) El , (5.5) |(∆q , ∆p )|2,ℓw ≤ exp (l − 1)h( √ + 2 β 16 β 3/2 where ã» Å 3 h4 C6 27 h2 C5 27 √ + √ C 5 h3 l + log(2d) . (5.6) El = 64 8 β(1 − η) 8 β(1 − η) 5.3. Assembling the bias bounds. Combining the discretization-error estimates above with the ℓ∞ propagator bounds of Section 4, under the coupling framework of Section 2, yields iterative inequalities for HMC and its underdamped counterpart. Iterating these to the stationary limit gives the bias bounds of Theorems 1.4 and 1.6. The full assembly is carried out in Appendix F. 6. Conclusions and Discussions We study the delocalization of bias phenomenon in unadjusted Hamiltonian Monte Carlo and underdamped Langevin algorithms. Although these algorithms are biased at finite step size h, we show that the W2 bias of low-dimensional marginals scales with the marginal dimension rather than the full dimension d, when the √ target potential has weak or sparse interactions, or is Gaussian. Consequently, O( K) integration steps suffice for a bounded W2 error across all K-marginals, up to logarithmic dependence on d. This offers a counterpoint to common practice in Bayesian statistics, where Metropolis–Hastings adjustment is incorporated to eliminate discretization bias. The adjusted variants instead require a step size scaling with the full dimension d, leading to a large number of iterations. More broadly, we introduce a general matrix-polynomial framework for the propagators of discretized Hamiltonian and Langevin dynamics. The polynomial representation is well-suited to the ℓ∞ norm, which is not amenable to continuous-time tools such as differentiation. The framework is also not specific to the leap-frog and BAOAB schemes studied here. It applies to other splitting integrators, including OBABO and ABOBA, where, possibly after regrouping the iterates, a single repeating core block satisfies a similar three-term recurrence and the boundary factors contribute only bounded prefactors that do not affect the coefficient estimates. We thus expect the framework to facilitate the study of delocalization of bias in other dynamics as well. Another interesting direction is randomized HMC [10, 8], where the integration time is randomized to avoid periodicity. Our proof strategy does not directly apply there, as the number of leap-frog steps per outer loop is random rather than fixed. We follow assumptions on the potential similar to those in [17, 30], and relaxing them toward conditions met by real physical systems is a significant future direction. One such setting is non-log-concave target distributions. Another is more general interaction patterns, such as the decaying but long-range interactions of physical systems rather than the sparse interactions considered here. These extensions are left for future work.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Acknowledgments We thank Nawaf Bou-Rabee, Aaron Dinner, Daniel Lacker, Peter A. Whalley, and Fuzhong Zhou for helpful discussions. Y. Chen acknowledges the support from the UCLA DataX Pilot Projects Grant Program. J. Niles-Weed is supported in part by National Science Foundation award DMS-2339829. J. Weare and X. Cheng are supported in part by National Science Foundation award DMS-2425899. X. Cheng is also supported in part by a Dean’s Dissertation Fellowship from the Graduate School of Arts and Science at New York University. AI tools were used only to edit language, flag typographical errors and minor inconsistencies, and generate the TikZ code for the figure.

References [1] Mark James Abraham, Teemu Murtola, Roland Schulz, Szilárd Páll, Jeremy C Smith, Berk Hess, and Erik Lindahl. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX, 1:19–25, 2015. [2] Milton Abramowitz and Irene A. Stegun. Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables. Dover, 1964. [3] Michael P. Allen and Dominic J. Tildesley. Computer Simulation of Liquids. Oxford University Press, 2017. [4] Hans C Andersen. Molecular dynamics simulations at constant pressure and/or temperature. The Journal of chemical physics, 72(4):2384–2393, 1980. [5] Simon Apers, Sander Gribling, and Dániel Szilágyi. Hamiltonian Monte Carlo for efficient Gaussian sampling: Long and random steps. Journal of Machine Learning Research, 25(348):1–30, 2024. [6] Alexandros Beskos, Natesh Pillai, Gareth Roberts, Jesus-Maria Sanz-Serna, and Andrew Stuart. Optimal tuning of the Hybrid Monte-Carlo Algorithm. Bernoulli, 19(5A):1501–1534, 2013. [7] Nawaf Bou-Rabee and Andreas Eberle. Mixing time guarantees for unadjusted Hamiltonian Monte Carlo. Bernoulli, 29(1):75–104, 2023. [8] Nawaf Bou-Rabee and Milo Marsden. Unadjusted Hamiltonian MCMC with stratified Monte Carlo time integration. The Annals of Applied Probability, 35(1):360–392, 2025. [9] Nawaf Bou-Rabee, Siddharth Mitra, and Andre Wibisono. Tail-sensitive KL and Rényi convergence of unadjusted Hamiltonian Monte Carlo via one-shot couplings. arXiv preprint arXiv:2601.09019, 2026. [10] Nawaf Bou-Rabee and Jesús Marı́a Sanz-Serna. Randomized Hamiltonian Monte Carlo. The Annals of Applied Probability, 27(4):2159–2194, 2017. [11] Nawaf Bou-Rabee and Katharina Schuh. Convergence of unadjusted Hamiltonian Monte Carlo for mean-field models. Electronic Journal of Probability, 28:1–40, 2023. [12] Evan Camrud, Alain Durmus, Pierre Monmarché, and Gabriel Stoltz. Second order quantitative bounds for unadjusted generalized Hamiltonian Monte Carlo. arXiv preprint arXiv:2306.09513, 2023. [13] Bob Carpenter, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus A. Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of Statistical Software, 76(1):1–32, 2017. [14] David A Case, David S Cerutti, Vinı́cius Wilian D Cruzeiro, Thomas A Darden, Robert E Duke, Mahdieh Ghazimirsaeed, George M Giambasu, Timothy J Giese, Andreas W Gotz, Julie A Harris, et al. Recent developments in Amber biomolecular simulations. Journal of Chemical Information and Modeling, 65(15):7835–7843, 2025. [15] Neil K. Chada, Benedict Leimkuhler, Daniel Paulin, and Peter A. Whalley. Unbiased kinetic Langevin Monte Carlo with inexact gradients. arXiv preprint arXiv:2311.05025, 2023. [16] Martin Chak and Pierre Monmarché. Reflection coupling for unadjusted generalized Hamiltonian Monte Carlo in the nonconvex stochastic gradient case. IMA Journal of Numerical Analysis, page draf045, 2025.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

[17] Yifan Chen, Xiaoou Cheng, Jonathan Niles-Weed, and Jonathan Weare. Convergence of unadjusted langevin in high dimensions: Delocalization of bias. Communications on Pure and Applied Mathematics, 79(6):1467–1491, 2026. [18] Yuansi Chen and Khashayar Gatmiry. When does Metropolized Hamiltonian Monte Carlo provably outperform Metropolis-adjusted Langevin algorithm? arXiv preprint arXiv:2304.04724, 2023. [19] 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. [20] Xiang Cheng, Niladri S. Chatterji, Peter L. Bartlett, and Michael I. Jordan. Underdamped Langevin MCMC: A non-asymptotic analysis. In Proceedings of the 31st Conference on Learning Theory (COLT), 2018. [21] Tiangang Cui, Shuigen Liu, and Xin T. Tong. Stein’s method for marginals on large graphical models. arXiv preprint arXiv:2410.11771, 2024. [22] Arnak S. Dalalyan and Lionel Riou-Durand. On sampling from a log-concave density using kinetic Langevin diffusions. Bernoulli, 26(3):1956 – 1988, 2020. [23] Simon Duane, Anthony D. Kennedy, Brian J. Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics letters B, 195(2):216–222, 1987. [24] Alain Durmus and Andreas Eberle. Asymptotic bias of inexact Markov chain Monte Carlo methods in high dimension. The Annals of Applied Probability, 34(4):3435–3468, 2024. [25] Peter Eastman, Jason Swails, John D. Chodera, Robert T. McGibbon, Yutong Zhao, Kyle A. Beauchamp, Lee-Ping Wang, Andrew C. Simmonett, Matthew P. Harrigan, Chaya D. Stern, Rafal P. Wiewiora, Bernard R. Brooks, and Vijay S. Pande. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. PLOS Computational Biology, 13(7):e1005659, 2017. [26] Andreas Eberle, Arnaud Guillin, and Raphael Zimmer. Couplings and quantitative contraction rates for Langevin dynamics. The Annals of Probability, 47(4):1982–2010, 2019. [27] Daan Frenkel and Berend Smit. Understanding Molecular Simulation: From Algorithms to Applications. Elsevier, 2023. [28] Nicolaı̈ Gouraud, Pierre Le Bris, Adrien Majka, and Pierre Monmarché. HMC and underdamped Langevin united in the unadjusted convex smooth case. SIAM/ASA Journal on Uncertainty Quantification, 13(1):278–303, 2025. [29] Stefanie Kieninger and Bettina G Keller. GROMACS stochastic dynamics and BAOAB are equivalent configurational sampling algorithms. Journal of Chemical Theory and Computation, 18(10):5792–5798, 2022. [30] Daniel Lacker and Fuzhong Zhou. A hierarchical entropy method for the delocalization of bias in high-dimensional Langevin Monte Carlo. arXiv preprint arXiv:2509.08619, 2025. [31] Ben Leimkuhler and Charles Matthews. Molecular Dynamics: With Deterministic and Stochastic Numerical Methods. Springer, 2015. [32] Benedict Leimkuhler and Charles Matthews. Rational construction of stochastic numerical methods for molecular sampling. Applied Mathematics Research eXpress, 2013(1):34–56, 2013. [33] Benedict Leimkuhler and Charles Matthews. Robust and efficient configurational molecular sampling via Langevin dynamics. The Journal of Chemical Physics, 138(17):174102, 2013. [34] Benedict J. Leimkuhler, Daniel Paulin, and Peter A. Whalley. Contraction and convergence rates for discretized kinetic Langevin dynamics. SIAM Journal on Numerical Analysis, 62(3):1226–1258, 2024. [35] Yi-An Ma, Niladri Chatterji, Xiang Cheng, Nicolas Flammarion, Peter Bartlett, and Michael I. Jordan. Is there an analog of Nesterov acceleration for MCMC? Bernoulli, 27(3):1942–1992, 2021. [36] Oren Mangoubi and Aaron Smith. Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions 2: Numerical integrators. In Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics (AISTATS), 2019. [37] 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. [38] Oren Mangoubi and Nisheeth Vishnoi. Dimensionally tight bounds for second-order Hamiltonian Monte Carlo. In Advances in Neural Information Processing Systems 31 (NeurIPS), 2018.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

[39] Pierre Monmarché. High-dimensional MCMC with a standard splitting scheme for the underdamped Langevin diffusion. Electronic Journal of Statistics, 15(2):4117–4166, 2021. [40] Pierre Monmarché. Almost sure contraction for diffusions on Rd . Application to generalised Langevin diffusions. Stochastic Processes and their Applications, 161:316–349, 2023. [41] Pierre Monmarché. An entropic approach for Hamiltonian Monte Carlo: The idealized case. The Annals of Applied Probability, 34(2):2243–2293, 2024. [42] Radford M. Neal. MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo, pages 113–162. Chapman and Hall/CRC, 2011. [43] James C Phillips, David J Hardy, Julio DC Maia, John E Stone, João V Ribeiro, Rafael C Bernardi, Ronak Buch, Giacomo Fiorin, Jérôme Hénin, Wei Jiang, et al. Scalable molecular dynamics on CPU and GPU architectures with NAMD. The Journal of Chemical Physics, 153(4):044130, 2020. [44] Joshua A Rackers, Zhi Wang, Chao Lu, Marie L Laury, Louis Lagardère, Michael J Schnieders, Jean-Philip Piquemal, Pengyu Ren, and Jay W Ponder. Tinker 8: Software tools for molecular design. Journal of Chemical Theory and Computation, 14(10):5273–5289, 2018. [45] Lionel Riou-Durand and Jure Vogrinc. Metropolis adjusted Langevin trajectories: a robust alternative to Hamiltonian Monte Carlo. arXiv preprint arXiv:2202.13230, 2022. [46] Jakob Robnik, Reuben Cohn-Gordon, and Uroš Seljak. Black-box unadjusted Hamiltonian Monte Carlo. arXiv preprint arXiv:2412.08876, 2024. [47] John Salvatier, Thomas V. Wiecki, and Christopher Fonnesbeck. Probabilistic programming in Python using PyMC3. PeerJ Computer Science, 2:e55, 2016. [48] Katharina Schuh and Peter A. Whalley. Convergence of kinetic Langevin samplers for non-convex potentials. arXiv preprint arXiv:2405.09992, 2024. [49] Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. In Advances in Neural Information Processing Systems 32 (NeurIPS), 2019. [50] Aidan P Thompson, H Metin Aktulga, Richard Berger, Dan S Bolintineanu, W Michael Brown, Paul S Crozier, Pieter J In’t Veld, Axel Kohlmeyer, Stan G Moore, Trung Dac Nguyen, et al. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer Physics Communications, 271:108171, 2022. [51] Matthew S. Zhang, Jason M. Altschuler, and Sinho Chewi. Algorithmic warm starts for Hamiltonian Monte Carlo. arXiv preprint arXiv:2603.22741, 2026.

The appendices contain all technical details. Section A gives the proof for Gaussian distributions. Section B derives the matrix polynomials. Sections C and D provide bounds for the propagators of the leap-frog and BAOAB integrators in ℓ2 and ℓ∞ . Discretization error analysis is laid out in Section E. The sampling bias bounds are concluded in Section F.

Appendix A. Proof for Gaussian Distributions Proof of Example 1.2. We assume m is chosen so that mh avoids the discrete set of resonant values at which some eigenmode fails to contract. Away from these, the leapfrog iteration contracts in every mode and the chain converges to πh . As these values are exceptional, we omit the condition from the example. According to [5], the invariant distribution for a target distribution π = N (µ, Σ) is given by πh = N (µ, Σ∞ ), where h2 −1 −1 1/2 Z + µ Σ−1 ∞ = Σ (I − 4 Σ ). Consider the coupling between πh and π as X = Σ 1/2

and Y = Σ∞ Z + µ with Z ∼ N (0, I). Then, (A.1)

2 2 1/2 2 − Σ1/2 W2,ℓ ∞ (πh , π) ≤ E[|X − Y |ℓ∞ ] = E[|(Σ ∞ )Z|ℓ∞ ].

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN 1/2

1/2

Entries of (Σ1/2 − Σ∞ )Z are |Σ1/2 − Σ∞ |2ℓ2 -subgaussian. And p 1 1 |Σ1/2 − Σ1/2 − = O( βh2 ) . ∞ |ℓ2 = max 1≤i≤d ωi ω̂i p √ 2p By√Lemma E.1, we have W2,ℓ∞ (πh , π) = O( βh log(2d)) = O(h log(2d)) as h < 1/ β. □ Appendix B. Derivation of Matrix Polynomials Throughout the appendices, when a multivariate damped Chebyshev polynomial is evaluated at the canonical argument determined by its index, with Hi as specified in the relevant statement, we suppress the argument and write 2 2 η η := Ũk−1 T̃kη := T̃kη {I − h2 Hi }k−1 Ũk−1 {I − h2 Hi }k−1 i=0 , i=1 . Arguments are written explicitly whenever they differ from this convention, as in the shifted windows of the outer-loop blocks and in the scalar reductions below. Proof of Proposition 3.5. Part (a). We derive the correspondence with multivariate damped matrix polynomials from three-term recurrence relations, similar to Proposition 3.1 but with varying Hessians. With the same noise applied to both trajectories, the BAOAB iteration (1.4) for the differences ∆qk := xk − yk and ∆pk := pk − p′k is 2

(B.1)

∆qk+1 = ∆qk + h2 (1 + η)∆pk − h4 (1 + η)Hk ∆qk , ∆pk+1 = η ∆pk − h2 Hk ∆qk − h2 Hk+1 ∆qk+1 .

The leap-frog iteration (1.3) is the special case η = 1, so it suffices to treat (B.1) for general η. We write down (B.1) also for (∆qk , ∆pk ) and eliminate the momentum ∆pk in two consecutive steps by matching its coefficients. We then obtain (B.2)

∆qk+1 = (1 + η)(I − h2 Hk )∆qk − η∆qk−1 ,

k ≥ 1.

At η = 1, the recurrence recovers the relation for the leap-frog scheme, ∆qk+1 = (2I − h2 Hk )∆qk − ∆qk−1 . We prove (3.8) by induction. It holds for k = 0, 1 by the initial conditions of T̃kη and η Ũk in Definition 3.4. For k ≥ 1, the induction hypothesis and (B.2) give 2 η ∆qk+1 = (1 + η)(I − h2 Hk )T̃kη − η T̃k−1 ∆q0 (B.3) 2 η η + (1 + η)(I − h2 Hk )hŨk−1 − ηhŨk−2 ∆p0 , and the three-term recurrences for T̃kη and Ũkη in (3.5) and (3.6) identify the right-hand side as (3.8), completing the induction. Part (b). We iterate the single-scheme identity at η = 1 across the HMC outer loops. For N ≥ 1, since both trajectories share the refreshment ξN −1 , the identity gives 2 m+k−1 1 ∆qN m+k = T̃m {I − h2 Hi }N i=(N −1)m+k ∆q(N −1)m+k (B.4)

←− Y

2 1 T̃m {I − h2 Hl }j+m−1 ∆qk , l=j

−1 j∈{k+im}N i=0

and the conclusion follows by substituting the expression for ∆qk from (3.8).

□

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

A property of these matrix polynomials, which we use when bounding the propagators, is that their coefficients have alternating signs depending only on the total order of each term. Lemma B.1. For any η ∈ [0, 1], the multivariate damped Chebyshev matrix polynomials T̃kη and Ũkη expand as (B.5)

(k) T̃kη = t0 I +

k X

(B.6)

k X

t̃i1 ,...,ij Hi1 · · · Hij ,

k−1≥i1 >···>ij ≥0

j=1

(k) Ũkη = u0 I +

(k)

(−1)j h2j

(k)

(−1)j h2j

ũi1 ,...,ij Hi1 · · · Hij ,

j=1

k≥i1 >···>ij ≥1

t0 = 1 ,

u0 = 12 (1 + η)

η with Ũ−1 ≡ 0. For k ≥ 0, (k)

(B.7)

(k)

η k+1 − 1 , η−1

and, for k ≥ i1 + 1, (B.8)

(k)

t̃i1 ,...,ij =

η k−i1 − 1 (i1 +1) t̃ , η − 1 i1 ,...,ij

(k−1)

ũi1 ,...ij =

η k−i1 − 1 (i1 ) ũi1 ,...,ij , η−1

where for η = 1 the quotients are read as the limit limη→1 (η k − 1)/(η − 1) = k. All of these coefficients are nonnegative: (B.9)

(k)

t0 , u0 , t̃i1 ,...,ij , ũi1 ,...,ij ≥ 0 .

Proof of Lemma B.1. The definitions of T̃kη and Ũkη in (3.5) and (3.6) give the matrix recurrences 2

(B.10)

η η , T̃k+1 = (1 + η)(I − h2 Hk )T̃kη − η T̃k−1

(B.11)

η η = (1 + η)(I − h2 Hk+1 )Ũkη − η Ũk−1 Ũk+1 .

We establish the expansions (B.5) and (B.6) by induction, the cases k = 0, 1 following from the initial conditions. (k+1) (k) (k−1) Constant terms. The constant coefficients satisfy t0 = (1 + η)t0 − ηt0 and (k+1) (k) (k−1) (0) (1) (−1) (0) u0 = (1 + η)u0 − ηu0 , with t0 = t0 = 1, u0 = 0, and u0 = (1 + η)/2. (k) (k) η k+1 −1 1 Solving these gives t0 = 1 and u0 = 2 (1 + η) η−1 , the latter read as k + 1 when η = 1. Higher-order terms. Assume (B.5) and (B.6) hold up to index k. Passing to k + 1 2 via (B.10) and (B.11), the only new matrix products arise from −(1 + η) h2 Hk T̃kη and 2 −(1 + η) h2 Hk+1 Ũkη , which prepend Hk (resp. Hk+1 ) to the existing terms. Applied to the order-j terms of T̃kη and Ũkη , this produces the order-(j + 1) contributions X (k) (B.12) (−1)j+1 h2(j+1) 1+η t̃i1 ,...,ij Hk Hi1 · · · Hij , 2 k−1≥i1 >···>ij ≥0

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN (k)

(−1)j+1 h2(j+1) 1+η 2

(B.13)

ũi1 ,...,ij Hk+1 Hi1 · · · Hij ,

k≥i1 >···>ij ≥1 2

(k)

together with the new order-one terms −(1+η) h2 t0 Hk and −(1+η) h2 u0 Hk+1 . In each new product the prepended factor carries a strictly larger index than all of i1 > · · · > ij , so the matrices remain in descending index order, and (B.5) and (B.6) hold at k + 1. Nonnegativity. We match matrix products across the recurrences and induct on j, the number of factors. The case j = 0 was settled above. A product Hi1 · · · Hij first appears in T̃iη1 +1 and Ũiη1 , with coefficients, for i1 ≥ 1, (i +1)

(i )

1 t̃i11,...,ij = 1+η 2 t̃i2 ,...,ij ,

(B.14)

(i −1)

(i )

1 ũi11,...,ij = 1+η 2 ũi2 ,...,ij ,

(i )

where for j = 1 the right-hand factors are read as the constant coefficients t0 1 and (i −1) (1) u0 1 . Besides the induction, we also note that the first-kind base case is t̃0 = 1+η 4 ≥ η 2 0, the coefficient of −h H0 in T̃1 . By the induction hypothesis on j − 1, these initial coefficients are nonnegative. Tracking the coefficient of this same product Hi1 · · · Hij in the higher-index polynomials, it satisfies (i +2)

(i +1)

t̃i11,...,ij = (1 + η)t̃i11,...,ij ,

(B.15)

(i +1)

(i )

ũi11,...,ij = (1 + η)ũi11,...,ij ,

and, for k ≥ i1 + 2, (B.16)

(k+1)

(k)

(k−1)

(k)

t̃i1 ,...,ij = (1 + η)t̃i1 ,...,ij − η t̃i1 ,...,ij ,

(k−1)

(k−2)

ũi1 ,...,ij = (1 + η)ũi1 ,...,ij − η ũi1 ,...,ij .

Solving these recurrences gives, for k ≥ i1 + 1, (B.17)

(k)

t̃i1 ,...,ij =

η k−i1 − 1 (i1 +1) t̃ , η − 1 i1 ,...,ij

(k−1)

ũi1 ,...,ij =

η k−i1 − 1 (i1 ) ũi1 ,...,ij , η−1

which are nonnegative, with the η = 1 values read as limits as before.

□

To bound the propagators we need the magnitude of these coefficients. Since all coefficients of a given order share one sign, it suffices to track their aggregate sums X X (k) (k) (k) (k) (B.18) tj = t̃i1 ,...,ij , uj = ũi1 ,...,ij . k−1≥i1 >···>ij ≥0

k≥i1 >···>ij ≥1

We characterize these aggregate coefficients in two regimes, exactly when η = 1, and by upper bounds when η ∈ [0, 1). Proposition B.2. When η = 1, (k)

(k + j − 1)! , (k − j)!(2j)!

k ≥ 1, 0 ≤ j ≤ k ,

(k + j + 1)! , (k − j)!(2j + 1)!

k ≥ 0, 0 ≤ j ≤ k ,

(B.19)

(B.20)

(k)

(0)

together with t0 = 1.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Proof of Proposition B.2. Taking Hi := I in Lemma B.1 collapses the expansions onto their aggregate coefficients, 2

k X (k) (−1)j h2j tj I ,

(k)

T̃kη ({(1 − h2 )I}k−1 i=0 ) = t0 I +

j=1

(B.21) 2 (k) Ũkη ({(1 − h2 )I}ki=1 ) = u0 I +

k X

(k)

(−1)j h2j uj I .

j=1

At η = 1 with equal arguments, the multivariate polynomials reduce to the classical Chebyshev polynomials Tk and Uk (as noted after Definition 3.4), so 2

Tk (1 − h2 ) I =

(B.22)

k X (k) (−1)j h2j tj I ,

Uk (1 − h2 ) I =

j=0

k X

(k)

(−1)j h2j uj I .

j=0

The hypergeometric forms of Tk and Uk [2] give explicit expressions for Tk (1 − x) and Uk (1 − x) at x = h2 /2, and matching coefficients yields the stated formulas. □ For η ∈ [0, 1), explicit formulas are no longer available, but the following upper bounds suffice to control the propagator norms. Proposition B.3. For η ∈ [0, 1) and the same ranges of j, k as in Proposition B.2, Å ãj Ç å ãj+1 Ç å Å 1 k k 1 (k) (k) , uj ≤ . (B.23) tj ≤ 1−η j 1−η j Proof of Proposition B.3. We derive iterative inequalities for the coefficients from their (k) (k) recurrences and conclude by induction. We treat tj , and the argument for uj is analogous. For j ≥ 1, Lemma B.1 gives X X η k−i1 − 1 (i1 +1) (k) (k) t̃ . (B.24) tj = t̃i1 ,...,ij = η − 1 i1 ,...,ij k−1≥i1 >···>ij ≥0

k−1≥i1 >···>ij ≥0

(1)

2 By (B.14) and the initial conditions t̃0 = 1+η 4 (which is the coefficient of −h H0 ) and (0)

(i +1)

(i )

t0 = 1, we have t̃i11,...,ij ≤ 21 (1 + η) t̃i21,...,ij , so (k)

≤ 12 (1 + η)

X k−1≥i1 ≥j−1

(B.25) = 21 (1 + η)

X k−1≥i1 ≥j−1

η k−i1 − 1 η−1

é (i )

t̃i21,...,ij

i1 −1≥i2 >···>ij ≥0

η k−i1 − 1 (i1 ) 1 t ≤ η − 1 j−1 1 − η

(i )

1 tj−1 ,

k−1≥i1 ≥j−1

(i )

1 using the definition of tj−1 in the equality. We now induct on j. The case j = 0 holds

(k)

since t0 = 1 for all k. Assuming (B.23) for some j and all k ≥ j, then for k ≥ j + 1, Ç å Ç å X 1 j i 1 1 1 j+1 k (k) (B.26) tj+1 ≤ = , 1−η 1−η j 1−η j+1 k−1≥i1 ≥j

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

where the last step uses Pascal’s rule

k−1≥i1 ≥j

k j+1 .

i1 j =

This proves the bound

(k) for tj . (k)

(k)

For uj , the same steps give the iterative inequality uj

1 ≤ 1−η

(i1 −1) k≥i1 ≥j uj−1 .

k+1

(k)

Combined with the initial condition u0 = 12 (1 + η) η η−1−1 from Lemma B.1 and the P 1 −1 □ = kj , this yields the stated bound. identity k≥i1 ≥j ij−1 Part (b) of Proposition 3.5 concerns the HMC propagators across multiple outer loops. The corresponding matrix polynomials again have coefficients of alternating (k) (k) signs, related to the inner-loop coefficients tj and uj . We only need that the outerloop aggregate coefficients inherit any uniform bound on the inner-loop ones. Proposition B.4. When η = 1, the propagators in (3.9) expand as ←− Y 2 1 T̃m ({I − h2 Hl }j+m−1 ) T̃k1 l=j −1 j∈{k+im}N i=0

(B.27)

(N m+k)

= a0

I+

NX m+k

(−1)j h2j

j=1 ←− Y

(B.28)

(N m+k)

ai1 ,...,ij Hi1 · · · Hij ,

i1 >···>ij

2 1 1 T̃m ({I − h2 Hl }j+m−1 ) Ũk−1 l=j

−1 j∈{k+im}N i=0

(N m+k−1) = b0 I+

N m+k−1 X

(−1)j h2j

j=1

(N m+k−1)

bi1 ,...,ij

Hi1 · · · Hij ,

i1 >···>ij

P (N m+k) (N m+k) (N m+k−1) with all coefficients nonnegative. Writing aj = i1 >···>ij ai1 ,...,ij and bj = P (N m+k−1) for the aggregate coefficients, we have i1 >···>ij bi1 ,...,ij (B.29)

(N m+k)

≤

max

j≤i≤N m+k

(i)

tj ,

(N m+k−1)

≤

max

j≤i≤N m+k−1

(i)

uj .

Proof of Proposition B.4. Alternating signs. In each factor T̃k1 and Ũk1 with argu2 ments I − h2 Hi , a product of j matrices carries the sign (−1)j . This sign is multiplicative across the product of factors in (B.27) and (B.28), so every order-j product carries (−1)j , and the coefficients are nonnegative. Reduction to Chebyshev polynomials. Taking Hi := I collapses the expansions onto their aggregate coefficients, and at η = 1 the multivariate polynomials reduce to the classical Tk and Uk , giving 2 N 2 (N m+k) Tm (1 − h2 ) Tk (1 − h2 ) I = a0 I+

(B.30)

NX m+k

(N m+k)

(−1)j h2j aj

j=1 2 N 2 (N m+k−1) Tm (1 − h2 ) Uk−1 (1 − h2 ) I = b0 I+

N m+k−1 X j=1

(N m+k−1)

(−1)j h2j bj

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN 2

29 2

Writing θ = arccos(1 − h2 /2), so that Tk (1 − h2 ) = cos(kθ) and Uk−1 (1 − h2 ) = sin(kθ)/ sin θ, and expanding N Ç å N 1 X N (B.31) cos(mθ) = N exp(i(2i − N )mθ) , 2 i i=0

we pair conjugate exponentials into cosines and apply 2 cos x cos y = cos(x−y)+cos(x+ y) and 2 cos x sin y = sin(x + y) − sin(x − y). Rewriting the resulting trigonometric sums N back in terms of Chebyshev polynomials expresses Tm Tk as a binomially weighted sum of T · at shifted indices 2im ± k (for N even) or (2i − 1)m ± k (for N odd), all with N positive weights. For Tm Uk−1 the sin(x + y) − sin(x − y) identity instead produces a signed sum of U · at indices 2im ± k − 1 (for N even) or (2i − 1)m ± k − 1 (for N odd), (·) (·) the x − y terms entering with a minus sign. Matching coefficients with tj and uj then (N m+k)

and bj

(i)

writes each aj coefficients tj

(N m+k−1)

as a 2−N -weighted sum of inner-loop aggregate

and uj , with binomial weights and indicator factors 1j≤i . The sum

(N m+k)

(N m+k−1)

for aj has all positive weights, while that for bj inherits the minus signs above. Explicit formulas can be read off in this way, but the main results need only the bound below. P (i) (N m+k) Bound. Each aj is a sum of the form 21N i wi tj 1j≤i with nonnegative binoP (i) mial weights wi satisfying i wi = 2N . Since the tj are nonnegative by Lemma B.1, (i)

(i)

bounding each term by tj 1j≤i ≤ maxj≤i≤N m+k tj gives 1 X (N m+k) (i) (i) wi (B.32) aj ≤ N max tj = max tj . j≤i≤N m+k j≤i≤N m+k 2 i

(N m+k−1)

The same argument bounds bj

(i)

. Here some uj enter with negative signs, but

(i)

since the uj are nonnegative, discarding the subtracted terms only increases the bound, (i)

and the remaining terms are controlled as above by maxj≤i≤N m+k−1 uj .

□

Appendix C. Bounds for Leap-frog Propagators C.1. ℓ2 bounds. Our ℓ2 bound for T̃k1 adapts the discrete-time HMC contraction result for the position component, Lemma 19 in [11], replacing the state-dependent Hessians there by general matrices. The underlying argument carries through unchanged. Ex1 . tending this to the momentum component gives the bound for Ũk−1 The adaptation rests on decoupling the matrix sequence from the trajectory. Because the leap-frog scheme is nonlinear in the state-dependent Hessians, a trajectory contraction bound |xk −yk |ℓ2 ≤ c|x0 −y0 |ℓ2 does not by itself give a matrix-norm bound. We therefore analyze a linear system driven by an arbitrary, state-independent matrix sequence {Hi }. The contraction proof of [11] does not use that the Hessians are evaluated along the trajectory, so the adaptation is direct. The resulting bounds hold for any initial condition and hence translate to matrix-norm bounds. Instantiating with the integrated Hessians yields Proposition 4.1.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Proposition C.1 (ℓ2 bounds for leap-frog, general matrices). Let {Hi } be symmetric positive definite matrices with αI ⪯ Hi ⪯ βI for α, β as in Assumption 1.1. For kh ≤ √ 1 , evaluated at the canonical arguments 1/ 20β, the matrix polynomials T̃k1 and Ũk−1 2 {I − h2 Hi } with these general Hi (recall the shorthand introduced in Section B), satisfy (C.1)

|T̃k1 |ℓ2 ≤ 1 −

α (kh)2 , 10

1 |hŨk−1 |ℓ2 ≤

2 . 5β

Proof of Proposition C.1. We interpret the polynomials with general {Hi } as the propagator of a leap-frog iteration with varying quadratic potentials, which lets us bound their ℓ2 norms through the contraction and boundedness of the discrete-time dynamics, analyzed via their continuous-time interpolation. Connection with dynamics. Consider the time-inhomogeneous linear dynamics driven by {Hi }, 2

(C.2)

qk+1 = qk + hpk − h2 Hk qk , pk+1 = pk − h2 Hk qk − h2 Hk+1 qk+1 ,

which is the leap-frog scheme (1.3) with the quadratic potential Vk (q) = 12 q ⊤ Hk q at step k. For two initial states (x0 , p0 ), (y0 , p′0 ), the differences ∆qk := xk − yk , ∆pk := pk − p′k satisfy 2

(C.3)

∆qk+1 = (I − h2 Hk )∆qk + h∆pk , ∆pk+1 = ∆pk − h2 Hk ∆qk − h2 Hk+1 ∆qk+1 .

This matches (B.1) in the proof of Proposition 3.5, now with general {Hi } in place of the integrated Hessians. As that proof does not use the specific Hessian structure, the same argument gives the three-term recurrence ∆qk+1 = (2I − h2 Hk )∆qk − ∆qk−1 for k ≥ 1 and, by induction, (C.4)

1 ∆qk = T̃k1 ∆q0 + hŨk−1 ∆p0 .

1 Thus T̃k1 and hŨk−1 are the propagators of (C.2). We bound their ℓ2 norms by choosing 1 . initial conditions that isolate each, setting ∆p0 = 0 for T̃k1 , and ∆q0 = 0 for hŨk−1 ℓ2 norm of T̃k1 . Adapting the proof of [11, Lemma 19] to a fixed matrix sequence {Hi } rather than the leap-frog scheme of the target potential, we show that the position difference still contracts when the potential is decoupled from the trajectory. We include it both for completeness and to set up the Ũk1 bound. Introduce the continuous-time interpolation of (C.2),

(C.5)

dqt = p⌊t⌋h − (t − ⌊t⌋h )H⌊t⌋h q⌊t⌋h , dt dpt = − 21 (H⌊t⌋h q⌊t⌋h + H⌈t⌉h q⌈t⌉h ) , dt

with ⌊t⌋h := ⌊t/h⌋h and ⌈t⌉h := ⌈t/h⌉h. Take two trajectories from (x0 , p0 ) and (y0 , p′0 ) with p0 = p′0 , set zt := xt − yt , wt := pt − p′t , and at each grid point tk = kh define Φkh := Hk zkh and φkh := Φkh · zkh , with · the inner product. Then w0 = 0 and

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

|Φt |2ℓ2 ≤ βφt , and (C.5) gives dzt = w⌊t⌋h − (t − ⌊t⌋h )Φ⌊t⌋h , dt (C.6) dwt = − 21 (Φ⌊t⌋h + Φ⌈t⌉h ) . dt Let at = |zt |2ℓ2 and bt = 2zt · wt , so b0 = 0 and αakh ≤ φkh ≤ βakh since αI ⪯ Hk ⪯ βI. We track at , bt to bound at . From (C.6), dat = bt + δt , dt (C.7) dbt = −C1 αat + ηt , dt where C1 > 0 is chosen later and the perturbations are δt = 2zt · (w⌊t⌋h − wt − (t − ⌊t⌋h )Φ⌊t⌋h ) = (t − ⌊t⌋h )zt · (Φ⌈t⌉h − Φ⌊t⌋h ) , (C.8)

ηt = C1 αat + 2|w⌊t⌋h − (t − ⌊t⌋h )Φ⌊t⌋h |2ℓ2 − zt · (Φ⌊t⌋h + Φ⌈t⌉h ) − (t − ⌊t⌋h )(w⌊t⌋h − (t − ⌊t⌋h )Φ⌊t⌋h ) · (Φ⌈t⌉h − Φ⌊t⌋h ) .

The term δt is piecewise smooth, with dδt = zt · (Φ⌈t⌉h − Φ⌊t⌋h ) + (t − ⌊t⌋h )(w⌊t⌋h − (t − ⌊t⌋h )Φ⌊t⌋h ) · (Φ⌈t⌉h − Φ⌊t⌋h ) dt in the interior of each step, and at grid points δkh+ = 0 and δkh− = h(φkh − φ(k−1)h ) − h(hw(k−1)h − 12 h2 Φ(k−1)h ) · Φ(k−1)h . √

Solving the perturbed system (C.7) by variation of parameters, with st = sin(√CC1ααt) and 1 √ ct = cos( C1 αt) (so st is increasing and nonnegative on 0 ≤ t ≤ 2√πC α ), gives 1 Z t Z t at = ct a0 + ct−r δr dr + st−r ηr dr 0 ⌈t/h⌉−1

(C.9)

= ct a0 +

+ r=t− [st−r δr ]|r=kh r=kh− − [st−r δr ]|r=0+ +

k=1 ⌈t/h⌉−1

= ct a0 −

X k=1

Z t 0

r st−r (ηr + dδ dr )dr

Z t st−kh δkh− + 0

st−r (−2φ⌊r⌋h + C1 αa⌊r⌋h + εr )dr ,

where we integrated by parts and set εt := ε1t + ε2t + ε3t , ε1t = 2|w⌊t⌋h − (t − ⌊t⌋h )Φ⌊t⌋h |2ℓ2 , ε2t = −2(zt − z⌊t⌋h ) · Φ⌊t⌋h , ε3t = C1 α(zt − z⌊t⌋h ) · (zt + z⌊t⌋h ) . We aim to show at ≤ ct a0 for a suitable C1 , which requires bounding the remaining terms in (C.9). Using |Φ⌊t⌋h |2ℓ2 ≤ βφ⌊t⌋h , (C.10)

ε1t ≤ 6|w⌊t⌋h |2ℓ2 + 3(t − ⌊t⌋h )2 |Φ⌊t⌋h |2ℓ2 ≤ 6|w⌊t⌋h |2ℓ2 + 3βh2 φ⌊t⌋h .

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Substituting zt − z⌊t⌋h from (C.6) into ε2t , ε2t = −2(t − ⌊t⌋h )w⌊t⌋h · Φ⌊t⌋h + (t − ⌊t⌋h )2 |Φ⌊t⌋h |2ℓ2

(C.11)

≤ 2|w⌊t⌋h |2ℓ2 + 32 h2 |Φ⌊t⌋h |2ℓ2 ≤ 2|w⌊t⌋h |2ℓ2 + 23 βh2 φ⌊t⌋h ,

and similarly for ε3t , ε3t = C1 α|(t − ⌊t⌋h )w⌊t⌋h − 12 (t − ⌊t⌋h )2 Φ⌊t⌋h |2ℓ2

(C.12)

+ 2C1 αz⌊t⌋h · ((t − ⌊t⌋h )w⌊t⌋h − 12 (t − ⌊t⌋h )2 Φ⌊t⌋h ) .

The first term is bounded by 2C1 αh2 |w⌊t⌋h |2ℓ2 + 12 C1 αh4 βφ⌊t⌋h . For the second, |z⌊t⌋h |2ℓ2 = a⌊t⌋h ≤ φ⌊t⌋h /α gives 2C1 αz⌊t⌋h · (t − ⌊t⌋h )w⌊t⌋h ≤ h2 C12 α2 |z⌊t⌋h |2ℓ2 + |w⌊t⌋h |2ℓ2 ≤ h2 C12 αφ⌊t⌋h + |w⌊t⌋h |2ℓ2 . Since z⌊t⌋h · Φ⌊t⌋h ≥ 0, ε3t ≤ (1 + 2C1 αh2 )|w⌊t⌋h |2ℓ2 + C1 α(C1 h2 + 21 βh4 )φ⌊t⌋h ,

(C.13)

and combining the three estimates, εt ≤ (9 + 2C1 αh2 )|w⌊t⌋h |2ℓ2 +

(C.14)

2 2 4 1 2 βh + C1 α(C1 h + 2 βh )

φ⌊t⌋h .

For the boundary sum in (C.9), when t/h ∈ Z+ , summation by parts gives t/h−1

−

t/h−1

st−kh δkh− =

k=1

(C.15)

h(st−kh − st−(k−1)h )φ(k−1)h − hsh φt−h + hst−h φ0

k=2 t/h−1

hst−kh (hw(k−1)h − 21 h2 Φ(k−1)h ) · Φ(k−1)h .

k=1

Since st−kh ≤ st−r for r ∈ [(k − 1)h, kh), φkh ≥ 0, and h2 w(k−1)h · Φ(k−1)h ≤ 12 h|w(k−1)h |2ℓ2 + 21 h3 |Φ(k−1)h |2ℓ2 ≤ 12 h|w(k−1)h |2ℓ2 + 12 βh3 φ(k−1)h , we obtain, for t/h ∈ Z+ , t/h−1

(C.16)

−

Z h st−kh δkh− ≤ 0

k=1

st−r φ⌊r⌋h dr + 12

Z t 0

st−r (|w⌊r⌋h |2ℓ2 + βh2 φ⌊r⌋h )dr .

From (C.6) with w0 = 0, Cauchy–Schwarz and |Φ⌊τ ⌋h |2ℓ2 ≤ βφ⌊τ ⌋h give (C.17) ÇZ ⌊t⌋ å 2 Z ⌊t⌋h Z ⌊t⌋h h 2 1 1 |w⌊t⌋h |ℓ2 = 2 (Φ⌊τ ⌋h + Φ⌈τ ⌉h )dτ ≤ 2 β⌊t⌋h φ⌊τ ⌋h dτ + φ⌈τ ⌉h dτ . 0

To bound

ℓ2

2 + 0 st−r |w⌊r⌋h |ℓ2 dr, we note that for t/h ∈ Z ,

ÇZ ⌊r⌋ h

Z t st−r ⌊r⌋h

(C.18) 0

å φ⌈τ ⌉h dτ

Z t ÇZ t dr = h

⌊τ ⌋h

å st−r ⌊r⌋h dr φ⌊τ ⌋h dτ .

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

As st is increasing, we decompose the integral in r into sub-intervals of length h, and direct calculation of the sum of ⌊r⌋h terms gives å ÇZ ⌊τ ⌋ +h Z t h t2 st−r dr st−r ⌊r⌋h dr ≤ (C.19) . 2h ⌊τ ⌋h ⌊τ ⌋h Plugging into (C.18) and switching the order of integration yield an upper bound å Z t ÇZ ⌊τ ⌋h +h Z t2 t t2 φ dτ = st−r dr st−r φ⌊r⌋h dr . (C.20) 2h ⌊τ ⌋h 2 h ⌊τ ⌋h h We can similarly derive for the other term in (C.17) that å ÇZ ⌊r⌋ Z t Z h t2 t−h φ⌊τ ⌋h dτ dr ≤ st−r ⌊r⌋h (C.21) st−r φ⌊r⌋h dr , 2 0 0 0 and combining these two bounds yields Z t Z t 2 1 2 st−r φ⌊r⌋h dr . (C.22) st−r |w⌊r⌋h |ℓ2 dr ≤ 2 βt 0

Hence the last two terms of (C.9) are bounded by ⌈t/h⌉−1

−

Z t 0

k=1

(C.23)

Z h ≤ 0

st−r (−2φ⌊r⌋h + C1 αa⌊r⌋h + εr )dr

st−kh δkh− + Z t

st−r φ⌊r⌋h dr +

Z t + 0

st−r (−2φ⌊r⌋h + C1 αa⌊r⌋h )dr

2 2 st−r ((5βh2 + C12 αh2 + 21 C1 αβh4 )φ⌊r⌋h + ( 19 2 + 2C1 αh )|w⌊r⌋h |ℓ2 )dr .

Applying (C.22) and αakh ≤ φkh , ⌈t/h⌉−1

−

X k=1

(C.24)

Z t ≤ 0

Z t st−kh δkh− + 0

st−r (−2φ⌊r⌋h + C1 αa⌊r⌋h + εr )dr

2 2 st−r −1 + C1 + 5βh2 + C12 αh2 + 12 C1 αβh4 + ( 19 4 + C1 αh )βt φ⌊r⌋h dr 2 2 2 3 st−r −1 + C1 + ( 39 + C + C βt )βt φ⌊r⌋h dr , 1 1 4 2

using t/h ∈ Z+ and α ≤ β in the last line. Taking C1 = 4/9 and βt2 ≤ 1/20 makes the bracketed coefficient negative, so (C.9) gives at ≤ ct a0 . The same choice gives √ t ≤ π/(2 C1 α), ensuring st is monotone on the interval, and 1 ct ≤ 1 − 21 C1 αt2 + 16 C12 α2 t4 ≤ 1 − 12 C1 αt2 + 120 C12 αt2 (α/β) 1 C12 )αt2 ≤ 1 − 51 αt2 . ≤ 1 − ( 12 C1 − 120

Thus for t = kh with p0 = p′0 , |xt − yt |2ℓ2 ≤ (1 − α5 (kh)2 )|x0 − y0 |2ℓ2 , and by (C.4), α (C.25) |T̃k1 |ℓ2 ≤ 1 − (kh)2 . 10

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

1 . By (C.4), bounding |hŨ 1 | ℓ2 norm of hŨk−1 k−1 ℓ2 again reduces to bounding at , now with initial conditions z0 = 0, w0 = p0 − p′0 , so a0 = b0 = 0 and Z t Z t at = ct−r δr dr + st−r ηr dr 0

(C.26)

0 ⌈t/h⌉−1

=−

Z t st−kh δkh− + 0

k=1

st−r (−2φ⌊r⌋h + C1 αa⌊r⌋h + εr )dr .

The boundary decomposition (C.15) still holds, and since now φ0 = 0, it refines (C.16) to Z t t/h−1 X 1 (C.27) − st−kh δkh− ≤ 2 st−r (|w⌊r⌋h |2ℓ2 + βh2 φ⌊r⌋h )dr . 0

k=1

The bound (C.14) on εr still applies, so combining with αa⌊r⌋h ≤ φ⌊r⌋h , Z t at ≤ st−r (−2 + C1 + 5βh2 + C12 αh2 + 12 C1 αβh4 )φ⌊r⌋h dr 0 (C.28) Z t 2 2 st−r ( 19 + 2 + 2C1 αh )|w⌊r⌋h |ℓ2 dr . 0

For w⌊r⌋h , repeating the derivation of (C.17) but with general w0 , for any C2 > 0, 2 Z ⌊t⌋h 2 1 (Φ⌊s⌋h + Φ⌈s⌉h )ds |w⌊t⌋h |ℓ2 = w0 − 2 0 ℓ2 (C.29) ÇZ å ≤ (1 + C2 )|w0 |2ℓ2 + (1 + C12 ) 21 β⌊t⌋h

⌊t⌋h

φ⌊τ ⌋h + φ⌈τ ⌉h dτ

With the derivation of (C.22), this bounds (C.28) by (C.30) Z t at ≤ st−r −2 + C1 + 5βh2 + C12 βh2 + 12 C1 β 2 h4 0

2 2 +(1 + C12 )( 19 4 + C1 βh )βt

2 φ⌊r⌋h dr + ( 19 2 + 2C1 αh )(1 + C2 )

Z t 0

st−r |w0 |2ℓ2 dr .

With C1 = 4/9, βt2 ≤ 1/20, and C2 = 1/4, the coefficient of φ⌊r⌋h is negative. Since 2 + αh2 ≤ βt2 gives 19 2 + 2C1 αh ≤ 10, for t/h ∈ Z , √ Z t 16(1 − cos( C1 αt)) 2 (C.31) at ≤ 16|w0 |2ℓ2 st−r dr = |w0 |2ℓ2 ≤ 8t2 |w0 |2ℓ2 ≤ |w0 |2ℓ2 , C α 5β 1 0 using cos t ≥ 1 − 21 t2 . By (C.4), (C.32)

1 |hŨk−1 |ℓ2 ≤

2 . 5β □

Proof of Proposition 4.1. This follows from Proposition C.1, since the integrated Hessian satisfies αI ⪯ Hi ⪯ βI by Assumption 1.1. □

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

C.2. ℓ∞ bounds under weak interactions. Proof of Proposition 4.3. Following the decomposition of ∇2 V in Assumption 1.3, set Z Z (D) (O) 2 (D) Hi := ∇ V dνi , Hi := ∇2 V (O) dνi , (D)

(O)

so that Hi = Hi + Hi . We begin with |T̃k1 |ℓ∞ , which we split into a diagonal contribution and a perturbation. By the sign structure of the coefficients in Lemma B.1, 2

(D) k−1 }i=0 )|ℓ∞

|T̃k1 |ℓ∞ ≤ |T̃k1 ({I − h2 Hi {z |

}

(a)

(C.33)

k X

(D)

(k)

h2j

(D)

t̃i1 ,...,ij |Hi1 · · · Hij − Hi1 · · · Hij |ℓ∞ .

k−1≥i1 >···>ij ≥0

j=1

}

(b)

The matrix in (a) is diagonal, so its ℓ∞ norm equals its ℓ2 norm. The ℓ2 bound of Proposition 4.1 gives α 2 (D) (C.34) (a) = |T̃k1 ({I − h2 Hi }k−1 (kh)2 . i=0 )|ℓ2 ≤ 1 − 10 (D)

For (b), expanding Hi1 · · · Hij = (Hi1 term yield

(C.35)

(O)

(D)

+ Hi1 ) · · · (Hij

(D) (D) |Hi1 · · · Hij − Hi1 · · · Hij |ℓ∞ ≤

(O)

+ Hij ) and bounding each

j Ç å X j (D) (O) (sup |Hk |ℓ∞ )j−l (sup |Hk |ℓ∞ )l l k k l=1

α j ) − βj , ≤ (β + 50

so that (b) ≤

k X j=0

(C.36) =

k X

h2j

(k) α j t̃i1 ,...,ij (β + 50 ) − βj

k−1≥i1 >...>ij ≥0 (k)

h2j tj

2 2 α α j (β + 50 ) − β j = Tk (1 + h2 (β + 50 )) − Tk (1 + h2 β) ,

j=0 (k)

using that the tj the estimates,

are the coefficients of Tk (1+ 12 x) =

(k) j j=0 tj x , by (B.22). Combining

α 2 2 α )) − Tk (1 + h2 β) . (kh)2 + Tk (1 + h2 (β + 50 10 It remains to show the perturbation in Tk is controlled, √ so the contraction factor stays of the same order. For x ≥ 1, Tk (x) = cosh(k log(x + x2 − 1)), with p d k Tk (x) = sinh(k log(x + x2 − 1)) √ . dx x2 − 1 (C.37)

|T̃k1 |ℓ∞ ≤ 1 −

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN 2

Substituting x = 1 + h2 z with z > 0, d Tk (x) dx

(C.38)

≤ 21 2 x=1+ h z 2

Å ãk q √ k h2 h2 z √ , 1+ 2 z+h z 1+ 4 h z

α so by the mean value theorem, for z ∈ [β, β + 50 ],

Ä Ä ä ä h2 α 2 2 α − Tk 1 + h2 β ≤ Tk 1 + h2 β + 50 100 (C.39)

sup 2

x=1+ h2 z, α z∈[β,β+ 50 ]

d Tk (x) dx

» » khα khα 51 51 2 √ exp kh 100 √ , hβ + 51 β 1 + h β ≤ 50 200 200 β 100 β √ where the last step uses h ≤ kh ≤ 1/ 20β. Substituting into (C.37), ≤

|T̃k1 |ℓ∞ ≤ 1 −

(C.40)

α khα √ . (kh)2 + 10 100 β

1 This remains a contraction as in the ℓ2 norm for large enough kh. Plugging in kh = √20β α gives the upper bound 1 − 400β . 1 . Lemma B.1 gives An analogous argument bounds hŨk−1 2

(D) k−1 }i=1 )|ℓ∞

1 1 |ℓ∞ ≤ |hŨk−1 ({I − h2 Hi |hŨk−1 {z |

}

(a)

(C.41)

k−1 X

h2j+1

(k)

(D)

k−1≥i1 >...>ij ≥1

j=1

| 2

(C.42)

k−1 X j=0

}

(b) (D) k−1 }i=1 )|ℓ2 ≤

1 ({I − h H where, as before, (a) = |hŨk−1 i 2

(b) ≤

(D)

ũi1 ,...,ij |Hi1 · · · Hij − Hi1 · · · Hij |ℓ∞ ,

2 5β by Proposition 4.1, and

(k) α j ũi1 ,...,ij (β + 50 ) − βj

h2j+1

k−1≥i1 >...>ij ≥1 2

α α = hUk−1 (1 + h2 (β + 50 )) − hUk−1 (1 + h2 β) ≤ hUk−1 (1 + h2 (β + 50 )) , √ √ since Uk−1 (x) ≥ 0 for x ≥ 1. Using Uk−1 (x) = sinh(k log(x + x2 − 1))/ x2 − 1 for x ≥ 1, » » 1 13 51 51 2 (C.43) (b) ≤ √ exp kh 100 hβ + 51 β 1 + h β ≤ √ . 50 200 2 β 20 β

Combining these in (C.41), (C.44)

3 1 |hŨk−1 |ℓ∞ ≤ √ . 2 β □

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

C.3. ℓ∞ bounds under sparse interactions. Here we use the sparsity of the Hessians to bound the propagators in ℓ∞ . The next lemma, adapted from Proposition C.1 of [17], shows that the integrated Hessian products inherit a sparsity property from Assumption 1.5. We include the proof for completeness. Lemma C.2. Suppose V satisfies Assumption 1.5, and let ν1 , . . . , νr be probability measures on Rd . Then the product ã ã ÅZ ã ÅZ ÅZ 2 2 2 ∇ V dνr ∇ V dν2 · · · (C.45) ∇ V dν1 is sr -sparse. More precisely, its (i, j) entry is nonzero only if j ∈ N2r (i). Proof of Lemma C.2. The (i, j) entry of the product is nonzero only if j ∈ Nr (i), or if i ∈ Nr (k) and j ∈ Nr (k) for some k. Since max1≤i≤d |N2r (i)| = sr , each row has at most sr nonzero entries, so the matrix is sr -sparse. □ We now prove an extended version of Proposition 4.5. Proposition C.3 (Extended version of Proposition 4.5). Suppose Assumptions 1.1 and 1.5 hold, √and let the integrated Hessians Hi be as in Proposition 4.1. Take ri = √ √ log d ⌈ih βe + log(5/3) ⌉, and let mh ≤ 1/ 20β with integers N ≥ 0 and k < m. Define (C.46)

AN m+k :=

←− Y

2 1 ) T̃k1 , T̃m ({I − h2 Hj+l }m−1 l=0

−1 j∈{k+im}N i=0

(C.47)

BN m+k−1 :=

←− Y

2 1 1 T̃m ({I − h2 Hj+l }m−1 l=0 ) hŨk−1 ,

−1 j∈{k+im}N i=0

with the product read as the identity when N = 0. Then 4 √ srN m+k . |BN m+k−1 |ℓ∞ ≤ √ β R R Moreover, for arbitrary probability measures ν0 , ν0′ and H := ∇2 V dν0 , H ′ := ∇2 V dν0′ , (C.48)

√ |AN m+k |ℓ∞ ≤ 4 srN m+k ,

(C.49)

√ |AN m+k H|ℓ∞ ≤ 4β srN m+k , p √ |BN m+k−1 H|ℓ∞ ≤ 4 β srN m+k ,

3√ |BN m+k−1 HH ′ |ℓ∞ ≤ 5β 2 srN m+k .

Proof of Proposition C.3. The bounds follow a common framework, so we give the details for AN m+k . The strategy is to split the matrix expansion by order in h at a threshold r that balances sparsity against magnitude, as for a tridiagonal matrix, loworder powers stay sparse (so their ℓ∞ and ℓ2 norms are comparable), while the growing density of higher-order powers is offset by the decay in h2j . The sparsity here is controlled by sr .

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

For any threshold r, Lemma C.2 gives that the entries of each order-r product are supported on N2r (i), and the absolute row sum decomposes as d X

|(AN m+k )ij | =

j=1

|(AN m+k )ij |

j ∈N / 2r (i)

j∈N2r (i)

(C.50)

|(AN m+k )ij | +

s X » ≤ |N2r (i)| |(AN m+k )ij |2 +

j∈N2r (i)

j ∈N / 2r (i)

√ ≤ sr |AN m+k |ℓ2 +

|(AN m+k )ij |

|(AN m+k )ij | .

j ∈N / 2r (i)

By Proposition B.4, (N m+k) AN m+k = a0 I+

r X

(−1)j h2j

j=1

(C.51) +

NX m+k

(−1)j h2j

j=r+1

(N m+k)

ai1 ,...,ij Hi1 · · · Hij

i1 >···>ij

(N m+k)

ai1 ,...,ij Hi1 · · · Hij .

i1 >···>ij

For j ∈ / N2r (i), the entry (AN m+k )ij can only come from the second line, so (C.52)

|(AN m+k )ij | ≤

NX m+k j=r+1

j ∈N / 2r (i)

(−1)j h2j

(N m+k)

ai1 ,...,ij Hi1 . . . Hij

i1 ,...,ij

ℓ∞

We control this ℓ∞ norm through the ℓ2 norm, (C.53)

NX m+k

(−1) h

j=r+1

(N m+k) ai1 ,...,ij Hi1 . . . Hij

j 2j

ℓ2

i1 ,...,ij

≤

(N m+k)

NX m+k

(N m+k)

h2j β j aj

j=r+1 (i)

By Proposition B.4, aj is bounded by maxj≤i≤N m+k tj , the Chebyshev coefficients whose explicit formulas are given in Proposition B.2. It therefore suffices to (i) bound h2j β j tj . For i > 0, (i + j − 1)! i(i + j − 1) · · · (i + 1) i(i − 1) · · · (i − j + 1) = (i − j)!(2j)! (2j) · · · (j + 1) j! (C.54) Ä Ä 2 ääj (i + j)j ij (i + j)j (ie)j ≤ ≤ = e ji 2 + ji , j 2j j j! j √ using j! ≥ j j / exp(j) (from exp(j) ≥ j j /j!). Hence for j ≥ (N m + k) βeh and i ≤ N m + k, Ä Ä 2 ääj Ä p äj j (i) (C.55) h2j β j tj ≤ h2 βe ji 2 + ji ≤ e−1 + βh ≤ 35 , √ where we use h ≤ mh ≤ 1/ 20β. Since this holds uniformly over i ≤ N m + k, j (N m+k) Proposition B.4 gives h2j β j aj ≤ 35 , so (C.53) yields (i)

tj = i

(C.56)

NX m+k j=r+1

j 2j

(−1) h

X i1 ,...,ij

(N m+k) ai1 ,...,ij Hi1 . . . Hij

≤ 2

ℓ

NX m+k j=r+1

3 j ≤3 5

3 r+1 . 5

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

√ Applying the norm inequality |M |ℓ∞ ≤ d |M |ℓ2 , valid for any matrix M ∈ Rd×d , to (C.52), X √ r+1 |(AN m+k )ij | ≤ 3 d 53 (C.57) . j ∈N / 2r (i)

NX m+k

(−1) h

j=r

(N m+k) ai1 ,...,ij Hi1 · · · Hij H

j 2j

i1 >···>ij

so taking rN m+k as before, √ |AN m+k H|ℓ∞ ≤ β srN m+k + 3β

(C.60)

ℓ2

≤β

NX m+k

3 j ≤ 3β 5

3 r , 5

j=r

√ √ 3 (N m+k)h βe ≤ 4β srN m+k . 5

For BN m+k−1 , the same decomposition gives, for any r, (C.61)

d X

|(BN m+k−1 )ij | ≤

√

sr |BN m+k−1 |ℓ2 +

X √ N m+k−1 (N m+k−1) d h2j+1 β j bj . j=r+1

j=1 (N m+k−1)

The difference from AN m+k is that bj

(i)

depends on uj . By Proposition B.2,

(i + j + 1)j (i + 1)j+1 ej+1 (i + j + 1)! ≤ , (i − j)!(2j + 1)! jj j j+1 √ so by Proposition B.4, for j ≥ (N m + k) βeh and i ≤ N m + k − 1, ääj h√β(i+1)e j+1 j 1 Ä p Ä 1 (N m+k−1) (C.63) h2j+1 β j bj ≤ √ h β i+1 + 1 ≤ √ 35 , j j β β and hence √ X √ N m+k−1 3 d 3 r+1 2j+1 j (N m+k−1) (C.64) d h β bj ≤ √ . 5 β j=r+1 » 2 (Proposition 4.1), (C.61) gives Taking r = rN m+k and using |BN m+k−1 |ℓ2 ≤ 5β q (N m+k)h√βe 4 √ 3 2 √ (C.65) |BN m+k−1 |ℓ∞ ≤ 5β srN m+k + √ 35 ≤√ srN m+k . β β The bounds for BN m+k−1 H and BN m+k−1 HH ′ follow from similar decompositions: » p 3 (N m+k)h√βe p √ √ (C.66) |BN m+k−1 H|ℓ∞ ≤ 2β s + 3 β 5 ≤ 4 β srN m+k , rN m+k 5 q 3√ 3√ 25 3 (N m+k)√βeh ′ ∞ (C.67) |BN m+k−1 HH |ℓ ≤ 2β5 srN m+k + β 2 53 ≤ 5β 2 srN m+k . 6 (C.62)

(i)

uj =

□

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Appendix D. Bounds for BAOAB Propagators In this section we bound the BAOAB propagators. The key difference from the leapfrog case is that we work with weighted ℓ2 and ℓ∞ norms, defined through the weight matrix in (2.9), which is what makes the propagators contractive. D.1. ℓ2 bounds. The ℓ2 bounds rest on the BAOAB contraction result, Theorem 5.1 of [34]. √ , then for Lemma D.1 (Theorem 5.1 in [34]). Let Assumption 1.1 hold. If h ≤ 21−η β

any initial states (x0 , p0 ), (y0 , p′0 ) ∈ R2d and any standard normal sequence {ξi }k−1 i=0 , the BAOAB chains driven by the same noise, (D.1) ξk−1 ξk−1 0 0 (y0 , p′0 ) , · · · Uξbaoab,h (x0 , p0 ) , (yk , p′k ) = Ubaoab,h · · · Uξbaoab,h (xk , pk ) = Ubaoab,h satisfy almost surely (D.2)

k−1

|(xk − yk , pk − p′k )|ℓ2w ≤ 7(1 − c(h)) 2 |(x0 − y0 , p0 − p′0 )|ℓ2w , 2

αh with c(h) = 4(1−η) .

As for HMC, we prove a general-matrix version, Proposition D.2 below, by reading the polynomials with general {Hi } as BAOAB propagators with varying quadratic potentials and turning the contraction of Lemma D.1 into operator-norm bounds. Proposition D.2 (ℓ2 bounds for BAOAB, general matrices). Let {Hi } be symmetric √ positive definite with αI ⪯ Hi ⪯ βI for α, β as in Assumption 1.1. For h ≤ 21−η β 2

η αh and c(h) = 4(1−η) , the polynomials T̃kη and Ũk−1 defined with {Hi }, and the full-space w weighted propagator Mbaoab,k , satisfy

(D.3)

k−1

w |ℓ2 ≤ 7(1 − c(h)) 2 , |Mbaoab,k

√ k−1 (D.4) |T̃kη |ℓ2 ≤ 7 2(1 − c(h)) 2 ,

√ k−1 1 η 2 √ . −bT̃kη + hŨk−1 ≤ 7 2 (1 − c(h)) ℓ2 a − b2 Proof of Proposition D.2. Consider the time-inhomogeneous linear dynamics driven by {Hi }, p 2 qk+1 = qk + h2 (1 + η)pk − h4 (1 + η)Hk qk + h2 1 − η 2 ξk , (D.5) p pk+1 = η(pk − h2 Hk qk ) + 1 − η 2 ξk − h2 Hk+1 qk+1 ,

which is the BAOAB scheme (1.4) with the varying potential Vk (q) = 12 q ⊤ Hk q at step k. For two trajectories coupled through the same noise {ξi }, the differences ∆qk , ∆pk obey 2

(D.6)

∆qk+1 = (I − h4 (1 + η)Hk )∆qk + h2 (1 + η)∆pk , ∆pk+1 = η(∆pk − h2 Hk ∆qk ) − h2 Hk+1 ∆qk+1 . 2

These yield the same three-term recurrence ∆qk+1 = (1 + η)(I − h2 Hk )∆qk − η∆qk−1 as in the proof of Proposition 3.5, so the same induction gives (D.7)

η ∆qk = T̃kη ∆q0 + hŨk−1 ∆p0 ,

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

η identifying T̃kη and hŨk−1 as the position-difference propagators of (D.5). Although Lemma D.1 is stated for a fixed potential, its proof uses only the spectral bounds αI ⪯ Hi ⪯ βI. [34] employs the decomposition

(D.8)

(BAOAB)k = BAO(ABAO)k−1 AB

and shows that the core operator ABAO is a weighted ℓ2 contraction (giving the (1−c(h)) factors) while the boundary terms BAO and AB stay bounded, all under those eigenvalue bounds. Since each Hi satisfies them, the result carries over to the varying-potential case, giving under (D.5) (D.9)

k−1

|(∆qk , ∆pk )|ℓ2w ≤ 7(1 − c(h)) 2 |(∆q0 , ∆p0 )|ℓ2w .

w We now derive the operator-norm bounds. Since Mbaoab,k = W Mbaoab,k W −1 maps W (∆q0 , ∆p0 ) to W (∆qk , ∆pk ), (D.9) is precisely (D.3). For the position blocks, the norm equivalence (2.10) gives » √ (D.10) |∆qk |ℓ2 ≤ |∆qk |2ℓ2 + a|∆pk |2ℓ2 ≤ 2|(∆qk , ∆pk )|ℓ2w ,

so √ k−1 |∆qk |ℓ2 ≤ 7 2(1 − c(h)) 2 |W (∆q0 , ∆p0 )|ℓ2 . ó î η W −1 W (∆q0 , ∆p0 ), the bound (D.11) gives Writing (D.7) as ∆qk = T̃kη hŨk−1 ó î η √ k−1 η W −1 2 ≤ 7 2(1 − c(h)) 2 , (D.12) T̃k hŨk−1

(D.11)

ℓ

and since (D.13)

ó î η î η ó η η b 1 T̃ η + √a−b hŨk−1 W −1 = T̃k − √a−b , T̃k hŨk−1 2 k 2

the block-wise bounds (D.4) follow.

□

We extract a corollary in ℓ∞ norms from this proof. Corollary D.3. Under the assumptions and notation of Proposition D.2, any state difference (∆qk , ∆pk ) under the varying-potential scheme (D.6) satisfies √ k−1 , (D.14) |(∆qk , ∆pk )|ℓ∞ ≤ 7 2d(1 − c(h)) 2 |(∆q0 , ∆p0 )|ℓ∞ w w (D.15)

√ k−1 . |∆qk |ℓ∞ ≤ 14 d(1 − c(h)) 2 |(∆q0 , ∆p0 )|ℓ∞ w

Proof of Corollary follows from (D.9) and the matrix norm equiv√ D.3. Inequality (D.14) 2d×2d alence |M |ℓ∞ ≤ 2d|M |ℓ2 for M ∈ R . For (D.15), (D.11) gives (D.16) √ √ k−1 k−1 |∆qk |ℓ∞ ≤ |∆qk |ℓ2 ≤ 7 2(1−c(h)) 2 |(∆q0 , ∆p0 )|ℓ2w ≤ 14 d(1−c(h)) 2 |(∆q0 , ∆p0 )|ℓ∞ . w □ Proof of Proposition 4.2. This follows from Proposition D.2, since the integrated Hessians satisfy the assumption. □

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

D.2. ℓ∞ bounds under weak interactions. We bound the BAOAB propagators in ℓ∞ following [34], decomposing the scheme into the operators ABAO, BAO, and AB as in (D.8), where ABAO drives the contraction and the boundary terms BAO, AB stay bounded. The contraction remains under the weak interactions. Unlike HMC, a diagonal Hessian does not give a diagonal propagator here. On the joint (q, p) space the√propagator is a 2 × 2 block matrix with diagonal blocks. The crude inequality |M |ℓ∞ ≤ 2|M |ℓ2 (from the two nonzero entries per row) would not preserve the ℓ2 contraction, so we bound each block in ℓ∞ directly, using the ℓ2 bounds as a guide to the contraction rate. w Proof of Proposition 4.4. Estimate for Mabao . The matrix is   1 + η 2 4 η (1−η) 1 2 2 3 h √a−b2 H0 I − 2 h H0 − 1−η h H0 w . (D.17) Mabao (H0 ) =  √ 1 −ηh a − b2 H0 η(I + (− 12 + 1−η )h2 H0 )

For the second block of rows, the ℓ1 norm of each row is at most p 2 21 1 1 (D.18) ηh a − b2 · 21 20 β + η(1 + (− 2 + 1−η )h · 20 β) . 2

h α This increases in η ∈ [0, 1], while the target bound 1 − 8(1−η) decreases in η. Comparing √ them at the largest admissible value η = 1 − 2 βh confirms the target bound for the second block. (D) (O) For the first block of rows, split off the diagonal and off-diagonal parts H0 , H0 (defined as in the proof of Proposition 4.3): ï ò η 1 (D)

(D.19)

(D)

η (1−η) I − 12 h2 H0 − 1−η h2 H0 h3 √a−b H0 2 ï ò 1 + η + − 1 h2 H0(O) − η h2 H0(O) h3 4 √ (1−η)2 2 H0(O) . 2

1−η

a−b

(D)

In the first matrix, H0 is diagonal, so each row reduces √ to a scalar λ ∈ [α, β]. Since η 1 2 2 β increases in η and is below 1 at η = 1 − 2 βh, it is below 1 throughout, h β + h 2 1−η and the ℓ1 norm of the row is (D.20)

η h2 λ + h3 1 − 21 h2 λ − 1−η

√ gives The condition h ≤ 21−η β

√

λ.

√

1 a − b2 ≥ 2√3β , so the last term is at most 4√ h2 λ + 3

√1 η h2 λ, and (D.20) is bounded by 1 − c(h). 3 1−η the ℓ∞ norm of the first block of rows is

(D.21)

1 + η 2 4 √ (1−η) a−b2

Adding the off-diagonal contribution,

ã Å 1 + η (O) η 2 3 4 √ (1−η)2 1 |H0 |ℓ∞ ≤ 1 − 12 c(h) . 1 − c(h) + ( 2 + 1−η )h + h a−b2

w . By the scheme definition in (2.8), Estimate for Mbao " # h3 (1+η) h2 η h2 h √ 1 I − H − H (− I + H ) 0 0 0 2 4 2 4(1−η) w a−b2 (D.22) Mbao (H0 ) = . √ 2(1−η) h2 η − h2 η a − b2 H0 ηI + 2(1−η) H0

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

We bound each 2×2 block in ℓ∞ . The (1, 1) block is bounded by 1+( h4 + h2 · 2√1 β ) 21 20 β ≤

√ 2√ β h h 1 21 101 383 √ , the (2, 1) block by 1 √1 h· 21 β ≤ 21 , , the (1, 2) block by ( + 4 ·2· 4β · 20 β) ≤ 160 320 2 β 20 80 3 2 3 181 and the (2, 2) block by 12 · 2√1 β h · 21 β + 1 ≤ . Summing the block bounds over each 20 160

row gives the result. w . By (2.8), Estimate for Mab " h2 H0 I − 2(1−η) w √ (D.23) Mab (H0 ) = h − 2 a − b2 H0

3 (1+η) √ 1 ( h I + h4(1−η) 2 H0 ) a−b2 2 h2 h2 I + ( 2(1−η) − 4 )H0

# .

h 21 181 By the same block-wise bounding, the (1, 1) block is at most 1 + 4√ β ≤ 160 , the β 20 √

1 21 h 21 101 21 √ √ (1, 2) block at most 2√3β ( h2 + h2 4β 20 β) ≤ 160 3 , the (2, 1) block at most 2 β 20 β ≤ 80 , 181 and the (2, 2) block at most 1 + h2 2√1 β 21 20 β ≤ 160 . A row sum concludes the proof.

□

D.3. ℓ∞ bounds under sparse interactions. This proof follows that of Proposition C.3, with two changes. The ℓ2 input is now the BAOAB contraction of Proposition 4.2, and the coefficient bounds come from the η < 1 estimates of Proposition B.3 rather than the exact η = 1 formulas. We therefore give only the parts that differ. η . As in (C.50), for each Proof of Proposition 4.6. Write Ak := T̃kη and Bk−1 := hŨk−1 threshold r, d X

(D.24)

|(Ak )ij | ≤

√

sr |Ak |ℓ2 +

j=1

|(Ak )ij | ,

j ∈N / 2r (i)

√ k−1 where |Ak |ℓ2 ≤ 7 2(1 − c(h)) 2 by Proposition 4.2, and the far-neighbor sum comes from the terms of order greater than r, (D.25)

|(Ak )ij | ≤

k X

whose ℓ2 norm is at most

(D.26)

(k) tj ≤

j=r+1 h

t̃i1 ,...,ij Hi1 . . . Hij

k−1≥i1 >···>ij ≥0

j=r+1

j ∈N / 2r (i)

(k)

(−1)j h2j

ℓ∞

2j β j t(k) . By Proposition B.3, j

äj 1 1−η

Ç å Ä äj k ke ≤ (1−η)j , j

β 2 using j! ≥ j j / exp(j). Hence for r > ke 1−η h , (D.27) k X X (k) (−1)j h2j t̃i1 ,...,ij Hi1 . . . Hij j=r+1

k−1≥i1 >···>ij ≥0

≤ 2

ℓ

k Ä äj X ke h2 (1−η)j β j=r+1

≤

k X j=r+1

exp(−j) ≤ exp(−r) .

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

√ √ 2β 2 Substituting into (D.24) with |M |ℓ∞ ≤ d|M |ℓ2 and taking rk := ⌈ ke 1−η h + log d⌉, Ä ä √ √ k−1 2β 2 |Ak |ℓ∞ ≤ 7 2 srk (1 − c(h)) 2 + exp − ke 1−η h (D.28) Ä ä √ √ αh2 (k − 1) . ≤ (7 2 + 1) srk exp − 8(1−η) For Bk−1 , Proposition B.3 gives Ç å äj+1 Ä (k−1)e äj Ä k−1 1 (D.29) , ≤ 1−η j j P 2j+1 β j u(k−1) ≤ so the terms of order greater than r have ℓ2 norm at most k−1 j=r+1 h j h exp(−r) = b exp(−r). Combining the tail estimates for A and B , k k−1 1−η (k−1) uj ≤

(D.30)

äj+1 1 1−η

ä Ä b √ 1 − √a−b A + B k k−1 2 a−b2

j ∈N / 2r (i)

√ b ≤ 2 d √a−b exp(−r) , 2

√ k−1 b 1 while Proposition 4.2 gives | − √a−b Ak + √a−b Bk−1 |ℓ2 ≤ 7 2(1 − c(h)) 2 . Hence, 2 2

b √ , simplifying √a−b via h ≤ 21−η 2 β

(D.31)

Ä ä √ k−1 2β √ b 1 2b − √a−b Ak + √a−b Bk−1 ∞ ≤ srk 7 2(1 − c(h)) 2 + √a−b exp − ke h2 2 2 2 1−η ℓ Ä √ ä√ Ä ä 2 αh2 srk exp − 8(1−η) ≤ 7 2 + √3 (k − 1) . □ Appendix E. Discretization Error Analysis

We now bound the discretization error by comparing the discrete and continuous dynamics. The following two bounds from [17], used throughout, produce the log(2d) dependence of the bias. Lemma E.1 (Lemma B.1 in [17]). Suppose Y = (Y (1) , . . . , Y (d) ) ∈ Rd with each Y (i) centered and sub-Gaussian with variance proxy σ 2 , i.e. E[exp(λY (i) )] ≤ exp( 21 λ2 σ 2 ). Then (E.1)

E[|Y |2ℓ∞ ] ≤ 4σ 2 log(2d) .

Lemma E.2 (Proposition 2.3 in [17]). Let Assumption 1.1 hold. Then » » (E.2) Eπ [|∇V |2ℓ∞ ] ≤ 2 β log(2d) . E.1. Discretization error in HMC. The discretization error decompositions (2.6) and (2.7) in Section 2.2 share the common structure (E.3)

1 h ∗ ∗ ∞ |M Π1 Um−i−1 hmc,h (Uhmc,h − Uhmc )(q , p )|2,ℓ ,

with a matrix M ∈ Rd×d and (q ∗ , p∗ ) ∼ π ⊗ N (0, I), where, as recalled in Section 5, (E.4)

(U1hmc,h − Uhhmc )(q ∗ , p∗ ) = ((Q1h − Qh ), (P1h − Ph ))(q ∗ , p∗ ) .

The following proposition bounds the error of each component.

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Proposition E.3 (Generalization of Proposition 5.1). For any matrix M ∈R Rd×d , suppose there existRconstants CR3 , C4 > 0 (possibly depending on β) such that |M ∇2 V dν1 |ℓ∞ 2 2 d ∞ < C3 and |M ( ∇ √ V dν1 )( ∇∗ V ∗dν2 )|ℓ < C4 for all probability measures ν1 , ν2 on R . Then for h ≤ 1/ 20β and (q , p ) ∼ π ⊗ N (0, I), » 1 |M (Q1h (q ∗ , p∗ ) − Qh (q ∗ , p∗ ))|2,ℓ∞ ≤ √ h3 C3 log(2d) , 2

(E.5)

|M (P1h (q ∗ , p∗ ) − Ph (q ∗ , p∗ ))|2,ℓ∞ ≤ 2h2 C3

(E.6)

» » 1 log(2d) + √ h4 C4 log(2d) . 2 5

Proof of Proposition E.3. Let (qt , pt ) = Uthmc (q ∗ , p∗ ). The continuous and discrete position updates have the integral representations Z hZ t h ∗ ∗ ∗ ∗ (E.7) Q (q , p ) = q +p h− ∇V (qs )dsdt , Q1h (q ∗ , p∗ ) = q ∗ +p∗ h− 12 h2 ∇V (q ∗ ) , 0

RhRt whose difference is Q1h (q ∗ , p∗ ) − Qh (q ∗ , p∗ ) = 0 0 (∇V (qs ) − ∇V (q ∗ ))dsdt. Writing R1 ∇V (qs ) − ∇V (q ∗ ) = ( 0 ∇2 V (τ qs + (1 − τ )q ∗ )dτ )(qs − q ∗ ) and using the hypothesis on M, Z hZ t h ∗ ∗ 2 1 ∗ ∗ 1 2 2 (E.8) |M (Qh (q , p ) − Q (q , p ))|ℓ∞ ≤ 2 h C3 |qs − q ∗ |2ℓ∞ dsdt . 0

Expanding qs bounds the remaining integral, (E.9) Z hZ t Z sZ τ Z hZ t Z hZ t 2 ∗ 2 ∗ 2 ∇V (qr )drdτ ∞ dsdt |p s|ℓ∞ dsdt + 2 |qs − q |ℓ∞ dsdt ≤ 2 0

≤ 61 h4 |p∗ |2ℓ∞ +

Z hZ t 0

Z sZ τ 0

ℓ

|∇V (qr )|2ℓ∞ drdτ dsdt .

Since (q ∗ , p∗ ) ∼ π ⊗ N (0, I), stationarity gives qt ∼ π for all t, so by Lemma E.1, Z hZ t (E.10)

E 0

1 6 |qs − q ∗ |2ℓ∞ dsdt ≤ 32 h4 log(2d) + 60 h Eπ [|∇V |2ℓ∞ ] .

Combining these with Lemma E.2 and h2 β ≤ 1/20, (E.11) » » 1 |M (Q1h (q ∗ , p∗ ) − Qh (q ∗ , p∗ ))|2,ℓ∞ ≤ √13 h3 C3 log(2d) + √120 h4 C3 Eπ [|∇V |2ℓ∞ ] » ≤ √12 h3 C3 log(2d) , which is (E.5). For the momentum component, (E.12) P1h (q ∗ , p∗ ) = p∗ − 12 h∇V (q ∗ ) − 12 h∇V (Q1h (q ∗ , p∗ )) ,

Ph (q ∗ , p∗ ) = p∗ −

Z h ∇V (qt )dt , 0

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

so |M (P1h (q ∗ , p∗ ) − Ph (q ∗ , p∗ ))|22,ℓ∞ Z h Z h 2 ≤ 12 hC32 E |qt − q ∗ |2ℓ∞ dt + 12 E M̃ (qt − Q1h (q ∗ , p∗ ))dt ∞ , ℓ {z } {z } | 0 | 0

(E.13)

(a)

where M̃ := give

(b)

2 1 ∗ ∗ 0 M ∇ V (τ qt + (1 − τ )Qh (q , p ))dτ . The bounds used in (E.9) and (E.10) 1 5 (a) ≤ 83 h3 log(2d) + 10 h Eπ [|∇V |2ℓ∞ ] .

(E.14)

Using the position differences in (E.7), (E.15) Z tZ s 1 ∗ ∗ ∗ (∇V (qτ ) − ∇V (q ∗ ))dτ ds + 21 (h2 − t2 )∇V (q ∗ ) qt − Qh (q , p ) = p (t − h) − 0 0 Z tZ sZ 1 =p∗ (t − h) − ∇2 V (rqτ + (1 − r)q ∗ )dr (qτ − q ∗ )dτ ds + 21 (h2 − t2 )∇V (q ∗ ) , 0

so (b) splits as (b) ≤ 3hC32 E |

Z h 0

|p∗ (t − h)|2ℓ∞ dt +3hC32 E | {z }

Z h 0

2 2 ∗ 1 2 2 (h − t )∇V (q ) ℓ∞ dt

(b1 )

(E.16) + 21 h3 C42 E |

Z hZ tZ s 0

(b2 )

}

|qτ − q ∗ |2ℓ∞ dτ dsdt . {z } (b3 )

Direct computation of (b1 ), (b2 ) and the approach of (E.10) for (b3 ) give (E.17) 2 5 2 5 1 7 (b1 ) ≤ 43 h3 log(2d) , (b2 ) ≤ 15 h Eπ [|∇V |2ℓ∞ ] , (b3 ) ≤ 15 h log(2d)+ 420 h Eπ [|∇V |2ℓ∞ ] . Combining these via (E.13) and applying Lemma E.2 gives (E.6).

□

E.2. Discretization error in UL. We adapt the framework of [34], in particular the analysis leading to their Proposition 8.3, to the W2,ℓ∞ norm. We use an HOH chain preserving the invariant distribution exactly, and measure the discrepancy between the BAOAB and HOH chains. Proposition E.4. Consider initial states (x, p) and (y, p′ ) with (x, p) ∼ π ⊗ N (0, I). Apply one HOH step to (x, p) and one BAOAB step to (y, p′ ), coupled through the same Gaussian noise ξ, and set (∆q , ∆p ) := Uξhoh,h (x, p) − Uξbaoab,h (y, p′ ). Suppose π satisfies R R R √ . If | ∇2 V dν1 |ℓ∞ ≤ C5 and |( ∇2 V dν1 )( ∇2 V dν2 )|ℓ∞ ≤ Assumption 1.1 and h ≤ 21−η β C6 for all probability measures ν1 , ν2 on Rd , then 3 (E.18) |∆q |2,ℓ∞ ≤ (1+ 41 h2 (η +1)C5 )|x−y|2,ℓ∞ + h2 (η +1)|p−p′ |2,ℓ∞ + 27 64 C5 h 2

» log(2d) ,

(E.19) ∆p = (ηI − h (η+1) Q2 )(p − p′ ) + (− h2 ηQ1 − h2 Q2 + h (η+1) Q2 Q1 )(x − y) + εp , 4 8

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

p p 3 4 d×d satisfy αI ⪯ 2 where |εp |2,ℓ∞ ≤ 27 8 h C5 log(2d) + 8 h C6 log(2d), and Q1 , Q2 ∈ R Q1 , Q2 ⪯ βI, |Q1 |ℓ∞ ≤ C5 , |Q2 |ℓ∞ ≤ C5 , and |Q2 Q1 |ℓ∞ ≤ C6 . Proof of Proposition E.4. We use integral representations of the updates. Let Uξho,h h/2

be the composition of a half-step of Hamiltonian dynamics H = Uhmc , followed by a full Ornstein–Uhlenbeck map O with noise ξ. We write (x̄, p̄) := Uξho,h (x, p). With (x(t), p(t)) = Uthmc (x, p) the continuous Hamiltonian trajectory (1.1), we have (x(t), p(t)) ∼ π ⊗ N (0, I) and (E.20) Z h/2 Z h/2 p h h x̄ = x + 2 p − ∇V (x(t))( 2 − t)dt , p̄ = η p − ∇V (x(t))dt + 1 − η 2 ξ . 0

Position component. Let (x̄(t), p̄(t)) = Uthmc (x̄, p̄), so (x̄(t), p̄(t)) ∼ π ⊗ N (0, I). The position updates are

qHOH :=x + h2 p(1 + η) − (E.21)

Z h/2 0

∇V (x(t))( h2 − t)dt −

Z h/2 0

∇V (x̄(t))( h2 − t)dt

Z h/2 p h h ∇2 V (x(t))p(t)( h2 − t)dt + 1 − η 2 ξ , + 2 η − 2 ∇V (x) − 0

qBAOAB := y + h2 p′ (1 + η) − h4 (1 + η)∇V (y) + h2

(E.22)

1 − η2ξ ,

and their difference satisfies (E.23) |∆q |2,ℓ∞ 2 ≤(1 + h4 ηC5 )|x − y|2,ℓ∞ + h2 (1 + η)|p − p′ |2,ℓ∞ +

+ η h2

Z h/2 0

∇ V (x(t))p(t)( h2 − t)dt +

Z h/2 0

(∇V (x(t)) − ∇V (y))( h2 − t)dt

(∇V (x̄(t)) − ∇V (y))( h2 − t)dt

2,ℓ∞

The last term is at most Z h/2 Z h/2 hηC5 h (E.24) C5 (|x(t) − y|2,ℓ∞ + |x̄(t) − y|2,ℓ∞ )( 2 − t)dt + 2 |p(t)|2,ℓ∞ ( h2 − t)dt . 0

Expanding the integral representation of x(t) gives |x(t) − y|2,ℓ∞ ≤ |x − y|2,ℓ∞ + |tp − Rt ∞ 0 ∇V (x(s))(t − s)ds|2,ℓ , and by Lemmas E.1 and E.2, Z t » » 2 (E.25) |tp − ∇V (x(s))(t − s)ds|2,ℓ∞ ≤ h log(2d) + h4 β log(2d) . 0

Similarly, with x̄ = x(h/2), |x̄(t) − y|2,ℓ∞ ≤ |x̄(t) − x̄|2,ℓ∞ + |x̄ − y|2,ℓ∞ , where » » 2 |x̄(t) − x̄|2,ℓ∞ ≤ h log(2d) + h4 β log(2d) , (E.26) » » 2 |x̄ − y|2,ℓ∞ ≤ |x − y|2,ℓ∞ + h log(2d) + h4 β log(2d) . √ Combining these and applying h ≤ 1/(2 β) gives (E.18).

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Momentum component. The momentum updates are Z h/2 p 2 h h := ∇ V (x(t))p(t)( 2 − t)dt + 1 − η 2 ξ pHOH η p − 2 ∇V (x) − 0 (E.27) Z h/2 − h2 ∇V (x̄) −

∇2 V (x̄(t))p̄(t)( h2 − t)dt ,

p pBAOAB := η p′ − h2 ∇V (y) + 1 − η 2 ξ − h2 ∇V (ŷ) , p 2 where ŷ = y + h2 (η + 1)p′ − h4 (η + 1)∇V (y) + h2 1 − η 2 ξ, so Z h/2 ′ h ∆p = η(p − p ) − 2 η(∇V (x) − ∇V (y)) − η ∇2 V (x(t))p(t)( h2 − t)dt 0 (E.29) Z h/2

(E.28)

− h2 (∇V (x̄) − ∇V (ŷ)) −

∇2 V (x̄(t))p̄(t)( h2 − t)dt .

Grouping the principal terms into αp := η(p − p′ ) − h2 η(∇V (x) − ∇V (y)) − h2 (∇V (x̂c ) − ∇V (ŷ)) , p 2 with x̂c = x + h2 (η + 1)p − h4 (η + 1)∇V (x) + h2 1 − η 2 ξ, (E.30)

αp = (ηI − h (η+1) Q2 )(p − p′ ) + (− h2 ηQ1 − h2 Q2 + h (η+1) Q2 Q1 )(x − y) , 4 8 R 1 where Q1 = ∇2 V ([y, x]), Q2 = ∇2 V ([x̂c , ŷ]), and ∇2 V ([v1 , v2 ]) := 0 ∇2 V (v1 + s(v2 − v1 ))ds. The remainder εp := ∆p − αp is I1 + I2 + I3 + I4 , with Z h/2 Z h/2 2 h h2 2 := I1 −η ∇ V (x(t))p( 2 − t)dt + η 4 ∇ V ([x̂c , x̄])p − η ∇2 V (x̄(t))p̃( h2 − t)dt (E.31)

Z h/2 p p 2 ∇2 V (x̄(t))ξ( h2 − t)dt , + h4 ∇2 V ([x̂c , x]) 1 − η 2 ξ − 1 − η 2 0

R h/2

where p̃ := p − h2 ∇V (x) − 0

∇2 V (x(t))p(t)( h2 − t)dt;

Z h/2 I2 := η 0

∇2 V (x(t))(p(t) − p)( h2 − t)dt ,

2 I3 := h2 ∇2 V ([x̂c , x̄])(− h4 (η + 1)∇V (x) +

Z h/2 I4 := − 0

∇2 V (x̄(t))(−t∇V (x̄) −

Z t 0

Z h/2 0

∇V (x(t))( h2 − t)dt) ,

∇2 V (x̄(s))p̄(s)(t − s)ds)( h2 − t)dt .

These terms are bounded by (E.32) » » 2 1 2 1 4 1 4 ∞ |I1 |2,ℓ∞ ≤ ( 21 h C + h C ) log(2d) , |I | ≤ ( h C + h C ) log(2d) , 5 6 2 2,ℓ 5 6 8 8 4 8 » » |I3 |2,ℓ∞ ≤ 38 h2 C5 log(2d) , |I4 |2,ℓ∞ ≤ ( 81 h2 C5 + 81 h4 C6 ) log(2d) , √ using h ≤ 1/(2 β), and combining them gives (E.19). □

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

|(∆lq , ∆lp )|2,ℓ∞ = (∆lq + b∆lp , w

p a − b2 ∆lp )

2,ℓ∞

≤ |∆lq |2,ℓ∞ +

√

a|∆lp |2,ℓ∞ .

To bound this surrogate, one could combine the single-step error of Proposition E.4 with the joint-space contraction of Corollary D.3. This is too crude, as it treats the position and momentum errors together. Instead, we use a Grönwall argument that tracks the two separately, mirroring the ℓ2 analysis of [34]. Proposition E.5 (Refined statement of Proposition 5.2). For the (∆lq , ∆lp ) process with √ , shared initial state (q ∗ , p∗ ) ∼ π ⊗ N (0, I), if π satisfies Assumption 1.1 and h ≤ 21−η β then the accumulated discretization error over l steps satisfies (E.35) h p i √ l l C5 3√ 1 C6 ∞ ≤ exp (l − 1)h ∞ + a|∆ | + + β El , |(∆lq , ∆lp )|2,ℓ∞ ≤ |∆ | 2,ℓ 2,ℓ 3/2 p q 2 β 16 β w where C5 , C6 are the constants in Proposition E.4 and » 2 4 27 √ h C5 + 3 √ h C6 C5 h3 l + 27 log(2d) . (E.36) El = 64 8 β(1−η) 8 β(1−η) Proof of Proposition E.5. By Proposition E.4, unrolling the position error recurrence to ∆0q = 0 gives » 3 l−1 l−1 27 ∞ ∞ ∞ |∆lq |2,ℓ∞ ≤ |∆l−1 | + h(hC |∆ | + |∆ | ) + C h log(2d) 5 2,ℓ 2,ℓ 2,ℓ q q p 64 5 (E.37) ≤

l−1 X

3 h(hC5 |∆iq |2,ℓ∞ + |∆ip |2,ℓ∞ ) + 27 64 C5 lh

» log(2d) .

i=1

Likewise, unrolling the momentum error recurrence and bounding the geometric sum of the constant term by (1 − η)−1 yields (E.38) 2 l−1 h3 |∆lp |2,ℓ∞ ≤ (η + h2 C5 )|∆l−1 p |2,ℓ∞ + (hC5 + 4 C6 )|∆q |2,ℓ∞ + |εp |2,ℓ∞ ≤

l−1 X

2 η i−1 h2 C5 |∆l−i p |2,ℓ∞ +

i=1

l−1 X

1 η i−1 (hC5 + h 4C6 )|∆l−i q |2,ℓ∞ + 1−η |εp |2,ℓ∞ .

i=1

√

Combining the two, with h ≤ 1/(2 β) and η ≤ 1, we obtain (E.39) l−1 X p √ √ l l l−i C5 1 C6 |∆q |2,ℓ∞ + a|∆p |2,ℓ∞ ≤ h 32 √ + + β (|∆l−i 3/2 q |2,ℓ∞ + a|∆p |2,ℓ∞ ) + El . 16 β β i=1

Applying Grönwall’s inequality leads to (E.35).

□

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

Appendix F. Sampling Bias Bounds Under the framework of Section 2, we now combine the propagator norm estimates with the discretization error to track the sampling error across steps, and obtain the bias bounds in the limit of infinitely many steps. F.1. HMC. We prove the iterative W2,ℓ∞ bounds of Propositions 2.1 and 2.2, from which the HMC cases of Theorems 1.4 and 1.6 follow by iterating to the stationary limit. We begin with the weak interaction case. Proof of Proposition 2.1. Following the coupling and one-step split (2.4) into the contraction term (a) and discretization term (b), we bound each. Contraction term. By Proposition 3.5, 2

1 (a) = |T̃m ({I − h2 Hi }m−1 i=0 )(Xk − Yk )|2,ℓ∞ ,

(F.1)

R1 with Hi = 0 ∇2 V (τ xi + (1 − τ )yi )dτ and xi , yi the internal positions after i leap-frog steps, xi := Π1 Uihmc,h (Xk , ξk ), yi := Π1 Uihmc,h (Yk , ξk ). The contractive ℓ∞ bound of Proposition 4.3 gives ã Å α |Xk − Yk |2,ℓ∞ . (F.2) (a) ≤ 1 − 400β Discretization error. The telescoping decomposition (2.6) and the component split (5.1) express (b) through single integration-step errors. Representing each propagated factor as a matrix polynomial via Proposition 3.5, with (q ∗ , p∗ ) ∼ π ⊗ N (0, I), (b) ≤

m−1 X

1 )(Q1h − Qh )(q ∗ , p∗ )|2,ℓ∞ |T̃m−i−1 ({I − h2 Hj }m−i−2 j=0 | {z } i=0 (b1 )

(F.3) +

m−1 X

1 |hŨm−i−2 ({I − h2 Hj }m−i−2 )(P1h − Ph )(q ∗ , p∗ )|2,ℓ∞ , j=1 | {z } i=0 (b2 )

where Hj = 0 ∇2 V (τ xj +(1−τ )yj )dτ and xj , yj are the j-th internal leap-frog positions started at U1hmc,h (q ∗ , p∗ ) and Uhhmc (q ∗ , p∗ ), i.e. we redefine xj := Π1 Ujhmc,h U1hmc,h (q ∗ , p∗ ), α 1 yj := Π1 Ujhmc,h Uhhmc (q ∗ , p∗ ). By Proposition 4.3, almost surely |T̃m−i−1 |ℓ∞ ≤ 1 + 400β ≤ 401 3 1 √ ∞ 400 and |hŨm−i−2 |ℓ ≤ 2 β .

To bound (b1 ), (b2 ) we apply Proposition E.3, instantiating its general matrix M to the specific polynomials appearing here and determining the corresponding constants 401 1 C3 , C4 . For (b1 ), we take M = T̃m−i−1 , which satisfies |M |ℓ∞ ≤ 400 . Under the weak R 2 α 51 interaction assumption the integrated Hessian obeys | ∇ V dν|ℓ∞ ≤ β + 50 ≤ 50 β, so √ √ R 2 1 |M ∇ V dν|ℓ∞ ≤ 2β and we may take C3 = 2β. For (b2 ), we take M = hŨm−i−2 153 √ with |M |ℓ∞ ≤ 2√3 β , giving C3 = 2√3 β · 51 50 β = 100 β and, for a product of two Hessians, √ 3 51 2 32 2 C4 = 2√3 β · ( 51 50 β) = 2 ( 50 ) β . With h ≤ 1/ 20β, Proposition E.3 yields » » » (F.4) (b1 ) ≤ h3 β log(2d) ≤ h2 β log(2d) , (b2 ) ≤ 4h2 β log(2d) ,

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

and summing over the m steps, (b) ≤ 5mh2

(F.5)

» » β log(2d) ≤ 2h log(2d) ,

√ using mh = 1/ 20β in the last step. Thus, Å ã » α (F.6) |Xk+1 − Yk+1 |2,ℓ∞ ≤ 1 − |Xk − Yk |2,ℓ∞ + 2h log(2d) . 400β Taking Xk , Yk optimally coupled so that W2,ℓ∞ (ρk , π) = |Xk − Yk |2,ℓ∞ , we obtain ã Å » α W2,ℓ∞ (ρk , π) + 2h log(2d) . (F.7) W2,ℓ∞ (ρk+1 , π) ≤ 1 − 400β □ √ Proof of the HMC case of Theorem 1.4. Under C (O) = 1/50 and mh = 1/ 20β, taking k → ∞ in the inequality of Proposition 2.1 gives ã Å » β (F.8) W2,ℓ∞ (πh , π) = O h log(2d) . α □ We turn to the sparse interaction case. Proof of Proposition 2.2. With Xk , Yk as before, the multi-step split (2.5) bounds |Xk+N − Yk+N |2,ℓ∞ by its contraction term (a) and discretization term (b). Contraction term. Proposition 3.5 represents the propagators as matrix polynomials, bounded in ℓ2 by Proposition 4.1 and in ℓ∞ by Proposition C.3, which extends √ 2 ∞ Proposition 4.5. The ℓ bound, converted to ℓ at a cost of d, gives √ α N ) (F.9) (a) ≤ (1 − 200β d |Xk − Yk |2,ℓ∞ . Discretization error. The telescoping decomposition (2.7) and the component split (5.1), with the matrix-polynomial representation of Proposition 3.5, give (b) ≤

N −1 m−1 X X

A(N −l−1)m+(m−i−1) (Q1h − Qh )(q ∗ , p∗ ) ∞ ℓ {z } l=0 i=0 | (b1 )

(F.10) +

N −1 m−1 X X

B(N −l−1)m+(m−i−2) (P1h − Ph )(q ∗ , p∗ ) ∞ , ℓ {z } l=0 i=0 | (b2 )

where, writing Pl,i :=

Q←−

−l−2 k∈{m−i−1+jm}N j=0

m−1 1 ({I − h H T̃m k+s }s=0 ) for the product of 2

full-loop propagators over the N − l − 1 subsequent outer loops, (F.11)

1 A(N −l−1)m+(m−i−1) = Pl,i T̃m−i−1 ,

1 B(N −l−1)m+(m−i−2) = Pl,i hŨm−i−2 .

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

The Hessians Hj = step count:

2 0 ∇ V (τ xj +(1−τ )yj )dτ are indexed by the total internal leap-frog

0≤j ≤m−i−1:

xj = Π1 Ujhmc,h U1hmc,h (q ∗ , p∗ ) , yj = Π1 Ujhmc,h Uhhmc (q ∗ , p∗ ) ;

0 ≤ ˆl ≤ N − l − 1, 0 ≤ j < m − 1 : xl̂m+j+(m−i−1) = Φjh,ξ

k+l+l̂+1

yl̂m+j+(m−i−1) = Φjh,ξ

k+l+l̂+1

k+l+l̂

m−i−1 1 ∗ ∗ · · · Φm h,ξk+l+1 Π1 Uhmc,h Uhmc,h (q , p ) ,

k+l+l̂

m−i−1 h ∗ ∗ · · · Φm h,ξk+l+1 Π1 Uhmc,h Uhmc (q , p ) ,

Φm h,ξ Φm h,ξ

with the convention that for ˆl = 0 only Φjh,ξ is applied. By Propositions C.3 k+l+l̂+1 and E.3, for r = r(N −l−1)m+(m−i−1) , » » (F.12) (b1 ) ≤ 3h2 sr β log(2d) , (b2 ) ≤ 9h2 sr β log(2d) , so summing over the N m steps, (b) ≤ 12h2

(F.13)

β log(2d)

NX m−1

√

sri .

i=0

Combining (F.9) and (F.13) with Xk , Yk optimally coupled completes the proof. □ The following corollary follows by specifying N and estimating polynomial growth condition.

P√

sri under the

√ log(2 d)⌉, so that the conProof of the HMC case of Theorem 1.6. Take N = ⌈ 200β α √ α N traction factor (1 − 200β ) d ≤ 12 . Under sk ≤ C(k + 1)n , an integral comparison gives (F.14) NX m−1 NX m−1 √ p n √ √ 2 log d sri ≤ C ih βe + log(5/3) +2 i=0

i=0

n +1 √ √ 2 log d Z Nm √ p n √ N mh βe + + 2 √ log(5/3) 2 log d √ ≤ C yh βe + log(5/3) + 2 dy ≤ C . (n/2 + 1)h βe 0 √ With this N and mh ≤ 1/ 20β, Proposition 2.2 gives the recursion » Ä ä n +1 (F.15) W2,ℓ∞ (ρk+N , π) ≤ 12 W2,ℓ∞ (ρk , π) + h log(2d) O( αβ log(2d)) 2 , and taking k → ∞ gives (F.16)

Å Å ãã n +1 » 2 β W2,ℓ∞ (πh , π) = h log(2d) O log(2d) . α □

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

F.2. UL. Proof of the UL case of Theorem 1.4. We follow the coupling framework (2.12) with the weighted matrix representation of Section 4.1. Denote the initial distribution of X0 by ρ0 . The BAOAB chain (Xk , Pk ) from (X0 , P0 ) ∼ ρ0 ⊗ N (0, I) and the HOH chain (Yk , Pk′ ) from (Y0 , P0′ ) ∼ π ⊗ N (0, I) are coupled through the same noise {ξi }. Introducξk−1 0 (Y0 , P0′ ), which starts at (Y0 , P0′ ) · · · Uξbaoab,h ing the auxiliary process (Ȳk , P̄k′ ) = Ubaoab,h but runs BAOAB with the same noise, the split (2.12) becomes |Xk − Yk |2,ℓ∞ ≤ |Xk − Ȳk |2,ℓ∞ + |Ȳk − Yk |2,ℓ∞ , | {z } | {z }

(F.17)

(a)

(b)

the contraction between two BAOAB chains (a) and the BAOAB-vs-HOH discretization error (b). ξk−1 ξk−1 0 ··· = Uab,h Uabao,h Contraction term. Using the decomposition Ubaoab,h · · · Uξbaoab,h

1 0 Uξabao,h Uξbao,h of Section 2.3 and the weighted position propagators (4.3), w w w w (F.18) Xk − Ȳk = I 0 W −1 Mab (Hk )Mabao (Hk−1 ) · · · Mabao (H1 )Mbao (H0 )W ∆Z0 ,

with ∆Z0 = (X0 − Y0 , P0 − P0′ ) and the matrices defined through the intermediate iterates 0 z1 := Uξbao,h (X0 , P0 ) ,

0 z1′ := Uξbao,h (Y0 , P0′ ) ,

i−1 1 · · · Uξabao,h z1 , zi := Uabao,h

zk+1 := Uab,h zk ,

i−1 1 · · · Uξabao,h z1′ , zi′ := Uabao,h

′ := Uab,h zk′ , zk+1

2 ≤ i ≤ k,

satisfying w W (z1 − z1′ ) = Mbao (H0 )W ∆Z0 ,

w ′ ) = Mab (Hk )W (zk − zk′ ) , W (zk+1 − zk+1

w ′ W (zi − zi′ ) = Mabao (Hi−1 )W (zi−1 − zi−1 ), 2 ≤ i ≤ k. −1 b By Proposition 4.4 and | I 0 W |ℓ∞ = 1 + √a−b , 2 b (1 − 12 c(h))k−1 |∆Z0 |2,ℓ∞ . (F.19) (a) ≤ 6 1 + √a−b 2 w

Discretization error. Write (∆lq , ∆lp ) = (Ulbaoab,h − Ulhoh,h )(q ∗ , p∗ ). The matrix k−(i+1)l̃

representation of Π1 Ubaoab,h (Ul̃baoab,h − Ul̃hoh,h )(q ∗ , p∗ ) parallels (F.18), with ∆Z0 = (∆l̃q , ∆l̃p ) and Hessians along the intermediate iterates. Each term of the telescoping decomposition (2.13) is thus bounded by k−(i+1)l̃

(F.20) (F.21)

|Π1 Ubaoab,h (Ul̃baoab,h − Ul̃hoh,h )(q ∗ , p∗ )|2,ℓ∞ b ≤ 6 1 + √a−b , (1 − 12 c(h))k−(i+1)l̃−1 |(∆l̃q , ∆l̃p )|2,ℓ∞ 2 w k−⌊k/l̃⌋l̃

k−⌊k/l̃⌋l̃

Π1 Ubaoab,h − Uhoh,h

(q ∗ , p∗ ) 2,ℓ

l̃⌋l̃ l̃⌋l̃ √ b ≤ 1 + |(∆k−⌊k/ , ∆k−⌊k/ )|2,ℓ∞ . q p w a−b2 ∞

In (F.20), when i = ⌊k/˜l⌋ − 1 there may be no internal ABAO steps and the exponent b can be negative, so we instead bound the operator by 6(1 + √a−b ), accounting for the 2

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

BAO and AB prefactors. Proposition E.5 gives for any l ≤ ˜l h p i C5 1 C6 ˜l − 1)h 3 √ (F.22) |(∆lq , ∆lp )|2,ℓ∞ ≤ exp ( β El̃ . + + 2 β 16 β 3/2 w Substituting these into (F.17), with X0 , Y0 optimally coupled so W2,ℓ∞ (ρ0 , π) = |X0 − Y0 |2,ℓ∞ and P0′ = P0 , (F.23) k−1 b 1 b ∞ (ρ0 , π) + 1 + √ |Xk − Yk |2,ℓ∞ ≤ 6 1 + √a−b 1 − c(h) W 2,ℓ 2 2 a−b2 Ñ é ⌊k/l̃⌋−2 h X p i C5 1 C6 + · 7+ 6(1 − 21 c(h))k−(i+1)l̃−1 exp (˜l − 1)h 23 √ + β El̃ . 16 β 3/2 β i=0 (O)

Under the weak interaction assumption |Hi |ℓ∞ ≤ α/20, the constants of Proposi21 2 2 2 tion E.4 satisfy C5 ≤ 21 20 β and C6 ≤ ( 20 ) β√, so that C5 /β and C6 /β are O(1). Since √ (˜l − 1)h β = O(1) by the choice ˜l = ⌈1/(2 βh)⌉, the exponential is O(1), and substituting El̃ gives h2 √β » h p i C5 3√ 1 C6 ˜ log(2d) . exp (l − 1)h 2 β + 16 β 3/2 + β El̃ = O 1−η 2

αh As k → ∞, the first term vanishes. For the geometric sum, with c(h) = 4(1−η) and √ ˜l ≥ 1/(2 βh), the inequality 1 − (1 − x)l̃ ≥ ˜lx − 1 (˜lx)2 for x > 0 gives 2

αh 63 √ 1 − (1 − 21 c(h))l̃ ≥ , 1024 β(1 − η) √ P⌊k/l̃⌋−2 so that i=0 6(1 − 12 c(h))k−(i+1)l̃−1 ≤ 6/(1 − (1 − 21 c(h))l̃ ) = O( β(1−η) ). Since αh b 1 + √a−b2 = O(1), (F.23) becomes (F.24) β » h2 √β » √β(1 − η) log(2d) = O h log(2d) . ·O W2,ℓ∞ (πh , π) ≤ O(1) · O αh 1−η α □

Proof of the UL case of Theorem 1.6. As in the HMC sparse case, we induce an ℓ∞ contraction from the ℓ2 contraction over multiple steps. We partition the k steps into blocks of B steps, with B specified below. We use the three processes (Xk , Pk ), (Yk , Pk′ ), (Ȳk , P̄k′ ) of the proof of the UL case of Theorem 1.4 with P0 = P0′ = P̄0′ , and the decomposition (F.25)

|Xk − Yk |2,ℓ∞ ≤ |Xk − Ȳk |2,ℓ∞ + |Ȳk − Yk |2,ℓ∞ . | {z } | {z } (a)

(b)

Contraction term. Corollary D.3 gives √ k−1 (F.26) (a) ≤ 14 d(1 − c(h)) 2 |X0 − Y0 |2,ℓ∞ , which vanishes as k → ∞. Discretization error. The telescoping decomposition (2.13) accumulates a local error every ˜l = ⌈ 2√1βh ⌉ steps, each propagated through the later steps. Write (∆lq , ∆lp ) = (Ulbaoab,h − Ulhoh,h )(q ∗ , p∗ ) for an l-step error, with (q ∗ , p∗ ) ∼ π ⊗ N (0, I), and recall

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

that Ulbaoab,h , Ulhoh,h are the l-fold compositions omitting the coupled noise. For the √ √ sparse case the constants of Proposition E.4 are C5 = β s1 and C6 = β 2 s2 , so Proposition E.5 bounds the ˜l-step error, (F.27) h2 √β » h p i C5 1 C6 ˜l − 1)h 3 √ β E = O log(2d) . + |(∆l̃q , ∆l̃p )|2,ℓ∞ + ≤ exp ( l̃ 2 β 16 β 3/2 w 1−η Setting the block length B := K1 ˜l with K1 specified later and grouping the trailing steps in (2.13) into full blocks of B steps followed by residual steps, ⌊k/l̃⌋ ⌉ K −1 1 1

⌈ K

(F.28)

(b) ≤

X X i=0

l̃ l̃ ∗ ∗ 1 UiB Π1 Ubaoab,h baoab,h (Ubaoab,h − Uhoh,h )(q , p )

k1 =0

k−⌊k/l̃⌋l̃ k−⌊k/l̃⌋l̃ + Π1 Ubaoab,h − Uhoh,h (q ∗ , p∗ )

2,ℓ∞

where i counts the full B-step blocks, k1 ∈ {0, . . . , K1 − 1} indexes the residual ˜l-step groups, and jk1 := k − ⌊k/˜l⌋˜l + k1 ˜l is the corresponding number of trailing steps. Each summand in the first sum is one ˜l-step error propagated through iB full-block steps and jk1 trailing steps. For two initial states (x, p), (y, p′ ), Proposition 3.5 and (4.7) give the weighted matrix-polynomial form of the propagated position difference over the residual steps, (F.29) î η ó ïx − y ò jk1 jk1 η η b 1 ′ √ √ Π1 Ubaoab,h (x, p) − Π1 Ubaoab,h (y, p ) = T̃jk − a−b2 T̃jk + a−b2 hŨjk −1 W , 1 1 1 p − p′ 2

the polynomials evaluated at matrices of the form I − h2 H with H the integrated Hessians along the intermediate iterates of the two trajectories. Applying this with ∗ ∗ l̃ ∗ ∗ ′ iB l̃ (x, p) = UiB baoab,h Ubaoab,h (q , p ) and (y, p ) = Ubaoab,h Uhoh,h (q , p ), the iB full-block steps contract the ˜l-step error by Corollary D.3, √ iB−1 (F.30) |(x − y, p − p′ )|2,ℓ∞ ≤ 7 2d(1 − c(h)) 2 |(∆l̃q , ∆l̃p )|2,ℓ∞ , w w while Proposition 4.6 bounds the ℓ∞ norm of the matrix-polynomial row in (F.29). Together they give jk

(F.31)

l̃ l̃ ∗ ∗ 1 Π1 Ubaoab,h UiB baoab,h (Ubaoab,h − Uhoh,h )(q , p ) 2,ℓ∞ » √ √ iB−1 srjk (1 − c(h)) 2 |(∆l̃q , ∆l̃p )|2,ℓ∞ . ≤ 7 2d 14 2 + 1 + √23 w 1

For i = 0 where there are no intermediate full blocks we instead employ (F.32) Ä √ ä» jk1 Π1 Ubaoab,h (Ul̃baoab,h − Ul̃hoh,h )(q ∗ , p∗ ) ∞ ≤ 14 2 + 1 + √23 srjk |(∆l̃q , ∆l̃p )|2,ℓ∞ . w

1 2,ℓ √ √ We now fix K1 = ⌈ 16 β(1−η)hαlog(14 2d) + 1⌉, so that the per-block factor contracts, √ B−1 7 2d(1 − c(h)) 2 ≤ 12 , making the contraction over i blocks at most 2−i . Substituting

DELOCALIZATION OF BIAS IN UNADJUSTED HMC AND UNDERDAMPED LANGEVIN

the ˜l-step error bound and these per-term estimates into (F.28), together with the last-term bound (F.21), and summing the geometric factor 2−i over i, we obtain Ñ é 1 −1 h2 √β » Ä √ ä KX » 2 √ (F.33) (b) ≤ 1 + 2 14 2 + 1 + 3 srjk O log(2d) . 1 1−η k1 =0

By an integral comparison as in (F.14), with sri increasing in i, 2 n +1 √ 2 l̃e β(K1 +1) 2 K KX 1 1 −1 h + log d + 2 X» √ » 1−η srjk ≤ srk l̃ ≤ C 2β 1 1 2 ( n2 + 1) l̃e (F.34) k1 =1 k1 =0 1−η h ä n +1 1−η Ä β = √ O( α log(2d)) 2 . βh Combining with (F.26) and (F.33) in (F.25) and taking k → ∞, » Ä ä n +1 (F.35) W2,ℓ∞ (πh , π) = h log(2d) O( αβ log(2d)) 2 . □

Record · ID 373375 · SHA-256 c6afece5926e0a46
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.