A proximal gradient algorithm for composite log-concave sampling Linghai Liu
Sinho Chewi
arXiv:2605.12461v1 [math.ST] 12 May 2026
Yale University
Abstract We propose an algorithm to sample from composite log-concave distributions over Rd , i.e., densities of the form π ∝ e−f −g , assuming access to gradient evaluations of f and a restricted Gaussian oracle (RGO) for g. The latter requirement means that we can easily sample from the 1 ∥y − x∥2 ), which is the sampling analogue of the proximal density RGOg,h,y (x) ∝ exp(−g(x) − 2h operator for g. If f + g is α-strongly and f is β-smooth, our sampler achieves ε error √ convex ̃ in total variation distance in O(κ d log4 (1/ε)) iterations where κ ∶= β/α, which matches prior state-of-the-art results for the case g = 0. We further extend our results to cases where (1) π is non-log-concave but satisfies a Poincaré or log-Sobolev inequality, and (2) f is non-smooth but Lipschitz.
1
1
Introduction
One of the cornerstones of convex optimization is the problem of composite optimization: min F (x) ∶= f (x) + g(x) ,
x∈Rd
where f is convex and smooth while g is convex but possibly non-smooth. The proximal gradient method approaches this problem by the iteration [PB14]: xk+1 = proxhg (xk − h∇f (xk )) ,
proxhg (x) ∶= arg min {g(y) + y∈Rd
1 ∥y − x∥2 } . 2h
The smooth part f admits a gradient step, whereas the non-smooth part g is handled by the proximal map proxhg . If g = ιC is the convex indicator of a closed convex set C, then the proximal step is the Euclidean projection onto C, and the method becomes projected gradient descent; if g = λ ∥⋅∥1 , then the method becomes the celebrated iterative shrinkage thresholding algorithm (ISTA) [BT09]. In our work, we are interested in the sampling counterpart of composite optimization, called composite log-concave sampling. Here, the goal is to sample from a distribution π on Rd with density π(x) ∝ exp(−F (x)) ∶= exp(−f (x) − g(x)) . Therefore, it is natural to ask: What is the analogue of the proximal gradient algorithm for composite log-concave sampling? Previous studies on log-concave sampling mainly focused on the (non-composite) case g = 0 with methods such as the unadjusted Langevin algorithm (ULA) [Che+25; Dal17; DM17; DM19; VW19] and the Metropolis-adjusted Langevin algorithm (MALA) [Che+21; Dwi+19; WSC22]. In contrast, our goal is to incorporate the non-smooth component g directly. To deal with non-trivial g, a natural first approach is to smooth g by replacing it with its Moreau–Yosida envelope 1 g h (y) ∶= inf {g(x) + ∥y − x∥2 } , 2h x∈Rd where h > 0, so that f + g h is smooth and Langevin algorithms can be applied to the smoothed potential. This idea underlies the Moreau–Yosida regularized unadjusted Langevin algorithm (MYULA) [Ber18; Bro+17; DMP18; Per16]. However, this strategy converges to a biased limit π h (usually not equal to π) and suffers from poor dimensional dependence, especially in the case when g = ιC in which case it reduces to the projected Langevin algorithm [BEL18]. A more direct approach is the proximal sampler introduced in [LST21]. Their seminal work puts forth that the sampling analogue of the proximal map of g is the restricted Gaussian oracle (RGO), which returns an exact sample from a Gaussian-tilted version of g: RGOg,h,y (x) ∝ exp(−g(x) −
1 ∥y − x∥2 ) . 2h
A similar idea appeared in the earlier work of [Mou+22], which used the RGO as a proposal distribution inside a Metropolis–Hastings algorithm; however, their algorithm requires computation of the normalizing constant of RGOg,h,y in addition to a sample. In [LST21], the proximal reduction framework considers the lifted distribution π̃(x, y) ∝ exp(−f (x) − g(y) − 2
1 ∥x − y∥2 ) . 2h
Their algorithm applies Gibbs sampling to π̃, which amounts to alternating between the RGOs for f and g, where the former is implemented via rejection sampling. After running Gibbs sampling and keeping the y sample, they still need to correct the resulting y-marginal π̃ Y to the original target distribution π by another rejection sampling step. Both algorithms in [LST21; Mou+22] achieve a ̃ complexity which scales polylogarithmically in 1/ε, where ε is the target accuracy, but incur O(d) dimension dependence. More recently, works in the non-composite setting √ g = 0 have developed faster RGO implementã d) [AC24; FYC23]. In particular, the rejection tions which only incur dimension dependence O( sampling schemes developed in [FYC23] also provide RGO implementations in the composite setting, but they require semi-smoothness assumptions on g (e.g., g is Lipschitz). This makes them unusable for our goal of making no assumptions on g, aside from convexity, in order to accommodate examples such as the convex indicator ιC . Contributions. In this paper, we propose an algorithm which is arguably closer to a true proximal gradient algorithm for composite log-concave sampling. Specifically, we consider applying the proximal sampler algorithm of [LST21], except that we lift to the joint density π(x, y) ∝ exp(−f (x) − g(x) −
1 ∥x − y∥2 ) . 2h
Again, we perform Gibbs sampling, and this time the x-marginal is precisely our target density π ∝ e−f −g , i.e., the algorithm is asymptotically unbiased. The non-trivial Gibbs sampling step requires sampling from the density RGOf +g,h,y (x) ∝ exp(−f (x) − g(x) −
1 ∥y − x∥2 ) . 2h
Here, our idea is to use as a proposal RGOg,h,y−h∇f (y) (x) ∝ exp(−f (y) − ⟨∇f (y), x − y⟩ − g(x) −
1 ∥y − x∥2 ) . 2h
Then, we apply a Metropolis–Hastings correction step. Note that our proposal directly takes inspiration from the proximal gradient algorithm by linearizing the smooth term, and our resulting algorithm uses the RGO for g alone. Our main result√ shows that our method yields an ε-accurate sample in total variation distance to ̃ the target π in O(κ d log4 (1/ε)) gradient evaluations, where κ is the condition number (smoothness of f divided by the strong convexity of f + g). This is the first high-accuracy sampler for composite √ log-concave sampling with d dimension dependence—matching the state-of-the-art for the noncomposite case g = 0—and which does not impose smoothness assumptions on g. Along the way, we provide an analysis of the independent Metropolis–Hastings algorithm (i.e., Metropolis–Hastings in which the proposal does not depend on the current state) based on easily checkable Rényi divergence conditions, which could be of independent interest. Finally, we extend our results to other standard settings, allowing for non-convex f + g or nonsmooth f , and we provide rudimentary experimental validation to demonstrate that the proposed algorithm is indeed implementable and samples reasonably well.
Other related works Composite sampling. Beyond the approaches discussed in the introduction, other works on composite log-concave sampling with a non-smooth component g include [DMM19; HHP24; SKR19; 3
SR20]. These works differ from ours in terms of the assumptions on the oracle. Specifically, [SKR19; SR20] assume access to proxhg instead of the RGO for g, while [DMM19; HHP24] study subgradientbased algorithms. Moreover, these methods use ULA-type discretizations or splittings with a constant ̃ step size, which incur dimensional dependence O(d) and polynomial dependence on 1/ε, as compared to high-accuracy samplers such as proximal sampler with polylog(1/ε) complexity dependence on the target accuracy ε. Proximal sampler. As discussed above, [Mou+22] proposed an algorithm based on the Metropolis– Hastings correction of the RGO, and the RGO was later incorporated into the proximal sampler framework in [LST21]. The work [Che+22] established convergence results for the ideal proximal sampler under weak log-concavity or functional inequalities, and we leverage these results in our analysis. Both papers [Che+22; LST21] considered simple rejection sampling implementations of the ̃ RGO, which incur dimension dependence O(d), and which were extended to composite/semi-smooth settings [LC23; Yua+23].√Subsequently, [AC24; FYC23] developed RGO implementations with ̃ d), with the former paper working in the composite setting but in which dimension dependence O( g is semi-smooth. More recently, [Che+26b] showed that the RGO implementation does not require function evaluations (only gradient evaluations), and [Che+26a] considered the role of stochastic gradient noise. Sampling from convex bodies. Sampling from a convex body C ⊆ Rd is a special case of composite log-concave sampling where we set f ≡ 0 and g = ιC . The RGO then reduces to sampling from a Gaussian restricted to C. Using this RGO, [KVZ24] pioneered the use of the proximal sampler for sampling convex bodies. A line of works [Koo25; KVZ24; KV25a; KV25b; KZ25] then studied the complexity of implementing the RGO using membership queries to C. Also, [DL25] studied RGO implementation given access to either a projection oracle or a separation oracle for C. The focus of our work is different, as we assume an RGO for C and show how this leads to an improved algorithm for composite log-concave sampling.
2
Preliminaries
2.1
Divergences between probability measures
We recall the divergences used throughout the paper to quantify distances between probability measures. Definition 2.1 (Rényi divergence). For two probability measures µ ≪ ν on Rd , the p-th order Rényi divergence (p > 1) is 1 dµ p Rp (µ∥ν) = log ∫ ( ) dν . p−1 dν Note that the Rényi divergence is also closely related to the KL and χ2 divergences: as p → 1 we have R1 (µ∥ν) ∶= KL(µ∥ν); for p = 2 we have R2 (µ∥ν) = log(1 + χ2 (µ∥ν)), where KL(µ∥ν) = ∫
2.2
dµ dµ log dν , dν dν
χ2 (µ∥ν) = ∫ (
dµ dµ 2 2 − 1) dν = ∫ ( ) dν − 1 . dν dν
Functional inequalities
We recall the definitions of the functional inequalities that we use. 4
Definition 2.2 (Poincaré inequality). We say that π satisfies a Poincaré inequality (PI) with constant 1/α if for all compactly supported and smooth test functions ϕ ∶ Rd → R, varπ (ϕ) ≤
1 Eπ [∥∇ϕ∥2 ] . α
Definition 2.3 (Log-Sobolev inequality). We say that π satisfies a log-Sobolev inequality (LSI) with constant 1/α if for all compactly supported and smooth test functions ϕ ∶ Rd → R, entπ (ϕ2 ) ∶= Eπ [ϕ2 log ϕ2 ] − Eπ [ϕ2 ] log Eπ [ϕ2 ] ≤
2 Eπ [∥∇ϕ∥2 ] . α
Under our notation, an LSI implies a PI with the same constant. In particular, by the Bakry– Émery criterion [BÉ85], any α-strongly log-concave distribution satisfies an LSI, and hence also a PI, with the constant 1/α.
2.3
Proximal sampler and restricted Gaussian oracle
Suppose our target distribution on Rd has density π ∝ exp(−V ), where V ∶ Rd → R ∪ {∞}. Fix a step size h > 0 and extend the distribution to Rd × Rd as: π(x, y) ∝ exp(−V (x) −
1 ∥x − y∥2 ) . 2h
The proximal sampler runs Gibbs sampling on this augmented distribution. Definition 2.4 (Proximal sampler). Initialize X0 ∼ µ0 . For k = 0, 1, 2, . . .: 1 • Sample Yk ∼ π Y ∣X=Xk (y) = N (Xk , hId ) ∝ exp (− 2h ∥y − Xk ∥2 ). 1 • Sample Xk+1 ∼ π X∣Y =Yk (x) ∝ exp (−V (x) − 2h ∥x − Yk ∥2 ).
As we can see from the definition, in order to implement the proximal sampler, we need to be able to sample from both conditional distributions, specifically from π X∣Y =Yk , whose density is a Gaussian tilt centered at Yk of the target distribution π. This idea is formulated in the following definition of the restricted Gaussian oracle [LST21]. Definition 2.5 (Restricted Gaussian oracle). A restricted Gaussian oracle for a convex function φ ∶ Rd → R ∪ {∞} with parameter h > 0 and center v ∈ Rd (RGOφ,h,v ) returns a sample from the density proportional to 1 exp(−φ(x) − ∥x − v∥2 ) . 2h The RGO is exactly the oracle needed to implement the second update of the proximal sampler with φ = V and v = Yk . In the composite setting V = f + g of this paper, we do not assume access to the RGO for V . Instead, we only assume access to RGO for the non-smooth term g, together with gradient access to the smooth term f .
3
Main results
Theorem 3.1. Let π be a distribution on Rd with density π ∝ e−f −g , where f is αf -convex and β-smooth and g is αg -strongly convex, with β ≥ αg ≥ 0. Let α ∶= αf + αg > 0, κ ∶= 1 ∨ β/α, and ε ∈ (0, 1). Assume access to x∗ = arg min{f + g}, proxg , and the RGO for g. Then, with a suitable √ ̃ choice of parameters, Algorithm 1 uses O(κ d log4 (1/ε)) evaluations of f , ∇f , and the RGO for g, and achieves ε total variation distance to π. 5
Remark 3.2. We make several clarifications on the assumptions of our main theorem. • The parameter αf is allowed to be negative, as long as αf ≥ −β and α > 0. In addition, we do not assume smoothness for g, so g can be non-smooth. • We can shift f and g so that x∗ becomes their shared minimizer, as established in [LST21, Proposition 23]. Throughout the paper, we implicitly assume that this preprocessing step has been done. • When applying Algorithm 2, we only need to compute the densities of the proposal and of the target (µk and νk respectively, in Algorithm 1) up to proportionality. In particular, we do not require the normalizing constant for the RGO for g, unlike [Mou+22]. The proof of this theorem is in Appendix A.7. The key ingredients of our theorem are the proximal sampler (Definition 2.4), stated in Algorithm 1, and the implementation of the conditional distribution π X∣Y using the independent Metropolis–Hastings algorithm (Algorithm 2) elaborated √ ̃ in §4. We show that a larger choice of h = Θ(1/(β d)) controls the Rényi divergence between the proposal and the target, which suffices for the analysis by Theorem 4.1. Since the iteration ̃ complexity of the proximal sampler typically scales as O(1/(αh)) (neglecting dependence on other parameters), our larger choice of h improves the overall dependence on the dimension. Algorithm 1 CompositeSampler(f, g, h, K, ε) Require: f αf -convex and β-smooth; g αg -convex; step size h > 0; K ∈ N∗ ; ε > 0. 2β−α Initialize x0 ∼ ρ0 ∝ exp(−g(⋅) − 2 g ∥ ⋅ −x∗ ∥2 ) = RGOg,(2β−αg )−1 ,x∗ . for k = 0 to K − 1 do Sample yk ∼ π Y ∣X=xk = N (xk , hId ). 1 Define νk (x) ∶= π X∣Y =yk (x) ∝ exp(−f (x) − g(x) − 2h ∥x − yk ∥2 ). 1 Define µk (x) ∶= RGOg,h,yk −h∇f (yk ) (x) ∝ exp(−⟨∇f (yk ), x − yk ⟩ − g(x) − 2h ∥x − yk ∥2 ). Sample xk+1 by IndependentMetropolis(νk , µk , ε). Return xK .
Extensions to other settings. We also extend our results to other standard settings. Here LSI(1/α) and PI(1/α) denote the log-Sobolev and Poincaré inequalities with constant 1/α, respectively. Theorem 3.3. Let π be a distribution on Rd with density π ∝ e−f −g and ε ∈ (0, 1). Assume access to x∗ = arg min{f + g}, proxg , prox(f +g) , and the RGO of g, where g is convex. Then, there are algorithms that sample from π up to ε error in total variation distance, with complexities given in the following table (up to hidden logarithmic factors).
f is β-smooth f is L-Lipschitz
π satisfies LSI(1/α) π satisfies PI(1/α) π log-concave √ √ √ β d W22 5 1 2 ) (log κ κ d log5 KL d χ + log ε ε ε2 4 KL L2 α log ε
L2 (log χ2 + log4 1ε ) α
L2 W22 ε2
Here, we use the shorthand KL ∶= KL(ρ0 ∥π), χ2 ∶= χ2 (ρ0 ∥π), W22 ∶= W22 (ρ0 , π), κ ∶= 1 ∨ β/α. Also, in the Lipschitz case, we assume that f is continuously differentiable. The proof of this theorem is in Appendix A.7. 6
4
Analysis of the independent Metropolis–Hastings algorithm
The analysis of our RGO implementation is based on an analysis of the independent Metropolis– Hastings (MH) chain, that is, the special case of the MH chain in which the proposal distribution does not depend on the current state. Since this could be of independent interest, we isolate this part of our proof into the following general result. Theorem 4.1. There is an absolute constant C > 0 such that the following holds. Let ε > 0 and assume that our target distribution π satisfies the Cheeger isoperimetric inequality with some constant Ch > 0 and that the independent proposal kernel µ satisfies Rq (µ∥π) ≤ CR q γ
and
Rq (π∥µ) ≤ CR q γ
for all 1 ≤ q ≤ C log(1/ε), an absolute constant γ > 0, and some constant CR > 0 such that CR ≪ (log(1/ε))−(γ+1) . Assume further that both µ and π are atomless. Then, the independent Metropolis–Hastings algorithm returns a sample whose marginal distribution is ε-close to π in total variation distance in N = O(log(1/ε)) iterations. Algorithm 2 IndependentMetropolis(π, µ, ε) Require: Target density π; proposal µ; ε > 0. Define N = O(log(1/ε)). Initialize x0 ∼ µ. for n = 0 to N − 1 do Sample U ∼ Unif(0, 1). Sample z ∼ µ. π(z) µ(xn ) ). Compute α(xn , z) = 12 min(1, π(xn ) µ(z) if U ≤ α(xn , z) then xn+1 ← z. else xn+1 ← xn . Return xN . The proof is in Appendix B.2. The crucial part of the analysis is the s-conductance cs of the MH chain with the proposal µ independent of the current state. By using the 12 -lazy version of the MH chain, as reflected in the algorithm, we show that cs ≳ 1. As the target distribution satisfies the Cheeger isoperimetric inequality, together with control of the Rényi divergence, we obtain the dimension-free complexity.
5
Examples
We give examples of functions g for which the RGO for g is easily implementable; see Appendix C for further details for some of the examples. • If g1 , . . . , gk admit easily implementable RGOs and g(x1 , . . . , xk ) ∶= ∑ki=1 gi (xi ), then RGOg,h,y = ⊗ki=1 RGOgi ,h,yi . In particular, since RGOs for one-dimensional functions are generally easy to implement (e.g., via rejection sampling), separable functions admit implementable RGOs. This covers numerous examples such as indicators of boxes and the ℓ1 penalty g(x) ∶= λ ∥x∥1 . 7
• If g(x) ∶= g̃(x − c) for some c ∈ Rd , then RGOg,h,y = (x ↦ x + c)# RGOg̃,h,y−c . • If g(x) ∶= g̃(x) + a2 ∥x∥2 + ⟨b, x⟩ + c, then RGOg,h,y = RGOg̃, h , y−hb . 1+ah 1+ah
• If g(x) ∶= ḡ(Bx) for a matrix B ∈ Rk×d , then a sample from RGOg,h,y is produced as follows. First, sample x̄ ∈ Rk from the density µ(x̄) ∝ exp(−ḡ(x̄) −
1 ⟨x̄ − By, (BB T )−1 (x̄ − By)⟩) . 2h
Then, sample x ∼ N (y + B T (BB T )−1 (x̄ − By), h(Id − B T (BB T )−1 B)). In particular, when k = 1 and B = bT , µ̄ = RGOḡ,h∥b∥2 ,⟨b,y⟩ is a one-dimensional RGO. This covers examples such as convex indicators of half-spaces {⟨b, ⋅⟩ ≤ c} and slabs {c ≤ ⟨b, ⋅⟩ ≤ c}. • For a quadratic g(x) ∶= 12 ⟨x, A x⟩ + ⟨b, x⟩ + c, RGOg,h,y = N (Ah (y/h − b), Ah ), where A−1 h = A + h−1 Id , provided that A + h−1 Id ≻ 0. • The ℓ∞ penalty g(x) ∶= λ ∥x∥∞ also admits a closed-form sampler (see Appendix C).
6
Numerical experiments
We test the effectiveness of Algorithm 1 (Composite Sampler) with two baseline algorithms [Per16], namely, the proximal Metropolis-adjusted Langevin algorithm (Prox–MALA), which adds a Metropolis– Hastings correction to the proposal Q(x, ⋅) = N (proxhg (x − h ∇f (x)), 2hId ) , and the proximal gradient Langevin algorithm (PGLA), whose iterates are xk+1 = proxhg (xk − h ∇f (xk ) + N (0, 2hId )) . We tune the step sizes for each method to optimize performance. For both experiments, we take f to correspond to the negative log-likelihood function for logistic regression: n
f (x) = ∑{log(1 + exp(⟨ai , x⟩)) − yi ⟨ai , x⟩} + i=1
τ ∥x∥22 , 2
where τ > 0, and ai ∈ Rd are the rows of the design matrix A ∈ Rn×d . As this is primarily a theoretical paper, the following three experiments are not meant to be comprehensive and only serve to illustrate a proof of concept. The Jupyter notebook containing the code for the numerical experiments is available in the supplementary material. The code chunks run in less than a minute or two in total by the default runtime type on Google Colab for the first two, and it takes around 10 minutes to run the third.
6.1
Sparse Bayesian logistic regression with ℓ1 prior
Here, we take the ℓ1 prior g = λ ∥⋅∥1 . We use the following settings: d = 36, n = 360, τ = 0.2, λ = 7; √ the design matrix A has rows drawn from N (0, Σ) rescaled to have norm n, where Σi,j = ρ∣i−j∣ and ρ = 0.65; the ground truth is chosen to be sparse: xtrue ∶= (1, −1, 0.8, 1.2, −0.9, 0, . . . , 0); and the labels {yi }i∈[n] are drawn from the logistic regression model with ground truth xtrue . 8
Sparse logistic 1: posterior mean error
10 1
Sparse logistic 1: marginal posterior intervals
method
posterior coefficient value
running mean RMSE
Composite sampler Prox-MALA PGLA
10 2 1
10
100 1000 gradient evaluations
0.5 0.0 0.5 1.0
10000
reference 90% interval composite 90% interval
1.0
0
5
10
15 20 coefficient index
25
30
35
Figure 1: (Left) RMSE comparison between the Composite Sampler, Prox-MALA, and PGLA for
sparse Bayesian logistic regression with an ℓ1 prior. (Right) Marginal posterior coverage of the Composite Sampler.
̂∥2 , where x ̂ In Figure 1, we plot the running root mean-squared error (RMSE) √1 ∥ k1 ∑kj=1 xj − x d is the posterior mean, as well as 90% coverage intervals for the coordinates. The ground truths are computed from a long, conservative run of Prox-MALA. For the sake of fair comparison, we plot the performance as a function of gradient evaluations; see Appendix D for more details. The experiment reveals that our composite sampler is competitive with Prox-MALA, and achieves good marginal coverage. In this experiment, PGLA noticeably converges to a smaller RMSE with fewer iterations; however, this is not the end of the story. Due to the application of proxhg , 11.3% of the coordinates of PGLA were exactly 0 at the end of the run, which does not accurately reflect the posterior distribution. Indeed, PGLA suffers from asymptotic bias for h > 0.
6.2
Bayesian logistic regression with box constraints
Here we take the convex indicator of a box, g = ι[−R,R]d and use the following settings: d = 24, n = 360, τ = 0.2, and R = 0.35. The data is generated similarly as before, except that ρ = 0.55 and the ground truth vector is xtrue = (3, −3, 3, −3, 0, . . . , 0). In this setting, the Composite Sampler clearly outperforms Prox-MALA, as shown in Figure 2; this appears to be because many of the proposals for Prox-MALA are infeasible, which is not a problem for the Composite Sampler since it samples each coordinate from the truncated Gaussian directly. PGLA is competitive in terms of RMSE but again suffers from bias: 10.9% of the final coordinates lie on the boundary.
6.3
Gaussian distribution with box constraints
In this experiment, we take the convex indicator of a box, g = ι[−R,R]d , and let κ = 1. Specifically, our target is π(x) ∝ exp(− 12 ∥x∥2 )1[−R,R]d with R = 1. We then vary d ∈ {4, 8, 16, . . . , 512}, and run each method with 3 random seeds. Convergence is measured by sliced W2 distance against the target π. For each method, we set the maximum number of iterations to be 106 . Figure 3 reports the results. It is worth noting√that the Composite Sampler reaches the threshold at every d tested, and the cost grows roughly as d — fitting a line through the eight points gives a log-log slope of approximately 0.5, which agrees with Theorem 3.1. Prox-MALA reaches the threshold within budget for d ≤ 32; PGLA also only reaches the threshold for some d’s, but the dependence of complexity on d is not clear, which may be the variability brought by random seeds. 9
Box-constrained logistic regression: posterior mean error
running mean RMSE
10 1
method
Composite sampler Prox-MALA PGLA
1
10
100 1000 gradient evaluations
10000
Figure 2: RMSE comparison between the Composite Sampler, Prox-MALA, and PGLA for sparse
Bayesian logistic regression with a box constraint.
10 1
1000
10000 gradient evaluations
0.05)
Composite sampler Prox-MALA PGLA threshold = 0.05
Cost-to-threshold vs dimension (Gaussian on box, = 1) log(gradient evaluations to reach W2
sliced Wasserstein-2 distance
Sliced W2 convergence on Gaussian-on-box, d = 32
100000
13 12 11 10 9
Composite sampler Prox-MALA PGLA
8 4 (1.39)
8 (2.08)
16 (2.77)
32 (3.47)
64 (4.16)
d (with logd below)
128 (4.85)
256 (5.55)
512 (6.24)
Figure 3: (Left) Sliced W2 distance to the exact target versus gradient evaluations, d = 64. (Right)
Gradient evaluations needed to first reach the threshold for each method, plotted against d on log-log axes. Hollow markers near the top of the right panel denote runs that did not reach the threshold within the per-method gradient budget.
7
Conclusion
In this paper, we have introduced an algorithm for sampling from composite log-concave target distributions of the form π ∝ e−f −g , inspired by the proximal sampler [Che+22; LST21]. With access to the restricted Gaussian oracle (RGO) of the potentially non-smooth component g, we have shown that √our method achieves high-accuracy sampling guarantees with improved dimensional dependence ̃ O( d), matching the current state-of-the-art results for g ≡ 0. The key insight is to linearize f in the RGO to control the Rényi divergence between the true √ RGO and our implementation via the ̃ tractable proposal (RGOg,h,y−h∇f (y) ) at step size h = Θ(1/(β d)). This larger step size is precisely what drives the improved dimensional dependence. Moreover, we provide verifiable Rényi divergence conditions for independent Metropolis–Hastings in our analysis of the implemented RGO. We have also generalized our results to settings in which π satisfies weaker conditions, such as isoperimetric inequalities, and to cases in which f is Lipschitz.
10
Limitations. Several questions still remain open. Our algorithms assume that the non-smooth component g admits an efficient RGO. While we have given many common examples, including ℓ1 and ℓ∞ penalties, box and half-space indicators, separable functions, and quadratics, it may still be difficult for more complicated composite structures. Our experimental results, while encouraging, are preliminary, and a more extensive study such as in high-dimensional Bayesian inverse problems and constrained sampling settings would help to fully understand the performance of our algorithms in practice at scale. A natural future direction could be an adaptive variant √ of the step size h, which would make the implementation easier, since our current choice h ≍ 1/(β d log(2κ) log2 (1/ζ)) depends on the prescribed choice of the accuracy ζ in the subroutine.
Acknowledgments We would like to acknowledge Andre Wibisono for his helpful suggestions and for pointing us to additional references.
11
A
Proofs of the main results
A.1
Preliminaries
We record a few standard lemmas which we use in our arguments. Lemma A.1 (Herbst’s argument, Ledoux [Led06, §2.3]). Suppose π is α-strongly log-concave and X ∼ π. Then for F that is L-Lipschitz, F (X) − Eπ [F (X)] is L2 /α-sub-Gaussian. Lemma A.2 (Basic lemma, Chewi [Che26, Ch. 4]). Let π ∝ exp(−V ), with V ⪰ αId ≻ 0. If V is minimized at x∗ , then Eπ [∥ ⋅ −x∗ ∥2 ] ≤ d/α. Lemma A.3 (Change of measure, Chewi [Che26, Ch. 6]). Suppose that a test function ϕ satisfies π(ϕ ≥ η) ≤ ψ(η) ,
∀η > 0 .
Then, for any measure µ and q > 1, q−1
µ(ϕ ≥ η) ≤ ψ(η) q exp(
A.2
q−1 Rq (µ∥π)) . q
Rényi divergence lemma
The log-Sobolev inequality bounds the KL divergence by the relative Fisher information, which is useful because the latter can be easier to bound. For instance, the relative Fisher information does not require computing normalizing constants. The following lemma provides a Rényi divergence analogue of this bound, which is helpful for our subsequent arguments. Definition A.4 (Relative Fisher information). For two probability measures µ ≪ ν on Rd , their relative Fisher information is defined as µ 2 FI(µ∥ν) = Eµ [∥∇ log ∥ ] . ν If ν satisfies LSI(1/α), by definition of log-Sobolev inequality 2.3, let ϕ = KL(µ∥ν) ≤
√
µ/ν, we have
1 FI(µ∥ν) . 2α
Lemma A.5 (Upper bound for Rényi divergence). Let π satisfy an LSI with constant 1/α. Then for any q > 1 and any λ > q 2 /(2α), Rq (µ∥π) ≤
q 2 /(2α) µ 2 log E exp(λ ∥∇ log ∥ ). π λ − q 2 /(2α) π
Also, if λ > q(q − 1)/(2α), Rq (µ∥π) ≤
q 2 /(2α) µ 2 log Eµ exp(λ ∥∇ log ∥ ) . λ − q(q − 1)/(2α) π
12
Proof. As in Lemma 5 of [VW19], by the log-Sobolev inequality, we have µ
Rq (µ∥π) ≤
µ
q 2 Eπ [( π )q ∥∇ log π ∥2 ] ⋅ . 2α Eπ [( πµ )q ]
Let γ = ( πµ )q π/Eπ [( πµ )q ]. The above reads Rq (µ∥π) ≤
µ 2 q2 Eγ [∥∇ log ∥ ] . 2α π
(A.1)
By Donsker–Varadhan duality, for any λ > 0, we have µ 2 1 µ 2 Eγ [∥∇ log ∥ ] ≤ {KL(γ∥π) + log Eπ exp(λ ∥∇ log ∥ )} . π λ π
(A.2)
By applying the log-Sobolev inequality again, we have KL(γ∥π) ≤
q2 µ 2 1 FI(γ∥π) = Eγ [∥∇ log ∥ ] . 2α 2α π
Substituting into (A.2) and rearranging terms, we have µ 2 µ 2 Eγ [∥∇ log ∥ ] ≤ (λ − q 2 /(2α))−1 log Eπ exp(λ ∥∇ log ∥ ) . π π Thus, (A.1) becomes Rq (µ∥π) ≤
µ 2 q 2 /(2α) log Eπ exp(λ ∥∇ log ∥ ) . 2 λ − q /(2α) π
The second statement follows similarly, except that we have KL(γ∥µ) ≤
A.3
q−1 KL(γ∥π) . q
RGO implementation
The key step in our analysis of Algorithm 1 is to implement the conditional probability π X∣Y , the RGO of f + g. Suppose we are given Y = y ∈ Rd and define 1 ∥x − y∥2 ) , 2h
(A.3)
1 ∥x − y∥2 ) 2h
(A.4)
ν(x) ∶= π X∣Y =y (x) ∝ exp(−f (x) − g(x) − and our proposal µ(x) ∝ exp(−⟨∇f (y), x − y⟩ − g(x) − with mode
y + = arg min{g(x) + x∈Rd
1 ∥x − (y − h∇f (y))∥2 } = proxhg (y − h∇f (y)) . 2h
The following key lemma controls the Rényi divergence between µ and ν, which will later allow us to invoke Theorem 4.1. Lemma A.6. Let µ and ν be those defined above in Eq. (A.3) and Eq. (A.4). Let g be convex. 13
1. Let p > 1. Assume that f is αf -convex and β-smooth, and that h ≲ 1/(βp) for a sufficiently small implied constant. Then, Rp (µ∥ν) ∨ Rp (ν∥µ) ≲ β 2 p2 h ∥y − y + ∥2 + β 2 p2 dh2 .
(A.5)
2. Let p > 1. Assume that f is L-Lipschitz and continuously differentiable, and that h ≲ 1/L2 . Then, (A.6)
Rp (µ∥ν) ∨ Rp (ν∥µ) ≲ p2 hL2 .
Proof. In the first case, note that ν is (1/h − β)-strongly log-concave, which is at least 1/(2h) when h ≤ 1/β. By Lemma A.5 with λ = 2p2 h, p2 h µ 2 log Eµ exp(2p2 h ∥∇ log ∥ ) 2 2p h − p(p − 1)h ν 2 2 ≤ log Ex∼µ exp(2p h ∥∇f (x) − ∇f (y)∥ )
Rp (µ∥ν) ≤
≤ log Ex∼µ exp(2β 2 p2 h ∥x − y∥2 )
(f is β-smooth)
+ 2
≤ log Ex∼µ exp(4β p h (∥x − y ∥ + ∥y − y + ∥2 )) 2 2
= 4β 2 p2 h ∥y − y + ∥2 + log Ex∼µ exp(4β 2 p2 h ∥x − y + ∥2 ) ,
(A.7)
where y + is the mode of µ. Since µ is 1/h-strongly log-concave and ∥ ⋅ −y + ∥ is 1-Lipschitz, by Herbst’s argument (Lemma A.1), Z = ∥x − y + ∥ − Ex∼µ [∥x − y + ∥]
is h-sub-Gaussian .
With this, we can further bound the latter term in Eq. (A.7): ∥x − y + ∥2 = (Z + Ex∼µ [∥x − y + ∥]) ≤ 2Z 2 + 2 Ex∼µ [∥x − y + ∥] , 2
2
which implies log Ex∼µ exp(4β 2 p2 h ∥x − y + ∥2 ) ≤ 8β 2 p2 h Ex∼µ [∥x − y + ∥2 ] + log Eµ exp(8β 2 p2 hZ 2 ) . For the first term, we apply Lemma A.2 and obtain Ex∼µ [∥x − y + ∥2 ] ≤ dh . For the second term, we use the equivalent definition of sub-Gaussian random variables [Ver18, 2 2 2 2 Proposition 2.5.2]. √ Specifically, h-sub-Gaussianity of Z implies E exp(λ̃ Z ) ≤ exp(C λ̃ ) with 1 ∣λ̃∣ ≤ C and C ≍ h. Hence, log Eµ exp(8β 2 p2 hZ 2 ) ≲ β 2 p2 h2 , where we used h ≲ 1/(βp) with a sufficiently small constant. Thus, Eq. (A.7) becomes: Rp (µ∥ν) ≲ β 2 p2 h ∥y − y + ∥2 + β 2 p2 dh2 + β 2 p2 h2 ≲ β 2 p2 h ∥y − y + ∥2 + β 2 p2 dh2 . The proof for the bound on Rp (ν∥µ) is the same as above: interchange the roles of µ and ν, and use the first inequality of Lemma A.5 instead. For the second statement, since f is differentiable and L-Lipschitz, ∥∇f ∥ ≤ L. Also, the negative log-density of ν is an L-Lipschitz perturbation of a 1/h-strongly convex function, so ν satisfies a log-Sobolev inequality with constant O(1/h) by Brigati and Pedrotti [BP25], provided that h ≲ 1/L2 . We apply Lemma A.5 with λ = Cp2 h for a universal constant C > 0: Rp (µ∥ν) ≤ log Ex∼µ exp(Cp2 h ∥∇f (x) − ∇f (y)∥2 ) ≲ p2 L2 h . The reverse bound also follows easily. 14
The first case of the lemma above suggests that we need a bound for ∥y − y + ∥2 along the iterates of Algorithm 1, which is the main focus of §A.4.
A.4
Concentration along the proximal sampler
In this subsection, we bound ∥y − y + ∥2 along Algorithm 1, where at each step k, the conditional distribution π Xk+1 ∣Yk is exact. In this subsection, we assume that f is αf -convex and β-smooth, and that g is αg -convex, where: α ∶= αf + αg > 0, αf , αg ≥ −β. Suppose Yk = y and recall that y + = proxhg (y − h∇f (y)) is the mode of our proposal µ. By first-step optimality, 1 −∇f (y) − (y + − y) ∈ ∂g(y + ) . h + By the definition of subgradient, ∀p ∈ ∂g(y ), we have g(z) ≥ g(y + ) + ⟨p, z − y + ⟩ ,
∀z ∈ Rd .
We evaluate this expression at z = y and p = −∇f (y) − h1 (y + − y): g(y) ≥ g(y + ) + ⟨−∇f (y) −
1 + (y − y), y − y + ⟩ h
1 ∥y − y + ∥2 − ⟨∇f (y), y − y + ⟩ h 1 Ô⇒ g(y + ) ≤ g(y) + ⟨∇f (y), y − y + ⟩ − ∥y − y + ∥2 . h = g(y + ) +
On the other hand, by β-smoothness of f , f (y + ) ≤ f (y) + ⟨∇f (y), y + − y⟩ +
β ∥y − y + ∥2 . 2
Add up the two inequalities and define F ∶= f + g. We obtain 1 β F (y + ) ≤ F (y) − ( − ) ∥y − y + ∥2 . h 2 Since h ≤ 1/β, we have ∥y − y + ∥2 ≤ 2h (F (y) − F (y + )) ≤ 2h (F (y) − F∗ ) ,
(A.8)
where F∗ ∶= inf F . Note that the argument above holds for any y. We first prove the concentration of F − F∗ under the stationary distribution π ∝ e−f −g . Then, we use Lemma A.3 to establish concentration at the initial distribution for the y variable. Finally, we apply the data-processing inequality to show that at each iteration k in Algorithm 1, we still have a suitable concentration bound. Lemma A.7. Let F be convex, F∗ ∶= inf F , and π ∝ exp(−F ). 1. For all λ ∈ (0, 1), Eπ exp(λ (F − F∗ )) ≤ (1 − λ)−d . 2. For all δ ∈ (0, 1), √ π(F − F∗ ≥ d +
2d log(1/δ) + log(1/δ)) ≤ δ . 15
Proof. For the first statement, exp (λ (F (x) − F∗ ) − F (x)) dx ∫ exp(−F (x)) dx exp(−(1 − λ) (F (x) − F∗ )) dx =∫ . ∫ exp(−(F (x) − F∗ )) dx
Eπ [exp(λ (F − F∗ ))] = ∫
By convexity of F , (1 − λ)F (x) + λF∗ ≥ F ((1 − λ)x + λx∗ ) , which implies (1 − λ) (F (x) − F∗ ) ≥ F ((1 − λ)x + λx∗ ) − F∗ . Thus, the ratio becomes ∫ exp(−(1 − λ) (F (x) − F∗ )) dx ≤ ∫ exp(−F ((1 − λ)x + λx∗ ) + F∗ ) dx ∫ exp(−(F (x) − F∗ )) dx ∫ exp(−(F (x) − F∗ )) dx exp(−F ((1 − λ)x + λx∗ )) dx =∫ ∫ exp(−F (x)) dx exp(−F (y)) ∣det[(1 − λ)Id ]∣−1 dy =∫ ∫ exp(−F (x)) dx = (1 − λ)−d . Then, log Eπ exp(λ (F − F∗ − d)) ≤ d log
1 dλ2 − dλ ≤ . 1−λ 2 (1 − λ)
According to the discussion in [BLM13, §2.4], F − F∗ − d has a sub-gamma right tail with variance factor d and scale parameter 1, and therefore satisfies the tail bound √ π(F − F∗ − d ≥ 2d log(1/δ) + log(1/δ)) ≤ δ . This proves the second statement. √ By the AM–GM inequality, 2d log(1/δ) ≤ 12 (2d + log(1/δ)). Based on Lemma A.7, for t > 0 and ζ ∈ (0, 1), letting log(1/δ) = log(1/ζ) + t, we obtain π(F − F∗ ≥ 2d +
3 3 3 log(1/ζ) + t) = π(F − F∗ ≥ 2d + log(1/δ)) 2 2 √2 ≤ π(F − F∗ ≥ d + 2d log(1/δ) + log(1/δ)) ≤ δ = ζe−t .
(A.9)
We next show how to transfer this concentration bound to the iterates by change of measure, beginning with a bound on the Rényi divergence at initialization. Recall that in Algorithm 1, the initial distribution ρ0 is defined as ρ0 (x) ∝ exp(−g(x) − where x∗ is the shared minimizer of f and g. 16
2β − αg ∥x − x∗ ∥2 ) , 2
Lemma A.8. Assume that f is αf -convex and β-smooth, and g is αg -convex, where αg ≥ −β and α ∶= αf + αg > 0. Then, for κ ∶= β/α, R∞ (ρ0 ∥π) ≤
d log(2κ) . 2
Proof. For any x ∈ Rd , 2β−α
g ρ0 (x) exp(−g(x) − 2 ∥x − x∗ ∥2 ) ∫ exp(−f (η) − g(η)) dη = ⋅ 2β−αg π(x) exp(−f (x) − g(x)) ∥η − x∗ ∥2 ) dη ∫ exp(−g(η) −
2 αf 2β−αg 2 ∥x − x ∥ ) exp(−g(η) − f − ∥η − x∗ ∥2 ) dη ∗ ∗ ∫ 2 2 ≤ ⋅ g exp(−f∗ − β2 ∥x − x∗ ∥2 ) ∫ exp(−g(η) − 2β−α ∥η − x∗ ∥2 ) dη 2 αf 2 2β − αg − αf d/2 ∫ exp(−g(η) − 2 ∥η − x∗ ∥ ) dη = ≤ (1 + ) = (2κ)d/2 , 2β−αg 2 α ∫ exp(−g(η) − 2 ∥η − x∗ ∥ ) dη
exp(−
where the last inequality follows from Lee, Shen, and Tian [LST21, Proposition 61] with f = α g + 2f ∥ ⋅ −x∗ ∥2 and λ1 = 2β − αg . Now we apply Lemma A.3 (change of measure) with q = ∞ and ϕ = F − F∗ − 2d − 3 log(1/ζ). By Eq. (A.9), we have: ρ0 (F − F∗ − 2d −
3 3 d log(1/ζ) ≥ t) ≤ ζ exp(−t + log(2κ)) ≤ ζ , 2 2 2
where we choose t = d2 log(2κ). By the data-processing inequality, for each k ∈ N∗ , we have R∞ (ρk ∥π) ≤ R∞ (ρ0 ∥π) . Therefore, by Lemma A.3, Eq. (A.9), Lemma A.8, and our choice of t, ρk (F − F∗ − 2d −
3 3 log(1/ζ) ≥ t) ≤ ζ exp(−t + R∞ (ρk ∥π)) 2 2 ≤ ζ exp(−t + R∞ (ρ0 ∥π)) ≤ ζ .
√ D Denote yk = xk + hξ, where ξ ∼ N (0, Id ) is independent of xk . Let x+k ∶= proxhg (xk − h∇f (xk )) and yk+ ∶= proxhg (yk − h∇f (yk )). Since proxhg is non-expansive, ∥x+k − yk+ ∥ = ∥proxhg (xk − h∇f (xk )) − proxhg (yk − h∇f (yk ))∥ ≤ ∥xk − h∇f (xk ) − (yk − h∇f (yk ))∥ ≤ ∥xk − yk ∥ + h∥∇f (xk ) − ∇f (yk )∥ ≤ (1 + βh)∥xk − yk ∥ D √ ≤ 2∥xk − yk ∥ = 2 h∥ξ∥ , since h ≤ 1/β, which implies ∥yk − yk+ ∥2 ≤ 3 (∥yk − xk ∥2 + ∥xk − x+k ∥2 + ∥x+k − yk+ ∥2 ) ≤ 3 ∥xk − x+k ∥2 + 9h ∥ξ∥2 . 17
(A.10)
√ By Gaussian concentration, we have that for any t1 = 2 log(1/ζ) > 0, √ P(∥ξ∥ ≥ d + t1 ) ≤ exp(−t21 /2) = ζ . Combined with Eq. A.8 and Eq. A.10, with probability ≥ 1 − 2ζ under the y-marginal at the k-th iteration, 3 3 ∥yk − yk+ ∥2 ≤ 3 ⋅ 2h (2d + d log(2κ) + log(1/ζ)) + 9h ⋅ (2d + 4 log(1/ζ)) 4 2 9 = 30dh + dh log(2κ) + 45h log(1/ζ) . 2
A.5
(A.11)
Analysis using the perfect RGO
We recall the following results from Chen, Chewi, Salim, and Wibisono [Che+22]. Lemma A.9 (Convergence of the ideal proximal sampler). Let (ρk )k∈N be the iterates of the proximal sampler with step size h and stationary distribution π. 1. Let π satisfy an LSI with constant 1/α. Then, KL(ρk ∥π) ≤
KL(ρ0 ∥π) . (1 + αh)2k
2. Let π satisfy a PI with constant 1/α. Then, χ2 (ρk ∥π) ≤
χ2 (ρ0 ∥π) . (1 + αh)2k
KL(ρk ∥π) ≤
W22 (ρ0 , π) . 2kh
3. Let π be log-concave. Then,
A.6
Error analysis for the proximal sampler
In this subsection, we denote ρk and ρ̂k as the marginal distributions for x along the proximal sampler using the perfect RGO π X∣Yk and our implemented RGO π̂ X∣Yk at the k-th iteration in Algorithm 1, respectively. Note that ρ0 = ρ̂0 as defined in the algorithm. We also denote the y-marginals along the perfect RGO and the implemented RGO at the k-th iteration as ρYk and ρ̂Yk , respectively. The tools that we use for the analysis are the results in the previous subsections that hold with high probability, as well as the following lemma on coupling: Lemma A.10 ([WLP09, Proposition 4.7]). Let µ and ν be two probability distributions on Ω. Then ∥µ − ν∥TV = inf{P(X ≠ Y ) ∣ (X, Y ) is a coupling of µ and ν} . Let ζ ∈ (0, 1). Suppose at the k-th iteration, we have iterates xk ∼ ρk and x̂k ∼ ρ̂k . By Lemma A.10, there exists a coupling of ρk and ρ̂k for which xk = x̂k with √ probability 1√− ∥ρk − ρ̂k ∥TV . Draw ξ ∼ N (0, Id ) independent of this coupling and set yk ∶= xk + hξ, ŷk ∶= x̂k + hξ. Then, we have yk = ŷk with the same probability under this coupling of ρYk and ρ̂Yk . At the (k + 1)-th step, xk+1 ∼ ρYk π X∣Yk and x̂k+1 ∼ ρ̂Yk π̂ X∣Yk . We are able to show that π X∣Yk =y and our implementation π̂ X∣Yk =y are close with the following corollary. 18
Corollary A.11 (RGO implementation). It holds that ∥π X∣Yk =yk − π̂ X∣Yk =yk ∥TV ≤ ζ
after N = O(log(1/ζ)) iterations
in the following two cases. 1. Let f be αf -convex and β-smooth and g be αg -convex, where αg ≥ 0 and α ∶= αf + αg > 0. Here, we assume that yk ∈ E, where E√is the event on which ∥yk − yk+ ∥2 enjoys the concentration bound (A.11), and that h ≍ 1/(β d log(2κ) log2 (1/ζ)) with a sufficiently small constant. 2. Let f be continuously differentiable and L-Lipschitz, and let g be convex. Here, we take h ≍ 1/(L2 log3 (1/ζ)) for a sufficiently small constant. Proof. This corollary directly follows from Theorem 4.1 on the independent Metropolis–Hastings algorithm. Let ζ ∈ (0, 1). In the first case, since yk ∈ E, the concentration bound in §A.4 reads ∥yk − yk+ ∥2 ≲ dh log(2κ) + h log(1/ζ) . Substituting into Eq. (A.5), we have Rp (µ∥ν) ∨ Rp (ν∥µ) ≲ β 2 p2 dh2 + β 2 p2 h (dh + dh log(2κ) + h log(1/ζ)) ≲ β 2 dh2 log(2κ) log(1/ζ) p2 . ´¹¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¸¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¶ CR
√ By choosing h ≍ 1/(β d log(2κ) log2 (1/ζ)) with a sufficiently small constant, we can apply Theorem 4.1 with γ = 2, thus proving the desired result. In the second case, we instead use Eq. (A.6), which states that Rp (µ∥ν) ∨ Rp (ν∥µ) ≲ hL2 p2 . ´¸¶ CR
Here, we choose h ≍ 1/(L2 log3 (1/ζ)) and apply Theorem 4.1 with γ = 2. By the triangle inequality and the data-processing inequality, ∥ρk+1 − ρ̂k+1 ∥TV = ∥ρYk π X∣Yk − ρ̂Yk π̂ X∣Yk ∥TV ≤ ∥ρYk π X∣Yk − ρYk π̂ X∣Yk ∥TV + ∥ρYk π̂ X∣Yk − ρ̂Yk π̂ X∣Yk ∥TV ≤ ∥ρYk π X∣Yk − ρYk π̂ X∣Yk ∥TV + ∥ρYk − ρ̂Yk ∥TV ≤ ∥ρYk π X∣Yk − ρYk π̂ X∣Yk ∥TV + ∥ρk − ρ̂k ∥TV .
(A.12)
Given Yk = y, if y ∈ E, draw (xk+1 , x̂k+1 ) from the maximal coupling as in Lemma A.10; otherwise, draw (xk+1 , x̂k+1 ) from the independent coupling π X∣Yk ⊗ π̂ X∣Yk . Then, ∥ρYk π X∣Yk − ρYk π̂ X∣Yk ∥TV ≤ ∥Law(y, xk+1 ) − Law(y, x̂k+1 )∥TV
(data-processing inequality)
≤ P(xk+1 ≠ x̂k+1 ) (by Lemma A.10) = P(xk+1 ≠ x̂k+1 , y ∉ E) + P(xk+1 ≠ x̂k+1 , y ∈ E) ≤ ρYk (E c ) + P(xk+1 ≠ x̂k+1 , y ∈ E) ≤ 2ζ + ζ = 3ζ , 19
where for the first term we used ρYk (E) ≥ 1 − 2ζ from §A.4 and for the second term we note that under our constructed coupling for y ∈ E, by Corollary A.11, P(xk+1 ≠ x̂k+1 , y ∈ E) = ∫ ∥π X∣Yk =y − π̂ X∣Yk =y ∥TV ρYk (dy) ≤ ζ . E
Therefore, the induction step in Eq. (A.12) becomes ∥ρk+1 − ρ̂k+1 ∥TV ≤ 3ζ + ∥ρk − ρ̂k ∥TV . Since ρ0 = ρ̂0 , by induction, we have ∥ρK − ρ̂K ∥TV ≤ 3Kζ .
A.7
Proof of Theorems 3.1 and 3.3
√ Proof of Theorem 3.1. Let ε, ζ ∈ (0, 1), and h ≍ 1/(β d log(2κ) log2 (1/ζ)) as in Corollary A.11. Since π ∝ e−f −g is α-strongly log-concave, it satisfies an LSI with constant 1/α. Hence, by Pinsker’s inequality, Lemma A.9, and monotonicity of Rényi divergence, we have ¿ ¿ √ Á KL(ρ0 ∥π) Á R∞ (ρ0 ∥π) 1 À À ∥ρk − π∥TV ≤ KL(ρk ∥π) ≤ Á ≤Á . (A.13) 2k 2 2 (1 + αh) 2 (1 + αh)2k log
d log(2κ)
4ε2 By Lemma A.8, R∞ (ρ0 ∥π) ≤ d2 log(2κ), so ∥ρK −π∥TV ≤ ε if K ≥ 2 log(1+αh) . Since log(1+x) ≥ x−x2 /2 for all x ≥ 0, log(1 + αh) ≥ αh (1 − αh/2) ≥ αh/2 because αh < 1. Then it suffices to choose
K = O(
d log(2κ) 1 log ). αh ε2
Moreover, by the bounds on total variation in §A.6, we have ∥ρK − ρ̂K ∥TV ≤ 3Kζ, where each iteration has complexity N = O(log ζ1 ). If both terms are bounded by ε/2, then the proof is complete. Furthermore, if the latter term is bounded by ε/2, so is the former because K ∈ N∗ . Hence, we can let ζ = ε/(6K). By plugging in h, K depends on ζ and itself as following: √ √ d log(2κ) 1 K = O(κ d log(2κ) log2 log ). ζ 6Kζ √ ̃ Neglecting logarithmic terms on d and κ, we obtain that K = O(κ d log3 ζ1 ). Since N = O(log ζ1 ) and ζ = ε/(6K), the final complexity is √ 1 ̃ KN = O(κ d log4 ) . ε Proof of Theorem 3.3. We use the proximal sampler algorithm. In all cases, the number of outer loop iteration K is given by Lemma A.9 to be K = O(
1 KL(ρ0 ∥π) log ), αh ε2
K = O(
1 χ2 (ρ0 ∥π) log ), αh ε2
K = O(
W22 (ρ0 , π) ) hε2
(A.14)
under LSI(1/α), PI(1/α), and log-concavity respectively. Then, we only need to consider the complexity N of the implementation of the RGO in each case, where we use ζ ∈ (0, 1) to denote the error of the implemented RGO to the true RGO. 20
Smooth case. Here, f is β-smooth. Choose h = 1/(2β) and for y ∈ Rd , define fh (x) ∶= f (x) +
1 ∥x − y∥2 . 2h
Note that ∇2 fh = ∇2 f + h1 Id , so fh is β-strongly convex and 3β-smooth. Thus, the RGO π X∣Y , whose density is proportional to exp(−f (x) − g(x) −
1 ∥x − y∥2 ) = exp(−fh (x) − g(x)) , 2h
is a composite log-concave density with condition number κ ≤ 3. By Theorem 3.1, we can implement √ 4 1 X∣Yk ̃ π in N = O( d log ζ ) function and gradient evaluations at each k given access to the mode of π X∣Yk , i.e., xk,∗ ∶= proxh(f +g) (Yk ). Combining ζ = ε/(6K) with K defined in Eq (A.14), we obtain the desired complexities under the three different conditions on π. Non-smooth case. Here, f is L-Lipschitz. In this case, according to Corollary A.11, we can choose h ≍ 1/(L2 log3 (1/ζ)), so N = O(log(1/ζ)). Combining this with the K defined in Eq (A.14) yields the results.
B
Analysis of independent Metropolis–Hastings
In this section, we study the independent Metropolis–Hastings algorithm. For any x ∈ Rd , let the proposal kernel be Q(x, ⋅) ∶= Qx ∶= µ. The transition kernel of the Metropolis-adjusted chain is defined to be P (x, dy) ∶= Q(x, dy) A(x, y) + (1 − ∫ Q(x, dy ′ ) A(x, y ′ )) δx (dy) , where δx is the Dirac mass at x and A(x, y) = min(1,
π(y) Q(y, x) π(y) µ(x) ) = min(1, ) π(x) Q(x, y) π(x) µ(y)
is the acceptance probability. Here, π is the stationary distribution. We also denote P (x, ⋅) by Px . Remark B.1. We study 12 -lazy chains in this section. A Markov chain is ℓ-lazy if at each step the chain transitions with the transition kernel P with probability 1 − ℓ and stays at its current position with probability ℓ. By studying lazy chains, we not only remove periodicity but also ensure that the Markov transition kernel is positive semidefinite at the cost of only a constant factor [LS93; Mou+22; Vem05].
B.1
Preliminaries on conductance analysis
We begin by recalling the standard toolbox for conductance analysis. Suppose we have a Markov chain on Rd with P as its transition kernel and π be the stationary distribution. Definition B.2 (Cheeger isoperimetry). We say that π satisfies a Cheeger isoperimetric inequality with constant Ch > 0 if lim inf + ε→0
π(Aε ) − π(A) 1 ≥ π(A) (1 − π(A)) ε Ch
for all measurable sets A ⊆ Rd , where Aε ∶= {x ∈ Rd ∶ dist(x, A) < ε} is the ε blow-up of A. 21
√ If π is α-strongly log-concave, then Ch ≲ 1/ α [Mil09]. Definition B.3 (s-conductance). For s ∈ [0, 12 ], the s-conductance of P is the largest cs > 0 such that for all events A ⊆ Rd , c c ∫ P (x, A ) π(dx) ≥ cs (π(A) − s) (π(A ) − s) . A
Note that for π(A) ≤ s, the inequality holds trivially. Lemma B.4 ([LS93], Corollary 1.6). For any s ∈ (0, 12 ], let ∆s ∶= sup{∣µ0 (A) − π(A)∣ ∶ A ⊆ Rd , π(A) ≤ s} . Then, the law of µN of the N -th iterate of a Markov chain with s-conductance cs , stationary distribution π, and initialized at µ0 satisfies ∥µN − π∥TV ≤ ∆s + In particular,
c2 N ∆s exp(− s ) . s 2 √
√ ∥µN − π∥TV ≤
sχ2 (µ0 ∥π) +
χ2 (µ0 ∥π) c2 N exp(− s ) . s 2 2
ε Remark B.5. Lemma B.4 implies that, by choosing s = 4χ2 (µ , we have ∥µN − π∥TV ≤ ε after 0 ∥π)
1 χ2 (µ0 ∥π) N = O( 2 log ) cs ε2
iterations.
To lower bound the s-conductance, we use the following geometric overlap lemma. Lemma B.6 (Overlap lemma for s-conductance, [Che26], Ch. 7). Assume: 1. π satisfies Cheeger isoperimetric inequality with constant Ch > 0. rs 2. There exists r ∈ [0, Ch] and E ⊆ Rd with probability π(E) ≥ 1 − 16Ch such that
∀x, y ∈ E , ∥x − y∥ ≤ r Ô⇒ ∥P (x, ⋅) − P (y, ⋅)∥TV ≤
1 . 2
Then cs ≳ r/Ch. In order to apply Lemma B.6, we need to control ∥Px − Py ∥TV ≤ 12 for x, y ∈ Rd on a set E with high probability under π. By the triangle inequality, ∥Px − Py ∥TV ≤ ∥Px − Qx ∥TV + ∥Qx − Qy ∥TV + ∥Qy − Py ∥TV = ∥Px − Qx ∥TV + ∥Py − Qy ∥TV ,
(B.1)
where the middle term is 0 for the independent MH algorithm. To control the other terms, we invoke the following lemma. Lemma B.7 (Pointwise projection property, [Che+21], Theorem 6). Let Q be an atomless proposal kernel and P be the corresponding Metropolis–Hastings kernel with target distribution π. Then, for any atomless kernel Q̄ which is reversible with respect to π, and for every x ∈ Rd , ∥Px − Qx ∥TV ≤ 2 ∥Qx − Q̄x ∥TV + ∫
π(y) Q̄(y, x) Q(y, x) ∣ − 1∣ dy . π(x) Q̄(y, x)
We next show how to control these terms under a Rényi divergence condition for the proposal. 22
B.2
Proof of Theorem 4.1
Proof of Theorem 4.1. Since both µ and π are atomless, we can apply Lemma B.7 and obtain that: ∀x ∈ Rd , ∥Px − Qx ∥TV ≤ 2 ∥µ − π∥TV + ∫
π(y) π(x) µ(x) µ(x) ∣ − 1∣ dy = 2 ∥µ − π∥TV + ∣ − 1∣ . π(x) π(x) π(x)
For the first term, by Pinsker’s inequality, √ √ √ 1 1 1 ∥µ − π∥TV ≤ KL(µ∥π) = R1 (µ∥π) ≤ CR . 2 2 2
(B.2)
(B.3)
For the second term, by Markov’s inequality: ∀t ∈ (0, 1), π(∣
µ µ π 1 − 1∣ ≥ t) = π( ≥ 1 + t) + π( ≥ ) π π µ 1−t ≤ (1 + t)−q Eπ [(µ/π)q ] + (1 − t)q−1 Eπ [(π/µ)q−1 ] = (1 + t)−q exp((q − 1) Rq (µ∥π)) + (1 − t)q−1 exp((q − 1) Rq (π∥µ)) ≤ ((1 + t)−q + (1 − t)q−1 ) exp(CR q γ+1 ) .
(B.4)
With t = 18 , Eq. (B.4) becomes: π(∣
1 8 q µ 7 q−1 8 q − 1∣ ≥ ) ≤ [( ) + ( ) ] exp(CR q γ+1 ) ≤ 3 ( ) exp(CR q γ+1 ) . π 8 9 8 9
Define the set E ∶= {∣µ/π − 1∣ ≤ 81 }. By the triangle inequality in Eq. (B.1), Eq. (B.2), Eq. (B.3), and Eq. (B.4), we have ∥Px − Py ∥TV ≤ 4 ∥µ − π∥TV + ∣
µ(x) µ(y) 1 − 1∣ + ∣ − 1∣ ≤ π(x) π(y) 2
1 for all x, y ∈ E, provided that CR ≤ 128 . 2 ε As in Remark B.5, let s = 4χ2 (µ∥π) , where we initialize the Markov chain with µ0 = µ. Then,
χ2 (µ∥π) = exp(R2 (µ∥π)) − 1 ≤ exp(CR 2γ ) = O(1) , rs s so s ≍ ε2 . In Lemma B.6, we can choose r = Ch. To bound π(E c ) by 16Ch = 16 , it suffices to take −(γ+1) q ≍ log(1/s) ≍ log(1/ε) and for CR ≪ q . Then, by Lemma B.6, the s-conductance satisfies cs ≳ r/Ch = 1. Finally, according to Remark B.5, the number of iterations we need to obtain ε-TV distance to π is
1 χ2 (µ∥π) N = O( 2 log ) = O(log(1/ε)) . cs ε2
C
RGO implementations
Here, we provide more details on concrete RGO implementations. Given a set C, we let ιC denote the convex indicator of C. • Indicator of a box. For the indicator of a box, g = ι⨉d [ai ,bi ] , sampling from the RGO i=1 reduces to sampling from one-dimensional truncated Gaussians N (y, h)∣[a,b] . 23
• ℓ1 penalty. For g(x) ∶= λ ∥x∥1 , by separability, it suffices to consider the one-dimensional case. One can see that ⎧ 1 2 ⎪ ⎪exp(− 2h (x − (y + λh)) ) , x ≤ 0 , RGOg,h,y (x) ∝ ⎨ 1 2 ⎪ ⎪ ⎩exp(− 2h (x − (y − λh)) ) , x ≥ 0 . Define the weights w− ∶= exp(λy +
1 2 y + λh λ h) Φ(− √ ) , 2 h
w+ ∶= exp(−λy +
1 2 y − λh λ h) Φ( √ ) , 2 h
where Φ(⋅) is the standard Gaussian CDF. Then, an exact sample from RGO(g, h, y) is produced as follows: with probability w− /(w− + w+ ), draw a sample from N (y + λh, h)∣R− ; otherwise, draw a sample from N (y − λh, h)∣R+ . • Other separable functions. Further examples of separable functions include ℓp norms g(x) ∶= λ ∥x∥pp for p ≥ 1; the elastic net g(x) ∶= λ1 ∥x∥1 + λ22 ∥x∥22 ; hinge penalties g(x) ∶= ∑di=1 λi (xi −bi )+ ; and the entropy map g(x) ∶= ∑di=1 xi log xi + ιRd+ (x). In each of these cases, the RGO reduces to sampling from one-dimensional densities, for which closed-form samplers or tailored rejection sampling schemes can be developed. • ℓ∞ penalty. For g(x) ∶= λ ∥x∥∞ , the density of the RGO is 1 d 2 ∑(xi − yi ) ) 2h i=1
RGOg,h,y (x) ∝ exp(−λ ∥x∥∞ − ∞
∝∫
0 ∞
∝∫
exp(−λt −
0
1 d 2 ∑(xi − yi ) ) 1t≥∥x∥∞ dt 2h i=1 d
exp(−λt) ∏[Zi (t)
1 exp(− 2h (xi − yi )2 ) 1t≥∣xi ∣
Zi (t)
i=1
] dt ,
√ √ i ) − Φ( −t−y √ i )]. This admits the following interpretation: first, where Zi (t) = 2πh [Φ( t−y h h sample t from the following density over R+ . d
p(t) ∝ exp(−λt) ∏ Zi (t) . i=1
Then, conditionally on t, draw X1 , . . . , Xd independently with Xi ∼ N (yi , h)∣[−t,t] . The max penalty g(x) = λ maxi∈[d] xi can be handled similarly.
D
Further experimental details
Measuring gradient evaluations. For each method, we roughly counted the number of gradient evaluations ∇f per iteration as follows. (The number of function evaluations of f is roughly the same.) • PGLA always uses one gradient evaluation per iteration. • Prox-MALA uses two gradient evaluations per iteration: one is used to generate the proposal, and another is needed to compute the acceptance ratio. 24
• In each outer iteration of the Composite Sampler, we use one initial computation of ∇f for the first proposal. Each inner iteration for the RGO implementation further requires one new gradient evaluation in order to evaluate the acceptance probability. In total, each outer iteration requires one evaluation of proxhf and N + 1 evaluations of ∇f , where N is the number of inner iterations. In practice, we capped the number of inner iterations to be Nmax and we count the cost of an outer iteration to be Nmax + 1, although this somewhat overcounts since the number of actual gradient evaluations can be smaller if the proposals are accepted in fewer than Nmax steps.
25
References [AC24]
J. M. Altschuler and S. Chewi. “Faster high-accuracy log-concave sampling via algorithmic warm starts”. In: Journal of the ACM 71.3 (2024), pp. 1–55.
[BÉ85]
D. Bakry and M. Émery. “Diffusions hypercontractives”. In: Séminaire de probabilités 19 (1985), pp. 177–206.
[BT09]
A. Beck and M. Teboulle. “A fast iterative shrinkage-thresholding algorithm for linear inverse problems”. In: SIAM Journal on Imaging Sciences 2.1 (2009), pp. 183–202.
[Ber18]
E. Bernton. “Langevin Monte Carlo and JKO splitting”. In: Conference on Learning Theory. PMLR. 2018, pp. 1777–1798.
[BLM13]
S. Boucheron, G. Lugosi, and P. Massart. Concentration inequalities. A nonasymptotic theory of independence, With a foreword by Michel Ledoux. Oxford University Press, Oxford, 2013, pp. x+481.
[BP25]
G. Brigati and F. Pedrotti. “Heat flow, log-concavity, and Lipschitz transport maps”. In: Electron. Commun. Probab. 30 (2025), Paper No. 71, 12.
[Bro+17]
N. Brosse, A. Durmus, É. Moulines, and M. Pereyra. “Sampling from a log-concave distribution with compact support with proximal Langevin Monte Carlo”. In: Conference on Learning Theory. PMLR. 2017, pp. 319–342.
[BEL18]
S. Bubeck, R. Eldan, and J. Lehec. “Sampling from a log-concave distribution with projected Langevin Monte Carlo”. In: Discrete & Computational Geometry 59.4 (2018), pp. 757–783.
[Che+26a]
F. Chen, S. Chewi, C. Daskalakis, and A. Rakhlin. “High-accuracy log-concave sampling with stochastic queries”. In: arXiv preprint 2602.14342 (2026).
[Che+26b]
F. Chen, S. Chewi, C. Daskalakis, and A. Rakhlin. “High-accuracy sampling for diffusion models and log-concave distributions”. In: arXiv preprint arXiv:2602.01338 (2026).
[Che+22]
Y. Chen, S. Chewi, A. Salim, and A. Wibisono. “Improved analysis for a proximal algorithm for sampling”. In: Conference on Learning Theory. PMLR. 2022, pp. 2984–3014.
[Che26]
S. Chewi. Log-concave sampling. Book draft. Available at https://chewisinho.github.io/ main.pdf. 2026+.
[Che+25]
S. Chewi, M. A. Erdogdu, M. Li, R. Shen, and M. S. Zhang. “Analysis of Langevin Monte Carlo from Poincaré to log-Sobolev”. In: Foundations of Computational Mathematics 25.4 (2025), pp. 1345–1395.
[Che+21]
S. Chewi, C. Lu, K. Ahn, X. Cheng, T. Le Gouic, and P. Rigollet. “Optimal dimension dependence of the Metropolis-adjusted Langevin algorithm”. In: Conference on Learning Theory. PMLR. 2021, pp. 1260–1300.
[Dal17]
A. S. Dalalyan. “Theoretical guarantees for approximate sampling from smooth and log-concave densities”. In: Journal of the Royal Statistical Society Series B: Statistical Methodology 79.3 (2017), pp. 651–676.
[DL25]
T. Dang and J. Liang. “Oracle-based Uniform Sampling from Convex Bodies”. In: arXiv preprint arXiv:2510.02983 (2025).
[DMM19]
A. Durmus, S. Majewski, and B. Miasojedow. “Analysis of Langevin Monte Carlo via convex optimization”. In: J. Mach. Learn. Res. 20 (2019), Paper No. 73, 46.
[DM17]
A. Durmus and É. Moulines. “Nonasymptotic convergence analysis for the unadjusted Langevin algorithm”. In: The Annals of Applied Probability 27.3 (2017), pp. 1551–1587.
[DM19]
A. Durmus and É. Moulines. “High-dimensional Bayesian inference via the unadjusted Langevin algorithm”. In: Bernoulli 25.4A (2019), pp. 2854–2882.
[DMP18]
A. Durmus, E. Moulines, and M. Pereyra. “Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau”. In: SIAM Journal on Imaging Sciences 11.1 (2018), pp. 473–506.
26
[Dwi+19]
R. Dwivedi, Y. Chen, M. J. Wainwright, and B. Yu. “Log-concave sampling: Metropolis–Hastings algorithms are fast”. In: Journal of Machine Learning Research 20.183 (2019), pp. 1–42.
[FYC23]
J. Fan, B. Yuan, and Y. Chen. “Improved dimension dependence of a proximal algorithm for sampling”. In: The Thirty Sixth Annual Conference on Learning Theory. PMLR. 2023, pp. 1473– 1521.
[HHP24]
A. Habring, M. Holler, and T. Pock. “Subgradient Langevin methods for sampling from nonsmooth potentials”. In: SIAM J. Math. Data Sci. 6.4 (2024), pp. 897–925.
[Koo25]
Y. Kook. “Zeroth-order log-concave sampling”. In: arXiv preprint 2507.18021 (2025).
[KVZ24]
Y. Kook, S. S. Vempala, and M. S. Zhang. “In-and-out: algorithmic diffusion for sampling convex bodies”. In: Advances in Neural Information Processing Systems 37 (2024), pp. 108354–108388.
[KV25a]
Y. Kook and S. S. Vempala. “Faster logconcave sampling from a cold start in high dimension”. In: arXiv preprint arXiv:2505.01937 (2025).
[KV25b]
Y. Kook and S. S. Vempala. “Sampling and integration of logconcave functions by algorithmic diffusion”. In: Proceedings of the 57th Annual ACM Symposium on Theory of Computing. STOC ’25. Prague, Czechia: Association for Computing Machinery, 2025, pp. 924–932.
[KZ25]
Y. Kook and M. S. Zhang. “Rényi-infinity constrained sampling with d3 membership queries”. In: Proceedings of the 2025 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA). SIAM. 2025, pp. 5278–5306.
[Led06]
M. Ledoux. “Concentration of measure and logarithmic Sobolev inequalities”. In: Seminaire de probabilites XXXIII. Springer, 2006, pp. 120–216.
[LST21]
Y. T. Lee, R. Shen, and K. Tian. “Structured logconcave sampling with a restricted Gaussian oracle”. In: Proceedings of Thirty Fourth Conference on Learning Theory. Ed. by M. Belkin and S. Kpotufe. Vol. 134. Proceedings of Machine Learning Research. PMLR, Aug. 2021, pp. 2993–3050.
[LC23]
J. Liang and Y. Chen. “A proximal algorithm for sampling”. In: Transactions on Machine Learning Research (2023).
[LS93]
L. Lovász and M. Simonovits. “Random walks in a convex body and an improved volume algorithm”. In: Random Structures & Algorithms 4.4 (1993), pp. 359–412.
[Mil09]
E. Milman. “On the role of convexity in isoperimetry, spectral gap and concentration”. In: Inventiones Mathematicae 177.1 (2009), pp. 1–43.
[Mou+22]
W. Mou, N. Flammarion, M. J. Wainwright, and P. L. Bartlett. “An efficient sampling algorithm for non-smooth composite potentials”. In: Journal of Machine Learning Research 23.233 (2022), pp. 1–50.
[PB14]
N. Parikh and S. Boyd. “Proximal algorithms”. In: Foundations and Trends in Optimization 1.3 (2014), pp. 127–239.
[Per16]
M. Pereyra. “Proximal Markov chain Monte Carlo algorithms”. In: Statistics and Computing 26.4 (2016), pp. 745–760.
[SKR19]
A. Salim, D. Kovalev, and P. Richtarik. “Stochastic proximal Langevin algorithm: potential splitting and nonasymptotic rates”. In: Advances in Neural Information Processing Systems. Ed. by H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché-Buc, E. Fox, and R. Garnett. Vol. 32. Curran Associates, Inc., 2019.
[SR20]
A. Salim and P. Richtarik. “Primal dual interpretation of the proximal stochastic gradient Langevin algorithm”. In: Advances in Neural Information Processing Systems. Ed. by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin. Vol. 33. Curran Associates, Inc., 2020, pp. 3786–3796.
[Vem05]
S. Vempala. “Geometric random walks: a survey”. In: Combinatorial and computational geometry 52.573-612 (2005), p. 2.
27
[VW19]
S. Vempala and A. Wibisono. “Rapid convergence of the unadjusted Langevin algorithm: isoperimetry suffices”. In: Advances in Neural Information Processing Systems 32 (2019).
[Ver18]
R. Vershynin. High-dimensional probability: an introduction with applications in data science. Vol. 47. Cambridge University Press, 2018.
[WLP09]
E. L. Wilmer, D. A. Levin, and Y. Peres. “Markov chains and mixing times”. In: American Mathematical Soc., Providence 107 (2009).
[WSC22]
K. Wu, S. Schmidler, and Y. Chen. “Minimax mixing time of the Metropolis-adjusted Langevin algorithm for log-concave sampling”. In: Journal of Machine Learning Research 23.270 (2022), pp. 1–63.
[Yua+23]
B. Yuan, J. Fan, J. Liang, A. Wibisono, and Y. Chen. “On a class of Gibbs sampling over networks”. In: Proceedings of Thirty Sixth Conference on Learning Theory. Ed. by G. Neu and L. Rosasco. Vol. 195. Proceedings of Machine Learning Research. PMLR, July 2023, pp. 5754–5780.
28