Conceptio › Archive › arXiv CS

arXiv CSopen access

Decentralized Proximal Stochastic Gradient Langevin Dynamics

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

Decentralized Proximal Stochastic Gradient Langevin Dynamics Mohammad Rafiqul Islam 1 , Lingjiong Zhu 2 May 4, 2026

arXiv:2605.00723v1 [stat.ML] 1 May 2026

Abstract We propose Decentralized Proximal Stochastic Gradient Langevin Dynamics (DE-PSGLD), a decentralized Markov chain Monte Carlo (MCMC) algorithm for sampling from a log-concave probability distribution constrained to a convex domain. Constraints are enforced through a shared proximal regularization based on the Moreau–Yosida envelope, enabling unconstrained updates while preserving consistency with the target constrained posterior. We establish nonasymptotic convergence guarantees in the 2-Wasserstein distance for both individual agent iterates and their network averages. Our analysis shows that DE-PSGLD converges to a regularized Gibbs distribution and quantifies the bias introduced by the proximal approximation. We evaluate DE-PSGLD for different sampling problems on synthetic and real datasets. As the first decentralized approach for constrained domains, our algorithm exhibits fast posterior concentration and high predictive accuracy. The code for the numerical experiments is publicly available at anonymous.4open.science/r/Decentralized-proximal-SGLD-74E8.

Introduction

Decentralized learning is a learning process in which data is distributed across computational agents or collected by individual agents, and model parameters are computed as the consensus of the agents. It has gained a lot of interest for applications where agents can collaboratively learn a predictive model without sharing their own data, but sharing only their local models with their immediate neighbors to generate a global model [He et al., 2018, Hendrikx et al., 2019, Arjevani et al., 2020]. We assume there are N agents who are connected over an undirected communication network G = (V, E) where V = {1, . . . , N } represents the agents and E ⊆ V × V denotes the set of edges; i.e., if agent i and j are connected then (i, j) ∈ E implies (j, i) ∈ E. Suppose we have a collection of n independent and identically distributed (i.i.d.) data pairs zi = (ai , yi ), where ai ∈ Rp is the feature vector and yi the label or response of the i-th observation. Let Z = [z1 , z2 , · · · , zn ] ∈ Rnp be sampled from the distribution p(Z|x) where the parameter x ∈ Rd has a common prior. The goal is to sample from the posterior distribution p(x|Z) ∝ p(Z|x)p(x) by distributing Z among N agents such that Zi = {z1i , z2i , · · · , zni i } is the subset of data exclusive to agent i. Thus, Z = ∪N i ̸= j. Therefore, by the additive property of i=1 Zi and Zi ∩ Zj = ∅ for PN Pni the log-likelihood function, we have log p(Z|x) = i=1 j=1 log p(zji |x). By defining the potential function f (x) as f (x) :=

N X

fi (x),

fi (x) := −

i=1

ni X j=1

log p(zji |x) −

1 log p(x), N

(1.1)

we set our goal to sample from the convex constrained set K ⊆ Rd , with the probability density function π(x) := p(x|Z) ∝ e−f (x) , where each component function fi is processed with the data that 1

Department of Mathematics, Florida State University, Tallahassee, Florida, United States of America; [email protected] 2 Department of Mathematics, Florida State University, Tallahassee, Florida, United States of America; [email protected]

are exclusive to agent i only. The choice of the log-likelihood function depends on the type of the learning problem, e.g., Bayesian linear regression [Hoff, 2009], Bayesian logistic regression [Hoff, 2009], Bayesian deep learning [Wang and Yeung, 2020, Polson and Sokolov, 2017] and Bayesian principal component analysis [Dubey et al., 2016].

Related Work

Decentralized optimization has been extensively studied in recent years to address large-scale learning problems where centralized computation is either infeasible or undesirable due to communication, privacy, or robustness constraints. A wide variety of decentralized first-order algorithms have been proposed, including decentralized gradient descent, EXTRA-type methods, decentralized stochastic gradient descent (DSGD), decentralized accelerated proximal stochastic gradient (DAPSGD), and accelerated variants; see, for example, [Ye et al., 2020, Arjevani et al., 2020, Shi et al., 2015, Lian et al., 2017]. Motivated by stochastic optimization literature, there has been a growing interest in decentralized sampling algorithms. Decentralized Langevin-type algorithms extend classical Langevin Monte Carlo by allowing each agent to access only its local data while communicating with neighbors to approximate the global posterior. Representative works include decentralized stochastic gradient Langevin dynamics (DE-SGLD), EXTRA-SGLD, DIGing-SGLD, and their accelerated variants [Gürbüzbalaban et al., 2021, Gurbuzbalaban et al., 2024b, Bajwa et al., 2025, Yao et al., 2025]. These methods are primarily designed for unconstrained sampling on Rd and rely on consensus mechanisms to ensure asymptotic agreement across agents. Sampling from probability distributions supported on constrained domains or involving non-smooth potentials has attracted significant attention in the centralized setting. One classical approach is projected Langevin Monte Carlo (PLMC), studied in Bubeck et al. [2015, 2018], where each Langevin step is followed by a projection onto the feasible set: p xk+1 = PK xk − η∇f (xk ) + 2ηξk+1 , (2.1) with PK denoting the Euclidean projection onto a convex set K, and the dynamics (2.1) is based on the discretization of the continuous-time overdamped Langevin stochastic differential equation (SDE) with reflected boundary: √ dXt = −∇f (Xt )dt + 2dWt + ν(Xt )L(dt), (2.2) the term ν(Xt )L(dt) ensures that Xt ∈ K for every t given that X0 ∈ K. In particular, Rwhere t 0 ν(Xs )L(ds) is a bounded variation reflection process and the measure L(dt) is such that L([0, t]) is finite, L(dt) is supported on {t|Xt ∈ ∂K}. Non-asymptotic convergence rates in total variation distance were established in Bubeck et al. [2018], though with a strong dependence on the underlying dimension. Projected stochastic gradient Langevin dynamics has also been studied under weaker assumptions on the objective function. In Lamperski [2021], the author analyzed projected SGLD for possibly non-convex smooth objectives with stochastic gradients satisfying sub-Gaussian noise conditions, obtaining convergence guarantees in the 1-Wasserstein distance; see also Zheng and Lamperski [2022] for related developments. Mirror descent-based Langevin algorithms (see e.g. Hsieh et al. [2018], Chewi et al. [2020], Zhang et al. [2020], Li et al. [2022], Ahn and Chewi [2021]) form another important class for constrained sampling. Mirrored Langevin dynamics was proposed in Hsieh et al. [2018], inspired by the classical mirror descent in optimization. An alternative strategy to handling constraints is through proximal and penalty-based methods. Proximal 2

Langevin Monte Carlo was introduced in Brosse et al. [2017], where the non-smooth component of the potential is handled via Moreau–Yosida regularization. This approach avoids explicit projections and yields polynomial complexity bounds under log-concavity assumptions. Subsequent work by Salim and Richtárik [2020] further examined proximal SGLD from a primal–dual perspective. More recently, penalized Langevin Monte Carlo methods have been proposed as an alternative to projection- or proximal-based approaches. Inspired by classical penalty methods in optimization, Gurbuzbalaban et al. [2024a] introduced penalized Langevin dynamics for constrained sampling with potentially non-convex objectives and demonstrated improved dimension dependence compared to earlier methods. Related work combining constraints with replica-exchange Langevin dynamics was studied in Zheng et al. [2024], showing that reducing the effective domain diameter can significantly accelerate mixing. In addition, motivated by the acceleration properties of non-reversible Langevin dynamics on Rd Hwang et al. [1993, 2005], Futami et al. [2020], Gao et al. [2020], Hu et al. [2020], recent studies have begun to explore non-reversible dynamics for constrained sampling. In particular, Du et al. [2025] proposed skew-reflected non-reversible Langevin dynamics and established non-asymptotic convergence guarantees, showing that breaking reversibility can yield acceleration even in the presence of constraints. Further theoretical insights based on large deviations and asymptotic variance, along with extensive numerical experiments, are provided in Wang et al. [2025]. Despite these advances, the literature on decentralized sampling and constrained Langevin methods has developed largely independently. To the best of our knowledge, existing decentralized Langevin algorithms are limited to smooth, unconstrained settings, while constrained and proximal Langevin methods rely on centralized computation. In particular, there is no prior work that systematically studies decentralized proximal or penalized stochastic gradient Langevin dynamics for sampling from constrained or composite posterior distributions. Addressing this gap is the main motivation of the present work. The contributions of our paper can be summarized as follows: • We propose a new sampling algorithm that can be used to generate samples from a constrained domain in a decentralized system. To the best of our knowledge, this is the first decentralized Langevin sampling algorithm in a constrained domain. • We provide a non-asymptotic convergence analysis in the 2-Wasserstein distance for both individual agent and network consensus chain (Theorem 5). Based on this, we provide an iteration complexity for a given accuracy level (Corollary 6). • Finally, we show the efficiency of our algorithm for sampling with synthetic and real data for different sampling problems, including Bayesian linear and logistic regressions.

Background and Problem Setup

3.1

Langevin Dynamics

Langevin algorithms are one of the most widely used Markov Chain Monte Carlo (MCMC) methods in statistical learning that allow sampling from a given density π(x) of interest. The classical Langevin algorithm is based on the discretization of overdamped Langevin SDE [Dalalyan, 2017, Durmus and Moulines, 2017]: √ dX(t) = −∇f (X(t))dt + 2dWt , (3.1) 3

where f : Rd → R and Wt is a standard d−dimensional Brownian motion with W0 = 0. Under some mild assumptions on f , the diffusion (3.1) admits a unique stationary distribution with the density π(x) ∝ e−f (x) , also known as the Gibbs distribution [Pavliotis, 2014]. The implementation of this algorithm requires the discretized version of the dynamics, and Euler-Maruyama discretization is the simplest one, known as Unadjusted Langevin Algorithm (ULA) [Durmus and Moulines, 2017, 2019]: p xk+1 = xk − η∇f (xk ) + 2ηwk+1 , (3.2) where η > 0 is the step-size (or learning rate), and wk ∈ Rd is a sequence of independent and identically distributed (i.i.d.) standard Gaussian random vectors N (0, Id ). However, the discretization (3.2) does not converge to the target distribution π(x), and it introduces a bias that needs to be properly characterized to ensure performance guarantees [Dalalyan and Karagulyan, 2019]. The Unadjusted Langevin Algorithm (ULA) (3.1) also requires the computation of ∇f at each iteration, which can be computationally expensive and often impractical when the data are large and multi-dimensional. This issue can be handled efficiently using stochastic gradient, instead of the full gradient [Bottou, 2010] which results in the algorithm Stochastic Gradient Langevin Dynamics (SGLD) [Welling and Teh, 2011, Raginsky et al., 2017] given as p ˜ (xk ) + 2ηwk+1 , xk+1 = xk − η ∇f (3.3) where η > 0 is the step-size (or learning rate), wk ∈ Rd is a sequence of i.i.d. standard Gaus˜ (xk ) is an unbiased estimator of the sian random vectors N (0, Id ), and the stochastic gradient ∇f deterministic gradient with a bounded variance. When the number of data points is large, these stochastic gradients are cheaper to estimate and can be computed from a mini-batch setup.

3.2

Decentralization

We have a collection of N computational agents connected over a communication network with the communication matrix W , where W = [Wi,j ] ∈ RN ×N is symmetric, doubly stochastic matrix with the properties, / E, P Wij = Wji > 0 for i ̸= j if {i, j} ∈ E, Wij = Wji = 0 for {i, j} ∈ and Wii = 1 − j̸=i Wij > 0 for every 1 ≤ i ≤ N . Moreover, the eigenvalues of W satisfy W W 1 = λW 1 > λ2 > · · · > λN > −1, with W 1 = 1 where 1 is a vector of length N with entries equal to 1. A common approach to compute W is taking W = IN − δL, where IN is an N × N identity matrix, L is the graph Laplacian, and δ > 0 is a small number satisfying 0 < δ < λ2L [Olfati-Saber N

W et al., 2007, Chung, 1997]. We define the spectral gap 1 − ρ ∈ (0, 1) where, ρ := max{|λW 2 |, |λN |} is the second largest eigenvalue of W . (k) Let xi be the local variable of the agent i at k-th iteration. The Decentralized Stochastic Gradient Langevin Dynamics (DE-SGLD) [Gürbüzbalaban et al., 2021] consists of a weighted av(k) eraging with local variables xj of node i’s immediate neighbors j ∈ Ωi := {j : (i, j) ∈ G} as well as a stochastic gradient step over the node’s component function fi (x), i.e. p X (k) ˜ i x(k) + 2ηw(k+1) , xk+1 = Wij xj − η ∇f (3.4) i i i j∈Ωi

˜ i x(k) is the stochastic gradient, and w(k) ∈ Rd are i.i.d. standard where η > 0 is the stepsize, ∇f i i Gaussian random vectors N (0, Id ). The DE-SGLD algorithm is an unconstrained sampling algorithm for any x ∈ Rd and can be used for each agent i to sample a common distribution e−f (x) . We 4

want to introduce regularization and constraints via a proximal function so that for any x ∈ Rd , the iterations are onto a convex compact set K. We propose and discuss this technique in the next section.

Main Results

Let K ∈ Rd be a compact convex set with non-empty interior. To sample from a distribution π restricted to a compact convex set K, we cannot use the ULA, and its discretized versions in either centralized (3.3) or decentralized (3.4) setting. For constraint sampling, we need to add some regularization in the form of a proximal function. The unconstrained potential function u : Rd → (−∞, +∞] of the form u(x) := f (x) + qK (x), is associated with the probability density function (PDF) π(x) ∝ e−u(x) , and qK (x) is the indicator function given as ( +∞ if x ̸= K, qK (x) = (4.1) 0 if x ∈ K. For the constraint sampling, we use the Moreau-Yosida envelope [Rockafellar and Wets, 1998], γ qK : Rd → R+ of the form 1 1 γ 2 qK (x) = inf qK (y) + ∥x − y∥ = ∥x − PK (x)∥2 , (4.2) 2γ 2γ y∈Rd where PK (x) := arg min ∥x − y∥ denotes the projection onto K, and γ > 0 is a regularization y∈K

γ parameter. Then the constrained potential uγ : Rd → R given by uγ (x) := f (x) + qK (x), is γ

1 associated with the stationary distribution π γ (x) ∝ e−u (x) = exp −f (x) − 2γ ∥x − PK (x)∥2 , for

any x ∈ Rd .

4.1

Distance Between π γ and π

In this section, we bound the distance between π γ and π in terms of the 2-Wasserstein distance. Let us define P2 (Rd ) as the space of all Borel probability measures µ on Rd with finite second moment (with respect to the Euclidean norm). For any two Borel probability measures µ, ν in P2 (Rd ), the 2-Wasserstein distance is given as [Villani, 2008]: Z 2 W2 (µ, ν) := inf ∥x − y∥2 dγ(x, y), γ∈Π(µ,ν) Rd ×Rd

where Π(µ, ν) is the set of couplings of (µ, ν), and µ is absolutely continuous with respect to ν. We have the following lemma that upper bounds the 2-Wasserstein distance between π γ and π. Lemma 1. Fix any γ0 ∈ (0, 1/e). Then, for any 0 < γ ≤ γ0 , W2 (π γ , π) ≤ Cγ 1/8 (log(1/γ))1/8 ,

(4.3)

where C is a constant that depends on γ0 , and defined in (C.9) in the proof given in Appendix C.1.

4.2

Decentralized Proximal SGLD

We are now ready to present our main algorithm that requires the following assumptions. Throughout this paper, we assume that the following assumptions hold. Assumption 2. We assume each component function fi : Rd → R, ∈ Sµ,L is µ−strongly convex and L-smooth with L > µ. That is for any fi ∈ Sµ,L , and for every x1 , x2 ∈ Rd , µ L ∥x1 − x2 ∥2 ≤ fi (x1 ) − fi (x2 ) − ∇fi (x2 )⊤ (x1 − x2 ) ≤ ∥x1 − x2 ∥2 . 2 2 Next, denote B(x, r) as the closed ball centered at x ∈ Rd with radius r > 0 such that B(x, r) := {y ∈ Rd : ∥y − x∥ ≤ r}. We impose the following assumption on K. Assumption 3. There exist r, R > 0, r ≤ R, such that, B(0, r) ⊂ K ⊂ B(0, R). Next, we impose the following assumption on the gradient noise. (k+1) ˜ i x(k) − ∇fi x(k) are unbiased Assumption 4. The gradient noises are defined as ξi := ∇f i i i h (k+1) (k+1) 2 2 with a finite second moment. That is, E ξi Fk = 0, and E ξi ≤ σ , where Fk is the (k)

natural filtration of the iterates xi

up to and including time k.

We also use two important properties from Rockafellar and Wets [1998] and Brosse et al. [2017]. λ is convex, continuously differentiable, and γ −1 -Lipschitz, i.e., The proximal function qK λ ∇qK (x) = γ −1 (x − PK (x)) ,

λ λ ∇qK (x) − ∇qK (y) ≤ γ −1 ∥x − y∥ ,

∀ x, y ∈ Rd .

Based on the setup and assumptions, we present the Decentralized Proximal Stochastic Gradient Langevin Dynamics (DE-PSGLD) as X p 1 (k) (k) (k+1) (k+1) (k+1) (k) (k) − ηξi + 2ηwi xi − P K xi xi = Wij xj − η ∇fi xi + , (4.4) Nγ j∈Ωi

(k)

where η > 0 is the stepsize, ξi are stochastic gradient noise satisfying Assumption 4, wi are i.i.d. Gaussian random vectors N (0, Id ) independent of the stochastic gradient noise, and Ωi = {j : (i, j) ∈ G} are the neighbors of the node i.

4.3

Convergence Analysis

Next, we concatenate the local decision variables into a single vector to facilitate convergence analysis. For K ⊂ Rd , we define the product set KN := K × · · · × K ⊂ RN d , and let x(k) define the decision vector from all agents x

(k)

⊤ (k) ⊤ (k) ⊤ (k) ⊤ , x2 , . . . , xN ∈ RN d . x1

(4.5)

P 2 Nd → R Define G : RN d → R with G(x) := G(x1 , x2 , · · · , xN ) = N i=1 ∥xi − PK (xi )∥ and F : R PN with F (x) := F (x1 , x2 · · · , xN ) = i=1 fi (xi ). Then, we further define U γ (x) := F (x) + 2N1 γ G(x) 6

γ PN γ qK (xi ) γ γ d such that U γ (x) = i=1 ui (xi ), where ui (xi ) := fi (xi ) + N . Let x∗ ∈ R be the unique γ ⊤ γ ⊤ ⊤ γ ∗ N d minimizer of u , and xγ = [(x∗ ) , · · · , (x∗ ) ] is a vector in R . Then the concatenated DEPSGLD iterations can be written as follows p x(k+1) = Wx(k) − η∇U γ x(k) − ηξ (k+1) + 2ηw(k+1) , (4.6)

where ∇U γ (x) = ∇F (x) + N1γ ∇G(x), W = W ⊗ Id , and ⊗ denotes the Kronecker product, and the noise terms are defined as ⊤ ⊤ (k+1) ⊤ (k+1) ⊤ (k+1) ⊤ (k+1) ⊤ (k+1) (k+1) w = w1 , . . . , wN , ξ = ξ1 , . . . , ξN , 2 with the properties E ξ (k+1) Fk = 0 and E ξ (k+1) ≤ σ 2 N that follows from Assumption 4. Let P (k) denote the mean iterate at the k-th iteration. Since W is doubly stochastic, x̄(k) := N1 N i=1 xi the average iterates satisfy: " !# N N 1 X 1 1 X (k) (k) (k+1) (k) (k) x̄ = x̄ − η ∇fi xi + PK xi x̄ − N Nγ N i=1 i=1 p (4.7) − η ξ¯(k+1) + 2η w̄(k+1) , P P (k+1) (k+1) , and w̄(k+1) := N1 N . where ξ¯(k+1) := N1 N i=1 ξi i=1 wi Now we state the main result that provides the finite-time convergence guarantees. n o 2 1+λW 1 N Theorem 5. Assume that E x(0) < ∞, and the stepsize 0 < η < ηmax := min 2N , , Lγ Lγ Lγ +µ (k)

with Lγ := L + N2γ . Then, for every k ∈ N, DE-PSGLD updates xi   W2 L x̄(k) , π ≤ (1 − µη)k C0 + 

ρ2k −

1 − µη N

ρ2 −

1 − µη N

ηL 1 − 2Nγ ηL 1 − 2Nγ

and their average x̄(k) satisfy

k  12 ρ √   √ C1 + η C2 N

+ Cγ 1/8 (log (1/γ))1/8 ,  k  12 2k − 1 − µη 1 − ηLγ N ρ N 2N 1 X ρ  (k)  √ C1 W2 L xi , π ≤ (1 − µη)k C0 +   ηLγ N N ρ2 − 1 − µη i=1 N 1 − 2N +

ρk √ ˜ η C2 + √ C3 + η C4 + Cγ 1/8 (log (1/γ))1/8 , N

where C, C0 , C1 , C2 , C˜2 , C3 , C4 > 0 are explicit constants defined in the Appendix that depend on µ, L, γ, ρ, σ, d, N, γ0 , and π ∝ e−f (x) is the Gibbs distribution supported on K. Next, by using Theorem 5, we obtain the iteration complexity in terms of the dependence on the accuracy level ε and the dimension d.

K, π Corollary 6. For any given ε > 0 that is sufficiently small, we have W L x̄ ≤ Õ(ε), and 2 P (K) N 1 ≤ Õ(ε), where Õ hides the logarithmic dependence on ε, provided that i=1 W2 L xi , π N √ γ = O(ε8 ), η = O(ε18 /d) and K = O η1 log ε9d . In particular, the iteration complexity is given √ by K = O εd18 log ε9d . The proof of Theorem 5 and Corollary 6 are given in Appendix A and Appendix B, respectively.

Numerical Experiments

This section presents numerical experiments illustrating how the proposed Decentralized Proximal Stochastic Gradient Langevin Dynamics (DE-PSGLD) can be used for sampling from constrained posteriors that arise in data science, machine learning, and sampling problems. In all of the experiments below, we use four types of network structures, such as fully connected, circular, star, and disconnected, as depicted in Figure 1.

(a) Fully connected

(b) Circular

(d) Disconnected

Figure 1: Different types of network structures Synthetic 1-Dimensional Sampling. We start our experiment by sampling from a nonGaussian target on a compact convex set K = [−1, 1] in one dimension, d = 1, with density π(x) ∝ 2 4 e−f (x) where f (x) = x2 + x8 − x, and x ∈ K. The target is smooth and log-concave near the origin. To evaluate sampling accuracy, we compute the 2-Wasserstein distance in 1 dimension [Panaretos and Zemel, 2019] based on the discretization below: Z 1 n 2 k − 12 1 X 2 −1 −1 , W2 (µ, ν) = Fµ (u) − Fν (u) du ≈ Q(uk ) − x(k) , uk := n n 0 k=1

where µ denotes the true distribution with CDF F , ν denotes the distribution of the sampler, Q = F −1 is the true quantile of the target on [−1, 1], and x(k) are stored samples at the k-th iteration. We use the step size η = 5 × 10−4 and proximal regularizer γ = 3.3 × 10−4 , and we iterate through 300 steps. In each iteration, we generate 100 samples of the variable x for both DE-PSGLD and PSGLD algorithms, and for decentralization, we used a 30-agent mixing matrix. The top row in Figure 2 shows the 2-Wasserstein distances between the true distribution and randomly chosen three representative agents and the network mean across four networks. These plots suggest that a stronger connection yields rapid agent agreement and a steady decay of the 2-Wasserstein distance. However, a weaker connection, such as a star network, yields slow convergence. From the bottom row, we see that the centralized PSGLD converges faster than the mean 8

Figure 2: Evolution of the 2-Wasserstein distance between the target distribution and the samples generated by decentralized proximal SGLD (DE-PSGLD) and centralized Moreau-Yosida SGLD (PSGLD) across different network styles. DE-PSGLD due to the absence of communication constraints. The DE-PSGLD samples concentrate most of their mass inside the feasible region, and the mean chain is more stable than a single agent. In contrast, centralized PSGLD exhibits a wider distribution with more mass outside [−1, 1] as shown in Figure 3. Synthetic 2-Dimension: Bayesian Linear Regression. The constrained sampling has great usefulness in the Bayesian regression problem. For instance, if we set the constraint that the model parameters are bounded in the ℓ1 -ball, then the model is referred to as the Lasso regression; if we set the constraint to be the ℓ2 -ball, then the regression is referred to as the ridge regression. As a toy experiment, we consider a two-dimensional problem with synthetic data generated using the following model yi = β∗⊤ Xi + ϵi ,

ϵi ∼ N (0, 0.25),

Xi ∼ N (0, I2 ),

β∗ = [1, 1]⊤ ,

(5.1)

where β∗ is the true value of the parameters. Given the data, our goal is to sample from the posterior distribution p(β | X, y) ∝ p(X | β)p(β). The prior p(β) is defined as a uniform distribution over the ℓp -ball of radius s. Specifically, we set s to be 80% of the norm of the Ordinary Least Squares (OLS) estimate, i.e., K = {β : ∥β∥2 ≤ s} and s = 54 ∥βOLS ∥2 . It is uniformly distributed in the ℓ2 -ball centered at the origin as shown in Figure 4. To run the sampling, we simulate 10,000 observations using the model (5.1) and distribute them among N = 20 agents. Under the assumption of Gaussian noise ϵi ∼ N (0, σ 2 ), the log-likelihood corresponds to the standard least-squares loss 1 PN Pni f (β) = 2σ2 i=1 j=1 (yj − β ⊤ Xj )2 , where ni = 500 is the number of data points in agent i. Real d-Dimension: Bayesian Logistic Regression. In this set of experiments, we use UCI ML Breast Cancer Wisconsin (Diagnostic) datasets [Wolberg et al., 1993]. The data is in the form 30 represents a feature vector containing 30 of an input-output pair (Xi , yi )569 i=1 where each Xi ∈ R predictive attributes, and the binary target yi ∈ {0, 1} represents the diagnostic classification: 212 malignant vs 357 benign. The prior p(β) is taken from the ℓ2 -ball of radius s, where s is the 80% of the norm of the β parameters from unconstrained logistic regression using the Maximum Likelihood Estimation (MLE) method. In our experiment, we consider a network of 5 agents and distribute 9

Figure 3: Comparison of the sampling from the target density π(x) ∝ e−f (x) 1[−1,1] with empirical densities obtained from DE-PSGLD and PSLGD. The mean and agent 1 are from a fully connected network style.

For this set of experiments, we take a batch size of 100 and run both DE-PSGLD and PSGLD algorithms for 500 iterations with a fixed proximal regularizer γ = 5 × 10−5 and stepsize η = 5 × 10−4 . At each iteration, we generate 300 samples of the parameters. The results from the PSGLD are shown in Figure 5: the PSGLD algorithm successfully concentrates the posterior density near the boundary of the ℓ2 -ball closest to the true parameter β∗ . Since the true parameter (the red star) lies outside the constraint set K, the posterior density is “pushed” against the boundary, illus1 trating how the proximal term 2γ ∥β − PK (β)∥2 allows the sampler to respect the hard constraint while targeting the likelihood. The result from the DE-PSGLD algorithm is shown in Figure 6. We observe the impact of the network structure on the sampling quality. Networks with some form of connections, such as fully connected, circular, and star, show strong consensus. The local posterior is nearly identical to the ensemble mean (bottom row). The fourth column, i.e., the disconnected network sets the baseline in sampling. Since the agents are not connected and do not communicate, the samples are nearly uniformly distributed within the constrained region, and the mean is far from the true parameters. It shows the importance of communication and decentralization.

Figure 4: The contour plot prior distribution constrained in the 2-norm ball around the origin.

Figure 5: The contour plot of the posterior distribution of the model parameter sampled using the PSGLD algorithm.

Figure 6: The top row shows the contour plot of samples of a randomly chosen agent out of 20 agents. The bottom row shows the mean of all agents across all 4 network structures. The red star is the true location of the parameter in 2D space. the data equally among the agents. Figure 7 illustrates the progress in classification accuracy as we increase the number of iterations up to 1000. Each iteration generates 1000 samples, and uses a mini-batch size of 10, stepsize η = 0.005, proximal parameter γ = 0.16.

Conclusion

In this paper, we studied the problem of constrained sampling in a distributed system where the agents learn to sample from a log-concave target distribution π(x) ∝ e−f (x) on a constrained domain K ⊂ Rd without sharing the local data. We proposed and studied the decentralized proximal stochastic gradient Langevin dynamics (DE-PSGLD). We provided a non-asymptotic 2-Wasserstein convergence guarantee and iteration complexity for a given accuracy level ε when the target f is smooth and strongly convex, and follow some assumptions. Finally, we provide numerical experiments to show the efficiency of our algorithm for sampling in a decentralized learning system.

Acknowledgments The authors are grateful to Mert Gürbüzbalaban for helpful discussions. Mohammad Rafiqul Islam is partially supported by the grant NSF DMS-2053454. Lingjiong Zhu is partially supported by the grants NSF DMS-2053454 and DMS-2208303.

Figure 7: The plots show the accuracy versus iterations using the DE-PSGLD and PSGLD algorithms in the classification of cancer types. The top row shows the relative performances of a randomly selected agent from the DE-PSGLD algorithm and the PSGLD algorithm across four different network structures. The bottom row displays the mean accuracy over all agents from the DE-PSGLD algorithm versus the PSGLD algorithm across four network structures.

References Kwangjun Ahn and Sinho Chewi. Efficient constrained sampling via the mirror-Langevin algorithm. In Advances in Neural Information Processing Systems (NeurIPS), volume 34, pages 28405–28418. Curran Associates, Inc., 2021. Yossi Arjevani, Joan Bruna, Bugra Can, Mert Gurbuzbalaban, Stefanie Jegelka, and Hongzhou Lin. IDEAL: Inexact DEcentralized Accelerated Augmented Lagrangian Method. In Advances in Neural Information Processing Systems, volume 33, pages 20648–20659. Curran Associates, Inc., 2020. Waheed U. Bajwa, Mert Gurbuzbalaban, Mustafa Ali Kutbay, Lingjiong Zhu, and Muhammad Zulqarnain. DIGing-SGLD: Decentralized and scalable Langevin sampling over time-varying networks. arXiv:2511.12836, 2025. François Bolley and Cédric Villani. Weighted Csiszár-Kullback-pinsker inequalities and applications to transportation inequalities. Annales-Faculté des sciences Toulouse Mathematiques, 14(3):331–352, 2005. Léon Bottou. Large-scale machine learning with stochastic gradient descent. In Proceedings of COMPSTAT’2010: 19th International Conference on Computational Statistics, Paris France, August 22-27, 2010 Keynote, Invited and Contributed Papers, pages 177–186. Springer, 2010. Nicolas Brosse, Alain Durmus, Éric Moulines, and Marcelo Pereyra. Sampling from a log-concave distribution with compact support with proximal Langevin Monte Carlo. In Conference on Learning Theory, volume 65, pages 319–342. PMLR, 2017.

Sebastien Bubeck, Ronen Eldan, and Joseph Lehec. Finite-time analysis of projected Langevin Monte Carlo. In Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015. Sébastien Bubeck, Ronen Eldan, and Joseph Lehec. Sampling from a log-concave distribution with projected Langevin Monte Carlo. Discrete & Computational Geometry, 59(4):757–783, 2018. Sinho Chewi, Thibaut Le Gouic, Cheng Lu, Tyler Maunu, Philippe Rigollet, and Austin Stromme. Exponential ergodicity of mirror-Langevin diffusions. In Advances in Neural Information Processing Systems (NeurIPS), volume 33, pages 19573–19585. Curran Associates, Inc., 2020. Fan RK Chung. Spectral Graph Theory, volume 92. American Mathematical Society, 1997. Arnak S Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society Series B: Statistical Methodology, 79(3):651–676, 2017. Arnak S Dalalyan and Avetik Karagulyan. User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient. Stochastic Processes and their Applications, 129(12):5278–5311, 2019. Hengrong Du, Qi Feng, Changwei Tu, Xiaoyu Wang, and Lingjiong Zhu. Non-reversible Langevin algorithms for constrained sampling. arXiv:2501.11743, 2025. Kumar Avinava Dubey, Sashank J Reddi, Sinead A Williamson, Barnabas Poczos, Alexander J Smola, and Eric P Xing. Variance reduction in stochastic gradient Langevin dynamics. In Advances in Neural Information Processing Systems, volume 29, pages 1154–1162. Curran Associates, Inc., 2016. Alain Durmus and Eric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. Annals of Applied Probability, 27(3):1551–1587, 2017. Alain Durmus and Eric Moulines. High-dimensional Bayesian inference via the unadjusted Langevin algorithm. Bernoulli, 25(4A):2854–2882, 2019. Futoshi Futami, Iseei Sato, and Masashi Sugiyama. Accelerating the diffusion-based ensemble sampling by non-reversible dynamics. In Proceedings of the 37th International Conference on Machine Learning, volume 119, pages 3337–3347. PMLR, 2020. Xuefeng Gao, Mert Gürbüzbalaban, and Lingjiong Zhu. Breaking reversibility accelerates Langevin dynamics for global non-convex optimization. In Advances in Neural Information Processing Systems (NeurIPS), volume 33, pages 17850–17862. Curran Associates, Inc., 2020. Mert Gürbüzbalaban, Xuefeng Gao, Yuanhan Hu, and Lingjiong Zhu. Decentralized stochastic gradient Langevin dynamics and Hamiltonian Monte Carlo. Journal of Machine Learning Research, 22(239):1–69, 2021. Mert Gurbuzbalaban, Yuanhan Hu, and Lingjiong Zhu. Penalized overdamped and underdamped Langevin Monte Carlo algorithms for constrained sampling. Journal of Machine Learning Research, 25(263):1–67, 2024a. Mert Gurbuzbalaban, Mohammad Rafiqul Islam, Xiaoyu Wang, and Lingjiong Zhu. Generalized EXTRA stochastic gradient Langevin dynamics. arXiv preprint arXiv:2412.01993, 2024b. Lie He, An Bian, and Martin Jaggi. COLA: Decentralized linear learning. In Advances in Neural Information Processing Systems, volume 31, pages 4536–4546. Curran Associates, Inc., 2018. Hadrien Hendrikx, Francis Bach, and Laurent Massoulié. An accelerated decentralized stochastic proximal algorithm for finite sums. In Advances in Neural Information Processing Systems, volume 32, pages 954–964. Curran Associates, Inc., 2019.

Peter D Hoff. A First Course in Bayesian Statistical Methods, volume 580. Springer, 2009. Ya-Ping Hsieh, Ali Kavis, Paul Rolland, and Volkan Cevher. Mirrored Langevin dynamics. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. Yuanhan Hu, Xiaoyu Wang, Xuefeng Gao, Gürbüzbalaban, and Lingjiong Zhu. Non-convex stochastic optimization via non-reversible stochastic gradient Langevin dynamics. arXiv:2004.02823, 2020. Chii-Ruey Hwang, Shu-Yin Hwang-Ma, and Shuenn-Jyi Sheu. Accelerating Gaussian diffusions. Annals of Applied Probability, 3:897–913, 1993. Chii-Ruey Hwang, Shu-Yin Hwang-Ma, and Shuenn-Jyi Sheu. Accelerating diffusions. Annals of Applied Probability, 15:1433–1444, 2005. Solomon Kullback and Richard A Leibler. On information and sufficiency. The Annals of Mathematical Statistics, 22(1):79–86, 1951. Andrew Lamperski. Projected stochastic gradient Langevin algorithms for constrained sampling and nonconvex learning. In Conference on Learning Theory, volume 134, pages 2891–2937. PMLR, 2021. Ruilin Li, Molei Tao, Santosh S. Vempala, and Andre Wibisono. The mirror Langevin algorithm converges with vanishing bias. In Sanjoy Dasgupta and Nika Haghtalab, editors, Proceedings of The 33rd International Conference on Algorithmic Learning Theory, volume 167, pages 718–742. PMLR, 2022. Xiangru Lian, Ce Zhang, Huan Zhang, Cho-Jui Hsieh, Wei Zhang, and Ji Liu. Can decentralized algorithms outperform centralized algorithms? a case study for decentralized parallel stochastic gradient descent. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. Reza Olfati-Saber, J Alex Fax, and Richard M Murray. Consensus and cooperation in networked multi-agent systems. Proceedings of the IEEE, 95(1):215–233, 2007. Victor M Panaretos and Yoav Zemel. Statistical aspects of Wasserstein distances. Annual Review of Statistics and Its Application, 6(1):405–431, 2019. Grigorios A Pavliotis. Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations, volume 60 of Texts in Applied Mathematics. Springer, 2014. Nicholas G. Polson and Vadim Sokolov. Deep learning: A Bayesian perspective. Bayesian Analysis, 12(4): 1275–1304, 12 2017. Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky. Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis. In Proceedings of the 2017 Conference on Learning Theory, volume 65, pages 1674–1703. PMLR, 2017. R Tyrrell Rockafellar and Roger JB Wets. Variational Analysis. Springer, 1998. Adil Salim and Peter Richtárik. Primal dual interpretation of the proximal stochastic gradient Langevin algorithm. In Advances in Neural Information Processing Systems (NeurIPS), volume 33, pages 3786–3796. Curran Associates, Inc., 2020. Wei Shi, Qing Ling, Gang Wu, and Wotao Yin. Extra: An exact first-order algorithm for decentralized consensus optimization. SIAM Journal on Optimization, 25(2):944–966, 2015. Cédric Villani. Optimal Transport: Old and New, volume 338. Springer, 2008. Hao Wang and Dit-Yan Yeung. A survey on Bayesian deep learning. ACM Computing Surveys, 53(5):1–37, 2020.

Yingli Wang, Changwei Tu, Xiaoyu Wang, and Lingjiong Zhu. Accelerated constrained sampling: A large deviations approach. arXiv:2506.07816, 2025. Max Welling and Yee W Teh. Bayesian learning via stochastic gradient Langevin dynamics. In Proceedings of the 28th International Conference on Machine Learning (ICML-11), pages 681–688. Citeseer, 2011. William Wolberg, Olvi Mangasarian, Nick Street, and W. Street. Breast cancer Wisconsin (Diagnostic). UCI Machine Learning Repository, 1993. https://doi.org/10.24432/C5DW2B. Nian Yao, Pervez Ali, Xihua Tao, and Lingjiong Zhu. Accelerating Langevin Monte Carlo sampling: A large deviations analysis. arXiv preprint arXiv:2503.19066, 2025. Haishan Ye, Ziang Zhou, Luo Luo, and Tong Zhang. Decentralized accelerated proximal gradient descent. In Advances in Neural Information Processing Systems, volume 33, pages 18308–18317. Curran Associates, Inc., 2020. Kun Yuan, Qing Ling, and Wotao Yin. On the convergence of decentralized gradient descent. SIAM Journal on Optimization, 26(3):1835–1854, 2016. Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, and Marcelo Pereyra. Wasserstein control of mirror Langevin Monte Carlo. In Conference on Learning Theory, volume 125, pages 3814–3841. PMLR, 2020. Haoyang Zheng, Hengrong Du, Qi Feng, Wei Deng, and Guang Lin. Constrained exploration via reflected replica exchange stochastic gradient Langevin dynamics. In Proceedings of the 41st International Conference on Machine Learning, volume 235, pages 61321–61348. PMLR, 2024. Yuping Zheng and Andrew Lamperski. Constrained Langevin algorithms with L-mixing external random variables. In Advances in Neural Information Processing Systems (NeurIPS), volume 35, pages 20511– 20521. Curran Associates, Inc., 2022.

Proof of Theorem 5

To prove Theorem 5, we split the arguments into the following steps: (k)

1. First, we compute an upper bound on the L2 distance between the iterate xi average iterate x̄(k) ;

and their

2. Next, we compute an upper bound on the L2 distance between the average iterate x̄(k) and the iterate of the centralized proximal Langevin algorithm (PLA), also known as Moreau–Yosida Unadjusted Langevin Algorithm (MYULA) in the literature [Salim and Richtárik, 2020]; 3. Finally, we compute an upper bound on the W2 distance between the law of x(k) and the Gibbs distribution π. The whole process breaks down into the following triangle inequality N

i=1

1 X 1 X (k) (k) (k) W2 L xi ,π ≤ W L xi , L x̄ + W L x̄(k) , π , N N

(A.1)

W2 L x̄(k) , π ≤ W2 L x̄(k) , L x̃(k) + W2 L x̃(k) , π γ + W2 (π γ , π) .

(A.2)

where

Here, equation (A.1) represents the average of the 2-Wasserstein distance between the distribution of the iterates of each agent and the target distribution, and equation (A.2) represents the 2Wasserstein distance between the distribution of the average iterates from all agents and the target distribution. The step-by-step proof is given in the subsequent sections.

A.1

(k)

Uniform L2 bounds between xi

and their average x̄k

We first start with computing a uniform L2 bound on the gradient ∇U γ x(k) using the following lemma. The idea is to consider the concatenated DE-PSGLD iterates in equation (4.6) as the decentralized gradient descent (DGD) subject to the proximal regularizer, stochastic gradient, and Gaussian noise. Lemma 7. Under the assumptions of Theorem 5, 2 E ∇U γ x(k) ≤ Dγ2 ,

for any k ∈ N,

where Dγ2 := 4L2γ 1 − µη 1 + λW N − ηLγ + 4 ∇U

x(0) − x∗γ

+ 8L2γ

2L2γ (ησ 2 N + 2dN ) 2 , (x∗γ ) + µ 1 + λW N − ηLγ

C12 η 2 N (1 − ρ)2

(A.3)

with Lγ := L + N2γ , and the constant C1 is defined in equation (C.13), 1 − ρ is the spectral gap, h i ⊤ ⊤ ⊤ x∗γ = (xγ∗ ) , · · · , (xγ∗ ) ∈ RN d with xγ∗ ∈ Rd being the unique minimizer of uγ . 16

The proof of Lemma 7 is given in Appendix C.2. This bound provides a useful tool to bound (k) the uniform L2 bound between the iterates xi and the mean iterates x̄(k) , which is given in the following lemma Lemma 8. Under the assumptions of Theorem 5, N 2 2 4η 2 Dγ2 4ρ2k 1 X 4η 2 σ 2 8ηd (k) (k) (0) ≤ + E xi − x̄ E x + + . N N N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) i=1

The proof of Lemma 8 is given in Appendix C.3. We use the Cauchy-Schwarz inequality in Lemma 8 to upper bound the first component of A.1, N

1 X (k) W2 L xi , L x̄(k) N i=1

≤

N 1 X

(k) 2 (k) W2 L xi , L x̄

i=1 N

≤

! 12

2 1 X (k) E xi − x̄(k) N

! 12

i=1

!1 2 2 4η 2 Dγ2 4ρ2k 4η 2 σ 2 8ηd (0) ≤ E x + + + N N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) √ 1 2 2 2ηDγ 2ησ 2ρk 2 2η d (0) √ +p E x ≤√ +p . + N (1 − ρ) N 1 − ρ2 1 − ρ2

A.2

(A.4)

L2 Distance between the mean x̄(k) and the centralized proximal Langevin (PLA) x̃(k)

When the data is distributed among the agents, the learning process introduces errors between the average of the gradients and the gradient of the average across all agents. Let us define the error term at the (k + 1)-th iteration as follows: " # " # N N X 1 X 1 (k) (k) Ek+1 := ∇fi x̄(k) − ∇fi xi + 2 PK x̄(k) − PK xi . (A.5) N N γ i=1

i=1

Since the gradients ∇fi ’s are Lipschitz based on our Assumption 2 and the projection is 1-Lipschitz, we can find an L2 bound for the error defined in (A.5) presented in the following lemma. Lemma 9. Under the assumptions of Theorem 5, for any k ∈ N, we have E ∥Ek+1 ∥2 ≤

2 8L2γ η 2 Dγ2 16L2γ ηd 8L2γ η 2 σ 2 8L2γ ρ2k E x(0) + + + . N N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 )

The proof of Lemma 9 is given in Appendix C.4. Next, let us define the mean-chain potential ûγ (x) :=

1 1 γ f (x) + qK (x), N N 17

such that ∇ûγ (x) =

1 1 ∇f (x) + (x − PK (x)) . N Nγ

Therefore, the mean iterates from (4.7) can be rewritten as follows: p x̄(k+1) = x̄(k) − η∇ûγ x̄(k) + ηEk+1 − η ξ¯(k+1) + 2η w̄(k+1) ,

(A.6)

P (k+1) (k+1) where Ek+1 is defined in (A.5), ξ¯(k+1) := N1 N with ξi being the gradient noise satisfyi=1 ξi (k+1) 1 PN (k+1) −1 ing Assumption 4, and w̄ := N i=1 wi ∼ N (0, N Id ). Now, let us define the centralized proximal Langevin algorithm (PLA) by p x̃(k+1) = x̃(k) − η∇ûγ x̃(k) + 2N −1 η z (k+1) , (A.7) √ where z (k+1) ∼ N (0, Id ), and coupled noise z (k+1) = N w̄(k+1) . This is the Euler-Maruyama discretization of the continuous-time overdamped Langevin diffusion √ dXt := −∇ûγ (Xt )dt + 2N −1 dWt , (A.8) where Wt is a standard d-dimensional Brownian motion. We are interested in upper-bounding the L2 distance between the mean iterate x̄(k) and the centralized proximal Langevin algorithm (PLA) iterate x̃(k) . We have the following lemma. Lemma 10. Under the assumptions of Theorem 5, for every k ∈ N, E x̄(k) − x̃(k) 

ηL (1 + Nγ )2



2 2 8L2γ ησ 2 16L2γ d η   8Lγ ηDγ + ≤ η + +  ηLγ µ N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) ηLγ 2 µ2 1 − 2 N 1 − 2N 2N N k 2k − 1 − µη 1 − ηLγ ρ 2 2 8L2γ ρ2 N 2N ησ + · + E x(0) . ηLγ ηL N µ 1 − 2Nγ ρ2 − 1 − µη N 1 − 2N

The proof of Lemma 10 is given in Appendix C.5.

A.3

W2 Distance between the centralized PLA and the Gibbs distribution π γ

The bounds in terms of the W2 distance between the law of the Euler-Maruyama discretization of x̃(k) from the centralized Moreau–Yosida regularized Langevin dynamics and the Gibbs distribution can be derived using the results from Dalalyan and Karagulyan [2019] with some adjustment in smoothness and convexity parameters. Recall the mean-chain potential ûγ (x) :=

1 1 γ f (x) + qK (x). N N

Note that ûγ (x) is µ/N -strongly convex and Lγ /N -smooth. Therefore, based on Theorem 4 from Dalalyan and Karagulyan [2019], we have the following lemma. 18

], we have Lemma 11. For any η ∈ (0, L2N γ +µ 1.65L p γ W2 L x̃(k) , π γ ≤ (1 − µη)k W2 L x̃(0) , π γ + ηd. Nµ Next, we need to bound the distance between the minimizer of ûγ (x) and the Gibbs distribution πγ . Lemma 12. Let x̂∗γ be the unique minimizer of ûγ (x). Then we have EX∼πγ X − x̂∗γ

≤

2d . µ

(A.9) N

The proof of Lemma 12 is given in Appendix C.6. Setting x̃(0) =

1 X (0) xi and using (A.9), N i=0

we have 2 1/2 2 1/2 W2 L x̃(0) , π γ ≤ E x̃(0) − x̂∗γ + EX∼πγ X − x̂∗γ p 2 1/2 (0) ∗ ≤ E x̃ − x̂γ + 2dµ−1 . Therefore, it follows from Lemma 11 that W2 L x̃(k) , π γ ≤ (1 − µη)k

A.4

! p 2 1/2 1.65Lγ p (0) ∗ ηd. E x̃ − x̂γ + 2dµ−1 + µ

(A.10)

Proof of Theorem 5

To complete the proof, we follow from Lemma 10 and obtain that W2 L x̄(k) , L x̃(k) 1 2 2 (k) (k) ≤ E x̄ − x̃ 

ηL (1 + Nγ )2

1 2

η √   + ≤ η 2  · ηLγ 2 µ ηL µ 1− γ 2 N 1 − 2N N

8L2γ ηDγ2 16L2γ d 8L2γ ησ 2 + + N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 )

k  12 √ 1 2k − 1 − µη 1 − ηLγ ρ 2 2 ησ N 2N 2 2ρLγ   (0)  · √ +r E x . + 2 ηLγ µη N ηLγ ρ − 1 − 1 − N 2N µ 1− 

√

(A.11)

Now we provide the final form of (A.2) using the results from (4.3), (A.10), and (A.11) as follows: W2 L x̄(k) , π  1 2 !1 ηLγ 2 8L2γ ηDγ2 8L2γ ησ 2 16L2γ d 2 (1 + N )  η √  + ≤ η + +  · ηLγ µ N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) ηLγ 2 µ2 1 − 1 − 2N N 2N N2  k  12 ηLγ √ µη 1 2k √ 2 2 ησ  ρ − 1 − N 1 − 2N  2 2ρLγ (0) E x +r  · √ + 2 ηLγ µη N ηLγ ρ − 1 − 1 − N 2N µ 1 − 2N ! 1 p 8 2 1/2 1 1.65Lγ p 1 k (0) ∗ −1 + 2dµ + (1 − µη) E x̃ − x̂γ + ηd + Cγ 8 log µ γ  k  12 2k − 1 − µη 1 − ηLγ ρ N 2N ρ √   √ C1 + ηC2 + Cγ 1/8 (log (1/γ))1/8 , (A.12) = (1 − µη)k C0 +   ηLγ N ρ2 − 1 − µη N 1 − 2N with C0 :=

1 √ 2 2 (0) C1 := 2 2Lγ E x ,

! p 2 1/2 E x̃(0) − x̂∗γ + 2dµ−1 ,

and C2 : =

1.65Lγ √ σ d+ r µ ηL µ 1 − 2Nγ 

ηL (1 + Nγ )2

1

η   + + 2  · ηLγ 2 µ ηL µ 1− γ 2 N 1 − 2N N

8L2γ ηDγ2 8L2γ ησ 2 16L2γ d + + N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 )

where x̂∗γ is the unique minimizer of ûγ , x̃(0) = N1 (k)

(0) (k) x , D is defined in (A.3), L x γ i=0 i i

denotes the distribution of xi , and π is the Gibbs distribution with probability density function proportional to e−f (x) . To obtain the final form of (A.1), we use (A.4) and (A.12), such that we

can compute that N 1 X (k) ,π W2 L xi N i=1 √ 1 2 2 2ηDγ 2 2η d 2ρk 2ησ (0) √ +p E x + +p ≤√ N (1 − ρ) N 1 − ρ2 1 − ρ2 1  k  2 2k − 1 − µη 1 − ηLγ ρ N 2N ρ √   √ C1 + ηC2 + Cγ 1/8 (log (1/γ))1/8 + (1 − µη)k C0 +   ηL γ N ρ2 − 1 − µη N 1 − 2N

  = (1 − µη)k C0 + 

ρ2k −

1 − µη N

ρ2 −

1 − µη N

k  12

ηL 1 − 2Nγ ηL 1 − 2Nγ

ρ   √ C1 N

ρk √ ˜ η C2 + √ C3 + ηC4 + Cγ 1/8 (log (1/γ))1/8 , N 1 √ 2 2 2D with C˜2 := C2 + √2 2 2 , C3 = 2 E x(0) , and C4 := (1−ρ)γ√N + √ 2σ 2 . This completes the proof +

1−ρ

of Theorem 5.

Proof of Corollary 6

Proof. From Theorem 5, for the average iterate x̄(k) , we have ρ2k − (Aγ (η))k ρ2 − Aγ (η)

!1/2

ρ √ C1 N √ 1/8 , + ηC2 + Cγ 1/8 log(1/γ)

W2 L x̄(k) , π ≤ (1 − µη)k C0 +

(B.1)

and N

1 X (k) W2 L xi , π ≤ (1 − µη)k C0 + N i=1

ρ2k − (Aγ (η))k ρ2 − Aγ (η)

!1/2

ρk ρ √ √ C1 + η C˜2 + √ C3 N N

+ η C4 + Cγ 1/8 (log (1/γ))1/8 ,

ηL

(B.2)

γ 1 1 where Aγ (η) := 1 − µη N 1 − 2N . Since η ≤ Lγ +µ ≤ Lγ , we have ηLγ ≤ 1, so that ηLγ µη Aγ (η) = 1 − 1− ≥ 1 − µη. N 2N 2

2 2 Now, choose η ≤ 1−ρ µ . Then we have 1 − µη ≥ ρ . This implies Aγ (η) ≥ ρ , and thus we can rewrite the second term in (B.1) as !1/2 (Aγ (η))k/2 (Aγ (η))k − ρ2k p ≤ . Aγ (η) − ρ2 Aγ (η) − ρ2

Therefore, 1/8 ρ (Aγ (η))k/2 √ C1 + η C2 + C γ 1/8 log(1/γ) , W2 L(x̄(k) ), π ≤ (1 − µη)k C0 + p N (Aγ (η) − ρ2 ) and N

ρ (Aγ (η))k/2 1 X ρk √ (k) k p W2 L xi , π ≤ (1 − µη) C0 + C1 + η C˜2 + √ C3 N N (Aγ (η) − ρ2 ) N i=1

+ η C4 + Cγ 1/8 (log (1/γ))1/8 . (B.3) n o 1+λW 1 N Next, let us recall that η < min 2N , , and Lγ = L + N2γ = O(γ −1 ) such that Lγ Lγ Lγ +µ η = O(γ −1 ).

(B.4)

Then, let us recall that C :=

1 2

(A∗ (γ)) +

A∗ (γ) 2

1 ! 4

sup Ĉ,

(B.5)

0<γ≤γ0

where Ĉ is defined in (C.8) and A∗ =

d VR

√ 2α̃ e−m∗ rZK

!d−1 p 2α̃γ0 log(1/γ0 ) J∗ 1+ + , r ZK

(B.6)

d/2

where VR := Γ(πd +1) Rd . We also recall that 2

C0 :=

! p 2 1/2 E x̃(0) − x̂∗γ + 2dµ−1 ,

1 √ 2 2 (0) , C1 := 2 2Lγ E x

and C2 : =

1.65Lγ √ σ d+ r µ ηL N µ 1 − 2γ 1



(1 + ηLγ )2

  + + 2  · ηLγ ηL γ µ 1− 2 µ2 1 −

8L2γ ηDγ2 8L2γ ησ 2 16L2γ d + + N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 )

!1 2

P (0) where x̂∗γ is the unique minimizer of ûγ , x̃(0) = N1 N i=0 xi , Dγ is defined in (A.3). Let us spell out the dependence of C, C0 , C1 , C2 on γ and d. First, by using Stirling’s approximation, we can derive from (B.6) that ! √ √ 2α̃ π d/2 Rd de r d γ0 log(1/γ0 ) = O(1), A∗ = O d log(d/2)− d2 2 e 22

(B.7)

as d → ∞. By plugging this into (B.5), we get C = O(1).

(B.8)

Next, since x̂∗γ is the unique minimizer of ûγ and as γ → 0, x̂∗γ will converge to the closure of the set K, and this implies that √ C0 = O d , as γ → 0. (B.9) After that, since Lγ = L + N2γ , it follows that C1 = O(1/γ),

as γ → 0.

(B.10)

Now, in order to find how C2 behaves as γ → 0, we first analyze Dγ where Dγ2 := 4L2γ 1 − µη 1 + λW N − ηLγ + 4 ∇U

x(0) − x∗η,γ

+ 8L2γ

C12 η 2 N (1 − ρ)2

2L2γ (ησ 2 N + 2dN ) 2 , (x∗γ ) + µ 1 + λW N − ηLγ

where we recall from (C.13) that v u N X 2(Lγ + µ) u t2L C1 := 1 + uγi (0) − uγ,∗ , γ i µ

with

uγ,∗ := min uγi (x). i

i=1

x∈Rd

(B.11)

(B.12)

Here, for the first term in (B.11), we have 4L2γ 1 − µη 1 + λW N − ηLγ

x(0) − x∗η,γ

= O(L2γ ) = O(γ −2 ).

For the second term in (B.11), we have 8L2γ

C12 η 2 N = O(L2γ L3γ η 2 ) = O(γ 2α−5 ), (1 − ρ)2

by taking η = O(γ α ) for some α ≥ 1 to be chosen later which satisfies the constraint (B.4), where 3/2 we used the fact that C1 = O(Lγ ) = O(1/γ 3/2 ) from (B.12) since 0 ≤ uγi (0) − uγ,∗ ≤ uγi (0) = fi (0) = O(1), i as γ → 0, where the first equality above is due to 0 ∈ K. For the last term in (B.11), we have 2L2γ (ησ 2 N + 2dN ) = O(dL2γ ) = O(dγ −2 ). µ 1 + λW − ηL γ N Therefore, Dγ2 = O dγ −2 + dγ

2α−5

 √ O dγ α−5/2 , √ Dγ = O dγ −1 , 23

1 ≤ α < 3/2, α ≥ 3/2.

For α ≥ 1, we have ηLγ = O γ α−1 . Therefore, L2γ η Dγ2 + L2γ η + dL2γ = O γ −2 γ α Dγ2 + γ −2 γ α + γ −2 = O γ α−2 Dγ2 + γ α−2 + dγ −2 . Now for 1 ≤ α < 3/2, we have Dγ2 = O(dγ 2α−5 ). Then, √ √ 3α−7 √ 3α−7 C2 = O dγ −1 + O(1) + O dγ 2 =O dγ 2 . If 3/2 ≤ α < 2, then C2 = O

√ √ α−4 √ α−4 dγ −1 + O(1) + O dγ 2 dγ 2 . =O

√ dγ −1 . Thus we can summarize the behavior of C2 as γ → 0 If α ≥ 2, then we have C2 = O depending on some α ≥ 1 as follows  √ 3α−7  dγ 2 , 1 ≤ α < 3/2, O    √ α−4 dγ 2 , 3/2 ≤ α < 2, C2 = O  √   O dγ −1 , α ≥ 2. By taking η = O(γ α /d), we get  4α−7  , 1 ≤ α < 3/2, O γ 2  √ ηC2 = O γ α−2 , 3/2 ≤ α < 2,   O γ α2 −1 , α ≥ 2.

(B.13)

Now for a given ε > 0, we choose γ and η such that √

ε η C2 ≤ , 6

ε Cγ 1/8 (log(1/γ))1/8 ≤ . 6 √ √ Thus, ε → 0, we require γ → 0 and ηC2 → 0. From (B.13), 2 → 0 as γ → 0 if and only if α ηC √ −1 α , and and we need to choose γ α > 2. For α > 2, with η = O(γ /d), we have η C2 = O γ 2 α and η such that η = O(γ /d) with α > 2 and √ α ε ε O dγ 2 −1 ≤ , and Cγ 1/8 (log(1/γ))1/8 ≤ . 6 6 α Then matching γ 2 −1 with γ 1/8 we obtain α = 94 . From these, we choose γ = O ε8 to obtain and

Cγ 1/8 (log(1/γ))1/8 ≤ Õ(ε),

(B.14)

where Õ hides the logarithmic dependence on ε. and we further choose η = O(γ 9/4 /d) = O(ε18 /d) such that √ ηC2 + Cγ 1/8 (log(1/γ))1/8 ≤ Õ(ε), where Õ hides the logarithmic dependence on ε. Also, recall that, √ 1 2 2 2Dγ 2 2 2σ (0) √ +p C˜2 := C2 + p , C3 := 2 E x , C4 := , (1 − ρ) N 1 − ρ2 1 − ρ2 24

√ √ where C˜2 = O( dγ −1 ), C3 = O(1), and C4 = O( dγ −1 ) as γ → 0. Accordingly, we get √ ˜ η C2 + ηC4 + Cγ 1/8 (log(1/γ))1/8 ≤ Õ(ε), where Õ hides the logarithmic dependence on ε. Now, given such η, γ > 0, it remains to choose K so that ρ (Aγ (η))K/2 ε p C1 ≤ , 2 3 N (Aγ (η) − ρ )

ε (1 − µη)K C0 ≤ , 3

(B.15)

which implies that W2 L x̄(k) , π ≤ Õ(ε), and

(B.16)

1 X (k) W2 L xi , π ≤ Õ(ε), N

(B.17)

i=1

where Õ√hides the logarithmic dependence on ε. From the first inequality in (B.15), by using C0 = O( d), we have √ !! log(3C0 /ε) log(3C0 /ε) 1 d K≥ ≍ =O log . | log(1 − µη)| µη η ε From the second inequality in (B.15), we have ρ (Aγ (η))K/2 ε p C1 ≤ , 3 N (Aγ (η) − ρ2 ) which implies that K/2

(Aγ (η))

p N (Aγ (η) − ρ2 ) ≤ . 3ρ C1 ε

Therefore, by using C1 = O(1/γ) = O(1/ε8 ), we get 3ρ C ! 1 2 log √ ε N (Aγ (η)−ρ2 ) 2 3ρ C1 1 1 p ≍ log =O log 9 . K≥ 2 | log(Aγ (η))| µη η ε ε N (Aγ (η) − ρ ) Finally, since η = O ε18 /d , we conclude that W2 L x̄(k) , π ≤ Õ(ε), and

1 X (k) W2 L xi , π ≤ Õ(ε), N i=1

provided that K=O

d log ε18

√ !! d , 9 ε

which gives the iteration complexity. This completes the proof. 25

Proofs of the Technical Lemmas

C.1

Proof of Lemma 1

Proof. We first compute an upper bound on the Kullback-Leibler (KL) divergence between π γ and π, and then provide an upper bound in terms of the 2-Wasserstein distance using the weighted Csiszár–Kullback–Pinsker inequality by Bolley and Villani [2005] (see e.g. Lemma B.1 in Gurbuzbalaban et al. [2024a]). First, we recall that if P and Q are probability distributions of continuous random variables with densities p and q on Rd such that P is absolutely continuous with respect to Q, the Kullback and Leibler (KL) divergence [Kullback and Leibler, 1951] between P and Q is defined as Z p(x) p(x) log D(P ∥ Q) = dx. q(x) Rd γ γ By definition qK (x) = 0 for any x ∈ K, and qK (x) > 0 for any x ∈ / K. Assume that e−f is integrable over K. Define the normalizing constants Z Z 1 −f (y)− 2γ ∥y−PK (y)∥2 −f (y) ZK := e dy, Zγ := e dy. Rd

Then, the corresponding densities are e−f (x) 1x∈K π(x) = R −f (y) ; dy Ke

π (x) = R

1 −f (y)− 2γ ∥y−PK (y)∥2

1 ∥y−PK (y)∥2 −f (y)− 2γ

. dy

Thus, the KL divergence between π and π γ can be computed as Z π(x) γ D(π∥π ) = log π(x)dx π γ (x) Rd ! R 1 Z −f (y)− 2γ ∥y−PK (y)∥2 e dy e−f (x) 1x∈K e−f (x) 1x∈K d R = log R −f (y) · R dx. 1 −f (y) dy −f (y)− 2γ ∥y−PK (y)∥2 dy K e Ke Ke Note that on K, we have PK (x) = x. Therefore, D(π∥π γ ) =

Rd e

Z log

1 −f (y)− 2γ ∥y−PK (y)∥2

R = log

−f (y) dy

1 ∥y−PK (y)∥2 −f (y)− 2γ

R R

= log

e−f (x) 1x∈K R dx −f (y) dy K Ke

−f (y) dy Ke

1 −f (y)− 2γ ∥y−PK (y)∥2 dy Rd e R −f (y) dy Ke

e−f (x) 1x∈K R dx −f (y) dy Ke

 R 1 −f (y)− 2γ ∥y−PK (y)∥2 −f (y) dy + e e dy d K R \K . R = log  −f (y) dy Ke R

Define T (x) = g(γK (x)) = (γK (x))2 , where γK (x) = ∥x − PK (x)∥ is the distance of the point x to the set K. Therefore R 1 T (y) −f (y)− 2γ dy Rd \K e γ R . (C.1) D(π∥π ) ≤ −f (y) dy Ke Notice that for any ω > 0, Z 1 −f (y)− 2γ T (y) e dy d R \K Z Z 1 T (y) −f (y)− 2γ dy + e = ≤ y ∈ R \ K : T (y) ≤ ω e

y∈Rd \K:T (y)>ω

y∈Rd \K:T (y)≤ω d

1 −f (y)− 2γ T (y)

− inf y∈Rd \K:T (y)≤ω f (y)

ω − 2γ

Z e

1 −f (y)− 2γ T (y)

dy.

(C.2)

Rd \K

We now state a result from Gurbuzbalaban et al. [2024a] (Lemma 2.4, without a proof) to provide an upper bound on the Lebesgue measure of the set K. Lemma 13. Assume the constraint set K is a bounded set containing an open ball with radius r > 0. Let T (x) = g(γK (x)) and g : R≥0 → R≥0 is a strictly increasing function with g(0) = 0 and g(x) → ∞ as x → ∞. Then for τ > 0, ! d g −1 (τ ) d x ∈ R \ K : T (x) ≤ τ ≤ 1+ − 1 |K|, r where | · | denotes the Lebesgue measure and g −1 is the inverse function of g. √ Using the result from Lemma 13 with g(x) = x2 such that g −1 (x) = x and taking ω = 2α̃γ log(1/γ) for any α̃ ≥ 12 , we obtain from (C.2): Z 1 − inf y∈Rd \K:T (y)≤2α̃γ log(1/γ) f (y) −f (y)− 2γ T (y) e dy ≤ y ∈ Rd \ K : T (y) ≤ 2α̃γ log(1/γ) e Rd \K Z 1 −f (y)− 2γ T (y) α̃ dy +γ e d R \K   !d p 2α̃γ log(1/γ) − inf y∈Rd \K:T (y)≤2α̃γ log(1/γ) f (y) − 1 |K|e ≤ 1+ r Z 1 −f (y)− 2γ T (y) α̃ +γ e dy d R \K   !d p 2α̃γ log(1/γ) − inf y∈Rd \K:T (y)≤2α̃γ log(1/γ) f (y) ≤ 1+ − 1  VR e r Z 1 −f (y)− 2γ T (y) + γ α̃ e dy, (C.3) Rd \K

d/2

where VR := Γ πd +1 Rd . Since K is contained in an Euclidean ball with radius R, we have |K| ≤ VR . (2 ) Thus, we obtain inequality (C.3). Therefore, the upper bound in (C.1) can be written as:   !d p 2α̃γ log(1/γ) J(γ) VR e−m(γ) + γ α̃ , (C.4) D(π∥π γ ) ≤  1 + − 1 r ZK ZK where

Z m(γ) :=

inf

y ∈K:T / (y)≤2α̃γ log(1/γ)

f (y),

1 −f (y)− 2γ T (y)

J(γ) :=

dy.

Rd \K

Take γ0 ∈ (0, 1/e), and set m∗ :=

inf

0<γ≤γ0

m(γ),

J∗ := sup J(γ). 0<γ≤γ0

Then for every 0 < γ ≤ γ0 , m(γ) ≥ m∗ ,

J(γ) ≤ J∗ .

Now, the first term of (C.4) can be bounded by using the inequality (1 + u)d − 1 ≤ du(1 + u)d−1 ,

u ≥ 0.

Therefore, 



2α̃γ log(1/γ) VR e−m(γ) − 1 r ZK ! !d−1 p p 2α̃γ log(1/γ) 2α̃γ log(1/γ) VR −m(γ) e ≤d 1+ r r ZK !d−1 p √ p 2α̃γ log(1/γ) d VR 2α̃ e−m(γ) γ log(1/γ) = 1+ rZK r p = A1 (γ) γ log(1/γ),

 1+

where A1 (γ) :=

d VR

√

2α̃ e−m(γ) rZK

p 1+

2α̃γ log(1/γ) r

!d−1 ,

d/2

where VR := Γ πd +1 Rd . (2 ) p √ Since we have α̃ ≥ 12 , for 0 < γ ≤ γ0 , we have γ α̃ ≤ γ ≤ γ log(1/γ). Thus, the second part of (C.4) can be written as: γ α̃

J(γ) J∗ J∗ p ≤ γ α̃ ≤ γ log(1/γ). ZK ZK ZK

J∗ A1 (γ) + ZK

Now, from (C.4), γ

D (π∥π ) ≤

p p γ log(1/γ) = A∗ γ log(1/γ), 28

(C.5)

where A∗ =

d VR

√ 2α̃ e−m∗ rZK

!d−1 p 2α̃γ0 log(1/γ0 ) J∗ 1+ , + r ZK

(C.6)

where we used the fact that the map x 7→ x log(1/x) is increasing on x ∈ (0, 1/e). Now we are ready to use the weighted Csiszár–Kullback–Pinsker inequality by Bolley and Villani [2005] (see e.g. Lemma B.1 in Gurbuzbalaban et al. [2024a]) to provide a bound in the 2-Wasserstein distance: 1 ! D (π∥π γ ) 4 γ γ 12 W2 (π, π ) ≤ Ĉ (D (π∥π )) + , (C.7) 2 where Ĉ = 2

inf

x̂∈Rd ,α̃≥ 21

1 α̃

3 + log 2

Z e

−α̃∥x−x̂∥2

dπ

(C.8)

By plugging the results from (C.5) into (C.7), we have 1 1 1 ! 4 4 8 1 1 A 1 ∗ + γ log W2 (π, π γ ) ≤ Ĉ (A∗ ) 2 γ log γ 2 γ 1 1 1 ! 8 8 1 1 A∗ 4 1 ≤ Ĉ (A∗ ) 2 γ log + γ log γ 2 γ ! 1 1 8 1 A∗ 4 1 2 = Ĉ (A∗ ) + γ log . 2 γ Let C :=

1 2

(A∗ ) +

A∗ 2

1 ! 4

sup Ĉ,

(C.9)

0<γ≤γ0

where A∗ is defined in (C.6) and Ĉ is defined in (C.8). Thus, 1 8 1 γ . W2 (π, π ) ≤ C γ log γ This completes the proof.

C.2

Proof of Lemma 7

Proof. Let us define γ UW,η (x) :=

1 ⊤ 1 ⊤ 1 x (IN d − W)x + U γ (x) = x (IN d − W)x + F (x) + G(x). 2η 2η 2N γ

(C.10)

Based on Lemma 2.6 in Gurbuzbalaban et al. [2024a], G(x) is ℓ−smooth with ℓ = 4. Therefore, 1−λW

γ UW,η is µ-strongly convex and Lη -smooth with Lη = L + η N + N2γ . Then equation (4.6) can be γ written in terms of UW,η as follows: p γ x(k+1) = x(k) − η∇UW,η x(k) − ηξ (k+1) + 2ηw(k+1) . (C.11)

⊤ Recall that we defined x∗γ = (xγ∗ )⊤ , (xγ∗ )⊤ , · · · , (xγ∗ )⊤ as a vector in the N d dimension, with γ γ d ∗ x∗ ∈ R . Let us also define xη,γ as the minimizer of UW,η . Since the gradients ∇U γ (x) are Lγ -Lipschitz with Lγ = L + N2γ , we have 2 2 = E ∇U γ x(k) − ∇U γ x∗η,γ + ∇U γ x∗η,γ E ∇U γ x(k) 2 2 ≤ 2E ∇U γ x(k) − ∇U γ x∗η,γ + 2 ∇U γ x∗η,γ ≤ 2L2γ E x(k) − x∗η,γ ≤ 2L2γ E x(k) − x∗η,γ ≤ 2L2γ E x(k) − x∗η,γ

2 2 2

+ 2 ∇U γ x∗η,γ

2 2 + 4 ∇U γ x∗γ + 4 ∇U γ x∗η,γ − ∇U γ x∗γ + 4L2γ x∗η,γ − x∗γ

+ 4 ∇U γ x∗γ

(C.12)

To control the term x∗η,γ − x∗γ in (C.12), we state the following result from Yuan et al. [2016]: n o 1+λW 1 N Lemma 14. For any η ≤ min , Lγ Lγ +µ , we have x∗η,γ − x∗γ

√ η N ≤ C1 , 1−ρ

with

ρ := max

W λW 2 , λN

where C1 :=

v u N X 2(Lγ + µ) u t 2Lγ uγi (0) − uγ,∗ , 1+ i µ

with

uγ,∗ := min uγi (x), i

i=1

x∈Rd

(C.13)

γ (x). where uγi (x) := fi (x) + qK

The proof of Lemma 14 can be found in Theorem 1 in Yuan et al. [2016] with some adjustments in notations. 2 Then, we control the term E x(k) − x∗η,γ in (C.12). From (C.11), we have p γ x(k) − ηξ (k+1) + 2ηw(k+1) , x(k+1) − x∗η,γ = x(k) − x∗η,γ − η∇UW,η which implies that E x(k+1) − x∗η,γ

p 2 γ = x(k) − x∗η,γ − η∇UW,η x(k) − ηξ (k+1) + 2ηw(k+1) D E 2 γ = E x(k) − x∗η,γ − 2ηE x(k) − x∗η,γ , ∇UW,η x(k) 2 2 γ + η 2 E ∇UW,η x(k) + η 2 E ξ (k+1) + 2η dN.

(C.14)

γ Since UW,η is µ-strongly convex and Lη -smooth, therefore

D E 2 γ γ γ γ Lη ∇UW,η (z) − ∇UW,η (y), z − y ≥ ∇UW,η (z) − ∇UW,η (y) D E γ γ z − y, ∇UW,η (z) − ∇UW,η (y) ≥ µ ∥z − y∥2 , ∀z, y ∈ Rd . 30

∀z, y ∈ Rd ,

(C.15) (C.16)

γ Also, from the first-optimality condition, we have ∇UW,η x∗η,γ = 0. Thus, from (C.14), we obtain E x(k+1) − x∗η,γ (k)

− x∗η,γ

2 2

(k)

γ − x∗η,γ , ∇UW,η

(k)

γ − ∇UW,η

x∗η,γ

− 2ηE x x 2 2 γ γ + η 2 E ∇UW,η x(k) − ∇UW,η x∗η,γ + η 2 E ξ (k+1) + 2η dN D 2 E γ γ ≤ E x(k) − x∗η,γ − 2ηE x(k) − x∗η,γ , ∇UW,η x(k) − ∇UW,η x∗η,γ D E γ γ + η 2 Lη E x(k) − x∗η,γ , ∇UW,η x(k) − ∇UW,η x∗η,γ + η 2 σ 2 N + 2η dN D 2 E γ γ = E x(k) − x∗η,γ − 2η (1 − ηLη /2) E x(k) − x∗η,γ , ∇UW,η (x(k) ) − ∇UW,η x∗η,γ

=E x

+ η 2 σ 2 N + 2η dN 2 ηLη ≤E x − 2ηµ 1 − E x(k) − x∗η,γ + η 2 σ 2 N + 2η dN 2 2 ηLη E x(k) − x∗η,γ + η 2 σ 2 N + 2η dN = 1 − 2µη 1 − 2 2 = 1 − µη 1 + λW E x(k) − x∗η,γ + η 2 σ 2 N + 2η dN, N − ηLγ (k)

where we used

− x∗η,γ

ηLη 2 < 1 and µη

2 1 + λW − η L + ∈ (0, 1) for any N Nγ 0<γ<

2N η 1 + λW N − ηL

1 − µη

Therefore, from a recursion and simplification, we obtain E x(k) − x∗η,γ

≤ 1 − µη 1 + λW N − ηLγ

x(0) − x∗η,γ

ησ 2 N + 2dN . µ 1 + λW N − ηLγ

Now we can obtain a complete bound for equation (C.12) as follows: 2 2 2 2 ≤ 2L2γ E x(k) − x∗η,γ + 4L2γ x∗η,γ − x∗γ + 4 ∇U γ (x∗γ ) E ∇U γ x(k) ≤ 2L2γ 1 − µη 1 + λW N − ηLγ + 4L2γ x∗η,γ − x∗γ ≤ 4L2γ

x(0) − x∗η,γ

2L2γ (ησ 2 N + 2dN ) µ 1 + λW N − ηLγ

+ 4 ∇U γ (x∗γ ) 2 k (0) 1 − µη 1 + λW x − x∗γ N − ηLγ k ∗ 2 + 4L2γ 1 − µη 1 + λW xη,γ − x∗γ N − ηLγ 2

+ 4L2γ x∗η,γ − x∗γ

+ 4 ∇U γ (x∗γ )

≤ 4L2γ 1 − µη 1 + λW N − ηLγ + 4 ∇U γ (x∗γ )

x(0) − x∗γ

2L2γ (ησ 2 N + 2dN ) µ 1 + λW N − ηLγ

+ 8L2γ x∗η,γ − x∗γ

2L2γ (ησ 2 N + 2dN ) . µ 1 + λW N − ηLγ

Finally, we use Lemma 14 to complete the proof: 2 2 k (0) C 2η2N ≤ 4L2γ 1 − µη 1 + λW E ∇U γ x(k) x − x∗γ + 8L2γ 1 N − ηLγ (1 − ρ)2 2L2γ (ησ 2 N + 2dN ) 2 . + 4 ∇U γ (x∗γ ) + µ 1 + λW N − ηLγ The proof is complete.

C.3

Proof of Lemma 8

Proof. Recall the concatenated version of the algorithm from equation (4.6): p x(k+1) = Wx(k) − η∇U γ x(k) − ηξ (k+1) + 2ηw(k+1) , where W = W ⊗ Id . From a telescoping computation, we obtain that k−1 X x(k) = W k ⊗ Id x(0) − η W k−1−r ⊗ Id ∇U γ x(r) r=0

−η

k−1 X

k−1−r

⊗ Id ξ

r+1

k−1 p X W k−1−r ⊗ Id w(s+1) . + 2η r=0

r=0

Next, define the concatenated average x̄ ∈ RN d as ⊤ ⊤ ⊤ i 1 h (k) (k) (k) x̄ := x̄ , · · · , x̄ = 1N 1⊤ N ⊗ Id x . N

(C.17)

Therefore, i 1 h (k) x(k) − 1N 1⊤ ⊗ I d x N N i 1 h k 1N 1⊤ W ⊗ I = W k ⊗ Id x(0) − x(0) d N N k−1 k−1 i X X 1 h k−1−r γ (r) k−1−r −η W ⊗ Id ∇U x +η 1N 1⊤ ⊗ Id U γ x(r) NW N −η

r=0 k−1 X

r=0

W k−1−r ⊗ Id ∇ξ (r+1) + η

r=0

p 2η

k−1 X r=0

k−1 X

W k−1−r ⊗ Id ∇w(r+1) −

r=0

i 1 h k−1−r (r+1) 1N 1⊤ W ⊗ I d ξ N N

k−1 i p X 1 h k−1−r 2η 1N 1⊤ W ⊗ I w(r+1) . d N N r=0

Since W is doubly stochastic, i i 1 h 1 h k ⊤ = ⊗ I 1N 1⊤ W 1 1 ⊗ I N N d d , N N N

k ≥ 1.

for all

Now applying the Cauchy-Schwarz inequality, we have 2 i 1 h (k) x 1N 1⊤ ⊗ I d N N 2 i 1 h (0) k ⊗ I x 1N 1⊤ W ≤ 4 W k ⊗ Id x(0) − d N N

x(k) −

+ 4 −η

k−1 X

k−1−r

⊗ Id ∇U

k−1−r

(r)

r=0

+ 4 −η

k−1 X

⊗ Id ∇ξ

(r+1)

+η

k−1 X r=0

r=0

p +4 2η

k−1 i X 1 h (r) γ k−1−r +η x U ⊗ I 1N 1⊤ W d N N

k−1 X

k−1−r

⊗ Id ∇w

(r+1)

−

2η

r=0

W k ⊗ Id x(0) −

+ 4 −η

k−1 X

k−1−r

1 h N

1N 1⊤ N ⊗ Id x

⊗ Id ∇U

(r)

r=0

+ 4 −η

k−1 X

k−1 X 1 h

r=0 2 (0)

2 k−1−r 1N 1⊤ NW

⊗ Id w

W k−1−r ⊗ Id ∇ξ (r+1) + η

k−1 X

k−1 X r=0

(r+1)

k−1 i 2 X 1 h ⊤ γ (r) +η 1N 1N ⊗ I d U x N r=0

r=0

p +4 2η

i 1 h (r+1) k−1−r ⊗ I 1N 1⊤ W d ξ N N

W k−1−r ⊗ Id ∇w(r+1) −

2 i 1 h (r+1) 1N 1⊤ N ⊗ Id ξ N

2η

r=0

k−1 X 1 h r=0

(r+1) 1N 1⊤ N ⊗ Id w

Factoring the common terms, 2 i 1 h (k) 1N 1⊤ ⊗ I x d N N " # 2 1 k ⊤ (0) ≤4 W − 1N 1N ⊗ I d x N " # k−1 2 X 1 2 k−1−r ⊤ γ (r) + 4η W − 1N 1N ⊗ Id ∇U x N r=0 " # 2 k−1 X 1 2 k−1−r ⊤ (r+1) + 4η W − 1N 1N ⊗ I d ξ N r=0 " # 2 k−1 X 1 k−1−r ⊤ (r+1) + 8η W − 1N 1N ⊗ I d w . N

x(k) −

(C.18)

r=0

W k Recall that W has the spectral gap 1 − ρ with ρ := max{|λW 2 |, |λN |}; therefore, W has the W k W W W eigenvalues (λi ) with 1 = λ1 > λ2 > · · · > λN > −1 for any k ≥ 1. This results

o n 1 W k−1−r W k−1−r = ρk−1−r 1N 1⊤ = max |λ | , |λ | N 2 N N o n 1 W k W k W k − 1N 1⊤ = ρk . N = max |λ2 | , |λN | N

W k−1−r −

Therefore, " # k−1 2 X 1 (r) k−1−r γ ⊤ 4η W − 1N 1N ⊗ Id ∇U x N r=0 " k−1 # 2 X 1 γ ≤ 4η 2 W k−1−r − 1N 1⊤ x(r) N ⊗ Id · ∇U N r=0 # " k−1 2 X 1 ≤ 4η 2 W k−1−r − 1N 1⊤ · ∇U γ x(r) N N r=0 " k−1 # 2 X 2 k−1−r γ (r) = 4η ρ · ∇U x 2

r=0

= 4η 2

k−1 X r=0

!2 ρk−1−r

Pk−1

k−1−r ∇U γ x(r) r=0 ρ Pk−1 k−1−r r=0 ρ

!2 .

(C.19)

For a real convex function Ψ : R → R, for any points (xi )ni=1 in the domain of Ψ, and for any ai > 0, Jensen’s inequality implies that Pn Pn ai xi a Ψ(xi ) i=1 Pn i Ψ Pn ≤ i=1 . i=1 ai i=1 ai 34

Using Jensen’s inequality in (C.19), we obtain " # k−1 2 X 1 k−1−r ⊤ γ (r) 4η W − 1N 1N ⊗ Id ∇U x N r=0 !2 k−1 k−1 2 X X ρk−1−r γ (r) 2 k−1−r . ≤ 4η ρ Pk−1 k−1−r ∇U x r=0 ρ r=0 r=0 2

Similarly, we compute that " 4

1 W k − 1N 1⊤ N N

⊗ Id x

(0)

2 2 2 1 k ⊤ ≤ 4 W − 1N 1N ⊗ I d x(0) ≤ 4ρ2k x(0) . N

Now we are ready to compute a bound between the individual iterates and their mean N X

(k)

E xi

− x̄(k)

i=1

=E x

(k)

2 i 1 h ⊤ (k) − 1N 1N ⊗ I d x N

≤ 4ρ2k E x(0)

k−1 X

+ 4η 2 E

ρk−1−r

!2 k−1 X

r=0

ρk−1−r Pk−1 r=0

ρk−1−r

∇U γ x(r)

" # 2 k−1 X 1 2 ⊤ k−1−r (r+1) + 4η E W − 1N 1N ⊗ I d ξ N r=0 # " 2 k−1 X 1 ⊤ (r+1) k−1−r . + 8ηE W − 1N 1N ⊗ I d w N

(C.20)

r=0

Recall from Lemma 7 that for every k = 1, 2, 3 · · · , 2 E ∇U γ x(k) ≤ Dγ2 , where Dγ2 is defined in equation (A.3): Dγ2 := 4L2γ 1 − µη 1 + λW N − ηLγ + 4 ∇U γ (x∗γ )

x(0) − x∗γ

+ 8L2γ

2L2γ (ησ 2 N + 2dN ) . µ 1 + λW N − ηLγ

C12 η 2 N (1 − ρ)2

Therefore, the second term on the right hand side of (C.20), k−1 X

4η 2 E

ρk−1−r

!2 k−1 X

ρk−1−r

∇U γ x(r) k−1−r ρ r=0 r=0 r=0 !2 k−1 k−1 X X ρk−1−r ≤ 4η 2 Dγ2 ρk−1−r Pk−1 k−1−r r=0 ρ r=0 r=0 2 k−1 k−1 k−1 X 1 − ρk−1 X ρk−1−r 2 2 1−ρ ρk−1−r = 4η 2 Dγ2 = 4η D γ 1−ρk−1 1−ρ 1−ρ r=0 r=0 1−ρ 2 k−1 1−ρ 1 = 4η 2 Dγ2 ≤ 4η 2 Dγ2 . 2 (1 − ρ) (1 − ρ)2 Pk−1

Using the result from (C.21) and Assumption 4 in (C.20), we obtain: N X

(k)

E xi

− x̄(k)

i=1

≤ 4ρ2k E x(0)

k−1

+ 4η 2 Dγ2

2 X 2 1 1 2 k−1−r ⊤ + 4η W − 1 1 E ξ (r+1) N N 2 (1 − ρ) N r=0

+ 8η

k−1 X

W k−1−r −

r=0

≤ 4ρ2k E x(0) ≤ 4ρ2k E

1 1N 1⊤ N N

E w(r+1)

k−1

+ 4η 2 Dγ2

k−1

X X 1 + 4η 2 N σ 2 ρ2(k−1−r) + 8ηdN ρ2(k−1−r) 2 (1 − ρ)

r=0 2 2 2 2 2 4η Dγ 4η N σ 8ηdN x(0) + + + . 2 2 (1 − ρ) (1 − ρ ) (1 − ρ2 )

r=0

Therefore, N 2 2 4η 2 Dγ2 1 X 4ρ2k 4η 2 σ 2 8ηd (k) E xi − x̄(k) ≤ E x(0) + + + . 2 2 N N N (1 − ρ) (1 − ρ ) (1 − ρ2 ) i=1

This completes the proof.

(C.21)

C.4

Proof of Lemma 9

Proof. By taking the L2 norm of the error term, we get E ∥Ek+1 ∥2 " # " # N N 2 X 1 X 1 (k) (k) =E ∇fi x̄(k) − ∇fi xi + 2 PK x̄(k) − PK xi N N γ i=1 i=1 # " " # 2 N N 2 X 1 X 1 (k) (k) ≤ 2E + 2E ∇fi x̄(k) − ∇fi xi PK x̄(k) − PK xi N N 2γ i=1

≤

2 N2

N X

2 (k) + N E ∇fi x̄(k) − ∇fi xi

i=1 N 2 2L X

N i=1 = 2 L2 + ≤2 L+ = 2L2γ

(k)

E xi 4

− x̄(k) N 1 X

N 2γ2

2 Nγ

N 1 X

i=1

N 1 N (k)

E xi

N 3γ2 (k)

E xi

i=1 N X

N X

(k)

E xi

2 4 N γ2

(k)

E xi

N X

2 (k) N E PK x̄(k) − PK xi

i=1

− x̄(k)

i=1

− x̄(k)

i=1

− x̄(k)

(C.22)

i=1

Using the results from Lemma 8, we obtain 2 8L2γ ρ2k 8L2γ η 2 Dγ2 8L2γ η 2 σ 2 16L2γ ηd (0) E ∥Ek+1 ∥ ≤ E x + + + . N N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) 2

This completes the proof.

C.5

Proof of Lemma 10

Proof. From (A.6) and (A.7) we obtain: h i x̄(k+1) − x̃(k+1) = x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) + ηEk+1 − η ξ¯(k+1) .

(C.23)

Therefore, x̄(k+1) − x̃(k+1)

h i 2 = x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) + ηEk+1 − η ξ¯(k+1) h i 2 2 + η 2 Ek+1 − ξ¯(k+1) = x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) D h i E + 2 x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) , ηEk+1 − η ξ¯(k+1) 2 2 2 + η 2 Ek+1 − ξ¯(k+1) = x̄(k) − x̃(k) + η 2 ∇ûγ x̄(k) − ∇ûγ x̃(k) D h iE − 2 x̄(k) − x̃(k) , η ∇ûγ x̄(k) − ∇ûγ x̃(k) D h i E + 2 x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) , ηEk+1 − η ξ¯(k+1) E 2 Lγ D (k) x̄ − x̃(k) , ∇ûγ x̄(k) − ∇ûγ x̃(k) ≤ x̄(k) − x̃(k) + η 2 N iE D h 2 (k) (k) + η 2 Ek+1 − ξ¯(k+1) − 2 x̄ − x̃ , η ∇ûγ x̄(k) − ∇ûγ x̃(k) E i D h + 2 x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) , ηEk+1 − η ξ¯(k+1) 2 2 ηLγ 2µη ≤ 1− 1− x̄(k) − x̃(k) + η 2 Ek+1 − ξ¯(k+1) N 2N E i D h + 2 x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) , ηEk+1 − η ξ¯(k+1) ,

(C.24)

where we used Lγ /N -smoothness and µ/N -strongly convexity of ûγ (x), and the assumption that 2N ¯(k+1) has mean zero η < 2N Lγ = L+ 2 in obtaining (C.24). Also, by Assumption 4, the term ξ Nγ

conditional on the natural filtration of the iterates till time k. Thus, by taking expectations, it yields that 2

E x̄(k+1) − x̃(k+1) 2 2 ηLγ 2µη ≤ 1− 1− E x̄(k) − x̃(k) + η 2 E Ek+1 − ξ¯(k+1) N 2N D h i E + 2E x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) , ηEk+1 − η ξ¯(k+1) 2 2 ηLγ 2µη = 1− 1− E x̄(k) − x̃(k) + η 2 E ∥Ek+1 ∥2 + η 2 E ξ¯(k+1) N 2N D h i E + 2E x̄(k) − x̃(k) − η ∇ûγ x̄(k) − ∇ûγ x̃(k) , ηEk+1 − η ξ¯(k+1) 2 2 ηLγ 2µη ≤ 1− 1− E x̄(k) − x̃(k) + η 2 E ∥Ek+1 ∥2 + η 2 E ξ¯(k+1) N 2N # " ηLγ +2 1+ ηE x̄(k) − x̃(k) · ∥Ek+1 ∥ . N

(C.25)

For any a, b ≥ 0 and ε > 0, we have the Young’s inequality 2ab ≤ εa2 +

a = x̄

(k)

− x̃

(k)

b = ∥Ek+1 ∥ ,

ε=

µ N

b2 . Choose ε

ηLγ 2 ηLγ 1+ N

1−

Then, from (C.25), we obtain 2

E x̄(k+1) − x̃(k+1) 2 ηLγ 2µη σ2 ≤ 1− 1− E x̄(k) − x̃(k) + η 2 E ∥Ek+1 ∥2 + η 2 N 2N N   ηL γ µ 1 − 2 2N N + ηLγ ηLγ E ∥Ek+1 ∥2  η + 1+ E x̄(k) − x̃(k) + ηLγ N N + ηLγ µ 1 − 2N   ηLγ 2 2 ) (1 + ηLγ σ2 µη 1− E x̄(k) − x̃(k) + η η + N  E ∥Ek+1 ∥2 + η 2 . = 1− ηL µ N 2N N 1− γ N

(C.26)

ηLγ Here, we assume that the leading term 1 − µη ∈ [0, 1). Now using the results from N 1 − 2N Lemma 9, we obtain 2

E x̄(k+1) − x̃(k+1) 2 ηLγ µη 1− E x̄(k) − x̃(k) ≤ 1− N 2N   ! ηLγ 2 2 8L2γ η 2 σ 2 16L2γ ηd 8L2γ ρ2k 8L2γ η 2 Dγ2 (1 + N ) (0)  + η η + E x + + + ηL µ N N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) 1− γ N

+ η2

σ2 N

Since E x̄(0) − x̃(0) = 0, by iterations, 2

E x̄(k) − x̃(k) r k−1 X ηLγ µη ≤ 1− 1− N 2N r=0  "  ηLγ 2 (1 + ) · η η + N  ηLγ µ N 1 − 2N

! # 2 8L2γ η 2 σ 2 8L2γ η 2 Dγ2 16L2γ ηd σ + + η2 + N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) N   r k−1 X 2 ηLγ (1 + ηLγ )2  8L2γ ρ2(k−r) µη 1− + 1− η η + E x(0) ηL N 2N N µ 1− γ r=0

ηL

γ 1 − 1 − µη N 1 − 2N = ηLγ 1 − 1 − µη 1 − N 2N  "  ! # ηL 2 8L2γ η 2 σ 2 16L2γ ηd 8L2γ η 2 Dγ2 (1 + Nγ )2 σ 2  · η η + + + +η ηLγ µ N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) N 1 − N 2N k ηLγ ρ2k − 1 − µη 1 − 2 8L2γ N 2N E x(0) . + ηLγ 1 N 1 − 1 − µη N 1 − 2N ρ2

Finally, for every k ∈ N, E x̄(k) − x̃(k)

ηL

η ≤

(1+ γ )2 η + µ NηLγ 1− 2N N

8L2γ η 2 Dγ2 8L2γ η 2 σ 2 16L2γ ηd + (1−ρ 2 ) + (1−ρ2 ) N (1−ρ)2 µη N

ηL

1 − 2Nγ k

+ η 2 σN

ηLγ ρ2k − 1 − µη 2 8L2γ ρ2 N 1 − 2N + · E x(0) ηLγ N ρ2 − 1 − µη N 1 − 2N ! 2 2 ηL 2 8Lγ ηDγ 8L2γ ησ 2 16L2γ d (1+ Nγ )2 η η + µ ηLγ + + + ησ 2 2 2 N N (1−ρ) (1−ρ ) (1−ρ ) N

1− 2N

µ N

ηL

1 − 2Nγ k

ηLγ ρ2k − 1 − µη 2 8L2γ ρ2 N 1 − 2N + E x(0) · ηLγ N ρ2 − 1 − µη N 1 − 2N  

ηL 8L2γ ησ 2 16L2γ d (1 + Nγ )2  8L2γ ηDγ2  η = η + + +  2 ηLγ µ N (1 − ρ)2 (1 − ρ2 ) (1 − ρ2 ) ηL µ2 1 − 2Nγ N 1 − 2N N2 k 2k − 1 − µη 1 − ηLγ ρ 2 2 8L2γ ρ2 N 2N ησ + · + E x(0) . ηL ηLγ N µ 1 − 2Nγ ρ2 − 1 − µη N 1 − 2N

This completes the proof.

C.6

Proof of Lemma 12

Proof. Define the process Yt := Xt − x̂∗γ , where Xt follows the overdamped Langevin diffusion: dXt := −∇ûγ (Xt )dt +

√

2N −1 dWt ,

(C.27)

where Wt is a standard d-dimensional Brownian motion, such that X∞ follows the unique stationary distribution π γ of (C.27). Then, we have √ dYt = −∇ûγ (Xt )dt + 2N −1 dWt . Now by Itô’s formula, √ d∥Yt ∥2 = −2 Yt , ∇ûγ (Xt ) − ∇ûγ (x̂∗γ ) dt + 2 2N −1 ⟨Yt , dWt ⟩ + (2N −1 d) dt.

Multiplying both sides by the integrating factor e N t and using the product rule yields Z t √ µ µ t 2 2 −1 e N ∥Yt ∥ = ∥Y0 ∥ + 2 2N e N s ⟨Ys , dWs ⟩ 0 Z t µ e N s Ys , ∇ûγ (Xs ) − ∇ûγ x̂∗γ ds −2 0 Z t Z µ µ t µs −1 s e N ds + + 2N d e N ∥Ys ∥2 ds. N 0 0 Since ûγ is µ/N -strongly convex, we have ⟨x − y, ∇ûγ (x) − ∇ûγ (y)⟩ ≥

µ ∥x − y∥2 . N

Applying this inequality with x = Xs and y = x̂∗γ implies Ys , ∇ûγ (Xs ) − ∇ûγ (x̂∗γ ) ≥

µ ∥Ys ∥2 . N

Therefore, Z t −2

µ s N

Ys , ∇û (Xs ) − ∇û

x̂∗γ

µ ds + N

Z t

e N s ∥Ys ∥2 ds ≤ 0.

Dropping this non-positive term yields e

µ t N

∥Yt ∥ ≤ ∥Y0 ∥ + 2

√

2N −1

Z t e

µ s N

⟨Ys , dWs ⟩ + 2N

−1

Z t d

e N s ds.

Taking expectations and using the fact that the stochastic integral has mean zero gives Z t i h µ µ 2N −1 d µ t 2 −1 t 2 N e N s ds = E∥Y0 ∥2 + eN − 1 . E e ∥Yt ∥ ≤ E∥Y0 ∥ + 2N d µ 0

Dividing by e

µ t N

yields µ

E∥Yt ∥2 ≤ e− N t E∥Y0 ∥2 +

2N −1 d µ N

µ µ 2N −1 d 1 − e− N t ≤ e− N t E∥Y0 ∥2 + . µ N

Letting t → ∞ and using convergence to the invariant distribution π γ gives E∥X∞ − x̂∗γ ∥2 ≤ This completes the proof.

2N −1 d µ N

2d . µ

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