A Convergence Theory for SVGD in the Population Limit
under Talagrand’s Inequality T1

Stein Variational Gradient Descent (SVGD) (Liu & Wang, 2016; Liu, 2017) is an alternative to the Langevin algorithm that has been applied in several contexts in Machine Learning, including Reinforcement Learning (Liu et al., 2017), sequential decision making (Zhang et al., 2018, 2019), Generative Adversarial Networks (Tao et al., 2019), Variational Auto Encoders (Pu et al., 2017), and Federated Learning (Kassab & Simeone, 2020).

The literature on theoretical properties of SVGD is scarce compared to that of Langevin algorithm, and limited to a few recent works (Korba et al., 2020; Lu et al., 2019; Duncan et al., 2019; Liu, 2017; Chewi et al., 2020; Gorham et al., 2020; Nüsken & Renger, 2021; Shi et al., 2021). In this paper, our goal is to provide a clean convergence theory for SVGD in the population limit, i.e., with an infinite number of particles.

The Machine Learning literature on the complexity of sampling from a non-logconcave target distribution has mainly focused on the Langevin algorithm. For instance,¹¹1The example of Vempala & Wibisono (2019) is taken only for illustration purpose. Many other results were obtained for Langevin algorithm, even in nonconvex cases, see above. Vempala & Wibisono (2019) showed that Langevin algorithm reaches $\varepsilon$ accuracy in terms of the Kullback-Leibler divergence after $\tilde{\Omega}(\frac{L^{2}d}{\lambda^{2}\varepsilon})$ iterations, assuming that the target distribution satisfies the logarithmic-Sobolev inequality (LSI) with constant $\lambda$. In this work, we will assume Talagrand’s inequality T1 with constant $\lambda$, which is milder than LSI with constant $\lambda$, and we will prove a complexity result in terms of another discrepancy called Kernelized Stein Discrepancy (KSD). Besides, a very recent work studies Langevin algorithm for a non-logconcave target distribution without assuming LSI and provides guarantees in terms of the Fisher information (Balasubramanian et al., 2022).

Most existing results on SVGD deal with the continuous time approximation of SVGD in the population limit, a Partial Differential Equation (PDE) representing SVGD with a vanishing step size and an infinite number of particles (Lu et al., 2019; Duncan et al., 2019; Liu, 2017; Nüsken & Renger, 2021; Chewi et al., 2020). In particular, Duncan et al. (2019) propose a Stein logarithmic Sobolev inequality that implies the linear convergence of this PDE. However, it is not yet understood when Stein logarithmic Sobolev inequality holds. Besides, Chewi et al. (2020) showed that the Wasserstein gradient flow of the chi-squared divergence can be seen as an approximation of that PDE, and showed linear convergence of the Wasserstein gradient flow of the chi-squared under Poincaré inequality. Other results, such as those of Lu et al. (2019); Liu (2017); Nüsken & Renger (2021), include asymptotic convergence properties of the PDE, but do not include convergence rates. In this paper, we will prove convergence rates for SVGD in discrete time.

We consider SVGD in the population limit, similarly to concurrent works such as Liu (2017); Korba et al. (2020); Gorham et al. (2020). Our paper intends to provide a clean analysis of SVGD, a problem stated in Liu (2017, Conclusion). To this end, we do not make any assumptions on the trajectory of the algorithm. Instead, our key assumption is that the target distribution $\pi$ satisfies T1, the mildest of the Talagrand’s inequalities, which holds under a mild assumption on the tails of the distribution; see Villani (2008, Theorem 22.10). Moreover, T1 is implied, for example, by the logarithmic Sobolev inequality (Villani, 2008, Theorem 22.17), with the same constant $\lambda$.

Although sampling algorithms are meant to be applied on high-dimensional problems, the question of the dependence of the complexity of SVGD in $d$ has not been studied in concurrent works, nor has been studied the generic weak convergence of SVGD under verifiable assumptions, to our knowledge. Assuming that the T1 inequality holds, we provide

• •

  a generic weak convergence result for SVGD (actually our result is a bit stronger: convergence holds in 1-Wasserstein distance),

• •

  a complexity bound for SVGD in terms of the dimension $d$ and the desired accuracy $\varepsilon$, under verifiable assumptions (i.e., assumptions that do not depend on the trajectory of the algorithm): $\tilde{\Omega}\left(\frac{L{d}^{3/2}}{\lambda^{1/2}\varepsilon}\right)$ iterations suffice to obtain a sample $\mu$ such that $\mathop{\mathrm{KSD}}\nolimits^{2}(\mu|\pi)<\varepsilon$, where $L$ is the smoothness constant of $F$ and $\lambda$ the constant in T1 inequality.

Note that these results hold without assuming $F$ convex. In particular, in the population limit, SVGD applied to non-logconcave target distributions satisfying T1 converges to the target distribution.

The remainder of the paper is organized as follows. In Section 2 we introduce the necessary mathematical and notational background on optimal transport, reproducing kernel Hilbert spaces and SVGD in order to be able to describe and explain our results. Section 3 is devoted to the development of our theory. Finally, in Section 4 we formulate three corollaries of our key result, capturing weak convergence and complexity estimates for SVGD. Technical proofs are postponed to the Appendix.

For any Hilbert space $H$, we denote by $\langle\cdot,\cdot\rangle_{H}$ the inner product of $H$ and by $\|\cdot\|_{H}$ its norm.

We denote by $C_{0}({\mathcal{X}})$ the set of continuous functions from ${\mathcal{X}}$ to ${\mathbb{R}}$ vanishing at infinity and by $C^{1}({\mathcal{X}},{\mathcal{Y}})$ the set of continuously differentiable functions from ${\mathcal{X}}$ to a Hilbert space ${\mathcal{Y}}$. Given $\phi\in C^{1}({\mathcal{X}},{\mathbb{R}})$, its gradient is denoted by $\nabla\phi$, and if $\phi\in C^{1}({\mathcal{X}},{\mathcal{X}})$, the Jacobian of $\phi$ is denoted by $J\phi$. For every $x\in{\mathcal{X}}$, $J\phi(x)$ can be seen as a $d\times d$ matrix. The trace of the Jacobian, also called divergence, is denoted by $\mathop{\mathrm{div}}\nolimits\phi$.

For any $d\times d$ matrix $A$, $\|A\|_{\mathop{\mathrm{HS}}\nolimits}$ denotes the Hilbert Schmidt norm of $A$ and $\|A\|_{\mathop{\mathrm{op}}\nolimits}$ the operator norm of $A$ viewed as a linear operator $A:{\mathcal{X}}\to{\mathcal{X}}$ (where ${\mathcal{X}}$ is endowed with the standard Euclidean inner product). Finally, $\delta_{x}$ is the Dirac measure at $x\in{\mathcal{X}}$.

Consider $p\geq 1$. We denote by ${\mathcal{P}}_{p}({\mathcal{X}})$ the set of Borel probability measures $\mu$ over ${\mathcal{X}}$ with finite $p^{\text{th}}$ moment: $\int\|x\|^{p}d\mu(x)<\infty$. We denote by $L^{p}(\mu)$ the set of measurable functions $f:{\mathcal{X}}\to{\mathcal{X}}$ such that $\int\|f\|^{p}d\mu<\infty$. Note that the identity map $I$ of ${\mathcal{X}}$ satisfies $I\in L^{p}(\mu)$ if $\mu\in{\mathcal{P}}_{p}({\mathcal{X}})$. Moreover, denoting the image (or pushforward) measure of $\mu$ by a map $T$ as $T\#\mu$, we have that if $\mu\in{\mathcal{P}}_{p}({\mathcal{X}})$ and $T\in L^{p}(\mu)$ then $T\#\mu\in{\mathcal{P}}_{p}({\mathcal{X}})$ using the transfer lemma.

For every $\mu,\nu\in{\mathcal{P}}_{p}({\mathcal{X}})$, the $p$-Wasserstein distance between $\mu$ and $\nu$ is defined by

where ${\mathcal{S}}(\mu,\nu)$ is the set of couplings between $\mu$ and $\nu$, i.e., the set of nonnegative measures over ${\mathcal{X}}^{2}$ such that $P\#s=\mu$ (resp. $Q\#s=\nu$) where $P:(x,y)\mapsto x$ (resp. $Q:(x,y)\mapsto y$) denotes the projection onto the first (resp. the second) component. The $p$-Wasserstein distance is a metric over ${\mathcal{P}}_{p}({\mathcal{X}})$. The metric space $({\mathcal{P}}_{2}({\mathcal{X}}),W_{2})$ is called the Wasserstein space.

In this paper, we consider a target probability distribution $\pi$ proportional to $\exp(-F)$, where $F$ satisfies the following.

Moreover, using $\int\exp(-F(x))dx<\infty$, $F$ admits a stationary point.

To specify the dependence in the dimension of our complexity bounds, we will initialize the algorithm from a Gaussian distribution centered at a stationary point. Such a stationary point can be found by gradient descent on $F$ for instance.

The task of sampling from $\pi$ can be formalized as an optimization problem. Indeed, define the Kullback-Leibler ($\mathop{\mathrm{KL}}\nolimits$) divergence as

if $\mu$ admits the density $\frac{d\mu}{d\pi}$ with respect to $\pi$, and $\mathop{\mathrm{KL}}\nolimits(\mu|\pi)\coloneqq+\infty$ else. Then, $\mathop{\mathrm{KL}}\nolimits(\mu|\pi)\geq 0$ and $\mathop{\mathrm{KL}}\nolimits(\mu|\pi)=0$ if and only if $\mu=\pi$. Therefore, assuming $\pi\in{\mathcal{P}}_{2}({\mathcal{X}})$, the optimization problem

where

admits a unique solution: the distribution $\pi$. We will see in Section 3 that SVGD can be seen as a gradient descent algorithm to solve (4).

Indeed, the Wasserstein space can be endowed with a differential structure. In particular, when it is well defined, the Wasserstein gradient of the functional ${\mathcal{F}}$ denoted by $\nabla_{W}{\mathcal{F}}(\mu)$ is an element of $L^{2}(\mu)$ and satisfies $\nabla_{W}{\mathcal{F}}(\mu)=\nabla\log\left(\frac{d\mu}{d\pi}\right)$.

We consider a kernel $k$ associated to a Reproducing Kernel Hilbert Space (RKHS) denoted by ${\mathcal{H}}_{0}$. We denote by $\Phi:{\mathcal{X}}\to{\mathcal{H}}_{0}$ the so-called feature map $\Phi:x\mapsto k(\cdot,x)$. The product space ${\mathcal{H}}_{0}^{d}$ is also a Hilbert space denoted ${\mathcal{H}}\coloneqq{\mathcal{H}}_{0}^{d}$. We make the following assumption on the kernel $k$.

2.6 is satisfied by the Gaussian kernel for example, with $B$ independent of $d$ using a scaling argument. 2.6 states that $\Phi:{\mathcal{X}}\to{\mathcal{H}}_{0}$ is bounded, Lipschitz and $C^{1}$. This is satisfied by many classical kernels used in practice. Note that $k(x,x)=\|\Phi(x)\|_{{\mathcal{H}}_{0}}^{2}$ and that $\mathop{\mathrm{div}}\nolimits_{1}\nabla_{2}k(x,a)=\langle\nabla\Phi(x),\nabla\Phi(a)\rangle_{{\mathcal{H}}}$ (in particular, $\mathop{\mathrm{div}}\nolimits_{1}\nabla_{2}k(x,x)=\|\nabla\Phi(x)\|_{{\mathcal{H}}}^{2}$). Hence, $\nabla\Phi$ is continuous iff $x\mapsto\mathop{\mathrm{div}}\nolimits_{1}\nabla_{2}k(x,x)$ and $x\mapsto\mathop{\mathrm{div}}\nolimits_{1}\nabla_{2}k(x,a)$ are continuous for every $a\in{\mathcal{X}}$.

Under Assumption 2.6, ${\mathcal{H}}\subset L^{2}(\mu)$ for every probability distribution on ${\mathcal{X}}$, and the inclusion map $\iota_{\mu}:{\mathcal{H}}\to L^{2}(\mu)$ is continuous. We denote by $P_{\mu}:L^{2}(\mu)\to{\mathcal{H}}$ its adjoint defined by the relation: for every $f\in L^{2}(\mu)$, $g\in{\mathcal{H}}$,

Then, $P_{\mu}$ can be expressed as a convolution with $k$ (Carmeli et al., 2010, Proposition 3):

or $P_{\mu}f=\int\Phi(y)f(y)d\mu(y)$ where the integral converges in norm.

We start by stating a fundamental inequality satisfied by ${\mathcal{F}}$ for any update of the form

where $g\in{\mathcal{H}}$.

Inequality (11) is a property of the functional ${\mathcal{F}}$, and not a property of the SVGD algorithm. Inequality (11) plays the role of a Taylor inequality for the functional ${\mathcal{F}}$, where $h_{\mu_{n}}$ is the Wasserstein gradient of ${\mathcal{F}}$ at $\mu_{n}$ under the metric induced by ${\mathcal{H}}$. Proposition 3.1 is a slight generalization of Korba et al. (2020, Proposition 5), and is not our main contribution, therefore we only sketch its proof in the Appendix.

Applying recursively the Taylor inequality Proposition 3.1 with $g=h_{\mu_{n}}$, we obtain the following descent property for SVGD, which is our main theoretical result. The proof of this result can be found in the Appendix.

If ${\mathcal{F}}(\mu_{0})<\infty$, then, using Theorem 3.2, $({\mathcal{F}}(\mu_{n}))_{n}$ is nonincreasing and $\mu_{n}$ has a density w.r.t. Lebesgue measure for every $n$ (since ${\mathcal{F}}(\mu_{n})<\infty$).

In the language of the gradient descent algorithms, Theorem 3.2 is called a descent property. It can be seen as a discrete time analogue of dissipation properties obtained for the PDE modeling SVGD in continuous time in the population limit (Duncan et al., 2019; Korba et al., 2020).

Unlike Korba et al. (2020, Proposition 5) and Liu (2017, Theorem 3.3), we do not assume that $\sup_{n}\mathop{\mathrm{KSD}}\nolimits(\mu_{n}|\pi)<\infty$ or that $\gamma\leq\mathop{\mathrm{KSD}}\nolimits(\mu_{n}|\pi)^{-1}$ to obtain our descent property. The step size $\gamma$ is bounded by a constant. Iterating Theorem 3.2, we obtain convergence results as corollaries in the next section.

We now show that Theorem 3.2 implies weak convergence and convergence in $W_{1}$.

In summary, under T1 and some smoothness assumptions but without convexity of the potential, SVGD in the population limit converges to the target distribution.

One can be surprised to see that SVGD converges without convexity assumption on $F$, but this is actually natural if one thinks about the gradient descent interpretation of SVGD. Indeed, SVGD in the population limit is a gradient descent on the KL divergence, which is

• •

  "smooth" if we restrict the descent directions to a RKHS (i.e., it satisfies a Taylor inequality Proposition 3.1),

• •

  coercive (i.e., sublevel sets are tight) Dupuis & Ellis (2011, Lemma 1.4.3),

• •

  and has a single stationary point which is its global minimizer (the KSD is the norm of the gradient of KL in our interpretation, and the KSD is equal to zero only at the optimum).

One can show that, over ${\mathbb{R}}^{d}$, gradient descent applied to a smooth coercive function with a single stationary point converges to the global minimizer. The situation here is similar.

Next, we provide a ${\mathcal{O}}(1/n)$ convergence rate for the empirical mean of the iterates $\mu_{n}$ in terms of the squared KSD. This result is obtained from our descent lemma (Theorem 3.2).

Note that this convergence rate is given in terms of the uniform mixture of $\mu_{0},\ldots,\mu_{n-1}$. Similar mixtures appear in the analysis of Langevin algorithm (see e.g. Durmus et al. (2019)). Note also that the convergence rate in Corollary 4.2 is similar to the convergence rate of the squared norm of the gradient in the gradient descent algorithm applied to a smooth function (Nesterov, 2013).

From the last result, we can characterize the iteration complexity of SVGD.

To our knowledge, Corollary 4.3 provides the first dimension-dependent complexity result for SVGD. Its proof can be found in the appendix. The dependence of the T1 constant $\lambda$ in the dimension $d$ is subject to active research in optimal transport theory (Villani, 2008, Remark 22.11) and is out of the scope of this paper. Yet, using Villani (2008, Theorem 22.10, Equation 22.16), $\lambda$ can be taken as

Note that the output $\mu$ of the algorithm is a mixture of the iterates: $\mu=\bar{\mu}_{n}$. Besides, optimizing the complexity over $\alpha$ leads to involved calculations that do not change the overall complexity. To see this, note that the larger the step size $\gamma$, the smaller the complexity. But, even if the step size $\gamma=\min(\frac{2}{BL},1/K)$ were allowed, the overall complexity would be the same.