# Provably Fast Finite Particle Variants of SVGD via Virtual Particle
Stochastic Approximation

Aniket Das
Google Research
Bangalore, India
<EMAIL_ADDRESS>
&Dheeraj Nagaraj
Google Research
Bangalore, India
<EMAIL_ADDRESS>

###### Abstract

Stein Variational Gradient Descent (SVGD) is a popular variational inference
algorithm which simulates an interacting particle system to approximately
sample from a target distribution, with impressive empirical performance
across various domains. Theoretically, its population (i.e, infinite-particle)
limit dynamics is well studied but the behavior of SVGD in the finite-particle
regime is much less understood. In this work, we design two computationally
efficient variants of SVGD, namely VP-SVGD (which is conceptually elegant) and
GB-SVGD (which is empirically effective), both of which exhibit provably fast
finite-particle convergence rates. By introducing the concept of _virtual
particles_ , we develop novel stochastic approximations of population-limit
SVGD dynamics in the space of probability measures, which are exactly
implementable using only a finite number of particles. Our algorithms can be
viewed as specific random-batch approximations of SVGD, which are
computationally more efficient than ordinary SVGD. We show that the $n$
particles output by VP-SVGD and GB-SVGD, run for $T$ steps with batch-size
$K$, are at-least as good as i.i.d samples from a distribution whose Kernel
Stein Discrepancy to the target is at most
$O(\nicefrac{{d^{1/3}}}{{(KT)^{1/6}}})$ under standard assumptions. Our
results also hold under a mild growth condition on the potential function,
which is much weaker than the isoperimetric (e.g. Poincare Inequality) or
information-transport conditions (e.g. Talagrand’s Inequality
$\mathsf{T}_{1}$) generally considered in prior works. As a corollary, we
consider the convergence of the empirical measure (of the particles output by
VP-SVGD and GB-SVGD) to the target distribution and demonstrate a _double
exponential improvement_ over the best known finite-particle analysis of SVGD.
Beyond this, our results present the _first known oracle complexities for this
setting with polynomial dimension dependence_ , thereby completely eliminating
the curse of dimensionality exhibited by previously known finite-particle
rates.

##  1 Introduction

Sampling from a distribution over $\mathbb{R}^{d}$ whose density
$\pi^{\star}(\mathbf{x})\propto\exp(-F(\mathbf{x}))$ is known only upto a
normalizing constant, is a fundamental problem in machine learning [51, 19],
statistics [40, 35], theoretical computer science [27, 16] and statistical
physics [38, 14]. A popular approach to this is the Stein Variational Gradient
Descent (SVGD) algorithm introduced by Liu and Wang [30], which uses a
positive definite kernel $k$ to evolve a set of $n$ interacting particles
$(\mathbf{x}^{(i)}_{t})_{i\in[n],t\in\mathbb{N}}$ as follows:

$\vspace{-0.05cm}\mathbf{x}^{(i)}_{t+1}\leftarrow\mathbf{x}^{(i)}_{t}-\frac{\gamma}{n}\sum_{j=1}^{n}\left[k(\mathbf{x}^{(i)}_{t},\mathbf{x}^{(j)}_{t})\nabla
F(\mathbf{x}^{(j)}_{t})-\nabla_{2}k(\mathbf{x}^{(i)}_{t},\mathbf{x}^{(j)}_{t})\right]$
(1)

SVGD exhibits remarkable empirical performance in a variety of Bayesian
inference, generative modeling and reinforcement learning tasks [30, 50, 53,
21, 49, 33, 44, 18] and usually converges rapidly to the target density while
using only a few particles, often outperforming Markov Chain Monte Carlo
(MCMC) methods. However, in contrast to its wide practical applicability,
theoretical analysis of the behavior SVGD is a relatively unexplored problem.
Prior works on the analysis of SVGD [25, 13, 29, 41, 6] mainly consider the
population limit, where the number of particles $n\to\infty$. These works
assume that the initial distribution of the (infinite number of) particles has
a finite KL divergence to the target $\pi^{\star}$ and subsequently, interpret
the dynamics of population-limit SVGD as projected gradient descent updates
for the KL divergence on the space of probability measures, equipped with the
Wasserstein geometry. Under appropriate assumptions on the target density, one
can then use the theory of Wasserstein Gradient Flows to establish non-
asymptotic (in time) convergence of population-limit SVGD to $\pi^{\star}$ in
the Kernel Stein Discrepancy (KSD) metric.

While the framework of Wasserstein Gradient Flows suffices to explain the
behavior of SVGD in the population limit, the same techniques are insufficient
to effectively analyze SVGD in the finite-particle regime. This is primarily
due to the fact that the empirical measure $\hat{\mu}^{(n)}$ of a finite
number of particles does not admit a density with respect to the Lebesgue
measure, and thus, its KL divergence to the target is always infinite (i.e.
$\mathsf{KL}\left(\hat{\mu}^{(n)}\boldsymbol{||}\pi^{\star}\right)=\infty$).
In such a setting, a direct analysis of the dynamics of finite-particle SVGD
becomes prohibitively difficult due to complex inter-particle dependencies. To
the best of our knowledge, the pioneering work of Shi and Mackey [42] is the
only result that obtains an explicit convergence rate of finite-particle SVGD
by tracking the deviation between the law of $n$-particle SVGD and that of its
population-limit. To this end, Shi and Mackey [42] show that for subgaussian
target densities, the empirical measure of $n$-particle SVGD converges to
$\pi^{\star}$ at a rate of $O(\sqrt{\tfrac{\mathsf{poly}(d)}{\log\log
n^{\Theta(1/d)}}})$ in KSD 111We explicate the dimension dependence in Shi and
Mackey [42] by closely following their analysis. The obtained convergence rate
is quite slow and fails to adequately explain the impressive practical
performance of SVGD.

Our work takes a starkly different approach to this problem and deliberately
deviates from tracking population-limit SVGD using a finite number of
particles. Instead, we directly analyze the dynamics of KL divergence along a
carefully constructed trajectory in the space of distributions. To this end,
our proposed algorithm, Virtual Particle SVGD (VP-SVGD) devises an _unbiased
stochastic approximation (in the space of measures) of the population-limit
dynamics of SVGD_. We achieve this by considering additional particles called
_virtual particles_ 222(roughly) analogous to virtual particles in quantum
field theory that enable interactions between real particles which evolve in
time but aren’t part of the output (i.e. _real particles_). These virtual
particles are used only to compute information about the current population-
level distribution of the real particles, and enable exact implementation of
our stochastic approximation to population-limit SVGD, while using only a
finite number of particles.

Our analysis is similar in spirit to non-asymptotic analyses of stochastic
gradient descent (SGD) that generally do not attempt to track gradient descent
(analogous to population-limit SVGD in this case), but instead directly track
the evolution of the objective function along the SGD trajectory using
appropriate stochastic descent lemmas [23, 20, 10]. The key feature of our
proposed stochastic approximation is the fact that it can be implemented using
only a finite number of particles. This allows us to design faster variants of
SVGD with provably fast finite-particle convergence.

###  1.1 Contributions

Result | Algorithm | Assumption | Rate |  Oracle Complexity
---|---|---|---|---
Korba et al. [25] |  Population Limit SVGD |  Uniformly Bounded $\mathsf{KSD}_{\pi^{\star}}(\bar{\mu}_{t}||\pi^{\star})$ | $\tfrac{\mathsf{poly}(d)}{\sqrt{T}}$ |  Not Implementable
Salim et al. [41] |  Population Limit SVGD | Sub-gaussian $\pi^{\star}$ | $\tfrac{d^{\nicefrac{{3}}{{2}}}}{\sqrt{T}}$ |  Not Implementable
Shi and Mackey [42] | SVGD | Sub-gaussian $\pi^{\star}$ | $\tfrac{\mathsf{poly}(d)}{\sqrt{\log\log n^{\Theta(\nicefrac{{1}}{{d}})}}}$ | $\tfrac{\mathsf{poly}(d)}{\epsilon^{2}}e^{{\Theta(de^{\nicefrac{{\mathsf{poly}(d)}}{{\epsilon^{2}}}})}}$
Ours, Corollary 1 | VP-SVGD |  Sub-gaussian $\pi^{\star}$ | $(\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{4}}}+(\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{2}}}$ | $\nicefrac{{d^{4}}}{{\epsilon^{12}}}$
Ours, Corollary 2 | GB-SVGD |  Sub-gaussian $\pi^{\star}$ | $\nicefrac{{d^{\nicefrac{{1}}{{3}}}}}{{n^{\nicefrac{{1}}{{12}}}}}+(\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{2}}}$ | $\nicefrac{{d^{6}}}{{\epsilon^{18}}}$
Ours, Corollary 1 | VP-SVGD |  Sub-exponential $\pi^{\star}$ | $\tfrac{d^{\nicefrac{{1}}{{3}}}}{n^{\nicefrac{{1}}{{6}}}}+\tfrac{d}{n^{\nicefrac{{1}}{{2}}}}$ | $\nicefrac{{d^{6}}}{{\epsilon^{16}}}$
Ours, Corollary 2 | GB-SVGD |  Sub-exponential $\pi^{\star}$ | $\tfrac{d^{\nicefrac{{3}}{{8}}}}{n^{\nicefrac{{1}}{{16}}}}+\tfrac{d}{n^{\nicefrac{{1}}{{2}}}}$ | $\nicefrac{{d^{9}}}{{\epsilon^{24}}}$

Table 1: Comparison of our results with prior works. $d$, $T$, and $n$ denote
the dimension, no. of iterations and no. of output particles respectively.
Oracle Complexity denotes number of evaluations of $\nabla F$ needed to
achieve $\mathsf{KSD}_{\pi^{\star}}(\cdot||\pi^{\star})\leq\epsilon$ and Rate
denotes convergence rate w.r.t KSD metric. Note that: 1. Population Limit SVGD
is not implementable as it requires infinite particles 2. The uniformly
bounded $\mathsf{KSD}_{\pi^{\star}}(\bar{\mu}_{t}||\pi^{\star})$ assumption
cannot be verified apriori and is much stronger than subgaussianity (see [41]
Lemma C.1)

VP-SVGD and GB-SVGD We propose two variants of SVGD that enjoy provably fast
finite-particle convergence guarantees to the target distribution: Virtual
Particle SVGD (VP-SVGD in Algorithm 1) and Global Batch SVGD (GB-SVGD in
Algorithm 2). VP-SVGD is a conceptually elegant stochastic approximation (in
the space of probability measures) of population-limit SVGD, and GB-SVGD is a
practically efficient version of SVGD which achieves good empirical
performance. Our analysis of GB-SVGD builds upon that of VP-SVGD. When the
potential $F$ is smooth and satisfies a quadratic growth condition (which
holds under subgaussianity of $\pi^{\star}$, a common assumption in prior
works [41, 42]), we show that the $n$ particles output by $T$ steps of our
algorithms, run with batch-size $K$, are at least as good as i.i.d draws from
a distribution whose Kernel Stein Discrepancy to $\pi^{\star}$ is at most
$O(\nicefrac{{d^{1/3}}}{{(KT)^{1/6}}})$. Our results also hold under a mild
subquadratic growth condition for $F$, which is much weaker than isoperimetric
(e.g. Poincare Inequality) or information-transport (e.g. Talagrand’s
Inequality $\mathsf{T}_{1}$) assumptions generally considered in the sampling
literature [46, 41, 42, 7, 2].

State-of-the-art Finite Particle Guarantees As corollaries of the above
result, we establish that _VP-SVGD and GB-SVGD exhibit the best known finite-
particle guarantees in the literature which significantly outperform that of
prior works_. Our results are summarized in Table 1. In particular, under
subgaussianity of the target distribution $\pi^{\star}$, we show that the
empirical measure of the $n$ particles output by VP-SVGD converges to
$\pi^{\star}$ in KSD at a rate of
$O((\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{4}}}+(\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{2}}})$.
Similarly, the empirical measure of the $n$ output particles of GB-SVGD
converges to $\pi^{\star}$ at a KSD rate of
$O(\nicefrac{{d^{\nicefrac{{1}}{{3}}}}}{{n^{\nicefrac{{1}}{{12}}}}}+(\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{2}}})$.
Both these results represent a double exponential improvement over the
$O(\tfrac{\mathsf{poly}(d)}{\sqrt{\log\log n^{\Theta(\nicefrac{{1}}{{d}})}}})$
KSD rate of $n$-particle SVGD obtained by Shi and Mackey [42], which, to our
knowledge, is the best known prior result for SVGD in the finite particle
regime. When benchmarked in terms of gradient oracle complexity, i.e., the
number of evaluations of $\nabla F$ required by an algorithm to achieve
$\mathsf{KSD}_{\pi^{\star}}(\cdot||\pi^{\star})\leq\epsilon$, we demonstrate
that for subgaussian $\pi^{\star}$, the oracle complexity of VP-SVGD is
$O(\nicefrac{{d^{4}}}{{\epsilon^{12}}})$ while that of GB-SVGD is
$O(\nicefrac{{d^{6}}}{{\epsilon^{18}}})$. To the best of our knowledge, our
result presents the _first known oracle complexity guarantee with polynomial
dimension dependence_ , and consequently, does not suffer from a curse of
dimensionality unlike prior works. Furthermore, as discused above, the
conditions under which our result holds is far weaker than subgaussianity of
$\pi^{\star}$, and as such, includes sub-exponential targets and beyond. In
particular, _our guarantees for sub-exponential target distributions are (to
the best of our knowledge) the first of its kind._

Computational Benefits: VP-SVGD and GB-SVGD can be viewed as specific random
batch approximations of SVGD. Our experiments (Section 8) show that GB-SVGD
obtains similar performance as SVGD but requires fewer computations. In this
context, a different kind of random batch method that divides the particles
into random subsets of interacting particles, has been proposed by Li et al.
[28]. However, the objective in Li et al. [28] is to approximate finite-
particle SVGD dynamics using the random batch method, instead of analyzing
convergence of the random batch method itself. As such, their guarantees also
suffer from an exponential dependence on the time $T$. As explained below,
their approach is also conceptually different from our method since we use the
same random batch to evolve every particle, allowing us to interpret this as a
stochastic approximation in the space of distributions instead of in the path
space.

###  1.2 Technical Challenges

We resolve the following important conceptual challenges, which may be of
independent interest.

Stochastic Approximation in the Space of Probability Measures Stochastic
approximations are widely used in optimization, control and sampling [26, 51,
22]. In the context of sampling, stochastic approximations are generally
implemented in path space, e.g., Stochastic Gradient Langevin Dynamics (SGLD)
[51] takes a random batch approximation of the drift term via the update
$\mathbf{x}_{t+1}=\mathbf{x}_{t}-\tfrac{\eta}{K}\sum_{j=0}^{K-1}\nabla
f(\mathbf{x}_{t},\xi_{j})+\sqrt{2\eta}\epsilon_{t},\
\epsilon_{t}\sim\mathcal{N}(0,\mathbf{I})$ where
$\mathbb{E}[f(\mathbf{x}_{t},\xi_{j})|\mathbf{x}_{t}]=F(\mathbf{x}_{t})$. Such
stochastic approximations are then analyzed using the theory of stochastic
processes over $\mathbb{R}^{d}$ [11, 39, 54, 24]. However, when viewed in the
space of probability measures (i.e, $\mu_{t}=\mathrm{Law}(\mathbf{x}_{t})$),
the time-evolution of these algorithms is deterministic. In contrast, our
approach designs _stochastic approximations in the space of probability
measures_. In particular, the time-evolution of the law of any particle in VP-
SVGD and GB-SVGD are a stochastic approximation of the dynamics of population-
limit SVGD. We ensure that requires only a finite number of particles for
exact implementation.

Tracking KL Divergence in the Finite-Particle Regime The population limit
($n\to\infty$) ensures that the initial empirical distribution ($\mu_{0}$) of
SVGD admits a density (w.r.t the Lebesgue measure). When the KL divergence
between $\mu_{0}$ and $\pi^{\star}$ is finite, prior works on population-limit
SVGD analyze the time-evolution of the KL divergence to $\pi^{\star}$.
However, this approach cannot be directly used to analyze finite-particle SVGD
since the empirical distribution of a finite number of particles does not
admit a density, and thus its KL divergence to $\pi^{\star}$ is infinite. Our
analysis of VP-SVGD and GB-SVGD circumvents this obstacle by considering the
dynamics of an infinite number of particles, whose empirical measure then
admits a density. However, the careful design ensures that the dynamics of $n$
of these particles can be computed exactly, using only a finite total number
of (real + virtual) particles. When conditioned on the virtual particles,
these particles are i.i.d. and their conditional law is close to the target
distribution with high probability.

##  2 Notation and Problem Setup

We use $\|\cdot\|,\langle\cdot,\cdot\rangle$ to denote the Euclidean norm and
inner product over $\mathbb{R}^{d}$ respectively. All other norms and inner
products are subscripted to indicate their underlying space.
$\mathcal{P}_{2}(\mathbb{R}^{d})$ denotes the space of probability measures on
$\mathbb{R}^{d}$ with finite second moment, with the Wasserstein-2 metric
denoted as $\mathcal{W}_{2}\left(\mu,\nu\right)$ for
$\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$. For any two probability measures
$\mu,\nu$, we denote their KL divergence as
$\mathsf{KL}\left(\mu\boldsymbol{||}\nu\right)$. For any function $f:X\to Y$
and any probability measure $\mu$ over $X$, we let $f_{\\#}\mu$ denote the law
of $f(\mathbf{x}):\mathbf{x}\sim\mu$. Given a sigma algebra $\mathcal{F}$ over
some space $\Omega$, and a measurable space $X$, $\mu(\cdot\
;\cdot):\mathcal{F}\times\mathcal{X}\to\mathbb{R}^{+}$ is a probability kernel
if for every $x\in\mathcal{X}$, $\mu(\cdot\ ;x)$ is a measure over
$\mathcal{F}$ and for every $A\in\mathcal{F}$, the map $x\to\mu(A\ ;x)$ is
measurable. We make use of probability measures $\mu(\cdot\ ;\mathbf{x})$
where $\mathbf{x}$ is a random element of some appropriate space
$\mathcal{X}$, resulting in random probability measures. We use $[m]$ and
$(m)$ to denote the sets $\\{1,\dots,m\\}$ and $\\{0,\dots,m-1\\}$
respectively. For any finite set $A$, $\mathbb{S}_{A}$ denotes the group of
all permutations of $A$. We use the $O$ notation to characterize the
dependence of our rates on the number of iterations $T$, dimension $d$ and
batch-size $K$, suppressing numerical and problem-dependent constants. We use
$\lesssim$ to denote $\leq$ upto universal constants.

We fix a symmetric positive definite reproducing kernel
$k:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ and let the corresponding
reproducing kernel Hilbert space (RKHS) [43] be denoted as $\mathcal{H}_{0}$.
We denote the product RKHS as $\mathcal{H}=\prod_{i=1}^{d}\mathcal{H}_{0}$,
equipped with the standard inner product for product spaces. We assume $k$ is
differentiable in both its arguments and let
$\nabla_{2}k(\mathbf{x},\mathbf{y})$ denote the gradient of $k(\cdot,\cdot)$
with respect to the second argument. For any
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$, we assume $\mathcal{H}\subset
L^{2}(\mu)$ and the inclusion map $i_{\mu}:\mathcal{H}\to L^{2}(\mu)$ is
continuous. We use $P_{\mu}:L^{2}(\mu)\to\mathcal{H}$ to denote the adjoint of
$i_{\mu}$, i.e., the unique operator which satisfies $\left\langle
f,i_{\mu}g\right\rangle_{L^{2}(\mu)}=\left\langle
P_{\mu}f,g\right\rangle_{\mathcal{H}}$ for any $f\in
L^{2}(\mu),g\in\mathcal{H}$. Carmeli et al. [5] shows that $P_{\mu}$ can be
expressed as a kernel convolution, i.e., $(P_{\mu}f)(\mathbf{x})=\int
k(\mathbf{x},\mathbf{y})f(\mathbf{y})\mathsf{d}\mu(\mathbf{y})$

We define the function $h:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ as
$h(\mathbf{x},\mathbf{y})=k(\mathbf{x},\mathbf{y})\nabla
F(\mathbf{y})-\nabla_{2}k(\mathbf{x},\mathbf{y})$ and $h_{\mu}\in\mathcal{H}$
as
$h_{\mu}=P_{\mu}(\nabla_{\mathbf{x}}\log(\tfrac{\mathsf{d}\mu}{\mathsf{d}\pi^{\star}}(\mathbf{x})))$
for any $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$. Integration by parts shows
that $h_{\mu}(\mathbf{x})=\int
h(\mathbf{x},\mathbf{y})\mathsf{d}\mu(\mathbf{y})$. The convergence metric we
use is the Kernel Stein Discrepancy (KSD) metric, which is widely used for
comparing probability distributions [31, 8] and analyzing SVGD [41, 25].

######  Definition 1 (Kernel Stein Discrepancy).

Define the Langevin Stein Operator of $\pi^{\star}$ acting on any
differentiable $g:\mathbb{R}^{d}\to\mathbb{R}^{d}$:

$(T_{\pi^{\star}}g)(\mathbf{x})=\nabla\cdot g(\mathbf{x})-\left\langle\nabla
F(\mathbf{x}),g(\mathbf{x})\right\rangle$

For any two probability measures $\mu,\nu$, the Kernel Stein Discrepancy
between $\mu$ and $\nu$ (with respect to $\pi^{\star}$), denoted as
$\mathsf{KSD}_{\pi^{\star}}(\mu||\nu)$ is defined as

$\mathsf{KSD}_{\pi^{\star}}(\mu||\nu)=\sup_{\|g\|_{\mathcal{H}}\leq
1}\mathbb{E}_{\mu}[T_{\pi^{\star}}g]-\mathbb{E}_{\nu}[T_{\pi^{\star}}g]$

Using integration by parts (see Chwialkowski et al. [8]), it follows that
$\mathsf{KSD}_{\pi^{\star}}(\mu||\nu)=\|h_{\mu}-h_{\nu}\|_{\mathcal{H}}$

##### Organization:

We review the population-limit SVGD in Section 3 and derive VP-SVGD and GB-
SVGD in Section 4. We state our technical assumptions in Section 5, main
results in Section 6 and provide a proof sketch in Section 7. We present
empirical evaluation in Section 8.

##  3 Background on Population-Limit SVGD

We briefly introduce the analysis of population-limit SVGD using the theory of
Wasserstein Gradient Flows and refer the readers to Korba et al. [25] and
Salim et al. [41] for a detailed treatment.

The space $\mathcal{P}_{2}(\mathbb{R}^{d})$ equipped with the 2-Wasserstein
metric $\mathcal{W}_{2}$ is known as the Wasserstein space, which admits the
following Riemannian structure : For any
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$, the tangent space
$T_{\mu}\mathcal{P}_{2}(\mathbb{R}^{d})$ can be identified with the Hilbert
space $L^{2}(\mu)$. We can then define differentiable functionals
$\mathcal{L}:\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}$ and compute their
Wasserstein gradients, denoted as $\nabla_{\mathcal{W}_{2}}\mathcal{L}$. Note
that the target $\pi^{\star}$ is the unique minimizer over
$\mathcal{P}_{2}(\mathbb{R}^{d})$ for the functional
$\mathcal{L}[\mu]=\mathsf{KL}\left(\mu\boldsymbol{||}\pi^{\star}\right)$. The
Wasserstein Gradient of $\mathcal{L}[\mu]$ is
$\nabla_{\mathcal{W}_{2}}\mathcal{L}[\mu]=\nabla_{\mathbf{x}}\log(\tfrac{\mathsf{d}\mu}{\mathsf{d}\pi^{\star}}(\mathbf{x}))$
[1]. This powerful machinery has served as a backbone for the analysis of
algorithms such as LMC [52, 3, 2] and population-limit SVGD [13, 25, 41, 44,
6].

The updates of population-limit SVGD can be viewed as Projected Gradient
Descent in the Wasserstein space. Recall from Section 2 that the function
$h_{\mu}(\mathbf{x})=P_{\mu}(\nabla\log(\tfrac{\mathsf{d}\mu}{\mathsf{d}\pi^{\star}}))(\mathbf{x})=\int
h(\mathbf{x},\mathbf{y})\mathsf{d}\mu(\mathbf{y})$. Let $\hat{\mu}^{n}_{t}$
denote the empirical measures of the SVGD particles
$(\mathbf{x}^{(i)}_{t})_{i\in[n]}$ at timestep $t$. We note that the SVGD
updates in (1) can be recast as $\hat{\mu}^{n}_{t+1}=(I-\gamma
h_{\hat{\mu}^{n}_{t}})_{\\#}\hat{\mu}^{t}_{n}$. In the limit of infinite
particles $n\to\infty$, suppose the empirical measure $\hat{\mu}^{n}_{t}$
converges to the population measure $\bar{\mu}_{t}$. In this population limit,
the updates of SVGD can be expressed as,

$\displaystyle\bar{\mu}_{t+1}=\left(I-h_{\bar{\mu}_{t}}\right)_{\\#}\bar{\mu}_{t}=\left(I-\gamma
P_{\bar{\mu}_{t}}\left(\nabla\log(\tfrac{d\bar{\mu}_{t}}{d\pi^{\star}})\right)\right)_{\\#}\bar{\mu}_{t}=\left(I-\gamma
P_{\bar{\mu}_{t}}\left(\nabla_{\mathcal{W}_{2}}\mathsf{KL}\left(\bar{\mu}_{t}\boldsymbol{||}\pi^{\star}\right)\right)\right)_{\\#}\bar{\mu}_{t}$

Recall from Section 2 that
$P_{\bar{\mu}_{t}}:L^{2}(\bar{\mu}_{t})\to\mathcal{H}$ is the Hilbert adjoint
of $i_{\bar{\mu}_{t}}$. Since $\mathcal{H}\subset L^{2}(\bar{\mu}_{t})$, the
updates of SVGD in the population limit can be seen as Projected Wasserstein
Gradient Descent for
$\mathcal{L}[\mu]=\mathsf{KL}\left(\mu\boldsymbol{||}\pi^{\star}\right)$, with
the Wasserstein Gradient at each step being projected onto the RKHS
$\mathcal{H}$. Assuming
$\mathsf{KL}\left(\bar{\mu}_{0}\boldsymbol{||}\pi^{\star}\right)<\infty$,
convergence of population limit SVGD is then established by tracking the
evolution of $\mathsf{KL}\left(\bar{\mu}_{t}\boldsymbol{||}\pi^{\star}\right)$
under appropriate structural assumptions (such as subgaussianity) on
$\pi^{\star}$.

##  4 Algorithm and Intuition

In this section, we derive VP-SVGD (Algorithm 1), and build upon it to obtain
GB-SVGD. Consider a countably infinite collection of particles
$\mathbf{x}_{0}^{(l)}\in\mathbb{R}^{d},\ l\in\mathbb{N}\cup\\{0\\}$, sampled
i.i.d from a measure $\mu_{0}$, having a density w.r.t. the Lebesgue measure.
By the strong law of large numbers, the empirical measure of
$\mathbf{x}_{0}^{(l)}$ is almost surely equal to $\mu_{0}$ (see Dudley [12,
Theorem 11.4.1]). Let batch size $K\in\mathbb{N}$ denote the batch size, and
$\mathcal{F}_{t}$ denote the filtration $\mathcal{F}_{t},\ t\geq 0$ as
$\mathcal{F}_{t}=\sigma(\\{\mathbf{x}_{0}^{(l)}\ |\ l\leq Kt-1\\}),\ \forall\
t\in\mathbb{N}$, with $\mathcal{F}_{0}$ being the trivial $\sigma$ algebra.
For ease of exposition, we discuss the case of $K=1$ in this section and
present a complete derivation for arbitrary $K\geq 1$ in Section C. Recall
from Section 3 that the updates of population-limit SVGD in
$\mathcal{P}_{2}(\mathbb{R}^{d})$ can be expressed as follows:

$\bar{\mu}_{t+1}=(I-\gamma h_{\bar{\mu}_{t}})_{\\#}\bar{\mu}_{t}$ (2)

We aim to design a stochastic approximation in
$\mathcal{P}_{2}(\mathbb{R}^{d})$ for the updates (2), such that it admits a
finite-particle realization. To this end, we propose the following dynamics in
$\mathbb{R}^{d}$

$\mathbf{x}_{t+1}^{(s)}=\mathbf{x}_{t}^{(s)}-\gamma
h(\mathbf{x}_{t}^{(s)},\mathbf{x}_{t}^{(t)}),\quad s\in\mathbb{N}\cup\\{0\\}$
(3)

Now, for each time-step $t$, we focus on the time evolution of the particles
$(\mathbf{x}^{(l)}_{t})_{l\geq t}$ (called the _lower triangular evolution_).
From (3), we observe that for any $t\in\mathbb{N}$ and $l\geq t$,
$\mathbf{x}^{(l)}_{t}$ depends only on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t-1)}_{0},\mathbf{x}^{(l)}_{0}$.
Therefore there exists a deterministic, measurable function $H_{t}$ such that:

$\displaystyle\mathbf{x}^{(l)}_{t}=H_{t}(\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t-1)}_{0},\mathbf{x}^{(l)}_{0})\,;\quad\text{for
every }l\geq t$ (4)

Since
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t-1)}_{0},\mathbf{x}^{(l)}_{0}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mu_{0}$, we conclude from (4) that
$(\mathbf{x}^{(l)}_{t})_{l\geq t}$ are i.i.d when conditioned on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t-1)}_{0}$. To this end, we define
the random measure $\mu_{t}|\mathcal{F}_{t}$ as the law of
$\mathbf{x}_{t}^{(t)}$ conditioned on $\mathcal{F}_{t}$, i.e.,
$\mu_{t}|\mathcal{F}_{t}$ is a probability kernel $\mu_{t}(\cdot\
;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t-1)}_{0})$, where
$\mu_{0}|\mathcal{F}_{0}:=\mu_{0}$. By the strong law of large numbers,
$\mu_{t}|\mathcal{F}_{t}$ is equal to the empirical measure of
$(\mathbf{x}^{(l)}_{t})_{l\geq t}$ conditioned on $\mathcal{F}_{t}$. We will
use $\mu_{t}|\mathcal{F}_{t}$ and $\mu_{t}(\cdot\
;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t-1)}_{0})$ interchangeably.

Define the random function $g_{t}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ as
$g_{t}(\mathbf{x}):=h(\mathbf{x},\mathbf{x}^{(t)}_{t})$. From (4), we note
that $g_{t}$ is $\mathcal{F}_{t+1}$ measurable. From (3), we infer that the
particles satisfy the following relation:

$\displaystyle\mathbf{x}^{(s)}_{t+1}=(I-\gamma
g_{t})(\mathbf{x}^{(s)}_{t}),\quad s\geq t+1$

Recall that
$\mathbf{x}_{t+1}^{(s)}|\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t)}_{0}\sim\mu_{t+1}|\mathcal{F}_{t+1}$
for any $s\geq t+1$. Furthermore, from Equation (4), we note that for $s\geq
t+1$, $\mathbf{x}^{(s)}_{t}$ depends only on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t-1)}_{0}$ and
$\mathbf{x}^{(s)}_{0}$. Hence, we conclude that
$\mathrm{Law}(\mathbf{x}^{(s)}_{t}|\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t)}_{0})=\mathrm{Law}(\mathbf{x}^{(s)}_{t}|\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(t-1)}_{0})=\mu_{t}|\mathcal{F}_{t}$.
With this insight, the dynamics of the lower-triangular evolution in
$\mathcal{P}_{2}(\mathbb{R}^{d})$ that the following holds almost surely:

$\mu_{t+1}|\mathcal{F}_{t+1}=(I-\gamma g_{t})_{\\#}\mu_{t}|\mathcal{F}_{t}$
(5)

$\mathbf{x}^{(t)}_{t}|\mathcal{F}_{t}\sim\mu_{t}|\mathcal{F}_{t}$ implies
$\mathbb{E}[g_{t}(\mathbf{x})|\mathcal{F}_{t}]=h_{\mu_{t}|\mathcal{F}_{t}}(\mathbf{x})$.
Thus _lower triangular dynamics_ (5) is a stochastic approximation in
$\mathcal{P}_{2}(\mathbb{R}^{d})$ to the population limit of SVGD (2). Setting
the batch size to general $K$ and tracking the evolution of the first $KT+n$
particles, we obtain VP-SVGD (Algorithm 1).

Algorithm 1 Virtual Particle SVGD (VP-SVGD)

Input: Number of steps $T$, number of output particles $n$, batch size $K$,
Initial positions
$\mathbf{x}_{0}^{(0)},\dots,\mathbf{x}_{0}^{(n+KT-1)}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mu_{0}$, Kernel $k$, step size $\gamma$.

1: for $t\in\\{0,\ldots,T-1\\}$ do

2: for $s\in\\{0,\ldots,KT+n-1\\}$ do

3:
$\mathbf{x}_{t+1}^{(s)}=\mathbf{x}_{t}^{(s)}-\frac{\gamma}{K}\sum_{l=0}^{K-1}[k(\mathbf{x}_{t}^{(s)},\mathbf{x}_{t}^{(tK+l)})\nabla
F(\mathbf{x}_{t}^{(tK+l)})-\nabla_{2}k(\mathbf{x}_{t}^{(s)},\mathbf{x}_{t}^{(tK+l)})]$

4: end for

5: end for

6: Draw $S$ uniformly at random from $\\{0,\dots,T-1\\}$

7: Output
$(\mathbf{y}^{(0)},\dots,\mathbf{y}^{(n-1)})=(\mathbf{x}_{S}^{(TK)},\dots,\mathbf{x}^{(TK+n-1)}_{S})$

Virtual Particles In Algorithm 1, $(\mathbf{x}_{t}^{(l)})_{KT\leq l\leq
KT+n-1}$ are the _real particles_ which constitute the output.
$(\mathbf{x}_{t}^{(l)})_{l<KT}$ are _virtual particles_ which propagate
information about the probability measure $\mu_{t}|\mathcal{F}_{t}$ to enable
computation of $g_{t}$, an unbiased estimate of the projected Wasserstein
gradient $h_{\mu_{t}|\mathcal{F}_{t}}$.

VP-SVGD as SVGD Without Replacement VP-SVGD is a without-replacement random-
batch approximation of SVGD (1), where a different batch is used across
timesteps, but the same batch is across particles given a fixed timestep. With
i.i.d. initialization, picking the ‘virtual particles’ in a fixed order or
from a random permutation does not change the evolution of the real particles.
With this insight, we design GB-SVGD (Algorithm 2) where we consider $n$
particles and output $n$ particles (instead of wasting $KT$ particles as
‘virtual particles’) via a random-batch approximation of SVGD.

Algorithm 2 Global Batch SVGD (GB-SVGD)

Input: $\\#$ of time steps $T$, $\\#$ of particles $n$,
$\mathbf{x}_{0}^{(0)},\dots,\mathbf{x}_{0}^{(n-1)}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mu_{0}$, Kernel $k$, step size $\gamma$, Batch size $K$,
Sampling method $\in\\{\text{with replacement},\text{without replacement}\\}$

1: for $t\in\\{0,\ldots,T-1\\}$ do

2: $\mathcal{K}_{t}\leftarrow$ random subset of $[n]$ of size $K$ (via.
sampling method)

3: for $s\in\\{0,\dots,n-1\\}$ do

4:
$\mathbf{x}_{t+1}^{(s)}=\mathbf{x}_{t}^{(s)}-\frac{\gamma}{K}\sum_{r\in\mathcal{K}_{t}}[k(\mathbf{x}_{t}^{(s)},\mathbf{x}_{t}^{(r)})\nabla
F(\mathbf{x}_{t}^{(r)})-\nabla_{2}k(\mathbf{x}_{t}^{(s)},\mathbf{x}_{t}^{(r)})]$

5: end for

6: end for

7: Draw $S$ uniformly at random from $\\{0,1,\dots,T-1\\}$

8: Output
$(\bar{\mathbf{y}}^{(0)},\dots,\bar{\mathbf{y}}^{(n-1)})=(\mathbf{x}_{S}^{(0)},\dots,\mathbf{x}^{(n-1)}_{S})$

In Algorithm 2, with replacement sampling means selecting a batch of $K$
particles i.i.d. from the uniform distribution over $[n]$. Without replacement
sampling means fixing a random permutation $\sigma$ over $\\{0,\dots,n-1\\}$
and selecting the batches in the order specified by the permutation.

##  5 Assumptions

In this section, we discuss the key assumptions required for our analysis of
VP-SVGD and GB-SVGD. Our first assumption is smoothness of $F$, which is
standard in optimization and sampling.

######  Assumption 1 (L-Smoothness).

$\nabla F$ exists and is $L$ Lipschitz. Moreover $\|\nabla
F(0)\|\leq\sqrt{L}$.

It is easy find a point such that $\|\nabla F(\mathbf{x}^{*})\|\leq\sqrt{L}$
(e.g., using $\Theta(1)$ gradient descent steps [36]) and center the
initialization at $\mu_{0}$ at $\mathbf{x}^{*}$. For clarity, we take
$\mathbf{x}^{*}=0$ without loss of generality. We now impose the following
growth condition on $F$.

######  Assumption 2 (Growth Condition).

There exist $\alpha,d_{1},d_{2}>0$ such that

$F(\mathbf{x})\geq
d_{1}\|\mathbf{x}\|^{\alpha}-d_{2}\quad\forall\mathbf{x}\in\mathbb{R}^{d}$

Note that Assumption 1 ensures $\alpha\leq 2$. Assumption 2 is essentially a
tail decay assumption on the target density $\pi^{\star}(\mathbf{x})\propto
e^{-F(\mathbf{x})}$. In fact, as we shall show in Appendix B, Assumption 1
ensures that the tails of $\pi^{\star}$ decay as $\propto
e^{-\|\mathbf{x}\|^{\alpha}}$. Consequently, Assumption 2 holds with
$\alpha=2$ when $\pi^{\star}$ is subgaussian and with $\alpha=1$ when
$\pi^{\star}$ is subexponential. Subgaussianity is equivalent to $\pi^{\star}$
satisfying the $\mathsf{T}_{1}$ inequality [48, Theorem 22.10], commonly
assumed in prior works on SVGD [41, 42]. We also note that subexponentiality
is implied when $\pi^{\star}$ satisfies the Poincare Inequality [4, Section
4], which is considered a mild condition in the sampling literature [46, 7, 2,
11, 6]. This makes Assumption 1 significantly weaker than the isoperimetric or
information-transport assumptions considered in prior works.

Next, we impose a mild assumption on the RKHS of the kernel $k$, which has
been used by several prior works [41, 25, 44, 42].

######  Assumption 3 (Bounded RKHS Norm).

For any $\mathbf{y}\in\mathbb{R}^{d}$, $k(\cdot,\mathbf{y})$ satisfies
$\|k(\cdot,\mathbf{y})\|_{\mathcal{H}_{0}}\leq B$. Furthermore,
$\nabla_{2}k(\cdot,\mathbf{y})\in\mathcal{H}$ and
$\|\nabla_{2}k(\cdot,\mathbf{y})\|_{\mathcal{H}}\leq B$

Assumption 3 ensures that the adjoint operator $P_{\mu}$, used in Sections 2
and 3, is well-defined. We also make the following assumptions on the kernel
$k$, which is satisfied by a large class of standard kernels such as Radial
Basis Function kernels and Matérn kernels of order $\geq\nicefrac{{3}}{{2}}$.

######  Assumption 4 (Kernel Decay).

The kernel $k$ satisfies the following for constants $A_{1},A_{2},A_{3}>0$.

$\displaystyle 0\leq
k(\mathbf{x},\mathbf{y})\leq\tfrac{A_{1}}{1+\|\mathbf{x}-\mathbf{y}\|^{2}},\quad\
\ \|\nabla_{2}k(\mathbf{x},\mathbf{y})\|\leq A_{2},\quad\ \
\|\nabla_{2}k(\mathbf{x},\mathbf{y})\|^{2}\leq A_{3}k(\mathbf{x},\mathbf{y})$

Finally, we make the following mild assumption on the initialization.

######  Assumption 5 (Initialization).

The initial distribution $\mu_{0}$ is such that
$\mathsf{KL}\left(\mu_{0}\boldsymbol{||}\pi^{\star}\right)<\infty$.
Furthermore, $\mu_{0}$ is supported in $\mathcal{B}(R)$, the $\ell_{2}$ ball
of radius $R$

Since prior works usually assume Gaussian initialization [41, 46], Assumption
5 may seem slightly non-standard. However, this is not a drawback. In fact,
whenever $R=\Theta(\sqrt{d}+\mathsf{polylog}(\nicefrac{{n}}{{\delta}}))$,
Gaussian initialization can be made indistinguishable from
$\mathsf{Uniform}(\mathcal{B}(R))$ initialization, with probability at least
$1-\delta$, via a coupling argument. To this end, we impose Assumption 5 for
ease of exposition and our results can be extended to consider Gaussian
initialization. In Appendix B we show that taking
$R=\sqrt{\nicefrac{{d}}{{L}}}$ and $\mu_{0}=\mathsf{Uniform}(\mathcal{B}(R))$
suffices to ensure
$\mathsf{KL}\left(\mu_{0}\boldsymbol{||}\pi^{\star}\right)=O(d)$.

##  6 Results

###  6.1 VP-SVGD

Our first result, proved in Appendix C, shows that the law of the _real
particles of_ VP-SVGD , when conditioned on the virtual particles, is close to
$\pi^{\star}$ in KSD. As a consequence, it shows that the particles output by
VP-SVGD are i.i.d. samples from a random probability measure
$\bar{\mu}(\cdot;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S)$ which
is close to $\pi^{\star}$ in KSD.

######  Theorem 1 (Convergence of VP-SVGD).

Let $\mu_{t}$ be as defined in Section 4. Let Assumptions 1 2, 3, 4, and 5 be
satisfied and let
$\gamma\leq\min\\{\nicefrac{{1}}{{2A_{1}L}},\nicefrac{{1}}{{(4+L)B}}\\}$.
There exist $(\zeta_{i})_{0\leq i\leq 3}$ depending polynomially on
$A_{1},A_{2},A_{3},B,L,d_{1},d_{2}$ for any fixed $\alpha\in(0,2]$, such that
whenever $\gamma\xi\leq\tfrac{1}{2B}$, with $\xi=\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$,
the following holds:

$\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[\mathsf{KSD}^{2}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})\right]\leq\frac{2\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{\gamma B(4+L)\xi^{2}}{K}$

Define the probability kernel $\bar{\mu}(\cdot\ ;\cdot)$ as follows: For any
$x_{\tau}\in\mathbb{R}^{d}$, $\tau\in(KT)$ and $s\in(T)$, $\bar{\mu}(\cdot\
;x_{0},\dots,x_{KT-1},s):=\mu_{s}(\cdot\ ;x_{0},\dots,x_{Ks-1})$ and
$\bar{\mu}(\cdot\ ;x_{0},\dots,x_{KT-1},s=0):=\mu_{0}(\cdot)$. Conditioned on
$\mathbf{x}^{(0)}_{\tau}=x_{\tau},\ S=s$ for every $\tau\in(KT)$, the outputs
$\mathbf{y}^{(0)},\dots,\mathbf{y}^{(n-1)}$ of VP-SVGD are i.i.d samples from
$\bar{\mu}(\cdot\ ;x_{0},\dots,x_{KT-1},s)$. Furthermore,

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\bar{\mu}(\cdot\
;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}_{0}^{(KT-1)},S)||\pi^{\star})]\leq\frac{2\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{\gamma B(4+L)\xi^{2}}{K}$

Convergence Rates Taking $\mu_{0}=\mathsf{Uniform}(\mathcal{B}(R))$ with
$R=\sqrt{\nicefrac{{d}}{{L}}}$ ensures
$\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)=O(d)$
(see Appendix B). Under this setting, choosing
$\gamma=O(\tfrac{(Kd)^{\eta}}{T^{1-\eta}})$ ensures that
$\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\bar{\mu}||\pi^{\star})]=O(\tfrac{d^{1-\eta}}{(KT)^{\eta}})$
where $\eta=\tfrac{\alpha}{2(1+\alpha)}$. Thus, for $\alpha=2$, (i.e, sub-
Gaussian $\pi^{\star}$),
$\mathsf{KSD}^{2}=O(\tfrac{d^{\nicefrac{{2}}{{3}}}}{(KT)^{\nicefrac{{1}}{{3}}}})$.
For $\alpha=1$ (i.e, sub-Exponential $\pi^{\star}$), the rate (in squared KSD)
becomes $O(\tfrac{d^{\nicefrac{{3}}{{4}}}}{(KT)^{\nicefrac{{1}}{{4}}}})$. To
the best of our knowledge, our convergence guarantee for sub-exponential
$\pi^{\star}$ is the first of its kind.

Comparison with Prior Works Salim et al. [41] analyzes population-limit SVGD
for subgaussian $\pi^{\star}$, obtaining
$\mathsf{KSD}^{2}=O(\nicefrac{{d^{\nicefrac{{3}}{{2}}}}}{{T}})$ rate. We note
that population-limit SVGD is not implementable whereas VP-SVGD is an
implementable algorithm whose outputs are samples from a distribution with
guaranteed convergence to $\pi^{\star}$.

We also establish a convergence guarantee for VP-SVGD which holds with high
probability, when conditioned on the virtual particles
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0}$. The proof of this result
is presented in Appendix C

######  Theorem 2 (VP-SVGD: High-Probability Rates).

Let the assumptions and parameter settings of Theorem 1 apply and let
$\delta\in(0,1)$. Then, the following holds with probability at least
$1-\delta$:

$\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}$
$\displaystyle\leq\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma(4+L)B\xi^{2}}{K}$
$\displaystyle+\frac{32\xi^{2}\log(\nicefrac{{2}}{{\delta}})}{KT}+12\gamma(4+L)B\xi^{2}\sqrt{\frac{\log(\nicefrac{{2}}{{\delta}})}{T}}$

Let $\bar{\mu}(\cdot;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S)$ be
the probability kernel defined in the statement of Theorem 1. Then,
conditioned on $\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S$, the $n$
particles output by VP-SVGD are i.i.d samples from
$\bar{\mu}(\cdot;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S)$.
Furthermore, with probability at least $1-\delta$

$\displaystyle\mathbb{E}_{S}[\mathsf{KSD}_{\pi^{\star}}(\bar{\mu}(\cdot;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S)||\pi^{\star})^{2}]$
$\displaystyle\leq\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma(4+L)B\xi^{2}}{K}$
$\displaystyle+\frac{32\xi^{2}\log(\nicefrac{{2}}{{\delta}})}{KT}+12\gamma(4+L)B\xi^{2}\sqrt{\frac{\log(\nicefrac{{2}}{{\delta}})}{T}}$

where $\mathbb{E}_{S}$ denotes that the expectation is being taken only with
respect to $S\sim\mathsf{Uniform}((T))$

####  6.1.1 VP-SVGD :Rapid Convergence of the Empirical Measure to the Target
Distribution

As a corollary of Theorem 1, we show that the empirical distribution of the
$n$ particles output by VP-SVGD rapidly converges to $\pi^{\star}$ in KSD. The
proof of this result is presented in Appendix E.1

######  Corollary 1 (VP-SVGD : Fast Finite-Particle Convergence).

Let the assumptions and parameter settings of Theorem 1 be satisfied. Let
$\hat{\mu}^{(n)}$ denote the empirical measure of the $n$ particles output by
VP-SVGD.

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq\frac{\xi^{2}}{n}+\frac{2\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{\gamma B(4+L)\xi^{2}}{K}$

where $\xi$ is as defined in Theorem 1. Setting
$R=\sqrt{\nicefrac{{d}}{{L}}},\gamma=O(\tfrac{(Kd)^{\eta}}{T^{1-\eta}})$ with
$\eta=\tfrac{\alpha}{2(1+\alpha)}$ and
$KT=d^{\tfrac{\alpha}{2+\alpha}}n^{\tfrac{2(1+\alpha)}{2+\alpha}}$ suffices to
ensure,

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq
O\left(\frac{d^{\tfrac{2}{2+\alpha}}}{n^{\tfrac{\alpha}{2+\alpha}}}+\frac{d^{\nicefrac{{2}}{{\alpha}}}}{n}\right)$

Comparison to Prior Work For subgaussian $\pi^{\star}$ (i.e., $\alpha=2$), VP-
SVGD exhibits a finite-particle rate of
$\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]=O((\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{4}}}+(\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{2}}})$.
This represents a _double exponential improvement_ over the SOTA
$\tilde{O}(\tfrac{\mathsf{poly}(d)}{\sqrt{\log\log
n^{\Theta(\nicefrac{{1}}{{d}})}}})$ KSD rate obtained by Shi and Mackey [42]
for SVGD for subgaussian targets. For subexponential $\pi^{\star}$ (i.e.
$\alpha=1$), the KSD rate of VP-SVGD is
$O(\tfrac{d^{\nicefrac{{1}}{{3}}}}{n^{\nicefrac{{1}}{{6}}}}+\tfrac{d}{n^{\nicefrac{{1}}{{2}}}})$.
To the best of our knowledge, this is the first known result of its kind.

Oracle Complexity As we shall show in Appendix E.3, the oracle complexity of
VP-SVGD to achieve $\epsilon$-convergence in KSD for subgaussian $\pi^{\star}$
is $O(\nicefrac{{d^{4}}}{{\epsilon^{12}}})$, which is a double exponential
improvement over the
$O(\tfrac{\mathsf{poly}(d)}{\epsilon^{2}}e^{{\Theta(de^{\nicefrac{{\mathsf{poly}(d)}}{{\epsilon^{2}}}})}})$
oracle complexity of SVGD as implied by Shi and Mackey [42]. Notably, the
obtained oracle complexity of VP-SVGD does not suffer from a curse of
dimensionality and to the best of our knowledge, is the _first known oracle
complexity for this problem with polynomial dimension dependence_. For
subexponential $\pi^{\star}$, the obtained oracle complexity is
$O(\nicefrac{{d^{6}}}{{\epsilon^{16}}})$ which is the first known result of
its kind.

###  6.2 GB-SVGD

We now use VP-SVGD as the basis to analyze GB-SVGD. Assume $n>KT$. Then, with
probability at least $1-\nicefrac{{K^{2}T^{2}}}{{n}}$ (for with-replacement
sampling) and $1$ (for without-replacement sampling), the random batches
$\mathcal{K}_{t}$ in GB-SVGD (Algorithm 2) are disjoint and contain distinct
elements. When conditioned on this event $\mathcal{E}$, we note that the
$n-KT$ particles that were not included in any random batch $\mathcal{K}_{t}$
evolve exactly like the $n$ real particles of VP-SVGD. With this insight, we
show that, conditioned on $\mathcal{E}$, the outputs of VP-SVGD and GB-SVGD
can be coupled such that the first $n-KT$ particles of both the algorithms are
exactly equal. This allows us to derive the following squared KSD bound
between the empirical measures of the outputs of VP-SVGD and GB-SVGD. The
proof of this result is presented in Appendix D

######  Theorem 3 (KSD Bounds for GB-SVGD).

Let $n>KT$ and let $\mathbf{Y}=(\mathbf{y}^{(0)},\dots,\mathbf{y}^{(n-1)})$
and $\bar{\mathbf{Y}}=(\bar{\mathbf{y}}^{(0)},\dots,\bar{\mathbf{y}}^{(n-1)})$
denote the outputs of VP-SVGD and GB-SVGD respectively. Moreover, let
$\hat{\mu}^{(n)}=\tfrac{1}{n}\sum_{i=0}^{n-1}\delta_{\mathbf{y}^{(i)}}$ and
$\hat{\nu}^{(n)}=\tfrac{1}{n}\sum_{i=0}^{n-1}\delta_{\bar{\mathbf{y}}^{(i)}}$
denote their respective empirical measures. Under the assumptions and
parameter settings of Theorem 1, there exists a coupling of $\mathbf{Y}$ and
$\bar{\mathbf{Y}}$ such that the following holds:

$\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\hat{\mu}^{(n)})]\leq\begin{cases}\frac{2K^{2}T^{2}\xi^{2}}{n^{2}}\quad&\text{(without
replacement sampling)}\\\
\frac{2K^{2}T^{2}\xi^{2}}{n^{2}}\left(1-\frac{K^{2}T^{2}}{n}\right)+\frac{2K^{2}T^{2}\xi^{2}}{n}\quad&\text{(with
replacement sampling)}\end{cases}$ (6)

####  6.2.1 GB-SVGD :Rapid Convergence of the Empirical Measure to the Target
Distribution

Theorem 3 ensures that the empirical distribution of the outputs of VP-SVGD
and GB-SVGD are close to each other in squared KSD (under an appropriate
coupling). To this end, we use Theorem 3 and the finite-particle convergence
guarantee of VP-SVGD (i.e., Corollary 1) to obtain the following finite-
particle guarantee for GB-SVGD. We present a proof of this result in Appendix
E.2

######  Corollary 2 (GB-SVGD : Fast Finite-Particle Convergence).

Let the assumptions and parameter settings of Theorem 1 be satisfied. Let
$\hat{\nu}^{(n)}$ denote the empirical measure of the $n$ particles output by
GB-SVGD. Then, under without-replacement sampling of the minibatches, the
following holds:

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]\leq\frac{4K^{2}T^{2}\xi^{2}}{n^{2}}+\frac{2\xi^{2}}{n}+\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma B(4+L)\xi^{2}}{K}$

and the following holds under with-replacement sampling of the minibatches

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]\leq\frac{4K^{2}T^{2}\xi^{2}}{n^{2}}(1-\frac{K^{2}T^{2}}{n})+\frac{4K^{2}T^{2}\xi^{2}}{n}+\frac{2\xi^{2}}{n}+\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma B(4+L)\xi^{2}}{K}$

where $\xi$ is as defined in Theorem 1. In particular, for GB-SVGD under
without-replacement sampling of the minibatches, setting
$R=\sqrt{\nicefrac{{d}}{{L}}},\gamma=O(\tfrac{(Kd)^{\eta}}{T^{1-\eta}})$ with
$\eta=\tfrac{\alpha}{2(1+\alpha)}$ and $KT=\sqrt{n}$ suffices to ensure the
following

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]\leq
O\left(\frac{d^{\nicefrac{{2}}{{\alpha}}}}{n}+\frac{d^{\tfrac{1}{1+\alpha}}}{n^{\tfrac{1+2\alpha}{2(1+\alpha)}}}+\frac{d^{\tfrac{2+\alpha}{2(1+\alpha)}}}{n^{\tfrac{\alpha}{4(1+\alpha)}}}\right)$

Comparison to Prior Work For subgaussian $\pi^{\star}$ (i.e., $\alpha=2$), GB-
SVGD exhibits a finite-particle rate of
$\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]=O(\nicefrac{{d^{\nicefrac{{1}}{{3}}}}}{{n^{\nicefrac{{1}}{{12}}}}}+(\nicefrac{{d}}{{n}})^{\nicefrac{{1}}{{2}}})$.
As before, this is a _double exponential improvement_ over the SOTA
$\tilde{O}(\tfrac{\mathsf{poly}(d)}{\sqrt{\log\log
n^{\Theta(\nicefrac{{1}}{{d}})}}})$ KSD rate obtained by Shi and Mackey [42]
for SVGD for subgaussian targets. For subexponential $\pi^{\star}$ (i.e.
$\alpha=1$), the KSD rate of GB-SVGD is
$O(\tfrac{d^{\nicefrac{{3}}{{8}}}}{n^{\nicefrac{{1}}{{16}}}}+\tfrac{d}{n^{\nicefrac{{1}}{{2}}}})$,
which, to the best of our knowledge, this is the first known result of its
kind.

Oracle Complexity As we shall show in Appendix E.3, the oracle complexity of
GB-SVGD to achieve $\epsilon$-convergence in KSD for subgaussian $\pi^{\star}$
is $O(\nicefrac{{d^{6}}}{{\epsilon^{18}}})$, which significantly improves upon
the oracle complexity of SVGD as implied by Shi and Mackey [42]. Notably, the
obtained oracle complexity of GB-SVGD does not suffer from a curse of
dimensionality and to the best of our knowledge, is the first known oracle
complexity for this problem with polynomial dimension dependence. For
subexponential $\pi^{\star}$, the obtained oracle complexity is
$O(\nicefrac{{d^{9}}}{{\epsilon^{24}}})$ which is the first known result of
its kind.

##  7 Proof Sketch

We now present a sketch of our analysis. As shown in Section 4, the particles
$(\mathbf{x}^{(l)}_{t})_{l\geq Kt}$ are i.i.d conditioned on the filtration
$\mathcal{F}_{t}$, and the random measure $\mu_{t}|\mathcal{F}_{t}$ is the law
of $(\mathbf{x}^{(Kt)}_{t})$ conditioned on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0}$. Moreover, from equation
(5), we know that $\mu_{t}|\mathcal{F}_{t}$ is a stochastic approximation of
population limit SVGD dynamics, i.e., $\mu_{t+1}|\mathcal{F}_{t+1}=(I-\gamma
g_{t})_{\\#}\mu_{t}|\mathcal{F}_{t}$. Lemma 1 (similar to Salim et al. [41,
Proposition 3.1] and Korba et al. [25, Proposition 5]) shows that under
appropriate conditions, the KL between $\mu_{t}|\mathcal{F}_{t}$ and
$\pi^{\star}$ satisfies a (stochastic) descent lemma . Hence
$\mu_{t}|\mathcal{F}_{t}$ admits a density and
$\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)$ is
almost surely finite.

######  Lemma 1 (Descent Lemma for $\mu_{t}|\mathcal{F}_{t}$).

Let Assumptions 1, 3 and 5 be satisfied and let $\beta>1$ be an arbitrary
constant. On the event $\gamma\|g_{t}\|_{\mathcal{H}}\leq\tfrac{\beta-1}{\beta
B}$, the following holds almost surely

$\displaystyle\mathsf{KL}\left(\mu_{t+1}|\mathcal{F}_{t+1}\boldsymbol{||}\pi^{\star}\right)\leq\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\gamma\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},g_{t}\right\rangle_{\mathcal{H}}+\frac{\gamma^{2}(\beta^{2}+L)B}{2}\|g_{t}\|^{2}_{\mathcal{H}}$

Lemma 1 is analogous to the noisy descent lemma which is used in the analysis
of SGD for smooth functions. Notice that
$\mathbb{E}[g_{t}|\mathcal{F}_{t}]=h_{\mu_{t}|\mathcal{F}_{t}}$ (when
interpreted as a Gelfand-Pettis integral [45], as discussed in Appendix B and
Appendix C) and hence in expectation, the KL divergence decreases in time. In
order to apply Lemma 1, we establish an almost-sure bound on
$\|g_{t}\|_{\mathcal{H}}$ below.

######  Lemma 2.

Let Assumptions 1, 2, 3, 4 and 5 hold. For
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$, the following holds almost surely,

$\displaystyle\|g_{t}\|_{\mathcal{H}}\leq\xi=\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$

where $\zeta_{0},\zeta_{1},\zeta_{2}$ and $\zeta_{3}$ which depend
polynomially on $A_{1},A_{2},A_{3},B,d_{1},d_{2},L$ for any fixed $\alpha$.

Let $K=1$ for clarity. To prove Lemma 2, we first note via smoothness of
$F(\cdot)$ and Assumption 3 that $\|g_{t}\|_{\mathcal{H}}\leq
C_{0}\|\mathbf{x}_{t}^{(t)}\|+C_{1}$, and then bound
$\|\mathbf{x}_{t}^{(t)}\|$. Now,
$g_{s}(\mathbf{x})=k(\mathbf{x},\mathbf{x}_{s}^{(s)})\nabla
F(\mathbf{x}_{s}^{(s)})-\nabla_{2}k(\mathbf{x},\mathbf{x}_{s}^{(s)})$. When
$\|\mathbf{x}_{s}^{(s)}-\mathbf{x}\|$ is large, $\|g_{s}(\mathbf{x})\|$ is
small due to decay assumptions on the kernel (Assumption 4) implying that the
particle does not move much. When $\mathbf{x}_{s}^{(s)}\approx\mathbf{x}$, we
have $g_{s}(\mathbf{x})\approx k(\mathbf{x},\mathbf{x}_{s}^{(s)})\nabla
F(\mathbf{x})-\nabla_{2}k(\mathbf{x},\mathbf{x}_{s}^{(s)})$ and
$k(\mathbf{x},\mathbf{x}_{s}^{(s)})\geq 0$. This is approximately a gradient
descent update on $F(\cdot)$ along with a bounded term
$\nabla_{2}k(\mathbf{x},\mathbf{x}_{s}^{(s)})$. Thus, the value of
$F(\mathbf{x}_{t}^{(l)})$ cannot grow too large after $T$ iterations. By
Assumption 2, $F(\mathbf{x}_{t}^{(l)})$ being small implies that
$\|\mathbf{x}_{t}^{(l)}\|$ is small.

Equipped with Lemma 2, we set the step-size $\gamma$ to ensure that the
descent lemma (Lemma 1) always holds. The remainder of the proof involves
unrolling through Lemma 1 by taking iterated expectations on both sides. To
this end we control $\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},g_{t}\right\rangle_{\mathcal{H}}$ and
$\|g_{t}\|^{2}_{\mathcal{H}}$ in expectation, in Lemma 3.

######  Lemma 3.

Let Assumptions 1,2,3,4,5 hold and $\xi$ be as defined in Lemma 2. Then, for
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$,

$\mathbb{E}\left[\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},g_{t}\right\rangle_{\mathcal{H}}|\mathcal{F}_{t}\right]=\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}\text{\quad
and
\quad}\mathbb{E}[\|g_{t}\|^{2}_{\mathcal{H}}]\leq\nicefrac{{\xi^{2}}}{{K}}+\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}$

##  8 Experiments

We compare the performance of GB-SVGD and SVGD. We take $n=100$ and use the
Laplace kernel with $h=1$ for both. We pick the stepsize $\gamma$ by a grid
search for each algorithm. Additional details are presented in Appendix G. We
observe that SVGD takes fewer iterations to converge, but the compute time for
GB-SVGD is lower. This is similar to the typical behavior of stochastic
optimization algorithms like SGD.

(a) MMD vs Compute Time

(b) MMD vs Iterations

Figure 1: Gaussian Experiment Comparing SVGD and GB-SVGD averaged over 10
experiments.

(a) Accuracy vs Compute Time

(b) Accuracy vs Iterations

Figure 2: Covertype Experiment, averaged over 50 runs. The error bars
represent 95% CI.

Sampling from Isotropic Gaussian (Figure 1): As a sanity check, we set
$\pi^{\star}=\mathcal{N}(0,\mathbf{I})$ with $d=5$. We pick $K=10$ for GB-
SVGD. The metric of convergence is MMD with respect to the empirical measure
of 1000 i.i.d. sampled Gaussians.

Bayesian Logistic Regression (Figure 2) We consider the Covertype dataset
which contains $\sim 580,000$ data points with $d=54$. We consider the same
priors suggested in Gershman et al. [15] and implemented in Liu and Wang [30].
We take $K=40$ for GB-SVGD. For both VP-SVGD and GB-SVGD, we use AdaGrad with
momentum to set the step-sizes as per Liu and Wang [30]

##  9 Conclusion

We develop two computationally efficient variants of SVGD with provably fast
convergence guarantees in the finite-particle regime, and present a wide range
of improvements over prior work. A promising avenue of future work could be to
establish convergence guarantees for SVGD with general non-logconcave targets,
as was considered in recent works on LMC and SGLD [2, 11]. Other important
avenues include establishing minimax lower bounds for SVGD and related
particle-based variational inference algorithms. Beyond this, we also
conjecture that the rates of GB-SVGD can be improved even in the regime $n\ll
KT$. However, we believe this requires new analytic tools.

## Acknowledgements

We thank Jiaxin Shi, Lester Mackey and the anonymous reviewers for their
helpful feedback. We are particularly grateful to Lester Mackey for providing
insightful pointers on the properties of Kernel Stein Discrepancy, which
greatly helped us in removing the curse of dimensionality from our finite-
particle convergence guarantees.

## References

  * Ambrosio et al. [2005] L. Ambrosio, N. Gigli, and G. Savaré.  _Gradient flows: in metric spaces and in the space of probability measures_.  Springer Science & Business Media, 2005.
  * Balasubramanian et al. [2022] K. Balasubramanian, S. Chewi, M. A. Erdogdu, A. Salim, and S. Zhang.  Towards a theory of non-log-concave sampling:first-order stationarity guarantees for langevin monte carlo.  In P.-L. Loh and M. Raginsky, editors, _Proceedings of Thirty Fifth Conference on Learning Theory_ , volume 178 of _Proceedings of Machine Learning Research_ , pages 2896–2923. PMLR, 02–05 Jul 2022.  URL https://proceedings.mlr.press/v178/balasubramanian22a.html.
  * Bernton [2018] E. Bernton.  Langevin monte carlo and jko splitting.  In _Conference On Learning Theory_ , pages 1777–1798. PMLR, 2018\.
  * Bobkov and Ledoux [1997] S. Bobkov and M. Ledoux.  Poincaré’s inequalities and talagrand’s concentration phenomenon for the exponential distribution.  _Probability Theory and Related Fields_ , 107:383–400, 1997\.
  * Carmeli et al. [2010] C. Carmeli, E. De Vito, A. Toigo, and V. Umanitá.  Vector valued reproducing kernel hilbert spaces and universality.  _Analysis and Applications_ , 8(01):19–61, 2010\.
  * Chewi et al. [2020] S. Chewi, T. Le Gouic, C. Lu, T. Maunu, and P. Rigollet.  Svgd as a kernelized wasserstein gradient flow of the chi-squared divergence.  _Advances in Neural Information Processing Systems_ , 33:2098–2109, 2020.
  * Chewi et al. [2022] S. Chewi, M. A. Erdogdu, M. Li, R. Shen, and S. Zhang.  Analysis of langevin monte carlo from poincare to log-sobolev.  In P.-L. Loh and M. Raginsky, editors, _Proceedings of Thirty Fifth Conference on Learning Theory_ , volume 178 of _Proceedings of Machine Learning Research_ , pages 1–2. PMLR, 02–05 Jul 2022.  URL https://proceedings.mlr.press/v178/chewi22a.html.
  * Chwialkowski et al. [2016] K. Chwialkowski, H. Strathmann, and A. Gretton.  A kernel test of goodness of fit.  In _International conference on machine learning_ , pages 2606–2615. PMLR, 2016.
  * Conway [2019] J. B. Conway.  _A course in functional analysis_ , volume 96.  Springer, 2019.
  * Das et al. [2022] A. Das, B. Schölkopf, and M. Muehlebach.  Sampling without replacement leads to faster rates in finite-sum minimax optimization.  _Advances in Neural Information Processing Systems_ , 2022.
  * Das et al. [2023] A. Das, D. Nagaraj, and A. Raj.  Utilising the clt structure in stochastic gradient based sampling: Improved analysis and faster algorithms.  In _Conference on Learning Theory_ , 2023.
  * Dudley [2018] R. M. Dudley.  _Real analysis and probability_.  CRC Press, 2018.
  * Duncan et al. [2019] A. Duncan, N. Nüsken, and L. Szpruch.  On the geometry of stein variational gradient descent.  _arXiv preprint arXiv:1912.00894_ , 2019.
  * El Alaoui et al. [2022] A. El Alaoui, A. Montanari, and M. Sellke.  Sampling from the sherrington-kirkpatrick gibbs measure via algorithmic stochastic localization.  In _2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 323–334. IEEE, 2022.
  * Gershman et al. [2012] S. Gershman, M. D. Hoffman, and D. M. Blei.  Nonparametric variational inference.  In _Proceedings of the 29th International Conference on International Conference on Machine Learning_ , 2012.
  * Gopi et al. [2022] S. Gopi, Y. T. Lee, and D. Liu.  Private convex optimization via exponential mechanism.  In _Conference on Learning Theory_ , pages 1948–1989. PMLR, 2022\.
  * Gorham et al. [2020] J. Gorham, A. Raj, and L. Mackey.  Stochastic stein discrepancies.  _Advances in Neural Information Processing Systems_ , 33:17931–17942, 2020.
  * Haarnoja et al. [2017] T. Haarnoja, H. Tang, P. Abbeel, and S. Levine.  Reinforcement learning with deep energy-based policies.  In _International conference on machine learning_ , pages 1352–1361. PMLR, 2017.
  * Ho et al. [2020] J. Ho, A. Jain, and P. Abbeel.  Denoising diffusion probabilistic models.  _Advances in Neural Information Processing Systems_ , 33:6840–6851, 2020.
  * Jain et al. [2021] P. Jain, D. M. Nagaraj, and P. Netrapalli.  Making the last iterate of sgd information theoretically optimal.  _SIAM Journal on Optimization_ , 31(2):1108–1130, 2021.
  * Jaini et al. [2021] P. Jaini, L. Holdijk, and M. Welling.  Learning equivariant energy based models with equivariant stein variational gradient descent.  _Advances in Neural Information Processing Systems_ , 34:16727–16737, 2021.
  * Jin et al. [2020] S. Jin, L. Li, and J.-G. Liu.  Random batch methods (rbm) for interacting particle systems.  _Journal of Computational Physics_ , 400:108877, 2020.
  * Khaled and Richtárik [2020] A. Khaled and P. Richtárik.  Better theory for sgd in the nonconvex world.  _arXiv preprint arXiv:2002.03329_ , 2020.
  * Kinoshita and Suzuki [2022] Y. Kinoshita and T. Suzuki.  Improved convergence rate of stochastic gradient langevin dynamics with variance reduction and its application to optimization.  _arXiv preprint arXiv:2203.16217_ , 2022.
  * Korba et al. [2020] A. Korba, A. Salim, M. Arbel, G. Luise, and A. Gretton.  A non-asymptotic analysis for stein variational gradient descent.  _Advances in Neural Information Processing Systems_ , 33:4672–4682, 2020.
  * Kushner and Clark [2012] H. J. Kushner and D. S. Clark.  _Stochastic approximation methods for constrained and unconstrained systems_ , volume 26.  Springer Science & Business Media, 2012.
  * Lee and Vempala [2022] Y. T. Lee and S. S. Vempala.  The manifold joys of sampling (invited talk).  In _49th International Colloquium on Automata, Languages, and Programming (ICALP 2022)_. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2022.
  * Li et al. [2020] L. Li, Y. Li, J.-G. Liu, Z. Liu, and J. Lu.  A stochastic version of stein variational gradient descent for efficient sampling.  _Communications in Applied Mathematics and Computational Science_ , 15(1):37–63, 2020.
  * Liu [2017] Q. Liu.  Stein variational gradient descent as gradient flow.  _Advances in neural information processing systems_ , 30, 2017.
  * Liu and Wang [2016] Q. Liu and D. Wang.  Stein variational gradient descent: A general purpose bayesian inference algorithm.  _Advances in neural information processing systems_ , 29, 2016.
  * Liu et al. [2016] Q. Liu, J. Lee, and M. Jordan.  A kernelized stein discrepancy for goodness-of-fit tests.  In _International conference on machine learning_ , pages 276–284. PMLR, 2016.
  * Liu et al. [2023] T. Liu, P. Ghosal, K. Balasubramanian, and N. S. Pillai.  Towards understanding the dynamics of gaussian-stein variational gradient descent, 2023.
  * Liu et al. [2017] Y. Liu, P. Ramachandran, Q. Liu, and J. Peng.  Stein variational policy gradient.  _arXiv preprint arXiv:1704.02399_ , 2017.
  * Lu et al. [2019] J. Lu, Y. Lu, and J. Nolen.  Scaling limit of the stein variational gradient descent: The mean field regime.  _SIAM Journal on Mathematical Analysis_ , 51(2):648–671, 2019.
  * Neal et al. [2011] R. M. Neal et al.  Mcmc using hamiltonian dynamics.  _Handbook of markov chain monte carlo_ , 2(11):2, 2011.
  * Nesterov [1998] Y. Nesterov.  Introductory lectures on convex programming volume i: Basic course.  _Lecture notes_ , 3(4):5, 1998.
  * Nüsken and Renger [2021] N. Nüsken and D. Renger.  Stein variational gradient descent: many-particle and long-time asymptotics.  _arXiv preprint arXiv:2102.12956_ , 2021.
  * Parisi [1981] G. Parisi.  Correlation functions and computer simulations.  _Nuclear Physics B_ , 180(3):378–384, 1981.
  * Raginsky et al. [2017] M. Raginsky, A. Rakhlin, and M. Telgarsky.  Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis.  In _Conference on Learning Theory_ , pages 1674–1703. PMLR, 2017\.
  * Roberts and Tweedie [1996] G. O. Roberts and R. L. Tweedie.  Exponential convergence of langevin distributions and their discrete approximations.  _Bernoulli_ , pages 341–363, 1996.
  * Salim et al. [2022] A. Salim, L. Sun, and P. Richtarik.  A convergence theory for svgd in the population limit under talagrand’s inequality t1.  In _International Conference on Machine Learning_ , pages 19139–19152. PMLR, 2022.
  * Shi and Mackey [2022] J. Shi and L. Mackey.  A finite-particle convergence rate for stein variational gradient descent.  _arXiv preprint arXiv:2211.09721_ , 2022.
  * Steinwart and Christmann [2008] I. Steinwart and A. Christmann.  _Support vector machines_.  Springer Science & Business Media, 2008.
  * Sun et al. [2023] L. Sun, A. Karagulyan, and P. Richtarik.  Convergence of stein variational gradient descent under a weaker smoothness condition.  In _International Conference on Artificial Intelligence and Statistics_ , pages 3693–3717. PMLR, 2023.
  * Talagrand [1984] M. Talagrand.  _Pettis integral and measure theory_.  American Mathematical Soc., 1984.
  * Vempala and Wibisono [2019] S. Vempala and A. Wibisono.  Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices.  _Advances in neural information processing systems_ , 32, 2019.
  * Vershynin [2018] R. Vershynin.  _High-dimensional probability: An introduction with applications in data science_ , volume 47.  Cambridge university press, 2018.
  * Villani et al. [2009] C. Villani et al.  _Optimal transport: old and new_ , volume 338.  Springer, 2009.
  * Wang and Liu [2016] D. Wang and Q. Liu.  Learning to draw samples: With application to amortized mle for generative adversarial learning.  _arXiv preprint arXiv:1611.01722_ , 2016.
  * Wang et al. [2018] D. Wang, Z. Zeng, and Q. Liu.  Stein variational message passing for continuous graphical models.  In _International Conference on Machine Learning_ , pages 5219–5227. PMLR, 2018.
  * Welling and Teh [2011] M. Welling and Y. W. Teh.  Bayesian learning via stochastic gradient langevin dynamics.  In _Proceedings of the 28th international conference on machine learning (ICML-11)_ , pages 681–688, 2011.
  * Wibisono [2018] A. Wibisono.  Sampling as optimization in the space of measures: The langevin dynamics as a composite optimization problem.  In _Conference on Learning Theory_ , pages 2093–3027. PMLR, 2018\.
  * Zhuo et al. [2018] J. Zhuo, C. Liu, J. Shi, J. Zhu, N. Chen, and B. Zhang.  Message passing stein variational gradient descent.  In _International Conference on Machine Learning_ , pages 6018–6027. PMLR, 2018.
  * Zou et al. [2021] D. Zou, P. Xu, and Q. Gu.  Faster convergence of stochastic gradient langevin dynamics for non-log-concave sampling.  In _Uncertainty in Artificial Intelligence_ , pages 1152–1162. PMLR, 2021.

###### Contents

  1. 1 Introduction
    1. 1.1 Contributions
    2. 1.2 Technical Challenges
  2. 2 Notation and Problem Setup
  3. 3 Background on Population-Limit SVGD
  4. 4 Algorithm and Intuition
  5. 5 Assumptions
  6. 6 Results
    1. 6.1 VP-SVGD
      1. 6.1.1 VP-SVGD :Rapid Convergence of the Empirical Measure to the Target Distribution
    2. 6.2 GB-SVGD
      1. 6.2.1 GB-SVGD :Rapid Convergence of the Empirical Measure to the Target Distribution
  7. 7 Proof Sketch
  8. 8 Experiments
  9. 9 Conclusion
  10. A Additional Notation and Organization
  11. B Preliminaries
    1. B.1 Gelfand-Pettis Integrals for Reproducing Kernel Hilbert Spaces
  12. C Analysis of VP-SVGD
    1. C.1 Population Level Dynamics : Proof of Lemma 1
    2. C.2 Iterate Bounds : Proof of Lemma 2
    3. C.3 Controlling $g_{t}$ in Expectation : Proof of Lemma 3
    4. C.4 Proof of Theorem 1
    5. C.5 Proof of Theorem 2
  13. D Analysis of GB-SVGD
    1. D.1 Proof of Theorem 3
  14. E Finite-Particle Convergence Guarantees for VP-SVGD and GB-SVGD
    1. E.1 VP-SVGD : Proof of Corollary 1
    2. E.2 GB-SVGD : Proof of Corollary 2
    3. E.3 Oracle Complexity of SVGD, VP-SVGD and GB-SVGD
      1. E.3.1 SVGD
      2. E.3.2 VP-SVGD
      3. E.3.3 GB-SVGD
  15. F Literature Review
  16. G Additional Experimental Details

##  Appendix A Additional Notation and Organization

We use $\Gamma$ to denote the Gamma function
$\Gamma(x)=\int_{0}^{\infty}t^{x-1}e^{-t}\mathsf{d}t$, and recall that for any
$n\in\mathbb{N}$, $\Gamma(n)=(n-1)!$. For any Lebesgue measurable
$A\subseteq\mathbb{R}^{d}$, we use $\mathsf{vol}(A)$ to denote it’s Lebesgue
Measure and $\mathsf{Uniform}(A)$ to denote the uniform distribution supported
on $A$. We use $\mathcal{B}(R)$ to denote the ball of radius $R$ centered at
the origin, and recall that
$\mathsf{vol}(\mathcal{B}(R))=\tfrac{\pi^{\nicefrac{{d}}{{2}}}}{\Gamma(\nicefrac{{d}}{{2}}+1)}R^{d}$.
For ease of exposition, we assume $d\geq 2$. We further assume
$\pi^{\star}(\mathbf{x})=e^{-F(\mathbf{x})}$. We note that this can be easily
ensured by absorbing the normalizing constant into $F(0)$, and does not affect
the dynamics of SVGD, VP-SVGD or GB-SVGD (since they only use the gradient
information of $F$). We highlight that both these assumptions are made purely
for the sake of clarity and are very easily removable with negligible changes
to our analysis.

In Appendix B, we discuss the technical lemmas used in our analysis, and
present a short exposition to the Gelfand-Pettis integral in Appendix B.1,
which we use to analyze VP-SVGD. We analyze VP-SVGD in Appendix C and GB-SVGD
in Appendix D. Convergence guarantees for the empirical measure of VP-SVGD and
GB-SVGD are presented in Appendix E. We give a brief review of the related
work in Section F. Additional details regarding our experimental setup are
stated in G.

##  Appendix B Preliminaries

The following lemma shows that setting the initial distribution
$\mu_{0}=\mathsf{Uniform}(\mathcal{B}(R))$ with $R=\sqrt{\nicefrac{{d}}{{L}}}$
suffices to ensure
$\mathsf{KL}\left(\mu_{0}\boldsymbol{||}\pi^{\star}\right)=O(d)$.The proof of
this result is similar to that of Vempala and Wibisono [46, Lemma 1] with the
Gaussian initialization replaced by $\mathsf{Uniform}(\mathcal{B}(R))$
initialization.

######  Lemma 4 (KL Upper Bound for Uniform Initialization).

Let Assumption 1 be satisfied and let
$\mu_{0}=\mathsf{Uniform}(\mathcal{B}(R))$ with
$R=\sqrt{\nicefrac{{d}}{{L}}}$. Then, the following holds:

$\displaystyle\mathsf{KL}\left(\mu_{0}\boldsymbol{||}\pi^{\star}\right)\leq\frac{d}{2}\log(\nicefrac{{L}}{{2\pi}})+d+F(0)+\nicefrac{{1}}{{2}}\leq
O(d)$

###### Proof.

For any $\mathbf{x}\in\mathbb{R}^{d}$, the following holds by Assumption 1

$\displaystyle F(\mathbf{x})$ $\displaystyle\leq F(0)+\left\langle\nabla
F(0),\mathbf{x}\right\rangle+\tfrac{L}{2}\|\mathbf{x}\|^{2}$
$\displaystyle\leq F(0)+\sqrt{L}\|\mathbf{x}\|+\tfrac{L}{2}\|\mathbf{x}\|^{2}$
$\displaystyle\leq F(0)+\nicefrac{{1}}{{2}}+L\|\mathbf{x}\|^{2}$

where the second inequality uses $\|\nabla F(0)\|\leq\sqrt{L}$ and the Cauchy
Schwarz inequality, and the last inequality uses the identity $ab\leq
a^{2}+\nicefrac{{b^{2}}}{{4}}$. It follows that,

$\displaystyle\mathbb{E}_{\mathbf{x}\sim\mu_{0}}[F(\mathbf{x})]\leq
F(0)+\nicefrac{{1}}{{2}}+LR^{2}$

By a slight abuse of notation, let $\mu_{0}$ denote the density of
$\mathsf{Uniform}(\mathcal{B}(R))$. Clearly.
$\mu_{0}(\mathbf{x})=\tfrac{1}{\mathsf{vol}(\mathcal{B}(R))}\mathbb{I}_{\mathbf{x}\in\mathcal{B}(R)}$.
It follows that,

$\displaystyle\int_{\mathbb{R}^{d}}\mu_{0}(\mathbf{x})\ln(\mu_{0}(\mathbf{x}))\mathsf{d}\mathbf{x}=\int_{\mathcal{B}(R)}\tfrac{1}{\mathsf{vol}(\mathcal{B}(R))}\log(\nicefrac{{1}}{{\mathsf{vol}(\mathcal{B}(R))}})\mathsf{d}\mathbf{x}=-\log(\mathsf{vol}(\mathcal{B}(R)))$

Now,
$\mathsf{vol}(\mathcal{B}(R))=\tfrac{\pi^{\nicefrac{{d}}{{2}}}}{\Gamma(\tfrac{d}{2}+1)}R^{d}$.
Furthermore, by Stirling’s Approximation,
$(\nicefrac{{x}}{{e}})^{x-1}\leq\Gamma(x)\leq(\nicefrac{{x}}{{2}})^{x-1}$.
Hence,

$\displaystyle\frac{d}{2}\log\left(\frac{2\pi
R^{2}}{\nicefrac{{d}}{{2}}+1}\right)\leq\log(\mathsf{vol}(\mathcal{B}(R)))\leq\frac{d}{2}\log\left(\frac{e\pi
R^{2}}{\nicefrac{{d}}{{2}}+1}\right)$

Without loss of generality, assume
$\pi^{\star}(\mathbf{x})=e^{-F(\mathbf{x})}$ (this can be easily ensured by
appropriately adjusting $F(0)$ upto constant factors). It follows that,

$\displaystyle\mathsf{KL}\left(\mu_{0}\boldsymbol{||}\pi^{\star}\right)$
$\displaystyle=\int_{\mathbb{R}^{d}}\mu_{0}(\mathbf{x})\log(\tfrac{\mu_{0}(\mathbf{x})}{\pi^{\star}(\mathbf{x})})\mathsf{d}\mathbf{x}=\int_{\mathbb{R}^{d}}\mu_{0}(\mathbf{x})\ln(\mu_{0}(\mathbf{x}))\mathsf{d}\mathbf{x}+\mathbb{E}_{\mathbf{x}\sim\mu_{0}}[F(\mathbf{x})]$
$\displaystyle\leq-\log(\mathsf{vol}(\mathcal{B}(R)))+F(0)+\nicefrac{{1}}{{2}}+LR^{2}$
$\displaystyle\leq\frac{d}{2}\log\left(\frac{\nicefrac{{d}}{{2}}+1}{2\pi
R^{2}}\right)+F(0)+\nicefrac{{1}}{{2}}+LR^{2}$

Setting $R=\sqrt{\nicefrac{{d}}{{L}}}$, we conclude that,

$\displaystyle\mathsf{KL}\left(\mu_{0}\boldsymbol{||}\pi^{\star}\right)\leq\frac{d}{2}\log(\nicefrac{{L}}{{2\pi}})+d+F(0)+\nicefrac{{1}}{{2}}\leq
O(d)$

∎

We now show that the growth condition on $F$, i.e. Assumption 2 is more
general than specific concentration assumptions on $\pi^{\star}$ (e.g.
subgaussianity, subexponentiality etc.). To this end, we define the notion of
$\alpha$-tail decay as follows:

######  Definition 2 ($\alpha$-Tail Decay).

A probability distribution $\nu$ on $\mathbb{R}^{d}$ is said to satisfy
$\alpha$-tail decay for some $\alpha>0$ if there exists some $C>0$ such that
$\mathbb{E}_{\mathbf{x}\sim\nu}\left[\exp\left(\|\tfrac{\mathbf{x}}{C}\|^{\alpha}\right)\right]<\infty$

The $\alpha$-tail decay condition essentially implies that the tails of
$\pi^{\star}$ decay as $\propto e^{-\|\mathbf{x}\|^{\alpha}}$. In particular,
Vershynin [47, Proposition 2.5.2 and Proposition 2.7.1] shows that
$\pi^{\star}$ satisfying the tail decay condition with $\alpha=2$ is
equivalent to $\pi^{\star}$ being subgaussian, whereas tail deay with
$\alpha=1$ is equivalent to $\pi^{\star}$ being subexponential.

In the following lemma, we establish that, under smoothness of $F$, the
$\alpha$-tail decay condition is equivalent to the growth condition on $F$
with the same exponent $\alpha$. Consequently, Assumption 2 is much weaker
than the standard isoperimetric and information transport assumptions
generally used in the literature.

######  Lemma 5 (Growth Condition and Tail Decay).

Let Assumption 2 be satisfied for some $\alpha>0$. Then, $\pi^{\star}$
satisfies the $\alpha$-tail decay condition. Conversely, let Assumption 1 be
satisfied and suppose $\pi^{\star}$ satisfies the $\alpha$-tail decay
condition. Then, $F$ satisfies Assumption 2 with the same exponent $\alpha$.

###### Proof.

Growth Condition Implies Tail Decay Since Assumption 2 is satisfied,
$F(\mathbf{x})\geq d_{1}\|\mathbf{x}\|^{\alpha}-d_{2}$ for some
$d_{1},d_{2},\alpha>0$. Let
$C=(\nicefrac{{2}}{{d_{1}}})^{\nicefrac{{1}}{{\alpha}}}$. It follows that,

$\displaystyle\mathbb{E}_{\mathbf{x}\sim\pi^{\star}}\left[e^{\|\nicefrac{{\mathbf{x}}}{{C}}\|^{\alpha}}\right]$
$\displaystyle=\int_{\mathbb{R}^{d}}e^{\tfrac{d_{1}}{2}\|\mathbf{x}\|^{\alpha}}\pi^{\star}(\mathbf{x})\mathsf{d}\mathbf{x}$
$\displaystyle\leq\int_{\mathbb{R}^{d}}e^{\tfrac{d_{1}}{2}\|\mathbf{x}\|^{\alpha}-d_{1}\|\mathbf{x}\|^{\alpha}+d_{2}}\mathsf{d}\mathbf{x}$
$\displaystyle=e^{d_{2}}\int_{\mathbb{R}^{d}}e^{-\tfrac{d_{1}}{2}\|\mathbf{x}\|^{\alpha}}\mathsf{d}\mathbf{x}<\infty$

From Definition 2, we conclude that $\pi^{\star}$ satisfies $\alpha$-tail
decay.

Smoothness and Tail Decay Imply the Growth Condition Since $F$ is smooth, it
suffices to consider $\alpha\in(0,2]$. By Assumption 1, the following
inequalities hold,

$\displaystyle F(\mathbf{y})-F(\mathbf{x})\leq\|\nabla
F(\mathbf{x})\|\|\mathbf{y}-\mathbf{x}\|+\frac{L}{2}\|\mathbf{y}-\mathbf{x}\|^{2}\leq(L\|\mathbf{x}\|+\sqrt{L})\|\mathbf{y}-\mathbf{x}\|+\frac{L}{2}\|\mathbf{y}-\mathbf{x}\|^{2}$
(7)

Ww now prove this result by contradiction. Since $\pi^{\star}$ satisfies
$\alpha$-tail decay, there exists a constant $C>0$ such that
$\mathbb{E}_{\mathbf{x}\sim\pi^{\star}}[e^{\|\nicefrac{{\mathbf{x}}}{{C}}\|^{\alpha}}]<\infty$.
Now, suppose $F$ does not satisfy the growth condition with exponent $\alpha$,
i.e., _assume there does not exist_ any $d_{1},d_{2}>0$ such that
$F(\mathbf{x})\geq d_{1}\|\mathbf{x}\|^{\alpha}-d_{2}\ \forall\
\mathbf{x}\in\mathbb{R}^{d}$. This implies that,
$\liminf_{\|\mathbf{x}\|\to\infty}\frac{F(\mathbf{x})}{\|\mathbf{x}\|^{\alpha}}=0$.
Thus, without loss of generality, we can assume there exists a diverging
sequence $a_{n}\in\mathbb{R}$ and a diverging sequence
$\mathbf{x}_{n}\in\mathbb{R}^{d}$ that satisfy the following for every
$n\in\mathbb{N}$ :

$\displaystyle\frac{F(\mathbf{x}_{n})}{\|\mathbf{x}_{n}\|^{\alpha}}\leq\frac{1}{a_{n}},\quad\|\mathbf{x}_{n}\|\geq
2n,\quad\|\mathbf{x}_{n+1}-\mathbf{x}_{n}\|\geq 1$ (8)

where, without loss of generality, we assume $a_{n},\|\mathbf{x}_{n}\|$ > 0\.
Now, let $r_{n}=\tfrac{1}{\|\mathbf{x}_{n}\|^{2}}$ and
$B_{n}\subseteq\mathbb{R}^{d}$ denote the ball of radius $r_{n}$ centered at
$\mathbf{x}_{n}$. Since $r_{n}\leq\nicefrac{{1}}{{4n^{2}}}$ and
$\|\mathbf{x}_{n+1}-\mathbf{x}_{n}\|\geq 1$, $B_{n}$ is a family of disjoint
subsets of $\mathbb{R}^{d}$. We shall now prove that there exists some
diverging sequence $b_{n}\in\mathbb{R}$ such that
$\tfrac{F(\mathbf{y})}{\|\mathbf{y}\|^{\alpha}}\leq\tfrac{1}{b_{n}}$ for every
$\mathbf{y}\in B_{n}$.

Consider any arbitrary $n\in\mathbb{N}$ and let $\mathbf{y}\in B_{n}$.
Applying (7) to $\mathbf{y}$ and $\mathbf{x}_{n}$, we obtain,

$\displaystyle\frac{F(\mathbf{y})}{\|\mathbf{x}_{n}\|^{\alpha}}$
$\displaystyle\leq\frac{F(\mathbf{x}_{n})}{\|\mathbf{x}_{n}\|^{\alpha}}+\frac{L\|\mathbf{x}_{n}\|r_{n}}{\|\mathbf{x}_{n}\|^{\alpha}}+\frac{r_{n}\sqrt{L}}{\|\mathbf{x}_{n}\|^{\alpha}}+\frac{Lr^{2}_{n}}{2\|\mathbf{x}_{n}\|^{\alpha}}$
$\displaystyle\leq\frac{1}{a_{n}}+\frac{L}{\|\mathbf{x}\|^{\alpha+1}}+\frac{\sqrt{L}}{\|\mathbf{x}_{n}\|^{\alpha+2}}+\frac{L}{2\|\mathbf{x}_{n}\|^{\alpha+4}}$
(9)

where we use (8) and $r_{n}=\nicefrac{{1}}{{\|\mathbf{x}_{n}\|^{2}}}$.
Moreover, we note that

$\displaystyle\|\mathbf{y}\|\geq\|\mathbf{x}_{n}\|-\|\mathbf{y}-\mathbf{x}_{n}\|\geq\|\mathbf{x}_{n}\|-r_{n}=\|\mathbf{x}_{n}\|-\frac{1}{\|\mathbf{x}_{n}\|^{2}}\geq\frac{\|\mathbf{x}_{n}\|}{2}$
(10)

where we use the fact that $\|\mathbf{x}_{n}\|\geq
2n>2^{\nicefrac{{1}}{{3}}}$. It follows that,

$\displaystyle\frac{F(\mathbf{y})}{\|\mathbf{y}\|^{\alpha}}$
$\displaystyle\leq\frac{2^{\alpha}F(\mathbf{y})}{\|\mathbf{x}_{n}\|^{\alpha}}$
$\displaystyle\leq\frac{4}{a_{n}}+\frac{4L}{\|\mathbf{x}\|^{\alpha+1}}+\frac{4\sqrt{L}}{\|\mathbf{x}_{n}\|^{\alpha+2}}+\frac{2L}{\|\mathbf{x}_{n}\|^{\alpha+4}}$
(11)

where we use (B) and the fact that $\alpha\in(0,2]$. We now define the
sequence $b_{n}\in\mathbb{R}$ as follows:

$\displaystyle b_{n}$
$\displaystyle=\left(\frac{4}{a_{n}}+\frac{4L}{\|\mathbf{x}\|^{\alpha+1}}+\frac{4\sqrt{L}}{\|\mathbf{x}_{n}\|^{\alpha+2}}+\frac{2L}{\|\mathbf{x}_{n}\|^{\alpha+4}}\right)^{-1}$

Since $\alpha>0$, and $a_{n},\|\mathbf{x}_{n}\|\to\infty$, it is clear that
$b_{n}$ is a diverging sequence. Furthermore, from (B), we conclude that
$\tfrac{F(\mathbf{y})}{\|\mathbf{y}\|^{\alpha}}\leq\tfrac{1}{b_{n}}\ \forall\
\mathbf{y}\in B_{n}$. Equipped with this construction, we note that

$\displaystyle\mathbb{E}_{\mathbf{x}\sim\pi^{\star}}\left[\exp\left(\frac{\|\mathbf{x}\|^{\alpha}}{C^{\alpha}}\right)\right]$
$\displaystyle=\int_{\mathbb{R}^{d}}\exp\left(\frac{\|\mathbf{y}\|^{\alpha}}{C^{\alpha}}\right)\exp(-F(\mathbf{y}))\mathsf{d}\mathbf{y}$
$\displaystyle\geq\sum_{n=1}^{\infty}\int_{B_{n}}\exp\left(\frac{\|\mathbf{y}\|^{\alpha}}{C^{\alpha}}\right)\exp(-F(\mathbf{y}))\mathsf{d}\mathbf{y}$
$\displaystyle\geq\sum_{n=1}^{\infty}\int_{B_{n}}\exp\left(\frac{\|\mathbf{y}\|^{\alpha}}{C^{\alpha}}-\frac{\|\mathbf{y}\|^{\alpha}}{b_{n}}\right)\mathsf{d}\mathbf{y}$

where the second inequality use the fact that $B_{n}$ is a disjoint family of
subsets of $\mathbb{R}^{d}$ and the third inequality uses the fact that
$\tfrac{F(\mathbf{y})}{\|\mathbf{y}\|^{\alpha}}\leq\tfrac{1}{b_{n}}\ \forall\
\mathbf{y}\in B_{n}$. Since $b_{n}$ is a diverging sequence, there exists some
$N_{0}\in\mathbb{N}$ such that $b_{n}\geq 2C^{\alpha}\ \forall\ n\geq N_{0}$.
It follows that,

$\displaystyle\mathbb{E}_{\mathbf{x}\sim\pi^{\star}}\left[\exp\left(\frac{\|\mathbf{x}\|^{\alpha}}{C^{\alpha}}\right)\right]$
$\displaystyle\geq\sum_{n=1}^{\infty}\int_{B_{n}}\exp\left(\frac{\|\mathbf{y}\|^{\alpha}}{C^{\alpha}}-\frac{\|\mathbf{y}\|^{\alpha}}{b_{n}}\right)\mathsf{d}\mathbf{y}$
$\displaystyle\geq\sum_{n=N_{0}}^{\infty}\int_{B_{n}}\exp\left(\frac{\|\mathbf{y}\|^{\alpha}}{2C^{\alpha}}\right)\mathsf{d}\mathbf{y}$
$\displaystyle=\sum_{n=N_{0}}^{\infty}\mathsf{vol}(B_{n})\mathbb{E}_{\mathbf{y}\sim\mathsf{Uniform}(B_{n})}\left[\exp\left(\frac{\|\mathbf{y}\|^{\alpha}}{2C^{\alpha}}\right)\right]$

Consider the function $g:[0,\infty)\to[0,\infty)$ defined as
$g(t)=e^{t^{\alpha}}$. We note that for $\alpha\geq 1$, $g$ is a convex
function for every $t\geq 0$, and for $\alpha\in(0,1)$, $g$ is convex for
every $t\geq(\nicefrac{{1}}{{\alpha}}-1)^{\nicefrac{{1}}{{\alpha}}}$. From
(10), we note that $\|\mathbf{y}\|\geq\nicefrac{{\|\mathbf{x}\|}}{{2}}\geq n$
for every $\mathbf{y}\in B_{n}$. Hence, there exists an $N_{1}\in\mathbb{N}$
such that $e^{t^{\alpha}}$ is a convex function for all
$t\geq\nicefrac{{\|\mathbf{y}\|}}{{2}},\ \forall\ \mathbf{y}\in B_{n},\ n\geq
N_{1}$. Let $N=\max\\{N_{0},N_{1}\\}+1$. Then,

$\displaystyle\mathbb{E}_{\mathbf{x}\sim\pi^{\star}}\left[\exp\left(\frac{\|\mathbf{x}\|^{\alpha}}{C^{\alpha}}\right)\right]$
$\displaystyle\geq\sum_{n=N_{0}}^{\infty}\mathsf{vol}(B_{n})\mathbb{E}_{\mathbf{y}\sim\mathsf{Uniform}(B_{n})}\left[\exp\left(\frac{\|\mathbf{y}\|^{\alpha}}{2C^{\alpha}}\right)\right]$
$\displaystyle\geq\sum_{n=N}^{\infty}\mathsf{vol}(B_{n})\exp\left(\tfrac{1}{2C^{\alpha}}\mathbb{E}_{\mathbf{y}\sim\mathsf{Uniform}(B_{n})}[\|\mathbf{y}\|]^{\alpha}\right)$
$\displaystyle\geq\sum_{n=N}^{\infty}\mathsf{vol}(B_{n})\exp\left(\tfrac{1}{2C^{\alpha}}\|\mathbb{E}_{\mathbf{y}\sim\mathsf{Uniform}(B_{n})}[\mathbf{y}]\|^{\alpha}\right)$
$\displaystyle\geq\sum_{n=N}^{\infty}\mathsf{vol}(B_{n})\exp\left(\tfrac{1}{2C^{\alpha}}\|\mathbf{x}_{n}\|^{\alpha}\right)$
$\displaystyle=\sum_{n=N}^{\infty}C_{d}(r_{n})^{d}\exp\left(\tfrac{1}{2C^{\alpha}}\|\mathbf{x}_{n}\|^{\alpha}\right)$
$\displaystyle=\sum_{n=N}^{\infty}\frac{C_{d}\exp\left(\tfrac{1}{2C^{\alpha}}\|\mathbf{x}_{n}\|^{\alpha}\right)}{\|\mathbf{x}_{n}\|^{2d}}$

where
$C_{d}=\tfrac{\pi^{\nicefrac{{d}}{{2}}}}{\Gamma(\nicefrac{{d}}{{2}}+1)}$. Let
$k$ be any positive integer such that $\alpha k\geq 2d+1$. It follows that,

$\displaystyle\frac{\exp\left(\tfrac{1}{2C^{\alpha}}\|\mathbf{x}_{n}\|^{\alpha}\right)}{\|\mathbf{x}_{n}\|^{2d}}\geq\frac{C_{d}}{2^{k}k!C^{\alpha
k}}\|\mathbf{x}_{n}\|^{\alpha k-2d}\geq\frac{C_{d}n}{2^{k-1}k!C^{\alpha k}}$

Thus, we infer that,

$\displaystyle\mathbb{E}_{\mathbf{x}\sim\pi^{\star}}\left[\exp\left(\frac{\|\mathbf{x}\|^{\alpha}}{C^{\alpha}}\right)\right]\geq\frac{C_{d}}{2^{k-1}k!C^{\alpha
k}}\sum_{n=N_{0}}^{\infty}n=\infty$

which is a contradiction. Thus, there exists some $d_{1},d_{2}>0$ such that
$F(\mathbf{x})\geq d_{1}\|\mathbf{x}\|^{\alpha}-d_{2}$, i.e., $F$ satisfies
the growth condition with exponent $\alpha$. ∎

The following lemma establishes boundedness and contractivity properties of
the function $h(\mathbf{x},\mathbf{y})=k(\mathbf{x},\mathbf{y})\nabla
F(\mathbf{y})-\nabla_{2}k(\cdot,\mathbf{y})$, that are vital for proving
almost-sure bounds such as Lemma 2.

######  Lemma 6 (Properties of $h$).

Let Assumptions 1, 3 and 4 be satisfied. Then, the following holds,

$\displaystyle\|h(.,\mathbf{y})\|_{\mathcal{H}}$ $\displaystyle\leq
BL\|\mathbf{y}\|+B\|\nabla F(0)\|+B$
$\displaystyle\|h(\mathbf{x},\mathbf{y})\|$
$\displaystyle\leq\frac{A_{1}L}{2}+A_{2}+k(\mathbf{x},\mathbf{y})\|\nabla
F(\mathbf{x})\|$ $\displaystyle-\left\langle\nabla
F(\mathbf{x}),h(\mathbf{x},\mathbf{y})\right\rangle$
$\displaystyle\leq-\tfrac{1}{2}k(\mathbf{x},\mathbf{y})\|\nabla
F(\mathbf{x})\|^{2}+L^{2}A_{1}+A_{3}$

###### Proof.

Recalling the definition of $h$ from Section 2, we observe that,

$\displaystyle h(\cdot,\mathbf{y})=k(\cdot,\mathbf{y})\nabla
F(\mathbf{y})-\nabla_{2}k(\cdot,\mathbf{y})$

Thus, by triangle inequality of $\|\cdot\|_{\mathcal{H}}$, Assumptions 1 and
3, we obtain

$\displaystyle\|h(\cdot,\mathbf{y})\|_{\mathcal{H}}$
$\displaystyle\leq\|\nabla
F(\mathbf{y})\|\|k(\cdot,\mathbf{y})\|_{\mathcal{H}_{0}}+\|\nabla_{2}k(\cdot,\mathbf{y})\|_{\mathcal{H}}$
$\displaystyle\leq BL\|\mathbf{y}\|+B\|\nabla F(0)\|+B$

To prove the remaining inequalities, we first note that,

$\displaystyle h(\mathbf{x},\mathbf{y})$
$\displaystyle=k(\mathbf{x},\mathbf{y})\nabla
F(\mathbf{y})-\nabla_{2}k(\mathbf{x},\mathbf{y})$
$\displaystyle=k(\mathbf{x},\mathbf{y})\nabla
F(\mathbf{x})+k(\mathbf{x},\mathbf{y})\left[\nabla F(\mathbf{y})-\nabla
F(\mathbf{x})\right]-\nabla_{2}k(\mathbf{x},\mathbf{y})$ (12)

Using Assumptions 1 and 4, we note that,

$\displaystyle\|h(\mathbf{x},\mathbf{y})\|$ $\displaystyle\leq
k(\mathbf{x},\mathbf{y})\|\nabla
F(\mathbf{x})\|+\frac{LA_{1}\|\mathbf{x}-\mathbf{y}\|}{1+\|\mathbf{x}-\mathbf{y}\|^{2}}+A_{2}$
$\displaystyle\leq\frac{A_{1}L}{2}+A_{2}+k(\mathbf{x},\mathbf{y})\|\nabla
F(\mathbf{x})\|$

where the second inequality uses the fact
$\tfrac{t}{1+t^{2}}\leq\nicefrac{{1}}{{2}}$

To prove the last inequality, we infer the following from (B)

$\displaystyle-\left\langle\nabla
F(\mathbf{x}),h(\mathbf{x},\mathbf{y})\right\rangle$
$\displaystyle\leq-k(\mathbf{x},\mathbf{y})\|\nabla
F(\mathbf{x})\|^{2}+k(\mathbf{x},\mathbf{y})\|\nabla F(\mathbf{x})-\nabla
F(\mathbf{y})\|\|\nabla F(\mathbf{x})\|$
$\displaystyle+\|\nabla_{2}k(\mathbf{x},\mathbf{y})\|\|\nabla F(\mathbf{x})\|$
$\displaystyle\leq-k(\mathbf{x},\mathbf{y})\|\nabla
F(\mathbf{x})\|^{2}+L\sqrt{k(\mathbf{x},\mathbf{y})}\sqrt{\frac{A_{1}\|\mathbf{x}-\mathbf{y}\|^{2}}{1+\|\mathbf{x}-\mathbf{y}\|^{2}}}\|\nabla
F(\mathbf{x})\|$ $\displaystyle+\sqrt{A_{3}k(\mathbf{x},\mathbf{y})}\|\nabla
F(\mathbf{x})\|$
$\displaystyle\leq-\frac{1}{2}k(\mathbf{x},\mathbf{y})\|\nabla
F(\mathbf{x})\|^{2}+L^{2}A_{1}+A_{3}$

where the second inequality uses Assumptions 1 and 4, and the last inequality
uses the identity $ab\leq a^{2}+\nicefrac{{b^{2}}}{{4}}$ ∎

To analyze the dynamics of VP-SVGD in the Wasserstein space, we use the
following lemma presented in Salim et al. [41]

######  Lemma 7 (Salim et al. [41], Proposition 3.1).

Let Assumptions 1 and 3 be satisfied. Consider any
$\nu_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ with
$\mathsf{KL}\left(\nu_{0}\boldsymbol{||}\pi^{\star}\right)<\infty$,
$f\in\mathcal{H}$ and let $\nu_{1}=(I-\eta f)_{\\#}\nu_{0}$ with
$\eta\|f\|_{\mathcal{H}}\leq\tfrac{\beta-1}{\beta B}$ for some $\beta>1$.
Then, the following holds,

$\displaystyle\mathsf{KL}\left(\mu_{1}\boldsymbol{||}\pi^{\star}\right)\leq\mathsf{KL}\left(\mu_{0}\boldsymbol{||}\pi^{\star}\right)-\eta\left\langle
h_{\mu_{0}},f\right\rangle+\frac{\eta^{2}(\beta^{2}+L)B}{2}\|f\|^{2}_{\mathcal{H}}$

###  B.1 Gelfand-Pettis Integrals for Reproducing Kernel Hilbert Spaces

The Gelfand-Pettis integral is a generalization of the Lebesgue integral to
functions that take values in an arbitrary topological vector space. In this
section, we describe the Gelfand-Pettis integral for an arbitrary Hilbert
space $(V,\left\langle\cdot,\cdot\right\rangle_{V})$ and refer the readers to
Talagrand [45] for a more general treatment.

Let $(X,\Sigma,\lambda)$ be a measure space and
$(V,\left\langle\cdot,\cdot\right\rangle_{V})$ be a Hilbert Space. A function
$g:X\to V$ is said to be Gelfand-Pettis integrable if there exists a vector
$w_{g}\in V$ such that $\left\langle
u,w_{g}\right\rangle_{V}=\int_{X}\left\langle
u,g(x)\right\rangle_{V}\mathsf{d}\lambda(x)\ \forall\ u\in V$. The vector
$w_{g}$ is called the Gelfand-Pettis integral of $g$

We now establish the following lemma for Gelfand-Pettis integrals with respect
to the RKHS $\mathcal{H}$, which is a key component of our analysis of VP-
SVGD.

######  Lemma 8.

Let $\mu$ be a probability measure on $\mathbb{R}^{d}$. Let
$G:\mathbb{R}^{d}\times\mathbb{R}^{d}$ be a function such that for every
$\mathbf{y}\in\mathbb{R}^{d}$, $G(.,\mathbf{y})\in\mathcal{H}$ with
$\|G(.,\mathbf{y})\|_{\mathcal{H}}\leq C$ holding $\mu$-almost surely. Let
$G_{\mu}(\mathbf{x})=\mathbb{E}_{\mathbf{y}\sim\mu}[G(\mathbf{x},\mathbf{y})]$.
Then, the map $\psi:\mathbb{R}^{d}\to\mathcal{H}$ defined as
$\psi(\mathbf{y})=G(\cdot,\mathbf{y})$ is Gelfand-Pettis integrable and
$G_{\mu}$ is the Gelfand-Pettis integral of $\psi$ with respect to $\mu$, i.e.
$G_{\mu}\in\mathcal{H}$ and for any $f\in\mathcal{H}$,
$\mathbb{E}_{\mathbf{y}\sim\mu}[\left\langle
f,G(.,\mathbf{y})\right\rangle_{\mathcal{H}}]=\left\langle
f,G_{\mu}\right\rangle_{\mathcal{H}}$

###### Proof.

Let $\Phi:\mathcal{H}\to\mathbb{R}$ denote the map
$\Phi(f)=\mathbb{E}_{\mathbf{y}\sim\mu}[\left\langle
f,G(.,\mathbf{y})\right\rangle_{\mathcal{H}}]\ \forall\ f\in\mathcal{H}$. By
linearity of expectations and inner products, we note that $\Phi$ is a linear
functional on $\mathcal{H}$. Furthermore, since
$\|G(.,\mathbf{y})\|_{\mathcal{H}}\leq C$ holds $\mu$-almost surely, we note
that for any $f\in\mathcal{H}$,
$|\Phi(f)|\leq\mathbb{E}_{\mathbf{y}\sim\mu}[|\left\langle
f,G(.,\mathbf{y})\right\rangle_{\mathcal{H}}|]\leq C\|f\|_{\mathcal{H}}$ by
Jensen’s inequality and Cauchy Schwarz inequality for $\mathcal{H}$. We
conclude that $\Phi$ is a bounded linear functional of $\mathcal{H}$. Thus, by
Reisz Representation Theorem [9], there exists $g\in\mathcal{H}$ such that for
any $f\in\mathcal{H}$, the following holds

$\displaystyle\mathbb{E}_{\mathbf{y}\in\mu}[\left\langle
f,G(.,\mathbf{y})\right\rangle_{\mathcal{H}}]=\left\langle
f,g\right\rangle_{\mathcal{H}}$

Hence, we conclude that the map $\psi$ is Gelfand-Pettis integrable. We now
use the reproducing property of $\mathcal{H}$ to show that $g=G_{\mu}$, i.e.,
$G_{\mu}$ is the Gelfand-Pettis integral of $\psi$. To this end, let
$\mathbf{x}\in\mathbb{R}^{d}$ be arbitrary. Setting $f=k(\mathbf{x},.)$ and
using the fact that $g\in\mathcal{H}$, $G(.,\mathbf{y})\in\mathcal{H}$ for any
$\mathbf{y}\in\mathbb{R}^{d}$,

$\displaystyle
g(\mathbf{x})=\mathbb{E}_{\mathbf{y}\in\mu}[G(\mathbf{x},\mathbf{y})]=G_{\mu}(\mathbf{x})$

Hence, $g=G_{\mu}$, i.e., $\mathbb{E}_{\mathbf{y}\sim\mu}[\left\langle
f,G(.,\mathbf{y})\right\rangle_{\mathcal{H}}]=\left\langle
f,G_{\mu}\right\rangle_{\mathcal{H}}$ ∎

##  Appendix C Analysis of VP-SVGD

In this section, we present our analysis of VP-SVGD. Throughout this section,
we define the random function $g_{t}:\mathbb{R}^{d}\times\mathbb{R}^{d}$ as
$g_{t}(\mathbf{x})=\tfrac{1}{K}\sum_{l=0}^{K-1}h(\mathbf{x},\mathbf{x}^{(Kt+l)}_{t})$
where $t\in\mathbb{N}\cup\\{0\\}$, $K$ is the batch-size of VP-SVGD, and
$h:\mathbb{R}^{d}\times\mathbb{R}^{d}$ is as defined in Section 2, i.e.,
$h(\mathbf{x},\mathbf{y})=k(\mathbf{x},\mathbf{y})\nabla
F(\mathbf{y})-\nabla_{2}k(\mathbf{x},\mathbf{y})$.

After proving the key lemmas required for our analysis of VP-SVGD, we present
the proof of Theorem 1 in Appendix C.4. We also present a high-probability
version of Theorem 1 in Appendix C.5

###  C.1 Population Level Dynamics : Proof of Lemma 1

###### Proof.

We now derive the population-limit dynamics of VP-SVGD for arbitrary batch-
size $K$, and subsequently prove the descent lemma (i.e. Lemma 1) for VP-SVGD.
The arguments of this section are a straightforward generalization of that
used in Section 4.

To this end, we recall from Section 4 that the countably infinite number of
particles $\mathbf{x}^{(l)}_{0}$, $l\in\mathbb{N}\cup\\{0\\}$ are i.i.d
samples from the measure $\mu_{0}$, which has a density w.r.t the Lebesgue
measure. Thus, by the strong law of large numbers (Dudley [12, Theorem
11.4.1]), the empirical measure of $(\mathbf{x}^{(l)}_{0})_{l\geq 0}$ is
almost surely equal to $\mu_{0}$. Furthermore, we recall the filtration
$\mathcal{F}_{t}$ defined in Section 4 as
$\mathcal{F}_{t}=\sigma(\mathbf{x}^{(l)}_{0}\ |\ l\leq Kt-1)$,
$t\in\mathbb{N}$ with $\mathcal{F}_{0}$ being the trivial $\sigma$ algebra. We
now consider the following dynamics in $\mathbb{R}^{d}$:

$\mathbf{x}^{(s)}_{t+1}=\mathbf{x}^{(s)}_{t}-\frac{\gamma}{K}\sum_{l=0}^{K-1}h(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(tK+l)}_{t}),\quad
s\in\mathbb{N}\cup\\{0\\}$ (13)

We note that the above updates are the same as that of VP-SVGD for
$s\in\\{0,\dots,KT+n-1\\}$. Now, for each time-step $t$, we focus on the lower
triangular evolution, i.e., the time evolution of the particles
$(\mathbf{x}^{(l)}_{t})_{l\geq Kt}$. From (13), we infer that for any
$t\in\mathbb{N}$ and $s\geq Kt$, $\mathbf{x}^{(s)}_{t}$ depends only on
$(\mathbf{x}^{(l)}_{0})_{l\leq Kt-1}$ and $\mathbf{x}^{(s)}_{0}$. Hence, there
exists a measurable function $H_{t}$ for every $t\in\mathbb{N}$ such that the
following holds almost surely:

$\mathbf{x}^{(s)}_{t}=H_{t}(\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0},\mathbf{x}^{(s)}_{0});\quad\forall\
s\geq Kt$ (14)

Since
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0},\mathbf{x}^{(s)}_{0}\
\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\ \mu_{0}$, we conclude from (14) that
$(\mathbf{x}^{(s)}_{t})_{s\geq Kt}$ are i.i.d when conditioned on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0}$. To this end, we define
the random measure $\mu_{t}|\mathcal{F}_{t}$ as the law of
$\mathbf{x}^{(Kt)}_{t}$ conditioned on $\mathcal{F}_{t}$, i.e.
$\mu_{t}|\mathcal{F}_{t}$ is a probability kernel
$\mu_{t}(\cdot;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0})$ with
$\mu_{0}|\mathcal{F}_{0}\coloneqq\mu_{0}$. By the strong law of large numbers,
$\mu_{t}|\mathcal{F}_{t}$ is equal to the empirical measure of
$(\mathbf{x}^{(l)}_{t})_{l\geq Kt}$ conditioned on $\mathcal{F}_{t}$.
Furthermore, we infer from (13) that the particles satisfy the following:

$\displaystyle\mathbf{x}^{(s)}_{t+1}=(I-\gamma
g_{t})(\mathbf{x}^{(s)}_{t}),\quad s\geq K(t+1)$

Recall that
$\mathbf{x}_{t+1}^{(s)}|\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(K(t+1)-1)}_{0}\sim\mu_{t+1}|\mathcal{F}_{t+1}$
for any $s\geq K(t+1)$. Furthermore, from Equation (14), we note that for
$s\geq K(t+1)$, $\mathbf{x}^{(s)}_{t}$ depends only on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0}$ and
$\mathbf{x}^{(s)}_{0}$, which implies that
$\mathrm{Law}(\mathbf{x}^{(s)}_{t}|\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(K(t+1)-1)}_{0})=\mathrm{Law}(\mathbf{x}^{(s)}_{t}|\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0})=\mu_{t}|\mathcal{F}_{t}$.
Finally, we note that $g_{t}$ is an $\mathcal{F}_{t+1}$-measurable random
function. With these insights, we conclude that the population-level dynamics
of the lower triangular evolution in $\mathcal{P}_{2}(\mathbb{R}^{d})$ is
almost surely described by the following update:

$\mu_{t+1}|\mathcal{F}_{t+1}=(I-\gamma g_{t})_{\\#}\mu_{t}|\mathcal{F}_{t}$
(15)

Setting $\gamma\|g_{t}\|_{\mathcal{H}}\leq\tfrac{\beta-1}{\beta B}$ for some
arbitrary $\beta>1$ and applying Lemma 7 to the population-level update (15),
we conclude that the following holds almost surely:

$\displaystyle\mathsf{KL}\left(\mu_{t+1}|\mathcal{F}_{t+1}\boldsymbol{||}\pi^{\star}\right)\leq\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\gamma\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},g_{t}\right\rangle_{\mathcal{H}}+\frac{\gamma^{2}(\beta^{2}+L)B}{2}\|g_{t}\|^{2}_{\mathcal{H}}$

∎

###  C.2 Iterate Bounds : Proof of Lemma 2

To establish almost sure bounds on $\|g_{t}\|_{\mathcal{H}}$, we prove the
following result which is stronger than Lemma 2.

######  Lemma 9 (Almost-Sure Iterate Bounds for VP-SVGD).

Let Assumptions 1, 2, 3, 4 and 5 be satisfied. Then, the following holds
almost surely for any $s\in\mathbb{N}\cup\\{0\\}$ and $t\in(T+1)$ whenever
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$

$\displaystyle\|\mathbf{x}^{(s)}_{t}\|$
$\displaystyle\leq\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$
$\displaystyle\|h(\cdot,\mathbf{x}^{(s)}_{t})\|_{\mathcal{H}}$
$\displaystyle\leq\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$
$\displaystyle\|g_{t}\|_{\mathcal{H}}$
$\displaystyle\leq\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$

where $\zeta_{0},\dots,\zeta_{3}$ are problem-dependent constants that depend
polynomially on $A_{1},A_{2},A_{3},B,d_{1},d_{2},L$ for any fixed $\alpha$.

###### Proof.

Let
$c^{(s)}_{t}=\tfrac{1}{K}\sum_{l=0}^{K-1}k(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(Kt+l)}_{t})$.
Note that by Assumption 4, $c^{(s)}_{t}\geq 0$ Since
$\mathbf{x}^{(s)}_{t+1}=\mathbf{x}^{(s)}_{t}-\gamma
g_{t}(\mathbf{x}^{(s)}_{t})$, it follows from the smoothness of $F$ that,

$F(\mathbf{x}^{(s)}_{t+1})-F(\mathbf{x}^{(s)})\leq-\gamma\left\langle\nabla
F(\mathbf{x}^{(s)}_{t}),g_{t}(\mathbf{x}^{(s)}_{t})\right\rangle+\frac{\gamma^{2}L}{2}\|g_{t}(\mathbf{x}^{(s)}_{t})\|^{2}$
(16)

By Lemma 6, we note that,

$\displaystyle-\gamma\left\langle\nabla
F(\mathbf{x}^{(s)}_{t}),g_{t}(\mathbf{x}^{(s)}_{t})\right\rangle$
$\displaystyle=-\frac{\gamma}{K}\sum_{l=0}^{K-1}\left\langle\nabla
F(\mathbf{x}^{(s)}_{t}),h(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(tK+l)}_{t})\right\rangle$
$\displaystyle\leq\frac{\gamma}{K}\sum_{l=0}^{L-1}\left[-\tfrac{1}{2}k(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(tK+l)}_{t})\|\nabla
F(\mathbf{x}^{(s)}_{t})\|^{2}+L^{2}A_{1}+A_{3}\right]$
$\displaystyle\leq-\frac{\gamma c^{(s)}_{t}}{2}\|\nabla
F(\mathbf{x}^{(s)}_{t})\|^{2}+\gamma L^{2}A_{1}+\gamma A_{3}$ (17)

Moreover, by Jensen’s Inequality and Lemma 6

$\displaystyle\|g_{t}(\mathbf{x}^{(s)}_{t})\|^{2}$
$\displaystyle\leq\frac{1}{K}\sum_{l=0}^{K-1}\|h(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(Kt+l)}_{t})\|^{2}$
$\displaystyle\leq\frac{1}{K}\sum_{l=0}^{K-1}2(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}+2k(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(tK+l)}_{t})^{2}\|F(\mathbf{x}^{(s)}_{t})\|^{2}$
$\displaystyle\leq\frac{1}{K}\sum_{l=0}^{K-1}2(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}+2A_{1}k(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(tK+l)}_{t})\|F(\mathbf{x}^{(s)}_{t})\|^{2}$
$\displaystyle\leq
2(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}+2A_{1}c^{(s)}_{t}\|F(\mathbf{x}^{(s)}_{t})\|^{2}$
(18)

Substituting (C.2) and (C.2) into (16), we obtain,

$\displaystyle F(\mathbf{x}^{(s)}_{t+1})-F(\mathbf{x}^{(s)}_{t})$
$\displaystyle\leq-\frac{\gamma c^{(s)}_{t}}{2}\|\nabla
F(\mathbf{x}^{(s)}_{t})\|^{2}+\gamma L^{2}A_{1}+\gamma A_{3}$
$\displaystyle+\gamma^{2}L(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}+\gamma^{2}LA_{1}c^{(s)}_{t}\|F(\mathbf{x}^{(s)}_{t})\|^{2}$
$\displaystyle\leq-\frac{\gamma c^{(s)}_{t}}{2}(1-2A_{1}L\gamma)\|\nabla
F(\mathbf{x}^{(s)}_{t})\|^{2}+\gamma A_{3}+\gamma
L^{2}A_{1}+\gamma^{2}L(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}$
$\displaystyle\leq\gamma A_{3}+\gamma
L^{2}A_{1}+\gamma^{2}L(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}$

where the last inequality uses the fact that $c^{(s)}_{t}\geq 0$ and
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$. Now, iterating through the above
inequality, we obtain the following for any $t\in[T],\
s\in\mathbb{N}\cup\\{0\\}$

$F(\mathbf{x}^{(s)}_{t})\leq F(\mathbf{x}^{(s)}_{0})+\gamma TL^{2}A_{1}+\gamma
TA_{3}+\gamma^{2}TL(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}$ (19)

Furthermore, by Assumption 1

$\displaystyle F(\mathbf{x}^{(s)}_{0})$ $\displaystyle\leq F(0)+\|\nabla
F(0)\|\|\mathbf{x}^{(s)}_{0}\|+\frac{L}{2}\|\mathbf{x}^{(s)}_{0}\|^{2}$
$\displaystyle\leq F(0)+\nicefrac{{1}}{{2}}+L\|\mathbf{x}^{(s)}_{0}\|^{2}$

Substituting the above inequality into (19), and using Assumption 2, we obtain
the following for any $t\in[T],\ s\in\mathbb{N}\cup\\{0\\}$

$\displaystyle d_{1}\|\mathbf{x}^{(s)}_{t}\|^{\alpha}-d_{2}\leq
F(\mathbf{x}^{s}_{t})$ $\displaystyle\leq
F(0)+\nicefrac{{1}}{{2}}+L\|\mathbf{x}^{(s)}_{0}\|^{2}+\gamma
TL^{2}A_{1}+\gamma TA_{3}$
$\displaystyle+\gamma^{2}TL(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}$

Rearranging and applying Assumption 5, we obtain

$\displaystyle\|\mathbf{x}^{(s)}_{t}\|$ $\displaystyle\leq
d^{-\nicefrac{{1}}{{\alpha}}}_{1}\left[F(0)+\nicefrac{{1}}{{2}}+LR^{2}+\gamma
TL^{2}A_{1}+\gamma
TA_{3}+\gamma^{2}TL(A_{1}L+A_{2})^{2}\right]^{\nicefrac{{1}}{{\alpha}}}$
$\displaystyle\leq\tilde{\zeta}_{0}+\tilde{\zeta}_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\zeta}_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\zeta}_{3}R^{\nicefrac{{2}}{{\alpha}}}$

where $\tilde{\zeta}_{0},\dots,\tilde{\zeta}_{3}$ are constants that depend
polynomially on $L,A_{1},A_{2},A_{3},R$. We note that, since $0<\alpha\leq 2$,
the above inequality also holds for $t=0$.

Using the above inequality Lemma 6 and Assumption 1, we conclude that the
following holds almost surely for any $t\in(T+1),s\in\mathbb{N}\cup\\{0\\}$

$\displaystyle\|h(\cdot,\mathbf{x}^{(s)}_{t})\|_{\mathcal{H}}$
$\displaystyle\leq BL\|\mathbf{x}^{(s)}_{t}\|+B\sqrt{L}+B$
$\displaystyle\leq\tilde{\eta_{0}}+\tilde{\eta}_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\eta}_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\eta}_{3}R^{\nicefrac{{2}}{{\alpha}}}$

where $\tilde{\eta}_{0},\dots,\tilde{\eta}_{3}$ are constants that depend
polynomially on $L,B,A_{1},A_{2},A_{3},R$. Using the above inequality, we
conclude that the following also holds for any $t\in(T+1)$.

$\displaystyle\|g_{t}\|_{\mathcal{H}}\leq\tilde{\eta_{0}}+\tilde{\eta}_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\eta}_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\eta}_{3}R^{\nicefrac{{2}}{{\alpha}}}$

Taking $\zeta_{i}=\max\\{\tilde{\zeta_{i}},\tilde{\eta}_{i}\\}$, the proof is
complete. ∎

###  C.3 Controlling $g_{t}$ in Expectation : Proof of Lemma 3

###### Proof.

Let $\xi=\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$
where $\zeta_{0},\dots,\zeta_{3}$ are as defined in Lemma 9. Recall that
$g_{t}=\tfrac{1}{K}\sum_{l=0}^{K-1}h(.,\mathbf{x}^{(Kt+l)}_{t})$. Since
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$,
$\|h(\cdot,\mathbf{x}^{(Kt+l)}_{t})\|_{\mathcal{H}}\leq\xi$ holds almost
surely y Lemma 9.

Consider any $l\in(K)$. Conditioned on the filtration $\mathcal{F}_{t}$,
$\mathrm{Law}(\mathbf{x}^{(Kt+l)}_{t}|\mathcal{F}_{t})=\mu_{t}|\mathcal{F}_{t}$.
Moreover, for any $\mathbf{x}\in\mathbb{R}^{d}$,
$\mathbb{E}_{\mathbf{x}^{(Kt+l)}_{t}}[h(\mathbf{x},\mathbf{x}^{(Kt+l)}_{t})|\mathcal{F}_{t}]=h_{\mu_{t}|\mathcal{F}_{t}}(\mathbf{x})$.
Thus, from Lemma 8, we conclude that $h_{\mu_{t}|\mathcal{F}_{t}}$ is the
Gelfand-Pettis Integral of the map $\mathbf{x}\to
h(\mathbf{x},\mathbf{x}^{(Kt+l)}_{t})$ with respect to
$\mu_{t}|\mathcal{F}_{t}$. Hence, the following holds

$\displaystyle\mathbb{E}_{\mathbf{x}^{(Kt+l)}_{t}}\left[\left\langle
h(\cdot,\mathbf{x}^{(Kt+l)}_{t}),f\right\rangle_{\mathcal{H}}\biggr{|}\mathcal{F}_{t}\right]=\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},f\right\rangle_{\mathcal{H}}$

In particular, setting $f=h_{\mu_{t}|\mathcal{F}_{t}}$ and using linearity of
expectation, we conclude,

$\displaystyle\mathbb{E}\left[\left\langle
g_{t},h_{\mu_{t}|\mathcal{F}_{t}}\right\rangle_{\mathcal{H}}|\mathcal{F}_{t}\right]$
$\displaystyle=\frac{1}{K}\sum_{l=0}^{K-1}\mathbb{E}_{\mathbf{x}^{(Kt+l)}_{t}}\left[\left\langle
h(\cdot,\mathbf{x}^{(Kt+l)}_{t}),h_{\mu_{t}|\mathcal{F}_{t}}\right\rangle_{\mathcal{H}}\biggr{|}\mathcal{F}_{t}\right]$
$\displaystyle=\|h_{\mu_{t}|\mathcal{F}_{t}}\|_{\mathcal{H}}^{2}$

To control $\mathbb{E}[\|g_{t}\|^{2}_{\mathcal{H}}|\mathcal{F}_{t}]$, we note
that,

$\displaystyle\|g_{t}\|^{2}_{\mathcal{H}}$
$\displaystyle=\frac{1}{K^{2}}\sum_{l_{1},l_{2}=0}^{K-1}\left\langle
h(\cdot,\mathbf{x}^{(Kt+l_{1})}_{t}),h(\cdot,\mathbf{x}^{(Kt+l_{2})}_{t})\right\rangle_{\mathcal{H}}$
$\displaystyle=\frac{1}{K^{2}}\sum_{l=0}^{K-1}\|h(\cdot,\mathbf{x}^{(Kt+l)}_{t})\|^{2}_{\mathcal{H}}+\sum_{0\leq
l_{1}\neq l_{2}\leq K-1}\left\langle
h(\cdot,\mathbf{x}^{(Kt+l_{1})}_{t}),h(\cdot,\mathbf{x}^{(Kt+l_{2})}_{t})\right\rangle_{\mathcal{H}}$
$\displaystyle\leq\frac{\xi^{2}}{K}+\sum_{0\leq l_{1}\neq l_{2}\leq
K-1}\left\langle
h(\cdot,\mathbf{x}^{(Kt+l_{1})}_{t}),h(\cdot,\mathbf{x}^{(Kt+l_{2})}_{t})\right\rangle_{\mathcal{H}}$

where the last inequality uses the fact that
$\|h(\cdot,\mathbf{x}^{(Kt+l)}_{t})\|_{\mathcal{H}}\leq\xi$ almost surely as
per Lemma 9.

To control the off-diagonal terms, let $i=Kt+l_{1}$ and $j=Kt+l_{2}$ for any
arbitrary $l_{1},l_{2}$ with $0\leq l_{1}\neq l_{2}\leq K-1$. Conditioned on
$\mathcal{F}_{t}$, $\mathbf{x}^{(i)}_{t}$ and $\mathbf{x}^{(j)}_{t}$ are i.i.d
samples from $\mu_{t}|\mathcal{F}_{t}$. Thus, by Lemma 8 and Fubini’s Theorem,

$\displaystyle\mathbb{E}_{\mathbf{x}^{(i)}_{t},\mathbf{x}^{(j)}_{t}}\left[\left\langle
h(\cdot,\mathbf{x}^{(i)}_{t}),h(\cdot,\mathbf{x}^{(j)}_{t})\right\rangle_{\mathcal{H}}\bigr{|}\mathcal{F}_{t}\right]$
$\displaystyle=\mathbb{E}_{\mathbf{x}^{(i)}_{t}}\left[\mathbb{E}_{\mathbf{x}^{(j)}_{t}}\left[\left\langle
h(\cdot,\mathbf{x}^{(i)}_{t}),h(\cdot,\mathbf{x}^{(j)}_{t})\right\rangle_{\mathcal{H}}\bigr{|}\right]\right]$
$\displaystyle=\mathbb{E}_{\mathbf{x}^{(i)}_{t}}\left[\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},h(\cdot,\mathbf{x}^{(i)}_{t})\right\rangle_{\mathcal{H}}\bigr{|}\mathcal{F}_{t}\right]$
$\displaystyle=\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}$

Thus, we conclude that,

$\displaystyle\mathbb{E}\left[\|g_{t}\|^{2}_{\mathcal{H}}|\mathcal{F}_{t}\right]\leq\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}+\frac{\xi^{2}}{K}$

∎

###  C.4 Proof of Theorem 1

###### Proof.

Let $\xi=\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$
where $\zeta_{0},\dots,\zeta_{3}$ are as defined in Lemma 9. Since
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$, $\|g_{t}\|_{\mathcal{H}}\leq\xi$ holds
almost surely as per Lemma 9.

Since $\gamma\xi\leq\nicefrac{{1}}{{2B}}$, Lemma 1 ensures that the following
holds almost surely

$\displaystyle\mathsf{KL}\left(\mu_{t+1}|\mathcal{F}_{t+1}\boldsymbol{||}\pi^{\star}\right)\leq\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\gamma\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},g_{t}\right\rangle_{\mathcal{H}}+\frac{\gamma^{2}(4+L)B}{2}\|g_{t}\|^{2}_{\mathcal{H}}$

Taking conditional expectations w.r.t $\mathcal{F}_{t}$ on both sides and
applying Lemma 3, we obtain,

$\displaystyle\mathbb{E}\left[\mathsf{KL}\left(\mu_{t+1}|\mathcal{F}_{t+1}\boldsymbol{||}\pi^{\star}\right)|\mathcal{F}_{t}\right]$
$\displaystyle\leq\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\gamma\left(1-\frac{\gamma(4+L)B}{2}\right)\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}+\frac{\gamma^{2}(4+L)B\xi^{2}}{2K}$
$\displaystyle\leq\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\frac{\gamma}{2}\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}+\frac{\gamma^{2}(4+L)B\xi^{2}}{2K}$
$\displaystyle=\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\frac{\gamma}{2}\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}+\frac{\gamma^{2}(4+L)B\xi^{2}}{2K}$

where the second inequality uses the fact that
$\gamma\leq\nicefrac{{1}}{{(4+L)B}}$. Taking expectations on both sides and
rearranging,

$\displaystyle\frac{\gamma}{2}\mathbb{E}\left[\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}\right]\leq\mathbb{E}\left[\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\mathsf{KL}\left(\mu_{t+1}|\mathcal{F}_{t+1}\boldsymbol{||}\pi^{\star}\right)\right]+\frac{\gamma^{2}(4+L)B\xi^{2}}{2K}$

Telescoping and averaging, we conclude,

$\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}\right]\leq\frac{2\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{\gamma(4+L)B\xi^{2}}{K}$ (20)

Now, recall from the proof of Lemma 1 in Section C.1 that for any $t\in[T]$
and $l\geq Kt$, $\mathbf{x}^{(l)}_{t}$ depends only on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0},\mathbf{x}^{(l)}_{0}$,
i.e., there exists a deterministic measurable function $H_{t}$ such that
$\mathbf{x}^{(l)}_{t}=H_{t}(\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0},\mathbf{x}^{(l)}_{0})$
holds almost surely. We note that the output
$\mathbf{Y}=(\mathbf{y}^{(0)},\dots,\mathbf{y}^{(n-1)})$ satisfies
$\mathbf{y}^{(l)}=\mathbf{x}^{(KT+l)}_{S}\ \forall\ l\in(n)$, where
$S\sim\mathsf{Uniform}((T))$ is sampled independently of everything else.

Thus, we infer that $\mathbf{y}^{(l)}$ depends only on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S,\mathbf{x}^{(KT+l)}_{0}$,
i.e., there exists a deterministic measurable function $G$ such that
$\mathbf{y}^{(l)}=G(\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S,\mathbf{x}^{(KT+l)}_{0})$
for every $l\in(n)$. Since
$\mathbf{x}^{(KT)}_{0},\dots,\mathbf{x}^{(KT+n-1)}_{0}\
\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\ \mu_{0}$, we infer that
$\mathbf{y}^{(0)},\dots,\mathbf{y}^{(n-1)}$ are i.i.d when conditioned on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S$.

We now show that, when conditioned on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S$, $\mathbf{y}^{(l)}$ is
distributed as $\bar{\mu}$, where $\bar{\mu}$ is the probability kernel
defined as
$\bar{\mu}(\cdot;\mathbf{x}^{(0)}_{0}=x_{0},\dots,\mathbf{x}^{(KT-1)}_{0}=\mathbf{x}_{KT-1},S=s)\coloneqq\mu_{s}(\cdot,\mathbf{x}^{(0)}_{0}=\mathbf{x}_{0},\dots,\mathbf{x}^{(Ks-1)}_{0}=\mathbf{x}_{Ks-1})$.
For any arbitrary fixed $s\in(T)$, note that, under the event $S=s$,
$\mathbf{y}^{(l)}=\mathbf{x}^{(KT+l)}_{s}$ for every $l\in(n)$. Thus, for any
Borel measurable set $A\subseteq\mathbb{R}^{d}$, $\\{\mathbf{y}^{(l)}\in
A\\}\cap\\{S=s\\}=\\{\mathbf{x}^{(KT+l)}_{s}\in A\\}\cap\\{S=s\\}$. For the
sake of clarity, we denote the conditioning
$\mathbf{x}^{(0)}_{0}=\mathbf{x}_{0},\mathbf{x}^{(KT-1)}_{0}=\mathbf{x}_{KT-1}$
as $\mathcal{C}$, only in this proof. Since $S$ is independent of
$\mathbf{x}^{(l)}_{t}$ for every $t\in(T+1),l\in(KT+n)$, we infer the
following:

$\displaystyle\mathbb{P}\left(\\{\mathbf{y}^{(l)}\in
A\\}\bigr{|}\mathcal{C},S=s\right)$
$\displaystyle=\frac{\mathbb{P}\left(\\{\mathbf{y}^{(l)}\in
A\\}\cap\\{S=s\\}\bigr{|}\mathcal{C}\right)}{\mathbb{P}\left(S=s\right)}$
$\displaystyle=T\mathbb{P}\left(\\{\mathbf{x}^{(KT+l)}_{s}\in
A\\}\cap\\{S=s\\}\bigr{|}\mathcal{C}\right)$
$\displaystyle=T\mathbb{P}\left(\\{S=s\\}\right)\mathbb{P}\left(\\{\mathbf{x}^{(KT+l)}\in
A\\}\bigr{|}\mathcal{C}\right)$
$\displaystyle=\mathbb{P}\left(\\{\mathbf{x}^{(KT+l)}_{s}\in
A\\}\bigr{|}\mathcal{C}\right)$

As discussed above, $\mathbf{x}^{(KT+l)}_{s}$ depends only on
$\mathbf{x}^{(0)}_{0},\mathbf{x}^{(Ks-1)}_{0},\mathbf{x}^{(KT+l)}_{0}$. It
follows that $\mathbb{P}\left(\\{\mathbf{x}^{(KT+l)}_{s}\in
A\\}|\mathcal{C}\right)=\mu_{s}(A;\mathbf{x}^{(0)}_{0}=x_{0},\dots,\mathbf{x}^{(Ks-1)}_{0}=x_{Ks-1})$
and,

$\displaystyle\mathbb{P}\left(\\{\mathbf{y}^{(l)}\in
A\\}\bigr{|}\mathcal{C},S=s\right)$
$\displaystyle=\mu_{s}(A;\mathbf{x}^{(0)}_{0}=x_{0},\dots,\mathbf{x}^{(Ks-1)}_{0}=x_{Ks-1})$
$\displaystyle=\bar{\mu}(A;\mathbf{x}^{(0)}_{0}=x_{0},\dots,\mathbf{x}^{(KT-1)}_{0}=x_{Kt-1},S=s)$

Thus, $\mathbf{y}^{(0)},\dots,\mathbf{y}^{(n-1)}$ are i.i.d samples from
$\bar{\mu}$ when conditioned on
$\mathbf{x}^{(0)}_{0},\dots\mathbf{x}^{(KT-1)}_{0},S$.

We now obtain an upper bound on the expected squared KSD between $\bar{\mu}$
and $\pi^{\star}$. We recall from the proof of Lemma 1 in Section C.1 that,
for any $t\in(T+1)$, conditioned on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0}$,
$(\mathbf{x}^{(l)}_{t})_{l\geq t}$ are i.i.d samples from
$\mu_{t}|\mathcal{F}_{t}$ where
$\mu_{t}|\mathcal{F}_{t}\coloneqq\mu_{t}(\cdot;\mathbf{x}^{(0)}_{0},\mathbf{x}^{(Kt-1)}_{0})$.
Hence, from (20), we conclude that,

$\displaystyle\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\bar{\mu}(\cdot;(\mathbf{x}^{(l)}_{0})_{l\in(KT)},S)||\pi^{\star})^{2}]$
$\displaystyle=\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[\mathbb{E}\left[\mathsf{KSD}_{\pi^{\star}}(\bar{\mu}(\cdot;(\mathbf{x}^{(l)}_{0})_{l\in(KT)},S=t)||\pi^{\star})^{2}\bigr{|}(\mathbf{x}^{(l)}_{0})_{l\in(KT)}\right]\right]$
$\displaystyle=\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[\mathsf{KSD}_{\pi^{\star}}(\mu_{t}(\cdot;\mathbf{x}^{(0)}_{0},\cdot,\mathbf{x}^{(Kt-1)}_{0})||\pi^{\star})^{2}\right]$
$\displaystyle=\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}\right]$
$\displaystyle\leq\frac{2\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{\gamma(4+L)B\xi^{2}}{K}$

where we use the fact that $S\sim\mathsf{Uniform}((T))$ is sampled independent
of everything else. ∎

###  C.5 Proof of Theorem 2

###### Proof.

Following the same steps as Theorem 1, we note that the following holds almost
surely.

$\displaystyle\mathsf{KL}\left(\mu_{t+1}|\mathcal{F}_{t+1}\boldsymbol{||}\pi^{\star}\right)$
$\displaystyle\leq\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\gamma\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},g_{t}\right\rangle_{\mathcal{H}}+\frac{\gamma^{2}(4+L)B}{2}\|g_{t}\|^{2}_{\mathcal{H}}$
$\displaystyle\leq\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\frac{\gamma}{2}\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}+\gamma\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},h_{\mu_{t}|\mathcal{F}_{t}}-g_{t}\right\rangle_{\mathcal{H}}$
$\displaystyle+\frac{\gamma^{2}(4+L)B\xi^{2}}{2K}+\frac{\gamma^{2}(4+L)B}{2}\left[\|g_{t}\|^{2}_{\mathcal{H}}-\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}-\frac{\xi^{2}}{K}\right]$
(21)

where the last inequality uses the fact that
$\gamma\leq\nicefrac{{1}}{{(4+L)B}}$. We now define $\Delta^{(l)}_{t}$,
$\Delta_{t}$ and $r_{t}$ for $l\in(K),\ t\in(T)$ as follows:

$\displaystyle\Delta^{(l)}_{t}$ $\displaystyle=\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},h_{\mu_{t}|\mathcal{F}_{t}}-h(\cdot,\mathbf{x}^{(Kt+l)}_{t})\right\rangle_{\mathcal{H}}$
$\displaystyle\Delta_{t}$
$\displaystyle=\frac{1}{K}\sum_{l=0}^{K-1}\Delta^{(l)}_{t}=\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},h_{\mu_{t}|\mathcal{F}_{t}}-g_{t}\right\rangle_{\mathcal{H}}$
$\displaystyle r_{t}$
$\displaystyle=\|g_{t}\|^{2}_{\mathcal{H}}-\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}-\frac{\xi^{2}}{K}$

Substituting the above into (C.5), we obtain the following:

$\displaystyle\mathsf{KL}\left(\mu_{t+1}|\mathcal{F}_{t+1}\boldsymbol{||}\pi^{\star}\right)$
$\displaystyle\leq\mathsf{KL}\left(\mu_{t}|\mathcal{F}_{t}\boldsymbol{||}\pi^{\star}\right)-\frac{\gamma}{2}\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}+\gamma\Delta_{t}+\frac{\gamma^{2}(4+L)B\xi^{2}}{2K}+\frac{\gamma^{2}(4+L)Br_{t}}{2}$

Telescoping and averaging both sides, and using
$\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}=\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}$,
we obtain the following:

$\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}$
$\displaystyle\leq\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma(4+L)B\xi^{2}}{K}$
$\displaystyle+\frac{4}{T}\sum_{t=0}^{T-1}\left(\Delta_{t}-\frac{\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}}{4}\right)+\frac{2\gamma(4+L)B}{T}\sum_{t=0}^{T-1}r_{t}$
(22)

We note that the first two terms are the same as that of the in-expectation
guarantee for VP-SVGD in Theorem 1. The third and fourth term are random
quantities that vanish in expectation. The remainder of our analysis upper
bounds them with high probability.

We begin by deriving a high probability upper bound for the fourth term in
(C.5). To this end, we note that, since $\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$,
$\|h(\cdot,\mathbf{x}^{(Kt+l)}_{t})\|_{\mathcal{H}}\leq\xi$ for any $t\in(T),\
l\in(K)$ as per Lemma 9. Furthermore, since
$\mathbb{E}[h(\cdot,\mathbf{x}^{(Kt+l)}_{t}|\mathcal{F}_{t}]=h_{\mu_{t}|\mathcal{F}_{t}}$
(both pointwise and in the sense of the Gelfand-Pettis integral, see proof of
Lemma 3 in Appendix C.3), it follows by Jensen’s inequality that
$\|h_{\mu_{t}|\mathcal{F}_{t}}\|_{\mathcal{H}}\leq\xi$. This further implies
that $|r_{t}|\leq 3\xi^{2}$. Moreover, $r_{t}$ is $\mathcal{F}_{t+1}$
measurable (as $g_{t}$ is an $\mathcal{F}_{t+1}$ measurable random function)
with $\mathbb{E}[r_{t}|\mathcal{F}_{t}]\leq 0$ (as per Lemma 3)

Thus, $S_{t}=\sum_{s=0}^{t-1}r_{t}$ is an $\mathcal{F}$-adapted
supermartingale difference sequence with bounded increments. Thus, by the
Hoeffding-Azuma inequality, we conclude that the following holds with
probability at least $1-\nicefrac{{\delta}}{{2}}$

$\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}r_{t}\leq
6\xi^{2}\sqrt{\frac{\log(\nicefrac{{2}}{{\delta}})}{T}}$ (23)

We now proceed to control the third term in (C.5). Recall from the proof of
Theorem 1 in Appendix C.4, that, for any fixed $t\in(T)$,
$(\mathbf{x}^{(l)}_{t})_{l\in(KT)}$ are i.i.d when conditioned on
$\mathcal{F}_{t}$. As discussed above,
$\mathbb{E}[h(\cdot,\mathbf{x}^{(Kt+l)}_{t})]=h_{\mu_{t}|\mathcal{F}_{t}}$ in
the sense of the Gelfand-Pettis integral, implying
$\mathbb{E}[\Delta^{(l)}_{t}]=0$. Moreover, $|\Delta^{(l)}\|_{t}\leq
2\xi\|h_{\mu_{t}|\mathcal{F}_{t}}\|$. Thus, when conditioned on
$\mathcal{F}_{t}$, $\Delta^{(l)}_{t}$ are independent zero-mean bounded random
variables. Hence, we conclude the following by Hoeffding’s Lemma

$\displaystyle\mathbb{E}\left[e^{\theta\Delta_{t}}|\mathcal{F}_{t}\right]$
$\displaystyle\leq\prod_{l=0}^{K-1}\mathbb{E}[e^{\tfrac{\theta\Delta^{(l)}_{t}}{K}}|\mathcal{F}_{t}]\leq
e^{\tfrac{2\theta^{2}\xi^{2}}{K}\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}},\quad\forall\
\theta\in\mathbb{R}$ (24)

We now define the sequence $M_{t}$ as follows, where
$\lambda=\nicefrac{{K}}{{8\xi^{2}}}$

$\displaystyle
M_{t}=\exp(\sum_{s=0}^{t-1}\lambda\Delta_{s}-\tfrac{\lambda}{4}\|h_{\mu_{s}|\mathcal{F}_{s}}\|^{2}_{\mathcal{H}})$

Since $g_{t}$ is $\mathcal{F}_{t+1}$ measurable, so is $\Delta_{t}$, which
implies $M_{t}$ is $\mathcal{F}_{t+1}$ measurable. Furthermore,

$\displaystyle\mathbb{E}[M_{t}|\mathcal{F}_{t}]$
$\displaystyle=M_{t-1}e^{-\tfrac{\lambda}{4}\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}}\mathbb{E}[e^{\lambda\Delta_{t}}|\mathcal{F}_{t}]$
$\displaystyle\leq
M_{t-1}e^{(-\tfrac{\lambda}{4}+\tfrac{2\lambda^{2}\xi^{2}}{K})\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}}\leq
M_{t-1}$

Thus, $M_{t}$ is an $\mathcal{F}$-adapted supermartingale sequence. Following
the same steps, we conclude $E[M_{1}]\leq 1$, which implies
$\mathbb{E}[M_{T}]\leq\mathbb{E}[M_{1}]\leq 1$. Thus, from Markov’s Inequality

$\displaystyle\mathbb{P}\left[\sum_{t=0}^{T-1}\Delta_{t}-\tfrac{1}{4}\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}>x\right]\leq
e^{-\lambda x}\mathbb{E}[M_{T}]\leq e^{-\lambda x}$

Hence, the following holds with probability at least
$1-\nicefrac{{\delta}}{{2}}$.

$\displaystyle\sum_{t=0}^{T-1}\Delta_{t}-\tfrac{1}{4}\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}\leq\frac{8\xi^{2}}{K}\log(\nicefrac{{2}}{{\delta}})$
(25)

Substituting (24) and (25) into (20) and taking a union bound, we conclude
that the following holds with probability at least $1-\delta$:

$\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}$
$\displaystyle\leq\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma(4+L)B\xi^{2}}{K}$
$\displaystyle+\frac{32\xi^{2}\log(\nicefrac{{2}}{{\delta}})}{KT}+12\gamma(4+L)B\xi^{2}\sqrt{\frac{\log(\nicefrac{{2}}{{\delta}})}{T}}$
(26)

Recall from the proof of Theorem 1 in Appendix C.4 that the outputs
$(\mathbf{y}^{(l)})_{l\in(n)}$ of VP-SVGD are i.i.d samples from the random
measure
$\bar{\mu}(\cdot;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(0)}_{KT-1},S)$ when
conditioned on $\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(0)}_{KT-1},S$.
Furthermore, when conditioned on $S=t$,
$\bar{\mu}(\cdot;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(0)}_{KT-1},S=t)=\mu_{t}|\mathcal{F}_{t}$.
Thus, from (C.5), we conclude that, upon taking an expectation over
$S\sim\mathsf{Uniform}((T))$ while conditioning on the virtual particles
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0}$, the following holds with
probability at least $1-\delta$:

$\displaystyle\mathbb{E}_{S}[\mathsf{KSD}_{\pi^{\star}}(\bar{\mu}(\cdot;\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0},S)||\pi^{\star})^{2}]$
$\displaystyle\leq\frac{1}{T}\sum_{t=0}^{T-1}\mathsf{KSD}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})^{2}$
$\displaystyle\leq\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma(4+L)B\xi^{2}}{K}$
$\displaystyle+\frac{32\xi^{2}\log(\nicefrac{{2}}{{\delta}})}{KT}+12\gamma(4+L)B\xi^{2}\sqrt{\frac{\log(\nicefrac{{2}}{{\delta}})}{T}}$

∎

##  Appendix D Analysis of GB-SVGD

In this section, we present our analysis of GB-SVGD. For any $t\in(T)$, we use
$\tilde{g}_{t}$ to denote the random function
$\tilde{g}_{t}(\mathbf{x})=\tfrac{1}{K}\sum_{r\in\mathcal{K}_{t}}h(\mathbf{x},\mathbf{x}^{(r)}_{t})$
where $\mathcal{K}_{t}$ is the random batch of size $K$ sampled at time-step
$t$ of GB-SVGD.

In order to prove Theorem 3, we first establish an almost-sure iterate bound
for GB-SVGD which is similar to that of Lemma 9 for VP-SVGD.

######  Lemma 10 (Almost-Sure Iterate Bounds).

Let Assumptions 1, 2, 3, 4 and 5 be satisfied. Then, the following holds
almost surely for any $s\in\mathbb{N}\cup\\{0\\}$ and $t\in(T+1)$ whenever
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$

$\displaystyle\|\mathbf{x}^{(s)}_{t}\|$
$\displaystyle\leq\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$
$\displaystyle\|h(\cdot,\mathbf{x}^{(s)}_{t})\|_{\mathcal{H}}$
$\displaystyle\leq\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$
$\displaystyle\|\tilde{g}_{t}\|_{\mathcal{H}}$
$\displaystyle\leq\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$

where $\zeta_{0},\dots,\zeta_{3}$ are problem-dependent constants that depend
polynomially on $A_{1},A_{2},A_{3},B,d_{1},d_{2},L$ for any fixed $\alpha$.

###### Proof.

The proof of this Lemma is identical to that of Lemma 9. To this end, let
$c^{(s)}_{t}=\tfrac{1}{K}\sum_{r\in\mathcal{K}_{t}}k(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(r)}_{t})$.
Note that by Assumption 4, $c^{(s)}_{t}\geq 0$ Since
$\mathbf{x}^{(s)}_{t+1}=\mathbf{x}^{(s)}_{t}-\gamma\tilde{g}_{t}(\mathbf{x}^{(s)}_{t})$,
it follows from the smoothness of $F$ that,

$F(\mathbf{x}^{(s)}_{t+1})-F(\mathbf{x}^{(s)})\leq-\gamma\left\langle\nabla
F(\mathbf{x}^{(s)}_{t}),\tilde{g}_{t}(\mathbf{x}^{(s)}_{t})\right\rangle+\frac{\gamma^{2}L}{2}\|\tilde{g}_{t}(\mathbf{x}^{(s)}_{t})\|^{2}$
(27)

By Lemma 6, we note that,

$\displaystyle-\gamma\left\langle\nabla
F(\mathbf{x}^{(s)}_{t}),\tilde{g}_{t}(\mathbf{x}^{(s)}_{t})\right\rangle$
$\displaystyle=-\frac{\gamma}{K}\sum_{r\in\mathcal{K}_{t}}\left\langle\nabla
F(\mathbf{x}^{(s)}_{t}),h(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(r)}_{t})\right\rangle$
$\displaystyle\leq\frac{\gamma}{K}\sum_{r\in\mathcal{K}_{t}}\left[-\tfrac{1}{2}k(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(r)}_{t})\|\nabla
F(\mathbf{x}^{(s)}_{t})\|^{2}+L^{2}A_{1}+A_{3}\right]$
$\displaystyle\leq-\frac{\gamma c^{(s)}_{t}}{2}\|\nabla
F(\mathbf{x}^{(s)}_{t})\|^{2}+\gamma L^{2}A_{1}+\gamma A_{3}$ (28)

Moreover, by Jensen’s Inequality and Lemma 6

$\displaystyle\|\tilde{g}_{t}(\mathbf{x}^{(s)}_{t})\|^{2}$
$\displaystyle\leq\frac{1}{K}\sum_{r\in\mathcal{K}_{t}}\|h(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(r)}_{t})\|^{2}$
$\displaystyle\leq\frac{1}{K}\sum_{r\in\mathcal{K}_{t}}2(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}+2k(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(r)}_{t})^{2}\|F(\mathbf{x}^{(s)}_{t})\|^{2}$
$\displaystyle\leq\frac{1}{K}\sum_{r\in\mathcal{K}_{t}}2(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}+2A_{1}k(\mathbf{x}^{(s)}_{t},\mathbf{x}^{(r)}_{t})\|F(\mathbf{x}^{(s)}_{t})\|^{2}$
$\displaystyle\leq
2(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}+2A_{1}c^{(s)}_{t}\|F(\mathbf{x}^{(s)}_{t})\|^{2}$
(29)

Substituting (D) and (D) into (27), we obtain,

$\displaystyle F(\mathbf{x}^{(s)}_{t+1})-F(\mathbf{x}^{(s)}_{t})$
$\displaystyle\leq-\frac{\gamma c^{(s)}_{t}}{2}\|\nabla
F(\mathbf{x}^{(s)}_{t})\|^{2}+\gamma L^{2}A_{1}+\gamma A_{3}$
$\displaystyle+\gamma^{2}L(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}+\gamma^{2}LA_{1}c^{(s)}_{t}\|F(\mathbf{x}^{(s)}_{t})\|^{2}$
$\displaystyle\leq-\frac{\gamma c^{(s)}_{t}}{2}(1-2A_{1}L\gamma)\|\nabla
F(\mathbf{x}^{(s)}_{t})\|^{2}+\gamma A_{3}+\gamma
L^{2}A_{1}+\gamma^{2}L(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}$
$\displaystyle\leq\gamma A_{3}+\gamma
L^{2}A_{1}+\gamma^{2}L(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}$

where the last inequality uses the fact that $c^{(s)}_{t}\geq 0$ and
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$. Now, iterating through the above
inequality, we obtain the following for any $t\in[T],\
s\in\mathbb{N}\cup\\{0\\}$

$F(\mathbf{x}^{(s)}_{t})\leq F(\mathbf{x}^{(s)}_{0})+\gamma TL^{2}A_{1}+\gamma
TA_{3}+\gamma^{2}TL(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}$ (30)

Furthermore, by Assumption 1

$\displaystyle F(\mathbf{x}^{(s)}_{0})$ $\displaystyle\leq F(0)+\|\nabla
F(0)\|\|\mathbf{x}^{(s)}_{0}\|+\frac{L}{2}\|\mathbf{x}^{(s)}_{0}\|^{2}$
$\displaystyle\leq F(0)+\nicefrac{{1}}{{2}}+L\|\mathbf{x}^{(s)}_{0}\|^{2}$

Substituting the above inequality into (30), and using Assumption 2, we obtain
the following for any $t\in[T],\ s\in\mathbb{N}\cup\\{0\\}$

$\displaystyle d_{1}\|\mathbf{x}^{(s)}_{t}\|^{\alpha}-d_{2}\leq
F(\mathbf{x}^{s}_{t})$ $\displaystyle\leq
F(0)+\nicefrac{{1}}{{2}}+L\|\mathbf{x}^{(s)}_{0}\|^{2}+\gamma
TL^{2}A_{1}+\gamma TA_{3}$
$\displaystyle+\gamma^{2}TL(\nicefrac{{A_{1}L}}{{2}}+A_{2})^{2}$

Rearranging and applying Assumption 5, we obtain

$\displaystyle\|\mathbf{x}^{(s)}_{t}\|$ $\displaystyle\leq
d^{-\nicefrac{{1}}{{\alpha}}}_{1}\left[F(0)+\nicefrac{{1}}{{2}}+LR^{2}+\gamma
TL^{2}A_{1}+\gamma
TA_{3}+\gamma^{2}TL(A_{1}L+A_{2})^{2}\right]^{\nicefrac{{1}}{{\alpha}}}$
$\displaystyle\leq\tilde{\zeta}_{0}+\tilde{\zeta}_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\zeta}_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\zeta}_{3}R^{\nicefrac{{2}}{{\alpha}}}$

where $\tilde{\zeta}_{0},\dots,\tilde{\zeta}_{3}$ are constants that depend
polynomially on $L,A_{1},A_{2},A_{3},R$. We note that, since $0<\alpha\leq 2$,
the above inequality also holds for $t=0$.

Using the above inequality Lemma 6 and Assumption 1, we conclude that the
following holds almost surely for any $t\in(T+1),s\in\mathbb{N}\cup\\{0\\}$

$\displaystyle\|h(\cdot,\mathbf{x}^{(s)}_{t})\|_{\mathcal{H}}$
$\displaystyle\leq BL\|\mathbf{x}^{(s)}_{t}\|+B\sqrt{L}+B$
$\displaystyle\leq\tilde{\eta_{0}}+\tilde{\eta}_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\eta}_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\eta}_{3}R^{\nicefrac{{2}}{{\alpha}}}$

where $\tilde{\eta}_{0},\dots,\tilde{\eta}_{3}$ are constants that depend
polynomially on $L,B,A_{1},A_{2},A_{3},R$. Using the above inequality, we
conclude that the following also holds for any $t\in(T+1)$.

$\displaystyle\|\tilde{g}_{t}\|_{\mathcal{H}}\leq\tilde{\eta_{0}}+\tilde{\eta}_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\eta}_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\tilde{\eta}_{3}R^{\nicefrac{{2}}{{\alpha}}}$

Taking $\zeta_{i}=\max\\{\tilde{\zeta_{i}},\tilde{\eta}_{i}\\}$, the proof is
complete. ∎

###  D.1 Proof of Theorem 3

###### Proof.

Let $\xi=\zeta_{0}+\zeta_{1}(\gamma
T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{2}(\gamma^{2}T)^{\nicefrac{{1}}{{\alpha}}}+\zeta_{3}R^{\nicefrac{{2}}{{\alpha}}}$
where $\zeta_{0},\dots,\zeta_{3}$ are constants as described in Lemma 9 and
Lemma 10. Since the assumptions and parameter settings of Theorem 1 holds,
$\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$ and thus, by Lemma 9 and Lemma 10, the
particles output by VP-SVGD and GB-SVGD are bounded as
$\|\mathbf{y}^{(l)}\|\leq\xi$ and $\|\bar{\mathbf{y}}^{(l)}\|\leq\xi$.

Let $\mathbf{Y}=(\mathbf{y}^{(0)},\dots,\mathbf{y}^{(n-1)})$ and
$\bar{\mathbf{Y}}=(\bar{\mathbf{y}}^{(0)},\dots,\bar{\mathbf{y}}^{(n-1)})$
denote the outputs of VP-SVGD and GB-SVGD. Let
$\hat{\mu}^{(n)}=\tfrac{1}{n}\sum_{i=0}^{n-1}\delta_{\mathbf{y}^{(i)}}$ and
$\hat{\nu}^{(n)}=\tfrac{1}{n}\sum_{i=0}^{n-1}\delta_{\bar{\mathbf{y}}^{(i)}}$
be their respective empirical distributions. We shall now explicitly construct
a coupling between the inputs of VP-SVGD and GB-SVGD such that the first
$n-KT$ particles of their respective outputs are equal. This in turn will