Title: TOWARDS MORE THEORETICALLY-GROUNDED PARTI-CLE OPTIMIZATION SAMPLING FOR DEEP LEARNING

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 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.

1. INTRODUCTION:

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 SG-MCMC 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.

2. PRELIMINARIES:

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.

2.1. STOCHASTIC GRADIENT MCMC:

In Bayesian sampling, we aim to generate random samples from a posterior distribution p(θ|X ) ∝ p(X |θ)p(θ), where θ ∈ R d represents the model parameter with a prior distribution p(θ), and X {x q } N q=1 represents the observed data with likelihood p(X |θ) = q p(x q |θ). Define the potential energy as: U (θ) -log p(X |θ) -log p(θ) = -N q=1 log p(x q |θ) + 1 N log p(θ) N q=1 U q (θ). 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:

EQUATION

where

F (θ) ∇ θ U (θ) = N q=1 ∇ θ U q (θ) N q=1

F q (θ); t is the time index; β > 0 is the temperature parameter; and W t ∈ R d 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) :

EQUATION

)

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: G

k N B k q∈I k F q (θ k ), where I k is a random subset of [1, 2, • • • , N ]

with size B k in each iteration. Based on the above settings, SGLD uses the Euler method with stepsize h k to numerically solve (1) and obtains the update equation:

θ k+1 = θ k -β -1 G k h k + 2β -1 h k ξ k , ξ k ∼ N (0, I).

2.2. STEIN VARIATIONAL GRADIENT DESCENT:

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 {θ

EQUATION

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, {θ

k } M i=1 , are used to approximate the target distribution.

2.3. PARTICLE-OPTIMIZATION BASED SAMPLING METHODS:

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 SG-MCMC 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.,

EQUATION

) 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:

EQUATION

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:

EQUATION

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 } M i=1 that evolve over time t, i.e. ν t ≈ 1 M 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.

3. STOCHASTIC PARTICLE-OPTIMIZATION SAMPLING ALGORITHMS:

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β -1 dW t , is equivalent to solving the corresponding FP equation: ∂ν t = β -1 ∇ θ • ∇ θ ν t . Consequently, we decompose RHS of (4) into two parts:

F 1 ∇ θ • ν t β -1 F (θ t ) + (K * ν t )ν t and F 2 β -1 ∇ θ • ∇ θ ν t .

Our idea is to solve F 1 deterministically under a PDE setting, and solve F 2 stochastically based on its corresponding SDE. When adopting particle approximation for the density ν t , both solutions of F 1 and F 2 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 ≈ 1 M 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 = F 1 reduces to the following system of differential equations describing evolutions of the particles over time: ∀i

EQUATION

On the other hand, by solving ∂ t ν t = F 2 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 = -β -1 F (θ (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β -1 dW (i) t (8)

Figure 1 : Comparison of SPOS (left) and SVGD (right) on a multi-mode distribution. The circles with different colors are the resulting 100 particles, which are able to spread over all modes for SPOS.

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 } M i=1 .

In our theory, we show this intuition is actually 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 B k q∈I k F q (θ (i)

k ) evaluated with a minibatch of data of size B k for computational feasibility. Based on the Euler method (Chen et al., 2015) with a stepsize h k , (8) leads to the following updates for the particles at the k-th iteration

EQUATION

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. 

4. NON-ASYMPTOTIC CONVERGENCE ANALYSIS: THE CONVEX CASE:

(i) T } M i=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 W 1 , a special case of p-Wasserstein metric de-

fined as Wp(µ, ν) = inf ζ∈Γ(µ,ν) R d ×R d Xµ -Xν p dζ(Xµ, Xν ) 1/p

, where Γ(µ, ν) is the set of joint distributions on R d × R d 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) .

4.1. BASIC SETUP AND ASSUMPTIONS:

Due to the exchangeability of the particle system {θ 8), if we initialize all the particles θ (i) t

(i) t } M i=1 in (

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 {θ 9), and thus we denote the distribution of each θ (i) k as µ k . To this end, our analysis aims at bounding W 1 (µ 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.

(i) k } M i=1 in (

Assumption 1 Assume F and K satisfy the following conditions:

• There exist positive m F and m K , such that

F (θ) -F (θ ), θ -θ ≥ m F θ -θ 2 and ∇K(θ) -∇K(θ ), θ -θ ≤ -m K θ -θ 2 .

• F is bounded by H F and L F -Lipschitz continuous i.e., F (θ) ≤ H F and F (θ 1 ) -

F (θ 2 ) ≤ L F θ 1 -θ 2 ; . • K is L K -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 W 1 (µ T , ν ∞ ) in this section is to decompose it as follows:

W 1 (µ T , ν ∞ ) ≤ W 1 (µ T , ρ T -1 k=0 h k ) + W 1 (ρ T -1 k=0 h k , ν T -1 k=0 h k ) + W 1 (ν T -1 k=0 h k , ν ∞ ) . (10) 4.2 BOUNDS WITH STOCHASTIC PARTICLE APPROXIMATION We firstly bound W 1 (ρ T -1 k=0 h k , ν T -1 k=0 h k ) and W 1 (ν T -1 k=0 h k , ν ∞ )

with the following theorems. Theorem 3 Under Assumption 1 and let ρ 0 = ν 0 , there exist some positive constants c 1 and c 2 independent of (M, t) and satisfying c 2 < β -1 such that

EQUATION

Remark 1 According to Theorem 3, we can bound the

W 1 (ρ T -1 k=0 h k , ν T -1 k=0 h k ) term as W 1 (ρ T -1 k=0 h k , ν T -1 k=0 h k ) ≤ c1 √ M (β -1 -c2) . Furthermore, by letting t → ∞, we have W 1 (ρ ∞ , ν ∞ ) ≤ c1 √ M (β -1 -c2)

, which is an important intermediate result to prove the following theorem. Theorem 4 Under Assumption 1, the following holds:

EQUATION

)

To ensure W 1 (ν T -1 k=0 h k , ν ∞ ) 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.

4.3. BOUNDS WITH A NUMERICAL SOLUTION:

To bound the W 1 (µ T , ρ T -1 k=0 h k ) 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 h k = h that is small enough, the correspond- ing W 1 (µ T , ρ T h ) is bounded as W 1 (µ T , ρ T h ) ≤ c 4 M d 3 2 β -3 (c 5 β 2 B -1 + c 6 h) 1 2 T 1 2 h 1 2 (13)

where B is the fixed size of the minibatch in each iteration and (c 4 , c 5 , c 6 ) 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 h k = h 0 /(k + 1), and let the minibatch size in each iteration k be

EQUATION

)

where B 0 is the initial minibatch size, and (c 4 , c 7 , c 8 , c 9 ) are some positive constants independent of (M, T, h 0 ).

Note B k increases at a very low speed, e.g., only by 15 after 10 5 iterations, thus it would not affect algorithm efficiency. According to Theorem 6, W 1 (µ T , ρ T -1 k=0 h k ) would approach zero when h 1/2 0 M → 0. By directly combining results from Theorem 3-6, one can easily bound the target W 1 (µ T , ν ∞ ). Detailed statements are given in Theorem 15-16 in Section H of the SM.

4.4. A PITFALL OF SVGD:

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 λ = m F β + m K -H F L K > 0. Then the EPD of SPOS is bounded as: EPD M i,j E θ (i) t -θ (j) t 2 ≤ C 1 e -2λt + 4 d β M λ , where C 1 = M (M -1)Γ -4 dβ -1 M λ .

Remark 2 There are two interesting cases: i) When C 1 > 0, the EPD would decrease to the bound 4 dβ -1 M/λ along time t. This represents the phenomenon of an attraction force between particles;

ii) When C 1 < 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 ≤ C 0 e -2λt

, where

C 0 = M i,j θ (i) 0 -θ (j) 0 2 and λ = m K -H F L K .

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.

5. NON-ASYMPTOTIC CONVERGENCE ANALYSIS: THE NONCONVEX CASE:

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 L f -continuous function f satisfying Assump- tion 2 below.

Note such metric has also been adopted in (Vollmer et al., 2016; Chen et al., 2015) . Our analysis considers (T, M, h k ) as variables in B. In addition, we use { θ(i) k } M i=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 ] ∈ R M d .

Consequently, the nonlinear SDE system (8) can be turned into a linear SDE ,which means Θ t is driven by the following linear SDE:

EQUATION

where

F Θ (Θ t ) [β -1 F (θ (1) t ) -1 M M j=1 ∇K(θ (1) t -θ (j) t ) + 1 M M j=1 K(θ (1) t - θ (j) t )F (θ (j) t ), • • • , β -1 F (θ (M ) t )-1 M M j=1 ∇K(θ (M ) t -θ (j) t )+ 1 M M j=1 K(θ (M ) t -θ (j) t )F (θ (j) t )] is a vector function R M d → R M d , and W (M d) t is Brownian motion of dimension M d. Similarly, we can define Θk [ θ(1) k , • • • , θ(M) k ] ∈ R M d

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, D k ψ, are bounded by a function V, i.e.,

D k ψ ≤ H k V p k for k = (0, 1, 2, 3, 4), H k , p k > 0. Furthermore, the expectation of V on {Θ t } is bounded: sup l EV p (Θ t ) < ∞, and V is smooth such that sup s∈(0,1) V p (sΘ + (1 -s) Θ ) ≤ H (V p (Θ) + V p (Θ )), ∀Θ, Θ , p ≤ max{2p k } for H > 0. Assumption 3 i) F , K and ∇K are L F , L K 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 h k = h, we can have the following results:

EQUATION

where

ς = 2L Θ (M bβ + m Θ β + M d)/m Θ , L Θ = √ 2β -1 L F + l , m Θ = β -1 m -m , and (σ, C 2 , C 3 , 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 B k = B and h k = h, B(µ T , μT ) is bounded as

B(µ T , μT ) ≤ C 5 T h(L Θ Γ + M C 4 ) ((6 + 2Γ )β/(BM )) 1/2 .

where

Γ = 2(1 + 1/m Θ )(M b + 2M 2 C 2 4 + M d/β

) and , (C 4 , C 5 ) 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

(C 2 , C 3 , C 4 , C 5 , C 6 ) such that: B(µ T , ν ∞ ) ≤ C 2 ςσ -M d/2 (1 + ςe m Θ h ) exp -2m Θ T hσ M d / log(ς) + C 3 h/β + C 5 T h(L Θ Γ + M C 4 ) ((6 + 2Γ )β/(BM )) 1/2 + C 6 M -1/2 ,

where σ, ς and Γ are the same as those in Theorem 9-10.

6. EXPERIMENTS:

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. We follow Chen et al. (2015) and consider a standard Gaussian model where x i ∼ N (θ, 1), θ ∼ N (0, 1). 1000 data samples {x i } 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.

6.1. BOUNDS ILLUSTRATION WITH A SIMPLE GAUSSIAN EXAMPLE:

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. 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(-θθ 2 2 /η 2 ), with the bandwidth set to η = med 2 / log M . 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.

6.2. BAYESIAN NEURAL NETWORKS FOR REGRESSION:



6.3. BAYESIAN EXPLORATION IN DEEP REINFORCEMENT LEARNING:

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 θ: 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.

EQUATION

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.

7. CONCLUSION:

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.

A DENSITY FUNCTION OF THE MULTI-MODE DISTRIBUTION IN SECTION 3:

The negative log-density function of the multi-mode distribution in Section 3 is defined as:

U (θ) e 3 4 θ 2 -3 2 10 i=1 ci sin( 1 4 π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, c i is the i-th element of c.

B GRONWALL LEMMA:

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)

C PROOF OF THEOREM 2:

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 Ω ⊂ R d , 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

EQUATION

for all differentiable function u(•).

In Theorem 2, we want to prove a particle representation of the following PDE:

∂ t µ t = F 1 = -∇ θ • (µ t F (θ t ) + (K * µ t )µ t ) -∂ θ (Gf )(θ, t) ⇒ ∂ t µ t u(θ)dθ = -∂ θ (Gf )(θ, t)u(θ)dθ ,

where f (θ, t) = µ t . Consequently, we have

EQUATION

By applying ( 20) in ( 21), we have

∂ t f (θ, t)u(θ)dθ = -∂ θ (Gf )(θ, t)u(θ)dθ = Gf (θ, t)∂ θ u(θ)dθ .

Since f dθ µ(dθ, t), we have

EQUATION

In particle approximation, we have

µ(t) = 1 M 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.

D PROOF OF THEOREM 3:

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 :

EQUATION

where L( θt ) denotes the probability law of θt , Wt ∈ R d 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 } M i=1 ,

d θ(i) t = -β -1 F ( θ(i) t ) -E Yi∼νt K( θ(i) t -Y i )F (Y i ) + ∇K * ν t ( θ(i) t ) + 2β -1 d W(i) t L( θ(i) t ) = ν t dθ (24)

where Y i 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

EQUATION

where

A ij (t) = -β -1 F (θ (i) t ) -F ( θ(i) t ) • θ (i) t - θ(i) t B ij (t) = ∇K(θ (i) t -θ (j) t ) -∇K( θ(i) t - θ(j) t ) • θ (i) t - θ(i) t C ij (t) = ∇K( θ(i) t - θ(j) t ) -∇K * ν t ( θ(i) t ) • θ (i) t - θ(i) t F ij (t) = -F (θ (j) t )K(θ (i) t -θ (j) t ) -F ( θ(j) t )K(θ (i) t -θ (j) t ) • θ (i) t - θ(i) t G ij (t) = -F ( θ(j) t )K(θ (i) t -θ (j) t ) -F ( θ(j) t )K( θ(i) t - θ(j) t ) • θ (i) t - θ(i) t H ij (t) = -F ( θ(j) t )K( θ(i) t - θ(j) t ) -E Yj ∼νt F (Y j )K( θ(i) t -Y j ) • θ (i) t - θ(i) t

For the A ij (t) term, according to the i) in Assumption 1 for F , we have

ij A ij (t) = - ij β -1 F (θ (i) t ) -F ( θ(i) t ) • θ (i) t - θ(i) t ≤ -β -1 m F M i θ (i) t - θ(i) t 2

For the B ij (t) term, applying the concave condition for K and the oddness of ∇K in Assumption 1, we have

ij B ij (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 m K M ij θ (i) t - θ(i) t -(θ (j) t - θ(j) t ) 2 ≤ 0

For the C ij (t) term, we have

E j C ij (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 ) 2   1/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 H ij (t) term,

E j H ij (t) ≤ E θ (i) t - θ(i) t 2 1/2   E j F ( θ(j) t )K( θ(i) t - θ(j) t ) -E Yj ∼νt F (Y j )K( θ(i) t -Y j ) 2    1/2 = E θ (i) t - θ(i) t 2 1/2   j E F ( θ(j) t )K( θ(i) t - θ(j) t ) -E Yj ∼νt F (Y j )K( θ(i) t -Y j ) 2   1/2 ≤H F √ 2M E θ (i) t - θ(i) t 2 1/2

For the F ij (t) and G ij (t) terms, we have:

ij F ij (t) = - ij F (θ (j) t )K(θ (i) t -θ (j) t ) -F ( θt (j) )K(θ (i) t -θ (j) t ) • θ (i) t - θ(i) t ≤ ij L F θ (j) t - θ(j) t θ (i) t - θ(i) t ≤ 2L F M i θ (i) t - θ(i) t 2 . ij G ij (t) = ij F ( θ(j) t )K(θ (i) t -θ (j) t ) -F ( θ(j) t )K( θ(i) t - θ(j) t ) • θ (i) t - θ(i) t ≤ H F L K ij θ (i) t - θ(i) t -(θ (j) t - θ(j) t ) θ (i) t - θ(i) t ≤ 3H F L K M 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 + H F √ 2 √ M γ(t) .

where

λ 1 = β -1 m F -3H F L K -2L F ⇒ ( γ(t) - (H ∇K + H F )/ √ 2 √ M (β -1 -3H F L K -2L F ) ) ≤ -λ 1 ( γ(t) - (H ∇K + H F )/ √ 2 √ M (β -1 -3H F L K -2L F ) ) 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 + H F )/ √ 2 √ M (β -1 -3H F L K -2L F )

Hence, there exist some positive constant (c 1 , c 2 ) such that:

EQUATION

where (1) holds due to the relationship between W 1 and W 2 metric Givens & Shortt (1984) , (2) due to the definition of W 2 and (3) due to the result from the previous proof.

E PROOF OF THEOREM 4:

Proof [Proof of Theorem 4] Firstly, what we aim at is W 1 (ν t , ν ∞ ) ≤ c 3 exp (-2λ 1 t) in this theorem.

According to the relationship between W 1 and W 2 metric Givens & Shortt (1984) , once we bound W 2 (ν t , ν ∞ ) as W 2 (ν t , ν ∞ ) ≤ c 3 exp (-2λ 1 t), 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 } M i=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 } M i=1 .

We denote the distribution of each θ

(i)

t,2 at t as ρ t,2 .

Since we need to bound W 2 (ν t , ν ∞ ), we make the following decomposition:

EQUATION

Note that ρ 0,1 = ν 0 and ρ 0,2 = ν ∞ . Then, according to (26), we have

W 2 (ν t , ρ t,1 ) ≤ c 1 √ M (β -1 -c 2 ) and W 2 (ρ t,2 , ν ∞ ) ≤ c 1 √ M (β -1 -c 2 )

Now we need to focus on the term

W 2 (ρ t,1 , ρ t,2 ). Since W 2 (ρ 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 (ξ 1 ij (t) + ξ 2 ij (t) + ξ 3 ij (t) + ξ 4 ij (t))dt where ξ 1 ij (t) = -β -1 F (θ (i) t,1 ) -F (θ (i) t,2 ) • θ (i) t,1 -θ (i) t,2 ξ 2 ij (t) = ∇K(θ (i) t,1 -θ (j) t,1 ) -∇K(θ (i) t,2 -θ (j) t,2 ) • θ (i) t,1 -θ (i) t,2 ξ 3 ij (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 ξ 4 ij (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 ξ 1 ij (t) terms, according to the i) in Assumption 1 for F , we have

ij ξ 1 ij (t) = - ij β -1 F (θ (i) t,1 ) -F (θ (i) t,2 ) • θ (i) t,1 -θ (i) t,2 ≤ -β -1 m F M i θ (i) t,1 -θ (i) t,2 2 .

For the ξ 2 ij (t) term, applying the concave condition for K and the oddness of ∇K in Assumption 1, we have

ij ξ 2 ij (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 m K M ij θ (i) t,1 -θ (i) t,2 -(θ (j) t,1 -θ (j) t,2 ) 2 ≤ 0 .

For the ξ 3 ij (t) terms, after applying the L F -Lipschitz property for F and K ≤ 1, we have

ij ξ 3 ij (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 L F θ (j) t,1 -θ (j) t,2 θ (i) t,1 -θ (i) t,2 ≤ 2L F M i θ (i) t,1 -θ (i) t,2 2 .

For the ξ 4 ij (t) terms :

EQUATION

t,2 ) • θ (i) t,1 -θ (i) t,2 ≤ H F L K ij θ (i) t,1 -θ (i) t,2 -(θ (j) t,1 -θ (j) t,2 ) θ (i) t,1 -θ (i) t,2 ≤ 3H F L K M i θ (i) t,1 -θ (i) t,2 2 . Now we have r (t) ≤ -2(β -1 m F -3H F L K -2L F )r(t) .

According to the Gronwall lemma,

r(t) ≤ r(0)e -2λ1t ,

where

λ 1 = β -1 m F -3H F L W -2L F .

Consequently, there exists some positive constant c 3 such that

W 2 (ρ t,1 , ρ t,2 ) ≤ c 3 e -2λ1t

Then we have

W 2 (ν t , ν ∞ ) ≤ c 3 e -2λ1t + c 1 √ M (β -1 -c 2 ) + c 1 √ M (β -1 -c 2 )

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 W 2 (ν t , ν ∞ ) ≤ c 3 e -2λ1t . Now we finish our proof.

F PROOF OF THEOREM 5:

To bound the W 1 (µ T , ρ T -1 k=0 h k ) term, note the original SDE driving the particles {θ 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 [θ

(i) t } in (

(1) t , • • • , θ (M ) t ] ∈ R M d .

Consequently, Θ t is driven by the following linear SDE:

EQUATION

where

F Θ (Θ t ) [β -1 F (θ (1) t ) -1 M M j=1 ∇K(θ (1) t -θ (j) t ) + 1 M M j=1 K(θ (1) t - θ (j) t )F (θ (j) t ), • • • , β -1 F (θ (M ) t )-1 M M j=1 ∇K(θ (M ) t -θ (j) t )+ 1 M M j=1 K(θ (M ) t -θ (j) t )F (θ (j) t )] is a vector function R M d → R M d , and W (M d) t

is Brownian motion of dimension M d.

Now we define the:

F (q)Θ (Θ t ) [β -1 F q (θ (1) t )-1 M N M j=1 ∇K(θ (1) t -θ (j) t )+ 1 M M j=1 K(θ (1) t - θ (j) t )F q (θ (j) t ), • • • , β -1 F q (θ (M ) t ) -1 M N M j=1 ∇K(θ (M ) t -θ (j) t ) + 1 M M j=1 K(θ (M ) t - θ (j) t )F q (θ (j) t )]. We can verify that F Θ (Θ t ) = N q=1 F (q)Θ (Θ t ). Then we define Θ k [θ (1) k , • • • , θ (M ) k ] and G Θ I k N B k q∈I k F (q)Θ (Θ k ).

We can verify that the following result holds:

EQUATION

where Ξ k ∼ N (0, I M d×M d ). 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 W 1 (µ k , ρ t ) ≤ 1 √ M W 1 (µ Θ k , ρ Θ t )

Proof [Proof of Lemma 13] Let us recall the definition of W 1 metric and its Kantorovich-Rubinstein duality Villani (2008) , i.e. W 1 (µ, ν)

sup g lip ≤1 |E θ∼µ [g(θ)] -E θ∼ν [g(θ)]|. We can prove the fact that if g(θ) : R d → R is a L g -Lipschitz function in R d , the g Θ (Θ), defined as g Θ (Θ) = 1 √ M M i g(θ (i) ), is a L g -Lipschitz function in R M d , where Θ [θ (1) , • • • , θ (M ) ]: g Θ (Θ 1 ) -g Θ (Θ 2 ) ≤ 1 √ M M i=1 g(θ (i) 1 ) -g(θ (i) 2 ) ≤ L g √ M M i=1 θ (i) 1 -θ (i) 2 ≤ L g √ M √ M M i=1 θ (i) 1 -θ (i) 2 2 = L g Θ 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 W 1 metric, we derive that

W 1 (µ 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 W 1 (µ Θ k , ρ Θ t )

Lemma 14 Assuming F (0) = 0. If F in ( 9) is Lipschitz with constant L F , and satisfies the dissipative property that

F (θ), θ ≥ m F θ 2 -b. Then F Θ in (28) is Lipschitz-continuous with constant √ 2β -1 L F + l and satisfies F Θ (Θ), Θ ≥ (β -1 m F -m ) Θ 2 -β -1 M b,

where l and m are some positive constants.

Proof [Proof of Lemma 14]

F Θ (Θ 1 ) -F Θ (Θ 2 ) = M i ω 1 i + ω 2 i + ω 3 i 2 ≤ M i ( ω 1 i + ω 2 i + ω 3 i ) 2

where

ω 1 i = β -1 F (θ (i) 1 ) -β -1 F (θ (i) 2 ) ≤ β -1 L F θ (i) 1 -θ (i) 2 ω 2 i = 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 ))

≤ (2L K H F + L F ) θ (i) 1 -θ (i) 2 ω 3 i = - 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 ω 1 i + ω 2 i + ω 3 i 2 ≤ M i 2(β -1 L F + 2L K H F + L F + L ∇K ) 2 θ (i) 1 -θ (i) 2 2 + 2 M j θ (j) 1 -θ (j) 2 2 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 Θ (Θ), Θ ≥ (β -1 m -2L F -3L ∇K ) Θ 2 -β -1 M b

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 (θ), θ ≥ m F θ 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 D KL (µ Θ k ρ Θ kh ) satisfies:

D KL (µ Θ k ρ Θ kh ) ≤ (A 0 β a B + A 1 h)kh with A 0 = 2( L F β + l ) 2 a 2 + 2(1 ∨ 1 β -1 m F -m )(2a 2 1 + M d β ) + a 2 1 A 1 =6( L F β + l ) 2 (βA 0 + M d)

and a 1 , a 2 are some positive constants. When the β is small enough, there exist some positive constants a 3 , a 4 such that

A 0 ≤ a 3 M d β 3 , A 1 ≤ a 4 M d β 4

Similar to the proof of Lemma 14, it is easy to verify that there exists some positive constant a 5 such that

F Θ (Θ 1 ) -F Θ (Θ 2 ), Θ 1 -Θ 2 ≥ (β -1 m F -a 5 ) Θ 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 W

1 (µ Θ k , ρ Θ kh ) ≤ C D KL (µ Θ k ρ Θ kh )

. According to Corollary 4 and Lemma 8 in Bolley & Villani (2005) , we can derive an explicit expression for C :

C ≤ a 6 β -1 M d ,

when β is small enough and a 6 is some positive constant.

Applying the Lemma 13, we have

W 1 (µ k , ρ kh ) ≤ 1 √ M W 1 (µ Θ k , ρ Θ kh ) ≤ a 6 M d 3 2 β -3 (a 3 a β 2 B -1 + a 4 h) 1 2 k 1 2 h 1 2

Let k = T and we can finish the poof.

G PROOF OF THEOREM 6:

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:

EQUATION

where Ξ k ∼ N (0, I M d×M d ) and h k = h0 k+1 . We assume E(G Θ

I k ) = F Θ (Θ k ), ∀Θ ∈ R M d

, which is a general assumption due to the way that we choose the minibatch I k . We need to define q(t), which will be used in the following proof:

q(t) = {k ∈ R| k-1 i=0 h i ≤ t < k i=0 h i } .

Furthermore, we define -1 i=0 h i 0 and 0 i=0 h i h 0 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 h i )   ds + 2 β t 0 W (M d) s , where I(s) ≡ I k for t ∈ k-1 i=0 h i , k i=0 h i , 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 h i ) 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 G s (Λ(s)) ds + 2 β t 0 W (M d) s where G t (x) := E   GΘ I(t)   Θ( q(t)-1 i=0 h i )   |Θ(t) = x   .

Let P t Λ := 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 P t Λ and P t Θ :

D KL (P t Λ P t Θ ) = -d P t Λ log d P t Λ d P t Θ = β 4 t 0 E F Θ (Λ(s)) -G s (Λ(s)) 2 ds = β 4 t 0 E F Θ (Θ(s)) -G s (Θ(s)) 2 ds

The last line follows that L(Θ(s)) = L(Λ(s)), ∀s.

In the following proof, we will let t = k-1 i=0 h i for some k ∈ R. Now we can use the martingale property of Itô integral to derive:

EQUATION

For the first part (31), we consider some s ∈

[ j-1 i=0 h i , j i=0 h i ).

The following equation holds:

Θ(s) -Θ( j-1 i=0 h i ) = -(s - j-1 i=0 h i )G Θ Ij + 2/β(W (M d) s -W (M d) j-1 i=0 hi ) = -(s - j-1 i=0 h i )G Θ Ij + (s - j-1 i=0 h i )(F Θ (Θ j ) -G Θ Ij ) + 2/β(W (M d) s -W (M d) 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 h i ) 2 ≤3 h 0 2 (j + 1) 2 E G Θ Ij 2 + 3 h 0 2 (j + 1) 2 E F Θ (Θ j ) -G Θ Ij 2 + 6h 0 M d β(j + 1) ≤12 h 0 2 (j + 1) 2 max 0≤j≤k-1 (L 2 F Θ E Θ j 2 + b 1 ) + 6h 0 M d β(j + 1)

where b 1 is some positive constant.

Consequently, the first part above (31) can be bounded as:

βL 2 F Θ 2 k-1 j=0 j i=0 hi j-1 i=0 hi E Θ(s) -Θ( q(s)-1 i=0 h i ) 2 ds ≤ βL 2 F Θ 2 k-1 j=0 12 h 0 3 (j + 1) 3 max 0≤j≤K-1 (L 2 F Θ E Θ j 2 + b 1 ) + 6h 0 2 M d β(j + 1) 2 ≤π 2 βL 2 F Θ h 0 3 max 0≤j≤K-1 (L 2 F Θ E Θ j 2 + b 1 ) + π 2 L 2 F Θ h 0 2 M d 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 h i )) -GΘ I(s)   Θ( q(s)-1 i=0 h i )   2 ds = k-1 j=0 βh 0 2(j + 1) E F Θ (Θ j ) -G Θ Ij 2 ≤βh 0 max 0≤j≤k-1 (L 2 F Θ E Θ j 2 + b 1 )   4 B 0 + k-1 j=1 4 (j + 1)(B 0 + log 100 99 (j + 1))   ≤βh 0 max 0≤j≤k-1 (L 2 F Θ E Θ j 2 + b 1 )   4 B 0 + k-1 j=1 4 (j + 1)log 100 99 (j + 1)   ≤(b 2 + 4 B 0 )βh 0 max 0≤j≤k-1 (L 2 F Θ E Θ j 2 + b 1 ) ,

where the last line follows from the fact that when r > 1,

k-1 j=1 4 (j + 1) log r (j + 1) ≤ ∞ j=1 4 (j + 1) log r (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

D KL (µ Θ k ρ Θ k-1 i=0 hi ) ≤ D KL (P k-1 i=0 hi Λ P k-1 i=0 hi Θ ) ≤ π 2 βL 2 F Θ h 0 3 max 0≤j≤k-1 (L 2 F Θ E Θ j 2 + b 1 ) + π 2 L 2 F Θ h 0 2 M d 2 + (b 2 + 4 B 0 )βh 0 max 0≤j≤k-1 (L 2 F Θ E Θj 2 + b 1 ) ≤ (π 2 βL 2 F Θ h 0 3 + b 2 βh 0 + 4 B 0 βh 0 ) max 0≤j≤k-1 (L 2 F Θ E Θ j 2 + b 1 ) + π 2 L 2 F Θ h 0 2 M d 2 .

Lemma 3.2 in Raginsky et al. (2017) has provided a uniform bound to max 0≤j≤k-1 (L 2 F Θ E Θ j 2 + b 1 ). Hence we can tell that D KL (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 L F Θ = √ 2β -1 L F + l , it is easy to verify that when β is small enough, there exists some positive constants b 3 , b 4 , b 5 and b 6 such that:

D KL (µ Θ k ρ Θ k-1 i=0 hi ) ≤(π 2 βL 2 F Θ h 0 3 + b 2 βh 0 + 4βh 0 B 0 ) max 0≤j≤K-1 (L 2 F Θ E Θj 2 + b 1 ) + π 2 L 2 F Θ h 0 2 M d 2 ≤(b 3 h 3 0 + b 4 β 3 h 0 B 0 + b 5 h 2 0 β 2 ) M d β 4 .

Similar to the proof of Theorem 5, we can bound the W 1 (µ Θ 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 a 6 such that:

W 1 (µ Θ k ρ Θ k-1 i=0 hi ) ≤ a 6 ( M d β ) D KL (µ Θ k ρ Θ k-1 i=0 hi ) ≤ a 6 β -3 M 3 2 d 3 2 (b 3 h 3 0 + b 4 β 3 h 0 B 0 + b 5 h 2 0 β 2 ) 1 2 .

According to Lemma 13, we have

W 1 (µ k , ρ kh ) ≤ 1 √ M W 1 (µ Θ k ρ Θ k-1 i=0 hi ) = a 6 β -3 M d 3 2 (b 3 h 3 0 + b 4 β 3 h 0 B 0 + b 5 h 2 0 β 2 ) 1 2

Let k = T and we can finish the proof.

H DETAILED STATEMENTS OF NONASYMPTOTIC CONVERGENCE UNDER THE CONVEX CASE:

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 h k = h 0 and B k = B 0 , we can bound the W 1 (µ T , ν ∞ ) as

W 1 (µ T , ν ∞ ) ≤ c 1 √ M (β -1 -c 2 ) + c 3 exp -2 β -1 m F -3H F L K -2L F T h + c 3 M d 3 2 β -3 (c 4 β 2 B -1 + c 5 h) 1 2 T 1 2 h 1 2 .

where (c 1 , c 2 , c 3 , c 4 , c 5 , c 6 , β) are some positive constants such that 1 β > c 2 and m F β > 3H F L K -2L F .

Theorem 16 (Decreasing Stepsize) Under Assumption 1, if we set h k = h 0 /(k + 1) and B k = B 0 + [log(k + 1)] 100/99 , we can bound the W 1 (µ T , ν ∞ ) as

W 1 (µ T , ν ∞ ) ≤ c 1 √ M (β -1 -c 2 ) + c 3 exp -2 β -1 m F -3H F L K -2L F ( T -1 k=0 h k ) + c 3 β -3 M d 3 2 (c 6 h 3 0 + c 7 β 3 h 0 B 0 + c 8 h 2 0 β 2 ) 1 2 .

where (c 1 , c 2 , c 3 , c 6 , c 7 , c 8 , β) are some positive constants such that 1 β > c 2 and m F β > 3H F L K -2L F , as As mentioned in Section 5, we firstly denote the distribution of Θ t as ρ Θ t . Then we define the

Θk [ θ(1) k , • • • , θ(M) k ] ∈ R M d

, 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 Θ (Θ) = C 0 + 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 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. 

EQUATION

:



Figures and Tables:

Figure 2: The estimation errors with different particles.

Figure 3: Policy learning with Bayesian exploration in policy-gradient methods on six scenarios with SVPG and SPOS-PG.

κ) ,where ρ ∈ (0, 1), C > 0 are positive constants, andκ = 2L Θ (b Θ β + m Θ β + M d)/m Θ . Now we adopt the f Θ : R M d → R defined as f Θ (Θ) = 1 M M i f (θ (i) ), where f : R d → R is a L f -Lipschitz function and satisfies our Assumption 2 and Θ [θ (1) , • • • , θ (M ) ]. Similar to the proof of Lemma 12, we can find that f Θ : R M d is a L f / √ M -Lipschitz function. Furthermore, according to the Lemma 14, it is easily check that F Θ is L Θ -Lipschitz where L Θ = √ 2β -1 L F + l . Hence, when the β is small enough, we have L f / √ M ≤ √ 2β -1 L F + l .Then we can set the C 0 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 , μ∞ ) asB(μ T , μ∞ ) ≤ E[f Θ ( ΘT ))] -E Θ∞∼μ Θ ∞ [f ( Θ∞ )] ≤ C 2 ςσ -M d/2 (1 + ςe m Θ h ) exp -2m Θ T hσ M d / log(ς)

Figure 4: Illustration of different algorithms on toy distributions. Dots are the final particles; the blue regions represent ground true densities. Each column is a distribution case. First row: standard SGLD; Second row: SVGD; Third row: SPOS.

Figure 5: Illustration of different algorithms on toy distributions. Dots are the final particles; the blue regions represent densities estimated by the particles. Each column is a distribution case. First row: ground true densities; Second row: standard SGLD; Third row: SVGD; Fourth row: SPOS. 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 Table2. 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.

Averaged predictions with standard deviations in terms of RMSE and log-likelihood on test sets.

Table content: <table><tr><td>We conduct experiments for</td><td/><td/><td/><td/></tr><tr><td>Bayesian learning of DNNs,</td><td/><td/><td/><td/></tr><tr><td>where we Bayesian DNNs are used to model weight uncer-tainty of neural networks, an</td><td colspan="5">Test RMSE SPOS Boston_Housing 2.961 ± 0.109 2.829 ± 0.126 -2.591 ± 0.029 -2.532 ± 0.082 Test Log likelihood Dataset SVGD SVGD SPOS</td></tr><tr><td>important topic that has been well explored (Hernández-</td><td>Concrete Energy Kin8nm (0.4)</td><td colspan="4">6.157 ± 0.082 5.071 ± 0.1495 -3.247 ± 0.01 -3.062 ± 0.037 1.291 ± 0.029 0.752 ± 0.0285 -1.534 ± 0.026 -1.158 ± 0.073 0.075 ± 0.001 0.079 ± 0.001 1.138 ± 0.004 1.092 ± 0.013</td></tr><tr><td>Lobato &amp; Adams, 2015; Blun-dell et al., 2015; Li et al.,</td><td>Naval (0.4) CCPP Winequality</td><td colspan="4">0.004 ± 0.000 0.004 ± 0.000 4.127 ± 0.027 3.939 ± 0.0495 -2.843 ± 0.006 -2.794 ± 0.025 4.032 ± 0.008 4.145 ± 0.02 0.604 ± 0.007 0.598 ± 0.014 -0.926 ± 0.009 -0.911 ± 0.041</td></tr><tr><td>2016; Louizos &amp; Welling,</td><td>Yacht (0.4) Protein</td><td colspan="4">1.597 ± 0.099 0.84 ± 0.0865 4.392 ± 0.015 4.254 ± 0.005 -2.905 ± 0.010 -2.876 ± 0.009 -1.818 ± 0.06 -1.446 ± 0.121</td></tr><tr><td/><td>YearPredict</td><td>8.684 ± NA</td><td>8.681 ± NA</td><td>-3.580 ± NA</td><td>-3.576 ± NA</td></tr></table>

Classification error of FNN on MNIST.

Table content: <table><tr><td>Method</td><td colspan="2">Test Error 400-400 800-800</td></tr><tr><td>SPOS</td><td>1.32%</td><td>1.24%</td></tr><tr><td>SVGD</td><td>1.56%</td><td>1.47%</td></tr><tr><td>SGLD</td><td>1.64%</td><td>1.41%</td></tr><tr><td>RMSprop</td><td>1.59%</td><td>1.43%</td></tr><tr><td>RMSspectral</td><td>1.65%</td><td>1.56%</td></tr><tr><td>SGD</td><td>1.72%</td><td>1.47%</td></tr><tr><td colspan="2">BPB, Gaussian 1.82%</td><td>1.99%</td></tr><tr><td>SGD, dropout</td><td>1.51%</td><td>1.33%</td></tr></table>

Peer review (iclr2019):
 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. This paper was a pleasant read until I decided to check the proof of Theorem 3. I was not able to understand transitions in some of the steps and certain statements in the proof seem wrong. Theorem 3: "Note that and are initialized with the same initial distribution µ0 = ν0 and we can also set to be independent of , we can have . ." - this doesn*t seem right to me. Expectation of squared difference of two independent and identically distributed random variables is not 0, assuming expectation is with respect to their joint density. "Then according to the Gronwall Lemma, we have" - I don*t see how the resulting inequality was obtained. When I tried applying Gronwall Lemma, it seems that authors forgot to multiply by and . Could you please elaborate how exactly Gronwall Lemma was used in this case. "... some positive constants c1 and c2 independent of (M, d)$ - in the proof authors introduce additional assumption "We can tune the bandwidth of the RBF kernel to make ∇K ≤ H_∇K, which is omitted in the Assumption due to the space limit." First, there is a missing norm, since ∇K is a vector and H_∇K is I believe a scalar constant. Second, c1 = H_∇K + H_F, which both bound norm of d-dimensional vector and hence depend on d. I also suggest that all assumptions are included in the theorem statements, especially since authors have another assumption requiring large bandwidth. Additionally, feasibility of these both assumptions being satisfied should be explored (it seems to me that they can hold together, but it doesn*t mean that part of assumptions can be moved to the supplement). I find using Wasserstein-1 metric misleading in the theorem statement . This is not what authors really bound - from the proof it can be seen that they bound W_1 with W_2 and then with just an expectation of l2 norm. Moreover I don*t understand the meaning of this bound. Theorem is concerned with W_1 distance between two atomic measures. What is the expectation over? Note that atom locations are supposed to be fixed for the W_1 to make sense in this context (and the expectation is over the coupling of discrete measures defined by weights of the atoms, not atom locations). "Note the first bullet indicates U to be a convex function and W to be ... " I think it should be K, not W. Theorems 3-6 could be lemmas, while there should be a unifying theorem for the bound. Finally, I think notation should be changed - same letter is used for Wasserstein distance and Wiener process. Other comments: Example in Figure 1 is somewhat contrived - clearly gradient based particle sampler will never escape the mode since all modes are disconnected by regions with 0 density. Proposed method on the other hand will eventually jump out due to noise, but it doesn*t necessarily mean it produces better posterior estimate. Something more realistic like a mixture of Gaussians, with density bounded away from zero across domain space, will be more informative. It is not sufficient to report RMSE and test log likelihood for BNNs. One of the key motivating points is posterior uncertainty estimation. Hence important metric, when comparing to other posterior inference techniques, is to show high uncertainty for out of distribution samples and low for training/test data.