Gradient-free Riemannian Langevin Sampler Ricardo Baptista∗, Olivier Zahm† July 9, 2026
arXiv:2607.07519v1 [cs.LG] 8 Jul 2026
Abstract We address the problem of efficiently sampling multimodal probability distributions, where standard Markov Chain Monte Carlo methods often suffer from poor mixing and mode trapping. To mitigate these issues, we propose Gradient-free Riemannian Langevin Sampler (GRiLS), a novel proposal that improves exploration without requiring gradient evaluations of the target density. Our approach introduces a Riemannian metric which reshapes the local geometry in order to facilitate transitions across modes. The resulting gradient-free MCMC algorithm is particularly suitable for complex, computationally expensive targets where derivatives are unavailable or impractical. The GRiLS proposal requires knowing the mean and covariance of the target density, which we estimate using an ensemble of interacting particles. Empirical results on multimodal benchmarks demonstrate that GRiLS achieves improved mixing compared to existing gradient-based and gradient-free MCMC approaches.
Keywords: Gradient-free MCMC, Riemannian Langevin dynamics, Lamperti transform, Interacting particle system.
1
Introduction
The accurate and efficient sampling of high-dimensional probability distributions is a central task in computational statistics. For sampling distributions on Rd whose density µ is known up to a normalization constant, Markov Chain Monte Carlo (MCMC) remains a predominant method [25, 54]. Given a proposal density q(·|·), MCMC algorithms build a Markov chain {x1 , x2 , . . .} by drawing a proposal candidate x† ∼ q(·|xk ) and directly accepting this as the next state xk+1 = x† (unadjusted algorithms) or rejecting it with a certain probability and setting xk+1 = xk (Metropolis adjusted algorithms). MCMC offers strong theoretical guarantees, especially when the target density µ is logconcave [10, 17, 18]. However, one major difficulty in practical applications is when the target is not log-concave, e.g., under multimodality. In such scenarios, the Markov chain may be trapped in a mode, leading to poor mixing, high autocorrelation, and large bias. The computational efficiency of various samplers, often quantified by the time required to traverse energy barriers between modes, deteriorates exponentially as the modes become more separated [34]. ∗ †
Department of Statistical Sciences, University of Toronto, Canada, [email protected] UGA, Inria, CNRS, Grenoble INP*, LJK, 38000 Grenoble, France, [email protected]
1
To overcome these mixing challenges, practitioners have explored a variety of techniques, including tempered MCMC (also known as Parallel Tempering [49]), transport map accelerated MCMC [15, 52], adaptive biasing force methods [12, 62], dimension reduction techniques [45, 67] and preconditioned Langevin dynamics [16, 28, 42, 66] to name just a few. We focus here on the latter approach, which defines the proposal density via a time discretization of a Langevin dynamic that has been preconditioned in order to improve the convergence of the continuous-time dynamic toward equilibrium. This preconditioning is achieved by equipping Rd with a suitable Riemannian metric that locally reduces the geodesic distances between modes, thereby facilitating the transitions of particles across different regions of high probability. Given an arbitrary field of symmetric definite matrices W : x 7→ W (x) ∈ Rd×d , we endow Rd with the Riemannian metric ⟨u, v⟩x = u⊤ W (x)−1 v. The Riemannian Langevin dynamic is given by p (1) dXt = divW (Xt ) + W (Xt )∇ ln µ(Xt ) dt + 2W (Xt )dBt , P where divW (x) = ( dj=1 ∂j Wi,j (x))1≤i≤d is the divergence of W and (Bt )t≥0 is the standard Brownian motion in Rd , see e.g. [11, 36]. Popular choices for W include the constant metric [30, 33, 64], the inverse negative Hessian of the log-density [28], or a suitable positive-definite approximation of it when the target is not strongly log-concave [6, 7]. More generally, the Hessian of any arbitrary strongly convex function may be used, giving rise to the so-called mirrored Langevin dynamics [35, 68]. It was proven in [16] that an optimal Riemannian metric exists (optimal in the sense of the fastest convergence of Xt towards equilibrium) and can be expressed as the Hessian of a strongly convex function; however, its exact computation remains intractable in general. In this work, we propose a gradient-free Riemannian Langevin sampler for targeting multimodal densities. To define the metric, we choose W (x) =
ν(x) Σ, µ(x)
(2)
where ν = N (m, Σ) is the Gaussian density with same mean m ∈ Rd and covariance Σ ∈ Rd×d as µ. Although it might not be optimal in the sense of [16], we show that this choice can greatly improve the convergence to equilibrium as compared to other choices, provided that ν “covers” the modes of µ by satisfying sup( µν ) < ∞. This choice has two key advantages. First, if a particle x lies between two modes of µ, we have W (x)−1 ≈ 0 which reduces the geodesic distances between the modes. Second, the choice (2) leads to a gradient-free Langevin dynamic which does not involve the term ∇ ln µ(Xt ) in (1). This offers a significant advantage when gradients of the target densities are computationally prohibitive or intractable to evaluate. Such gradient-free dynamics are receiving growing attention, see [16] when W is the optimal metric, [19, 20] for densities defined on the torus, and [9, 24, 39] in the setting of ensemble-based Langevin samplers. To define the proposal density based on the Riemannian Langevin dynamics (1) with (2), we will consider a Lamperti transform [40] and an appropriate time discretization with stepsize ∆t > 0. The resulting proposal sample x† ∼ q(·|θk ) is given by: q ξk ∼ N (0, Σ), ν(xk ) † 2 , θk = exp −∆t µ(x (3) x = θk xk + (1 − θk )m + 1 − θk ξk , where k) ν(x) ∝ exp − 12 ∥x − m∥2Σ−1 . 2
We denote the MCMC algorithm based on this proposal by Gradient-free Riemannian Langevin Sampler (GRiLS). For the proposal in (3), the target density µ can be known up to a normalizing constant, since the latter can be absorbed in the stepsize ∆t. Interestingly, letting θk = θ0 ∈ (0, 1) be a constant parameter recovers the preconditioned Crank-Nicolson proposal (pCN) [14]. In this sense, the proposal in (3) can be interpreted as an extension of the pCN scheme, which is consistent with the Langevin dynamics in (1) as ∆t → 0. In practice, we approximate the mean m and covariance Σ in (3) using approximate samples from µ. These samples can either be obtained from previous MCMC iterations {x1 , . . . , xk }, following adaptive MCMC strategies [3, 30, 33], or from an ensemble of N ≥ 1 particles {x1k , . . . , xN k } at step k, as in more recent approaches [8, 23, 41, 61]. We focus on the ensemble-based approach and introduce a block-ensemble version of GRiLS (BE-GRiLS) in order to improve the computational efficiency. Given a partition ∪Pℓ=1 B ℓ = {1, . . . , N }, each particle xik within a block i ∈ B ℓ is updated using the sample mean and covariance computed from all other blocks, that is: ! X X 1 1 xik and Σℓk = (xik )(xik )⊤ − (mℓk )(mℓk )⊤ . (4) mℓk = ℓ N − #B ℓ N − #B ℓ ℓ i∈B /
i∈B /
As pointed out in [61], this strategy enables the parallel update of all particles within a block, thereby improving the computational efficiency of the algorithm. It is worth mentioning that the resulting sampler BE-GRiLS is similar to the Consensus-based sampler (CBS) [8] in several aspects, as they both use the ensemble mean and covariance to update each particle as in (3). The main difference is that, in CBS, the parameter θk = θ0 is taken to be constant during the iterations and hence, as for pCN, CBS is not consistent with a Langevin dynamic targeting µ when ∆t → 0. Moreover, instead of using a Metropolis correction, CBS employs importance sampling to estimate the mean and covariance of µ. The rest of the paper is structured as follow. In Section 2 we motivate the choice of metric (2) for preconditioning the Riemannian Langevin dynamic (1). In Section 3 we propose a time-integration scheme which leads to the proposal in (3). In Section 4, a spectral analysis of MCMC algorithms relates the convergence rate of MCMC algorithms with the one of the continuous Langevin dynamic (1). Finally, in Sections 5 and 6 we illustrate the performance of the resulting MCMC algorithm on several one-dimensional and multivariate benchmark problems.
2
Riemannian Langevin Dynamics
Let µ be a probability density on Rd defined by µ(x) ∝ exp(−V (x)), where V : Rd → R is a smooth potential function. A classic approach for sampling µ is to consider the overdamped Langevin dynamics √ (5) dXt = −∇V (Xt )dt + 2dBt , where Bt is the d-dimensional Brownian motion. Under standard regularity assumptions on V , the invariant density of this dynamic is µ, i.e., Xt converges in law to µ as t → ∞. 3
constant
(inverse) Metric W (x) = Id W (x) = Covµ
Comments Standard Langevin dynamic (5) Covariance-based preconditionning [30], see [24] for an ensemblebased covariance estimation. Inverse of Fisher matrix of µ [64]
W (x) = Eµ [∇ ln µ ∇ ln µ⊤ ]−1
x-dependent
−1 W (x) = (I(x) + Σ−1 0 )
W (x) = SoftAbs[−Hess ln µ(x)] W (x) = (Hess ψ(x))−1 W (x) = µ(x)−1 Id W ∈ arg minW :x7→W (x) C(µ, W ) ν(x) Covµ W (x) = µ(x)
When dµ(x) ∝ Ly (x)dµ0 (x) where µ0 = N (0, Σ0 ), use the Fisher information matrix I(x) = EY |x [∇ ln LY (x)∇ ln LY (x)⊤ ] of the likelihood Ly (x) [28], see also [37] for generalization to a class of latent variable models, and [5, 7] for efficient approximate computation in high-dimension. Soft-absolute value of local Hessian of the log-density [6] Mirrored Langevin dynamic, for ad-hoc strictly convex ψ [35, 68] Derivative-free dynamics (for µ defined on the torus Td ) [20, 42] Optimizing the Poincaré constant [16, 42] For ν = N (m, Σ) with m = Eµ [X] and Σ = Covµ (present paper)
Table 1: Possible choices of metric W An alternative dynamic to (5) that does not change the invariant density is p dXt = divW (Xt ) − W (Xt )∇V (Xt ) dt + 2W (Xt )dBt ,
(6)
where W : Rd → S+d is any smooth field taking value in S+d ⊂ Rd×d , the set of symmetric P semi-definite positive matrices, see [11, 36]. Here, divW = ( dj=1 ∂j Wi,j )1≤i≤d is the √ √ √ divergence of W and W is any square root of W such that W (x) W (x)⊤ = W (x). This dynamic corresponds to a Langevin diffusion process on the Riemannian manifold Rd endowed with the metric induced by W −1 , meaning ⟨u, v⟩X := u⊤ W (X)−1 v for any v, u ∈ Rd . Over the past two decades, several choices have been proposed for W , ranging from constant matrix fields to location-dependent fields; see Table 1 for a non-exhaustive overview of possible metrics.
2.1
Optimal Metric
The idea proposed in [16, 42] is to identify W by optimizing the convergence rate of the preconditioned dynamics in (6). Denoting the density of Xt that solves (6) by µt , the chi-square divergence χ2 (µt ||µ) := Varµ (µt /µ) satisfies χ2 (µt ||µ) ≤ e−2t/C(µ,W ) χ2 (µ0 ||µ),
(7)
for any t ≥ 0 and any initial condition µ0 ≫ µ. Here, C(µ, W ) ≥ 0 denotes the Poincaré constant, which is defined as the smallest constant such that the Riemannian Poincaré inequality Z Varµ (f ) ≤ C(µ, W ) ∥∇f (x)∥2W (x) µ(x)dx, (8) holds for any smooth function f : Rd → R; see [4] for more details. Here we use the notation ∥v∥2A = v ⊤ Av. It turns out that C(µ, W ) is also the smallest constant so that the inequality (7) holds for any t ≥ 0 and any µ0 ≫ µ, therefore it characterizes exactly the exponential convergence of the dynamic (6). Based on this observation, [16, 42] propose to identify the matrix field W which yields the best convergence rate by solving min
d W :Rd →S+ Eµ [tr(W )]=tr(Covµ )
4
C(µ, W ).
(9)
The constraint Eµ [tr(W )] = tr(Covµ ) is introduced in [16] to fix the scaling of W in order to prevent arising from C(µ, αW ) = α1 C(µ, W ) → 0 as α → ∞. Here, R a trivial solution R Covµ = (x − m)(x − m)⊤ µ(x)dx and m = xµ(x)dx are the covariance and the mean of µ. While [42] considers alternative normalization constraints, the choice Eµ [tr(W )] = tr(Covµ ) offers significant advantages. First, by testing the Poincaré inequality (8) with affine functions f , this constraint yields the lower bound C(µ, W ) ≥ 1.
(10)
Second, under some assumptions on µ, it is shown in [16, Section 2] that an optimal solution Wµopt to (9) exists with C(µ, Wµopt ) = 1. This optimal field Wµopt is given by Wµopt (x) = Hess φ⋆ (x)−1 ,
(11)
where φ⋆ : supp(µ) → R is the strictly convex and smooth function corresponding to the convex conjugate of the moment map φ of µ. We refer to [13] for more details on moment maps. Because Wµopt is the inverse Hessian of a strictly convex function, the resulting dynamic (6) is a Mirrored Langevin dynamic [35, 68] of the form Xt = ∇ψ ⋆ (Yt ) p with dYt = −∇V (Xt )dt + 2Hess ψ(x)dBt , with an optimal choice for the mirror map ψ = φ⋆ . Third, [16, Theorem 3.2.] shows that Wµopt is actually a Stein Kernel 1 , and therefore the dynamic (6) can be shown to further simplify as q dXt = −(Xt − m)dt + 2Wµopt (Xt )dBt , where we recall that m is the mean of µ. Remarkably, this dynamic is gradient-free in the sense that it no longer involves the gradient of the potential V of µ ∝ exp(−V ). Computing the optimal metric Wµopt (or the associated moment map φ) for a general measure µ, however, is a difficult task. While [16] proposes a gradient-descent algorithm for solving (9), computing the gradient of W 7→ C(µ, W ) requires the eigendecomposition of a diffusion operator on Rd , which is not computationally tractable for dimensions d ≥ 3. Nonetheless, there are two cases where Wµopt can be computed in closed form: • For d = 1, if µ is supported on a convex domain, Proposition 2 in [16] shows Z ∞ Z 1 opt (t − m)µ(t)dt, m = xµ(x)dx. (12) Wµ (x) = µ(x) x See Figure 1 for the representation of Wµopt when µ is a mixture of two Gaussian. This object has received growing attention in the probability literature [21, 26, 58] and in sensitivity analysis [31, 32, 57, 60]. • For a probability density µ defined on the torus Td = (R/Z)d , i.e., µ(x + k) = µ(x) for all x ∈ Rd and k ∈ Zd , the optimal metric cannot be expressed as in (11). In fact, there is no moment map φ for such µ because convex functions on Td are necessarily constant. Nonetheless, [42, Proposition 10] states that Wµopt,hom (x) = 1
1 Id , µ(x)
d d A Stein Kernel for µ is any matrix field R W : R → S+ such that d smooth function f : R → R, where m = xdµ.
5
(13) R
(x − m)f dµ =
R
W ∇f dµ for all
(x) W opt(x) (x) W (x) (x)
10.0
7.5
5.0
2.5
0.0 x
2.5
5.0
7.5
10.0
Figure 1: Optimal metric Wµopt as in (12) (red curve) for the Gaussian mixture µ = 1 N (−4, 1) + 12 N (+4, 1) (bashed black curve). The blue curve represents the sub-optimal 2 metric Wµ as in (15). For better visualisation, the metrics are multiplied by µ. solves (9) in the homogenized limit. That is, when replacing µ and W in (9) with µ♯,k (x) = µ(kx) and W♯,k (x) = W (kx), respectively, and letting k ∈ N tend to +∞ yields this solution; see [1] for more details. The metric Wµopt,hom is used in [20] to p derive a gradient-free Langevin dynamic dXt = 2/µ(Xt )dBt from (6). The two analytical solution (12) and (13) will guide our construction of a suboptimal, but computationally practical, metric in the next subsection.
2.2
Gradient-Free Dynamics
For simplicity, we assume the probability density µ is fully supported on Rd , so that µ(x) > 0 for all x ∈ Rd . With µ(x) ∝ e−V (x) , the drift term in (6) can then be written as divW (x) − W (x)∇V (x) =
1 div(W (x)µ(x)). µ(x)
(14)
Thus, any field of the form W (x) = A(x)/µ(x) for some field A : Rd → S+d results in a Langevin dynamic (6) in which the gradient of V does not appear: only the divergence of A is required. This is the case for the optimal metrics (12) and (13) mentioned earlier. We now construct a simple and tractable field A that yields a suitable preconditioner for the Langevin dynamic. Based on Figure 1, two important features of Wµopt are that Wµopt (x)µ(x) is unimodal and that it encompass the support of µ. Hence, we propose exp(− 21 ∥x−m∥2 −1 ) Σ √ ν(x) = , ν(x) 2π det Σ Σ, where (15) Wµ (x) = m = Eµ [X], µ(x) Σ = Covµ . Figure 1 presents Wµ and the optimal Wµopt given by (15) and (12), respectively. While the choice in (15) might not be optimal with respect to optimization problem (9), it is a natural and convenient choice in practice. In addition, setting Σ = Covµ ensures that Wµ satisfies the normalization constraint of (9) by construction, meaning Eµ [tr(Wµ )] = Eν [tr(Σ)] = tr(Covµ ). Combining (14) with (15), the preconditioned dynamic in (6)
6
simplifies to ν(Xt ) dXt = − (Xt − m)dt + µ(Xt )
s 2
ν(Xt ) ΣdBt . µ(Xt )
(16)
This Langevin dynamic does not involve the gradient of the potential V , while still admitting µ as its invariant measure by construction. The following proposition provides a simple upper bound for the Poincaré constant C(µ, Wµ ). Proposition 1. Let µ be a probability density that is fully supported on Rd , and let Wµ take the form in (15). Then µ(x) . x∈Rd ν(x)
1 ≤ C (µ, Wµ ) ≤ sup
(17)
Proof. The left inequality 1 ≤ C (µ, Wµ ) is a direct consequence of the normalization Eµ [tr(Wµ )] = tr(Covµ ); see (10). For the upper bound, we use that Z Z µ 2 Varµ (f ) = min (f − α) dµ ≤ (sup ν ) min (f − α)2 dν = (sup µν ) Varν (f ), α∈R
α∈R
for any f . Then, the Gaussian Poincaré inequality Varν (f ) ≤ Eν [∥∇f ∥2Σ ] (see [4]) yields Varµ (f ) ≤ (sup µν )Eν [∥∇f ∥2Σ ] = (sup µν )Eµ [∥∇f ∥2Wµ ], which proves (17). We now illustrate Proposition 1 on the one-dimensional Gaussian mixture 1 1 µh = N (−h, 1) + N (+h, 1), 2 2 for some h ≥ 0. It is well known that, for such multimodal densities, the Poincaré constant associated with the standard Langevin dynamic (5) grows exponentially with the height of the energy barrier between modes (here, h2 ); see [34, 48]. We show in Appendix B that the lower bound 2
eh /3 C(µh , Covµh ) ≥ , 1 + h2
(18)
holds for any h ≥ 1/2. The sharpest upper bound for C(µh , Covµh ) to the best of our 4h2
3+e knowledge is derived in [59, Section 4.1] and yields C(µh , Covµh ) ≤ 4(1+h 2 ) . Instead, by considering the metric Wµh in (15), the Poincaré constant C(µh , Wµh ) associated with the preconditioned Langevin dynamic (16) satisfies (17)
√ µh (x) ≤ exp(1) 1 + h2 , x∈Rd νh (x)
C(µh , Wµh ) ≤ sup
(19)
for any h ≥ 0, where νh = N (Eµ [X], Covµh ) = N (0, 1 + h2 ). We give a proof of inequality (19) in Appendix A. Both Wµh and Covµh satisfy the normalization constraint R tr(Wµh )dµh = tr(Covµh ). Therefore, they can be compared as suboptimal solutions to (9): while the Poincaré constant associated with the standard Langevin dynamic (5) grows at least exponentially in h2 , the one of the preconditioned Langevin dynamic (16) growths at most linearly in h for h ≫ 1. 7
Remark 2.1 (Generalization to high-dimensional mixtures). The bound (19) generalizes to P Gaussian mixtures with N ≥ 2 components in dimension d ≥ 1. Consider µN h (x) = N 1 Letting h ≥ 0 be the smallest constant such that i=1 µi (x) where µi = N (mi , Id ). P N 1 N ∥mi − m∥ ≤ h for all i, where m = N N i=1 mi is the mean of µh , we have C(µ, WµNh ) ≤ exp(N/2)(1 + h2 )d/2 .
(20)
The proof is given in Appendix A. While this bound still grows polynomially in h, it depends exponentially on d and N . We believe this behavior is an artifact of the proof technique, which relies on the bound (17). This bound is likely overly pessimistic, particularly in high-dimensional settings. Remark 2.2 (Generalization to non-Gaussian reference densities). The previous development can be generalized to a non-Gaussian reference density ν. Given an approximation ν to µ, the choice Wµν (x) = Ω
ν(x) Covν , µ(x)
where Ω =
tr(Covµ ) , tr(Covν )
yields the following preconditioned Langevin dynamic s ν(Xt ) ν(Xt ) dXt = −Ω Covν ∇ log ν(Xt )dt + 2Ω Covν dBt . µ(Xt ) µ(Xt ) This dynamic does not involve the gradient of the target µ, but only the one of the approximation ν. By construction we have Eµ [tr(Wµν )] = tr(Covµ ) so that Wµν satisfies the normalization constraint in (9). Thus, similarly to Proposition 1, we can show that 1 ≤ µ) C(µ, Wµν ) ≤ tr(Cov (sup µν )C(ν, Covν ). This analysis suggests to finding ν by minimizing tr(Covν ) µ) ν → tr(Cov (sup µν ). We leave this for future work. tr(Covν )
3
Sampling Algorithms
In this section we consider different numerical time integration methods for the dynamics (16) in order to construct a MCMC proposal. Given a stepsize ∆t > 0, the proposal approximately draws a sample xk+1 from X∆t , where Xt is the solution to (16) initialized at X0 = xk . A first approach is to apply the Euler-Maruyama scheme directly to (16), leading to s ν(xk ) ν(xk ) (xk − m)∆t + 2∆t ξk , ξk ∼ N (0, Σ). (21) xk+1 = xk − µ(xk ) µ(xk ) This naive discretization is problematic when the ratio ν/µ exhibits strong spatial variations. Indeed, if xk is between two modes of µ, we have µ(xk ) ≪ ∆tν(xk ). Hence, the mean and covariance of the Gaussian vector xk+1 |xk defined above will explode. The following Lamperti transformation avoids this behavior. Proposition 2 (Lamperti transformation). Let µ be a probability density that is fully supported on Rd , and let ν = N (m, Σ) be the Gaussian density with mean m and covariance Σ. Let Yt be the solution of the Ornstein-Uhlenbeck process √ dYt = −(Yt − m)dt + 2ΣdBt , (22) 8
with initial condition X0 = Y0 , and consider the random variable Z t ν(Ys ) τ (t) = ds. 0 µ(Ys )
(23)
d
Then, for all t ≥ 0, we have Yτ (t) = Xt , where (Xt )t≥0 is the solution to (16).
Proof. Applying the stochastic time change in [51, Theorem 8.5.1] yields the result; see also [50] for the proof. Proposition 2 provides an alternative way to draw a sample xk+1 from X∆t |X0 = xk by following the steps: (i) draw a trajectory (Yt )t≥0 of the Ornstein-Uhlenbeck (OU) process (22) initialized at Y0 = xk and (ii) compute the (deterministic) integral τ (∆t) in (23) conditioned on (Yt )t≥0 . Since Proposition 2 ensures Yτ (∆t) has the same law as X∆t |X0 = xk , we can set xk+1 = Yτ (∆t) . In practice, we need to numerically approximate τ (∆t) in step (ii). We propose to use the first-order quadrature scheme: Z ∆t ν(Y0 ) ν(xk ) ν(Ys ) ds ≈ ∆t = ∆t . (24) τ (∆t) = µ(Ys ) µ(Y0 ) µ(xk ) 0 Notably, the quadrature scheme does not depend on the values of Yt for t > 0. Next, using the closed form solution of the OU process Ys |Y0 ∼ N (m + e−s (Y0 − m), (1 − e−2s )Σ), we can draw a sample xk+1 with law (Y∆tν(xk )/µ(xk ) |Y0 = xk ) as follows ( q ξk ∼ N (0, Σ), xk+1 = m + θk (xk − m) + 1 − θk2 ξk , where (25) ν(xk ) θk = exp(−∆t µ(x ). k) Contrarily to the Euler-Maruyama scheme (21), µ(xk ) ≪ ∆tν(xk ) implies θk ≪ 1 so that the Gaussian vector xk+1 |xk defined in (25) will be approximately drawn from N (m, Σ), which has finite mean and covariance. Notably, as ∆t → 0, the Taylor expansion θk = ν(xk ) + O(∆t2 ) shows that (25) recovers the Euler-Maruyama scheme (21) up to 1 − ∆t µ(x k) second order terms O(∆t2 ). Remark 3.1. Classical error analysis of quadrature scheme (24) yields " 2 #1/2 ν(Y0 ) E τ (∆t) − ∆t = O(∆t3/2 ). µ(Y0 ) Higher order quadrature methods can be used to improve the estimation of τ (∆t) at the cost of evaluating the gradient of the integrand ν/µ. We refer to [38, Chapter 5] for more details. Remark 3.2. The Lamperti transform can be related to an importance sampling strategy. Indeed, for any integrable function f , we can write Z Z µ(y) ν(y)dy f (x)µ(x)dx = f (y) ν(y) Z µ(Yt ) 1 T (For Yt solving (22)) = lim f (Yt ) dt T →∞ T t=0 ν(Yt ) Z −1 1 τ (T ) (t ← τ (t) as in (23)) = lim f (Yτ (t) )dt T →∞ T t=0 Z 1 T = lim f (Yτ (t) )dt T →∞ T t=0 9
where, for the last step, we Rused the fact that τ −1 (T ) ∼ T when T → ∞ (i.e., by letting f (x) = 1, we deduce 1 = f (x)µ(x)dx = limT →∞ τ −1 (T )/T and thus τ −1 (T ) ∼ T ). Instead of using the ratio µ(y)/ν(y) as an importance weight, as classically done by the importance sampling method, the Lamperti transform stretches time via t ← τ (t) in d order to correct for the bias introduced by using Yt in place of Xt = Yτ (t) .
3.1
MCMC Proposal
Given a target probability density µ (possibly unnormalized) and a proposal density q(·|·), our MCMC algorithm builds a Markov chain {x0 , x1 , . . .} ⊂ Rd by iteratively proposing a candidate x† ∼ q(·|xk ). For an unadjusted algorithm, one always accepts this candidate xk+1 = x† . For Metropolis adjusted algorithms, one accepts the candidate xk+1 = x† with probability µ(x† )q(xk |x† ) † . (26) α(x |xk ) = min 1, µ(xk )q(x† |xk ) If the candidate is rejected, then xk+1 = xk . This accept/reject step ensures that the Markov chain {xk }k≥0 admits µ as an invariant measure [25, 54]. Motivated by the scheme in (25), we introduce the proposal density ∥x† − m − θk (xk − m)∥2Σ−1 1 † exp − q(x |xk ) = p , (27) 2(1 − θk2 ) 2πdet((1 − θk2 )Σ) ν(xk ) ). In this definition, m and Σ are not necessarily the mean and where θk = exp(−∆t µ(x k) the covariance of the density µ, which are typically unknown. In practice, m and Σ can be defined as the sample mean and sample covariance as described later in Section 3.2. Alternatively, as done in [14], when µ ∝ Lµ0 is the posterior density of a Bayesian inverse problem with prior µ0 and likelihood L, then m and Σ can be chosen to be the mean and covariance of µ0 . With this proposal density, the acceptance rate simplifies to ∥x† −m−θ (x −m)∥2 k k Σ−1 d/2 exp − †2 † 2 2(1−θ µ(x ) 1 − θ k) † , (28) α(x |xk ) = min 1 , ∥xk −m−θ† (x† −m)∥2 −1 µ(xk ) 1 − θk2 Σ exp − 2(1−θ†2 )
†
ν(x ) where θ† = exp(−∆t µ(x † ) ). The resulting MCMC algorithm is described in Algorithm 1. Let us make a few comments on this algorithm. First, in order to compute θk and θ† , there is no need to know the normalizing constant of µ nor of ν as they can both be absorbed in the parameter ∆t. Second, as mentioned earlier, the proposal (27) is similar to the preconditioned Crank-Nicolson proposal (pCN) introduced in [14]. Indeed, replacing θk with a constant value θ ∈ (0, 1) yields the pCN proposal for which the acceptance µ(x† )ν(xk ) probability simplifies to α(x† |xk ) = min{1; µ(x † }. Third, the acceptance rate exhibits k )ν(x ) an exponential dependence on the dimension d. This behavior is the counterpart of using a location-dependent proposal variance. As the following proposition demonstrates, the effect of this dimensional dependence can be mitigated when ∆t is chosen sufficiently large.
Proposition 3. Let µ and ν be two probability densities on Rd (possibly unnormalized) d such that Ω = sup µν < ∞. Then, for any β > 1, taking ∆t ≥ Ω2 ln(1 + 2 ln(β) ) ensures d/2 1 − θ(x)2 −1 β ≤ ≤ β, 1 − θ(y)2 10
for any x, y ∈ Rd , where θ(x) = exp(−∆tν(x)/µ(x)). Proof. By definition of Ω = sup µν we have θ(x) ≤ exp(−∆t/Ω) for all x ∈ Rd . Then 1 − exp(− 2∆t ) Ω 1
!d/2
≤
1 − θ(x)2 1 − θ(y)2
d/2 ≤
1 1 − exp(− 2∆t ) Ω
!d/2 .
))−d/2 ≤ β, which is equivalent to ∆t ≥ − Ω2 ln(1 − It remains to show that (1 − exp(− 2∆t Ω −2/d β ). Using the inequality 1 − exp(−u) ≥ u/(1 + u), we have ! 2 ln(β) d 2 ln(β) −2/d d ≤ − ln − ln(1 − β ) = − ln 1 − exp − d = ln 1 + . 2 ln(β) 1 + 2 ln(β) d
d Then, the assumption ∆t ≥ Ω2 ln(1 + 2 ln(β) ) is sufficient to ensure ∆t ≥ − Ω2 ln(1 − β −2/d ), which concludes the proof.
Remark 3.3 (Independence Sampler). While Proposition 3 shows that for large enough †2 ∆t, the ratio ( 1−θ )d/2 in (28) can be uniformly bounded, it is also worth mentioning 1−θk2 that ∆t ≫ Ω = sup µν simplifies drastically the proposal density. Indeed, we have θk ≤ exp(−∆t/Ω) → 0 uniformly on k when ∆t → ∞, and therefore d
q(·|xk ) −→ ν(·) = N (m, Σ). ∆t→∞
(29)
In this case, GRiLS becomes an Independence Sampler (IS) in which the proposal candidate is drawn from ν independently of the current state. This type of sampler is efficient if the proposal density ν is close to µ. Using the formula min{a, b} = (a + b − |a − b|)/2, the mean acceptance rate is given by ( X ∼ µ(·) α = E[α(Y |X)], Y ∼ q(·|X) Z (26) = min{µ(x)q(y|x); µ(y)q(x|y)}dxdy Z (29) −→ min{µ(x)ν(y); µ(y)ν(x)}dxdy ∆t→∞ Z 1 µ(x)ν(y) − µ(y)ν(x) dxdy =1− 2 = 1 − ∥µ ⊗ ν − ν ⊗ µ∥TV , (30) where ∥ · ∥TV denotes the total variation distance. We observe numerically in Section 5 that the mean acceptance rate of GRiLS is always greater than the above limit.
3.2
Block Ensemble GRiLS
In this section, we propose an ensemble version of GRiLS in which we estimate the mean m and the covariance Σ of µ using an ensemble {x1 , . . . , xN } of N particles in Rd . For simplicity, we represent this ensemble by the matrix X = (x1 , . . . , xN ) ∈ Rd×N . 11
Algorithm 1: GRiLS (Gradient-free Riemannian Langevin Sampler) Input: Initial state x0 ∈ Rd , target probability density µ (possibly unnormalized), mean m ∈ Rd , covariance Σ ∈ S+d , step size ∆t > 0, number of steps K ∈ N, Output: MCMC chain {x1 , . . . , xK } ⊂ Rd Algorithm GRiLS(x0 , µ, m, Σ, ∆t, K) 2 for k = 0 to K − 1 do ν(xk ) 3 Compute θk = exp(−∆t µ(x ), where ν(x) = exp(−∥x − m∥2Σ−1 ) k) 4 Draw a sample ξk ∼ N (0, Σ) and build the candidate proposal 1
x† = m + θk (xk − m) +
q
1 − θk2 ξk ,
†
5
ν(x ) Compute θ† = exp(−∆t µ(x † ) ) and, with probability
α(x† |xk ) = min
1;
µ(x† ) µ(xk )
∥x† −m−θ (x −m)∥2 k k Σ−1 d/2 exp − 2(1−θk2 ) 1 − θ†2 ∥xk −m−θ† (x† −m)∥2 −1 1 − θk2 Σ exp − 2(1−θ†2 )
accept the candidate by setting xk+1 = x† or reject by setting xk+1 = xk . 7 end 8 return {x1 , . . . , xK }; 6
We aim to build a MCMC chain {X1 , X2 , . . .} ⊂ Rd×N that targets the probability density in which the members are independent and drawn from the tensor product measure of the target µ⊗N (X) = µ(x1 ) · · · µ(xN ). (31) To do this, we consider a Metropolis-within-Gibbs algorithm, where Xk+1 is constructed by updating each column of Xk one after the other. One Gibbs iteration starts by initializing Xk→k+1 = Xk and, for ℓ = 1 . . . N , updates the ℓ-th column of Xk→k+1 to ℓ+1 ℓ N Xk→k+1 ← (x1k+1 , . . . , xℓ−1 k+1 , xk+1 , xk , . . . , xk ).
(32)
At this step, xℓk+1 can be generated by any Metropolis algorithm targeting the conditional ℓ+1 N ℓ density µ⊗N (xℓ |x1k+1 , . . . , xℓ−1 k+1 , xk , . . . , xk ) which, by (31), simplifies to µ(x ). After the N columns are updated, we set Xk+1 = Xk→k+1 and we move on to the next Gibbs iteration k ← k + 1. We now detail the Metropolis step which generates the update xℓk+1 in (32). For simplicity, we denote the current ensemble Xk→k+1 excluding its ℓ-th element by ℓ+1 −ℓ N d×(N −1) = (x1k+1 , . . . , xℓ−1 . Xk→k+1 k+1 , xk , . . . , xk ) ∈ R
(33)
−ℓ ), we draw a candidate sample Given a proposal density of the form q(x† |xℓk , Xk→k+1 −ℓ x† ∼ q(·|xℓk , Xk→k+1 ) and accept it in the ensemble xℓk+1 = x† with probability ) ( † q(xℓ |x† , X −ℓ ) µ(x ) k k→k+1 −ℓ α(x† |xℓk , Xk→k+1 ) = min 1; . (34) −ℓ µ(xℓk ) q(x† |xℓk , Xk→k+1 )
12
Otherwise, we reject x† and set xℓk+1 = xℓk . Defined that way, this Metropolis-withinGibbs algorithm guarantees detailed balance, which ensures µ⊗N is the stationary density of the chain {Xk }k≥1 ; see [56] for more details. We propose to employ the GRiLS proposal (27) in which m and Σ are defined as the −ℓ sample mean and the sample covariance of the ensemble Xk→k+1 defined in (33). The resulting proposal is given by ! ∥x† − mℓk − θkℓ (xℓk − mℓk )∥2(Σℓ )−1 1 −ℓ k exp − , q(x† |xℓk , Xk→k+1 )= p 2(1 − (θkℓ )2 ) 2πdet((1 − (θkℓ )2 )Σℓk ) (35) ℓ ℓ ℓ N −1 where θk = exp(−∆tν(xk )/µ(xk )). Using 1N −1 ∈ R to denote the vector of ones, the mean and covariance of the proposal can be computed as ℓ−1 (Xk→k+1 )1N −1 ℓ mk =
N −1
and
ℓ−1 ℓ−1 (Xk→k+1 )(Xk→k+1 )⊤ ℓ Σk = − (mℓk )(mℓk )⊤ .
N −1
(36)
By excluding the ℓ-th member in the computation of mℓk and Σℓk , the acceptance proba−ℓ ) is the same as the one in (28) with m = mℓk and Σ = Σℓk . bility α(x† |xℓk , Xk→k+1 So far, we have updated one particle at a time xℓk → xℓk+1 . In the same way, we now derive an algorithm which updates multiple particles in parallel {xik }ℓ∈B → {xik+1 }ℓ∈B for some subset B ⊂ {1, . . . , N }, thereby reducing the number of iterations by a factor #B and improving the overall computational efficiency. To preserve detailed balance for the resulting MCMC algorithm, we estimate the mean and the covariance using the particles that are not in B. We refer to this algorithm as Block-Ensemble GRiLS (BE-GRiLS). The steps of the procedure are summarized in Algorithm 2. Remark 3.4 (Avoiding covariance matrix factorizations). The ensemble covariance matrix Σℓk in (36) can be expressed equivalently in matrix form as 1
1⊤
ℓ−1 ℓ−1 N −1 (Xk→k+1 )(IN −1 − N −1 )(Xk→k+1 )⊤ N −1 ℓ Σk = , N −1
where IN −1 is the identity matrix of size N −1. This expression allows for straightforward sampling of a Gaussian vector ξkℓ ∼pN (0, Σℓk ), which is required to generate the candidate proposal x† = mℓk + θk (xk − mℓk ) + 1 − θk2 ξkℓ . In particular, for Z ∼ N (0, IN −1 ), we have 1⊤
1
N −1 (IN −1 − N −1 ) N −1 ℓ−1 ℓ √ ξk = (Xk→k+1 ) Z N −1
∼
N (0, Σℓk ).
Thus, we can sample from N (0, Σℓk ) without computing any factorization of Σℓk .
13
Algorithm 2: BE-GRiLS (Block-Ensemble GRiLS) d×N Input: Initial ensemble X0 = (x10 , . . . , xN , target probability density µ 0 ) ∈ R (possibly unnormalized), partition 0 = N0 ≤ . . . ≤ NP = N , step size ∆t > 0, number of steps K ∈ N, Output: MCMC chain {X1 , . . . , XK } ⊂ Rd×N Algorithm BE-GRiLS(X0 , µ, (Nℓ )Pℓ=0 , ∆t, K) 2 For i ≤ d do ## Gibbs iteration 3 Initialize Xk→k+1 := Xk 4 For ℓ = 1 to P do ## Update the ℓ-th block 5 Remove the ℓ-th block of ensemble from Xk→k+1
1
6 N
(N )+1
−ℓ ℓ−1 Xk→k+1 = (x1k+1 , . . . , xk+1 , xk ℓ 7
Compute N (ℓ) = N − (Nℓ − Nℓ−1 ) and the ensemble mean and covariance mℓk =
8 9 10 11
, . . . , xN k )
ℓ−1 (Xk→k+1 )1N (ℓ) N (ℓ)
and Σℓk =
ℓ−1 ℓ−1 (Xk→k+1 )(Xk→k+1 )⊤ − (mℓk )(mℓk )⊤ N (ℓ)
For i = (Nℓ−1 ) + 1 to Nℓ do ## Metropolis (parallelizable) xik+1 = GRiLS(xik , µ, mℓk , Σℓk , ∆t, 1) end Update the ℓ-th block of ensemble (N )+1
ℓ ℓ Xk→k+1 = (x1k+1 , . . . , xN k+1 , xk
, . . . , xN k )
end Update the ensemble Xk+1 = Xk→k+1 14 end 15 return {X1 , . . . , XK }; 12 13
4
Spectral Analysis of MCMC Algorithms
In this section we compare the efficiency of GRiLS against other MCMC algorithms via a spectral analysis of their respective transition kernels. The transition kernel P (·, ·) of a homogeneous Markov chain {x0 , x1 , x2 , . . .} ⊂ Rd is a function defined on Rd × B(Rd ) such that P (x, ·) is the probability measure of xk+1 |xk = x for all x ∈ Rd . Let µk denote the probability measure of xk ∼ µk . Then, we can write µk+1 = µk P,
(37)
R where we use the notation ρP (·) = x∈Rd P (x, ·)ρ(dx) for any probability measure ρ on B(Rd ). We assume that there exists a unique invariant measure µ∞ such that µ∞ = µ∞ P . The transition kernel of unadjusted MCMC algorithms with proposal density q(·|·) is then P (x, dy) = q(y|x)dy,
(38)
whose invariant measure is not necessarily equal to µ in general. For Metropolis adjusted MCMC algorithms with proposal density q(·|·) targeting a probability density µ, the 14
transition kernel is P (x, dy) = q(y|x)α(y|x)dy + δx (dy)
Z
1 − α(y ′ |x) q(y ′ |x)dy ′ ,
(39)
µ(y)q(x|y) where δx is the Dirac measure centered at x and α(y|x) = min{1; µ(x)q(y|x) }. This definition of α implies detailed balance µ(x)P (x, dy)dx = µ(y)P (y, dx)dy, which ensures µ∞ = µ; see for instance [63]. To analyze the convergence of µk towards µ∞ , we consider the linear operator f 7→ P f on L2µ∞ defined by Z P f (x) = f (y)P (x, dy). y∈Rd
This operator has many interesting properties that are detailed in [44, Chapter 12]. First, together Rwith the property µ∞ = µ∞ P , Jensen’s inequality allows us to write R R (P f )2 dµ∞ ≤ P (f 2 )dµ∞ = f 2 dµ∞ so that f 7→ P f is a contraction in L2µ∞ . Second, R detailed balance ensures the operator is symmetric in L2µ∞ , meaning g(P f )dµ∞ = R (P g)f dµ∞ . Third, since P (·|x) is a probability measure for all x, we have P 1 = 1. the density of µk with respect to µ∞ , the detailed Finally, denoting by fk (x) = µµ∞k (x) (x) balance (again) permits us to write Z Z P (x, dy) P (y, dx)dy µk+1 (dx) µk (y) P fk (x) = µk (y) = = = fk+1 (x). µ∞ (y) µ∞ (x)dx µ∞ (x)dx y∈Rd y∈Rd To analyze the convergence of µk towards µ∞ , we consider the chi-square divergence χ2 (µk+1 ||µ∞ ) = Varµ (µk+1 /µ∞ ) of µk+1 from µ∞ . Because P (fk − 1) = fk+1 − 1, we can write Z 2 χ (µk+1 ||µ∞ ) = (P (fk − 1))2 dµ∞ ! R (P (fk − 1))2 dµ∞ χ2 (µk ||µ∞ ) = R (fk − 1)2 dµ∞ R (P f )2 dµ∞ χ2 (µk ||µ∞ ) R ≤ sup 2 f dµ∞ f ∈L2µ∞ s.t. R
f dµ∞ =0
= (1 − Gap(P ))2 χ2 (µk ||µ∞ ),
hence the geometric convergence χ2 (µk ||µ∞ ) ≤ (1 − Gap(P ))2k χ2 (µ0 ||µ∞ ). Here, Gap(P ) denotes the spectral gap of P which is defined by Gap(P ) = 1 − λ2 where λ2 is the 2nd largest eigenvalue of P seen as an operator in L2µ∞ (the largest eigenvalue begin λ1 = 1). An equivalent expression for the spectral gap is E(f, f ) Gap(P ) = inf2 , f ∈Lµ∞ Varµ∞ (f ) where E(f, f ) is the Dirichlet form associated with P defined by Z E(f, f ) := f 2 − f (P f )dµ∞ Z 1 = f 2 − 2f (P f ) + P (f 2 )dµ∞ 2 Z Z 2 1 = f (x) − f (y) P (x, dy)µ∞ (x)dx. 2 15
(40)
We show in Section 5.1 how to numerically estimate Gap(P ) in dimension d = 1. Remark 4.1 (Spectral Gap and Poincaré constant). The above expression (40) permits to link the spectral gap with the Poincaré constant. Indeed, a Taylor expansion on f permits to write (f (x) − f (y))2 = (∇f (x)⊤ (y − x))2 + O(∥x − y∥3 ) so that Z 1 ∇f (x)⊤ Σ(x)∇f (x)µ∞ (x)dx + O(δ), E(f, f ) = 2 R RR where Σ(x) = (y − x)(y − x)⊤ P (x, dy) and δ = ∥x − y∥3 P (x, dy)µ∞ (x)dx. For instance, for the unadjusted algorithm (38) with proposal density ( m∆t (x) = x + ∆tb(x) + O(∆t3/2 ) q∆t (·|x) = N (m∆t (x), Σ∆t (x)) where (41) Σ∆t (x) = 2∆tW (x) + O(∆t3/2 ), with b : Rd → Rd and W : Rd → S+d two arbitrary functions, we have Σ(x) = 2∆tW (x) + d
O(∆t3/2 ) and δ = O(∆t3/2 ). Thus, if µ∞ → µ as ∆t → 0 (which is the case when b(x) = divW (x) − W (x)∇V (x), see e.g. [47, Theorem 5.1]), we deduce E(f, f ) = ∆tEµ [∥∇f ∥2W ] + O(∆t3/2 ).
(42)
A similar expression can be obtain with the Metropolis adjusted algorithm (39) with the same proposal density q∆t (·|·). Let us emphasis that the constants hidden in the term O(∆t3/2 ) depend on f : if we had the stronger statement E(f, f ) = ∆tEµ [∥∇f ∥2W ] + ∆t3/2 R(f ) with R(f ) ≤ C(1 + Eµ [∥∇f ∥2W ]) for all f ∈ L2µ , we would be able to conclude that ∆t + O(∆t3/2 ), (43) Gap(P ) = C(µ, W ) where C(µ, W ) is the Poincaré constant such that Varµ (f ) ≤ C(µ, W )Eµ [∥∇f ∥2W ] for all smooth function f . A precise derivation of (43) is beyond the scope of this paper and is left for future work. Remark 4.2. The spectral gap relates to the integrated autocorrelation time (IACT) as Gap(P ) =
2 1 + supf ∈L2µ∞ IACT(f )
where IACT(f ) = 1 + 2
,
∞ X Cov(f (x0 ), f (xk )) k=1
Var(f (x0 ))
,
see for instance [27]. In principle, a sample-based estimator of the IACT allows one to \ ) for any (fixed) function f ∈ L2µ . However, estimating approximate IACT(f ) ≈ IACT(f ∞ supf ∈L2µ∞ IACT(f ) is more challenging, since one must estimate IACT(f ) simultaneously for several functions f . For instance, consider the finite-dimensional subspace VN = span{1Ωi , . . . , 1ΩN } of piecewise constant functions where Ω1 , . . . , ΩN is a partition of Rd . Accurate estimation of supf ∈VN IACT(f ) requires a sufficient number of samples in each subdomain Ωi , which becomes impractical when N is large. We show in Section 5.1 an alternative approach for computing Gap(P ) directly in dimension d = 1.
16
Acronym IS RW MALA (ULA) pCN GRiLS (U-GRiLS)
iid
Candidate x† ∼ q(·|xk ) (with ξk ∼ N (0, Σ)) x† = m + ξ√ k x† = xk + 2∆tξk √ x† = xk + ∆tΣ∇ ln µ(xk√) + 2∆tξk x† = θxk + (1 − θ)m + 1 − θ2 ξk , where θ = exp(−∆t) q x† = θk xk + (1 − θk )m +
1 − θk2 ξk ,
ν(xk ) where θk = exp(−∆t µ(x ) and ν = N (m, Σ) k)
Parameters Σ ≻ 0, m ∈ Rd Σ ≻ 0, ∆t > 0 Σ ≻ 0, ∆t > 0 Σ ≻ 0, m ∈ Rd , ∆t > 0 Σ ≻ 0, m ∈ Rd , ∆t > 0
Table 2: Metropolis-adjusted MCMC algorithms: Independence Sampler (IS), Random walk (RW), Metropolis Adjusted Langevin Algorithm (MALA), preconditionned Crank-Nicolson (pCN) and Gradient-free Riemannian Langevin Sampler (GRiLS). Unadjusted algorithms: Unadjusted Langevin Algorithm (ULA) and unadjusted GRiLS (U-GRiLS).
5
Numerical Experiments in Dimension One
We compare numerically the the MCMC algorithms listed in Table 2 on two target measures in dimension one2 . We consider a Gaussian mixture µh and a smoothed piecewise constant density µε defined respectively by (x − h/2)2 (x + h/2)2 µh (x) ∝ 1[−10,10] (x) exp − + exp − 2 2 4 X µε (x) ∝ 1[−2,2] (x) αi (σε (x − ri ) − σε (x − ri+1 )) . i=1
1 is a sigmoid function (with the convention σε=0 (t) = 1+ 12 sign(t)), Here σε (t) = 1+exp(−t/ε) and α = (0.5, 4, 1, 3) and r = (−2, 1, 0, 1, 2). As illustrated on Figure 2, the parameters h and ε control the difficulty of the problem: large values of h ≫ 1 enforce stronger multimodality in the Gaussian mixture and small values of ε ≪ 1 make µε closer to a piecewise-constant measure. In these regimes, standard Langevin-based algorithms are expected to perform poorly. In these experiments, the mean m and the covariance Σ are computed analytically (see Section 6 for the use of Algorithm 2 to estimate both m and Σ from interacting particles).
5.1
Ulam’s Method
In order to numerically estimate µ∞ (for unadjusted algorithms), the expected acceptance rate (for adjusted algorithms) and Gap(P ), we employ a piecewise constant discretization scheme. Let [a, b] = supp(µ) be the support of the target measure and let ci = a + (b − a) Ni−1 for 1 ≤ i ≤ N . In all our experiments, we take N = 2000. We consider the −1 ci−1 +ci c1 +ci+1 partition R = ∪N , 2 ] with the convention c0 = −∞ and i=1 Ωi where Ωi = [ 2 2
The implementation to reproduce the numerical results can be found at: https://gitlab.inria. fr/ozahm/grils_dim1.
17
=1 = 0.1 = 0.01 =0
h=1 h=2 h=4 h=8
3
10.0 7.5 5.0 2.5 0.0 2.5 5.0 7.5 10.0 (a) Mixture of two Gaussian µh
2
1
0
1
2
3
(b) Smoothed piecewise constant density µε
Figure 2: Two target densities in dimension d = 1. Left: Mixture of two Gaussian µh = 12 N (−h/2, 1) + 12 N (+h/2, 1) cN +1 = +∞. We then assemble the matrix Z N Qij = q(y|ci )dy, Ωj
by using the closed form expression of the CDF of the Gaussian density q(·|ci ). By construction, QN = (QN ij ) is a stochastic matrix (all its entries are nonnegative and each row sums to one) which approximates the unadjusted algorithm. By the FrobeniusPerron theorem, there is a nonnegative vector v ∈ RN whose components sum to one, such that v = (QN )⊤ v, and the corresponding density function µN ∞ (x) =
N X vi
|Ωi | i=1
1Ωi (x),
is a piecewise constant approximation to µ∞ called the Ulam approximation [65]. Although this method is simple, its convergence is usually slow since it only uses piecewise −1 ln(N )), see [22]. In addition, the constant functions, typically ∥µ∞ − µN ∞ ∥L1 = O(N N N spectral gap Gap(Q ) = 1−λ2 can be computed by finding the second largest eigenvalue N ⊤ N λN 2 of (Q ) w = λ2 w. The stochastic PN µ(ci ) matrix representing the Metropolis adjusted algorithm targeting N µ (x) = i=1 |Ωi | 1Ωi (x) is assembled via ( N −1 min{µ(ci )QN i ̸= j, ij , µ(cj )Qji }µ(ci ) PijN = P N 1 − k̸=i Pik i = j. The spectral gap Gap(P N ) = 1 − λN 2 is also be computed by finding the second largest N N ⊤ N eigenvalue λ2 of (P ) w = λ2 w, and the expected acceptance rate is given by N
E[α ] =
N X i,j=1
5.2
N min{µ(ci )QN ij , µ(cj )Qji }.
Unadjusted Algorithms
We compare here ULA and U-GRiLS. Because both algorithms are based on a timeN discretization of a Langevin dynamic, their stationary measure µN ∞ converges to µ (x) = 18
h=1
100
h=2
h=4
h=8
10 1
Gap(Q N)
t=1 10 2 10 3 ULA U-GRiLS
10 4
2
10 5 10 7
4
10 2
100
Gap(Q N)
10 1
t = 10 4
( || )
10 1
10 7
10 4
10 1
( || N )
2
= 0.1
= 0.01
10 7
( || N )
ULA U-GRiLS 2( || )
10 4 2
10 1
( || N )
= 0.001
10 7
10 4 2
10 1
( || N )
= 0.0
t=1 t = 10 4
10 2 10 3 10 4 4 10 10 3 10 2 10 1 2( || N )
100 10 4 10 3 10 2 10 1 2( || N )
100 10 4 10 3 10 2 10 1 2( || N )
100 10 4 10 3 10 2 10 1 2( || N )
100
Figure 3: Unadjusted algorithms: spectral gap V.S. bias computed for different values of ∆t ∈ [10−4 , 100 ]. Top row: Gaussian mixture µh for h ∈ {1, 2, 4, 8}. Bottom row: smoothed piecewise constant density µε for ε ∈ {10−1 , 10−2 , 10−3 , 0}. √ = O( ∆t), see [47, Theorem 5.1]. In addition, their spectral gap decreases with ∆t as Gap(QN ) = O(∆t) according to (43). Thus we expect 2 Gap(QN ) = O χ2 (µN ||µN ) . (44) ∞ PN
µ(ci ) 2 N N i=1 |Ωi | 1Ωi (x) as χ (µ ||µ∞ )
Figure 3 represents the spectral gap as a function of the bias, both quantities being computed for a range of ∆t ∈ [10−4 , 100 ]. Note that relation (44) is observed for UGRiLS on µh for h = 8 and on µε with ε ∈ {0.01, 0.001, 0}. In addition, in these plots, a larger spectral gap Gap(QN ) indicates faster convergence toward the stationary measure µN ∞ , assuming a comparable bias. Overall, we observe that GRiLS often outperforms ULA. Note that the worth performance of ULA is for the smoothed piecewise constant density µε=0 for which the bias do not decrease. This is not surprinsing, since ∇ ln µ(ε=0) (x) = 0 almost surely. In contrast, the stationary measure of U-GRiLS is consistently approaching µ(ε=0) as ∆t → 0, as shown on Figure 4.
5.3
Metropolis Adjusted Algorithms
Figure 5 shows the performance of the Metropolis adjusted algorithms by plotting Gap(P N ) as a function of the rejection rate 1 − E[αN ]. Both quantities are being computed for a range of ∆t ∈ [10−4 , 100 ]. For all √ considered algorithms, the rejection rate is known to behave like 1 − E[αN ] = O( ∆t) with the expection of MALA for which
19
N, N, N, N,
t = 1e t = 1e t = 1e t = 1e
01 02 03 04
∆t 10−1 10−2 10−3 10−4
( = 0)
3
2
1
0
1
χ2 (µε=0 ||µN ∞) −1 1.44 × 10 5.21 × 10−2 1.99 × 10−2 6.55 × 10−3
2
(a) Stationary measures µN ∞
(b) Bias
Figure 4: U-GRiLS targeting the piecewise constant density µ(ε=0) : stationary measures −1 −4 2 N µN ∞ and bias χ (µ||µ∞ ) for different ∆t ∈ [10 , 10 ]. 1 − E[αN ] = O(∆t3/2 ), see e.g. [55]. Thus we expect ( O (1 − E[αN ])2/3 for MALA, N Gap(Q ) = O (1 − E[αN ])2 for other algorithms. This asymptotic behaviors are clearly observed for µh for h ∈ {1, 2, 4}. Is it worth to note that, albeit being based on the discretization of a Langevin dynamic, GRiLS do not heritate from the favourable rejection rate 1 − E[αN ] = O(∆t3/2 ) and has a spectral gap behaving as Gap(QN ) ∼ (1 − E[αN ])2 . Among all algorithms, the independence sampler (IS) consistently yields the smallest spectral gap, making it the most effective Metropolis-adjusted algorithm for the onedimensional problems we consider here. This indicates that, for the Gaussian mixture µh (top row), the most efficient mechanism to jump across the modes is to draw a independent sample. Note also that the expected rejection rate 1 − E[αN ] of IS increases rapidly with h, corroborating formula (30). This means that, albeit having the largest spectral gap, IS accepts less frequently when h ≫ 1. We also observe that GRiLS consistently performs better than pCN and, for large ∆t, both GRiLS and pCN recover IS, which is consistent with their construction (remember Remark 3.3). It is worth noting that, for h = 8, GRiLS achieves significantly largest spectral gap over a wider range of small acceptance probability. This indicates that GRiLS is able to explore efficiently each modes (small rejection rate) while still being able to jump across modes (large spectral gap). When considering MALA, we observe on both target densities that its performance deteriorates significantly for large h and small ε. Note also that, for the piecewise constant density µε=0 , MALA reduces to RW.
6
Higher-Dimensional Numerical Experiments
In this section, we demonstrate the performance of BE-GRiLS on multimodal densities3 . Subsection 6.1 considers two-dimensional benchmark densities, and Subsection 6.2 exhibits the performance with a ten-dimensional mixture density. 3
The implementation to reproduce the numerical results can be found at: https://github.com/ baptistar/GradientFreeRiemannianLangevinSampler.
20
h=1
100
h=2
h=4 t=1
10 1 Gap(P N)
h=8
10 2 MALA GRiLS RW pCN IS
10 3 10 4
t = 10 4
10 5 4 10 10 3 10 2 10 1 100 10 4 10 3 10 2 10 1 100 10 4 10 3 10 2 10 1 100 10 4 10 3 10 2 10 1 100 1 [ N] 1 [ N] 1 [ N] 1 [ N]
= 0.1
100
= 0.01
= 0.001 t=1
Gap(P N)
10 1
= 0.0
t = 10 4
10 2 MALA GRiLS RW pCN IS
10 3 10 4 10 5 2 10
1
10 1 [ N]
100 10 2
1
10 1 [ N]
100 10 2
1
10 1 [ N]
100 10 2
1
10 1 [ N]
100
Figure 5: Adjusted algorithms: spectral gap V.S. expected rejection rate for different values of ∆t ∈ [10−4 , 100 ]. Top row: Gaussian mixture µh for h ∈ {1, 2, 4, 8}. Bottom row: smoothed piecewise constant density µε for ε ∈ {10−1 , 10−2 , 10−3 , 0}.
6.1
Benchmark Two-Dimensional Distributions
First, we assess the performance of BE-GRiLS on benchmark non-Gaussian densities µ ∈ P(R2 ) including a mixture of Gaussians, the two moons and the two rings densities defined by: 3 1X ∥x − mk ∥2 µ(x) ∝ exp − 3 k=1 σ2 (∥x∥ − 1)2 (x1 − 2)2 4x1 /a − µ(x) ∝ (1 + exp ) exp − b 2a 2 2 1 (∥x∥ − 1) (∥x∥ − 2) µ(x) ∝ − − , ∥x∥ 2σR2 2σR2 √ √ respectively. For the mixture of Gaussians example, we set m1 = [− 3, 1], m2 = [ 3, 1] and m3 = [0, 2] so that the modes are centered at the vertices of an equilateral triangle, and σ 2 = 0.1. For the two moons example, we set a = 0.08, b = 0.08. For the two rings example, we set σR = 0.1. We evaluate the block ensemble GRiLS in comparison to two gradient-free methods: adaptive Metropolis [30], which updates the proposal covariance adaptively based on the history of the Markov chain, and the Affine Independent Ensemble Sampler (AIES, [29]) which, similarly to BE-GRiLS, is based on an ensemble of particles. We also compare to the gold-standard Metropolis Adjusted Langevin Algorithm (MALA), which has access to gradients from the target density and whose covariance is also adapted based on the history of the chain. For AIES and BE-GRiLS, we initialize the samplers using N = 10 21
Adaptive Metropolis
Adaptive MALA
AIES
Block Ensemble GRiLS
3 2
x2
1 0 −1 −2 0 x1
−2
2
Adaptive Metropolis
−1
−2
0
2
−1
x2
0 −1
x2
3
−2 10000 20000 Iteration
30000
10000 20000 Iteration
30000
2
Block Ensemble GRiLS 0.0 −2.5 2
−2 0
0 x1
2.5
0.0 −2.5
0
−2
2
AIES
−3
1
0 x1
2.5
−1
x1
x1
0
x1
−2
2
Adaptive MALA
1
x2
0 x1
x1
−2
x2
−3
0
10000 20000 Iteration
30000
0 −2
0
10000 20000 Iteration
30000
Figure 6: Samples and trace plots for the mixture of Gaussians density across four algorithms. The BE-GRiLS algorithm shows fast mixing behavior across and coverage of all three modes. particles drawn from a standard Gaussian density, while the Adaptive Metropolis and Adaptive MALA algorithms are initialized at the MAP of the target density with an identity covariance matrix for the Gaussian proposal. We run BE-GRiLS with P = 2 blocks of equal size. Each method is run for 21,000 steps and we consider 1,000 steps as burn-in. We select ∆t = 0.3 in these experiments to minimize the probability of missing modes due to initialization of the chain. Figures 6-8 plot the sample history relative to the contours of each target density and the trace plots for the two coordinates. Overall, we observe that only BE-GRiLS consistently captures all modes of the target density. Moreover, the trace plots demonstrate fast mixing between the modes as a result of the preconditioning. We note that while adaptive Metropolis can adapt the covariance structure to capture some of these two-dimensional densities with a sufficiently large step-size, the next subsection will demonstrate that it can generally miss modes in higher-dimensional problems.
6.2
Mixture of Gaussians in Ten Dimension
Lastly, we assess the performance of BE-GRiLS on approximating the first two moments of a d = 10 dimensional Gaussian mixture target density with three equal-weight components. The first two coordinates match the target density of the Gaussian mixture in Subsection 6.1, while the remaining d − 2 coordinates follow the law of independent standard normal variables. We compare the algorithm’s performance to adaptive random-walk Metropolis (AM), MALA with an adaptive covariance, and the affine-invariant ensemble sampler (AIES). For each algorithm we run K = 105 total iterations with the first 103 steps discarded as burn-in. We initialize the single-chain methods (AM and MALA) from the target’s MAP and the ensemble methods (AIES and BE-GRiLS) from N = 40 walkers drawn from N (µMAP , Σπ ), where Σπ is the true global covariance in order to maximize the success of 22
Adaptive Metropolis
Adaptive MALA
AIES
Block Ensemble GRiLS
1
x2
0
−1
−1
0 x1
1
−2
2
0 x1
1
−2
2
10000 15000 Iteration
0 x1
1
2
Block Ensemble GRiLS
0 −2
1
1
0
0
x2
x2
1.0
1
0
5000
10000 15000 Iteration
20000
0 −1
−1
−1
20000
−1
2
0.5
x2 5000
−2
2
−2 1
−1
1
1.5
−1
−1.5
0
0 x1
AIES x1
−1.0
0
−1
Adaptive MALA x1
x1
Adaptive Metropolis
−1
x1
−2
x2
−2
0
5000
10000 15000 Iteration
20000
0
5000
10000 15000 Iteration
20000
Figure 7: Samples and trace plots for the two moons density across four algorithms. Adaptive Metropolis
Adaptive MALA
AIES
Block Ensemble GRiLS
2
x2
1 0 −1 −2 −2
2
0 −2
0
10000 20000 Iteration
30000
0 x1
−2
2
AIES
0
0
0
−2
−2
2
2
2
−2
0
10000 20000 Iteration
30000
0 −2
2
2
−2
0
0 x1
Block Ensemble GRiLS
2 x1
x1 x2
x1 x2
2
Adaptive MALA
2
0 −2
−2
2
x2
Adaptive Metropolis 2
0 x1
x1
0 x1
x2
−2
0
10000 20000 Iteration
30000
0 −2
0
10000 20000 Iteration
30000
Figure 8: Samples and trace plots for the two rings density across four algorithms. all methods. The step sizes ∆t (or the AIES stretch scale factor) are independently tuned for each algorithm by a pilot run of length 103 steps by maximizing the expected squared jump distance per dimension, ESJD/d = d1 E∥xt+1 − xt ∥22 . This criterion avoids select an overly small step, which would miss certain modes of the multimodal density [53]. For each iteration k, we compute the errors between the running empirical mean and covariance (m̂k , Σ̂k ) and the true mean and covariance (m, Σ) of the target density. The errors normalized by an appropriate dimension scaling are given by: ∥Σ̂k − Σ∥F . d
∥m̂k − m∥2 √ , d
Figure 9 plots the mean and covariance error as a function of increasing iterations for the four algorithms. The results are averaged over 5 independent repetitions with the mean 23
100
4 × 10−1
km̂k − mk2 /
√
kΣ̂k − ΣkF / d
d
3 × 10−1
2 × 10−1
10−1
Adaptive Metropolis Adaptive Covariance MALA AIES BE-GRiLS 101
102
103
10−1
104
105
101
Number of steps, k
102
103
104
105
Number of steps, k
Figure 9: Convergence of the empirical mean (left) and covariance (right) for the d = 10 dimensonal Gaussian mixture across four algorithms. reported along with one standard deviation in the figure. Across the tested settings, BE-GRiLS generally attains the lowest normalized mean and covariance errors at a given sample budget than the single-chain AM and MALA baselines or the ensemble method AIES. As in Figure 6, we observe fast mixing of BE-GRiLS across the three modes of the density, unlike the alternative methods which are more likely to become trapped within a single mixture component and hence result in lower global accuracy.
7
Discussions and Conclusions
In this work, we have introduced GRiLS, a gradient-free sampler based on the time discretization of a Riemannian Langevin dynamics, and demonstrated its connections to several established MCMC methods. Specifically, GRiLS unifies a broad family of proposals: it recovers an independence sampler (IS) in the limit ∆t → ∞, it simplifies to preconditioned Crank-Nicolson (pCN) scheme when the parameter θk = θ is held constant, and shares the structural character of MALA/ULA through its Langevin foundation as ∆t → 0. Moreover, the ensemble variant BE-GRiLS has a close resemblance to Consensus-Based Sampling (CBS), with each proposal combining a contraction toward the ensemble mean with a Gaussian perturbation. This unifying perspective suggests that GRiLS provides a principled and flexible framework that interpolates between and generalizes several state-of-the-art MCMC proposals. Several directions for future work remain open. An important question in practice concerns the systematic selection of the step size ∆t: as in gradient-based optimization, this choice is critical for algorithmic performance, and principled heuristics or adaptive strategies remain to be developed. On the theoretical side, it is natural to ask whether the classical optimal scaling analysis for MALA as in [55] can be adapted to GRiLS, which, despite being gradient-free, is fundamentally a Langevin sampler. Also, following [2, 46], a convergence analysis via Cheeger inequalities, offers another promising theoretical avenue. Beyond the Gaussian setting, generalizing GRiLS to non-Gaussian reference measures ν would broaden its applicability, as outlined in Remark 2.2. Improving the accuracy of the time integration through higher-order quadrature schemes represents a further refinement of the method. In addition, as done in CBS [8], a rigorous analysis of the mean-field 24
limit of BE-GRiLS as N → ∞ would provide a deeper understanding of the ensemble dynamics in a non-stationary regime, that is, when the target mean and covariance are poorly estimated by the ensemble. Finally, a critical question that this paper does not fully resolve concerns the robustness of GRiLS in high dimensions d ≫ 1. As shown in Equation (28), the acceptance rate of the Metropolis-corrected GRiLS algorithm degenerates as the dimension grows, revealing an intrinsic limitation of proposals with location-dependent covariance. A promising remedy, to be pursued in future work, is to confine the spatial variations of the proposal covariance to a low-dimensional subspace, in the spirit of [43], thereby mitigating this adverse effect of dimensionality. Acknowledgments. The authors appreciate discussion of some of the results with Björn Sprungk and Simon Barthelmé and Bryan Toto.
A
Proof for (19) and (20)
Inequality (19) can be deduce from (20) by letting d = 1, N = 2, mi = ±h. To proof (20), we use the following proposition. P Proposition 4. Let µ(x) = N i=1 ωi µi (x) be a d-dimensional Gaussian mixture, where µi = N (mi , Σi ), Σi ≻ 0, for some weights ωi ≥ 0 that sum to one. Its mean and covariance are given by m=
N X i=1
ωi mi ,
and
Covµ =
N X i=1
ωi Σi + (mi − m)(mi − m)⊤ .
For any matrix Σ such that Σ ≻ Σi for all i ≤ N , the Gaussian measure ν = N (m, Σ) satisfies s ! N X ∥mi − m∥2(Σ−Σi )−1 |Σ| µ(x) ≤ ωi exp . (45) sup |Σi | 2 x∈Rd ν(x) i=1 −1 −1 Proof. Let Ai = Σ−1 and ai = (Σ−1 i −Σ i mi − Σ m). Because Σ ≻ Σi , we have Ai ≻ 0 and then s ! N ∥x − mi ∥2Σ−1 ∥x − m∥2 −1 µ(x) X |Σ| i Σ = ωi exp − + ν(x) |Σ | 2 2 i i=1 s ! N −1 ∥mi ∥2Σ−1 ∥m∥2 −1 X |Σ| x⊤ (Σ−1 i − Σ )x −1 ⊤ −1 i Σ ωi exp − = + x (Σi mi − Σ m) − + |Σ | 2 2 2 i i=1 s ! N 2 −1 2 ∥mi ∥2Σ−1 ∥m∥2 −1 X ∥x − A−1 a ∥ ∥A a ∥ |Σ| i i Ai Ai i i Σ = ωi exp − + − + |Σ | 2 2 2 2 i i=1 s ! N 2 2 X ∥m − m∥ −1 ∥x − A−1 a ∥ i |Σ| i Ai (Σ−Σi ) i = ωi exp − + . |Σ | 2 2 i i=1
25
For the last step, we used −1 −1 −1 2 −1 ⊤ −1 −1 −1 ∥A−1 i ai ∥Ai = (Σi mi − Σ m) (Σi − Σ ) (Σi mi − Σ m) −1 −1 ⊤ −1 = m⊤ i (Σi − Σi Σ Σi ) mi − 2mi (Σ − Σi ) m
−1 + m⊤ (ΣΣ−1 i Σ − Σ) m
−1 = ∥mi − m∥2(Σ−Σi )−1 − m⊤ (Σ − Σi )−1 m − m⊤ i (Σ − Σi ) mi
−1 −1 −1 ⊤ −1 + m⊤ i (Σi − Σi Σ Σi ) mi + m (ΣΣi Σ − Σ) m −1 = ∥mi − m∥2(Σ−Σi )−1 − m⊤ Σ−1 m + m⊤ i Σi mi ,
where we employed the Woodbury formula −1 (Σi − Σi Σ−1 Σi )−1 = Σ−1 i + (Σ − Σi ) −1 (ΣΣ−1 = (Σ − Σi )−1 − Σ−1 . i Σ − Σ)
This yields N µ(x) X ≤ ωi ν(x) i=1
s
|Σ| exp |Σi |
∥mi − m∥2(Σ−Σi )−1
!
2
,
which is (45) We now prove Inequality (20). By Propositions 1 and 4, the Gaussian mixture N
µN h (x) =
1 X µi (x), N i=1
µi = N (mi , Id ),
(46)
satisfies µ(x) x∈Rd ν(x) s ! N (45) X ∥mi − m∥2(Σ−Σi )−1 |Σ| ≤ ωi exp |Σ | 2 i i=1 s ! N ∥mi − m∥2(Covµ −Id )−1 | Covµ | (46) 1 X exp = N i=1 |Id | 2 (17)
C(µ, WµNh ) ≤ sup
Next, because ∥mi − m∥ ≤ h for all i, we have Covµ = Id + N1 (1 + h2 )Id so that | Covµ | ≤ (1 + h2 )d . Furthermore, we have
PN
⊤ i=1 (mi − m)(mi − m) ⪯
∥mi − m∥2(Covµ −Id )−1 = N ∥mi − m∥2(PN (mi −m)(mi −m)⊤ )−1 ≤ N, i=1
for all i ≤ N , so that the above inequality becomes C(µ, WµNh ) ≤ (1 + h2 )d/2 exp(N/2) which is (20).
26
B
Proof for inequality (18)
Consider the Gaussian mixture µ = 12 N (−h, 1) + 21 N (+h, 1) in dimension d = 1. In order to show that 2 eh /3 C(µ, Covµ ) ≥ , 1 + h2 holds for any h ≥ 1/2, one simply test the Poincaré inequality (8) with a smooth function Rx 2 /2) √ f that is approximately ±1 on each of the two modes of µ. Choose f (x) = 2 0 exp(−t dt. 2π ′ 2 2 Because f (x) = 2/π exp(−x ) is even, we have Z exp(−(x − h)2 /2) ′2 √ dx Eµ [f ] = f ′ (x)2 2π Z 2 (x − h)2 2 = √ exp −x − dx 2 π 2π Z (x − 13 h)2 h2 h2 2 = √ exp − + − dx 2/3 6 2 π 2π h2 2 = √ e− 3 . π 3 Because Covµ = 1 + h2 we deduce √ h2 Varµ (f ) Varµ (f )π 3 e 3 C(µ, Covµ ) ≥ = . E[(f ′ )2 (1 + h2 )] 2 1 + h2 (8)
Furthermore, since f is odd we have Eµ [f ] = 0 and Z Z 2 exp(−x2 /2) 2 exp(−(x − h) /2) √ √ dx = f (x + h)2 dx. Varµ (f ) = f (x) 2π 2π Thus, the function h 7→ Varµ (f ) is monotonically increasing (because f is) and one can 2 check numerically that Varµ (f ) > π√ as soon as h ≥ 1/2. This gives the result. 3
References [1] G. Allaire, Shape Optimization by the Homogenization Method, Springer, 2002. [2] C. Andrieu, A. Lee, S. Power, and A. Q. Wang, Explicit convergence bounds for Metropolis Markov chains: Isoperimetry, spectral gaps and profiles, The Annals of Applied Probability, 34 (2024), pp. 4022–4071. [3] C. Andrieu and J. Thoms, A tutorial on adaptive MCMC, Statistics and computing, 18 (2008), pp. 343–373. [4] D. Bakry, I. Gentil, and M. Ledoux, Analysis and geometry of Markov diffusion operators, vol. 348, Springer Science & Business Media, 2013. [5] A. Beskos, M. Girolami, S. Lan, P. E. Farrell, and A. M. Stuart, Geometric MCMC for infinite-dimensional inverse problems, Journal of Computational Physics, 335 (2017), pp. 327–351. 27
[6] M. Betancourt, A general metric for Riemannian manifold Hamiltonian Monte Carlo, in International Conference on Geometric Science of Information, Springer, 2013, pp. 327–334. [7] L. Cao, T. O’Leary-Roseberry, and O. Ghattas, Derivative-informed neural operator acceleration of geometric MCMC for infinite-dimensional Bayesian inverse problems, Journal of Machine Learning Research, 26 (2025), pp. 1–68. [8] J. A. Carrillo, F. Hoffmann, A. M. Stuart, and U. Vaes, Consensus-based sampling, Studies in Applied Mathematics, 148 (2022), pp. 1069–1140. [9] D. Chakraborty, R. Harris, R. Klein, G. Olicón-Méndez, S. Reich, and C. Schillings, Affine invariant langevin dynamics for rare-event sampling, (2025). [10] S. Chewi, Log-concave sampling, Book draft available at https://chewisinho. github. io, 9 (2023), pp. 17–18. [11] K. L. Chung, Lectures from Markov processes to Brownian motion, vol. 249, Springer Science & Business Media, 2013. [12] J. Comer, J. C. Gumbart, J. Hénin, T. Lelièvre, A. Pohorille, and C. Chipot, The adaptive biasing force method: Everything you always wanted to know but were afraid to ask, The Journal of Physical Chemistry B, 119 (2015), pp. 1129–1151. [13] D. Cordero-Erausquin and B. Klartag, Moment measures, Journal of Functional Analysis, 268 (2015), pp. 3834–3866. [14] S. L. Cotter, G. O. Roberts, A. M. Stuart, and D. White, MCMC methods for functions: modifying old algorithms to make them faster, Statistical Science, (2013), pp. 424–446. [15] T. Cui, S. Dolgov, and O. Zahm, Scalable conditional deep inverse Rosenblatt transports using tensor trains and gradient-based dimension reduction, Journal of Computational Physics, 485 (2023), p. 112103. [16] T. Cui, X. Tong, and O. Zahm, Optimal Riemannian metric for Poincaré inequalities and how to ideally precondition Langevin dynamics, arXiv preprint arXiv:2404.02554, (2024). [17] A. Durmus and E. Moulines, Nonasymptotic convergence analysis for the unadjusted Langevin algorithm, The Annals of Applied Probability, 27 (2017), pp. 1551– 1587. [18] R. Dwivedi, Y. Chen, M. J. Wainwright, and B. Yu, Log-concave sampling: Metropolis-Hastings algorithms are fast, Journal of Machine Learning Research, 20 (2019), pp. 1–42. [19] B. Engquist, K. Ren, and Y. Yang, Adaptive state-dependent diffusion for derivative-free optimization, Communications on Applied Mathematics and Computation, 6 (2024), pp. 1241–1269.
28
[20]
, Sampling with adaptive variance for multimodal distributions, arXiv preprint arXiv:2411.15220, (2024).
[21] M. Ernst, G. Reinert, and Y. Swan, First-order covariance inequalities via Stein’s method, Bernoulli, 26 (2020), pp. 2051 – 2081. [22] G. Froyland, Approximating physical invariant measures of mixing dynamical systems in higher dimensions, Nonlinear Analysis: Theory, Methods & Applications, 32 (1998), pp. 831–860. [23] A. Garbuno-Inigo, F. Hoffmann, W. Li, and A. M. Stuart, Interacting Langevin diffusions: Gradient structure and ensemble Kalman sampler, SIAM Journal on Applied Dynamical Systems, 19 (2020), pp. 412–441. [24] A. Garbuno-Inigo, N. Nüsken, and S. Reich, Affine invariant interacting langevin dynamics for bayesian inference, SIAM Journal on Applied Dynamical Systems, 19 (2020), pp. 1633–1658. [25] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin, Bayesian data analysis, Chapman and Hall/CRC, 1995. [26] G. Germain and Y. Swan, A note on one-dimensional Poincaré inequalities by Stein-type integration, Bernoulli, 29 (2023), pp. 1714–1740. [27] C. J. Geyer, Practical Markov chain Monte Carlo, Statistical science, (1992), pp. 473–483. [28] M. Girolami and B. Calderhead, Riemann manifold Langevin and Hamiltonian Monte Carlo methods, Journal of the Royal Statistical Society Series B: Statistical Methodology, 73 (2011), pp. 123–214. [29] J. Goodman and J. Weare, Ensemble samplers with affine invariance, Communications in applied mathematics and computational science, 5 (2010), pp. 65–80. [30] H. Haario, E. Saksman, and J. Tamminen, An adaptive Metropolis algorithm, Bernoulli, (2001), pp. 223–242. [31] D. Heredia, On weighted poincaré inequalities for multivariate liouville distributions–application to global sensitivity analysis, arXiv preprint arXiv:2605.30979, (2026). [32] D. Heredia, A. Joulin, and O. Roustant, On one dimensional weighted Poincaré inequalities for Global Sensitivity Analysis, Journal of Mathematical Analysis and Applications, (2025), p. 129992. [33] M. Hird, F. Maire, and J. Negrea, A non-asymptotic analysis for learning and applying a preconditioner in MCMC, arXiv preprint arXiv:2602.10714, (2026). [34] R. A. Holley, S. Kusuoka, and D. W. Stroock, Asymptotics of the spectral gap with applications to the theory of simulated annealing, J. Funct. Anal, 83 (1989), pp. 333–347.
29
[35] Y.-P. Hsieh, A. Kavis, P. Rolland, and V. Cevher, Mirrored Langevin dynamics, Advances in Neural Information Processing Systems, 31 (2018). [36] J. Kent, Time-reversible diffusions, Advances in Applied Probability, 10 (1978), pp. 819–835. [37] T. S. Kleppe, Log-density gradient covariance and automatic metric tensors for Riemann manifold Monte Carlo methods, Scandinavian Journal of Statistics, 51 (2024), pp. 1206–1229. [38] P. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Stochastic Modelling and Applied Probability, Springer Berlin Heidelberg, 2013. [39] R. Kutri and R. Scheichl, Fast-mixing markov chains without gradients, arXiv preprint arXiv:2606.27564, (2026). [40] J. Lamperti, Semi-stable Markov processes. I, Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 22 (1972), pp. 205–225. [41] B. Leimkuhler, C. Matthews, and J. Weare, Ensemble preconditioning for Markov chain Monte Carlo simulation, Statistics and Computing, 28 (2018), pp. 277– 290. [42] T. Lelièvre, G. Pavliotis, G. Robin, R. Santet, and G. Stoltz, Optimizing the diffusion coefficient of overdamped Langevin dynamics, Mathematics of Computation, (2025). [43] T. Lelièvre, R. Santet, and G. Stoltz, Improving sampling by modifying the effective diffusion, Journal of Computational Physics, (2025), p. 114313. [44] D. A. Levin and Y. Peres, Markov chains and mixing times, vol. 107, American Mathematical Soc., 2017. [45] M. T. Li, T. Cui, F. Li, Y. Marzouk, and O. Zahm, Sharp detection of lowdimensional structure in probability measures via dimensional logarithmic Sobolev inequalities, Information and Inference: A Journal of the IMA, 14 (2025), p. iaaf021. [46] S. Liu and X. T. Tong, Spectral gap of Metropolis algorithms for non-smooth distributions under isoperimetry, arXiv preprint arXiv:2601.21763, (2026). [47] J. C. Mattingly, A. M. Stuart, and M. V. Tretyakov, Convergence of numerical time-averaging and stationary measures via Poisson equations, SIAM Journal on Numerical Analysis, 48 (2010), pp. 552–577. [48] G. Menz and A. Schlichting, Poincaré and logarithmic Sobolev inequalities by decomposition of the energy landscape, The Annals of Probability, 42 (2014), p. 1809. [49] R. M. Neal, Sampling from multimodal distributions using tempered transitions, Statistics and computing, 6 (1996), pp. 353–366. [50] B. Øksendal, When is a stochastic integral a time change of a diffusion?, Journal of theoretical probability, 3 (1990), pp. 207–226. 30
[51]
, Stochastic differential equations, in Stochastic differential equations: an introduction with applications, Springer, 2003, pp. 38–50.
[52] M. D. Parno and Y. M. Marzouk, Transport map accelerated Markov chain Monte Carlo, SIAM/ASA Journal on Uncertainty Quantification, 6 (2018), pp. 645– 682. [53] C. Pasarica and A. Gelman, Adaptively scaling the metropolis algorithm using expected squared jumped distance, Statistica Sinica, (2010), pp. 343–364. [54] C. P. Robert, G. Casella, and G. Casella, Monte Carlo statistical methods, vol. 2, Springer, 1999. [55] G. O. Roberts and J. S. Rosenthal, Optimal scaling of discrete approximations to Langevin diffusions, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60 (1998), pp. 255–268. [56] G. O. Roberts and J. S. Rosenthal, Harris recurrence of Metropolis-withinGibbs and trans-dimensional Markov chains, The Annals of Applied Probability, 16 (2006), pp. 2123–2139. [57] O. Roustant, N. Lüthen, D. Heredia, and B. Sudret, Gradientenhanced global sensitivity analysis with Poincaré chaos expansions, arXiv preprint arXiv:2510.03056, (2025). [58] A. Saumard, Weighted Poincaré inequalities, concentration inequalities and tail bounds related to Stein kernels in dimension one, Bernoulli, 25 (2019), pp. 3978– 4006. [59] A. Schlichting, Poincaré and log-Sobolev inequalities for mixtures, Entropy, 21 (2019), p. 89. [60] S. Song, T. Zhou, L. Wang, S. Kucherenko, and Z. Lu, Derivative-based new upper bound of Sobol’ sensitivity measure, Reliability Engineering & System Safety, 187 (2019), pp. 142–148. [61] B. Sprungk, S. Weissmann, and J. Zech, Metropolis-adjusted interacting particle sampling, Statistics and Computing, 35 (2025), p. 64. [62] G. Stoltz, M. Rousset, et al., Free energy computations: A mathematical perspective, World Scientific, 2010. [63] L. Tierney, A note on Metropolis-Hastings kernels for general state spaces, Annals of applied probability, (1998), pp. 1–9. [64] M. Titsias, Optimal preconditioning and Fisher adaptive Langevin sampling, Advances in Neural Information Processing Systems, 36 (2023), pp. 29449–29460. [65] S. M. Ulam, A collection of mathematical problems, (No Title), (1960). [66] T. Xifara, C. Sherlock, S. Livingstone, S. Byrne, and M. Girolami, Langevin diffusions and the Metropolis-adjusted Langevin algorithm, Statistics & Probability Letters, 91 (2014), pp. 14–19. 31
[67] O. Zahm, T. Cui, K. Law, A. Spantini, and Y. Marzouk, Certified dimension reduction in nonlinear Bayesian inverse problems, Mathematics of Computation, 91 (2022), pp. 1789–1835. [68] K. S. Zhang, G. Peyré, J. Fadili, and M. Pereyra, Wasserstein control of mirror Langevin Monte Carlo, in Conference on learning theory, PMLR, 2020, pp. 3814–3841.
32