Year: 2018

Abstract:
Many deep-learning based methods such as Bayesian deep learning (DL) and deep reinforcement learning (RL) have heavily relied on the ability of a model being able to efficiently explore via Bayesian sampling. Particle-optimization sampling (POS) is a recently developed technique to generate high-quality samples from a target distribution by iteratively updating a set of interactive particles, with a representative algorithm the Stein variational gradient descent (SVGD). Though obtaining significant empirical success, the non-asymptotic convergence behavior of SVGD remains unknown. In this paper, we generalize POS to a stochasticity setting by injecting random noise in particle updates, called stochastic particle-optimization sampling (SPOS). Notably, for the first time, we develop non-asymptotic convergence theory for the SPOS framework, characterizing convergence of a sample approximation w.r.t. the number of particles and iterations under both convexand noncovexenergy-function settings. Interestingly, we provide theoretical understanding of a pitfall of SVGD that can be avoided in the proposed SPOS framework, i.e., particles tend to collapse to a local mode in SVGD under some particular conditions. Our theory is based on the analysis of nonlinear stochastic differential equations, which serves as an extension and a complementary development to the asymptotic convergence theory for SVGD such as (Liu, 2017). With such theoretical guarantees, SPOS can be safely and effectively applied on both Bayesian DL and deep RL tasks. Extensive results demonstrate the effectiveness of our proposed framework.

text: Many deep-learning based methods such as Bayesian deep learning (DL) and deep reinforcement learning (RL) have heavily relied on the ability of a model being able to efficiently explore via Bayesian sampling. Particle-optimization sampling (POS) is a recently developed technique to generate high-quality samples from a target distribution by iteratively updating a set of interactive particles, with a representative algorithm the Stein variational gradient descent (SVGD). Though obtaining significant empirical success, the non-asymptotic convergence behavior of SVGD remains unknown. In this paper, we generalize POS to a stochasticity setting by injecting random noise in particle updates, called stochastic particle-optimization sampling (SPOS). Notably, for the first time, we develop non-asymptotic convergence theory for the SPOS framework, characterizing convergence of a sample approximation w.r.t. the number of particles and iterations under both convex- and noncovexenergy-function settings. Interestingly, we provide theoretical understanding of a pitfall of SVGD that can be avoided in the proposed SPOS framework, i.e., particles tend to collapse to a local mode in SVGD under some particular conditions. Our theory is based on the analysis of nonlinear stochastic differential equations, which serves as an extension and a complementary development to the asymptotic convergence theory for SVGD such as (Liu, 2017). With such theoretical guarantees, SPOS can be safely and effectively applied on both Bayesian DL and deep RL tasks. Extensive results demonstrate the effectiveness of our proposed framework.

heading: 1 INTRODUCTION
text: Recent development of deep-learning techniques have required the ability of an algorithm to efficiently explore some particular space (e.g., parameter space) via Bayesian sampling, due to high model complexity of modern deep models. Meanwhile, recent years have seen the development of scalable sampling methods such as stochastic gradient MCMC (SG-MCMC) (Welling & Teh, 2011; Chen et al., 2014; Ding et al., 2014; Chen et al., 2015) and Stein variational gradient descent (SVGD) (Liu & Wang, 2016) to facilitate big-data analysis. SG-MCMC is a family of scalable Bayesian sampling algorithms relying on Itó diffusions, linear stochastic differential equations (SDE) with appropriately designed coefficients such that the corresponding stationary distributions match a target distribution. One potential issue is that samples might be highly correlated partially due to the nature of Markov chains, leading to low sample efficiency, an undesirable property of SG-MCMC.
SVGD, on the other hand, belongs to the family of particle-based sampling methods, which optimize a set of interacting particles to minimize some distance metric between the target distribution and the distribution induced by the particles. By optimizing the distance measure, one maintains an optimal set of particles at each time. Recent development of SVGD has shown that the underlying mathematical principle is based on a family of nonlinear SDEs, in the sense that coefficients of the SDE depend on the current density of the particles. Though achieving numerous practical successes (Liu & Wang, 2016; Feng et al., 2017; Liu et al., 2017; Haarnoja et al., 2017; Liu & Zhu, 2018), little theory has been developed to understand the convergence property of the algorithm. A recent theoretical development has interpreted SVGD as a special type of gradient flow in the space of probability measures, and developed theory to disclose its asymptotic convergence behavior (Liu, 2017).
Recently, Chen et al. (2018) unified SG-MCMC and SVGD by proposing a particle-optimization sampling framework that interprets both as Wasserstein gradient flows (WGFs). Generally speaking, a WGF is a partial differential equation (PDE) defined on the space of probability measures, describing the evolution of a density function over time. In (Chen et al., 2018), the authors define a WGF by combining the corresponding Fokker-Planck equations for both SG-MCMC and SVGD, and solve it with deterministic particle approximations. However, due to the diffusion nature, deterministicparticle approximation leads to a hard-to-control error, challenging for theoretical analysis.
Based on the unified framework in (Chen et al., 2018), we propose to solve WGFs with stochastic particle approximation, leading to stochastic particle-optimization sampling (SPOS). The idea is instead of solving the WGF with an uncontrollable deterministic approximation for a diffusion term, we solve it stochastically by injecting random noise to the particle-update equations. Remarkably, for the first time, we develop nonasymptotic convergence theory for the family of SPOS algorithms, considering both convex- and nonconvex-energy functions. Different from existing theory for SGMCMC based algorithms (Teh et al., 2016; Vollmer et al., 2016; Chen et al., 2015; Raginsky et al., 2017; Zhang et al., 2017; Xu et al., 2018), our development relies on the theory of nonlinear SDEs, which is more involved and less explored in literature. Particularly, we adopt tools from granular media equations (Malrieu, 2003; Cattiaux et al., 2008) to develop nonasymptotic error bounds in terms of the 1-Wasserstein distance. Please refer to Section M in the Supplementary Material (SM) for detailed distinctions of our work to existing work. Within our theoretical framework, we provide a formal theoretical understanding of a pitfall of SVGD, i.e., particles tend to collapse to one point under some particular conditions; whereas this can be avoided in the proposed SPOS framework due to the injected random noise. Our theory is general enough to be applied for various deep-learning tasks, including Bayesian deep learning and Bayesian exploration in reinforcement learning. Our extensive experimental results well suggest advantages of our framework compared to existing methods.

heading: 2 PRELIMINARIES
text: This section introduces necessary preliminaries, along with notations used in this paper. For the sake of clarity, through out the paper, we use bold letters to denote variables in continuous-time diffusions and model definitions, e.g., θt in (1) defined below (indexed by “time” t). By contrast, normal unbold letters are used to denote parameters in algorithms (discrete solutions of continuous-time diffusions), e.g., θ(i)k in (3) below (indexed by “iteration” k). For conciseness, all the proofs as well as some extra experimental results are presented in the SM. Discussion on the complexity of our algorithm is also included in Section L of the SM.

heading: 2.1 STOCHASTIC GRADIENT MCMC
text: In Bayesian sampling, we aim to generate random samples from a posterior distribution p(θ|X ) ∝ p(X|θ)p(θ), where θ ∈ Rd represents the model parameter with a prior distribution p(θ), and X , {xq}Nq=1 represents the observed data with likelihood p(X|θ) = ∏ q p(xq |θ). Define the potential en-
ergy as: U(θ) , − log p(X|θ)− log p(θ) = − ∑N q=1 ( log p(xq |θ) + 1N log p(θ) ) , ∑N q=1 Uq(θ). SG-MCMC algorithms belong to diffusion-based sampling methods, where a continuous-time diffusion process is designed such that its stationary distribution matches the target posterior distribution. The diffusion process is driven by a specific stochastic differential equation (SDE). For example, in stochastic gradient Langevin dynamic (SGLD), the SDE endows the following form:
dθt = −β−1F (θt)dt+ √ 2/βdWt , (1)
where F (θ) , ∇θU(θ) = ∑N q=1∇θUq(θ) , ∑N q=1 Fq(θ); t is the time index; β > 0 is the temperature parameter; andWt ∈ Rd is a d-dimensional Brownian motion. More instances of SDEs corresponding to other SG-MCMC algorithms can be defined by specifying different forms of F and potentially other diffusion coefficients. We focus on SGLD and (1) in this paper, and refer interested readers to (Ma et al., 2015) for more detailed description of general SG-MCMC algorithms. Denote the probability density function of θt in (1) as νt, and a ·b , a> b for two vectors a and b. It is known that νt is characterized by the following Fokker-Planck (FP) equation (Risken, 1989): ∂tνt = ∇θ · (β−1νtF (θ) + β−1∇θνt) . (2) where the stationary distribution ν∞ equals to our target distribution p(θ|X ) according to Chiang & Hwang (1987). SGLD generates samples from p(θ|X ) by numerically solving the SDE (1). For
scalability, it replaces F (θk) in each iteration with an unbiased evaluation by randomly sampling a subset of X , i.e. F (θk) is approximated by: Gk , NBk ∑ q∈Ik Fq(θk), where Ik is a random subset of [1, 2, · · · , N ] with size Bk in each iteration. Based on the above settings, SGLD uses the Euler method with stepsize hk to numerically solve (1) and obtains the update equation: θk+1 = θk − β−1Gkhk + √ 2β−1hkξk, ξk ∼ N (0, I).

heading: 2.2 STEIN VARIATIONAL GRADIENT DESCENT
text: Different from SG-MCMC, SVGD is a deterministic particle-optimization algorithm that is able to generate samples from a target distribution. In the algorithm, a set of particles interact with each other, driving them to high density locations in the parameter space while keeping them far away from each other with repulsive force. The update equations of the particles follow the fastest descent direction of the KL-divergence between current empirical distribution of the particles and the target distribution, on an RKHS induced by a kernel function κ(·, ·) (Liu & Wang, 2016). Formally, Liu & Wang (2016) derived the following updating rules for the particles {θ(i)k }Mi=1 at the k-th iteration with stepsize hk and G (i) k , N Bk ∑ q∈Ik Fq(θ (i) k ):
θ (i) k+1 = θ (i) k + hk M M∑ j=1 [ κ(θ (j) k , θ (i) k )G (i) k +∇θ(j)k κ(θ (j) k , θ (i) k ) ] , ∀i (3)
where the first term in the bracket encourages particles to locate on high density modes, and the second term serves as repulsive force that pushes away different particles. Different from SG-MCMC, only particles at the current iteration, {θ(i)k }Mi=1, are used to approximate the target distribution.

heading: 2.3 PARTICLE-OPTIMIZATION BASED SAMPLING METHODS
text: SG-MCMC and SVGD, though look closely related, behave very differently in terms of algorithms, e.g., stochastic and noninteractive versus deterministic and interactive particle updates. Recently, Chen et al. (2018) proposed a deterministic particle-optimization framework that unifies both SGMCMC and SVGD. Specifically, the authors viewed both SG-MCMC and SVGD as Wasserstein gradient flows (WGFs) on the space of probabilistic measures, and derived several deterministic particle-optimization techniques for particle evolutions, like what SVGD does. For SG-MCMC, the FP equation (2) for SGLD is a special type of WGFs. Together with an interpretation of SVGD as a special case of the Vlasov equation in nonlinear PDE literature, Chen et al. (2018) proposed a general form of PDE to characterize the evolution of the density for the model parameter θ, denoted as νt at time t with ν∞ matching our target (posterior) distribution, i.e.,
∂tνt = ∇θ · ( νtβ −1F (θ) + νt (K ∗ νt(θ)) + β−1∇θνt ) , (4)
where K is a function controlling the interaction of particles in the PDE system. For example, in SVGD, Chen et al. (2018) showed that K and K ∗ νt(θ) endow the following forms:
K(θ,θ′) , F (θ′)κ(θ′,θ)−∇θ′κ(θ′,θ) and K ∗ νt(θ) , ∫ K(θ,θ′)νt(θ′)dθ′ (5)
where κ(·, ·) is a kernel function such as the RBF kernel. In the following, we introduce a new unary function K(θ) = exp(−‖θ‖ 2
η2 ), thus κ(θ,θ ′) can be rewritten as κ(θ,θ′) = K(θ − θ′). Hence, (4)
with K defined in (5) can be equivalently rewritten as: ∂tνt = ∇θ · ( νtβ −1F (θ) + νt (EY∼νtK(θ − Y )F (Y )−∇K ∗ νt(θ)) + β−1∇θνt ) , (6) where Y is a random sample from νt but independent of θ. Importantly,
Proposition 1 (Chen et al. (2018)) The stationary distribution of (6) equals to our target distribution, which means ν∞(θ) = p(θ|X ).
Chen et al. (2018) proposed to solve (4) numerically with deterministic particle-optimization algorithms such as the blob method. Specifically, the continuous density νt is approximated by a set of M particles {θ(i)t }Mi=1 that evolve over time t, i.e. νt ≈ 1M ∑M i=1 δθ(i)t (θ), where δ θ (i) t (θ) = 1 if θ = θ (i) t and 0 otherwise. Note ∇θνt in (4) is no longer a valid definition when adopting particle approximation for νt. Consequently,∇θνt needs nontrivial approximations, e.g., by discrete gradient flows or blob methods proposed in (Chen et al., 2018). We omit the details here for simplicity.

heading: 3 STOCHASTIC PARTICLE-OPTIMIZATION SAMPLING ALGORITHMS
text: The deterministic particle-approximation methods proposed by Chen et al. (2018) to approximately solve the WGF problem (4) introduce approximation errors for ∇θνt that are hard to control analytically. To overcome this problem, we propose to solve (4) stochastically to replace the ∇θνt term with a Brownian motion. Specifically, first note that the term β−1∇θ · ∇θνt is contributed from Brownian motion, i.e., solving the SDE, dθt = √ 2β−1dWt, is equivalent to solving the corresponding FP equation: ∂νt = β−1∇θ · ∇θνt. Consequently, we decompose RHS of (4) into two parts: F1 , ∇θ · ( νtβ −1F (θt) + (K ∗ νt)νt ) and F2 , β−1∇θ · ∇θνt. Our idea is to solve F1 deterministically under a PDE setting, and solve F2 stochastically based on its corresponding SDE. When adopting particle approximation for the density νt, both solutions of F1 and F2 are represented in terms of particles {θ(i)t }. Thus we can combine the solutions from the two parts directly to approximate the original exact solution of (4). Similar to the results of SVGD in Section 3.3 in (Liu, 2017), we first formally show in Theorem 2 that when approximating νt with particles, i.e., νt ≈ 1M ∑M i=1 δθ(i)t
(θ), the PDE can be transformed into a system of deterministic differential equations with interacting particles.
Theorem 2 When approximating νt in (4) with particles {θ(i)t }i, the PDE ∂tνt = F1 reduces to the following system of differential equations describing evolutions of the particles over time: ∀i
dθ (i) t = −β −1F (θ (i) t )dt−
1
M M∑ j=1 K(θ (i) t − θ (j) t )F (θ (j) t )dt+ 1 M M∑ j=1 ∇K(θ(i)t − θ (j) t )dt (7)
On the other hand, by solving ∂tνt = F2 stochastically with its equivalent SDE counterpart, we arrive at the following differential equation system, describing evolutions of the particles {θ(i)t } over time t: ∀i dθ
(i) t = ( −β−1F (θ(i)t )− 1
M M∑ j=1 K(θ (i) t − θ (j) t )F (θ (j) t ) + 1 M M∑ j=1 ∇K(θ(i)t − θ (j) t )
) dt+ √ 2β−1dW(i)t
Our intuition is that if the particle evolution (8) can be solved exactly, the solution of (6) νt will be wellapproximated by the particles {θ(i)t }Mi=1. In our theory, we show this intuition is actu-
ally true. In practice, however, solving (8) is typically infeasible, and thus numerical methods are adopted. Furthermore, in the case of big data, following SG-MCMC, F (θ(i)k ) is typically replaced by a stochastic version G(i)k , N Bk ∑ q∈Ik Fq(θ (i) k ) evaluated with a minibatch of data of size Bk for computational feasibility. Based on the Euler method (Chen et al., 2015) with a stepsize hk, (8) leads to the following updates for the particles at the k-th iteration
θ (i) k+1 = θ (i) k −hkβ −1G (i) k − hk M M∑ j=1 K(θ (i) k − θ (j) k )G (j) k
+ hk M M∑ j=1 ∇K(θ(i)k − θ (j) k ) + √ 2β−1hkξ (i) k , ξ (i) k ∼ N (0, I) ∀i . (9)
We called the algorithm with particle update equations (9) stochastic particle-optimization sampling (described in Algorithm 3), in the sense that particles are optimized stochastically with extra random Gaussian noise. Intuitively, the added Gaussian noise would enhance the ability of the algorithm to jump out of local modes, leading to better ergodic properties compared to standard SVGD. This serves as one of our motivations to generalize SVGD to SPOS. To illustrate the advantage of introducing the noise term, we compare SPOS and SVGD on sampling a difficult multi-mode distribution, with the density function given in Section A of the SM. The particles are initialized on a local mode close to zero. Figure 1 plots the final locations of the particles along with the true density, which shows that particles in SPOS are able to reach different modes, while they are all trapped at one mode in SVGD. This pitfall of SVGD will be studied formally in Section 4.4.

heading: 4 NON-ASYMPTOTIC CONVERGENCE ANALYSIS: THE CONVEX CASE
text: Algorithm 1 Stochastic Particle-Optimization Sampling
Input: Initial particles {θ(i)0 }Mi=1 with θ (i) 0 ∈ Rd, step size hk, batch size Bk 1: for iteration k= 0,1,...,T do 2: Update θ(i)k+1 with (9) for ∀i. 3: end for Output:{θ(i)T }Mi=1 In this section, we develop nonasymptotic convergence theory for the proposed SPOS when the energy function U(θ) is convex. The nonconvex case is discussed in Section 5. We prove nonasymptotic convergence rates for SPOS algorithm under the 1- Wasserstein metric W1, a special case of p-Wasserstein metric de-
fined as Wp(µ, ν) = ( infζ∈Γ(µ,ν) ∫ Rd×Rd ‖Xµ −Xν‖ pdζ(Xµ, Xν) )1/p , where Γ(µ, ν) is the set of joint distributions on Rd × Rd with marginal distribution µ and ν. Note that SPOS reduces to SVGD when β → ∞, thus our theory also sheds light on the convergence behavior of SVGD, where non-asymptotic theory is currently missing, despite the asymptotic theory developed recently (Liu, 2017; Lu et al., 2018). It is worth noting that part of our proofs are generalization of techniques for analyzing granular media equations in (Malrieu, 2003; Cattiaux et al., 2008).

heading: 4.1 BASIC SETUP AND ASSUMPTIONS
text: Due to the exchangeability of the particle system {θ(i)t }Mi=1 in (8), if we initialize all the particles θ (i) t with the same distribution ρ0, they would endow the same distribution for each time t. We denote the distribution of each θ(i)t as ρt. Similar arguments hold for the particle system {θ (i) k }Mi=1 in (9), and thus we denote the distribution of each θ(i)k as µk. To this end, our analysis aims at bounding W1(µT , ν∞) since ν∞ equals to our target distribution p(θ|X ) according to Proposition 1. Before proceeding to our theoretical results, we first present the following basic assumptions.
Assumption 1 Assume F and K satisfy the following conditions:
• There exist positive mF and mK , such that 〈F (θ)− F (θ′),θ − θ′〉 ≥ mF ‖θ − θ′‖2 and 〈∇K(θ)−∇K(θ′),θ − θ′〉 ≤ −mK‖θ − θ′‖2.
• F is bounded by HF and LF -Lipschitz continuous i.e., ‖F (θ)‖ ≤ HF and ‖F (θ1) − F (θ2)‖ ≤ LF ‖θ1 − θ2‖; .
• K is LK-Lipschitz continuous;∇K is bounded by H∇K and L∇K-Lipschitz continuous
• F (0) = 0 and K is an even function, i.e., K(−θ) = K(θ).
Note the first bullet indicates U to be a convex function and K to be a concave function. For an RBF kernel, the later could be achieved by setting the bandwidth large enough and only considering the concave region for simplicity. This assumption is used for revealing some undesired property of SVGD developed below. We do not need such an assumption when analyzing under a nonconvex energy function U in Section 5. Then "F (0) = 0" in the second bullet is reasonable, as F in our setting corresponds to an unnormalized log-posterior, which can be shifted such that F (0) = 0 for a specific problem. Since we often care about bounded space in practice, we can realize the third bullet due to the continuity of K and∇K. The high-level idea of bounding W1(µT , ν∞) in this section is to decompose it as follows:
W1(µT , ν∞) ≤W1(µT , ρ∑T−1 k=0 hk ) +W1(ρ∑T−1 k=0 hk , ν∑T−1 k=0 hk ) +W1(ν∑T−1 k=0 hk , ν∞) . (10)

heading: 4.2 BOUNDS WITH STOCHASTIC PARTICLE APPROXIMATION
text: We firstly bound W1(ρ∑T−1 k=0 hk , ν∑T−1 k=0 hk ) and W1(ν∑T−1 k=0 hk , ν∞) with the following theorems.
Theorem 3 Under Assumption 1 and let ρ0 = ν0, there exist some positive constants c1 and c2 independent of (M, t) and satisfying c2 < β−1 such that
W1(ρt, νt) ≤ c1(β−1 − c2)−1M−1/2, ∀t. (11)
Remark 1 According to Theorem 3, we can bound the W1(ρ∑T−1 k=0 hk , ν∑T−1 k=0 hk
) term as W1(ρ∑T−1 k=0 hk , ν∑T−1 k=0 hk ) ≤ c1√ M(β−1−c2) . Furthermore, by letting t→∞, we haveW1(ρ∞, ν∞) ≤ c1√ M(β−1−c2) , which is an important intermediate result to prove the following theorem.
Theorem 4 Under Assumption 1, the following holds: W1(νt, ν∞) ≤ c3e−2λ1t, where λ1 = β−1mF − 3HFLK − 2LF and c3 is some positive constant independent of (M, t). Hence, the W1(ν∑T
k=0 hk , ν∞) term in (10) can be bounded as:
W1(ν∑T−1 k=0 hk , ν∞) ≤ c3 exp ( −2λ1( T−1∑ k=0 hk) ) . (12)
To ensure W1(ν∑T−1 k=0 hk
, ν∞) to decrease over time, one needs to choose β small enough such that λ1 > 0. This also sheds light on a failure case of SVGD (where β →∞) discussed in Section 4.4.

heading: 4.3 BOUNDS WITH A NUMERICAL SOLUTION
text: To bound the W1(µT , ρ∑T−1 k=0 hk
) term, we adopt techniques from (Raginsky et al., 2017; Xu et al., 2018) on analyzing the behaviors of SGLD, and derive the following results for our SPOS algorithm:
Theorem 5 Under Assumptions 1, for a fixed step size hk = h that is small enough, the corresponding W1(µT , ρTh) is bounded as
W1(µT , ρTh) ≤ c4Md 3 2 β−3(c5β 2B−1 + c6h) 1 2T 1 2h 1 2 (13)
whereB is the fixed size of the minibatch in each iteration and (c4, c5, c6) are some positive constants independent of (M,T, h).
The dependence of T in the bound above makes the bound relatively loose. Fortunately, we can improve the bound to make it independent of T by considering a decreasing-stepsize SPOS algorithm, stated in Theorem 6.
Theorem 6 Under Assumptions 1, for a decreasing step size hk = h0/(k+ 1), and let the minibatch size in each iteration k be Bk = B0 + [log(k + 1)]100/99, the corresponding W1(µT , ρ∑T−1
k=0 hk )
term is bounded, for some β small enough, as W1(µT , ρ∑T−1
k=0 hk ) ≤ c4β−3Md 3 2 (c7h 3 0 + c8β 3h0/B0 + c9h 2 0β 2)1/2 , (14)
where B0 is the initial minibatch size, and (c4, c7, c8, c9) are some positive constants independent of (M,T, h0).
Note Bk increases at a very low speed, e.g., only by 15 after 105 iterations, thus it would not affect algorithm efficiency. According to Theorem 6, W1(µT , ρ∑T−1
k=0 hk ) would approach zero when
h 1/2 0 M → 0. By directly combining results from Theorem 3–6, one can easily bound the target W1(µT , ν∞). Detailed statements are given in Theorem 15–16 in Section H of the SM.

heading: 4.4 A PITFALL OF SVGD
text: Based on the above analysis, we now formally show a pitfall of SVGD, i.e., particles in SVGD tend to collapse to a local mode under some particular conditions. Inspired by the work on analyzing the granular media equations by Malrieu (2003); Cattiaux et al. (2008), we measure this by calculating the expected distance between particles, called expected particle distance (EPD). Firstly, we bound the EPD for the proposed SPOS algorithm in Theorem 7.
Theorem 7 Under Assumption 1, further assuming every {θ(i)t } with the same initial probability law ρ0 and Γ , Eθ∼ρ0,θ′∼ρ0 [‖θ− θ′‖2] <∞. Choose a β such that λ = mFβ +mK −HFLK > 0.
Then the EPD of SPOS is bounded as: EPD , √∑M
i,j E‖θ (i) t − θ (j) t ‖2 ≤ C1e−2λt + 4 √ d β M λ ,
where C1 = M(M − 1)Γ− 4 √ dβ−1Mλ . Remark 2 There are two interesting cases: i) When C1 > 0, the EPD would decrease to the bound 4 √ dβ−1M/λ along time t. This represents the phenomenon of an attraction force between particles;
ii) When C1 < 0, the EPD would increase to the same bound, which represents the phenomenon of a repulsive force between particles, e.g., when particles are initialized with the same value (Γ = 0), they would be pushed away from each other until the EPD increases to the aforementioned bound.
Intuitively, the EPD for SVGD can be obtained by taking the β → ∞ limit. Corollary 8 formally characterizes the particle-degeneracy phenomenon of SVGD, which has been empirically studied in (Zhuo et al., 2018).
Corollary 8 Under the same conditions of Theorem 7, the EPD in SVGD is bounded as: EPD ,√∑M i,j ‖θ (i) t − θ (j) t ‖2 ≤ C0e−2λt, where C0 = √∑M i,j ‖θ (i) 0 − θ (j) 0 ‖2 and λ = mK −HFLK .
Remark 3 We would like to emphasize two points: 1) In the case of λ ≥ 0, Corollary 8 indicates that particles in SVGD would collapse to a point when t → ∞. In practice, we usually find that particles are trapped in a local mode instead of collapsing. This is due to two reasons: i) numerical errors inject noise into the particles; ii) some particles are out of the concave region of K stated in Assumption 1 in SVGD, which is required for the theory to hold. All these make the empirical EPD not exactly the same as the true particle distance. 2) Corollary 8 also applies when the energy function is nonconvex. Our proof in the SM considers the nonconvex case as well. Consequently, this serves as a strong theoretical motivation to apply SPOS instead of SVGD in deep learning.

heading: 5 NON-ASYMPTOTIC CONVERGENCE ANALYSIS: THE NONCONVEX CASE
text: Since the non-convex case is much more complicated than the convex case, we reply on different assumptions and adopt another distance metric, denoted as B̃, to characterize the convergence behavior of SPOS under the non-convex case. Specifically, B̃(µ, ν) is defined as B̃(µ, ν) , |Eθ∼µ[f(θ)]− Eθ∼ν [f(θ)]| for a known Lf -continuous function f satisfying Assumption 2 below. Note such metric has also been adopted in (Vollmer et al., 2016; Chen et al., 2015). Our analysis considers (T,M, hk) as variables in B̃. In addition, we use {θ̂(i)k }Mi=1 to denote the particles when full gradients are adopted in (9). The distribution of the particles is denoted as µ̂k.
Our high-level idea of bounding B̃(µT , ν∞) is to decompose it as follows:
B̃(µT , ν∞) ≤ B̃(µT , µ̂T ) + B̃(µ̂T , µ̂∞) + B̃(µ̂∞, ρ∞) + B̃(ρ∞, ν∞) (15)
Our second idea is to concatenate the particles at each time into a single vector representation, i.e. defining the new parameter at time t as Θt , [θ (1) t , · · · ,θ (M) t ] ∈ RMd. Consequently, the nonlinear SDE system (8) can be turned into a linear SDE ,which means Θt is driven by the following linear SDE:
dΘt = −FΘ(Θt)dt+ √ 2β−1dW(Md)t , (16)
where FΘ(Θt) , [β−1F (θ (1) t ) − 1M ∑M j=1∇K(θ (1) t − θ (j) t ) + 1 M ∑M j=1K(θ (1) t − θ (j) t )F (θ (j) t ), · · · , β−1F (θ (M) t )− 1M ∑M j=1∇K(θ (M) t −θ (j) t )+ 1 M ∑M j=1K(θ (M) t −θ (j) t )F (θ (j) t )] is a vector function RMd → RMd, andW(Md)t is Brownian motion of dimension Md. Similarly, we can define Θ̂k , [θ̂ (1) k , · · · , θ̂ (M) k ] ∈ RMd for the full-gradient case. Hence, it can be seen that through such a decomposition in (15), the bound related to a nonlinear SDE system (8) reduces to that of a linear SDE. The second term B̃(µ̂T , µ̂∞) reflexes the geometric ergodicity of a linear dynamic system with a numerical method. It is known that even if a dynamic system has an exponential convergence rate to its equilibrium, its corresponding numerical method might not. Our bound for B̃(µ̂T , µ̂∞) is essentially a specification of the result of Mattingly et al. (2002), which has also been applied by Xu et al. (2018). The third term B̃(µ̂∞, ρ∞) reflects the numerical error of a linear SDE, which has been studied in related literature such as (Chen et al., 2015). To this end, we adopt standard assumptions used in the analysis of SDEs (Vollmer et al., 2016; Chen et al., 2015), rephrased in Assumption 2.
Assumption 2 (Assumption in (Vollmer et al., 2016; Chen et al., 2015)) For the linear SDE (16) and a Lipschitz function f , let ψ be the solution functional of the Poisson equation: Gψ(Θ̂k) =
1 M ∑M i=1 f(θ̂ (i) k ) − Eθ∼p(θ|D)[f(θ)], where G denotes the infinite generator of the linear SDE (16). Assume ψ and its up to 4th-order derivatives, Dkψ, are bounded by a function V , i.e., ‖Dkψ‖ ≤ HkVpk for k = (0, 1, 2, 3, 4), Hk, pk > 0. Furthermore, the expectation of V on {Θt} is bounded: supl EVp(Θt) <∞, and V is smooth such that sups∈(0,1) Vp (sΘ + (1− s)Θ′) ≤ H (Vp (Θ) + Vp (Θ′)), ∀Θ,Θ′, p ≤ max{2pk} for H > 0.
Assumption 3 i) F , K and ∇K are LF , LK and L∇k Lipschitz; ii) F satisfies the dissipative property, i.e., 〈F (θ),θ〉 ≥ m ‖θ‖2 − b for some m, b > 0; iii) Remark 1 applies to the nonconvex setting, i.e. sup‖f‖Lip≤1 |Eθ∼µ∞ [f(θ)]− Eθ∼ν∞ [f(θ)]| =W1(ρ∞, ν∞) = O(M −1/2).
Remark 4 Assumption 2 is necessary to control the gap between a numerical solution and the exact solution of an SDE. Specifically, it is used to bound the B̃(µ̂∞, ρ∞) term and the B̃(µT , µ̂T ) term above. Purely relying on the dissipative assumption in Assumption 3 as in non-convex optimization with SG-MCMC (Raginsky et al., 2017; Xu et al., 2018) would induce a bound increasing linearly w.r.t. time t. Thus it is not suitable for our goal. Finally, iii) in Assumption 3 is a mild condition and reasonable because we expect particles to be able to approximate all distributions equally well in the asymptotic limit of t→∞ by ergodicity due to the injected noise. How to remove/replace this assumption is an interesting future work.
Based on the assumptions above, the bounds for B̃(µ̂T , µ̂∞) and B̃(µ̂∞, ρ∞) are summarized below. Theorem 9 Under Assumption 2–3, if we set the stepsize hk = h, we can have the following results:
B̃(µ̂T , µ̂∞) ≤ C2ςσ−Md/2(1 + ςemΘh) exp ( −2mΘThσMd/ log(ς) ) (17)
B̃(µ̂∞, ρ∞) ≤ C3h/β, (18)
where ς = 2LΘ(Mbβ + mΘβ + Md)/mΘ, LΘ = √ 2β−1LF + l ′, mΘ = β−1m − m′, and (σ,C2, C3, l ′,m′) are some positive constants independent of (T, M, h) and σ ∈ (0, 1)
Remark 5 It is seen that in order to make the B̃(µ̂T , µ̂∞) term asymptotically decrease to zero, the number of running iteration T should increase at a rate faster enough to compensate the effect of increasing M . We believe there is room for improving this bound, which is an interesting future work.
Next we bound the B̃(µT , µ̂T ) term related to stochastic gradients. By adapting results from linear SDE (Xu et al., 2018), B̃(µT , µ̂T ) can be bounded with Theorem 10.
Theorem 10 Under Assumptions 2–3, if we set Bk = B and hk = h, B̃(µT , µ̂T ) is bounded as
B̃(µT , µ̂T ) ≤ C5Th(LΘΓ′ +MC4) ((6 + 2Γ′)β/(BM)) 1/2 .
where Γ′ = 2(1 + 1/mΘ)(Mb + 2M2C24 + Md/β) and , (C4, C5) is some positive constant independent of (T, M, h)
Finally, by combining the results from Theorem 9, 10 and iii) in Assumption 3, we arrive at a bound for our target B̃(µT , ν∞), summarized in Theorem 11.
Theorem 11 Under Assumptions 2–3, there exist some positive constants (C2, C3, C4, C5, C6) such that: B̃(µT , ν∞) ≤ C2ςσ−Md/2(1 + ςemΘh) exp ( −2mΘThσMd/ log(ς) ) + C3h/β + C5Th(LΘΓ ′ +MC4) ((6 + 2Γ ′)β/(BM)) 1/2 + C6M −1/2 ,
where σ, ς and Γ′ are the same as those in Theorem 9–10.

heading: 6 EXPERIMENTS
text: In this section, we illustrate the effectiveness of our proposed Bayesian sampling framework on several deep-learning models, including Bayesian learning of deep neural network and Bayesian exploration in deep reinforcement learning. We start with a simple illustrative example.
6.1 BOUNDS ILLUSTRATION WITH A SIMPLE GAUSSIAN EXAMPLE
We follow Chen et al. (2015) and consider a standard Gaussian model where xi ∼ N (θ, 1), θ ∼ N (0, 1). 1000 data samples {xi} are generated, and every minibatch in the stochastic gradient is of size 10. The test function is defined as f(θ) , θ2, with explicit expression for the posterior average. To evaluate the expectations in the bias and MSE, we average over 200 runs with random initializations. The estimation errors are plotted in Figure 2. It is observed that at the beginning, the errors for the ones with less particles decrease faster than those with more particles. This is reflected in the bounds in
Theorem 15, which are dominated by the bound in Theorem 5 (indicating larger M results in larger errors at the beginning). When more running time/iterations are given, the increase of error in Theorem 5 by increasing M essentially cancels out the exponentially-decay term in Theorem 4. According to Theorem 3, the error would eventually decrease with increasing number of particles.

heading: 6.2 BAYESIAN NEURAL NETWORKS FOR REGRESSION
text: We conduct experiments for Bayesian learning of DNNs, where we Bayesian DNNs are used to model weight uncertainty of neural networks, an important topic that has been well explored (HernándezLobato & Adams, 2015; Blundell et al., 2015; Li et al., 2016; Louizos & Welling, 2016). We assign simple isotropic Gaussian priors to
the weights, and perform posterior sampling with different methods. For all methods, we use a RBF kernel K(θ,θ′) = exp(−‖θ − θ′‖22/η2), with the bandwidth set to η = med2/ logM . Here med is the median of the pairwise distance between particles. We use a single-layer BNN for regression tasks. Following Li et al. (2015), 10 UCI public datasets are considered: 100 hidden units for 2 large datasets (Protein and YearPredict), and 50 hidden units for the other 8 small datasets. Following Zhang et al. (2018b), we repeat the experiments 20 times with batchsize 100 for all datasets except for Protein and YearPredict, which we repeat 5 times and once with batchsize 1000. The datasets are randomly split into 90% training and 10% testing. For a fair comparison, we use the same split of data (train, val and test) for SVGD and SPOS. The test results are reported on the best model on the validation set. We adopt the root mean squared error (RMSE) and test log-likelihood as the evaluation criteria. The experimental results are shown in Table 1, from where we can see the proposed SPOS outperforms other methods, achieving state-of-the-art results.

heading: 6.3 BAYESIAN EXPLORATION IN DEEP REINFORCEMENT LEARNING
text: it is well-accepted that RL performance directly measures how well the uncertainty is learned , reflected by the exploration stage. As a result, we apply our method for RL. Following Liu et al. (2017); Zhang et al. (2018a), we define policies in RL with Bayesian neural networks. This naturally introduces uncertainty into action selections under a specific state-action pair, rendering Bayesian explorations to make policy learning more effective.
Specifically, denote the policy as πθ(a | s) parameterized by θ with prior distribution p(θ), where a represent the action variable, and s the state variable. According to Liu et al. (2017), learning the optimal policy corresponds to calculating the following posterior distribution for θ:
q(θ) ∝ exp ( 1
α J(θ)
) p(θ), (19)
where J(θ) denotes the expected cumulative reward under the policy with parameter θ and α a hyperparameter. Consequently, θ could be updated by drawing samples from (19) with the proposed SPOS. We denote this method as SPOS-PG. In addition, when drawing samples with SVGD, the resulting algorithm is called Stein variational policy gradient (SVPG) (Liu et al., 2017). Note in implementation, the term J(θ) can be approximated with REINFORCE (Williams, 1992) or advantage actor critic (Schulman et al., 2015), which we will investigate in our experiments.
We follow the same setting as in (Liu et al., 2017) except using simpler policy-network architectures as Houthooft et al. (2016). We conduct experiments on three classical continuous control tasks are considered: Cartpole Swing-Up, Double Pendulum, and Cartpole. Specifically, the policy is parameterized as a two-layer (25-10 hidden units) neural network with tanh as the activation function. The maximal length of horizon is set to 500. We use a sample size of 10000 for policy gradient estimation, and M = 16, α = 10. For the simplest task, Cartpole, all agents are trained for 100 episodes. For the other two complex tasks, all agents are trained up to 1000 episodes. The average reward versus number of episodes are plotted in Figure 3. It is observed that our proposed SPOS-PG obtains much larger average rewards as well as smaller variance compared to SVPG, though the convergence behaviors are similar in the simplest task Carpole.

heading: 7 CONCLUSION
text: Motivated by the need of effective and efficient Bayesian sampling techniques in modern deep learning, we propose a probability approach for particle-optimization-based sampling that unifies SG-MCMC and SVGD. Notably, for the first time, by analyzing the corresponding nonlinear SDE, we develop non-asymptotic convergence theory for the proposed SPOS framework, a missing yet important theoretical result since the development of SVGD. Within our theoretical framework, a pitfall of SVGD, which has been studied empirically (Zhuo et al., 2018), is formally analyzed. Our theory also indicates the convergence of SPOS to the true posterior distribution in the asymptotic limit of infinite particles and iterations under appropriate conditions. Our theory is of great practice value, as for the first time it provides nonasymptotic theoretical guarantees for the recently proposed particleoptimization-based algorithms such as the SVGD, whose advantages have also been extensively evaluated by experiments on Bayesian learning of DNNs and Bayesian exploration of DRL. There are a number of interesting future works. For example, one might explore more recently developed techniques such as (Cheng et al., 2018) to improve the convergence bound; one can also adopt the SPOS framework for non-convex optimization like what SG-MCMC is used for, and develop corresponding theory to study the convergence property of the algorithm to the global optimum.

heading: A DENSITY FUNCTION OF THE MULTI-MODE DISTRIBUTION IN SECTION 3
text: The negative log-density function of the multi-mode distribution in Section 3 is defined as:
U(θ) , e 3 4θ 2− 32 ∑10 i=1 ci sin( 14πi(θ+4)) ,
where c = (−0.47,−0.83,−0.71,−0.02, 0.24, 0.01, 0.27,−0.37, 0.87,−0.37) is a vector, ci is the i-th element of c.

heading: B GRONWALL LEMMA
text: The Gronwall Lemma plays an important role in parts of our proofs, which is stated in Lemma 12.
Lemma 12 (Gronwall Lemma) Let I denotes an interval of the form [a,+∞) for some a ∈ R. If v(t), defined on I, is differentiable in I and satisfies the following inequality:
v′(t) ≤ β(t)v(t) ,
where β(t) is a real-value continuous function defined on I. Then v(t) can be bounded as: v(t) ≤ v(a) exp( ∫ t a β(s)ds)

heading: C PROOF OF THEOREM 2
text: To prove Theorem 2, we rely on the definition of generalized derivative in Definition 1.
Definition 1 (Generalized Derivative) Let g and φ be locally integrable functions on an open set Ω ⊂ Rd, that is, Lebesgue integrable on any closed bounded set F ⊂ ω. Then φ is the generalized derivative of g with respect to θ on Ω, written as φ = ∂θg, if for any infinitely-differentiable function u with compact support in Ω, we have∫
Ω
g(θ)∂θu(θ)dθ = − ∫
Ω
φ(θ)u(θ)dθ .
Proof The proof relies on further expansions on the definition of generalized derivative on specific functions. Specifically, let the function g in Definition 1 be in a form of g , Gf for two functions G and f . The generalized derivative of (Gf)(v, t), written as ∂θ(Gf)(θ, t), satisfies∫
∂θ(Gf)(θ, t)u(θ)dθ = − ∫ Gf(θ, t)∂θu(θ)dθ (20)
for all differentiable function u(·). In Theorem 2, we want to prove a particle representation of the following PDE:
∂tµt = F1 = −∇θ · (µtF (θt) + (K ∗ µt)µt) , −∂θ(Gf)(θ, t) ⇒ ∫ ∂tµtu(θ)dθ = − ∫ ∂θ(Gf)(θ, t)u(θ)dθ ,
where f(θ, t) = µt. Consequently, we have∫ ∂tf(θ, t)u(θ)dθ = − ∫ ∂θ(Gf)(θ, t)u(θ)dθ (21)
By applying (20) in (21), we have∫ ∂tf(θ, t)u(θ)dθ = − ∫ ∂θ(Gf)(θ, t)u(θ)dθ = ∫ Gf(θ, t)∂θu(θ)dθ .
Since fdθ , µ(dθ, t), we have∫ ∂tµ(dθ, t)u(θ) = ∫ G · µ(dθ, t)∂θu(θ)
⇒ d dt
∫ µ(dθ, t)u(θ) = ∫ G · µ(dθ, t)∂θu(θ)
⇒ d dt Eµ(t)[u(θ)] = Eµ(t)[G · ∂θu(θ)] . (22)
In particle approximation, we have µ(t) = 1M ∑M i=1 δ(θ(i)t )
(θ). For each particle, let u(θ) = θ, (22) reduces to the following equation:
dθ (i) t = G(θ (i) t )dt .
This completes the proof.

heading: D PROOF OF THEOREM 3
text: Note that one challenge in our analysis compared with the analysis for linear SDEs, such as those for SG-MCMC Vollmer et al. (2016); Chen et al. (2015), is how to bound the gap between the original nonlinear PDE (4) and the reduced nonlinear SDE (8). Based on the techniques on analyzing granular media equations in Malrieu (2003); Cattiaux et al. (2008); Durmus et al. (2018), we introduce a nonlinear SDE as an element in-between (6) and (8) like :{
dθ̄t = −β−1F (θ̄t)− EY∼νtK(θ̄t − Y )F (Y ) +∇K ∗ νt(θ̄t) + √
2β−1dW̄t L(θ̄t) = νtdθ
(23)
where L(θ̄t) denotes the probability law of θ̄t, W̄t ∈ Rd is a d-dimensional Brownian motion and Y is a random variable independent of θ̄t and just used here for the sake of clarity. In order to match the SDE system (8) of the particles {θ(i)t }Mi=1, we duplicate (23) M times, each endowing with an exact solution θ̄(i)t indexed by i. The distribution of each particles {θ̄ (i) t }Mi=1 is νt and the corresponding W̄(i)t can be set exactly the same as theW (i) t in (8):{
dθ̄ (i) t = −β−1F (θ̄ (i) t )− EYi∼νtK(θ̄ (i) t − Yi)F (Yi) +∇K ∗ νt(θ̄ (i) t ) + √ 2β−1dW̄(i)t L(θ̄(i)t ) = νtdθ (24)
where Yi is a random variable independent of θ̄ (i) t and just used here for the convenience of the proof.
Proof [Proof of Theorem 3] Firstly we have d ( θ
(i) t − θ̄ (i) t
) =− β−1 ( F (θ
(i) t )− F (θ̄ (i) t ) ) dt+ 1
M M∑ j [ ∇K(θ(i)t − θ (j) t )−∇K ∗ νt(θ̄ (i) t ) ] dt
− 1 M M∑ j ( F (θ (j) t )W (θ (i) t − θ (j) t )− EYj∼νtF (Yj)W (θ̄ (i) t − Yj) ) dt
⇒ d ( M∑ i ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2 ) = 2 M M∑ i,j (Aij(t) +Bij(t) + Cij(t) + Fij(t) +Gij(t) +Hij(t))dt
(25) where Aij(t) = −β−1 ( F (θ (i) t )− F (θ̄ (i) t ) ) · ( θ (i) t − θ̄ (i) t ) Bij(t) = ( ∇K(θ(i)t − θ (j) t )−∇K(θ̄ (i) t − θ̄ (j) t ) ) · ( θ (i) t − θ̄ (i) t )
Cij(t) = ( ∇K(θ̄(i)t − θ̄ (j) t )−∇K ∗ νt(θ̄ (i) t ) ) · ( θ (i) t − θ̄ (i) t ) Fij(t) = − ( F (θ (j) t )K(θ (i) t − θ (j) t )− F (θ̄ (j) t )K(θ (i) t − θ (j) t ) ) · ( θ (i) t − θ̄ (i) t
) Gij(t) = − ( F (θ̄ (j) t )K(θ (i) t − θ (j) t )− F (θ̄ (j) t )K(θ̄ (i) t − θ̄ (j) t ) ) · ( θ (i) t − θ̄ (i) t
) Hij(t) = − ( F (θ̄ (j) t )K(θ̄ (i) t − θ̄ (j) t )− EYj∼νtF (Yj)K(θ̄ (i) t − Yj) ) · ( θ (i) t − θ̄ (i) t
) For the Aij(t) term, according to the i) in Assumption 1 for F , we have∑
ij Aij(t) = − ∑ ij β−1 ( F (θ (i) t )− F (θ̄ (i) t ) ) · ( θ (i) t − θ̄ (i) t ) ≤ −β−1mFM
∑ i ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2 For the Bij(t) term, applying the concave condition for K and the oddness of ∇K in Assumption 1, we have∑
ij
Bij(t) = M∑ ij ( ∇K(θ(i)t − θ (j) t )−∇K(θ̄ (i) t − θ̄ (j) t ) ) · ( θ (i) t − θ̄ (i) t )
= 1
2 M∑ ij ( ∇K(θ(i)t − θ (j) t )−∇K(θ̄ (i) t − θ̄ (j) t ) ) · ( θ (i) t − θ̄ (i) t − (θ (j) t − θ̄ (j) t ) )
≤ −1 2 mK M∑ ij ∥∥∥θ(i)t − θ̄(i)t − (θ(j)t − θ̄(j)t )∥∥∥2 ≤ 0 For the Cij(t) term, we have
E ∑ j Cij(t) (1) ≤ ( E ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2)1/2 E ∥∥∥∥∥∥ ∑ j ( ∇K(θ̄(i)t − θ̄ (j) t )−∇K ∗ νt(θ̄ (i) t ) )∥∥∥∥∥∥ 2  1/2
(2) = ( E ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2)1/2 ∑ j E ( ∇K(θ̄(i)t − θ̄ (j) t )−∇K ∗ νt(θ̄ (i) t ) )21/2
(3) ≤ H∇K √ 2M ( E ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2)1/2
where the (1) is obtained by applying the Cauchy-Schwarz inequality and (2) by the fact that E ( K(θ̄
(i) t − θ̄ (j) t )−K ∗ νt(θ̄ (i) t ) )
= 0. We can tune the bandwidth of the RBF kernel to make ‖∇K‖ ≤ H∇K . Hence (3) is obtain by the boundedness of∇K(θ).
Similarly, since K ≤ 1, we have the following result for Hij(t) term, E ∑ j Hij(t) ≤ ( E ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2)1/2 E ∥∥∥∥∥∥ ∑ j ( F (θ̄ (j) t )K(θ̄ (i) t − θ̄ (j) t )− EYj∼νtF (Yj)K(θ̄ (i) t − Yj) )∥∥∥∥∥∥ 2  1/2
= ( E ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2)1/2 ∑ j E ( F (θ̄ (j) t )K(θ̄ (i) t − θ̄ (j) t )− EYj∼νtF (Yj)K(θ̄ (i) t − Yj) )21/2
≤HF √ 2M ( E ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2)1/2
For the Fij(t) and Gij(t) terms, we have:∑ ij Fij(t) = − ∑ ij ( F (θ (j) t )K(θ (i) t − θ (j) t )− F (θ̄t (j) )K(θ (i) t − θ (j) t ) ) · ( θ (i) t − θ̄ (i) t ) ≤ ∑ ij LF
∥∥∥θ(j)t − θ̄(j)t ∥∥∥ ∥∥∥θ(i)t − θ̄(i)t ∥∥∥ ≤ 2LFM
∑ i ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2 . ∑ ij Gij(t) = ∑ ij ( F (θ̄ (j) t )K(θ (i) t − θ (j) t )− F (θ̄ (j) t )K(θ̄ (i) t − θ̄ (j) t ) ) · ( θ (i) t − θ̄ (i) t
) ≤ HFLK
∑ ij ∥∥∥θ(i)t − θ̄(i)t − (θ(j)t − θ̄(j)t )∥∥∥∥∥∥θ(i)t − θ̄(i)t ∥∥∥ ≤ 3HFLKM
∑ i ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2
We denote γi(t) , E ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2. Due to the exchangeability of the particles, γi(t) are the same
for all the particles, denoted as γ(t). Then according to (25), we have
γ′(t) ≤ −2λ1γ(t) + H∇K
√ 2 +HF √ 2√
M
√ γ(t) .
where λ1 = β−1mF − 3HFLK − 2LF
⇒ ( √ γ(t)− (H∇K +HF )/ √ 2√
M(β−1 − 3HFLK − 2LF ) )′ ≤ −λ1(
√ γ(t)− (H∇K +HF )/ √ 2√
M(β−1 − 3HFLK − 2LF ) )
Note that θ(i)t and θ̄ (i) t are initialized with the same initial distribution µ0 = ν0. In the definition of θ̄ (i) t , there is no restriction on how the initial value is set. As a result, we can set θ (i) 0 to be identical to θ̄(i)0 , leading to γ(0) = 0. Then according to the Gronwall Lemma, we have√ γ(t) ≤ (H∇K +HF )/ √ 2√
M(β−1 − 3HFLK − 2LF )
Hence, there exist some positive constant (c1, c2) such that:
W1(ρt, νt) (1)
≤ W2(ρt, νt) (2) ≤ √ E ∥∥∥θ(i)t − θ̄(i)t ∥∥∥2 (3)≤ c1√
M(β−1 − c2) , (26)
where (1) holds due to the relationship between W1 and W2 metric Givens & Shortt (1984), (2) due to the definition of W2 and (3) due to the result from the previous proof.

heading: E PROOF OF THEOREM 4
text: Proof [Proof of Theorem 4] Firstly, what we aim at is W1(νt, ν∞) ≤ c3 exp (−2λ1t) in this theorem. According to the relationship between W1 and W2 metric Givens & Shortt (1984), once we bound W2(νt, ν∞) as W2(νt, ν∞) ≤ c3 exp (−2λ1t), we will finish our proof.
Next, look at the equation (8): If we set the initial distribution of each particle to be ν0, which means ρ0 = L(θ(i)0 ) = ν0, we will derive M particles denoted as {θ(i)t,1}Mi=1. We denote the distribution of each θ (i) t,1 at t as ρt,1. If we set the initial distribution of each particle to be ν∞, which means ρ0 = L(θ(i)0 ) = ν∞, we will derive M particles denoted as {θ(i)t,2}Mi=1. We denote the distribution of each θ (i) t,2 at t as ρt,2.
Since we need to bound W2(νt, ν∞), we make the following decomposition:
W2(νt, ν∞) ≤W2(νt, ρt,1) +W2(ρt,1, ρt,2) +W2(ρt,2, ν∞) . (27)
Note that ρ0,1 = ν0 and ρ0,2 = ν∞. Then, according to (26), we have
W2(νt, ρt,1) ≤ c1√
M(β−1 − c2) and W2(ρt,2, ν∞) ≤ c1√ M(β−1 − c2)
Now we need to focus on the term W2(ρt,1, ρt,2). Since W2(ρt,1, ρt,2) ≤ E (∥∥∥θ(i)t,1 − θ(i)t,2∥∥∥2) ,
r(t), we will derive a bound for E (∥∥∥θ(i)t,1 − θ(i)t,2∥∥∥2) in the following. We have
d ( θ
(i) t,1 − θ (i) t,2
) =− β−1 ( F (θ
(i) t,1)− F (θ (i) t,2) ) dt
+ 1
M M∑ j [ ∇K(θ(i)t,1 − θ (j) t,1 )−∇K(θ (i) t,2 − θ (j) t,2 ) ] dt
− 1 M M∑ j ( F (θ (j) t,1 )K(θ (i) t,1 − θ (j) t,1 )− F (θ (j) t,2 )K(θ (i) t,2 − θ (j) t,2 ) ) dt
⇒ d ( M∑ i ∥∥∥θ(i)t,1 − θ(i)t,2∥∥∥2 ) = 2 M M∑ i,j (ξ1ij(t) + ξ 2 ij(t) + ξ 3 ij(t) + ξ 4 ij(t))dt
where ξ1ij(t) = −β−1 ( F (θ (i) t,1)− F (θ (i) t,2) ) · ( θ (i) t,1 − θ (i) t,2 ) ξ2ij(t) = ( ∇K(θ(i)t,1 − θ (j) t,1 )−∇K(θ (i) t,2 − θ (j) t,2 ) ) · ( θ (i) t,1 − θ (i) t,2
) ξ3ij(t) = − ( F (θ (j) t,1 )K(θ (i) t,1 − θ (j) t,1 )− F (θ (j) t,2 )K(θ (i) t,1 − θ (j) t,1 ) ) · ( θ (i) t,1 − θ (i) t,2
) ξ4ij(t) = − ( F (θ (j) t,2 )K(θ (i) t,1 − θ (j) t,1 )− F (θ (j) t,2 )K(θ (i) t,2 − θ (j) t,2 ) ) · ( θ (i) t,1 − θ (i) t,2 )
For the ξ1ij(t) terms, according to the i) in Assumption 1 for F , we have∑ ij ξ1ij(t) = − ∑ ij β−1 ( F (θ (i) t,1)− F (θ (i) t,2) ) · ( θ (i) t,1 − θ (i) t,2 ) ≤ −β−1mFM
∑ i ∥∥∥θ(i)t,1 − θ(i)t,2∥∥∥2 .
For the ξ2ij(t) term, applying the concave condition for K and the oddness of∇K in Assumption 1, we have∑
ij
ξ2ij(t) = M∑ ij ( ∇K(θ(i)t,1 − θ (j) t,1 )−∇K(θ (i) t,2 − θ (j) t,2 ) ) · ( θ (i) t,1 − θ (i) t,2 )
= 1
2 M∑ ij ( ∇K(θ(i)t,1 − θ (j) t,1 )−∇K(θ (i) t,2 − θ (j) t,2 ) ) · ( θ (i) t,1 − θ (i) t,2 − (θ (j) t,1 − θ (j) t,2 ) )
≤ −1 2 mK M∑ ij ∥∥∥θ(i)t,1 − θ(i)t,2 − (θ(j)t,1 − θ(j)t,2 )∥∥∥2 ≤ 0 . For the ξ3ij(t) terms, after applying the LF -Lipschitz property for F and K ≤ 1, we have∑
ij ξ3ij(t) = ∑ ij − ( F (θ (j) t,1 )K(θ (i) t,1 − θ (j) t,1 )− F (θ (j) t,2 )K(θ (i) t,1 − θ (j) t,1 ) ) · ( θ (i) t,1 − θ (i) t,2 ) ≤ ∑ ij LF
∥∥∥θ(j)t,1 − θ(j)t,2∥∥∥∥∥∥θ(i)t,1 − θ(i)t,2∥∥∥ ≤ 2LFM
∑ i ∥∥∥θ(i)t,1 − θ(i)t,2∥∥∥2 . For the ξ4ij(t) terms :∑
ij ξ4ij(t) = − ∑ ij ( F (θ (j) t,2 )K(θ (i) t,1 − θ (j) t,1 )− F (θ (j) t,2 )K(θ (i) t,2 − θ (j) t,2 ) ) · ( θ (i) t,1 − θ (i) t,2 ) ≤ HFLK
∑ ij ∥∥∥θ(i)t,1 − θ(i)t,2 − (θ(j)t,1 − θ(j)t,2 )∥∥∥∥∥∥θ(i)t,1 − θ(i)t,2∥∥∥ ≤ 3HFLKM
∑ i ∥∥∥θ(i)t,1 − θ(i)t,2∥∥∥2 .
Now we have
r′(t) ≤ −2(β−1mF − 3HFLK − 2LF )r(t) .
According to the Gronwall lemma,
r(t) ≤ r(0)e−2λ1t,
where λ1 = β−1mF − 3HFLW − 2LF .
Consequently, there exists some positive constant c3 such that
W2(ρt,1, ρt,2) ≤ c3e−2λ1t
Then we have
W2(νt, ν∞) ≤ c3e−2λ1t + c1√
M(β−1 − c2) + c1√ M(β−1 − c2)
However, it worth noting that νt is the solution of (6) which has nothing to do with the number of particles, M . Then let M → ∞, we can derive that W2(νt, ν∞) ≤ c3e−2λ1t. Now we finish our proof.

heading: F PROOF OF THEOREM 5
text: To bound the W1(µT , ρ∑T−1 k=0 hk
) term, note the original SDE driving the particles {θ(i)t } in (8) is a nonlinear SDE, which is hard to deal with. Fortunately, (8) can be turned into a linear SDE by concatenating the particles at each time into a single vector representation, i.e., by defining the new parameter at time t as Θt , [θ (1) t , · · · ,θ (M) t ] ∈ RMd. Consequently, Θt is driven by the following linear SDE:
dΘt = −FΘ(Θt)dt+ √ 2β−1dW(Md)t , (28)
where FΘ(Θt) , [β−1F (θ (1) t ) − 1M ∑M j=1∇K(θ (1) t − θ (j) t ) + 1 M ∑M j=1K(θ (1) t − θ (j) t )F (θ (j) t ), · · · , β−1F (θ (M) t )− 1M ∑M j=1∇K(θ (M) t −θ (j) t )+ 1 M ∑M j=1K(θ (M) t −θ (j) t )F (θ (j) t )] is a vector function RMd → RMd, andW(Md)t is Brownian motion of dimension Md.
Now we define the F(q)Θ(Θt) , [β−1Fq(θ (1) t )− 1MN ∑M j=1∇K(θ (1) t −θ (j) t )+ 1 M ∑M j=1K(θ (1) t − θ (j) t )Fq(θ (j) t ), · · · , β−1Fq(θ (M) t ) − 1MN ∑M j=1∇K(θ (M) t − θ (j) t ) + 1 M ∑M j=1K(θ (M) t − θ (j) t )Fq(θ (j) t )]. We can verify that FΘ(Θt) = ∑N q=1 F(q)Θ(Θt).
Then we define Θk , [θ (1) k , · · · , θ (M) k ] and G Θ Ik , N Bk ∑ q∈Ik F(q)Θ(Θk). We can verify that the following result holds:
Θk+1 = Θk − β−1GΘIkhk + √ 2β−1hkΞk , (29)
where Ξk ∼ N (0, IMd×Md). Now we reach the conclusion that Θk of (29) is accutually the numerical solution of the SDE (28) via stochastic gradients.
We denote the distribution of Θk as µΘk and the distribution of Θt as ρ Θ t . Before proceeding to our theoretical results, we need to present the following Lemmas which is very important in our proof.
Lemma 13 W1(µk, ρt) ≤ 1√MW1(µ Θ k , ρ Θ t )
Proof [Proof of Lemma 13] Let us recall the definition of W1 metric and its Kantorovich-Rubinstein duality Villani (2008), i.e. W1(µ, ν) , sup‖g‖lip≤1 |Eθ∼µ[g(θ)]− Eθ∼ν [g(θ)]|. We can prove the fact that if g(θ) : Rd → R is a Lg-Lipschitz function in Rd, the gΘ(Θ), defined as gΘ(Θ) = 1√ M ∑M i g(θ (i)), is a Lg-Lipschitz function in RMd, where Θ , [θ(1), · · · ,θ(M)]:
‖gΘ(Θ1)− gΘ(Θ2)‖ ≤ 1√ M M∑ i=1 ‖g(θ(i)1 )− g(θ (i) 2 )‖ ≤ Lg√ M M∑ i=1 ‖θ(i)1 − θ (i) 2 ‖
≤ Lg√ M
√ M √√√√ M∑ i=1 ‖θ(i)1 − θ (i) 2 ‖2 = Lg‖Θ1 −Θ2‖
Then we have:
1
M M∑ i=1 ∣∣∣Eθ(i)k ∼µk [g(θ(i)k )]− Eθ(i)t ∼ρt [g(θ(i)t )]∣∣∣ (1)= 1√M ∣∣∣∣∣ 1√M M∑ i=1 (E θ (i) k ∼µk [g(θ (i) k )]− Eθ(i)t ∼ρt [g(θ (i) t )]) ∣∣∣∣∣ =
1√ M |EΘk∼µk [gΘ(Θk)]− EΘt∼ρt [gΘ(Θt)]|
The (1) holds since E θ (1) k ∼µk [g(θ (1) k )] = · · · = Eθ(M)k ∼µk [g(θ (M) k )] for all the particles θ (i) k and E θ (1) t ∼ρt [g(θ (1) t )] = · · · = Eθ(M)t ∼ρt [g(θ (M) t )] for all the particles θ (i) t . Then according to the
definition of W1 metric, we derive that
W1(µk, ρt) = sup ‖g‖lip≤1
1
M M∑ i=1 ∣∣∣Eθ(i)k ∼µk [g(θ(i)k )]− Eθ(i)t ∼ρt [g(θ(i)t )]∣∣∣ = 1√ M sup ‖g‖lip≤1 |EΘk∼µk [gΘ(Θk)]− EΘt∼ρt [gΘ(Θt)]| = 1√ M sup ‖gΘ‖lip≤1 |EΘk∼µk [gΘ(Θk)]− EΘt∼ρt [gΘ(Θt)]|
≤ 1√ M W1(µ Θ k , ρ Θ t )
Lemma 14 Assuming F (0) = 0. If F in (9) is Lipschitz with constant LF , and satisfies the dissipative property that 〈F (θ),θ〉 ≥ mF ‖θ‖2 − b. Then FΘ in (28) is Lipschitz-continuous with constant √ 2β−1LF + l
′ and satisfies 〈FΘ(Θ),Θ〉 ≥ (β−1mF −m′) ‖Θ‖2 − β−1Mb, where l′ and m′ are some positive constants.
Proof [Proof of Lemma 14]
‖FΘ(Θ1)− FΘ(Θ2)‖ = √√√√ M∑ i ‖ω1i + ω2i + ω3i ‖ 2 ≤ √√√√ M∑ i (‖ω1i ‖+ ‖ω2i ‖+ ‖ω3i ‖) 2
where
‖ω1i ‖ = ‖β−1F (θ (i) 1 )− β−1F (θ (i) 2 )‖
≤ β−1LF ‖θ(i)1 − θ (i) 2 ‖
‖ω2i ‖ = ‖ 1
M ( M∑ j K(θ (i) 1 − θ (j) 1 )F (θ (j) 1 )− M∑ j K(θ (i) 2 − θ (j) 2 )F (θ (j) 2 ))‖
≤ ‖ 1 M ( M∑ j K(θ (i) 1 − θ (j) 1 )F (θ (j) 1 )− M∑ j K(θ (i) 2 − θ (j) 2 )F (θ (j) 1 ))‖
+ ‖ 1 M ( M∑ j K(θ (i) 2 − θ (j) 2 )F (θ (j) 1 )− M∑ j K(θ (i) 2 − θ (j) 2 )F (θ (j) 2 ))‖
≤ (2LKHF + LF )‖θ(i)1 − θ (i) 2 ‖
‖ω3i ‖ = ‖ − 1
M ( M∑ j ∇K(θ(i)1 − θ (j) 1 )− M∑ j ∇K(θ(i)2 − θ (j) 2 ))‖
≤ L∇K M ( M∑ j ‖θ(i)1 − θ (i) 2 − (θ (j) 1 − θ (j) 2 )‖)
≤ L∇K(‖θ(i)1 − θ (i) 2 ‖+
1
M M∑ j ‖θ(j)1 − θ (j) 2 ‖)
It is easy to verify that there exits some positive constant l′ such that
‖FΘ(Θ1)− FΘ(Θ2)‖ = √√√√ M∑ i ‖ω1i + ω2i + ω3i ‖ 2
≤ √√√√ M∑ i 2(β−1LF + 2LKHF + LF + L∇K)2‖θ(i)1 − θ (i) 2 ‖2 + 2 M∑ j ‖θ(j)1 − θ (j) 2 ‖2
≤ √ 2(β−1LF + 2LKHF + LF + L∇K)2 + 2 √√√√ M∑ i ‖θ(i)1 − θ (i) 2 ‖2
≤ ( √ 2β−1LF + l ′) √√√√ M∑ i ‖θ(i)1 − θ (i) 2 ‖2 = ( √ 2β−1LF + l ′)‖Θ1 −Θ2‖
Next, we have 〈FΘ(Θ),Θ〉
= M∑ i β−1F (θ(i))θ(i) + 1 M M∑ j K(θ(i) − θ(j))F (θ(j))θ(i) − 1 M M∑ j ∇K(θ(i) − θ(j))θ(i) 
Notice : M∑ i β−1F (θ(i))θ(i) ≥ β−1mF M∑ i ‖θ(i)‖2 − β−1β−1Mb = β−1mF ‖Θ‖2 − β−1Mb Since we have assumed F (0) = 0, we have:
M∑ i 1 M M∑ j K(θ(i) − θ(j))F (θ(j))θ(i) ≥− 1 M M∑ i M∑ j LF ‖θ(i)‖‖θ(j)‖
≥ − 2LF M∑ i=1 ‖θ(i)‖2 = −2LF ‖Θ‖2
Since∇K is an odd function, we have: M∑ i 1 M M∑ j ∇K(θ(i) − θ(j))θ(i) ≥ − M∑ i 1 M M∑ j L∇k‖θ(i) − θ(j)‖‖θ(i)‖
≥ −3L∇K M∑ i ‖θ(i)‖2 = −3L∇K‖Θ‖2
As a result, we can derive the following result: 〈FΘ(Θ),Θ〉 ≥ (β−1m− 2LF − 3L∇K)‖Θ‖2 − β−1Mb
Now it is ready to prove Theorem 5. It worth noting that after assuming F (0) = 0, the first bullet in Assumption 1 recovers the dissipative assumption as 〈F (θ),θ〉 ≥ mF ‖θ‖2. Proof Next we use Lemma C.5 in Xu et al. (2018) to verify that FΘ satisfies the assumptions in Raginsky et al. (2017) by setting δ = a′B with a
′ a positive constant and B the size of the random set I. Let µΘk := L(Θk) and ρΘt := L(Θt). Now we can borrow the result of Lemma 3.6 in Raginsky et al. (2017). The relative entropy DKL(µΘk ‖ρΘkh) satisfies:
DKL(µ Θ k ‖ρΘkh) ≤ (A0β a′ B +A1h)kh
with
A0 = ( 2( LF β + l′)2 ( a2 + 2(1 ∨
1
β−1mF −m′ )(2a21 +
Md
β )
) + a21 ) A1 =6(
LF β + l′)2(βA0 +Md)
and a1, a2 are some positive constants. When the β is small enough, there exist some positive constants a3, a4 such that
A0 ≤ a3 Md
β3 , A1 ≤ a4
Md
β4
Similar to the proof of Lemma 14, it is easy to verify that there exists some positive constant a5 such that 〈FΘ(Θ1)− FΘ(Θ2),Θ1 −Θ2〉 ≥ (β−1mF − a5)‖Θ1 −Θ2‖2. Notice, when β is small enough, (28) satisfies the conditions of Proposition 4.2 in Cattiaux et al. (2008). Hence, there exits
some positive constant C such that W1(µΘk , ρΘkh) ≤ C √ DKL(µΘk ‖ρΘkh) .
According to Corollary 4 and Lemma 8 in Bolley & Villani (2005), we can derive an explicit expression for C :
C ≤ a6β−1Md ,
when β is small enough and a6 is some positive constant.
Applying the Lemma 13, we have
W1(µk, ρkh) ≤ 1√ M W1(µ Θ k , ρ Θ kh) ≤ a6Md 3 2 β−3(a3a ′β2B−1 + a4h) 1 2 k 1 2h 1 2
Let k = T and we can finish the poof.

heading: G PROOF OF THEOREM 6
text: Proof Our proof is based on the proof of Lemma 3.6 in Raginsky et al. (2017) with some modifications. Firstly, adopt the same notations in the Section F and we get the following update:
Θk+1 = Θk − β−1GΘIkhk + √ 2β−1hkΞk , (30)
where Ξk ∼ N (0, IMd×Md) and hk = h0k+1 . We assume E(G Θ Ik) = FΘ(Θk), ∀Θ ∈ R Md , which is a general assumption due to the way that we choose the minibatch Ik. We need to define q(t), which will be used in the following proof:
q(t) = {k ∈ R| k−1∑ i=0 hi ≤ t < k∑ i=0 hi} .
Furthermore, we define ∑−1 i=0 hi , 0 and ∑0 i=0 hi , h0 here just used for the convenience of statement in the following.
Now we focus on the following continuous-time interpolation of Θk:
Θ(t) =Θ0 − ∫ t
0
G̃ΘI(s) Θ(q(s)−1∑ i=0 hi) ds+√ 2 β ∫ t 0 W(Md)s ,
where I(s) ≡ Ik for t ∈ [∑k−1 i=0 hi, ∑k i=0 hi ) , G̃ΘI(s)(Θ) , N B(s) ∑ q∈I(s) F(q)Θ(Θ) and B(s) is
the size of the minibatch I(s). And for each k, Θ( ∑k−1 i=0 hi) and Θk have the same probability law ρΘk . Since Θ(t) is not a Markov process, we define the following Itô process which has the same one-time marginals as Θ(t)
Λ(t) = Θ0 − ∫ t
0
Gs (Λ(s)) ds+
√ 2
β ∫ t 0 W(Md)s
where Gt(x) := E G̃ΘI(t) Θ(q(t)−1∑
i=0
hi)
 |Θ(t) = x  .
Let PtΛ := L (Λ(s) : 0 ≤ s ≤ t) and P t Θ := L (Θ(s) : 0 ≤ s ≤ t). According to the proof of lemma 3.6 in Raginsky et al. (2017), we can derive a similar result for the relative entropy of PtΛ and PtΘ:
DKL(P t Λ ‖P t Θ) =− ∫ dPtΛ log
dPtΛ dPtΘ
= β
4 ∫ t 0 E‖FΘ(Λ(s))−Gs (Λ(s)) ‖2ds
= β
4 ∫ t 0 E‖FΘ(Θ(s))−Gs (Θ(s)) ‖2ds
The last line follows that L(Θ(s)) = L(Λ(s)), ∀s. In the following proof, we will let t = ∑k−1 i=0 hi for some k ∈ R. Now we can use the martingale property of Itô integral to derive:
DKL(P ∑k−1 i=0 hi Λ ‖P ∑k−1 i=0 hi Θ )
= β
4 k−1∑ j=0 ∫ ∑j i=0 hi∑j−1 i=0 hi E‖FΘ(Θ(s))−Gs (Θ(s)) ‖2ds
≤ β 2 k−1∑ j=0 ∫ ∑j i=0 hi∑j−1 i=0 hi E‖FΘ(Θ(s))− FΘ(Θ( q(s)−1∑ i=0 hi))‖2ds
+ β
2 k−1∑ j=0 ∫ ∑j i=0 hi∑j−1 i=0 hi E‖FΘ(Θ( q(s)−1∑ i=0 hi))− G̃ΘI(s) Θ(q(s)−1∑ i=0 hi)  ‖2ds ≤ βL2FΘ
2 k−1∑ j=0 ∫ ∑j i=0 hi∑j−1 i=0 hi E‖Θ(s)−Θ( q(s)−1∑ i=0 hi)‖2ds (31)
+ β
2 k−1∑ j=0 ∫ ∑j i=0 hi∑j−1 i=0 hi E‖FΘ(Θ( q(s)−1∑ i=0 hi))− G̃ΘI(s) Θ(q(s)−1∑ i=0 hi)  ‖2ds (32) For the first part (31), we consider some s ∈ [ ∑j−1 i=0 hi, ∑j i=0 hi). The following equation holds:
Θ(s)−Θ( j−1∑ i=0 hi) =− (s− j−1∑ i=0 hi)G Θ Ij + √ 2/β(W(Md)s −W (Md)∑j−1 i=0 hi )
= −(s− j−1∑ i=0 hi)G Θ Ij + (s− j−1∑ i=0 hi)(FΘ(Θj)−GΘIj )
+ √
2/β(W(Md)s −W (Md)∑j−1 i=0 hi )
Thus, we can use Lemma 3.1 and 3.2 in Raginsky et al. (2017), and Lemma C.5 in Xu et al. (2018) to get the following result:
E‖Θ(s)−Θ( j−1∑ i=0 hi)‖2
≤3 h0 2
(j + 1)2 E‖GΘIj‖
2 + 3 h0
2
(j + 1)2 E‖FΘ(Θj)−GΘIj‖
2 + 6h0Md
β(j + 1)
≤12 h0 2
(j + 1)2 max 0≤j≤k−1 (L2FΘE‖Θj‖
2 + b1) + 6h0Md
β(j + 1)
where b1 is some positive constant.
Consequently, the first part above (31) can be bounded as:
βL2FΘ 2 k−1∑ j=0 ∫ ∑j i=0 hi∑j−1 i=0 hi E‖Θ(s)−Θ( q(s)−1∑ i=0 hi)‖2ds
≤ βL2FΘ
2 k−1∑ j=0 [ 12 h0 3 (j + 1)3 max 0≤j≤K−1 (L2FΘE‖Θj‖ 2 + b1) + 6h0 2Md β(j + 1)2 ]
≤π2βL2FΘh0 3 max 0≤j≤K−1 (L2FΘE‖Θj‖
2 + b1) + π2L2FΘh0 2Md
2 ,
where the last line follows from the fact that
k−1∑ j=0
1 (j + 1)3 ≤ k−1∑ j=0
1 (j + 1)2 ≤ ∞∑ j=0
1 (j + 1)2 = π2 6 .
According to Lemma C.5 in Xu et al. (2018), the second part (32) can be bounded as follows:
β
2 k−1∑ j=0 ∫ ∑j i=0 hi∑j−1 i=0 hi E‖FΘ(Θ( q(s)−1∑ i=0 hi))− G̃ΘI(s) Θ(q(s)−1∑ i=0 hi)  ‖2ds =
k−1∑ j=0 βh0 2(j + 1) E‖FΘ(Θj)−GΘIj‖ 2
≤βh0 max 0≤j≤k−1 (L2FΘE‖Θj‖ 2 + b1)  4 B0 + k−1∑ j=1
4
(j + 1)(B0 + log 100 99 (j + 1))  ≤βh0 max
0≤j≤k−1 (L2FΘE‖Θj‖ 2 + b1)  4 B0 + k−1∑ j=1
4
(j + 1)log 100 99 (j + 1)  ≤(b2 + 4
B0 )βh0 max 0≤j≤k−1 (L2FΘE‖Θj‖ 2 + b1) ,
where the last line follows from the fact that when r > 1,
k−1∑ j=1
4 (j + 1) logr(j + 1) ≤ ∞∑ j=1
4 (j + 1) logr(j + 1) ≤ 4 log
1−r 2
r − 1 .
Denote µΘk := L(Θk) and ρΘt := L(Θt). Due to the data-processing inequality for the relative entropy, we have
DKL(µ Θ k ‖ρΘ∑k−1
i=0 hi ) ≤ DKL(P
∑k−1 i=0 hi Λ ‖P ∑k−1 i=0 hi Θ )
≤ π2βL2FΘh0 3 max 0≤j≤k−1 (L2FΘE‖Θj‖
2 + b1) + π2L2FΘh0 2Md
2
+ (b2 + 4
B0 )βh0 max 0≤j≤k−1 (L2FΘE‖Θ̃j‖ 2 + b1)
≤ (π2βL2FΘh0 3 + b2βh0 +
4
B0 βh0) max 0≤j≤k−1 (L2FΘE‖Θj‖
2 + b1) + π2L2FΘh0 2Md
2 .
Lemma 3.2 in Raginsky et al. (2017) has provided a uniform bound to max0≤j≤k−1(L2FΘE‖Θj‖ 2 + b1). Hence we can tell that DKL(P ∑k−1 i=0 hi Λ ‖P ∑k−1 i=0 hi
Θ ) would not increase w.r.t. k. This is a nice property that the fixed-step-size SPOS does not endow. Since LFΘ = √ 2β−1LF + l ′, it is easy to
verify that when β is small enough, there exists some positive constants b3, b4, b5 and b6 such that:
DKL(µ Θ k ‖ρΘ∑k−1
i=0 hi )
≤(π2βL2FΘh0 3 + b2βh0 + 4βh0 B0 ) max 0≤j≤K−1 (L2FΘE‖Θ̃j‖ 2 + b1) +
π2L2FΘh0 2Md
2
≤(b3h30 + b4β 3h0 B0 + b5h 2 0β 2) Md β4 .
Similar to the proof of Theorem 5, we can bound the W1(µΘk ‖ρΘ∑k−1 i=0 hi
) term with Corollary 4, Lemma 8 in Bolley & Villani (2005) and Proposition 4.2 in Cattiaux et al. (2008). Specifically, when β is small enough, there exist some positive constant a6 such that:
W1(µ Θ k ‖ρΘ∑k−1
i=0 hi ) ≤ a6(
Md β ) √ DKL(µΘk ‖ρΘ∑k−1 i=0 hi )
≤ a6β−3M 3 2 d 3 2 (b3h 3 0 + b4β 3h0 B0 + b5h 2 0β 2) 1 2 .
According to Lemma 13, we have
W1(µk, ρkh) ≤ 1√ M W1(µ Θ k ‖ρΘ∑k−1 i=0 hi ) = a6β −3Md 3 2 (b3h 3 0 + b4β 3h0 B0 + b5h 2 0β 2) 1 2
Let k = T and we can finish the proof.

heading: H DETAILED STATEMENTS OF NONASYMPTOTIC CONVERGENCE UNDER THE CONVEX CASE
text: We give detailed statements of non-asymptotic convergence of SPOS under the convex case, which can be proved by directly combining results from Theorem 3–6, thus they are omitted for simplicity. We consider a fixed-stepsize and a decreasing-stepsize cases in Theorem 15 and Theorem 16, respectively.
Theorem 15 (Fixed Stepsize) Under Assumption 1, if we set hk = h0 and Bk = B0, we can bound the W1(µT , ν∞) as
W1(µT , ν∞) ≤ c1√
M(β−1 − c2) + c3 exp
{ −2 ( β−1mF − 3HFLK − 2LF ) Th }
+ c3Md 3 2 β−3(c4β 2B−1 + c5h) 1 2T 1 2h 1 2 .
where (c1, c2, c3, c4, c5, c6, β) are some positive constants such that 1β > c2 and mF β > 3HFLK − 2LF .
Theorem 16 (Decreasing Stepsize) Under Assumption 1, if we set hk = h0/(k + 1) and Bk = B0 + [log(k + 1)] 100/99, we can bound the W1(µT , ν∞) as
W1(µT , ν∞) ≤ c1√
M(β−1 − c2) + c3 exp
{ −2 ( β−1mF − 3HFLK − 2LF ) (
T−1∑ k=0 hk)
}
+ c3β −3Md 3 2 (c6h 3 0 + c7β 3h0 B0 + c8h 2 0β 2) 1 2 .
where (c1, c2, c3, c6, c7, c8, β) are some positive constants such that 1β > c2 and mF β > 3HFLK − 2LF , as

heading: I PROOF OF THEOREM 7 AND COROLLARY 8
text: Proof [Proof of Theorem 7] d ( θ
(i) t − θ (j) t
) = −β−1 ( F (θ
(i) t )− F (θ (j) t ) ) dt
+ 1
M M∑ q [ ∇K(θ(i)t − θ (q) t )−∇K(θ (j) t − θ (q) t ) ] dt
− 1 M M∑ q ( F (θ (q) t )K(θ (i) t − θ (q) t )− F (θ (q) t )K(θ (j) t − θ (q) t )) ) dt
+
√ 2
β (dW(i)t − dW (j) t )
⇒d E M∑ ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2  = −E2 M∑ ij β−1 ( F (θ (i) t )− F (θ (j) t ) )( θ (i) t − θ (j) t ) dt
+ E M∑ ij 2 M M∑ q [ ∇K(θ(i)t − θ (q) t )−∇K(θ (j) t − θ (q) t ) ] ( θ (i) t − θ (j) t ) dt − E M∑ ij 2 M M∑ q ( F (θ (q) t )K(θ (i) t − θ (q) t )− F (θ (q) t )K(θ (j) t − θ (q) t )) )( θ (i) t − θ (j) t ) dt + E2 M∑ ij √ 2 β (dW(i)t − dW (j) t ) ( θ (i) t − θ (j) t ) ≤ −2β−1mFE M∑ ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2 dt− 2mKE M∑ ij
∥∥∥θ(i)t − θ(j)t ∥∥∥2 dt + 2HFLKE
M∑ ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2 dt + 2 √ 2
β E M∑ ij (dW(i)t − dW (j) t ) 2 1/2E M∑ ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2 1/2 .
Denote z(t)= E ∑M ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2. We have z(t)′ ≤ −(2β−1mF + 2mK − 2HFLK)z(t) + 4M √ d
β z(t) (33)
We can finish our proof by applying Gronwall Lemma on (33).
Proof [Proof for Corollary 8] d ( θ
(i) t − θ (j) t
) =
+ 1
M M∑ q [ ∇K(θ(i)t − θ (q) t )−∇K(θ (j) t − θ (q) t ) ] dt
− 1 M M∑ q ( F (θ (q) t )K(θ (i) t − θ (q) t )− F (θ (q) t )K(θ (j) t − θ (q) t )) ) dt
⇒d E M∑ ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2  =
+ E M∑ ij 2 M M∑ q [ ∇K(θ(i)t − θ (q) t )−∇K(θ (j) t − θ (q) t ) ] ( θ (i) t − θ (j) t ) dt − E M∑ ij 2 M M∑ q ( F (θ (q) t )K(θ (i) t − θ (q) t )− F (θ (q) t )K(θ (j) t − θ (q) t )) )( θ (i) t − θ (j) t ) dt ≤ −2mKE M∑ ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2 dt+ 2HFLKE M∑ ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2 dt
Denote z(t)= E ∑M ij ∥∥∥θ(i)t − θ(j)t ∥∥∥2. We have z(t)′ ≤ −(2mK − 2HFLK)z(t) (34)
We can finish our proof by applying Gronwall Lemma on (34).

heading: J PROOF OF THEOREM 9
text: Proof [Proof of Theorem 9] Our conclusion for B̃(µ̂T , µ̂∞) is essentially a specification of the result in Mattingly et al. (2002), which has also been applied in Xu et al. (2018).
Specifically, we rely on the following lemma, which is essentially Theorem 7.3 in Mattingly et al. (2002) and Lemma C.3 in Xu et al. (2018), considering the SDE (16):
dΘt = −FΘ(Θt)dt+ √
2β−1dW(Md)t As mentioned in Section 5, we firstly denote the distribution of Θt as ρΘt . Then we define the Θ̂k , [θ̂ (1) k , · · · , θ̂ (M) k ] ∈ RMd, which is actually the numerical solution of (16) using full gradient with Euler method. And we denote the distribution of Θ̂k as µ̂Θk .
Lemma 17 Let FΘ be Lipschitz-continuous with constant LΘ, and satisfy the dissipative property that 〈FΘ(Θ),Θ〉 ≥ mΘ ‖Θ‖2 − bΘ. Define VΘ(Θ) = C0 + LΘ/2‖Θ‖2. The Euler method for (16) has a unique invariant measure µ̂Θ∞, and for all test function fΘ such that |fΘ| ≤ VΘ(Θ), we have ∣∣∣E[fΘ(Θ̂k))]− EΘ̂∞∼µ̂Θ∞ [f(Θ̂∞)]∣∣∣ ≤ Cκρ−Md/2(1 + κemΘh) exp ( −2mΘkhρ Md
log(κ)
) ,
where ρ ∈ (0, 1), C > 0 are positive constants, and κ = 2LΘ(bΘβ +mΘβ +Md)/mΘ.
Now we adopt the fΘ : RMd → R defined as fΘ(Θ) = 1M ∑M i f(θ
(i)), where f : Rd → R is a Lf -Lipschitz function and satisfies our Assumption 2 and Θ , [θ(1), · · · ,θ(M)]. Similar to the proof of Lemma 12, we can find that fΘ : RMd is a Lf/ √ M -Lipschitz function. Furthermore, according to the Lemma 14, it is easily check that FΘ is LΘ-Lipschitz where LΘ = √ 2β−1LF + l ′.
Hence, when the β is small enough, we have Lf/ √ M ≤ √ 2β−1LF + l
′. Then we can set the C0 large enough to force fΘ to satisfy the condition in Lemma 17 that |fΘ| ≤ VΘ(Θ). Then, according to the exchangeability of the particle system {θ̂(i)k } and Lemma 14, we can bound the B̃(µ̂T , µ̂∞) as
B̃(µ̂T , µ̂∞) ≤ ∣∣∣E[fΘ(Θ̂T ))]− EΘ̂∞∼µ̂Θ∞ [f(Θ̂∞)]∣∣∣
≤ C2ςσ−Md/2(1 + ςemΘh) exp ( −2mΘThσMd/ log(ς) )
where ς = 2LΘ(Mbβ + mΘβ + Md)/mΘ, LΘ = √ 2β−1LF + l ′, mΘ = β−1m − m′, and (σ,C2, l ′,m′) are some positive constants independent of (T, M, h) and σ ∈ (0, 1).
To prove the bound for B̃(µ̂∞, ρ∞) Since Θ̂k = (θ̂(1)k , · · · , θ̂ (M) k ) can be considered as a solution to the SDE (16), standard results from linear FP equation can be applied. Specifically, for the B̃(µ̂∞, ρ∞) term, we rely on the following lemma adapted from Lemma C.4 in Xu et al. (2018)Chen et al. (2015), and is essentially the work of Chen et al. (2015) when taking K →∞.
Lemma 18 Under the same assumption as in Lemma 17, for the Lipschitz-continuous function fΘ(Θ) = 1 M ∑M i f(θ
(i)) mentioned above, the following bound is satisfied for some positive constant C: ∣∣∣∣∣ 1K K−1∑ k=1 E[fΘ(Θ̂k)]− EΘ∞∼ρΘ∞ [f(Θ∞)] ∣∣∣∣∣ ≤ C(hβ + βKh ) . The uniqueness of invariant measure of the Euler method from Lemma 17 implies the numerical solution Θ̂k to be ergodic. Then similar to the proof of Lemma 4.2 in Xu et al. (2018), we consider the case where K →∞. Take average over the {Θ̂k}K=1k=0 , we have
EΘ̂∞∼µ̂Θ∞ [fΘ(Θ̂∞)] = limK→∞ 1 K K−1∑ k=1 E[fΘ(Θ̂k)]
Now according to the exchangeability of the particle system {θ̂(i)k } and {θ (i) t }, we can bound the B̃(µ̂∞, ρ∞) as :
B̃(µ̂∞, ρ∞) ≤ ∣∣∣EΘ̂∞∼µ̂Θ∞ [fΘ(Θ̂∞)]− EΘ∞∼ρΘ∞ [fΘ(Θ∞)]∣∣∣ ≤ C3h/β
where C3 are some positive constant.

heading: K PROOF OF THEOREM 10
text: Proof [Proof of Theorem 10] Adopting the same notation used in the proof of the Theorem 5, we define Θk , [θ (1) k , · · · , θ (M) k ] and G Θ Ik , N Bk ∑ q∈Ik F(q)Θ(Θk). We denote the distribution of Θk as µΘk .
Θk+1 = Θk − β−1GΘIkhk + √ 2β−1hkΞk ,
We firstly give a bound to the W2(µΘk , µ̂ Θ k ) (the definition of the µ̂ Θ k has been mention in the last section). According to the proof of the Lemma 4.4 in Xu et al. (2018)
W2(µ Θ k , µ̂ Θ k ) ≤ kh(LΘΓ′ +MC4) ((6 + 2Γ′)β/B) 1/2
where Γ′ = 2(1 + 1/mΘ)(Mb+ 2M2C24 +Md/β) and C4 is some positive constant independent of (T, M, h). Note the fact that W1(µΘk , µ̂ Θ k ) ≤ W2(µΘk , µ̂Θk ) and W1(µk, µ̂k) ≤ 1√MW1(µ Θ k , µ̂ Θ k ) (see the proof of Lemma 13, similar result holds here), We get
W1(µT , µ̂T ) ≤ Th(LΘΓ′ +MC4) ((6 + 2Γ′)β/(BM)) 1/2 .
Let us compare the definition ofW1(µ, ν) and B̃(µ, ν)
W1(µ, ν) , sup ‖g‖lip≤1 |Eθ∼µ[g(θ)]− Eθ∼ν [g(θ)]|
B̃(µ, ν) , |Eθ∼µ[f(θ)]− Eθ∼ν [f(θ)]|
and we can derive the result that B̃(µT , µ̂K) ≤ LfW1(µT , µ̂T ). Now we finish our proof.

heading: L DISCUSSION ON THE COMPLEXITY OF OUR METHOD
text: The complexity of an algorithm mainly refers to its time complexity (corresponding to the number of iterations in our method i.e. T) and space complexity (corresponding to the number of particles used in our method i.e. M). Hence the complexity of our method can be well explored with our work, since our non-asymptotic convergence theory was developed w.r.t. the number of particles i.e. M and iterations i.e. T. Their relationship (tradeoff) was even discussed further in Experiment 6.1. Moreover, comparing (9) with (3) , one can easily find that our space complexity is exactly the same as SVGD and our computational time in each iteration is almost the same as SVGD with an extra addition operation. However, it worth noting that our method have much better performance in practice and no "pitfall" verified by both our theory and our experiments. .

heading: M COMPARISON WITH RELATED WORK
text: Firstly, our proposed framework SPOS is different from the recently proposed particle-optimization sampling framework (Chen et al., 2018), in the sense that we solve the nonlinear PDE (6) stochastically. For example they deterministically solve the equation in (6) ∂νt = β−1∇θ ·∇θνt approximately using blob method adopted from (Carrillo et al., 2017). Secondly, our method is also distinguishable to existing work on granular media equations such as (Durmus et al., 2018). Their work about the Granular media equations focuses on the following PDE:
∂tνt = ∇θ · ( νtβ −1F (θ) + νt (∇K ∗ νt(θ)) + β−1∇θνt ) , (35)
whereas our framework focuses on the following one: ∂tνt = ∇θ · ( νtβ −1F (θ) + νt (EY∼νtK(θ − Y )F (Y )−∇K ∗ νt(θ)) + β−1∇θνt ) . (36)
The extra term νt (EY∼νtK(θ − Y )F (Y )) in our framework makes the analysis much more challenging. The main differences are summarized below:
• Formulations are different. The extra term EY∼µtK(θ−Y )F (Y ) cannot be combined with the F (θ) term in their (35). This is because function F (θ) itself is a function independent of t; while EY∼µtK(θ− Y )F (Y ) depends on both θ and t. This makes our problem much more difficult.
• Assumptions are different. For example, the analysis on granular media equations in (Cattiaux et al., 2008) requires that F satisfies a special condition C(A, α), which is a strong condition impractical to be satisfied in our case; And Durmus et al. (2018) adopts different assumptions from ours with a different goal.
• For the Euler integrator, Durmus et al. (2018) does not consider an Euler solution. Furthermore, our sampling method needs "stochastic gradient" i.e. G(i)k , N Bk ∑ q∈Ik Fq(θ (i) k ) in
(9) for computational feasibility, which is quite different from the former work on particleSDE such as (Malrieu, 2003; Cattiaux et al., 2008). Few of the former work on particle-SDE considered the stochastic gradient issue.
To sum up, the main purpose of our paper is to provide a non-asymptotic analysis of our method instead of improving the former work on a certain type of PDE. This is also the reason why we said that part of our proof techniques are based on those for analyzing granular media equations.

heading: N EXTRA EXPERIMENTS
text: N.1 TOY EXPERIMENTS
We compare the proposed SPOS with other popular methods such as SVGD and standard SGLD on four mutil-mode toy examples. We aim to sample from four unnormalized 2D densities p(z)/exp{U(z)}, with the functional form provided in Rezende & Mohamed (2015). We optimize/sample 50 and 2000 particles to approximate the target distributions. The results are illustrated in Figure 4 and Figure 5, respectively.
N.1.1 BAYESIAN NEURAL NETWORKS FOR MNIST CLASSIFICATION
We perform the classification tasks on the standard MNIST dataset. A two-layer MLP 784-X-X-10 with ReLU activation function is used, with X being the number of hidden units for each layer. The training epoch is set to 100. The test errors are reported in Table 2. Surprisingly, the proposed SPOS outperforms other algorithms such as SVGD at a significant level, though it is just a simple modification of SVGD by adding in random Gaussian noise. This is partly due to the fact that our SPOS algorithm can jump out of local modes efficiently, as explained in Section 4.4.

Review summaries:
This paper considers the problem of Bayesian inference using particle optimization sampler. Similarly to SGLD, authors propose Stochastic Particle Optimization Sampler (SPOS), augmenting Stein Variational Gradient Descent (SVGD) with diminishing Gaussian noise, replacing the hard-to-compute term of the Chen et al. (2018) formulation. Various theoretical results are given.