0.0.1 Introduction↩︎

Bayesian methods are widely implemented in various inference tasks that emerge in computational statistics [2], [3], machine learning [4], [5], inverse problems [6], [7] and model selection [8], [9]. Often Bayesian methods boil down to approximate integration problems which are solved using Markov Chain Monte-Carlo algorithms [10]. In practice, the target distribution \(\pi\) defined on \(\mathbb{R}^d\) is often absolutely continuous w.r.t. the Lebesgue measure and we have access to its density up to a normalization constant. Throughout the paper we will use the same notation for probability distributions and their corresponding densities. In Bayesian statistics the target distribution is generally given by \[\label{eq:target} \pi(x)=\frac{1}{Z} e^{-V(x)},\tag{1}\] where \(Z\) is the normalization constant and \(V:{\mathbb{R}}^d\to \mathbb{R}\) is called potential function. Here, we usually do not know the value of \(Z\). The approximate sampling methods are proceed as follows. For every \(\varepsilon > 0\), the goal is to construct a distribution \(\mu\) that approximates the target \(\pi\) in some probability distance dist: \[{\sf dist}(\mu,\pi) < \varepsilon.\] In this paper, we will study the case when \({\sf dist}\) is the Kullback-Leibler divergence: \[\label{eq:KL} \mathop{\mathrm{KL}}\nolimits(\mu\mid\pi):=\int \log\left(\frac{\mu(x)}{\pi(x)}\right)\mu({\mathrm d}x) = \mathbb{E}_{X\sim\mu}\big[V(X)\big] - H(\mu).\tag{2}\] Here, \(H(\mu) := - \int \log(\mu(x)) \mu({\mathrm d}x)\) is the negative entropy. The sampling problem can also be seen as a minimization problem in the space of probability measures (see e.g. [11]–[13]). Indeed, consider the functional \({\mathcal{F}}(\cdot) := \mathop{\mathrm{KL}}\nolimits(\cdot \mid \pi)\). This functional is non-negative and it admits its minimum value \(0\) only when for the target distribution: \(\pi = \arg\min_{\mu} {\mathcal{F}}(\mu)\). The connection of sampling and optimization has been repeatedly leveraged in the previous literature of Langevin sampling. Various algorithms such as Langevin Monte-Carlo (e.g. [13]–[15]), Underdamped Langevin Algorithm (e.g. [16], [17]) have often been influenced by the known optimization methods. Conversely, another line of research has studied the application of sampling algorithms to solve optimization problems (see e.g. [18]). In this paper, we will study another algorithm called Stein Variational Gradient Descent (SVGD), which is designed as a gradient descent algorithm for the \(\mathop{\mathrm{KL}}\nolimits\) divergence in the space of probability measures.

0.0.1.1 SVGD

The LMC algorithm, treats the functional \(\mathop{\mathrm{KL}}\nolimits(\cdot|\pi)\) as a composite functional described in 2 . For details, we refer the reader to 0.0.6 (see also [12]). Unlike the LMC, the Stein Variational Gradient Descent (SVGD) algorithm applies gradient descent directly to \(\mathop{\mathrm{KL}}\nolimits(\cdot|\pi)\) (see 0.0.2.2 for the complete definition). SVGD is an important alternative to the Langevin algorithm and already has been used extensively in different settings of machine learning, such as variational auto-encoders [19], reinforcement learning [20], sequential decision making [21], [22], generative adversarial networks [23] and federated learning [24].

The seminal work of [11] introduced SVGD as a sampling method. Since then several variants of SVGD have been proposed. Here is a non-exclusive list: random batch method SVGD [25], matrix kernel SVGD [26], Newton version SVGD [27]. However, the theoretical understanding of SVGD is still limited to the continuous time approximation of SVGD or the so called infinite particle regime. The work by [28] was the first that proved a convergence result of the SVGD in the population limit. Later, [29] studied the geometry related to the SVGD and proposed a scheme to choose the kernel. The clean analysis of the finite particle regime remains an open question (see also [30], [31]). In [32], a descent lemma was established for the SVGD in population limit in terms of Kullback-Leibler divergence. The drawback of this result was that the analysis relied on the path information of the SVGD which is unknown beforehand. [33] improved the work of  [32] and provided a clean analysis for the convergence. They assumed \(\pi\) satisfies \(\text{T}_1\) (see 0.0.2) inequality which essentially replaced the initial trajectory condition. This implied a complexity bound for the SVGD in terms of the desired accuracy \(\epsilon\) and dimension \(d\). To the best of our knowledge, this is the state of the art result on SVGD.

0.0.1.2 Contributions

The main contribution of the paper relies on its weaker set of assumptions, that allow to treat a larger class of probability distributions which includes densities with polynomials. We enlarge the class of probability distributions two-fold.

0.0.1.2.0.1 Smoothness:

The gradient smoothness assumption is very common in the sampling literature (see e.g., [13], [32]–[38]). It is formulated as the Lipschitz continuity of the gradient. That is the Hessian \(\nabla^2 V(\cdot)\) is well-defined on \(\mathbb{R}^d\) and \[\label{eq:L-smooth} \left\| {\nabla^2 V}(x) \right\|_{o p} \leq L, \quad \forall x\in\mathbb{R}^d.\tag{3}\]

Though people have made great progress towards the understanding of these algorithms, the \(L\)-smoothness condition is quite restrictive. In fact from \(L\)-smoothness condition we can easily deduce that \(V\) has at most quadratic growth rate. In particular, the large class of polynomial potentials whose order is higher than 2 does not satisfy this condition. For the LMC algorithm, several papers have proposed different methods to relax or remove this assumption [39]–[44].

However, to the best of our knowledge there is no result for the SVGD under a relaxed smoothness assumption. In this paper we introduce the concept of \((L_0,L_1)\) smoothness (see [A-l0l1] for the definition). The latter was initially proposed by [1], for the gradient clipping algorithm. The new condition with parameters \((L,0)\) is equivalent to the \(L\)-smoothness and therefore, it is indeed a weaker assumption. Important example of such functions are the higher order polynomials (see 0.0.4 for details).

0.0.1.2.0.2 Functional inequalities:

The analysis of Langevin algorithms often relies on the strong convexity of the potential function \(V\). A line of work has proposed different methods to relax or bypass this assumption. One possible solution is to slightly modify the original algorithm (e.g. [45], [46]). Another approach relies on functional inequalities such as Poincaré or logarithmic Sobolev inequalities (e.g. [37], [47], [48]). However, the verification of these inequalities is not straightforward.

Unlike the Langevin algorithm, the SVGD algorithm does not require convexity in any form. However, functional inequalities also serve for the analysis of the SVGD but with a different purpose. [33] have proved a complexity result in dimension and precision error for the SVGD algorithm performed on the targets that satisfy Talagrand’s \({\mathrm T}_1\) inequality (see 8 for details). This assumption replaces the bound on the trajectory that was initially proposed by [32]. Despite the major improvement, they still assumed the \(L\)-smoothness of the potential, which as mentioned previously does not cover the polynomials.

In this regard, we propose to replace the Talagrand’s inequality with the generalized \({\mathrm T}_p\) (see [Tp-assumption] for the definition). Due its general form, this new condition includes inequalities that are easy to verify (see 0.0.10). See also the [Table:contribution] for a visual representation of the literature review.

0.0.1.3 Paper structure

The paper is organized as follows. In 0.0.2, we present the mathematical setting of our problem. We introduce the basic notions, describe the main algorithm and list the necessary assumptions. The first part of 0.0.3 is devoted to the main result, which is a descent lemma for \(\mathop{\mathrm{KL}}\nolimits(\cdot,\pi)\) functional on the Wasserstein space. The second part contains the complexity results and discussions. 0.0.4 presents two examples. 0.0.5 concludes the results of the paper and discusses possible future work. Finally in the Appendix, the reader may find the postponed proofs, intermediary technical lemmas and supplementary discussions.

0.0.2 Mathematical setting of the problem↩︎

0.0.2.1 Notations

We will denote the random variables with uppercase (resp. lowercase) Latin letters the random (resp. deterministic) vectors, unless specified. The space of \(d\)-dimensional real vectors is denoted by \(\mathbb{R}^d\), while the set of non-negative real numbers is denoted by \(\mathbb{R}^+\). Our target distribution \(\pi\) is defined on \(\mathbb{R}^d\), which is enhanced with the Euclidean norm. The notation \(\|\cdot\|\) will correspond to the \({\ell}_2\) norm on \(\mathbb{R}^d\) unless specified. Also, we will assume that \(\pi\) has the \(p\)-th moment, that is \(\pi \in \mathcal{P}_{p} (\mathbb{R}^d)\). The Jacobian of a vector valued function \(h(\cdot) = (h_1(\cdot),\ldots,h_d(\cdot))^\top\) is a \(d\times d\) matrix defined as \[J h(x) := \big({ \partial_{x_i} h_j}\big)_{i=1,j=1}^{d,d}.\] The divergence of the vector valued function \(h\) is the trace of its Jacobian: \[\operatorname{div} h(x) := \sum_{i=1}^{d} \partial_{x_i} h_i (x).\] The Hessian of a real valued function \(U : \mathbb{R}^d \rightarrow \mathbb{R}\) is defined as the following square matrix: \[\nabla^2 U (x) := \left(\frac{\partial^2 h}{\partial_{x_i} \partial_{x_j}} (x) \right)_{i=1,j=1}^{d,d}.\]

For any Hilbert space \(\mathcal{H}\), we denote by \(\langle\cdot, \cdot\rangle_{\mathcal{H}}\) the inner product of \(\mathcal{H}\) and by \(\|\cdot\|_{\mathcal{H}}\) its norm. Moreover, \(\|\cdot\|_{\mathrm{op}}\) denotes the operator norm on the set of matrices. Let \(\mu,\nu \in \mathcal{P}_{p} (\mathbb{R}^d)\). The set of couplings \(\Gamma(\mu,\nu)\) is the set of all joint distributions defined on \(\mathbb{R}^d \times \mathbb{R}^d\) having \(\mu\) and \(\nu\) as its marginals. The Wasserstein-\(p\) distance between two probability measures is defined as \[W_p(\mu,\nu) := \inf _{\eta \in {\Gamma}(\mu, \nu)} \left[\int\|x-y\|^{p} \eta({\mathrm d}x, {\mathrm d}y)\right]^{1/p}.\] The Kullback-Leibler divergence is defined as \[\mathop{\mathrm{KL}}\nolimits(\mu\mid\nu) = \begin{cases} \int_{\mathbb{R}^d} \log\left(\frac{\mu(x)}{\nu(x)}\right)\mu({\mathrm d}x), \text{ if } \mu \ll \nu;\\ +\infty, \text{ otherwise}. \end{cases}\] We will use the spectral and the Hilbert-Schmidt norms for matrices. For \(M \in \mathbb{R}^{d\times d}\) they are respectively defined as \[\begin{align} \|M\|_{op} &:= \sqrt{\lambda_{\max} (M^{\top}M)}, \qquad \text{and} \qquad \|M\|_{HS} &:= \sqrt{\sum_{i=1}^{d} \sum_{j=1}^{d} M_{i,j}^2}. \end{align}\] Here \(\lambda_{\max}\) corresponds to the largest eigenvalue.

0.0.2.2 The definition of the SVGD

Let us present briefly Reproducing Kernel Hilbert Spaces (RKHS) and some of its essential properties. We refer the reader to [50][Chapter 4] for a detailed introduction. Let the map \(k : \mathbb{R}^d \times \mathbb{R}^d \rightarrow \mathbb{R}\) be a reproducing kernel and let \({\mathcal{H}}_0\) be its corresponding RKHS. This means that \({\mathcal{H}}_0\) consists of real-valued maps from \(\mathbb{R}^d\) to \(\mathbb{R}\), including the feature maps \(\Phi(x) := k(x,\cdot) \in {\mathcal{H}}_0\), and the reproducing property is satisfied: \[f(x) = \langle f, k(\cdot,x) \rangle_{{\mathcal{H}}_0}.\] See the discussion on the reproducing property in 0.0.10.1. Let \({\mathcal{H}}\) be the space of the \(d\)-dimensional maps \(\{(f_1,\ldots,f_d)^{\top} \mid f_i \in {\mathcal{H}}_0, i = 1,\ldots,d\}\). For two vector functions \(f = (f_1,\ldots,f_d)^{\top}\) and \(g = (g_1,\ldots,g_d)^{\top}\) from \({\mathcal{H}}\), we define the scalar product as \[\langle f,g \rangle_{{\mathcal{H}}} := \sum_{i=1}^{d} \langle f_i,g_i \rangle_{{\mathcal{H}}_0}.\]

Suppose that we have a kernel \(k : \mathbb{R}^d\times\mathbb{R}^d \rightarrow \mathbb{R}^+\) and \({\mathcal{H}}_0\) is its corresponding RKHS. As described above, we construct the Hilbert space \({\mathcal{H}}= {\mathcal{H}}_0^d\). Our goal is to construct an iterative algorithm, that has its iterates in the space of probability measures \({\mathcal{P}}_p(\mathbb{R}^d)\). Each iterate is defined as a pushforward measure from the previous one in a way that it minimizes the \(\mathop{\mathrm{KL}}\nolimits\) distance the most. To do so, for every \(\psi \in {\mathcal{H}}\) and \(\gamma > 0\) let us define the operator \[T_{\gamma}(x) := x - \gamma \psi(x).\] The operator \(\psi\) will serve us as the direction or the perturbation, while as \(\gamma\) is the step-size. The goal is to choose the direction in which the \(\mathop{\mathrm{KL}}\nolimits\) error descends the most. Thus, for every \(\mu\), the optimal choice of \(\psi\) is the solution of the following problem: \[\label{eq:opt-dir-problem} g_\mu := \arg\max_{\psi \in \mathcal{H}}\left\{- \frac{\mathrm{d}}{\mathrm{d} \gamma} \mathrm{KL}\left(T_{\gamma}{\#} \mu \| \pi \right)\Big|_{\gamma=0} \quad \text{ s.t. } \quad\|\psi\|_{\mathcal{H}} \leq 1\right\} .\tag{4}\] The functional objective depends linearly on the perturbation function \(\psi\). Indeed, [11] have proved that \[\label{eq:linear-obj} - \frac{\mathrm{d}}{\mathrm{d} \gamma} \mathrm{KL}\left(T_{\gamma}{\#} \mu \| \pi \right)\Big|_{\gamma=0} = \int_{\mathbb{R}^d} \big(- V(x)\psi(x) + \operatorname{div} \psi(x)\big) \mu({\mathrm d}x).\tag{5}\] If \(\psi\) satisfies mild conditions, then the right-hand side of 5 is equal to zero if and only if the measures \(\mu\) and \(\pi\) coincide [51]. This property motivates to define a discrepancy measure called Stein discrepancy as follows: \[I_{\rm Stein}(\mu \mid \pi) := \max_{\psi \in \mathcal{H}} \left\{ \int_{\mathbb{R}^d} \big(- V(x)\psi(x) + \operatorname{div} \psi(x)\big) \mu({\mathrm d}x) \quad \text{ s.t. } \quad\|\psi\|_{\mathcal{H}} \leq 1 \right\}^2.\] [52] have proved that the solution of the optimization problem 4 is given explicitly by \[g_\mu (\cdot) = - \int_{\mathbb{R}^d} \big[\nabla\log \pi(x) k(x,\cdot) + \nabla_x k(x,\cdot)\big]\mu ({\mathrm d}x)\] and \(\|g_{\mu}\|_{{\mathcal{H}}} = \sqrt{I_{\rm Stein} (\mu \mid \pi)}\). For the proof of this inequality, we refer the reader to the remark in 0.0.7.1. In our case, \(I_{\rm Stein}\) is often referred to as the squared Kernelized Stein Discrepancy (KSD) or the Stein-Fisher information. Under certain conditions (see 0.0.10.4) on the target distribution \(\pi\) and the kernel \(k\), the convergence in KSD implies weak convergence: \[I_{\rm Stein} (\mu_n \mid \pi) \rightarrow 0 \implies \mu_n \Rightarrow \pi.\]

This latter equality is the alternative definition of the optimal direction \(g_\mu\) given by [32]. The proof of the remark can be found in 0.0.10.5. Thus, we have determined the descent direction. In order to formulate the algorithm let us initialize our algorithm at \(\mu_0 \in {\mathcal{P}}_{p}(\mathbb{R}^d)\). The \((i+1)\)-th iterate \(\mu_{i+1}\) is obtained by transferring \(\mu_i\) with the map \(I - \gamma g_{\mu_i}\): \[\label{eq:SVGD} \mu_{i+1} = (I - \gamma g_{\mu_i})\# \mu_i.\tag{6}\] In case when \(k \in L_p(\mu)\), then \({\mathcal{H}}_0 \subset L_p(\mu)\) (see [50][Theorem 4.26]). This guarantees that every iteration of the SVGD has a finite \(p\)-th moment, that is for every \(i\in \mathbb{N}\), \(\mu_i \in {\mathcal{P}}_p(\mathbb{R}^d)\). The rest of the section is devoted to the assumptions that we will use later in the analysis.

0.0.2.3 Assumptions

In order to perform SVGD we need the target \(\pi\), in particular the gradient of its potential potential \(V\) and the kernel \(k(\cdot,\cdot)\). Below we present the four main assumptions on \(\pi\),\(V\) and \(k\), that we use in the analysis. The first assumptions is a relaxed form of smoothness condition for \(\nabla V\).

We study the SVGD algorithm under the so called \((L_0,L_1)\) smoothness, which has been borrowed from optimization literature. In [1], the authors studied the convergence of the clipped gradient descent under this condition. The assumption goes as follows. The Hessian \({\nabla^2 V}\) of \(V=-\log \pi\) is well-defined and \(\exists L_0, L_1\geq 0\) s.t. \[\label{eq:l0l1} \left\| {\nabla^2 V}(x) \right\|_{op} \leq L_0+L_1 \| \nabla V(x) \|,\tag{7}\] for any \(x\in\mathbb{R}^n\). First let us notice that when \(L_1 = 0\) the inequality 7 becomes the usual smoothness condition 3 . On the other hand, let us consider the generalized normal distribution (also known as exponential power distribution) which is given by its density \(\pi(x) \propto \exp(-\|x-a\|^{\beta})\), where \(\beta \geq 2\). In this case, Hessian \(\nabla^2 V\) is not bounded by a constant in terms of the operator norm. However, it is easy to verify that \(V\) satisfies [A-l0l1], which yields that this assumption is indeed a relaxation of the \(L\)-smoothness.

The next two assumptions are intended to replace the trajectory condition [29], [32] for polynomial potentials. However, this condition limits the applicability of the algorithm as these results cannot guarantee the convergence of the algorithm before actually implementing it. This condition was replaced by Talagrand’s \({\mathrm T}_1\) inequality in the recent work by [33], where complexity bounds are proved for \(L\)-smooth potentials. Talagrand’s \({\mathrm T}_p\) inequality goes as follows: \[\label{eq:Tp} W_p(\mu,\pi)\leq \sqrt{\frac{2\mathop{\mathrm{KL}}\nolimits(\mu\mid\pi)}{\lambda_{{\mathrm T}}}}, \quad \forall \mu\in{\mathcal{P}}_p({\mathbb{R}}^d).\tag{8}\] Since the Wasserstein-\(p\) distance is increasing with respect to the order \(p\), then \({\mathrm T}_p\) implies \({\mathrm T}_1\). In general, it is hard to verify if the distribution satisfies the \({\mathrm T}_p\) inequality. However, one may check that \({\mathrm T}_2\) is true for \(m\)-strongly log-concave distributions with \(\lambda_{{\mathrm T}} = m\). According to Bakry-Émery theory, sharp estimates on this constant are available due to perturbation arguments, such as the well-known Holley-Stroock method (see [53]). We propose a more general form of the \({\mathrm T}_p\) inequalities which will help us to analyze the convergence of the SVGD under the \((L_0,L_1)\) smoothness condition.

The target distribution \(\pi\) satisfies the generalized \({\mathrm T}_p\) inequality for some increasing function \({\mathrm J}: \mathbb{R}^+ \rightarrow \mathbb{R}^+\). If we take \(J(r) \equiv \sqrt{2r/\lambda_{{\mathrm T}}}\),for every \(r \in \mathbb{R}^+\), then we retrieve the classical \({\mathrm T}_p\) inequality. Thus, we indeed generalize \({\mathrm T}_p\) inequality with this assumption. An important example of [Tp-assumption] is a consequence of [54][Corollary 2.3]. They prove that if \[\int_{\mathbb{R}^d} \exp(\|x - x_0\|^p) \pi({\mathrm d}x) < +\infty, \text{ for some }x_0 \in \mathbb{R}^d,\] then [Tp-assumption] is satisfied for \({\mathrm J}(r) = \lambda_{BV}(r^{\frac{1}{p}} + (r/2)^{\frac{1}{2p}})\), where \(\lambda_{BV}\) is a constant that may depend on the dimension (see 0.0.10.2 for details). Therefore, in this case the assumption is essentially a condition on the tails of the target distribution \(\pi\). In particular, for the generalized normal distribution the verification of this bound is straightforward.

The third assumption is chosen to essentially restrict ourselves for potentials with (at most) polynomial growth. Mathematically, it is expressed as follows:

For some \(p>0\), there exists a polynomial with positive coefficients such that \(\mathrm{ord} ({{\mathrm Q}}) = p\) and the following inequality is true: \[\big\|\nabla V(x)\big\| \leq {{\mathrm Q}}( \|x \|).\] Using Taylor formula, one may check that \(L\)-smooth potentials satisfy this condition with \({\mathrm Q}(r) = L r + \|\nabla V(0)\|\). We want highlight that the constant \(p\) is the same for the assumptions [Tp-assumption] and [poly-assumption]. This will allow us to treat the polynomials of order \(p\). For detailed examples of distributions that satisfy our set of assumptions please refer to 0.0.4.

We conclude our list of assumptions with a bound on the kernel \(k(\cdot,\cdot)\). There exists \(B>0\) such that \(\|k(x, .)\|_{\mathcal{H}_{0}} \leq B\) and \[\left\|\nabla_{x} k(x, .)\right\|_{\mathcal{H}}=\left(\sum_{i=1}^{d}\left\|\partial_{x_{i}} k\left(x, .\right)\right\|_{\mathcal{H}_{0}}^{2}\right)^{\frac{1}{2}} \leq B,\] for all \(x \in \mathbb{R}^d\). Due to the reproducing property, these conditions are equivalent to \(k(x,x)\leq B\) and \(\partial^2_{x_i,x_i} k(x,x) \leq B\). Here, the first and the second partial derivatives are operated, respectively, on the first and the second variables of \(k(\cdot,\cdot)\). Based on this criterion, one may check that the multiquadratic kernel \(k(x,y) := (c^2 + \|x-y\|^2)^{\beta}\) for some \(c>0>\beta>-1\) satisfies the [A-ker-bound]. See appendices 0.0.10.1 for details on the reproductive property.

0.0.3 Main result↩︎

In this section we present our main results. We start with a proposition that bounds the difference between the value of the objective at two consecutive iterations. The proof of the proposition can be found in 0.0.7.1

1 may be seen as a descent lemma on the \(\mathop{\mathrm{KL}}\nolimits\) divergence. Let us develop the condition ?? . The following lemma provides us with an upper bound on \(\|g_{\mu_n}\|_{{\mathcal{H}}}\).

Thus ?? may be replaced by \[\label{eq:gamma-cond-new} \gamma \leq (\alpha-1)\min\{ 1 , {1}/{L_1} \}\left[\alpha B^2\left(\mathbb{E}_{X \sim \mu_n}\big[\|\nabla V(X)\|\big]+1\right)\right]^{-1}.\tag{9}\] On the other hand, we would like the right-hand side of ?? to be negative and it yields the following condition on the step-size: \[\label{eq:gamma-cond2} \begin{align} \gamma &< 2B^{-2}\big({\alpha^2 + (e-1) A_n}\big)^{-1} = { 2B^{-2}\big(\alpha^2 + (e-1) \big(L_{0}+L_{1}\mathbb{E}_{X \sim \mu_n}[\|\nabla V(X)\|]\big)\big)}^{-1}. \end{align}\tag{10}\] Both in 9 and 10 we have dependence on \(\mathbb{E}_{X \sim \mu_n}[\|\nabla V(X)\|]\). The lemma below establishes an upper bound on this quantity.

The important implication of 2 is that both upper bounds on the step-size are inversely proportional to \(\mathop{\mathrm{KL}}\nolimits(\mu_n\mid \pi)\). Let us now state the main theorem.

The proof is postponed to 0.0.7.2.

This yields in particular that the series with the general term \(I_{\rm Stein}(\mu_n \mid \pi)\) is convergent, that is \(I_{\rm Stein}(\mu_n \mid \pi) \rightarrow 0\). When \(\pi\) is distant dissipative and the \(k\) is the inverse multiquadratic kernel, then convergence in \(I_{\rm Stein}\) yields weak convergence of the sequence \(\mu_n\) to the target \(\pi\) (see 0.0.10.4). Proofs of the corollary can be found in 0.0.7.3. 1 yields that if \(n > \frac{\mathop{\mathrm{KL}}\nolimits(\mu_0\mid\pi)}{\gamma\varepsilon}\), then \(\frac{1}{n} \sum_{i=0}^{n-1} I_{\rm Stein}\left(\mu_{i} \mid \pi\right)<\varepsilon\). Based on this inequality, the next theorem estimates the complexity of the algorithm for some particular choices of \({\mathrm J}\) and \(\mu_0\).

The proof is postponed to 0.0.8.1. We observe that in both cases we achieve polynomial convergence of the algorithm in terms of the dimension and the precision. As mentioned in 0.0.2, in the first setting, [Tp-assumption] essentially becomes a tail condition, which can be easily verified for a wide class of log-polynomial densities. The Talagrand’s \({\mathrm T}_p\) is harder to verify in the general case. However, unlike \(\lambda_{{\mathrm T}}\), the constant \(\lambda_{BV}\) may be dimension dependent (see 0.0.10.2). Therefore, when \(p \leq 2\) Talagrand’s inequality yields better convergence. In particular, when \(p = 1\), we recover the complexity \({\mathcal{O}}(d^{3/2})\) which coincides with the result from [33]. On the other hand, when \(p\) is large, the algorithm is likely to perform significantly faster under the generalized \({\mathrm T}_p\).

0.0.4 Examples↩︎

In this section, we discuss a few cases that our method can be applied on. As mentioned in 0.0.2, the family of generalized Gaussian distributions satisfy our set of assumptions: \[\pi(x)\propto \exp\left(-\frac{\|x - a\|^p}{2\sigma^p}\right).\] This family of distributions has received considerable attention from the engineering community, due to the flexible parametric form of its probability density function in modeling many physical phenomena. For instance values of \(p=2.2\) and \(p=1.6\) have been found to model the ship transit noise and the sea surface agitation noise respectively [55]. Due to its simple explicit form of the density function, many important probabilistic quantities, such as moment, entropy etc. are easy to compute for the generalized Gaussians. It is known that that this class minimizes the entropy under under a \(p\)-th absolute moment constraint [56]. An overview of the analytic properties of this family of distributions can be found in [57]. As seen previously, when \(p\geq 2\), this family of distributions meets the assumptions of this work, so we can apply SVGD on them.

The second example has its roots in the Bayesian LASSO problem . The LASSO estimator (proposed by[58]) can be seen as the \(\ell_1\)-penalized least square estimator: This can be interpreted also as a Bayesian posterior mode (also known as Bayesian LASSO), with Laplacian priors (see [59]–[61]).

From the sampling perspective, the Bayesian LASSO boils down to \[\pi(x)\propto \exp(-f(x)-\tau\|x\|_p^p).\] Here \(f\) is the likelihood function of the prior. In [59], the authors consider the simpler case where \(f(x)=\|{A x}-{b}\|^{2}\), with \({A}\) is the measurement matrix, while \({b}\) is the vector of labels. In the case when \(p > 1\) the classical sampling methods fail to provide guarantees as the potential does not satisfy the smoothness condition. When \(p\in(1,2)\), [49] applied the Gaussian smoothing technique on the potential function \(V(x)=f(x)+\tau\|x\|_p^p\). This allowed to perform the usual Langevin Monte-Carlo algorithm on the modified potential. However, their method does not deal with the case \(p\geq 2\). Moreover, they assume the strong convexity of the potential. Our analysis for the SVGD can deal with the case when \(p\geq 2\). In particular, we do not need to put any restriction on the matrix \({A}\).

0.0.5 Conclusion and Future work↩︎

In this work, we studied the convergence of the SVGD under a relaxed smoothness condition. The latter was first proposed in the context of gradient clipping for the optimization problem. This relaxed smoothness condition allows to treat distributions with polynomial potentials. Main result consists of a descent lemma for the \(\mathop{\mathrm{KL}}\nolimits\) error of the SVGD algorithm. The result implies polynomial convergence for the average Stein-Fisher information error for two types of functional inequalities. However, at the moment this setting is of purely theoretical interest as the algorithm operates on the space of probability distributions. Despite numerous applications, the discretized algorithm does not have practical convergence guarantees in the general case. This problem remains open for future work.

Appendix↩︎

0.0.6 Langevin Monte-Carlo from the optimization perspective↩︎

Langevin algorithms have extensively been studied by statistics and machine learning communities during the last decades. The algorithms rely on an SDE called (vanilla) Langevin dynamics: \[\label{eq:LangevinDynamic} dX_t=-\nabla V(X_t)dt+\sqrt{2}dB_t,\tag{11}\] where \(B_t\) is a \(d\)-dimensional Brownian motion on \(\mathbb{R}^d\). The important property of this SDE is that under technical conditions the target \(\pi(\cdot) \propto \exp(-V(\cdot))\) is its unique invariant distribution. Moreover, it is ergodic and \(\mathop{\mathrm{KL}}\nolimits(\mu_t,\pi)\) converges linearly to zero, where \(\mu_t\) is the distribution of \(X_t\) [62]. The Euler-Maruyama discretization of this SDE over the time axis results in the Langevin Monte-Carlo (LMC) algorithm: \[\label{eq:LMC} \theta_{i+1} = \theta - \gamma \nabla V (\theta_i) + \sqrt{2h} \xi_{i+1}, \text{ for } i\in \mathbb{N}.\tag{12}\] Here \(\gamma>0\) is the discretization step-size and \(\xi_{1},\xi_{2},\ldots\) is a sequence of independent standard Gaussians independent from the initial point \(\theta_0\). This algorithm initially was studied by [3], [63]–[65] who proposed to apply Metropolis-Hastings step (MALA). The reason for this adjustment is that the constant step-size \(\gamma\) induces a bias and the target distribution \(\pi\) is no longer the invariant distribution of the discrete-time process. Later, [35] suggested to remove the Metropolis-Hastings step with a result that essentially controls the bias depending on the step-size \(\gamma\). This instigated a new line of research which studied the convergence properties of the LMC in different settings (see [34], [36], [37], [66]–[68]).

The LMC algorithm essentially performs a Forward-Flow iteration in the space of distributions (see [12], [13]). Indeed, if we define by \(\mu_i\) the distribution of the \(i\)-th iterate \(\theta_i\), then 12 is equivalent to \[\begin{align} &\mu_{i+\frac{1}{2}} = \big( I - \gamma \nabla V(\mu_i)\big){\#}\mu_i\\ &\mu_{i+1} = \mathcal{N}(0,{\rm I}_d)\textasteriskcentered\mu_{i+\frac{1}{2}}, \end{align}\] where the \(\#\) is the pushforward measure, while \(\textasteriskcentered\) is the convolution. The first equation corresponds to the gradient descent for \(\mathbb{E}_{X\sim\mu}[V(X)]\) where the argument is the measure \(\mu\). The second equation can be interpreted as the exact gradient flow of the negative entropy. The combination of this two steps results in a biased algorithm (as noticed in [3]) since the flow step is not the adjoint of the gradient step.

0.0.7 Convergence results↩︎

0.0.7.1 Proofs of the 1

For the sake of brevity we will omit the index when referring to \(g_{\mu_n}\) and will write simply \(g\). Let us define \(\phi_{t} := I - t g\) and \(\rho_t:= {\phi_{t}}_{\#} \mu_n\) for every \(t \in[0, \gamma]\). Then, applying Taylor formula to the function \[\label{eq:g9y980fd8} \varphi(t):=\mathrm{KL}\left(\rho_{t} \mid \pi\right)\tag{13}\] we have the following \[\label{eq:phi-taylor} \varphi(\gamma)=\varphi(0)+\gamma \varphi^{\prime}(0)+\int_{0}^{\gamma}(\gamma-t) \varphi^{\prime \prime}(t) {{\mathrm d}} t.\tag{14}\] By the definition of the SVGD iteration, we have that \(\varphi(0)={\mathop{\mathrm{KL}}\nolimits}\left(\mu_{n} \mid \pi\right)\) and \(\varphi(\gamma)={\mathop{\mathrm{KL}}\nolimits}\left(\mu_{n+1} \mid \pi\right)\). Let us now compute the term of 14 corresponding to the first order derivative.

The proof of the lemma can be found in 0.0.9.3. Applying the lemma to the function \(g\), we obtain the following: \[\label{eq:ineq-op-HS} \|t J g(x)\|_{op} \leq \|t J g(x)\|_{{\rm HS}} \leq t B\|g\|_{\mathcal{H}} < 1.\tag{15}\] The latter inequality is due to the condition on the step-size \(\gamma\). The bound 15 implies that \(\phi_{t}\) is a diffeomorphism. Therefore, \(\rho_{t}\) admits a density given by the change of variables formula: \[\rho_{t}(x)=\left|J \phi_{t}\left(\phi_{t}^{-1}(x)\right)\right|^{-1} \mu_{n}\left(\phi_{t}^{-1}(x)\right).\] Changing the variable of integration and applying the transfer lemma we get the following formula for \(\varphi(t)\): \[\begin{align} \varphi(t) &=\int \log \left(\frac{\rho_{t}(y)}{\pi(y)}\right) \rho_{t}({\mathrm d}y) \\ &=\int \log \left(\frac{\mu_{n}(x)\left|J \phi_{t}(x)\right|^{-1}}{\pi\left(\phi_{t}(x)\right)}\right) \mu_{n}({\mathrm d}x)\\ &=\int \left[ \log \left(\mu_{n}(x)\right) + \log \left(\left|J \phi_{t}(x)\right|^{-1}\right) - \log \left({\pi\left(\phi_{t}(x)\right)}\right) \right] \mu_{n}({\mathrm d}x). \end{align}\] Let us then compute the time derivative of \(\varphi(t)\). Taking the derivative inside and applying Jacobi’s formula for matrix determinant differentiation we obtain the following equality: \[\varphi^{\prime}(t)=-\int \operatorname{tr}\left((J \phi_{t}(x))^{-1} \frac{d J \phi_{t}(x)}{d t}\right) \mu_{n}({\mathrm d}x) - \int\left\langle\nabla \log \pi\left(\phi_{t}(x)\right), \frac{d \phi_{t}(x)}{d t}\right\rangle \mu_{n}({\mathrm d}x).\] By definition, \(d\phi_t / dt = g\). Therefore, we can use the explicit expression of \(\phi_{t}\) to write: \[\label{eq:phi-prime} \varphi^{\prime}(t)=\int \operatorname{tr}\left((J \phi_{t}(x))^{-1} J g(x)\right) \mu_{n}({\mathrm d}x) + \int \left \langle \nabla V\left(\phi_{t}(x)\right), g(x)\right\rangle \mu_{n}({\mathrm d}x) .\tag{16}\]

The Jacobian at time \(t=0\) is simply equal to the identity matrix since \(\phi_{0}={\rm I}_d\). It follows that \(\operatorname{tr}\left((J \phi_{0}(x))^{-1} J g(x)\right)=\) \(\operatorname{tr}(J g(x))=\operatorname{div}(g)(x)\) by the definition of the divergence operator. Using integration by parts: \[\begin{align} \varphi^{\prime}(0) &=-\int[-\operatorname{div}(g)(x)-\langle\nabla \log \pi(x), g(x)\rangle] \mu_{n}({\mathrm d}x)\\ &=-\int\left\langle\nabla \log \left(\frac{\mu_{n}}{\pi}\right)(x), g(x)\right\rangle \mu_{n}({\mathrm d}x). \end{align}\] Based on the alternative definition of \(g_{\mu}\) (see 1) and the reproducing property, we have \[\begin{align} \int\Big\langle\nabla\log\left(\frac{\mu_n}{\pi}\right)(x)&,g(x)\Big\rangle\mu_n({\mathrm d}x)\\ &=\iint k(x,y)\left\langle\nabla\log\left(\frac{\mu_n}{\pi}\right)(x), \nabla\log\left(\frac{\mu_n}{\pi}\right)(y)\right\rangle \mu_n({\mathrm d}x)\mu_n({\mathrm d}y)\\ &=\iint \langle k(x,\cdot),k(y,\cdot)\rangle_{{\mathcal{H}}_0}\left\langle \nabla\log\left(\frac{\mu_n}{\pi}\right)(x), \nabla\log\left(\frac{\mu_n}{\pi}\right)(y)\right\rangle \mu_n({\mathrm d}x)\mu_n({\mathrm d}y)\\ &=\left\langle \int\nabla\log\left(\frac{\mu_n}{\pi}\right)(x)k(x,\cdot)\mu_n({\mathrm d}x), \int\nabla\log\left(\frac{\mu_n}{\pi}\right)(y)k(y,\cdot)\mu_n({\mathrm d}y)\right\rangle_{{\mathcal{H}}}\\ &=\left\| g \right\|_{{\mathcal{H}}}^2. \end{align}\] Therefore, \[\label{eq:phi3940041} \varphi^{\prime}(0) = - \left\| g \right\|_{{\mathcal{H}}}^2.\tag{17}\]

Next, we calculate the term of 14 that contains the second derivative using 16 . First, \[\begin{align} \frac{d}{dt}\operatorname{tr}\left((J \phi_{t}(x))^{-1} J g(x)\right) &=\operatorname{tr}\left(\frac{d}{dt}(J \phi_{t}(x))^{-1} J g(x)\right) \\ &=-\operatorname{tr}\left((J \phi_{t}(x))^{-2} \frac{d}{dt}J\phi_t(x)J g(x)\right) \\ &=\operatorname{tr}\left(J \phi_{t}(x)^{-2} J g(x)J g(x)\right) \\ &=\operatorname{tr}\left(\left(J \phi_{t}(x)^{-1} J g(x)\right)^2\right). \end{align}\] From the definition of the function \(\phi_t\), we know that \(J \phi_{t}(x) = ({\rm I}_d + t Jg)(x)\). Thus, \((J \phi_{t})^{-1}\) and \(Jg\) commute. The latter yields \(\operatorname{tr}\left(\left(J g(x)\left(J \phi_{t}(x)\right)^{-1}\right)^{2}\right) =\left\|J g(x)\left(J \phi_{t}(x)\right)^{-1}\right\|_{\rm H S}^{2}\). On the other hand, applying the chain rule on second term of 16 yields \[\partial_t \big(\left\langle \nabla V\left(\phi_{t}(x)\right), g(x)\right\rangle\big) = \left\langle g(x), {\nabla^2 V}\left(\phi_{t}(x)\right) g(x)\right\rangle.\] Summing up, we have the following:

\[\varphi^{\prime \prime}(t) =\underbrace{\int \left\|J g(x)\left(J \phi_{t}(x)\right)^{-1}\right\|_{\rm H S}^{2} \mu_{n}({\mathrm d}x)}_ {:=\psi_1(t)} +\underbrace{\int\left\langle g(x), {\nabla^2 V}\left(\phi_{t}(x)\right) g(x)\right\rangle\mu_{n}({\mathrm d}x)}_{:=\psi_2(t)} .\] First, we bound \(\psi_{1}(t)\). Cauchy-Schwarz implies that \[\left\|J g(x)\left(J \phi_{t}(x)\right)^{-1}\right\|_{\rm H S}^{2} \leq \|J g(x)\|_{\rm H S}^{2}\left\|\left(J \phi_{t}(x)\right)^{-1}\right\|_{o p}^{2}.\] From 3, we have \(\|J g(x)\|_{\rm H S} \leq B\|g\|_{{\mathcal{H}}}\). To bound the second term, let us recall that \(\phi_t = I - tg\) and that \(t \leq \gamma\). Thus, the following bound is true: \[\begin{align} \left\|\left(J \phi_{t}(x)\right)^{-1}\right\|_{o p} & = \left\|\big( ( {\rm I}_d - t J g)(x)\big)^{-1}\right\|_{o p} \leq \sum_{i=0}^{\infty}\|t J g(x)\|_{o p}^{i} \leq \sum_{i=0}^{\infty}\|\gamma J g(x)\|_{\rm H S}^{i}. \end{align}\] Recalling ?? and combining it with 3 we obtain \[\begin{align} \left\|\left(J \phi_{t}(x)\right)^{-1}\right\|_{o p} \leq \sum_{i=0}^{\infty}(\gamma B \|g\|_{{\mathcal{H}}})^{i} \leq \sum_{i=0}^{\infty} \left(\frac{\alpha-1}{\alpha}\right)^i = \alpha. \end{align}\] Summing up, we have that \[\label{eq:psi1-bound} \psi_1(t) \leq \alpha^2 B^2\|g\|_{{\mathcal{H}}}^2.\tag{18}\] Next, we bound \(\psi_{2}(t)\). By definition, \[\begin{align} \psi_2(t) &= \mathbb{E}_{X \sim \mu_n} \left[\left\langle g(X), {\nabla^2 V}\left(\phi_{t}(X)\right) g(X)\right\rangle \right] \leq \mathbb{E}_{X \sim \mu_n} \left[ \|{\nabla^2 V}\left(\phi_t(X)\right)\|_{op} \|g(X)\|_2^2 \right]. \end{align}\] Let us bound the norm of \(g(x)\). The reproduction property of the RKHS yields the following: \[\label{eq:norm2-H} \|g(x)\|_2^{2}=\sum_{i=1}^{d}\left\langle k(x, .), g_{i}\right\rangle_{\mathcal{H}_{0}}^{2} \leq\|k(x, .)\|_{\mathcal{H}_{0}}^{2}\|g\|_{\mathcal{H}}^{2} \leq B^{2}\|g\|_{\mathcal{H}}^{2}.\tag{19}\] Therefore, \[\label{eq:psi2-initial-bound} \psi_2(t) \leq B^{2}\|g\|_{\mathcal{H}}^{2}\mathbb{E}_{X \sim \mu_n} \left[ \|{\nabla^2 V}\left(\phi_t(X)\right)\|_{op} \right].\tag{20}\] Let us bound \(\mathbb{E}_{X \sim \mu_n} \left[ \|{\nabla^2 V}\left(\phi_t(X)\right)\|_{op} \right]\). [A-l0l1] implies the following inequality:

\[\left\| \nabla^2 V (\phi_t(x))\right\|_{op} \leq L_{0}+L_{1}\|\nabla V(\phi_t(x))\|,\] for every \(x \in \mathbb{R}^d\). To bound the term \(\|\nabla V(\phi_t(x))\|\), we introduce the following lemma.

We will apply 4 to \(\phi_t(x)\) and \(\phi_0(x)\). By definition, \(\phi_t(x) - \phi_0(x) = t g(x)\) and according to inequality 19 , \[\|\phi_t(x) - \phi_0(x)\|_2\leq t B\| g\|_{{\mathcal{H}}}.\] Thus, using 4 for \(x = \phi_0(x)\), \(x^+ = \phi_t(x)\) and \(\Delta = tB\|g\|_{{\mathcal{H}}}\), we obtain the following: \[\label{eq:bound-Hess-phi} \begin{align} \left\| \nabla^2 V (\phi_t(x))\right\|_{op} &\leq L_0 + L_1 \left(\frac{L_0}{L_1}\big(\exp(t B \|g\|_{{\mathcal{H}}} L_1) - 1\big) + \left\|\nabla V\left(\phi_0(x) \right)\right\| \exp(t B \|g\|_{{\mathcal{H}}}L_1)\right) \\ &= \big(L_{0} + L_1\left\|\nabla V\left(x\right)\right\|\big)\exp (t B L_1 \|g\|_{{\mathcal{H}}} ), \end{align}\tag{21}\] Combining 20 and 21 we obtain \[\label{eq:psi2-final-bound} \psi_2(t)\leq B^2 \|g\|_{{\mathcal{H}}}^2 \left(L_{0} + L_1 \mathbb{E}_{X \sim \mu_n}\big[\left\|\nabla V (X) \right\|\big]\right) \exp (t B L_1 \|g\|_{{\mathcal{H}}} ).\tag{22}\] Summing up, the bounds on \(\psi_1\) and \(\psi_2\) yield the following inequality: \[\varphi^{\prime \prime}(t) \leq B^2 \|g\|_{{\mathcal{H}}}^2 \Big[ \alpha^2 + \Big(L_{0} + L_1 \mathbb{E}_{X \sim \mu_n}\big[\left\|\nabla V (X) \right\|\big]\Big) \exp (t B L_1\|g\|_{{\mathcal{H}}} )\Big].\] Recall that by definition \(A_n = \big(L_{0} + L_1 \mathbb{E}_{X \sim \mu_n}\big[\left\|\nabla V (X) \right\|\big]\big).\) Inserting the previous inequality along with 17 to 14 , we get the following bound: \[\begin{align} \varphi(\gamma) - \varphi(0) &\leq - \gamma \|g\|_{{\mathcal{H}}}^2 + \Big[ \frac{1}{2}\gamma^2\alpha^2 B^2 \|g\|_{{\mathcal{H}}}^2 + A_n \big( \exp (\gamma B L_1 \|g\|_{{\mathcal{H}}}) - \gamma B L_1\|g\|_{{\mathcal{H}}} - 1 \big) \Big]. \end{align}\] One can check that \(\exp(t) - t - 1 \leq (e-1)t^2/2\), when \(t \in [0,1]\). Since \(\gamma B L_1\|g\|_{{\mathcal{H}}} < 1\), we deduce \[\label{eq:final-prop} \varphi(\gamma) - \varphi(0) \leq \Big[ - \gamma + \frac{1}{2} (\gamma^2 \alpha^2 B^2 + (e-1) A_n \gamma^2 B^2)\Big] \|g\|_{{\mathcal{H}}}^2 .\tag{23}\] Finally, by the definition of Stein’s information \(\|g\|_{{\mathcal{H}}}^2 = I_{\rm Stein}(\mu_n \mid \pi)\). This concludes the proof.

0.0.7.2 Proof of 1

0.0.7.3 Proof of 1

Let us start with the first statement. Summing the descent bounds of 1 for \(i=0,1,\ldots,n-1\), we obtain the following: \[\mathrm{KL}\left(\mu_{n} \mid \pi\right)-\mathrm{KL}\left(\mu_{0} \mid \pi\right) \leq - \frac{\gamma}{2} \sum_{i=0}^{n-1}I_{\rm Stein}(\mu_i \mid \pi).\] Rearranging the terms we get \[\begin{align} \sum_{i=0}^{n-1}I_{\rm Stein}(\mu_i \mid \pi) &\leq \frac{2}{\gamma}(\mathrm{KL}\left(\mu_{0} \mid \pi\right)-\mathrm{KL}\left(\mu_{n} \mid \pi\right))\\ &\leq \frac{2}{\gamma}\mathrm{KL}\left(\mu_{0} \mid \pi\right). \end{align}\] This means that the series on the left-hand side is convergent and thus, its general term converges to zero. Thus the first point is proved. Dividing both sides of the previous inequality on \(n\), we deduce the second statement.

0.0.8 Complexity results↩︎

0.0.8.1 Proof of 2

0.0.9 Proofs of the lemmas↩︎

0.0.9.1 Proof of 1

0.0.9.2 Proof of 2

Let the polynomial \({\mathrm Q}\) have the following explicit form: \[\label{eq:q-explicit} {\mathrm Q}(r) = \sum_{i=0}^{m} a_i r^{p_i}\tag{28}\] for every \(r \in \mathbb{R}\). Here \(p = p_0 > p_1 > \ldots > p_m\) and \(a_i > 0\). Then [poly-assumption] yields \[\begin{align} \mathbb{E}_{X \sim \mu_n}\big[ \|\nabla V(X)\| \big] &\leq \mathbb{E}_{X \sim \mu_n}\big[ {\mathrm P}(\|X\|)\big] \\ &= \mathbb{E}_{X \sim \mu_n}\left[ \sum_{i=0}^{m} a_i \|X\|^{p_i}\right]\\ &= \sum_{i=0}^{m} a_i W_{p_i}\big(\mu_n,\delta_0 \big)^{p_i}, \end{align}\] where \(\delta_0\) is Dirac measure on \(\mathbb{R}^d\) at point \(0\). Using the fact that \(W_r\) is monotonically increasing w.r.t. to \(r\), we have the following \[\begin{align} \mathbb{E}_{X \sim \mu_n}\big[ \|\nabla V(X)\| \big] &\leq \sum_{i=0}^{m} a_i W_{p}\big(\mu_n,\delta_0 \big)^{p_i}\\ &\leq \sum_{i=1}^{m} a_i \big(W_p(\mu_n,\pi) + W_p(\pi,\delta_0)\big)^{p_i}\\ &\leq \sum_{i=1}^{m} a_i \left( {\mathrm J}\big(\mathop{\mathrm{KL}}\nolimits(\mu_n\mid \pi)\big) + W_p(\pi,\delta_0)\right)^{p_i}.\\ &\leq {\mathrm Q}\left( {\mathrm J}\big(\mathop{\mathrm{KL}}\nolimits(\mu_n\mid \pi)\big) + W_p(\pi,\delta_0)\right). \end{align}\] The third inequality is due to [Tp-assumption]. Recalling 28 we conclude the proof.

0.0.9.3 Proof of 3

The proof is based on the reproducing property and Cauchy-Schwarz inequality in the RKHS space. Indeed, \[\begin{align} \|J h(x)\|_{\rm H S}^{2} &= \sum_{i, j=1}^{d}\left|\frac{\partial h_{i}(x)}{\partial x_{j}}\right|^{2}\\ &=\sum_{i, j=1}^{d}\left\langle\partial_{x_{j}} k(x, .), h_{i}\right\rangle_{\mathcal{H}_{0}} \\ &\leq \sum_{i, j=1}^{d}\left\|\partial_{x_{j}} k(x, .)\right\|_{\mathcal{H}_{0}}^{2}\left\|h_{i}\right\|_{\mathcal{H}_{0}}^{2} \\ &=\|\nabla k(x, .)\|_{\mathcal{H}}^{2}\left\|h\right\|_{\mathcal{H}}^{2} \\ &\leq B^{2}\left\|h\right\|_{\mathcal{H}}^{2} . \end{align}\] This concludes the proof.

0.0.9.4 Proof of 4

The rest of the section contains the proofs of the lemmas appearing in the proof of 2. Without loss of generality, we may assume that the normalizing constant of the density of \(\pi\) is equal to 1. Thus, \(\pi(x)=e^{-V(x)}\) and the following lemma is true.

0.0.9.5 Proof of 5

0.0.9.6 Proof of 6

0.0.10 Miscellaneous↩︎

In this section, we present previously known auxiliary results that were mentioned in the paper. Some of the results are proved, the others refer to their papers of origin.

0.0.10.1 Reproducing property

Here we remind you of the reproducing property: let \(f\in{\mathcal{H}}_0\), then we have \(f(x)=\langle f(\cdot),k(x,\cdot)\rangle_{{\mathcal{H}}_0}\), if we choose \(f(\cdot)=k(y,\cdot)\in{\mathcal{H}}_0\), then \(k(y,x)=\langle k(y,\cdot),k(x,\cdot)\rangle_{{\mathcal{H}}_0}\).

Let us first use the reproducing property to calculate \(\left\| k(x,\cdot) \right\|_{{\mathcal{H}}_0}\). The definition of the norm in Hilbert spaces yields: \[\left\| k(x,\cdot) \right\|_{{\mathcal{H}}_0}^2=\langle k(x,\cdot),k(x,\cdot)\rangle_{{\mathcal{H}}_0}=k(x,x).\]

Now, let us proceed to \(\left\| \partial_{x_{i}}k(x,\cdot) \right\|_{{\mathcal{H}}_0}\). Let \(e_i\) be the \(i\)-th standard basis of \({\mathbb{R}}^d\). Then \[\begin{align} \langle \partial_{x_{i}}k(x,\cdot)&,\partial_{y_{i}}k(y,\cdot)\rangle_{{\mathcal{H}}_0}=\lim_{\epsilon_1\to 0}\lim_{\epsilon_2\to 0} \left\langle \frac{k(x+\epsilon_1e_i,\cdot)-k(x,\cdot)}{\epsilon_1}, \frac{k(y+\epsilon_2e_i,\cdot)-k(y,\cdot)}{\epsilon_2}\right\rangle_{{\mathcal{H}}_0}\\ &=\lim_{\epsilon_1\to 0}\lim_{\epsilon_2\to 0}\frac{k(x+\epsilon_1e_i,y+\epsilon_2e_i)-k(x+\epsilon_1e_i,y)-k(x,y+\epsilon_2e_i)+k(x,y)}{\epsilon_1\epsilon_2}\\ &=\lim_{\epsilon_1\to 0}\lim_{\epsilon_2\to 0}\frac{\left(k(x+\epsilon_1e_i,y+\epsilon_2e_i)-k(x+\epsilon_1e_i,y)\right)-\left(k(x,y+\epsilon_2e_i)-k(x,y)\right)}{\epsilon_1\epsilon_2}\\ &=\lim_{\epsilon_1\to 0}\frac{\partial_{y_i}k(x+\epsilon_1e_i,y)-\partial_{y_{i}}k(x,y)}{\epsilon_1}\\ &=\partial_{x_{i}}\partial_{y_{i}}k(x,y). \end{align}\] Setting \(x=y\), we obtain \(\left\| \partial_{x_{i}}k(x,\cdot) \right\|_{{\mathcal{H}}_0}^2=\partial_{x_{i}}\partial_{x_{i}}k(x,x)=:\partial_{x_{i},x_i}^2k(x,x)\), where the first \(\partial_{x_{i}}\) is operated on the first variable of \(k(\cdot,\cdot)\), the second \(\partial_{x_{i}}\) is operated on the second variable of \(k(\cdot,\cdot)\).

0.0.10.2 The generalized \({\mathrm T}_p\) inequality

The following lemma is from [54]. The proof is omitted.

We want to underline the fact that the constant here may depend on the dimension. Below, we explicit the dimension-dependence of \(\lambda_{BV}\) for the general class of distributions \(\pi(x)\propto \exp(-\left\| x \right\|^p)\) with \(p\geq 1\). It is straightforward the condition of the lemma is satisfied for \(x_0 = 0\). Now, let us compute the objective function in ?? . We start with the integral term: \[\begin{align} \int_{x\in{\mathbb{R}}^d}\exp(s\left\| x \right\|^p)\pi(x)dx &=\frac{1}{Z}\int_{x\in{\mathbb{R}}^d}\exp(-(1-s)\left\| x \right\|^p)dx\\ &=\frac{1}{Z}(1-s)^{-\frac{d}{p}}\int_{x\in{\mathbb{R}}^d}\exp(-\left\| x \right\|^p)dx\\ &=(1-s)^{-\frac{d}{p}}. \end{align}\] Inserting the previous formula into ?? , we obtain \[\label{eq:cpiporder} \begin{align} \lambda_{BV}:&=2 \inf _{x_{0} \in{\mathbb{R}}^d, s>0}\left(\frac{1}{s}\left(\frac{3}{2}+\log \int \exp({s \left\| x-x_0 \right\|^{p}} )d \pi(x)\right)\right)^{\frac{1}{p}}\\ &=2\inf_{s\in (0,1)}\left(\frac{3}{2s}-\frac{d\log(1-s)}{ps}\right)^{\frac{1}{p}}\\ &\geq 2\inf_{s\in (0,1)}\left(\frac{3}{2s}+\frac{d}{p}\right)^{\frac{1}{p}}\\ &=2 d^{\frac{1}{p}}. \end{align}\tag{30}\] Thus, indeed the constant \(\lambda_{BV}\) can be dimension dependent.

0.0.10.3 Grönwall’s integral inequality

The following lemma is the integral form of Grönwall inequality from [70].

0.0.10.4 Stein-Fisher information and weak convergence

We provide a sufficient condition on which \(\lim_{n\to\infty}I_{Stein}(\mu_n\mid\pi)\) implies \(\mu_n\to\pi\) weakly. This condition can be found in [71]. Since Stein Fisher information \(I_{Stein}(\cdot\mid\pi)\) depends on the target distribution \(\pi\) and the kernel \(k(\cdot,\cdot)\), we need the following two properties, respectively:

1. \(\pi\) is distant dissipative, that is \(\kappa_{0} \triangleq \liminf _{r \rightarrow \infty} \kappa(r)>0\) with \[\kappa(r)=\inf \left\{2 \frac{\langle \nabla V(x)-\nabla V(y), x-y\rangle}{\|x-y\|^{2}}:\|x-y\|=r\right\}.\] If \(V\) is strongly convex outside a compact set, then \(\pi\) is distant dissipative, for instance \(V(x)=\left\| x \right\|^{2+\delta}\) with \(\delta\geq 0\).

2. \(k(\cdot,\cdot)\) is an inverse multiquadratic kernel, i.e., \(k(x, y)=\left(c^{2}+\|x-y\|^{2}\right)^{\beta}\) for some \(c>0\) and \(\beta \in(-1,0)\). It is easy to check that [A-ker-bound] is satisfied.

0.0.10.5 Proof of 1

Chain rule implies \[\begin{align} \int_{\mathbb{R}^d} \nabla &\log\left(\frac{\mu(x)}{\pi(x)}\right) k(x,\cdot) \mu({\mathrm d}x) \\ &= \int_{\mathbb{R}^d} \nabla\log({\mu(x)}) k(x,\cdot) \mu({\mathrm d}x) - \int_{\mathbb{R}^d} \nabla\log(\pi(x))k(x,\cdot)\mu({\mathrm d}x)\\ &= \int_{\mathbb{R}^d} \nabla {\mu(x)} k(x,\cdot) {\mathrm d}x - \int_{\mathbb{R}^d} \log(\pi(x)) \nabla_x k(x,\cdot)\mu({\mathrm d}x).\\ \end{align}\] Since \(\mu(x) \rightarrow 0\), when \(\|x\|\rightarrow +\infty\), the integration by parts yields \[\begin{align} \int_{\mathbb{R}^d} \nabla {\mu(x)} k(x,\cdot) {\mathrm d}x &- \int_{\mathbb{R}^d} \log(\pi(x)) \nabla_x k(x,\cdot)\mu({\mathrm d}x) \\&= - \int_{\mathbb{R}^d} \left\{{\mu(x)} \nabla k(x,\cdot) + \log(\pi(x)) \nabla_x k(x,\cdot)\right\}\mu({\mathrm d}x) \\ &=g_{\mu}(\cdot). \end{align}\] This concludes the proof.

References↩︎