Source-linked AI summary

High-dimensional Bayesian inference via the Unadjusted Langevin Algorithm

Alain Durmus, Eric Moulines

arXiv:1605.01559v4math.STstat.MEstat.ML

TL;DR

The paper studies how to sample high-dimensional distributions, including Bayesian posteriors, using ULA, the Euler discretization of the Langevin stochastic differential equation. Under smoothness and strong-convexity assumptions, it gives iteration bounds for Wasserstein and total-variation accuracy, with explicit dependence on dimension.

  • Problem

    Sampling high-dimensional Bayesian posteriors matters because the full posterior represents parameter uncertainty and can help prevent overfitting.

    Source: §1 · pp. 1–2
  • Method

    The paper analyzes ULA, the Euler discretization of the Langevin stochastic differential equation, with constant and decreasing step sizes.

    Source: §1 · p. 2
  • Results

    For fixed precision, ULA’s iteration bounds for Wasserstein or total-variation accuracy are O(dε^-2) or O(dε^-1), depending on smoothness, up to logarithmic terms.

    Source: §1 · p. 3
  • Takeaways & Limitations

    For d-dimensional Gaussian targets, the fixed-precision iteration bound is optimal up to logarithmic factors.

    Source: §1 · p. 3

    The analysis is scoped to potentials U that are continuously differentiable and have globally Lipschitz gradients.

    Source: §2 · p. 5

Abstract

from arXiv · show

We consider in this paper the problem of sampling a high-dimensional probability distribution $π$ having a density with respect to the Lebesgue measure on $\mathbb{R}^d$, known up to a normalization constant $x \mapsto π(x)= \mathrm{e}^{-U(x)}/\int_{\mathbb{R}^d} \mathrm{e}^{-U(y)} \mathrm{d} y$. Such problem naturally occurs for example in Bayesian inference and machine learning. Under the assumption that $U$ is continuously differentiable, $\nabla U$ is globally Lipschitz and $U$ is strongly convex, we obtain non-asymptotic bounds for the convergence to stationarity in Wasserstein distance of order $2$ and total variation distance of the sampling method based on the Euler discretization of the Langevin stochastic differential equation, for both constant and decreasing step sizes. The dependence on the dimension of the state space of these bounds is explicit. The convergence of an appropriately weighted empirical measure is also investigated and bounds for the mean square error and exponential deviation inequality are reported for functions which are measurable and bounded. An illustration to Bayesian inference for binary regression is presented to support our claims.

1 Introduction

Motivated by Bayesian methods that sample full posteriors to represent uncertainty and limit overfitting [33], [42], this work gives explicit, dimension-dependent non-asymptotic guarantees for Euler-discretized Langevin sampling.

1 Introduction

  • For strongly convex U, explicit bounds for ULA in Wasserstein and total variation cover fixed and non-increasing steps, improving earlier total-variation bounds and.The bounds concern the distribution of the n-th iterate relative to the target π.Source: §1 · p. 3
  • At fixed precision ε, ULA needs O(dε^-2) Wasserstein or O(dε^-1) total-variation iterations, up to logarithmic terms, depending on U's smoothness.For a fixed horizon n, the Wasserstein error can be bounded by O(n^-1/2) or O(n^-1), depending on smoothness; the precision result is optimal up to logarithmic factors for d-dimensional Gaussians.Source: §1 · p. 3
  • For decreasing steps γk = γ1k^-α, the marginal distribution converges to π, with the optimal reported Wasserstein and total-variation rates at α = 1 and sufficiently large γ1.The reported rates improve those in.Source: §1 · p. 3
  • The stationary-distribution bias between π and πγ is bounded in Wasserstein and total variation, with total variation scaling as γ for sufficiently smooth U.The total-variation bound improves the one reported in.Source: §1 · p. 3
  • Weighted empirical measures are analyzed using mean-square-error and exponential-deviation bounds for Lipschitz or bounded measurable functions.For Lipschitz functions under sufficient smoothness, the optimal any-time rate with decreasing γk = γ1/k^α occurs at α = 1/3. With constant steps, the MSE target ε requires O(dε^-4) or O(dε^-3) iterations depending on U's smoothness; for bounded measurable functions, the requirement is the same up to logarithmic terms.Source: §1 · pp. 3–4

2 Non-asymptotic bounds in Wasserstein distance of order 2 for ULA

Under global Lipschitz-gradient and strong-convexity assumptions, the section establishes Wasserstein-2 convergence bounds for Langevin diffusion and its Euler discretization. It also examines how step-size schedules and additional smoothness affect convergence rates and dimension dependence.

2 Non-asymptotic bounds in Wasserstein distance of order 2 for ULA

  • The diffusion contracts exponentially in W2 at rate m, while the Euler kernels are strict contractions and converge geometrically, with both stated rates independent of dimension.These results assume H1 and H2; the Euler contraction applies under the stated step-size conditions.Source: §2 · pp. 6–7
  • For constant steps γ ≤ 2/(m + L), the Euler kernel has a unique stationary distribution πγ; decreasing steps with γ1 ≤ 1/(m + L), γk → 0, and cumulative step size diverging yield convergence to π.The convergence-to-π condition is stated under H1 and H2.Source: §2 · pp. 7–8
  • When ˜L = 0, including non-degenerate Gaussian targets, the iteration-complexity dependence on dimension is O(d^1/2 log(d)); under H3, the precision dependence improves.Under H1 and H2, the reported complexity matches results for total variation distance.Source: §2 · p. 9
  • For finite-horizon sampling, a step size can be chosen for each n to give W2^2 error at most O((log(n)/n)^2) when H3 holds.The passage states this bound alongside the H1 and H2 finite-horizon result.Source: §2 · p. 9
  • With γk = γ1/k and γ1 > 2κ^-1, W2^2 error is O(n^-1) under H1 and H2; under H3, the same rate is reported for γ1 > 4κ^-1.The constants are exponential in γ1, and the text notes that the rates are explicit.Source: §2 · p. 10

3 Quantitative bounds in total variation distance

The section derives total-variation bounds for continuous Langevin dynamics and its Euler discretization, motivated in part by computing highest posterior density credible regions and intervals in Bayesian inference.

3 Quantitative bounds in total variation distance

  • Reflection coupling yields an upper bound on total variation between semigroup evolutions under H1 and H2.The section uses the semigroup’s regularizing property to obtain total-variation bounds.Source: §3 · pp. 10–11
  • The discrete-time bounds parallel the continuous results, using a different coupling because reflection coupling does not apply in discrete time.Theorem 11 assumes a non-increasing step sequence with γ1 ≤ 2/(m + L); its coupling construction is inspired by a method for Gaussian random walks.Source: §3 · p. 12
  • The passages report a total-variation error of d^1/2 O(ℓ^-3α/4), with γ1 > 4κ^-1 required when α = 1.Source: §3 · p. 14
  • For constant-step sampling under H1–H2, the iteration bound to reach TV precision ε is d log(d) O(|log(ε)| ε^-2).The text states this is the same order as the result in.Source: §3 · p. 15
  • For constant steps, the stationary TV error is at most C1 γ^1/2, while H3 improves its order to d O(γ |log(γ)|).C1 is of order d^1/2; when ˜L = 0, the H3 bound is d^1/2 O(γ |log(γ)|) and is sharp up to a logarithmic factor.Source: §3 · pp. 15–16
  • Under H1–H3, a well-chosen constant step needs d log^2(d) O(ε^-1 log^2(ε)) iterations, improving the H1–H2 bound.For the standard Gaussian case, the reported bound is sharp up to logarithmic factors in dimension and precision.Source: §3 · pp. 15–16

4 Mean square error and concentration for bounded measurable functions

The section bounds the mean-square error of a weighted Langevin estimator for Lipschitz and bounded measurable functions, and establishes exponential deviation bounds under non-increasing step sizes. It also identifies optimal polynomial step-size decay and fixed-horizon MSE rates under the paper’s assumptions.

4 Mean square error and concentration for bounded measurable functions

  • The weighted-average estimator’s MSE decomposes into squared bias and variance, and the section derives bias bounds for Lipschitz or bounded functions.The estimator approximates π(f) using the Euler-discretized Langevin chain with non-increasing step sizes.Source: §4 · pp. 17–18
  • For polynomial step sizes gamma_k=gamma1/k^alpha, the infinite-horizon MSE bound is optimized at alpha=1/2 under H1 and H2, and alpha=1/3 under H1, H2, and H3.These conclusions combine the squared-bias and variance bounds; their orders are summarized in Table 6.Source: §4 · p. 18
  • For Lipschitz functions, the variance bound is independent of dimension under H1 and H2 for non-increasing step sizes with gamma1≤2/(m+L).The variance result applies for all N≥0 and n≥1.Source: §4 · p. 18
  • With constant step size γ and fixed total iteration count n + N, optimizing γ and burn-in period N yields an upper bound of order n^-1/2 under H1 and H2, and n^-2/3 under H1, H2, and H3.Source: §4 · p. 18
  • Theorem 16 gives an MSE bound for bounded measurable functions under H1 and H2, with non-increasing step sizes satisfying gamma1≤2/(m+L).It applies for all N≥0, n≥1, and starting points x in Rd.Source: §4 · p. 19

5 Numerical experiments

The numerical experiments apply ULA to Bayesian binary logistic regression, comparing it with Polya–Gamma Gibbs sampling and MALA; the reported results support ULA as an alternative.

5 Numerical experiments

  • The experiments consider Bayesian binary logistic regression with d-dimensional covariates and coefficients and a zero-mean Gaussian prior with covariance Σβ.Source: §5 · p. 21
  • The Polya–Gamma data-augmentation algorithm [36] has a cubic-dimensional augmentation step and an ergodicity constant exponentially small in parameter dimension and sample count [7, Proposition 1].The cited result establishes uniform ergodicity in total variation and gives an explicit expression for the constant.Source: §5 · p. 21
  • ULA and the Polya–Gamma Gibbs sampler [36] are compared using first-component histograms with constant and decreasing step sizes.The comparison uses synthetic data with d = 5 and p = 100; the Gibbs sampler produces 10^8 samples, while ULA uses 10^3 runs of 10^6 effective iterations.Source: §5 · p. 21, Figure 2 · p. 22
  • The real-data comparison also uses preconditioned MALA and ULA, transforming draws back to the original posterior parameterization.The samplers target a density transformed by the covariance matrix before post-multiplication of draws.Source: §5 · p. 22
  • For five real data sets, 100 runs each of MALA and ULA use the same constant step size, with MALA tuned to about 0.5 stationarity acceptance; marginal accuracy is averaged across dimensions.Each run uses 10^6 iterations, and posterior marginals are estimated using 2 · 10^7 Polya–Gamma Gibbs iterations.Source: §5 · pp. 21–22
  • Results across five real data sets imply ULA is an alternative to the Polya–Gamma Gibbs sampler and MALA.Figure 3 shows marginal-accuracy boxplots for the data sets.Source: §5 · p. 22, Figure 3 · p. 23

6 Contraction in total variation for functional autoregressive models

For functional autoregressive models with Lipschitz update functions and Gaussian noise, the section bounds total-variation distance between marginals using a coupling construction. The bound does not depend on the state-space dimension.

6 Contraction in total variation for functional autoregressive models

  • The model uses measurable update functions h_k, positive noise scales σ_k, and independent standard Gaussian inputs; each h_k satisfies the AR1 Lipschitz condition.The resulting process is an inhomogeneous Markov chain.Source: §6 · p. 23
  • A paired Markov kernel couples the two marginal chains, so Lindvall’s inequality bounds total variation by the probability that the coupled states differ.For initial states x and y, the paired process at time k couples δxQk and δyQk.Source: §6 · p. 24
  • Under AR1, Theorem 19 gives an upper bound on ∥δxQn − δyQn∥TV for all initial states and n ≥ 1 that is independent of dimension d.The supplied passage states the dimension-free scope but omits the displayed bound’s formula.Source: §6 · pp. 24–25

A Proofs of Section 2

This section collects the proofs postponed from Section 2 and begins an argument under H1 using results from [34].

A Proofs of Section 2

  • Under H1, the section invokes [34, Theorem 2.1.12, Theorem 2.1.9] for a claim stated to hold for all x, y ∈ R^d.The passage does not include the claim that follows.Source: A Proofs of Section 2 · p. 30

A.1 Proof of Proposition 1

The proof uses a quadratic Lyapunov function to control the diffusion and a coupled SDE to compare trajectories in Wasserstein distance. These arguments yield a stationary second-moment bound and support the remaining proposition claim.

A.1 Proof of Proposition 1

  • The quadratic Lyapunov function V(x)=||x−x⋆||^2 satisfies the drift bound AV(x)≤2(−mV(x)+d), yielding a differential inequality handled by Grönwall.The proof uses strong convexity, expressed as ⟨∇U(x)−∇U(x⋆),x−x⋆⟩≥m||x−x⋆||^2.Source: §A.1 · p. 31
  • The invariant distribution satisfies π(||x−x⋆||^2)≤d/m, obtained by truncation, dominated convergence, and monotone convergence.The argument applies invariance of π and Jensen’s inequality to the bound from Proposition 1-(i).Source: §A.1 · p. 31
  • A coupled SDE has a unique strong solution, and H2 with Grönwall bounds the coupled trajectories; their joint law then provides a coupling for the W2 distance.The proof concludes by bounding W2(δxPt,δyPt) by the square root of the expected squared distance between the coupled trajectories.Source: §A.1 · p. 31
  • The proof of Proposition 1-(iv) follows directly from parts (ii) and (iii).Source: §A.1 · p. 31

A.2 Proof of Proposition 2

The proof treats ℓ<n as trivial and handles ℓ≥n using inequalities for a restricted step size and induction. It also establishes uniqueness of the stationary distribution for R_γ.

  • The ℓ<n case is trivial, so the proof focuses on ℓ≥n.Source: §A.2 · p. 32
  • For ℓ≥n, the proof uses ∇U(x⋆)=0 and γ≤2/(m+L) in deriving the Q_n,ℓ relation, then concludes by induction.The stated step-size range is γ∈(0, 2/(m+L)).Source: §A.2 · p. 32
  • Because every compact set is accessible and small for R_γ, [32, Theorem 15.0.1] yields a unique stationary distribution π_γ.The proof then follows the lines of Proposition 1-(ii) using (46).Source: §A.2 · p. 32

A.3 Proof of Proposition 3 … A.6 Proofs of Theorem 8

These proof sections establish uniqueness of the ULA invariant distribution under a step-size bound and develop coupling-based arguments for Theorems 5 and 8 and Corollary 6.

A.3 Proof of Proposition 3

  • For γ ≤ 2/(m + L), Proposition 3’s transition kernel is a strict W2 contraction and has a unique invariant distribution.The proof invokes this contraction property on P2(R^d).Source: §A.3 · p. 33

A.4 Proof of Theorem 5

  • Theorem 5’s proof develops a technical diffusion estimate using the martingale for squared displacement and the strong-convexity assumption H2.The argument defines the squared displacement from the starting point and uses H2 to bound its evolution.Source: §A.4 · p. 33
  • Theorem 5 is proved by coupling the diffusion with the Euler interpolation and bounding their Wasserstein distance, using invariance of π for the diffusion.The proof takes the diffusion’s initial law to be π and applies a coupling estimate to the Euler iterates.Source: §A.4 · pp. 34–35

A.5 Proof of Corollary 6

  • The Corollary 6 proof shows its two error terms vanish when step sizes are non-increasing, tend to zero, and cumulative step time diverges.It uses the limits γ_k → 0 and Γ_n → +∞ to establish the required vanishing.Source: §A.5 · pp. 36–37

A.6 Proofs of Theorem 8

  • Lemma 24 supplies an auxiliary estimate under H1, H2, and H3 for non-increasing step sizes with γ_1 ≤ 1/(m + L).The lemma is stated for a coupling initialized with a joint law in P2(R^d × R^d).Source: §A.6 · p. 37
  • Theorem 8’s proof follows the proof of Theorem 5, using Lemma 24 in place of Lemma 22.Theorem 5’s argument likewise uses an auxiliary coupling lemma.Source: §A.4 · p. 34, §A.6 · p. 38

B Proofs of Section 3 … C.2 Bounds on u(4)

The appendix proofs derive convergence and empirical-measure bounds using contraction, Pinsker, and Gaussian Poincaré and variance controls, and establish sequence-independent bounds on u(4).

B Proofs of Section 3 / B.1 Proof of Theorem 12

  • The postponed Section 3 proofs derive Theorem 12 by combining a triangle-inequality decomposition with a Pinsker bound and H1-based estimates.The total-variation step uses [12, Equation 15] with the paper’s stated convention for total variation.Source: B Proofs of Section 3 · p. 38, §B.1 · p. 39

B.2 Proof of (24)

  • For constant step size γ in (0, 1/(m + L)], the proof derives bound (24) from Theorem 12 using exponential contraction.It uses κ ≥ 2m and 1 − t ≤ e^−t for t ≥ 0.Source: §B.2 · pp. 39–40

B.3 Proof of Theorem 14

  • Theorem 14 is proved under H1, H2, and H3 by bounding terms in a preliminary lemma and taking a limit using Theorems 10 and 11.The proof uses contraction bounds and Pinsker’s inequality in the preliminary lemma.Source: §B.3 · pp. 40–42

C Proof of Section 4 / C.1 Proof of Theorem 15 and Theorem 16

  • The postponed Section 4 proofs of Theorems 15 and 16 use Gaussian Poincaré inequality for the Gaussian transition and express empirical sums as martingale increments.The transition has mean y − γ∇U(y) and covariance 2γ Id; the cited inequality is [3, Theorem 3.20].Source: C Proof of Section 4 · p. 42, §C.1 · pp. 42–43

C.1.1 Proof of Theorem 15

  • Under H1 and H2 with non-increasing step sizes, Theorem 15’s proof controls martingale increments through Lipschitz and variance bounds.The step sizes satisfy γ1 ≤ 2/(m + L); the proof invokes Lemma 26 and Corollary 28, with Lemma 27 supplying a variance bound.Source: §C.1.1 · pp. 43–44

C.1.2 Proof of Theorem 16

  • For Theorem 16, the proof handles a non-Lipschitz function by decomposing it into Lipschitz terms and applying a variance bound.It uses Theorem 11-(ii) and the counterpart of Corollary 28 derived from Lemma 27.Source: §C.1.2 · p. 45

C.2 Bounds on u(4)

  • For non-increasing step sizes tending to zero with divergent cumulative sum, the proof establishes constants C1 and C2 independent of the step-size sequence for bounds on u(4)_0,n.The conditions are lim γ_k = 0 and lim Γ_k = +∞.Source: §C.2 · pp. 45–46

C.3 Proof of Theorem 17 … D.1 Coupling

The proofs derive concentration bounds for weighted empirical measures of Lipschitz and bounded measurable functions using Laplace-transform estimates. The technical appendix also constructs a coupling kernel that preserves the marginals of two Euler transitions.

C.3 Proof of Theorem 17

  • The proof derives Laplace-transform bounds by decomposing the empirical-measure functional into martingale increments and applying the Gaussian log-Sobolev inequality to Euler transitions.The transitions have Gaussian noise with covariance 2γ Id.Source: §C.3 · pp. 47–49
  • Lemma 30 bounds the Laplace transform of a Lipschitz function of the Euler approximation over indices n through p under the same step-size conditions.Its proof uses the Gaussian log-Sobolev inequality.Source: §C.3 · pp. 48–49
  • Under H1 and H2 with non-increasing step sizes satisfying γ1 ≤ 2/(m + L), Proposition 32 bounds the Laplace transform used to prove Theorem 17’s concentration result.The proof applies Markov's inequality to the proposition's bound.Source: §C.3 · p. 49

C.4 Proof of Theorem 18

  • Theorem 18’s proof decomposes the bounded-measurable-function terms into Lipschitz and bounded measurable components before combining their Laplace-transform bounds.The decomposition is combined with the empirical-measure decomposition used in the preceding proof.Source: §C.4 · pp. 49–50
  • The proof combines the component bounds to control the empirical measure’s Laplace transform, then uses Markov’s inequality to obtain a deviation bound.The passages identify these steps but do not show the bound’s formula.Source: §C.4 · p. 50

D.1 Coupling

  • For every x, y, and k ≥ 1, Kk((x, y), ·) is a transference plan between the Euler-transition distributions Pk(x, ·) and Pk(y, ·).The proof verifies the second marginal when hk(x) differs from hk(y) using an orthogonal change of variables.Source: §D.1 · pp. 50–51

D.2 Distribution of hitting time of 0 for Ornstein-Ulhenbeck processes

This section defines the first time a one-dimensional Ornstein–Uhlenbeck process reaches 0 and cites a proposition giving its distribution.

D.2 Distribution of hitting time of 0 for Ornstein-Ulhenbeck processes

  • The hitting time is ˜T0 = inf{t ≥ 0 : ˜Ut = 0} for the process started at ˜U0 = a.The process follows d˜Ut = −θ˜Utdt + σdB1_t, with θ, σ > 0.Source: §D.2 · p. 51
  • Proposition 34 gives the distribution of this hitting time for all a ∈ R, θ, σ > 0, and t > 0 [2].The passage identifies Φ as the cumulative distribution function of the standard normal distribution.Source: §D.2 · p. 52
Loading 1605.01559v4…