The empirical measure with samples from some probability measure (which might be known up to a multiplicative factor) has many applications in Bayesian inference [5; 3] and data assimilation [17]. A class of widely used sampling methods is the Markov chain Monte Carlo (MCMC) methods, where the trajectory of a particle is given by some constructed Markov chain with the desired distribution invariant. The trajectory of the particle is clearly stochastic, and the Monte Carlo methods take effect slowly for small number of samples. Unlike MCMC, the Stein variational gradient method (proposed by Liu and Wang in [20]) belongs to particlebased variational inference sampling methods (see also [22; 9]). These methods update particles by solving optimization problems, and each iteration is expected to make progress. As a nonparametric variational inference method, SVGD gives a deterministic way to generate points that approximate the desired probability distribution by solving an ODE system. Suppose that we are interested in some target probability distribution with density π(x) ∝ exp(-V (x)) (x ∈ ޒ d ). In SVGD, one sets V = -log π, chooses some symmetric positive definite kernel (x, y), 38 LEI LI, YINGZHOU LI, JIAN-GUO LIU, ZIBU LIU AND JIANFENG LU and solves the following ODE system for given initial points {X i (0)} N i=1 [20; 19]:
(X i , X j )∇V (X j ), i = 1, . . . , N , (
where N is the number of particles for the sampling purpose. The subindex "y" in ∇ y means that the gradient is taken with respect to the second variable in ( • , • ); i.e., ∇ y (X i , X j ) := ∇ y (x, y)| (x,y)=(X i ,X j ) . When t is large enough, the empirical measure constructed using {X i (t)} N i=1 is expected to be close to π, i.e.,
1
Below, in Section 2, we will explain why this is expected to be true. Theoretic understanding of (1-1) is limited. For example, the convergence of the particle system (1-1) is still open. Recently, there have been a few attempts at understanding the limiting mean field PDE [19; 21]. In particular, Lu et al. [21] showed the convergence of the mean field PDE to the desired measure π.
In practice, SVGD seems to perform quite well, better compared with some typical Monte Carlo methods in some examples [19; 10]. It provides consistent estimation for generic distributions as Monte Carlo methods do, but with fewer samples. SVGD seems to be more efficient than some Monte Carlo methods in the particle level for approximating the desired measure, when the number of particles is small. Interestingly, it reduces to the maximum a posterior (MAP) method when N = 1 [20].
Though (1-1) behaves well when the particle number N is not very big, one sometimes still needs an efficient algorithm to simulate (1-1). For example, when the dimension of the problem is not very high, in a typical MCMC method, the number of particles is several millions, or N ≈ 10 6 , while in SVGD, one may have N ≈ 10 3 . Simulating (1-1) needs O(N 2 ) work to compute the interactions for each iteration, especially for interaction kernels that are not superlocalized or particles that are not sparse. In fact, for such situations, to compute the interaction force for one particle, one must consider all the other N -1 particles to have enough accuracy. There are N particles, so one must consider O(N 2 ) interactions, which yields the O(N 2 ) complexity for one iteration. Though N ≈ 10 2 -10 3 is not large, the O(N 2 ) complexity makes the cost of SVGD for these cases comparable with MCMC with larger number of particles. Hence, it is highly motivated to develop a cheap version of SVGD.
In this work, we propose RBM-SVGD, a stochastic version of SVGD for sampling from a given probability measure. The idea is very natural: we apply the random batch method in [16] to the interacting particle system (1-1). Note that in the random batch method, the "batch" refers to the set for computing the interaction forces, not to be confused with the "batch" of samples for computing gradient as in stochastic gradient descent (SGD). Of course, if V is the loss function corresponding to many samples, or the probability density in Bayesian inference corresponding to many observed data, the data-mini-batch idea can be used to compute ∇V in SVGD as well [20]. With the random batch idea for computing interaction, the complexity for each iteration now is only O(N ). Moreover, it inherits the advantages of SVGD (i.e., efficient for sampling when the number of particles is not large) since the random batch method is designed to approximate the particle system directly. In fact, we will prove that the one marginal of the random batch method converges to the one marginal of the interacting particle systems under Wasserstein-2 distance on a fixed time interval [0, T ]. Note that the behavior of randomness in RBM-SVGD is different from that in MCMC. In MCMC, the randomness is required to ensure that the desired probability is invariant under the transition. The randomness in RBM-SVGD is simply due to the batch for computing the interaction forces, which is mainly for speeding up the computation. Though this randomness is not essential for sampling from the invariant measure, it may have other benefits. For example, it may lead to better ergodic properties for the particle system.
We now give a brief introduction to the SVGD proposed in [20] and provide some discussions. The derivation here is a continuous counterpart of that in [20].
Assume that random variable X ∈ ޒ d has density p 0 (x). Consider some mapping -: ޒ d → ޒ d , and we denote the distribution of -(X ) by p := -# p 0 , which is called the push-forward of p 0 under -. The goal is to make -# p 0 closer to π(x) in some sense. The way to measure the closeness of measures in [20] is taken to be the Kullback-Leibler (KL) divergence, which is also known as the relative entropy, defined by
where dµ dν is the well known Radon-Nikodym derivative. In [20,Theorem 3.1], it is shown that the Gateaux differential of -→ G(-) := KL( p π) is given by
where S q associated with a probability density q is called the Stein operator given by
In fact, using the formula
(2-5)
The quantity δG δ-, φ can be understood as the directional derivative of G( • ) in the direction given by φ. The paring in the second term above is in L 2 ޒ( d ) sense.
Based on this calculation, we now consider a continuously varying family of mappingsτ with τ ≥ 0 and
Here, "•" means composition, i.e., for any given x, d dττ (x) = φ τ (τ (x)). In this sense x → X (τ ; x) :=τ (x) is the trajectory of x under this mapping; x can be viewed as the so-called Lagrangian coordinate as in fluid mechanics while φ τ is the flow field. We denote p τ := (τ ) # p 0 .
(2-6)
The idea is then to choose φ τ such that the functional τ → G(τ ) decays as fast as possible. Note that to optimize the direction, we must require the field to have bounded magnitude φ τ H ≤ 1, where H is some subspace of the functions defined on ޒ d . The optimized curve τ →τ is a constant-speed curve (in some manifold).
Hence, the problem is reduced to the optimization problem
It is observed in [20] that this optimization problem can be solved by a convenient closed formula if H is the so-called (vector) reproducing kernel Hilbert space (RKHS) [1; 2]. A (scalar) RKHS is a Hilbert space, denoted by Ᏼ, consisting of functions defined on some space (in our case = ޒ d ) such that the evaluation function f → E x ( f ) := f (x) is continuous for all x ∈ . There thus exists k x ∈ Ᏼ such that E x ( f ) = f, k x Ᏼ . Then the kernel (x, y) := k x , k y Ᏼ is symmetric and positive definite, meaning that n i=1 n j=1 (x i , x j )c i c j ≥ 0 for any x i ∈ and c i ∈ .ޒ Reversely, given any positive definite kernel, one can construct a RKHS consisting of functions f (x) of the form f (x) = (x, y)ψ(y) dµ(y) where µ is some suitably given measure on . For example, if µ is the counting measure, choosing ψ(y) = ∞ j=1 a j 1 x j (y) (a j ∈ )ޒ can recover the form of RKHS in [20]. All such constructions yield isomorphic RKHS as guaranteed by the Moore-Aronszajn theorem [1]. Now, consider a given µ and H = Ᏼ d to be the vector RKHS:
The inner product is defined as
This inner product therefore induces a norm f H = √ f, f H . Clearly, H consists of functions with • H to be finite. The optimization problem (2-7) can be solved by the Lagrange multiplier method
where dy means Lebesgue measure and φ(x) = ޒ d (x, y)ψ(y) dµ(y). Using
where S * π is given by
(2-10)
The ODE flow
gives the constant-speed optimal curve, so that the velocity is the unit vector in H along the gradient of G. Reparametrizing the curve t = t (τ ) so that dτ dt = 2λ and denoting ρ t := p τ (t) , then
Clearly, the curve oft is not changed by this reparametrization. Using (2-4), one finds that ρ satisfies the equation
Here, * f (x) := (x, y) f (y) dy. It is easy to see that exp(-V ) is invariant under this PDE. According to the explanation here, the right-hand side gives the optimal decreasing direction of KL divergence if the transport flow is measured by RKHS. Hence, one expects it to be the negation of gradient of KL divergence in the manifold of probability densities with metric defined through RKHS. Indeed, Liu made the first attempt to justify this in [19, §3.4].
The above theory has a little trouble for empirical measures because the KL divergence is simply infinity. For empirical measure, ∇ρ must be in the distributional sense. The good thing for RKHS is that we can move the gradient from ∇ρ onto the kernel (x, y) so that the flow (2)(3)(4)(5)(6)(7)(8)(9)(10)(11) becomes (1-1), which makes perfect sense. In fact, if (1-1) holds, the empirical measure is a measure solution to (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12) (by testing on smooth function ϕ) [21, Proposition 2.5]. Hence, the ODE system is justified in this level, and one expects that (1-1) will give an approximation for the desired density. The numerical tests in [20] indeed justify this expectation. In this sense, the ODE system is formally a gradient flow of KL divergence, though the KL divergence functional is infinity for empirical measures.
Typical examples of (x, y) include (x, y) = (αx • y + 1) m , Gaussian kernel (x, y) = e -|x-y| 2 /(2σ 2 ) for ޒ d , and (x, y) = (sin a(xy))/(π(xy)) for 1D space .ޒ By Bochner's theorem [25], if a function K has a positive Fourier transform, then
is a positive definite kernel. With this kernel, (1-1) becomes
as used in [21]. Both Gaussians and 1/|x| α with α ∈ (0, d) have positive Fourier transforms. The difference is that the Gaussian has a short range of interaction while the latter has a long range of interaction. One can smoothen 1/|x| α out by mollifying with Gaussian kernels, resulting in positive definite smooth kernels but with long-range interaction. Choosing localized kernels like Gaussians may have some issues in very high-dimensional spaces [12; 10]. Due to its simplicity, when the dimension is not very high, we choose Gaussian kernels in Section 4. As a further comment, one may consider other metrics to gauge the closeness of probability measures, such as Wasserstein distances. Also, one can consider other norms for φ and get gradient flows in different spaces. These variants have been explored by some authors already [18; 8]. In general, computing the Frechét derivatives in closed form for these variants seems not that easy.
Remark. If we optimize (2-7) for φ in L 2 ޒ( d ; ޒ d ) spaces, the flow is then given by
The corresponding PDE is
. This is in fact the case when we choose (x, y) = δ(xy). This PDE, however, will not make sense for empirical measures since ρ∇ρ is hard to justify (clearly, the equivalent ODE system has the same trouble). By using RKHS, the derivative on ∇ρ can be moved onto the kernel and then the ODE system makes sense.
In this section, we introduce the "random batch" or "mini-batch" idea, which has already appeared in many places, and recall the random batch method for simulating interacting particle systems in [16]. By applying the random batch method to (1-1), we obtain a new algorithm, called RBM-SVGD. The proof that RBM-SVGD is close to SVGD on finite time interval is given in Section 3.2. (g(z j ; ω), y j ),
where (z j , y j ) M j=1 are some given data set, g( • ; ω) is a model that takes z j as an input and gives some prediction to y j , and ( • , • ) is some function to gauge the discrepancy between g(z j ; ω) and y j . Hence, the problem is to find ω such that the discrepancy is small enough. Often ( • , • ) is a neural network so that computing the gradient is not easy. Hence, if one aims to find the minimizer using gradient descent, the computation cost is high. The idea of "mini-batch" or "random batch" is to choose a small random subset ξ of {1, 2, . . . , N }, and consider the unbiased random estimate
with B = |ξ |, the size of ξ . Using this unbiased estimation L ξ to replace the original true gradient ∇ L ξ ≈ ∇ L, one can form the so-called stochastic gradient descent (SGD) [4; 6]. Using a similar idea for Langevin dynamics, Welling and Teh obtained a Markov chain Monte Carlo method, called the stochastic gradient Langevin dynamics (SGLD), useful for Bayesian inference [29].
for m in 1 : N T do Divide {1, 2, . . . , pn} into n batches randomly.
for each batch Ꮿ q do Update X i (i ∈ Ꮿ q ) by solving the equation for t ∈ [t m-1 , t m ):
(3-4)
end for end for Algorithm 1. Random batch method without replacement.
Consider in general the interacting particle system of the form
Here, F(x, y) does not have to be symmetric, and also F(x, x) is not necessarily zero. It is desirable to develop some cheap random approximation to the interacting forces so that the one-step O(N 2 ) complexity can be reduced. One idea is to use the "random batch" idea, but how to develop the concrete "random batch" algorithm depends on the concrete applications. Regarding the interacting particle systems, Jin et al. proposed some random grouping approach to achieve this goal in [16].
Here, we adopt the random batch method in [16] to and then obtain a stochastic version method for the SVGD ODE system (1-1). For this reason, we explain the random batch method a little bit. Choose a time step η. We define time grid points t m = mη.
(3-3)
At t m , one divides the particles into groups randomly, and each group is called a "batch", and then turns on interactions inside batches only. As indicated in [16], the random division of the particles into n batches takes O(N ) operations (one can for example use random permutation). Depending on whether one does batches without or with replacement, one can have different versions (see Algorithms 1 and 2). For the ODEs in the algorithms, one can apply any suitable ODE solver. For example, one can use the forward Euler discretization if F is smooth like Gaussian kernels. If K is singular, one may take p = 2 and apply the splitting strategy in [16].
For the SVGD ODE system (1-1), the kernel F takes the form
Applying the random batch method to this special kernel and using any suitable ODE solvers, we get a class of sampling algorithms, which we will call RBM-SVGD. In for m in 1 : N T * (N / p) do Pick a set Ꮿ of size p randomly. Update X i (i ∈ Ꮿ) by solving the following with pseudotime s ∈ [s m-1 , s m ):
(3-5)
end for Algorithm 2. Random batch method with replacement.
for k in 0 : N T -1 do Divide {1, 2, . . . , pn} into n batches randomly.
for each batch Ꮿ q do For all i ∈ Ꮿ q ,
where
end for end for Algorithm 3. RBM-SVGD.
this work, we will focus on the ones without replacement. The one with forward Euler discretization (with possible variant step size) is shown in Algorithm 3. Clearly, the complexity is O( pN ) for each iteration.
Here, N T is the number of iterations and {η k } is the sequence of time steps, which play the same role as learning rate in SGD [4; 6]. For some applications, one may simply set η k = η 1 to be a constant and get relatively good results. However, in many high-dimensional problems, choosing η k to be constant may yield divergent sequences [23]. One may decrease η k to obtain convergent data sequences. For example, one may simply choose η k = 1/k as in SGD. Another frequently used strategy is the AdaGrad approach [11; 28].
. We now give convergence analysis regarding the timecontinuous version of RBM-SVGD on torus ޔ d (i.e., choosing the particular force (3-6) for Algorithm 1 and X i ∈ ޔ d ). The analysis in this section justifies the expectation that RBM-SVGD should give similar performance as the original SVGD, as confirmed by the numerical experiments in Section 4.
By "torus", we mean the domain is equipped with periodic boundary conditions. The derivation of SVGD clearly stays unchanged for the torus. The reason we consider a torus is that (1-1) is challenging to analyze in ޒ d because of the nonlocal effect of the external force. On the torus, all functions are smooth and bounded. Moreover, using bounded domains with periodic boundary condition can always approximate the problem in ޒ d in practice.
Consider the random force for z = (x 1 , . . . , x N ) ∈ ޔ N d defined by
where Ꮿ is the random batch that contains i in the random batch method. Correspondingly, the exact force is given by
Define the "noise" by
We have the following consistency result regarding the random batch.
Lemma 1. For given z = (x 1 , . . . , x N ) ∈ ޔ N d (or ޒ N d ), it holds that
Moreover, the second moment is given by
where
(3-12)
The proof is similar to that in [16], but we also attach it in Appendix A for convenience.
We recall that the Wasserstein-2 distance is given by [26]
where (µ, ν) is called the transport plan, consisting of all the joint distributions whose marginal distributions are µ and ν, respectively: i.e., for any Borel set We now state the convergence result for the time-continuous version of RBM-SVGD, where we recall that F(x, y) given by (3)(4)(5)(6). We use X to denote the process generated by the random algorithm while X is the process by (1-1). The particles are exchangeable if the initial values are sampled i.i.d. from the same distribution. Hence, the distributions of X i are the same, and we call this one-particle distribution the one marginal distribution, which is a probability measure in ޔ d or ޒ d . We denote it by µ (1) N for convenience. Similarly, we introduce the one marginal distribution for the particles generated by the random algorithm, denoted by μ(1)
N .
Theorem 2. Assume V and K are smooth on torus ޔ d . The initial data X 0 i are drawn independently from the same initial distribution. Given T > 0, there exists C(T ) > 0, such that
Consequently, the one marginals µ (1) N and μ(1) N are close under Wasserstein-2 distance:
Proof. In the proof below, the constant C will represent a general constant independent of N and p, but its concrete meaning can change for every occurrence.
Consider the corresponding two processes and t ∈ [t m-1 , t m ]:
and
Taking the difference and dotting with X i -X i , one has
where
Hence, introducing
Due to the smoothness of K and V on the torus, we easily find
where C is independent of N . Note that Ꮿ is not independent of X j (t) for t > t m-1 , so to continue we must consider conditional expectation. Let Ᏺ m-1 be the σalgebra generated by X i (τ ), X i (τ ) for τ ≤ t m-1 (including the initial data drawn independently) and the random division of the batches at t m-1 . Then (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14) directly implies almost surely that
Thus, defining the error process
Note that
The inside of the parentheses can be estimated as
The first term on the right-hand side then becomes u(t m-1 ) by Lemma 1. By (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16), it is clear that
Hence,
where C is independent of N . Since u(t m-1 ) ≤ Cu(t) + Cη 2 by (3-18), then
Letting Z = (X 1 , . . . , X N ), one sees easily that I 2 = χ i (Z (t)). Then, we find
In J 2 , Y i (t m-1 ) is independent of the random batch division at t m-1 . Then, Lemma 1 tells us that ޅJ 2 = 0.
Using (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14), we have
where C is related to the infinity norm of the variance of χ i (t). This is the main term in the local truncation error. Just as we did for I 1 ,
this is controlled by C u(t m-1 ) + Cη. Hence,
Using the fact that u(t m-1 ) ≤ u(t) + Cη, one eventually has that
Applying Grönwall's inequality, we find
The last claim for W 2 distance follows from the definition of W 2 .
Note that the one marginal µ (1) N (t) is the distribution of X i (t) for any i, which is deterministic. This should be distinguished from the empirical measure µ N = (1/N ) i δ(x -X i (t)) which is random. As can be seen from the proof, the main contribution in the local truncation error comes from the variance of the noise χ i .
As can be seen, the error bound is given by the square root of variance of the random force times √ η = √ T /N T with N T being the number of steps. Hence, the result is a type of law of large number convergence result (see [16] for more details). The bigness of the variance on one hand depends on the batch size as 1/( p -1) -1/(N -1), while on the other hand depends on the bigness of the interaction. As long as the variance is bounded, the convergence of random batch method is ensured.
One crucial part is that the bigness of the variance depends on the bigness of the interaction, instead of the range of the interaction. This means that the random batch version of the algorithm is particularly useful when the interaction has long range or when the particles are not sparse. In fact, if the interaction has short range and the particles are sparse, one can use some data structure like cell-list [13, Appendix F] to reduce the computation of the interactions from O(N 2 ) to O(N ). However, when the interaction has long range or is not sparse (like the case in the example in Section 4.2), those data structures cannot be used any more, and RBM-SVGD becomes useful: it can still reduce the cost from O(N 2 ) to O(N ).
As another observation, according to (3)(4)(5)(6)(7)(8)(9)(10)(11) and (3)(4)(5)(6)(7)(8)(9)(10)(11)(12), the bigness of the variance depends on the bigness of the interaction kernel. As long as the variance stays controlled, the convergence of RBM-SVGD to SVGD is guaranteed. In this sense, the range of the interaction kernel is not sensitive to RBM-SVGD, so it can intrinsically be used for kernels that have long range. The choice of kernels clearly affects the performance of SVGD, but it seems not so significant for RBM-SVGD to approximate SVGD. In other words, we expect RBM-SVGD to work well when the kernel is chosen such that SVGD behaves well. In fact, our experience in Section 4 confirmed this.
Remark. We believe the error bound in Theorem 2 can be made independent of T due to the intrinsic structure of SVGD discussed above in Section 2. Then RBM-SVGD can be used as the efficient sampling algorithm from the desired distribution π . Such long time estimates are often established by some contracting properties of the ODE flows, so one may want to find the intrinsic converging structure of (1-1). However, rigorously establishing such results seems nontrivial due to the nonlocal effects of the external forces (∇V terms).
We consider some test examples in [19] to validate RBM-SVGD algorithm and compare with the original SVGD algorithm. In particular, in a toy example for 1D Gaussian mixture, RBM-SVGD is proved to be effective in the sense that the particle system converges to the expected distribution with less running time than the original SVGD method. A more practical example, namely Bayesian logistic regression, is also considered to verify the effectiveness of RBM-SVGD on large data sets in high dimension. Competitive prediction accuracy is presented by RBM-SVGD, and less time is needed. Hence, RBM-SVGD seems to be a more efficient method.
All numerical results in this section are implemented with Matlab R2018a and performed on a machine with Intel Xeon CPU E5-1650v2 at 3.50 GHz with 64 GB memory.
4.1. 1D Gaussian mixture. As a first example, we use the Gaussian mixture probability in [20] for RBM-SVGD. The initial distribution is ᏺ(-10, 1), Gaussian with mean -10 and variance 1. The target density is given by the Gaussian mixture
The kernel for the RKHS is the Gaussian kernel where h is the bandwidth parameter. For a fair comparison with the numerical results in [20], we first reproduce their results using N = 100 particles and dynamic bandwidth parameter h = med 2 /(2 log N ), where med is the median of the pairwise distance between the current points. Since dynamic bandwidth is infeasible for RBM-SVGD, we produce the results with fixed bandwidth h = 2 for the comparison between SVGD and RBM-SVGD. The RBM-SVGD uses Algorithm 3 with initial step size 0.2 and the following step sizes generated from AdaGrad. Different batch sizes are tested to demonstrate the efficiency of RBM-SVGD. Numerical results are illustrated in Figure 1 with the same initial random positions of particles following an ᏺ(-10, 1) distribution.
As stated in [20], the difficulty lies in the strong disagreement between the initial density function and the target density π(x). According to the first and second rows in Figure 1, SVGD with and without the fixed bandwidth parameter capture the target density efficiently and the corresponding convergence behaviors are similar to each other. Reading from the last column of Figure 1, we observe that RBM-SVGD inherits the advantage of SVGD in the sense that it can conquer the challenge and also show compelling result with SVGD. When the batch size is small, e.g., p = 2 or p = 5, the estimated densities differ from that of SVGD, and according to our experience, the estimated densities are not very stable across several executions while, in theory, RBM-SVGD runs N / p times faster than SVGD. Hence, RBM-SVGD with p = 5 at the 500-th iteration costs the same as 50 iterations of SVGD. According to Figure 1, RBM-SVGD(2) at the 500-th iteration significantly outperforms the 50-th iteration of SVGD. As we increase the batch size, as in the last two rows of Figure 1, more stable and similar behavior to SVGD is observed.
Provided the good performance of RBM-SVGD, we also check the sampling power and its computational cost. We conduct the following simulations with N = 256 particles for 500 iterations with the Gaussian kernel (4-2). For RBM-SVGD, we use fixed bandwidth h = 2 whereas SVGD uses the aforementioned dynamic bandwidth strategy. When we apply SVGD or RBM-SVGD with different batch sizes, the same initial random positions of particles is used. For a given test function h(x), we compute the estimated expectation h = (1/N ) N i=1 h(X i (T )) and the sampling accuracy is measured via the minimum square error (MSE) over 100 random initializations following the same distribution as before:
where ޅ X ∼π h(X ) denotes the underlying truth. Three test functions are explored, h 1 (x) = x, h 2 (x) = x 2 , and h 3 (x) = cos 2x, with their corresponding true expectations being 2 3 , 5, and (cos 4)/e 2 . The reported run time is also averaged over 100 random initializations.
Figure 2 shows the MSE against different batch sizes for h 1 (x), h 2 (x), and h 3 (x), respectively. The results of RBM-SVGD with different batch sizes are connected by lines, whereas the results of SVGD are the isolated points with batch size p = 256. In general, the estimations of h 1 (x) and h 2 (x) are better than that of h 3 (x), which agrees with the difficulty of the problems. Table 1 shows the averaged run time of RBM-SVGD and SVGD for different batch sizes under two different implementations in Matlab. RBM-SVGD is faster than SVGD for all choices of batch size. With respect to the two implementations in Matlab, for the first block row, within each batch, a matrix operation is adopted in computing the kernel matrix whereas for the second block row, the kernel matrix is computed row by row. Matlab naturally is more favorable in the first implementation, which hence achieves fastest run time for all different batch sizes. For other programming languages, e.g., C++, Fortran, etc., the speedup of the second block row is excepted, which is close to ideal case as we predicted earlier.
In this section, we will compare RBM-SVGD with MCMC, specifically Metropolis-Hastings (MH) [15]. The algorithmic detail of Metropolis-Hastings is available in Appendix B. The performance of RBM-SVGD and MH on a Bayesian inference task is compared to illustrate the advantage of RBM-SVGD. When the number of particles is not very large and desired accuracy is not high, RBM-SVGD can be more efficient.
We run MH and RBM-SVGD on a Bayesian inference task which is exactly the experiment in [10]. In this inference problem, our unknown parameter x is in ޒ 2 . The observational data y is a real number which is determined by the forward map Ᏺ(x) and the observational noise, i.e., y = Ᏺ(x) + ξ , where the forward map is a scalar logarithmic Rosenbrock function [24] Ᏺ(x) = log((1 -x 1 ) 2 + 100(x 2 -x 2 1 ) 2 ) for x = (x 1 , x 2 ) and the Gaussian noise ξ satisfies ξ ∼ ᏺ(0, σ 2 ) for σ = 0.3. The relationship between parameter x and observation y implies that the likelihood function is p(y | x) = ᏺ(F(x), σ 2 ). Finally, we set the prior distribution for x to be Gaussian, i.e., π 0 (x) = ᏺ(0, τ 2 I 2 ), where I 2 is the identity matrix and τ will be specified later. Thus, the unnormalized posterior density is given by
N = 512 particles are sampled in RBM-SVGD, and the maximum iteration number is 800. Different batch sizes are tested for performance, and the bandwidth parameter is fixed to be h = 0.1. To make MH comparable with RBM-SVGD regarding the number of sampling points, we viewed MH as a method with batch size 1, so the total number of iterations we performed for MH was N • 200. We apply burn-in technique by only considering the second-half iterations. To reduce correlation, only 1 sample is drawn from every 100 iterations. Therefore, a total number of N samples are selected from MH, which agrees with the number of particles we employ in RBM-SVGD. According to the performance test in Appendix B, we compare RBM-SVGD with MH by choosing τ = 5 • 10 -3 , which is tested to be convergent and presents the best visual performance among different choices of τ . For both RBM-SVGD and MH, the initial points are sampled from a Gaussian distribution ᏺ(0, 0.4 2 ). The target distribution is double banana with centers near (0, 0.5) and (0, -0.5). Hence, we adopt two test functions as
)). In Figure 3, we plotted the position of each particle after RBM-SVGD iteration or MH together with the contour map of the target distribution. From the picture we can tell that both MH and RBM-SVGD can recover the shape of the target density and produce persuasive samplings. Although RBM-SVGD slightly harmed the aggregation of particles around the true distribution (which also paid off with a much shorter running time) compared to the original SVGD (RBM-SVGD with batch size = 512), it can still provide a convincing sampling by almost recovering the shape of the target density. In Table 2, we give further quantitative comparison. All numbers in the table are averaged over 100 different initializations. The run time for any RBM-SVGD with different batch sizes is faster than that of MH, and RBM-SVGD with batch size 2 is more than 20× faster while, regarding the MSE for both h 1 and h 2 , RBM-SVGD is much better than MH for h 1 and better than MH for h 2 . Hence, we conclude, for this example, SVGD outperforms MH both in run time and accuracy. RBM-SVGD further significantly reduces the run time of regular SVGD without loss of accuracy.
In this experiment, we apply RBM-SVGD to conduct Bayesian logistic regression for binary classification for the Covertype data set with 581012 data points and 54 features [14]. Under the same setting as Gershman [14; 20], the regression weights w of dimension 54 are assigned with a Gaussian prior p 0 (ω | α) = ᏺ(w, α -1 ), and the variance satisfies p 0 (α, 1, 0.01), where represents the density of Gamma distribution. The inference is applied on posterior p(x | D) with x = [w, log α] of dimension 55. The kernel K ( • ) is taken again to be the same Gaussian kernel as (4-2).
Since the problem is in high dimension, we adopt N = 512 particles in this experiment, which also create more space for the selection of batch sizes. The training is done on 80% of the data set, and the other 20% is used as the test data set. For particle system (1-1), the computation of -∇V = ∇ log p(x) is expensive. Hence, we use the same strategy as mentioned in [20, §3.2], i.e., using data-mini-batch 1 of the data to form a stochastic approximation of p(x) with the data-mini-batch size being 100. Since ∇ log p depends only on x as in Algorithm 3, 1 To avoid confusion with our batch of particles, we call it data-mini-batch instead. As shown in Figure 4, RBM-SVGD is almost as efficient as SVGD even for small batch sizes. When p = 2, the test accuracy converges to a value slightly off that of SVGD. RBM-SVGD with p = 4 converges to the same accuracy as SVGD but at a slower convergent rate. For RBM-SVGD with batch size greater than 4, we observe similar convergence behavior as that of SVGD. The run time of RBM-SVGD, as shown in Table 3, is lower than that of SVGD, where the run time of 6000 iterations is reported. Comparing to the similar run time table for the 1D Gaussian mixture example (Table 1), the acceleration of RBM-SVGD is not as significant as before. This is due to the linear but expensive evaluation of ∇ log p, where RBM-SVGD and SVGD spend the same amount of time in the evaluation each iteration. Although the evaluation of ∇ log p is expensive, it is linear in N . As N increases, the advantage of RBM-SVGD would be more significant. In Table 4, we list the mean and standard deviation of RBM-SVGD with p = 2 and p = 8 and SVGD of different iterations. Based on the statistics, we conclude that RBM-SVGD and SVGD are of similar prediction power and RBM-SVGD is efficient also in high-dimensional particle systems as well.
We have applied the random batch method for interacting particle systems to SVGD, resulting in RBM-SVGD, which turns out to be a cheap sampling algorithm and inherits the efficiency of the original SVGD algorithm. Theory and numerical experiments have validated the algorithm, and hence, it can potentially have many applications, like Bayesian inference. Moreover, as a hybrid strategy, one may increase the batch size as time goes on to increase the accuracy, or apply some variance reduction approach. Given the prior and target distribution as in (4-3), we first test the performance of MH among different values of parameter τ , which represents the variance of the proposal distribution.
Figure 5 illustrates the samples together with the contour map of target distribution for τ = 5 • 10 -1 , 5 • 10 -2 , 5 • 10 -3 , 5 • 10 -4 . Clearly, the best performance was attained when τ = 5 • 10 -3 . For τ greater than 5 • 10 -3 , samples are still wandering around the true distribution without accumulating due to a high variance, whereas for smaller τ , samples are merely aggregating around the "lower" banana rather than the "upper" banana. This phenomenon can be explained by the small variance of proposal distribution, which confines the particles around upper banana.
Moreover, a convergence diagnosis for MCMC was also conducted by computing the auto-correlation [7]. A lower auto-correlation always implies better convergence because a higher auto-correlation indicates that effective sampling size is smaller and more iteration is necessary [27]. Figure 6