A stochastic version of Stein Variational Gradient Descent for efficient sampling

TL;DR

This paper introduces RBM-SVGD, a stochastic variant of Stein Variational Gradient Descent that leverages the Random Batch Method to reduce computational costs in sampling tasks, especially with long-range kernels.

Contribution

The paper develops RBM-SVGD, integrating RBM into SVGD to improve efficiency for Bayesian sampling without sacrificing core behaviors.

Findings

RBM-SVGD reduces computational cost compared to standard SVGD.

Numerical experiments demonstrate the efficiency of RBM-SVGD.

RBM-SVGD maintains the effectiveness of SVGD in sampling tasks.

Abstract

We propose in this work RBM-SVGD, a stochastic version of Stein Variational Gradient Descent (SVGD) method for efficiently sampling from a given probability measure and thus useful for Bayesian inference. The method is to apply the Random Batch Method (RBM) for interacting particle systems proposed by Jin et al to the interacting particle systems in SVGD. While keeping the behaviors of SVGD, it reduces the computational cost, especially when the interacting kernel has long range. Numerical examples verify the efficiency of this new version of SVGD.

Figures

Tables

Table 1. Table 1 : Averaged runtime for different batch sizes.

	RBM-SVGD							SVGD
Batch size	2	4	8	16	32	64	128	256
Runtime(s)	0.055	0.095	0.178	0.341	0.270	0.238	0.314	0.733
Speedup	13.3x	7.7x	4.1x	2.1x	2.7x	3.1x	2.3x

Table 2. Table 2 : Average runtime of 6000 6000 6000 iterations

	RBM-SVGD						SVGD
Batch size	2	4	8	16	32	128	512
Runtime(s)	8.59	11.24	16.28	26.15	21.66	19.42	47.01
Speedup	5.5x	4.2x	2.9x	1.8x	2.2x	2.4x

Table 3. Table 3 : Statistics of RBM-SVGD and SVGD.

	Iteration	1000	2000	3000	4000	5000	6000
RBM-SVGD $p=2$	Mean	0.7090	0.7349	0.7409	0.7446	0.7457	0.7471
	Std	0.0045	0.0040	0.0040	0.0034	0.0034	0.0038
RBM-SVGD $p=8$	Mean	0.7342	0.7470	0.7508	0.7518	0.7527	0.7534
	Std	0.0073	0.0056	0.0041	0.0045	0.0039	0.0033
SVGD	Mean	0.7347	0.7530	0.7523	0.7529	0.7504	0.7511
	Std	0.0068	0.0048	0.0071	0.0048	0.0061	0.0062

Equations

$${\dot{X}}_i=\frac{1}{N}\sum _{j=1}^N{\nabla }_y\mathcal{K}(X_i,X_j)-\frac{1}{N}\sum _{j=1}^N\mathcal{K}(X_i,X_j)\nabla V(X_j),$$

$$\mathrm{KL}(\mu ||\nu )={\mathbb{E}}_{Y\sim \mu }\log \left(\frac{d\mu }{d\nu }(Y)\right),$$

$$⟨\frac{\delta G}{\delta T},\varphi ⟩=-{\mathbb{E}}_{Y\sim p}{S}_{\pi }\varphi (Y),\forall \varphi \in C_c^{\infty }({\mathbb{R}}^d;{\mathbb{R}}^d)$$

$$S_q\varphi (x)=\nabla (\log q(x))\cdot \varphi (x)+\nabla \cdot \varphi (x).$$

$$\frac{d}{dϵ}(T+ϵ\varphi \circ T{)}_{\#}{p}_0{|}_{ϵ=0}=\frac{d}{dϵ}(I+ϵ\varphi {)}_{\#}p{|}_{ϵ=0}=-p{S}_p\varphi =-\nabla \cdot (p\varphi ),$$

$$⟨\frac{\delta G}{\delta T},\varphi ⟩=⟨\frac{\delta \mathrm{KL}(p||\pi )}{\delta p},-\nabla \cdot (p\varphi )⟩=-{\int }_{{\mathbb{R}}^d}p{S}_{\pi }\varphi dx.$$

$$\frac{d}{d\tau }{T}_{\tau }={\varphi }_{\tau }\circ T_{\tau }.$$

$$p_{\tau }:=(T_{\tau }{)}_{\#}{p}_0.$$

$$\sup \{{\mathbb{E}}_{Y\sim p}{S}_{\pi }\varphi (Y)|\Vert \varphi {\Vert }_H\le 1\}.$$

$$H=\left\{f=\int_{\mathbb{R}^{d}}\mathcal{K}(\cdot,y)\psi(y)\,d\mu(y)\Big{|}\psi:\mathbb{R}^{d}\to\mathbb{R}^{d},\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\mathcal{K}(x,y)\psi(x)\cdot\psi(y)d\mu(x)d\mu(y)<\infty\right\}.$$

$$\begin{array}{rl}⟨f^{(1)},f^{(2)}{⟩}_H& ={\iint }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\mathcal{K}(x,y){\psi }^{(1)}(x)\cdot {\psi }^{(2)}(y)d\mu (x)d\mu (y)\\ & =\sum _{j=1}^d{\iint }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\mathcal{K}(x,y){\psi }_j^{(1)}(x){\psi }_j^{(2)}(y)d\mu (x)d\mu (y).\end{array}$$

$$\mathcal{L}={\int }_{{\mathbb{R}}^d}(S_{\pi }\varphi )p_{\tau }(y)dy-\lambda {\iint }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\mathcal{K}(x,y)\psi (x)\cdot \psi (y)d\mu (x)d\mu (y),$$

$$2\lambda \varphi ={\int }_{{\mathbb{R}}^d}\mathcal{K}(x,y)(S_{\pi }^{\ast }{p}_t)(y)dy=:\mathcal{V}(p_t),$$

$$S_{\pi }^{\ast }(f)=f(y)\nabla (\log \pi )-\nabla f(y)=-f(y)\nabla V(y)-\nabla f(y).$$

$$\frac{d}{d\tau }{T}_{\tau }=\frac{1}{2\lambda (\tau )}\mathcal{V}(p_{\tau })\circ T_{\tau },$$

$$\frac{d}{dt}{T}_t=\mathcal{V}({\rho }_t)\circ T_t.$$

$${\partial }_t\rho =-\nabla \cdot (\mathcal{V}(\rho )\rho )=\nabla \cdot (\rho \mathcal{K}\ast (\rho \nabla V+\nabla \rho )).$$

$$\mathcal{K}(x,y)=K(x-y)$$

$${\dot{X}}_i=-\frac{1}{N}\sum _{j=1}^N\nabla K(X_i-X_j)-\frac{1}{N}\sum _{j=1}^{N}K(X_i-X_j)\nabla V(X_j),$$

$$\frac{d}{dt}T=(S_{\pi }^{\ast }\rho )\circ T.$$

$${\dot{X}}_i=\frac{1}{N}\sum _{j=1}^{N}F(X_i,X_j)=\frac{1}{N}F(X_i,X_i)+\frac{1}{N}\sum _{j:j\ne i}F(X_i,X_j).$$

$$t_m=m\eta .$$

$${\dot{X}}_i=\frac{1}{N}F(X_i,X_i)+(1-\frac{1}{N})\frac{1}{p-1}\sum _{j\in {\mathcal{C}}_q,j\ne i}F(X_i,X_j).$$

$${\dot{X}}_i=\frac{1}{N}F(X_i,X_i)+(1-\frac{1}{N})\frac{1}{p-1}\sum _{j\in \mathcal{C},j\ne i}F(X_i,X_j).$$

$$F(x,y)={\nabla }_y\mathcal{K}(x,y)-\mathcal{K}(x,y)\nabla V(y).$$

$$X_{i}^{(k+1)}\leftarrow X_{i}^{(k)}+\frac{1}{N}\Big{(}\nabla_{y}\mathcal{K}(X_{i}^{(k)},X_{i}^{(k)})-\mathcal{K}(X_{i}^{(k)},X_{i}^{(k)})\nabla V(X_{i}^{(k)})\Big{)}\eta_{k}+\Phi_{k,i}\eta_{k},$$

$${\Phi }_{k,i}=\frac{N-1}{N(p-1)}\sum _{j\in {\mathcal{C}}_q,j\ne i}\left({\nabla }_y\mathcal{K}(X_i^{(k)},X_j^{(k)})-\mathcal{K}(X_i^{(k)},X_j^{(k)})\nabla V(X_j^{(k)})\right).$$

$$f_i(z):=\left(1-\frac{1}{N}\right)\frac{1}{p-1}\sum _{j:j\in \mathcal{C}}F(x_i,x_j),$$

$${\chi }_i(z):=\frac{1}{N}\sum _{j:j\ne i}F(x_i,x_j)-f_i(z).$$

$$\mathbb{E}{\chi }_i(z)=0$$

Full text

1 Introduction

The empirical measure with samples from some probability measure (which might be known up to a multiplicative factor) has many applications in Bayesian inference [1, 2] and data assimilation [3]. 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 [4]) belongs to particle based variational inference sampling methods (see also [5, 6]). These methods update particles by solving optimization problems, and each iteration is expected to make progress. As a non-parametric 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 $\pi(x)\propto\exp(-V(x))$ ($x\in\mathbb{R}^{d}$). In SVGD, one sets $V=-\log\pi$ and solve the following ODE system for given initial points $\{X_{i}(0)\}$ (see [4, 7]):

[TABLE]

where $\mathcal{K}(x,y)$ is a symmetric positive definite kernel. When $t$ is large enough, the empirical measures constructed using $\{X_{i}(t)\}$ is expected to be close to $\pi$.

SVGD seems to be more efficient in the particle level for approximating the desired measure and interestingly, it reduces to the maximum a posterior (MAP) method when $N=1$ [4]. It provides consistent estimation for generic distributions as Monte Carlo methods do, but with fewer samples. Theoretic understanding of (1.1) is limited. For example, the convergence of the particle system (1.1) is still open. Recently, there are a few attempts for the understanding of the limiting mean field PDE [7, 8]. In particular, Lu et al [8] showed the convergence of the PDE to the desired measure.

Though (1.1) behaves well when the particle number $N$ is not very big, one sometimes still needs efficient algorithm to simulate (1.1). For example, in a typical MCMC method $N=10^{4}\sim 10^{6}$ while in SVGD, one may have $N=10^{2}\sim 10^{3}$. Though $N=10^{2}\sim 10^{3}$ is not large, simulating (1.1) needs $O(N^{2})$ work to compute the interactions for each iteration, especially for interaction kernels that are not super localized (such as kernels with algebraic decaying rate, like $K(x)\sim|x|^{-\alpha}$). The computation cost of SVGD for these cases is therefore 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 [9] 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 $\nabla V$ in SVGD as well (see [4]). 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 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 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 particle system.

2 Mathematical background of SVGD

We now give a brief introduction to the SVGD proposed in [4] and make some discussion. The derivation here is a continuous counterpart of that in [4].

Assume that random variable $X\in\mathbb{R}^{d}$ has density $p_{0}(x)$. Consider some mapping $T:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and we denote the distribution of $T(X)$ by $p:=T_{\#}p_{0}$, which is called the push-forward of $p_{0}$ under $T$. The goal is to make $T_{\#}p_{0}$ closer to $\pi(x)$ in some sense. The way to measure the closeness of measures in [4] is taken to be the Kullback-Leibler (KL) divergence, which is also known as the relative entropy, defined by

[TABLE]

where $\frac{d\mu}{d\nu}$ is the well-known Radon-Nikodym derivative. In [4, Theorem 3.1], it is shown that the Frechet differential of $T\mapsto G(T):=\operatorname{KL}(p||\pi)$ is given by

[TABLE]

where $S_{q}$ associated with a probability density $q$ is called the Stein operator given by

[TABLE]

In fact, using the formula

[TABLE]

and $\frac{\delta\operatorname{KL}(p||\pi)}{\delta p}=\log p-\log\pi$, one finds

[TABLE]

The quantity $\langle\frac{\delta G}{\delta T},\phi\rangle$ can be understood as the directional derivative of $G(\cdot)$ in the direction given by $\phi$.

Based on this calculation, we now consider a continuously varying family of mappings $T_{\tau},\tau\geq 0$ and

[TABLE]

Here, ’$\circ$’ means composition, i.e. for any given $x$, $\frac{d}{d\tau}T_{\tau}(x)=\phi_{\tau}(T_{\tau}(x))$. In this sense $x\mapsto X(\tau;x):=T_{\tau}(x)$ is the trajectory of $x$ under this mapping; $x$ can be viewed as the so-called Lagrangian coordinate as in fluid mechanics while $\phi_{\tau}$ is the flow field. We denote

[TABLE]

The idea is then to choose $\phi_{\tau}$ such that the functional $\tau\mapsto G(T_{\tau})$ decays as fast as possible. Note that to optimize the direction, we must impose the field to have bounded magnitude $\|\phi_{\tau}\|_{H}\leq 1$, where $H$ is some subspace of the functions defined on $\mathbb{R}^{d}$. The optimized curve $\tau\mapsto T_{\tau}$ is a constant speed curve (in some manifold). Hence, the problem is reduced to the following optimization problem

[TABLE]

It is observed in [4] that this optimization problem can be solved by a convenient closed formula if $H$ is the so-called (vector) reproducing kernel Hilbert space (RKHS) [10, 11]. A (scalar) RKHS is a Hilbert space, denoted by $\mathcal{H}$, consisting of functions defined on some space $\Omega$ (in our case $\Omega=\mathbb{R}^{d}$) such that the evaluation function $f\mapsto E_{x}(f):=f(x)$ is continuous for all $x\in\Omega$. There thus exists $k_{x}\in\mathcal{H}$ such that $E_{x}(f)=\langle f,k_{x}\rangle_{\mathcal{H}}$. Then the kernel $\mathcal{K}(x,y):=\langle k_{x},k_{y}\rangle_{\mathcal{H}}$ is symmetric and positive definite, meaning that $\sum_{i=1}^{n}\sum_{j=1}^{n}\mathcal{K}(x_{i},x_{j})c_{i}c_{j}\geq 0$ for any $x_{i}\in\Omega$ and $c_{i}\in\mathbb{R}$. Reversely, given any positive definite kernel, one can construct a RKHS consisting of functions $f(x)$ of the form $f(x)=\int\mathcal{K}(x,y)\psi(y)\,d\mu(y)$ where $\mu$ is some suitably given measure on $\Omega$. For example, if $\mu$ is the counting measure, choosing $\psi(y)=\sum_{j=1}^{\infty}a_{j}1_{x_{j}}(y)$ ($a_{j}\in\mathbb{R}$) can recover the form of RKHS in [4]. All such constructions yield isomorphic RKHS as guaranteed by Moore-Aronszajn theorem [10]. Now, consider a given $\mu$ and $H=\mathcal{H}^{d}$ to be the vector RKHS:

[TABLE]

The inner product is defined as

[TABLE]

This inner product therefore induces a norm $\|f\|_{H}=\sqrt{\langle f,f\rangle_{H}}$. Clearly, $H$ consists of functions with $\|\cdot\|_{H}$ to be finite. The optimization problem (2.7) can be solved by the Lagrange multiplier method

[TABLE]

where $dy$ means Lebesgue measure and $\phi(x)=\int_{\mathbb{R}^{d}}\mathcal{K}(x,y)\psi(y)\,d\mu(y)$. Using $\frac{\delta\mathcal{L}}{\delta\phi}=0$, we find

[TABLE]

where $S_{\pi}^{*}$ is given by

[TABLE]

The ODE flow

[TABLE]

gives the constant speed optimal curve, so that the velocity is the unit vector in $H$ along the gradient of $G$. Re-parametrizing the curve $t=t(\tau)$ so that $\frac{d\tau}{dt}=2\lambda$, and we denote $\rho_{t}:=p_{\tau(t)}$, then

[TABLE]

Clearly, the curve of $T_{t}$ is not changed by this reparametrization. Using (2.4), one finds that $\rho$ satisfies the following equation

[TABLE]

Here, $\mathcal{K}*f(x):=\int\mathcal{K}(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 [7, Sec. 3.4].

While everything looks great for continuous probability densities, the above theory does not work for empirical measures because the KL divergence is simply infinity. For empirical measure, $\nabla\rho$ must be in the distributional sense. However, the good thing for RKHS is that we can move the gradient from $\nabla\rho$ onto the kernel $\mathcal{K}(x,y)$ so that the flow (2.11) becomes (1.1), which makes perfect sense. In fact, if (1.1), holds, the empirical measure is a measure solution to (2.12) (by testing on smooth function $\varphi$) [8, Proposition 2.5]. Hence, one expects that (1.1) will give approximation for the desired density. The numerical tests in [4] 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 $\mathcal{K}(x,y)$ include $\mathcal{K}(x,y)=(\alpha x\cdot y+1)^{m}$, Gaussian kernel $\mathcal{K}(x,y)=e^{-|x-y|^{2}/(2\sigma^{2})}$ for $\mathbb{R}^{d}$, and $\mathcal{K}(x,y)=\frac{\sin a(x-y)}{\pi(x-y)}$ for 1D space $\mathbb{R}$. By Bochner’s theorem [12], if a function $K$ has a positive Fourier transform, then

[TABLE]

is a positive definite kernel. With this kernel, (1.1) becomes

[TABLE]

as used in [8]. Both Gaussians and $1/|x|^{\alpha}$ with $\alpha\in(0,d)$ have positive Fourier transforms. The difference is that Gaussian has short range of interaction while the latter has long range of interaction. One can smoothen $1/|x|^{\alpha}$ 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 [13, 14]. 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 metric to gauge the closeness of probability measures, such as Wasserstein distances. Also, one can consider other norms for $\phi$ and get gradient flows in different spaces. These variants have been explored by some authors already [15, 16]. In general, computing the Frechet derivatives in closed form for these variants seems not that easy.

Remark 1**.**

If we optimize (2.7) for $\phi$ in $L^{2}(\mathbb{R}^{d};\mathbb{R}^{d})$ spaces, the flow is then given by

[TABLE]

The corresponding PDE is $\partial_{t}\rho=\nabla\cdot(\rho(\rho\nabla V+\nabla\rho))=\nabla\cdot(\rho^{2}\nabla\log\frac{\rho}{\pi})$. This is in fact the case when we choose $\mathcal{K}(x,y)=\delta(x-y)$. This PDE, however, will not make sense for empirical measures since $\rho\nabla\rho$ is hard to justify (Clearly, the equivalent ODE system has the same trouble.) By using RKHS, the derivative on $\nabla\rho$ can be moved onto the kernel and then the ODE system makes sense.

3 The new sampling algorithm: RBM-SVGD

We consider in general the particle system of the following form.

[TABLE]

Here, $F(x,y)$ does not have to be symmetric, and also $F(x,x)$ is not necessarily zero.

3.1 The algorithms

We apply the random batch method in [9] to this particle system. In particular, choose a time step $\eta$. We define time grid points

[TABLE]

The idea of random batch method is to form some random batches at $t_{m}$ and then turn on interactions inside batches only. As indicated in [9], the random division of the particles into $n$ batches takes $O(N)$ operations (we can for example use random permutation). Depending on whether we do batches without or with replacement, we can have different versions (see Algorithm 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 [9].

For the Stein Variational Gradient Descent (1.1), the kernel $F$ takes the following form.

[TABLE]

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 this work, we will mainly focus on the ones without replacement. Some discussion for RBM-SVGD with or without replacement will be made in section 4.2. 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 $\{\eta_{k}\}$ is the sequence of time steps, which play the same role as learning rate in stochastic gradient descent (SGD). For some applications, one may simply set $\eta_{k}=\eta\ll 1$ to be a constant and gets reasonably good results. However, in many high dimensional problems, choosing $\eta_{k}$ to be constant may yield divergent sequences [17]. One may decreas $\eta_{k}$ to obtain convergent data sequences. For example, one may simply choose $\eta_{k}=\eta\ll 1$ as in SGD. Another frequently used strategy is the Adagrad approach [18, 19].

3.2 Theoretic results

We now give convergence analysis regarding the time continuous version of RBM-SVGD on torus $\mathbb{T}^{d}$ (i.e. choosing the particular force (3.5) for Algorithm 1 and $X_{i}\in\mathbb{T}^{d}$). The derivation of SVGD clearly stays unchanged for torus. The reason we consider torus is that (1.1) is challenging to analyze in $\mathbb{R}^{d}$ because of the nonlocal effect of the external force. On torus, all functions are smooth and bounded. Moreover, using bounded domains with periodic boundary condition can always approximate the problem in $\mathbb{R}^{d}$ in practice.

Consider the random force for $z=(x_{1},\ldots,x_{N})\in\mathbb{T}^{Nd}$ defined by

[TABLE]

where $\mathcal{C}$ is the random batch that contains $i$ in the random batch method. Correspondingly, the exact force is given by $F_{i}(z)=\frac{1}{N}\sum_{j:j\neq i}F(x_{i},x_{j})$. Define the ’noise’ by

[TABLE]

We have the following consistency result regarding the random batch.

Lemma 1**.**

For given $z=(x_{1},\ldots,x_{N})\in\mathbb{T}^{Nd}$ (or $\mathbb{R}^{Nd}$), it holds that

[TABLE]

Moreover, the second moment is given by

[TABLE]

where

[TABLE]

The proof is similar as in [9], but we also attach it in the Appendix A for convenience.

We recall that the Wasserstein-$2$ distance is given by [20]

[TABLE]

where $\Pi(\mu,\nu)$ is called the transport plan, consisting of all the joint distributions whose marginal distributions are $\mu$ and $\nu$ respectively: i.e. for any Borel set $E\subset\mathbb{T}^{d}$, $\mu(E)=\iint_{x\in E,y\in\mathbb{T}^{d}}\gamma(dx,dy)$ and $\nu(E)=\int_{x\in\mathbb{T}^{d},y\in E}\gamma(dx,dy)$.

We now state the convergence result for the time continuous version of RBM-SVGD, where $F(x,y)=\nabla_{y}\mathcal{K}(x,y)-\mathcal{K}(x,y)\nabla V(y)$. We use $\tilde{X}$ to denote the process generated by the random algorithm while $X$ is the process by (1.1).

Theorem 1**.**

Assume $V$ and $K$ are smooth on torus $\mathbb{T}^{d}$. The initial data $X_{i}^{0}$ are drawn independently from the same initial distribution. Given $T>0$, there exists $C(T)>0$, such that $\mathbb{E}|X_{i}-\tilde{X}_{i}|^{2}\leq C(T)\eta$. Consequently, the one marginals $\mu_{N}^{(1)}$ and $\tilde{\mu}_{N}^{(1)}$ are close under Wasserstein-$2$ distance:

[TABLE]

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\in[t_{m-1},t_{m}]$.

[TABLE]

and

[TABLE]

Taking the difference and dotting with $\tilde{X}_{i}-X_{i}$, one has

[TABLE]

where

[TABLE]

[TABLE]

Hence, introducing

[TABLE]

we have

[TABLE]

Due to the smoothness of $K$ and $V$ on torus, we easily find

[TABLE]

where $C$ is independent of $N$. Note that $\mathcal{C}$ is not independent of $X_{j}$ for $t>t_{m-1}$, so to continue we must consider conditional expectation. Let $\mathcal{F}_{m-1}$ be the $\sigma$-algebra generated by $X_{i}(\tau),\tilde{X}_{i}(\tau)$ for $\tau\leq t_{m-1}$ (including the initial data drawn independently) and the random division of the batches at $t_{m-1}$. Then, (3.13) directly implies almost surely it holds that

[TABLE]

Thus, defining the error process

[TABLE]

we have

[TABLE]

yielding

[TABLE]

Note that

[TABLE]

The inside of the parenthesis can be estimated as

[TABLE]

The first term on the right hand side then becomes $u(t_{m-1})$ by Lemma 1. By (3.15), it is clear that

[TABLE]

Hence,

[TABLE]

where $C$ is independent of $N$. Since $u(t_{m-1})\leq Cu(t)+C\eta^{2}$ by (3.18), then

[TABLE]

Letting $Z=(X_{1},\ldots,X_{N})$, one sees easily that $I_{2}=\chi_{i}(Z(t))$. Then, we find

[TABLE]

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

[TABLE]

Using (3.13), we have

[TABLE]

Since $\chi_{i}$ is bounded,

[TABLE]

where $C$ is related to the infinity norm of the variance of $\chi_{i}(t)$. This is the main term in the local truncation error. Just as we did for $I_{1}$,

[TABLE]

Since

[TABLE]

This is controlled by $C\sqrt{u(t_{m-1})}+C\eta$. Hence,

[TABLE]

where the $\eta$ term is from the variance term.

Eventually,

[TABLE]

Applying Grönwall’s inequality, we find

[TABLE]

The last claim for $W_{2}$ distance follows from the definition of $W_{2}$. ∎

Note that the one marginal $\mu_{N}^{(1)}$ is the distribution of $X_{i}$ for any $i$, which is deterministic. This should be distringuished from the empirical measure $\mu_{N}=\frac{1}{N}\sum_{i}\delta(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 the noise $\chi_{i}$. We believe the error bound here can be made independent of $T$ due to the intrinsic structure of SVGD discussed above in section 2. Often, such long time estimates are established by some contracting properties, 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 (the $\nabla V$ terms).

4 Numerical Experiments

We consider some test examples in [7] 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 datasets 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-1650 v2 @ 3.50GHz and 64GB memory.

4.1 1D Gaussian Mixture

As a first example, we use the Gaussian mixture probability in [4] for RBM-SVGD. The initial distribution is $\mathcal{N}(-10,1)$, Gaussian with mean $-10$ and variance 1. The target density is given by the following Gaussian mixture

[TABLE]

The kernel for the RKHS is the following Gaussian kernel

[TABLE]

where $h$ is the bandwidth parameter. For a fair comparison with the numerical results in [4], we first reproduce their results using $N=100$ particles and dynamic bandwidth parameter $h=\frac{\mathrm{med}^{2}}{2\log N}$, where $\mathrm{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 stepsize being 0.2 and the following stepsizes are 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 $\mathcal{N}(-10,1)$ distribution.

As stated in [4], the difficulty lies in the strong disagreement between the initial density function and the target density $\pi(x)$. According to the first and second row in Figure 1, SVGD with and without fixed bandwidth parameter capture the target density efficiently and the corresponding convergent 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 500th iteration costs the same as 50 iterations of SVGD. According to Figure 1, RBM-SVGD(2) at 500th iteration significantly outperform the 50th iteration of SVGD. As we increase the batch size, as the last two rows of Figure 1, more stable and similar behavior as 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=0.35$ whereas SVGD use 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 $\bar{h}=\frac{1}{N}\sum_{i=1}^{N}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:

[TABLE]

where $\mathbb{E}_{X\sim\pi}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 $\frac{2}{3}$, $5$, and $\frac{\cos 4}{e^{2}}$. The reported runtime is also averaged over $100$ random initializations.

Figure 2 (a), (b), and (c) show 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 through 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. However, in all three figures, we observe that the MSE decays first as $p$ increases and then increases for $p\geq 64$. Such a behavior is due to the choice of bandwidth parameter. Table 1 shows the averaged runtime of RBM-SVGD and SVGD for different batch sizes. RBM-SVGD is faster than SVGD for all choices of batch sizes. Ideally, RBM-SVGD with $p=2$ should be 128 times faster than SVGD, which turns out to be $13.3$ times speedup in runtime. This is due to the nature of Matlab, since Matlab is better optimized for block matrix operations. We expect that if the code is implemented with other programming languages, e.g., C++, Fortran, etc., close-to-optimal speedup should be observed.

4.2 Bayesian Logistic Regression

In this experiment, we apply RBM-SVGD to conduct Bayesian logistic regression for binary classification for the Covertype dataset with 581012 data points and 54 features [21]. Under the same setting as Gershman [21, 4], the regression weights $w$ are assigned with a Gaussian prior $p_{0}(\omega|\alpha)=\mathcal{N}(w,\alpha^{-1})$, and the variance satisfies $p_{0}(\alpha)=\Gamma(\alpha,1,0.01)$, where $\Gamma$ represents the density of Gamma distribution. The inference is applied on posterior $p(x|D)$ with $x=[w,\log\alpha]$. The kernel $K(\cdot)$ is taken again to be the same Gaussian kernel as (4.2).

Since the problem is in high dimension, we adopt $N=512$ particles $N=512$ in this experiment, which also create more space for the selection of batch sizes. The training is done on 80% of the dataset, and the other 20% is used as the test dataset. For particle system (1.1), the computation of $-\nabla V=\nabla\log p(x)$ is expensive. Hence, we use the same strategy as mentioned in [4, section 3.2], i.e. using data-mini-batch 111To avoid confusion with our batch of particles, we call it data-mini-batch instead. of the data to form a stochastic approximation of $p(x)$ with the data-mini-batch size being $100$. Since $\nabla\log p$ depends only on $x$ as in Algorithm 3, at each time step, we call this function only once and compute $\nabla\log p$ for all particles, which means the same mini-batches are used for $\nabla\log p$ of all particles. In this experiment, we use fixed bandwidth $h=256$ for RBM-SVGD and dynamic bandwidth strategy for SVGD. The RBM-SVGD uses Algorithm 3 with initial stepsize being 0.05 and the following stepsizes are generated from AdaGrad. Different batch sizes are tested to demonstrate the efficiency of RBM-SVGD. Each configuration is executed on $50$ random initializations. The averaged test accuracies for different batch sizes are illustrated in Figure 3.

As shown in Figure 3, 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. RBM-SVGD with batch size greater than $4$, we observe similar convergent behavior as that of SVGD. The runtime of RBM-SVGD, as shown in Table 2, is faster than SVGD, where the runtime of $6000$ iterations is reported. Comparing to the similar runtime table for 1D Gaussian mixture example, as Table 1, the acceleration of RBM-SVGD is not as significant as before. This is due to the linear but expensive evaluation of $\nabla\log p$, where RBM-SVGD and SVGD spend the same amount time in the evaluation each iteration. Although the evaluation of $\nabla\log p$ is expensive, it is linear in $N$. As $N$ increases, the advantage of RBM-SVGD would be more significant. In Table 3, we list the mean and standard deviation of RBM-SVGD with $p=2$, $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.

5 Conclusion

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.

Acknowledgement

This work is supported by KI-Net NSF RNMS11-07444. The work of L. Li was partially sponsored by Shanghai Sailing Program 19YF1421300, the work of Y. Li was partially supported by OAC-1450280, the work of J.-G. Liu was partially supported by NSF DMS-1812573, and the work of J. Lu was supported in part by NSF DMS-1454939.

Appendix A Proof of Lemma 1

Proof of Lemma 1.

The proof is pretty like the one in [9]. We use the random variable $I(i,j)$ to indicate whether $i$ and $j$ are in a common batch. In particular, $I(i,j)=1$ if $i$ and $j$ are in a common batch while $I(i,j)=0$ if otherwise. Then, it is not hard to compute (see [9])

[TABLE]

We note

[TABLE]

The first equation in (A.1) clearly implies that $\mathbb{E}\chi_{i}(x)=0$. Using (A.1), we can compute directly that

[TABLE]

Rearranging this, we get the claimed expression. ∎