Overdamped Langevin dynamics (1.1) dX t = -∇V (X t )dt + √ 2dB t , has been used extensively to sample from high dimensional probability distributions in applications ranging from molecular dynamics to Bayesian inverse problems and data assimilation [33]. Here V : R d → R is a function in R d referred to as the potential and B t is the d-dimensional Brownian motion. The target distribution π is proportional to exp(-V ), which is the stationary distribution of (1.1). The overdamped Langevin Monte Carlo algorithm, also known as the unadjusted Langevin algorithm, is obtained by applying the Euler-Maruyama scheme to (1.1):
(1.2)
where h is the step size. The distribution of X kh is denoted by ρ kh .
1.1. Motivation. The unadjusted Langevin algorithm (1.2) is biased, meaning that even when k → ∞, the distribution ρ kh will not converge to the exact target distribution π for any finite step size h. For strongly log-concave distributions, ρ kh converges in the W 2 metric to the stationary distribution of (1.2), π h , at a rate that is independent of dimension (see, e.g., [17]). Consequently, the dependence of the error W 2 (ρ kh , π) on dimension, d, is completely determined by the bias W 2 (π h , π). This bias depends on the step size h, which, according to state-of-the-art analyses, needs to scale proportionally to 1/d or 1/ √ d to achieve a bounded bias for any d. Since the iteration complexity typically scales as 1/h (up to logarithmic factors), this small step size implies a computational cost of order √ d or d. Unbiased samplers based on (1.1), such as the Metropolis-adjusted Langevin algorithm [37,36] and proximal samplers [26,5], are available, but these algorithms also require a small step size when d is large in order to ensure that the acceptance rate is non-negligible; see the review in Section 1.2. In short, existing theoretical analyses of both the unadjusted Langevin algorithm and unbiased variants require step sizes scaling as d -c for some c > 0. In the former case, small step sizes are required to ensure small bias; in the latter, small step sizes are required to maintain a reasonable acceptance rate.
In the case of the unadjusted Langevin algorithm, however, these theoretical predictions seem at odds with abundant empirical evidence that the scheme samples efficiently in extremely high dimensions. As an example, molecular dynamics simulations using integrators closely related to (1.2) with up to billions of atoms are not uncommon [22]. In these simulations, the step size is typically set to several femtoseconds, irrespective of the system size [28].
In this paper, we argue that the mismatch between theory and practice is due, in large part, to the metric used to measure convergence. Standard metrics such as W 2 measure accuracy of the entire distribution, but practitioners are often interested in averages of functions involving a relatively small set of variables. Based on this observation, we aim to investigate the convergence behavior of the algorithm under alternative metrics designed to characterize the accuracy of low-dimensional marginals. We identify a new "delocalization" phenomenon for the bias of the unadjusted Langevin algorithm to show that, even in high-dimensional settings, the step size may not need to be very small if the quantities of interest depend on low-dimensional marginals only. Consequently, the iteration complexity required to achieve a bounded error exhibits a benign dependence on dimension d.
1.2. Literature review. We begin by reviewing existing works on the analysis of Langevin algorithms, with a particular focus on the dependence on dimension d. Analysis of (1.2) dates back to the work of Roberts and Tweedie [36], where asymptotic properties such as ergodicity of the discrete Markov chain are studied. Asymptotic bias of the discrete SDEs can be investigated using Taylor's expansion and the Poisson equation [39,40,1,31]. Non-asymptotic analysis of the algorithm (1.2) aims to characterize the step size and iteration complexity to achieve bounded error under certain metrics; again, the iteration complexity typically scales inversely with the step size (up to logarithmic factors of d). Here, we focus on the scaling of the step size with dimension. A large body of work considers strongly log-concave distributions: Assumption 1.1. Let π ∝ exp(-V ) and V ∈ C 2 (R d ). Assume V is α-strongly convex and β-smooth such that αI ≼ ∇ 2 V (x) ≼ βI at any x ∈ R d . Here 0 < α ≤ β < ∞.
Under Assumption 1.1, existing non-asymptotic bounds for the unadjusted Langevin algorithm (1.2) mainly focus on metrics such as the total variation distance, the W 2 distance, and the Kullback-Leibler (KL) divergence; see [13,12,16,17,8,15]. In these studies, to maintain a bounded error or bias as the dimension d increases, the step size h must scale inversely with d, i.e., h ∼ 1/d. With additional smoothness assumptions on the Hessian of V , the bound on the step size can be improved to h ∼ 1/ √ d; see [29] and the mean squared analysis framework developed in [30]. Other convergence results exist, such as those using the χ 2 and Rényi divergences as metrics [41,20], under weaker assumptions like the log-Sobolev inequality [41,9] or other conditions on π [4,19,32,27], and alternative discretization schemes [38,23]. All these results lead to power law scaling of the step size with d (up to logarithmic terms) to attain a desired error for any d.
Metropolis-adjusted Langevin algorithms (MALA) [36] and proximal samplers [26] are two common unbiased schemes for the Langevin dynamics (1.1). For MALA, the asymptotic analysis in [35] shows that, when h ∼ 1/d 1/3 , the algorithm admits a nondegenerate diffusion limit as d → ∞, provided π is a smooth product measure and the algorithm starts at stationarity; see also further generalizations beyond product measures and stationarity [3,11,34,24] and note that out of stationarity, the results suggest the scaling h ∼ 1/ √ d. Non-asymptotic bounds on the mixing time of MALA have also been investigated in the literature [18,6,10,43,7]. These bounds require a specific choice of step size, scaling as h ∼ 1/ √ d for MALA initialized at a warm start and h ∼ 1/d for a feasible start, up to log d terms; see also related studies on the lower bounds of the mixing time [10,25,43]. For proximal samplers [26], a key component is the implementation of the restricted Gaussian oracle (RGO), and the step size h needs to be small to ensure efficient implementation. Existing analyses suggest choosing h ∼ 1/d if the RGO is implemented via rejection sampling [5], and h ∼ 1/ √ d if it is implemented via approximate rejection sampling [21].
In summary, for both MALA and proximal samplers, choosing the step size inversely proportional to d or √ d is necessary for efficient implementation and to achieve the desired mixing time bounds. For the unadjusted Langevin algorithm, similar scaling is required to attain a bounded error on the bias under the aforementioned metrics.
Under this context, the goal of this paper is to show that the requirement on the step size in the unadjusted Langevin algorithm can be significantly improved if we use an alternative W 2,ℓ ∞ metric to measure convergence. The convergence bounds we prove in this metric indicate that, in some situations, the bias of the unajusted Langevin algorithm is delocalized, in the sense that the bias in individual coordinates is nearly dimension-free. This is precisely the behavior observed when π is a product measure. Our results show that this phenomenon holds in wider generality.
References [2] and [14] studied target distributions for which larger stepsizes are possible in high-dimensions if the observable of interest has Lipschitz constant that scales as 1/ √ d, such as the averaged observable
In these cases, a constant step size is sufficient to achieve a bounded error regardless of dimensionality. Our results go beyond these studies by establishing bounds that are not limited to averaged observables. Our analysis implies that the algorithm could achieve a bounded error for low-dimensional marginals in any dimension using a step size that is dimension-independent up to logarithmic terms. 1.3. Main results. In this section, we outline our main results, including a novel metric W 2,ℓ ∞ to measure convergence, two positive and one negative example that illustrate the delocalization effects of bias, and our main theorem regarding the convergence result for strongly log-concave distributions with sparse interactions.
A new metric to measure convergence. We introduce the following W p,ℓ ∞ metric:
where | • | ∞ is the ℓ ∞ norm of a vector and Π(µ, ν) represents the set of measures in the joint space R d × R d that have marginals µ, ν.
We note that |x -y| ∞ ≥ |x (j)y (j) | for any 1 ≤ j ≤ d, where we use the superscript (j) to denote the j-th component of a vector. This observation implies that W p,ℓ ∞ (µ, ν) serves as an upper bound for the W p distance between one-dimensional marginals of µ and ν. Moreover, since K|x -y| p ∞ ≥ K t=1 |x (jt)y (jt) | p for any 1 ≤ j t ≤ d, we have that K 1/p • W p,ℓ ∞ (µ, ν) serves as the upper bound for the W p distance between any K-dimensional marginals of µ and ν. In summary, the metric W p,ℓ ∞ is capable of measuring the accuracy of low-dimensional marginals. In this paper, we specifically focus on the case p = 2, i.e., the W 2,ℓ ∞ metric.
Another important observation is that
which, combined with the dimension-independent contraction result in the W 2 metric, implies that the W 2,ℓ ∞ bias governs the dependence of the W 2,ℓ ∞ error on d.
Positive and negative examples: delocalization effect. How should we expect the W 2,ℓ ∞ bias of (1.2) behave? To motivate the discussion, we first consider the examples of product measures and Gaussian measures. We begin with the product measure case for which a bound on W 2,ℓ ∞ (π h , π) can be obtained by a contraction argument similar to the W 2 analysis (e.g., [12]). We include a sketch of the argument for later reference in the sketch of the proof of our main result, Theorem 2.2 in Section 2.2.
Sketch of Proof. Given k ∈ N, we couple the continuous-time Langevin dynamics at stationarity Y t , kh ≤ t ≤ (k + 1)h and the discrete-time iterates X kh in (1.2) so that they share the same Brownian motion. Introducing an auxiliary random variable
) and using the triangle inequality then lead to
.
Here (b) is the one-step discretization error which is bounded by
where we apply the strong convexity of V i to bound each coordinate of the vector yielding the final ℓ ∞ norm bound. This shows that there is a one-step contraction in the ℓ ∞ norm. We then couple the distribution of X kh and
Iterating this inequality implies W 2,ℓ ∞ (π h , π) = O( β α h log(2d)). A detailed proof can be found in Appendix A.
Example 1.2 shows that the W 2,ℓ ∞ bias scales only as the square root of log(2d); that is, the bias is nearly independent of the dimension. The key in the proof for Example 1.2 is the one-step contraction property in the ℓ ∞ norm, which relies on the structure of the product measures. This property does not hold for general π. On the other hand, we know that for Gaussian distributions, we have an explicit formula for the law of iterates in the algorithm and thus the biased distribution π h , so we can investigate the W 2,ℓ ∞ bias directly without concerning the one-step iteration property; see the following Example 1.3 with proof in Appendix B.
Again, Example 1.3 shows that the W 2,ℓ ∞ bias is nearly independent of the dimension. Using the property of the W 2,ℓ ∞ metric, this further implies that the W 2 distance between K-marginals of π h and π is bounded by O( Kh log(2d)), nearly independent of d. We can interpret this as the overall bias being nearly delocalized accross all onedimensional marginals. A step size of O(1/K), up to logarithmic terms, suffices for a bounded bias and error in K-marginals for any dimension; the iteration complexity scales with K and is also nearly independent of d.
On the other hand, a simple example shows that delocalization of bias does not always hold:
Example 1.4. Consider π = ρ ⊗d where ρ is a one-dimensional centered log concave distribution for which the biased distribution ρ h obtained by the corresponding one dimensional unadjusted Langevin has nonzero mean 1 , so that their mean differs by δ > 0.
Now, consider the rotation matrix Q which satisfies (Qx
1 For example, Proposition 3.1 can be used to show that ρ can be taken to be a univariate mixture of Gaussians pN (µ1, 1) + (1p)N (µ2, 1) with pµ1 + (1p)µ2 = 0, 0 < p < 1/2, and |µ1 -µ2| sufficiently small.
We then find
Example 1.4 shows that for the rotated product measure π, the W 2,ℓ ∞ bias is of order √ d. This indicates that the bias is not declocalized, but concentrated on one specific dimension. Thus, the delocalization effect of the bias over marginals is a delicate phenomenon and does not hold universally. We note that Example 1.4 is characterized by strong, dense interactions between coordinates after the rotation. 1.3.3. Strongly log-concave target distributions with sparse interactions. In this paper, our main result is to show that the delocalization effect holds for distributions with sparse interactions; see the following Figure 1 and Theorem 1.5. Some i-th variable:
Theorem 1.5 (informal). Let Assumption 1.1 hold. Suppose V further satisfies the sparsity conditions illustrated in Figure 1 and rigorously formulated in Section 2.1 below with the sparsity parameter
The precise statement of the theorem is in Section 2. The key technical challenge in extending the delocalization of bias result beyond product measures and Gaussian measures is the lack of one-step contraction in the W 2,ℓ ∞ metric and the lack of explicit formulas of the solutions. We employ a novel multi-step coupling argument to derive a multi-step contraction result and crucially use the sparsity of the potential to control the accumulated errors over these steps to prove the theorem.
The theorem implies that for such a sparse potential, the W 2,ℓ ∞ bias is nearly independent of d, up to logarithmic terms and assuming the condition number β/α is also nearly independent of d. As a consequence, by taking h ∼ 1/K, the iteration complexity scales with K to ensure a desired W 2 error for all K-marginals, up to logarithmic terms on d. Thus, for these potentials, the unadjusted Langevin algorithm proves more scalable than its unbiased variants when the quantities of interest are low-dimensional marginals. We note that these sparse potential illustrated in Figure 1 can arise in physical applications and Bayesian inference; see Section 2.3.
Generalizing the delocalization of bias effect beyond sparse potentials can provide further insights on the behavior of the algorithm in high dimensions. We also approach this through an alternative asymptotic perspective, in the spirit of methods using Taylor's expansion and the Poisson equation [39,40,1,31]. We derive explicit formulas for the bias of observables in first orders and demonstrate how the delocalization of bias, in the context of observables, can hold in a wider generality.
1.4. Organization of this paper. Section 2 presents our main result on sparse potentials. In Section 3, we discuss generalizations of our result through asymptotic arguments. We conclude the paper in Section 4. All proofs are deferred to the appendices. 1.5. Notations. We write A = O(B) or A B to mean that there exists a constant C independent of α, β, h, d such that A ≤ CB. On the other hand, A = Ω(B) or A B means that there exists a constant
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.
In this section, we begin by mathematically defining the potential with sparse interactions (as illustrated in Figure 1) in subsection 2.1. Subsection 2.2 presents the convergence bound. Finally, in subsection 2.3, we discuss examples of sparse potentials that satisfy the assumptions of our theorem.
2.1. Sparse graphical models. Consider an undirected graph G with d nodes, labeled by 1 ≤ i ≤ d. If there is an edge between node i and j, we write i ∼ j, meaning that the two nodes are connected. We denote by N(i) the neighborhood of the node i, which is a set of nodes connected to i. Here the neighborhood relationship is symmetric, and without loss of generality, we assume i ∼ i for all 1 ≤ i ≤ d. Furthermore, we define recursively that
where X i = {x (j) : j ∈ N(i)} and V i is a function depending on variables in the set X i .
We can also understand G as a factor graph, where F = {X i , 1 ≤ i ≤ d} represents the factor vertices. By definition, the potential V described above exhibits sparse interaction when the graph is sparse.
Convergence bound for unadjusted Langevin. Below we present the bias and convergence bounds in terms of the W 2,ℓ ∞ metric. The sparsity parameters s k will play important roles in our bound. Theorem 2.2. Let Assumptions 1.1 and 2.1 hold. Assume 0 ≤ h ≤ 1/β and denote q = exp(-hα). Let r i = ⌈e 2 ihβ + log √ d⌉ for i ∈ N. Then, for any N ∈ N such that 2q N √ s r N < 1, we have the bound
Moreover, the following convergence bound holds for any k ∈ N:
Sketch of proof. Using the same coupling argument and notations in Example 1.2, we have
,
. Now π is not a product measure and we no longer have the contraction (a)
. In fact, using Taylor's expansion, we can write
where
When π is a product measure, H k is a diagonal matrix and thus its ℓ ∞ norm equals ℓ 2 norm so we can get the contraction property. However, in general, |H k | ∞ can be much larger than
on the fact that H k has only s 1 nonzero entries in each column. As a result, we must employ a multi-step coupling argument to analyze the iterations. We now have
where error(1) = O(βh 3/2 log(2d)). Applying the bound again within the first term on the right hand size, yields
where error(2) represents the two-step error. More generally, for any k, N ∈ N, we get
Now by choosing a large N , we will get
Thus N ∼ log d h will lead to a contraction. In our proof, we leverage the sparsity of the potential to produce a potentially tighter bound
The primary challenge lies in controlling the growth of error(N ). A naïve approach would lead to polynomial dependence on d since the bound may grow exponentially fast in N without contraction. In contrast, our analysis hinges on a sparsity bound of the propagator of the unadjusted Langevin algorithm over multiple steps; see Appendix C.1. This sparsity enables us to obtain tighter bounds on these ℓ ∞ errors across multiple iterations. In fact, we can see that in (2.2), r i scales with the physical time ih, which justifies that s r i can characterize the sparsity of interactions after i steps of the algorithm. Using the sparsity bound, we finally control the accumulated discretization errors in N steps by
With the N -step contraction and the bound on the accumulated discretization errors, we will get
we can let k → ∞ and get the final bound in (2.2). The complete proof for Theorem 2.2 can be found in Appendix C. Overall, Theorem 2.2 implies that the bias depends on the sparsity growth, the expectation of |∇V | 2 ∞ , and log(2d). In the following, we state a result that provides an upper bound on the expected |∇V | 2 ∞ . Proposition 2.3. Let Assumptions 1.1 hold. Then, it holds that (2.4)
Sketch of Proof. When V is strongly convex, by the Bakry-Émery criterion, π satisfies the log-Sobolev inequality with constant 1/α. Using this fact and Herbst's argument leads to the bound (2.4). The complete proof can be found in Appendix C.3. Now, we are ready to provide one instantiation of Theorem 2.2 and Proposition 2.3 for the case that the sparsity of the neighborhood of the graphical model grows at most polynomially. We need to carefully calculate an explicit bound on the terms in (2.2); the detailed proof is in Appendix C. Section 2.3 will provide concrete examples illustrating this case.
Theorem 2.4. Let Assumptions 1.1 and 2.1 hold with s k ≤ C(k + 1) n for any k ∈ N, where C > 0 is a generic constant and n ∈ N. Then for h ≤ α β 2 , we have
, O
Moreover, it holds that (2.6)
In Theorem 2.4, we express
in this manner to ensure that the hidden constant is also independent of n. In the worst-case scenario, where there is no sparsity present at all, n scales with log d. Theorem 2.4 implies that the W 2,ℓ ∞ bias is nearly independent of d, up to logarithmic terms, if we assume n is independent of d and the condition number β/α is also nearly independent of d. In such case, by taking h ∼ 1/K, the iteration complexity scales with K to ensure a bounded W 2 error for all K-marginals, up to logarithmic terms on d, per the discussion in Section 1.3.
Examples of sparse graphical models. In this subsection, we provide examples that satisfy Assumptions 1.1 and 2.1. For the bound in Theorem 2.4 to be non-trivial, the ratio β/α should not grow too rapidly with the dimension. The simplest case would be when α and β are fixed constants, independent of the dimension. It is not immediately apparent that given constant values of α and β such that α ≤ β, there exists a potential function V for any dimension d ∈ N with αI ≼ ∇ 2 V (x) ≼ βI at any point x ∈ R d , and where s k grows at most polynomially in k, as required in Theorem 2.4. As the following example shows, a natural class of examples of such potentials arises from Laplacians of sparse graphs.
Example 2.6. We provide an example of V : R d → R for which αI ≼ ∇ 2 V (x) ≼ βI at any x ∈ R d with α, β independent of d, and s k = min{2k + 1, d}.
Consider the following matrix
This matrix can be seen as -∆ + λ(x)I, where ∆ is the discretization of the 1D Laplace operator with homogeneous Dirichlet boundary conditions. This operator satisfies 0 More generally, consider a sparse, bounded-degree graph G on the variables x (1) , . . . , x (d) and potentials of the form
, where V 1 , . . . , V d are arbitrary smooth and strongly convex functions and E(G) is the edge set of G. The resulting probability measures are log-concave perturbations of a Gaussian free field. Then V will satisfy Assumptions 1.1 and 2.1 with α and β independent of the dimension, analogous to Example 2.6. An important special case is when G is a subset of the lattice Z n , in which case the sparsity parameter satisfies
These examples can arise in physical systems where there are local interactions and a restoring force for each state which plays the role of V i (x (i) ), or in Bayesian inverse problems where such V i (x (i) ) terms can come from priors.
In this section, we discuss generalizations of the results beyond log-concave distributions with sparse interactions. While our analysis requires strict sparsity, where most interactions are zero, we anticipate that the analysis can be extended to cases where all but a small number of interactions are weak. Moreover, our sparsity assumptions treat different coordinates equally, and it is of interest to consider heterogeneous sparse models where some coordinates have more, but potentially weaker, connections with others. The study of these scenarios is left for future endeavors.
Alternatively, we can approach the question from a different perspective, employing asymptotic analysis. More precisely, we adopt the Poisson argument and utilize the generator of the Markov process to derive an asymptotic formula for the bias of certain observables. By doing so, we can gain additional insights into the algorithm's bias.
Recall our definition V (x) =log π(x). Let L and L h be the generators of the continuous Langevin dynamics and the unadjusted Langevin algorithm, respectively. By definition, Lu = ∇ log π • ∇u + ∆u and
Below we calculate the first order formula for the bias of an observable f . We assume that f is sufficiently regular so that the Taylor expansions and integration by parts used in the Proof of Proposition 3.1 below are valid. In this section, we use the notation A = o(h) to mean that the term A is a higher order term than h as h → 0.
Proposition 3.1. Without loss of generality, assume f π = 0. Then it holds that
Moreover, we also have the following formula:
The proof of the formula can be found in Appendix D. Let us make some observations regarding Proposition 3.1. If π is Gaussian, then f (∆ log π)π = 0 since f π = 0. We have
This means that if f is a linear observable, then the first-order term of the bias is zero. Moreover, if f depends only on a small number of coordinates of x ∈ R d , then the integral (∆f )π will scale with that number, rather than d, because the entire integral will rely only on the marginal distribution of π at these coordinates. Thus, to leading order, the bias of the observable is delocalized across dimensions.
In fact, the above argument can be further generalized to π that is a perturbation of Gaussians; see the following Proposition 3.2, and its proof is in Appendix D. Proposition 3.2. Let π ∝ exp(-V (x)) ∝ N (x; m, Σ) exp(-U (x)) where N (x; m, Σ) is the density of a Gaussian with mean m and covariance Σ. Assume π is centered, and (∆U ) 2 π ≤ C 0 for some constant C 0 independent of d. Suppose f depends only on K coordinates of x ∈ R d , then the first-order term of the bias depends on K, not d. In particular, for the observable f (x) = K k=1 x (k) , under the additional assumption that π satisfies a Poincaré inequality with a positive dimension-independent constant, we have
where C 1 is a constant independent of d.
Proposition 3.2 shows that the delocalization of bias for observables holds if π is a perturbation of a Gaussian distribution, suggesting that the delocalization effect may hold in significant generality. Note, however, the asymptotic arguments in this section do not provide insight into the behavior of the higher order contributions to the bias. On the other hand, the asymptotic formula is applicable beyond log-concave distributions and can be applied to study the bias of any observable of interest. Proposition 3.2 also suggests that the Gaussian part may not have a significant effect on the bias of the observable f . Recall that our W 2,ℓ ∞ bounds apply to Gaussian distributions (which can have dense interactions) and log-concave distributions with sparse interactions. It is natural to inquire whether our W 2,ℓ ∞ analysis can be extended to the case of a distribution that is the product of a dense Gaussian and a sparse log-concave distribution.
In this article, we studied the convergence of unadjusted Langevin algorithms in high dimensions. For strongly log-concave distributions, existing results showed that the iteration complexity scales proportionally to d or √ d in order to achieve a desired error in the W 2 metric for any dimension d. We demonstrate that for Gaussian distributions or distributions with certain sparsity structures, a constant number of iterations, up to some logarithmic terms in d, suffices to achieve a bounded error in the W 2,ℓ ∞ metric. Consequently, a number of iterations proportional to K (up to some log d terms) can achieve a bounded W 2 error for all K-marginals. This result implies that even in extremely high-dimensional settings, unadjusted Langevin algorithms can still be highly scalable if the quantities of interest depend only on low-dimensional marginals. We note that this desirable property is not satisfied for unbiased schemes such as the MALA or proximal samplers, where the constraint on the step size necessitates poor scaling in d.
Our discussion hinges on a delocalization effect of the bias: the overall bias of the algorithm is delocalized to a comparable order across each one-dimensional marginal distribution. Thus, to obtain accurate low-dimensional marginals, one does not suffer from a poor dependence on the full dimension d. The delocalization effect is rigorously shown in this paper for Gaussian distributions and distributions with certain sparsity structures. Our error bounds have a favorable dependence on d; however, the dependence on the condition number β/α may potentially be improved. We also provide counterexamples of a rotated product measure for which the bias is not delocalized.
Extending the delocalization analysis to other broader classes of distributions is of future interest. In fact, our asymptotic arguments based on Poisson equations provided further insights into how the first-order terms of the bias of general observables scale with dimension.
More broadly, it is an important question to understand how the complexity of other unadjusted MCMC algorithms scales with dimension when the quantities of interest depend solely on low-dimensional marginals. Addressing this question can provide valuable insights to help practitioners select appropriate algorithms and understand their computational complexity when dealing with high-dimensional sampling problems.
Proof. First, for the continuous-time Langevin dynamics Y t and the discrete-time iterates X kh in the unadjusted Langevin algorithm, we have:
We couple the two processes using the same Brownian motion. We aim to estimate
. For this purpose we introduce an auxiliary random variable
Using the triangle inequality then leads to
.
For (a), we have
where we used the facts that
where in the last inequality, we used the fact that ∇ 2 V (uY t + (1u)Y kh ) is a diagonal matrix with each diagonal entry bounded by β in magnitude. We further have (A.7)
where we used the fact that all Y t ∼ π, and the bound E|B u | 2 ∞ ≤ 4u log(2d) holds for any u ≥ 0 due to Lemma B.1.
Then by Proposition 2.3, we have
. Combining all the bounds above, we get
We can now couple the distribution of X kh and Y kh such that
With this and using the definition of the W 2,ℓ ∞ norm, we get
Iterating this inequality leads to the bound on the W 2,∞ bias. Furthermore under the assumption h ≤ α β 2 , we get
The proof is complete.
We first state a lemma for the expected squared maximal norm of a random vector whose entries are sub-Gaussians. This lemma will also be used in the proof of our main theorem.
Lemma B.1. Suppose Y = (Y (1) , Y (2) , ..., Y (d) ) ∈ R d and each Y (i) is centered and sub-Gaussian with variance proxy σ 2 , namely
Then, it holds that
By the property of sub-Gaussian random variables [42], we have that for 0 ≤ λ < 1/(2σ 2 ),
Using the convexity of the exponential function, we can derive
Then, we can bound the right hand side as follows:
Thus, combining the above two inequalities and taking logarithms, we find
With Lemma B.1, we present the proof for the statement in Example 1.3 as follows.
Proof for Example 1.3. When π is Gaussian, the potential V = 1 2 (xm) T Σ -1 (xm) is quadratic. In this case, the iteration takes the form
For Σ k , we have the identity:
Letting k → ∞, we get
∞ Z + m where Z ∼ N (0, I). Then, for this specific coupling, we have the bound
Then we have the formula
where
C.1. Sparsity of the propagator of unadjusted Langevin. First, we present a proposition concerning the ℓ ∞ norm of matrices that include the propagator of the Langevin Monte Carlo algorithm as a special case. We will use this proposition when analyzing the convergence of the algorithm in the next subsection.
Proposition C.1. Assume 0 ≤ h ≤ 1/β and denote q = exp(-hα). For the potential function in Assumption 2.1, the following facts hold:
, which is the multiplication of r Hessian matrices, is s r sparse, for any x 1 , x 2 , ..., x r ∈ R d , r ∈ N. Here we say a matrix is s r -sparse if each row of the matrix contains at most s r nonzero entries. (ii) Let ν 1 , ..., ν N be any probability measures in R d . Define the matrix (C.1)
We have the inequality
for any r ≥ e 2 Nhβ. In particular, taking r N = ⌈e 2 Nhβ + log √ d⌉, we get that
(iii) In the context of (ii), consider additionally a probability measure ν 0 and the matrix (C.4)
We have the inequality
for any r ≥ e 2 Nhβ. In particular, taking r N = ⌈e 2 Nhβ + log √ d⌉, we get that
Proof of Proposition C.1. We prove the above facts one by one.
Proof for (i). We know that the ij-th entry of the matrix
is nonzero if and only if j ∈ N r (i). Since max 1≤i≤d |N r (i)| = s r , we will have at most s r nonzero entries in each row. This implies that the matrix is s r -sparse.
Expanding the product, we have (C.7)
where in the last inequality, we used the fact that k! ≥ k k / exp(k). This fact can be seen by taking x = k in the inequality exp(x) ≥ x k /k!.
When r ≥ e 2 Nhβ, we can bound (C.9)
Furthermore, we obtain r) by the Cauchy-Schwarz inequality.
where we used the fact that h ≤ 1/β. Then, for each row i of P N , (C.10)
where in the second inequality, we used the fact that j / ∈Nr(i)
since a nonzero term (P N ) ij for j / ∈ N r (i) can only be produced by product of more than r matrices of the kind of A i k , as a consequence of the argument in (i). In summary, we get the result that if r ≥ e 2 Nhβ, then
We note that the above inequality will also hold true when r ≥ d.
Taking r N = ⌈e 2 Nhβ + log √ d⌉, we get (C.11)
Proof for (iii). Note that J N = ∇ 2 V dν 0 H N . Following the proof for (ii), we have that when r ≥ e 2 Nhβ, (C.12)
Thus the ∞-norm of the above matrix is bounded by √ dβ exp(-r). Then, for any
The proof is complete. C.2. A multistep coupling argument. To prove Theorem 2.2, we first show the following Proposition C.2, which bounds the W 2,ℓ ∞ metric through a multistep coupling argument. Note that the one-step coupling argument that is commonly used in bounding W 2 distance is not enough here, as the one-step contraction property is lost in the W 2,ℓ ∞ metric. Proposition C.2. Let Assumptions 1.1 and 2.1 hold. Assume 0 ≤ h ≤ 1/β and denote q = exp(-hα). Consider the Langevin dynamics dY t = -∇V (Y t )dt + √ 2dB t with Y 0 ∼ π. Let ρ kh represent the law of X kh from the iterations of Langevin Monte Carlo with stepsize h. Then, the following estimates hold for any k, N ∈ N:
where r i , 1 ≤ i ≤ N is defined in Proposition C.1, and
Proof of Proposition C.2. Let dY t = -∇V (Y t )dt+ √ 2dB t . We write down the following identity:
We couple the two processes using the same Brownian motion. The goal is to estimate
To do so we introduce an auxiliary random variable (C. 16)
Using the triangle inequality then leads to (C.17
.
For (a), we have
where
We can view H k as a random matrix depending on X kh and Y kh . In particular, each realization of this random matrix can be written in the form Ih ∇ 2 V dν where ν is a probability measure; here ν has its mass on a line. This form allows us to use Proposition C.1 to analyze the iterations.
For (b), we have (C. 19)
which can be understood as the one-step discretization error of the overdamped Langevin dynamics. In the above, we used the fact that
Now, combining the above estimates, we arrive at
We move one step back further and get
,
where, similar as before, we define
). For (c), we use the same argument as earlier to get
where we applied (iii) of Proposition C.1 to the the matrix
As a summary, we get
Iterating the above arguments N times, and for simplicity of notations writing k + N in place of k + 1, we get
where in the last step, we applied (ii) of Proposition C.1 to the matrix
We can now couple the distribution of X kh and
With this and using the definition of the W 2,ℓ ∞ norm, we get
With this expression, we can ensure contraction in the W 2,ℓ ∞ metric by iterating sufficiently large number of steps N , so that the factor q N is small enough to offset the increasing sparsity parameter s r N . In particular, in the dense case where s r N = d, such a contraction can be achieved by taking N ∼ log d h . We can further bound ϵ 2 j using the same approach as in (A.7), which implies that (C.24)
The proof is complete.
Proof of Theorem 2.2. With the bound on ϵ j in Proposition C.2, we take the limit k → ∞ in (C.14) to obtain the bias bound.
Once we have the bound on the bias, we can utilize the convergence bound for W 2 to establish the following:
where we used the contraction in W 2 to get W 2 (ρ kh , π h ) ≤ q k W 2 (ρ 0 , π h ). The proof is complete.
Proof of Proposition 2.3. By the Bakry-Émery criterion, for distribution π that is αstrongly log concave, π satisfies the log-Sobolev inequality with constant 1/α. Then, by Herbst's argument [42], we have
) .
Furthermore, since ∇V is β-Lipschitz, Herbst's argument leads to
where ∇ i V is the i-th component of the vector ∇V . This implies that the distribution of ∇ i V (Y ) has a sub-Gaussian tail. Thus, by Lemma B.1, we get (C.28)
Note that E π [∇V (Y )] = -π∇ log π = ∇π = 0. The proof is complete.
Proof of Proposition 2.4. We will use Theorem 2.2. Recall the definition r N = ⌈e 2 Nhβ+ log √ d⌉ and q = exp(-hα). There is a universal upper bound on s i , given by s i ≤ d.
Choose N = ⌈ log(4 √ d) hα ⌉, which leads to 2q N √ s r N ≤ 1/2. Based on Theorem 2.2, it remains to calculate the bound on N i=1 q i-1 √ s r i . Since q ≤ 1, we have the bound (C.29)
√ C(e 2 ihβ + log √ d + 2) n 2 ≤ √ C N +1 1 (e 2 yhβ + log √ d + 2) n 2 dy ≤ √ C Ä e 2 (N + 1)hβ + log √ d + 2 ä n 2 +1 (n/2 + 1)(e 2 hβ) ≤ √ C Ä 4e 2 log(4 √ d) β α + log √ d + 2 ä n 2 +1
(n/2 + 1)(e 2 hβ) .
Therefore, we get (C.30)
(n/2 + 1)(e 2 h)
.
On the other hand, the trivial bound s i ≤ d leads to (C.31)
Thus, we have another bound
Moreover, by Proposition 2.3, we have
Using the fact that h ≤ α β 2 and combining the above two inequalities leads to the final result. i ] = 3 for ξ i ∼ N (0, 1). In the third identity, we performed the Taylor expansion at x, based on the fact x = x + h∇ log π(x). In the last identity, we noticed the fact that the terms in the big bracket equals ∆ 2 u.
Moreover, u(x)u(x) = h∇u(x) • ∇ log π(x) + 1 2 h 2 (∇ log π(x)) T ∇ 2 u(x)∇ log π(x) + o(h 2 ). Therefore, we get When f only depends on K number of coordinates of x ∈ R d , then the integral (f ) 2 π will only scale with that number, rather than d, because the whole integral will only rely on the marginal distribution of π at these coordinates. As such argument also applies