Towards Understanding the Dynamics of Gaussian-Stein Variational Gradient Descent
Authors: Tianle Liu
Promit Ghosal
Krishna Balasubramanian
Natesh S. Pillai
Venue: NeurIPS 2023 poster (2023)
Keywords: Stein variational gradient descent
Gaussian variational inference
Rates of Convergence

Abstract: Stein Variational Gradient Descent (SVGD) is a nonparametric particle-based deterministic sampling algorithm. Despite its wide usage, understanding the theoretical properties of SVGD has remained a challenging problem. For sampling from a Gaussian target, the SVGD dynamics with a bilinear kernel will remain Gaussian as long as the initializer is Gaussian. Inspired by this fact, we undertake a detailed theoretical study of the Gaussian-SVGD, i.e., SVGD projected to the family of Gaussian distributions via the bilinear kernel, or equivalently Gaussian variational inference (GVI) with SVGD. We present a complete picture by considering both the mean-field PDE and discrete particle systems. When the target is strongly log-concave, the mean-field Gaussian-SVGD dynamics is proven to converge linearly to the Gaussian distribution closest to the target in KL divergence. In the finite-particle setting, there is both uniform in time convergence to the mean-field limit and linear convergence in time to the equilibrium if the target is Gaussian. In the general case, we propose a density-based and a particle-based implementation of the Gaussian-SVGD, and show that several recent algorithms for GVI, proposed from different perspectives, emerge as special cases of our unified framework. Interestingly, one of the new particle-based instance from this framework empirically outperforms existing approaches. Our results make concrete contributions towards obtaining a deeper understanding of both SVGD and GVI.

Towards Understanding the Dynamics of Gaussian-Stein Variational Gradient Descent
Tianle Liu
Promit Ghosal
Department of Statistics Harvard University Cambridge, MA 02138
Department of Mathematics
Massachusetts Institute of Technology
Waltham, MA 02453 promit promitemit... edu
tianleliu@fas. www.eliiu@faa..aaaa..e.eee harvard.
Krishnakumar Balasubramanian
Natesh Pillai
Department of Statistics University of California, Davis
Department of Statistics Harvard University
Davis, CA 95616 kbala@ucdavis edu
Cambridge, MA 02138 pillai falliefaa..o www.aardddd. edu
Abstract
Stein Variational Gradient Descent (SVGD) is nonparametric particle-based deter- ministic sampling algorithm. Despite its wide usage, understanding the theoretical properties of SVGD has remained challenging problem. For sampling from a Gaussian target, the SVGD dynamics with a bilinear kernel will remain Gaussian long the initializer is Gaussian. Inspired by this fact, undertake detailed theoretical study of the Mannssian-VGii i.e., SVGD projected the family Gaussian distributions via the bilinear kernel, equivalently Gaussian variational inference (GVI) with SVGD. We present a complete picture by considering both the mean-field PDE and discrete particle systems. When the target is strongly log-concave, the mean-field linearly to the Gaussian distribution closest to the target in KL divergence. In the finite-particle setting, there is both uniform in time convergence to the mean-field limit and linear convergence in time to the equilibrium if the target is Gaussian. In the general case, we propose a density-based and particle-based implementation of the Gaussian-SVe SVGD, and show that several recent algorithms for GVL, proposed from different perspectives, emerge as special cases of our unifying framework. In- terestingly, one the new particle-based instance from this framework empirically outperforms existing approaches. Our results make concrete contributions towards obtaining a deeper understanding of both SVGD and GVI.
Coussian-SVVCI dynamics is proven to converge
Introduction
Sampling from given target density arises frequently in Bayesian statistics, machine learning and applied mathematies. Specifically, given potential V : Bd - R. the target density is given by
p(x) -V(o), where Z= V(x) is the normalizing constant.
Traditionally-used Markov Chain Monte Carlo (MCMC) sampling algorithms are invariably not scalable large-scale datasets [8] 61]. Variational inference and particle-based methods are two related alternatives proposed in the literature, both motivated by viewing sampling as optimization over the space of densities. We refer to 33 for additional details related to this line of works.
In the literature on variational inference, recent efforts have focused on the Gaussian Variational Inference (GVD) problem. On the theoretical side, this is statistically motivated by the Bernstein-von
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Mises theorem, which posits that in the limit of large samples posterior distributions tend to be Gaussian distributed under certain regularity assumptions. We refer to [T Chapter 10] for details of the classical results, and to 34
73 for some recent non-asymptotic analysis. On the algorithmic side, efficient algorithms with both statistical and computational guarantees are developed for GVI
[14 20 35, 45 19.. From a practical point-of-view, several works [64] 78 80, 67] have shown
superior performance of GVI, especially in the presence of large datasets.
Turning to particle-based methods, [53] proposed the Stein Variational Gradient Descent (SVGD) algorithm, kernel- deterministic approach for sampling. has gained sigmificant attention in the machine learning and applied mathematics communities due to its intriguing theoretical
properties and wide applicability [23] 31 [56] [87. Researchers have also developed variations of
SVGD motivated by algorithrmc and applied challenges [89] 49 [8] (28 84 12 7T. In its original
form, SVGD could be viewed as a nonparametric variational inference method with kernel-based practical implementation.
The flexibility offered by the nonparametric aspect of SVGD also leads unintended consequences. On one hand, from a practical perspective, the question of how to pick the right kernel for imple- menting the SVGD algorithm unclear. Existing approaches are mostly ad-hoc and do not provide clear instructions on the selection of kernels. On the other hand, developing a deeper theoretical understanding of SVGD dynamics is challenging due to its nonparametric formulation. Notably [58] derived the continuous- time PDE for the evolving density that emerges as the mean-fiee -ield limit of the finite-particle SVGD systems, and shows the well- sell-posedness of the PDE solutions. In general, the following different types of convergences could be examined regarding SVGD:
Initial particles 1-11
Evolving panieles -5 Σ
(d)
Equilibrium
8 3
(0)
>
(0)
(a)
(e)
Initial density Po
Evolving density
(b)
Target
Ps
p*
(a) Unified convergence of the empirical measure for N finite particles to the continuous target as time and N jointly grow to infinity;
(b) Convergence of mean-field SVGD to the target distribution over time;
(c) Convergence of the empirical measure for finite particles to the mean-field distribution at any finite given time t E [0, oo);
(d) Convergence of finite-partiche SVGD to the equilibrium over time;
(e) Convergence of the empirical measure for finite particles to the continuous target at time xc.
From a practical point of view (a) is the ideal type of result that fully characterizes the algorithmic
hehavior of SVGD, which could be obtained by combining either [b) [c] [d] and (e) Regarding
(b) [51) showed the convergence of mean-field SVGD in kernel Stein discrepancy (KSD), 17 52 29
which is known to imply weak convergence under appropriate assumptions. [40 15 70 75 22
sharpened the results with weaker conditions or explicit rates. [32] extended the above result to
the stronger Fisher information metric and Kullback-Leibler divergence based on a regularization
technique. [58
obtained -ime-dependent mean-field convergence (c)
[c] of N particles under
various assumptions using techniques from the literature of propagation of chaos. 72 obtained even
stronger results for [c] and combined (b) to get the first unified convergence (a) terms of KSD. However, they have a rather slow rate 1/ log log N, resulting from the fact that their bounds for (c) still depends on the time (sum step sizes) double-exponentially. Moreover, there has not been any work that studies the convergence (d) and (e) for SVGD, which illustrate a new way to characterize the unified convergence [a)
In an attempt to overcome the drawbacks of the nonparametric formulation of SVGD and also taking cue from the GVI literature, in this work we study the dynamics of Gaussian-SVGGI a parametric formulation of SVGD. Our contributions in this work are three-fold:
Mean-field results: We study the dynamics of Gaussian-SVGD in the mean-field setting and establish linear convergence for both Gaussian and strongly log- log-coneave targets. As an example of the obtained results, Table shows the convergence rates of covariance for centered Gaussian

---

Table Convergence rates of SVGD with different bilinear kernels for Gaussian families.
K1-SVGD
K2-SVGD
WGF (K3-SVGD)
R-SVGD (K--VVGD)
Centered
O ?!)
3.2
O(e 21)
[22
O(e $)
F.1
Ole (1
3.4
Gaussian
General
Theorem 3.1
[F2
O
F.
O - - -
3.4
Gaussian
families for several revelant algorithms. All of them will be introduced later in Sections and 3 We also establish the well-posedness of the solutions for the mean-field PDE and discrete particle systems that govern SVGD with bilinear kernels (see Appendix C). Prior work [58] requires that the kernel be radial which rules out the important class of bilinear kernels that we consider. [2] relaxed the radial kernel assumption, however, they required boundedness assumptions which avoid in this work for the case of bilinear kernels.
Finite-particle results: We study the finite-particle SVGD systems in both continuous and discrete time for Gaussian targets. We show that for SVGD with bilinear kernel if the target and initializer are both Gaussian, the mean-field convergence can be uniform in time (See Theorem 3.7] To
the best of our knowledge, this is the first uniform in time result for SVGD dynamics and should contrasted with the double exponential dependency on for nonparametric SVGD [58] 72).
Our numerical simulations suggest that similar results should hold for certain classes of non- Gaussian target as well for Gaussian- SVGD. We also study the convergence [d] directly solving the finite-part particle systems (See Theorems 3.6 and 3.8). Moreover, in Theorem 3.10 assuming centered Gaussian targets, we obtain a linear rate for covariance convergence in the finite- satte-patticles discrete-time setting, precisely characterizing the step size choice for the practical algorithm.
Unifying algorithm frameworks: We propose two unifying algorithm frameworks for finite- particle, discrete-time implementations of the CouussianSVVGi dynamics. The first approach assumes access to samples from Gaussian densities with the mean and covariance depending the current time instance of the dynamics. The second is purely particle-based approach, in that, it assumes access only to an initial set of samples from Gaussian density. In particular, show that three previously proposed methods from (27 43 for GVI emerge as special cases of the proposed frameworks, by picking different bilinear kernels, thereby further strengthening the connections between GVI and the kernel choice in SVGD. Furthermore, we conduct experiments for eight algorithms that can be implied from our framework, and observe that the particle- based algorithms are invariably more stable than density- ones. Notably one of the new particle- based algorithms emerging from our analysis outperforms existing approaches.
2 Preliminaries
fp Denote the space of probability densities on Rd by P(R")= == {o F(R*): : pdz = 1, 0). where F(Rd)) is the set of smooth functions. As studied in [42]. P(IR)) forms Fréchet manifold called the density manifold. For any "point" p P(Ra). we denote the tangent space and cotangent space at p by (TR(( and (RR)) respectively. A Riemannian metric tensor assigns to each p E P(RRd) positive definite inner product y, T (RR)) T, ((Rd) + R and uniquely corresponds to an isomorphism G, (called the canonical isomorphism) between the tangent and cotangent bundles [66]. i.e., we have G, (RR)) - (TE))
Definition 2.1 (Wasserstein metric). The Wasserstein metric is induced by the following canonical isomorphism GWaas ((RR) + P(Rd)) such that
GG Wasss 1 V (DVΦ), $ TTR(())))
The Wasserstein gradient flow (WGF) can be seen as the natural gradient flow for KL divergence on the density manifold with respect to the Wasserstein metric. In specific, the mean-field PDE is given by the linear BokkerrPPanncl equation [33]
P - Wass 1
(p: p²) = V (PLV & KL(P, 1l p0)) =V (V (1) op:
opi
3

---

where Sps denotes the variational derivative with respect to Pt. KL(p 11 p²) is the so-called energy function, and (ax) is the potential function that satisfies p²(x) x exp( exp((V(ax)).
Interestingly, the mean field flow of SVGD can also seen as a gradient flow for the KL divergence but with respect to the Stein metric, where we perform kernelization in the cotangent space before taking the divergence [51 222.
Definition 2.2 (Stein metric). The Stein metric is induced by the following canonical isomorphism
G T(P(( + ((R))) such that
(Giici)) 14 - (p(-) K(, V)ll)
& E ETT(()))
where : Fd F R. a positive-definite kernel.
In particular, the mean field PDE of the SVGD algorithm can be written as
P: = (G
Apc KL |I p²) = V (pc(-), ((.)))) + (r(g) dy). (2)
Gaussian Famílies as Submanifolds. We consider the family of (multivariate) Gaussian densities
Pa E P(()) where 0 (H, 2) E ② _dd x Sym" (d, R). Note that Sym (d, IR) is the set of (symmetric) positive definite x matrices. In this way, can be identified as a Riemannian submanifold of the density manifold P(R4) with the induced Riemannian structure. If we further restrict the Gaussian family to have zero mean, it further induces the submanifold Og Sym+ RR). Notably the Wasserstein metric on the density manifold induces the Bures Wasseesassnn metric for the Gaussian families 66 82
In paper, however, consider the induced Stein metric on θ Oo (we call it the Gaussian- metric; for details see Appendix B
Different Bilinear Kernels and Induced Metrics. There are several different bilinear kernels that y + while (2 [3
appear in literature. [44] considers the simple bilinear kernel K₁ (x. y)
suggest the of an affine-invariant bilinear kernel K2(2,, y) (xc (1) y 1 [13 further points out that with the rescaled affine- invariant kernel (3
(x [1) -y (4) the
Gaussian- Gaussian-Stein metric magically coincides with the Bures- Wasserstein metric on è (not true on the whole density manifold). Note that here μa and Σ are the mean and covariance of the current mean- field Gaussian distribution which could change with time. Moreover, (32 proposed a regularized version of SVGD (R-SVGD) that interpolates the dynamics between WGF and SVGD. Interestingly R-SVGD with the affine-invariant kernel K₂ for Gaussian families can be reformulated as Gaussian- SVGD with a new kernel Ka y) (ar 11) ((1 v)Σ + -/)--((p (y - 1z), which interpolates between K₂ and Kg (see Theorem 3.4). For clarity we present the results for K in the main article while leave the analogues for K2-Kı in the appendix. We also point out that Kj and K2 are the
same on Og.
Gaussian-Stein Variational Gradient Descent. With bilinear kernel and Gaussian targets, will prove in the next subsection that the SVGD dynamics would remain Gaussian as long as the initializer is Gaussian. However, this is not true in more general situations especially when the target is non-Gaasiaan. Fortunately for Gaussian variational inference we can still consider the gradient flow restricted the Gaussian submanifold. In general, we denote by Gbeei the canonical isomorphism on θ induced by GStcin and define the Gaussian-Stenn variational gradient descent as
à, = -G KTT(po, || p*),
where p" might not be Gaussian density. Notably Gaussian-SVGr solves the following optimization problem
min KL(ps p²), where pa is the density of N(,, E),
via gradient descent under the Caassian-Steii metric.
Dynamics of Gaussian-SVGD for Gaussian Targets
3.1 Mean-Field Analysis from WGF to SVGD
The Wasserstein gradient flow (WGF) restricted to the general Gaussian family θ known as the Wassersesnnn gradient flow 43 19]. For consistency in this subsection always set

---

the initializer to be (uo; Eo) with density Po and target target.W(b. Q) with density p" for the general Gaussian family. Then the WGF at any time remains Gaussian, and can be fully characterized by the following dynamics of the mean la and covariance matrix Eg:
P. = -- (14 b),
Σ = 21 - E2Q - 100---Q-110 2
(3)
For SVGD with bilinear kernels we have similar results:
Theorem 3.1 (SVGD). For any t. 2 O the solution of SVGD 2 with the bilinear kernel K remains Gaussian density with mean 14, and covariance matrix En given by
j 1-1---1221 E = -3, + -
(1 u1
(4)
(4:-b0) b) En
which has unique global solution on [0, 20) given any 1o E Rd and Eo E Sym* (d. IR). And P converges weakly to p* as + xo at the following rates
|| b ()- 1 - QI1 Ole ((e-2(r-)))) Ve>0,
where y is the smallest eigenvalue of the matrix
La
1 8 Q-112 v2
1
with a lower bound 7 1 - b b + 2X
b & QQ1122 1/2
(1+b
V22
where X is the largest eigenvalue of Q.
Note that for any vector æ. || denotes its Euclidean norm and for any matrix A we use || for its spectral norm, ||A]], for the nuclear norm and ||A|| for the Frobenius norm. All matrix convergence are considered under the spectral norm in default for technical simplicity (though all matrix norms are equivalent in finite dimensions).
If restrict to the centered Gaussian family where both the initializer and target have zero mean (setting N = b = 0), the dynamics can further be simplified.
Theorem 3.2 (SVGD for centered Gaussian). Let po and p" be two centered Gaussian densities. Then for any t 2 0 the solution Pt of SVGD ② with the bilinear kernel K₁ or K₂ remains a centered Gaussian density with the covariance matrix E, given by the following Riccati rype equation:
Σ 2 EQ- 1-1
(5)
which has a unique global solution on (0, 00) given any Eo, Q Sym (d. 1f QEa, we have the closed-form solution:
= + (1 (1- c 2 20=
(6)
In particular, if we let Eo I and Q = I + Laak for some v >> O and v E Rd such that v then En can be rewrilten as
1,
Σ, = I 111---2) 1+ne-2 ve
13 shows that for WGF if E. = QTag. then we have I|, - b|= ((2) and | QII
O(e-21/X) For the centered Gaussian family SVGD converges faster X > For the general Gaussian family WGF and SVGD have rather comparable rates (e.g.. take A >>>>>> ||b]] then the lower bound here is roughly ((e-(A))) Another observation is that the WGF rates depend on Q alone but the SVGD rates here sometimes also depend on b, which breaks the affine invariance of the system. This is problem originated from the choice of kernels as addressed in [13.. where they propose to use K2 instead of K1. Such approach has both advantages and disadvantages. The convention in SVGD is that the kernel should not depend on the mean-field density because the density is usually unknown and changes with time. But for GVI the affine-invarian bilinear kernel K₂ only requires estimating the means from Gaussian distributions and is not a big issue
Regularized Stein Variational Gradient Descent. Section show that SVGD can be regarded as WGF kernelized in the cotangent space (RRd). The regularized Stein variational gradient descent (R-SVGD) (32 interpolates WGF and SVGD by pulling back part of the kernelized gradient of the cotangent vector 4 which is also seen as gradient flows under the regularized Stein metric:
5

---

Definition 3.3 (Regularized Stein metric). The regularized Stein metric is induced by the following canonical map (G, -
(p ((1- vI) TK.
where Tk.p is the kernelization operator given by T == [- (()))))()))()) dy.
The R-SVGD is defined
P -Gp))
& (L. ll (²)) - (p::p(ppp ()Tk. + +UT KK.p. V log opc
f p'
(8)
Theorem 3.4 (R-SVGD). Let po and p¹ be two Gaussian densities. Then the solution Pr. of R-SVGD 8 with the bilineur kernel K2 converges to p² as l - oo, and pr. is the density of N(, Ex) with
i = -Q= (44,-b bb)
(1-1-)))1 - -1000-00000000000000000000000000 Σ = 2((1-))z + +11 (11-
+ +
(9)
If VEQQ= -
-Eo, QE,, we have | - Q11 = (--2/(0)>)tt)) where X is the largest eigenvalue of Q-
From this theorem we see that R-SVGD can take the advantage of both regimes by choosing wisely. Another interesting connection is that on θ the induced regularized Stein metric coincides with the Stein metric with a different kernel (4 y) = (x 1s) T ((1 - ك)Σ + uI) Theorem B.7p
14) + (see
Stein AIG Flow. Accelerating methods are widely used in first-order optimization algorithms and have attracted considerable interest in particle-bass based variational inference [50]. [76 86] study the accelcrated information gradient (AIG) flows as the analogue of Nesterov's accelerated gradient method [63] on the density manifold. Given probability space P(RR)) with a metric tensor 9p(---), -),
let Go TPR( + TFP((R be the corresponding isomorphism. The Hamiltonian flow in probability space [16 follows from
11
$ H (pt, 2c) B0z
0:
Pt 42
-
-1
0
H (pt, p.) 34.
where H(pi, k.) := 1. G.G. G.G.G./... 1 d. + KL(P |I p*)
pr
is called the Hamiltonian function, which consists of a kinetic energy 1 [ dæ and a potential energy KL(p || p*). Following (86) we introduce the accelerated information gradient flow in probability space. Let a 2 0 be scalar function of time t. We add damping term a b2 to the Hamiltonian flow:
[p 0
0
[
0
1
H (pe. 2+)
Φ+
at $
+
0
H (pr. 22)
(10)
By adopting the Stein metric we obtain the Stein AIG flow (S-AIGF):
Pe - (pt (()
(11)
- - / (V)))) ()) - opt KL p*)
Again we characterize the dynamics of S-AIGF with the linear kernel for the Gaussian family: Theorem 3.5 (S-AIGF). Let po and p² be two centered Gaussian densities. Then the solution Pe of S-AIGF II with the bilinear kernel K₁ or K2 is the density of N Ex) where E, satisfies
ΣS = 2(S - E.S.)
(12)
= - - 2(5, + EyS?) +
where St E Sym(d, IR) with initial value So 0.
Note that the convergence properties of S-AIGF still remains open in contrast to the Wasserstein AIG flow (W-AIGF) as [86] shows that the W-AIGF for the centered Gaussian family is
Y 2(S E,
š, = - - + | -
and that if or well-chosen, the KL divergence in W-AIGF converges at the rate of Ole Thus, when À is large it converges faster than WGF. It is also interesting to point out that the acceleration effect of Nesterov's scheme also comes from time discretization of the ODE system (see 774 as moves roughly Ve rather than along the gradient path when the step size is

---

3.2 Finite-Particle Systems
In this subsection, we consider the case where N 30 particles evolve time t. We set a Gaussian target Q) (i.e., the potential is V(x) ( b) with a bilinear kernel, and obtain the dynamics of æ SVGD corresponds to the following deterministic interactive system in RR:
Q [ b)) and run the SVGD algorithm
I N
The continuous- time particle-based
a ( = Σ Va, K
(t) (1)
N N Σ K
VV
(t) ()))
(13)
with initial particles given by zi (0)
(i 1,... ... ,). Now denoting the sample mean and covariance
matrix time by M4 N m== ( and C. Eド= (1) (() following theorem.
2 -.tt, we have the
Theorem 3.6 (SVGD). Suppose the initial particles satisfy that Co non-singular. Then SVGD with the bilinear kernel K₁ and Gaussian potential V has unique solution given by
(3)
x Ac(ar, - PP) + Pno
(14)
where As is the unique (matrix) solution of the linear system
A. (I-Q--AC: + p. + Q Q- Ac: Ao I,
(15)
and pe and Ct are the unique solution of the ODE system
L4 = Q - (1 p4 (2)a)QQ1 (1+-b)
(16)
- Co (4, -)) Q
-0=220 -
- - -2
Q-1 102-0(( (H, C.
The ODE system (T6) is exactly the same as that in the density flow 4.. Thus, we have the the same convergence rates in Theorem 3.1 Theorem c.6 be interpreted as: At each time the particle positions are a linear transformation of the initialization. On one hand, if we initialize i.i.d. from Gaussian, there is uniform in time convergence as shown in the theorem below.
On the other hand, if we initialize i.i.d. from a non-Gaussian distribution po. At each time the mean field limit Pr should be a linear transformation of po and cannot converge to the Gaussian target p² as I ++ 0o. Note that in general SVGD with the bilinear kernel might not always converge to the target distribution for nonparametric sampling but for GVI there is such issue. This will be discussed in detail in Appendix together with general results of well-pooe polsededsssss and mean-fien -feld convergence of SVGD with the bilinear kernel, which has not yet been studied in literature.
Theorem 3.7 (Uniform in time convergence). Given the same setling as Theorem 3.6 further suppose the initial particles are drawn i.i.d. from N(Ho: Ea). Then there exists constant Ca F.F... such that for all & E [0, o for all N 2 2, with the empirical measure C/ N E 8 l0 the second
moment of W distance between eQ and Pt converges:
{
lol log N N N N)2
if = 1 if 2
W (<1)) pi)] S Ca. x
(17)
2/2
if 23 3
Similar to Theorem 3.2 we also provide the finite-particle result for centered Gaussian target.
Theorem 3.8 (SVGD for centered Gaussian). Suppose the SVGD particle system (3) with the
bilinear kernel K₁ or K₂ is targeting a centered Gaussian distribution and initialized by (a
(0) N
such that 1o= = 0 and CoQ = QC2. Then we have the following closed-form solution
(c-"II 11-=
-1/2 (0) x
(18)
Analogous Result for R-SVGD. Next we consider the particle dynamics of R-SVGD. As shown in 33], the finite-pann particle system of R-SVGD is
x-= v) K. + vv) ( C INN LVK (a)²) :))
where X, : (1) (1)
(1) ,.., x
T
L.f = (f(a()) ... 1m
for all f Rd. + Rd and
(Ki)g == K [ (111) for all
S i,j S N. Similar to Theorem 3.6 we have the following result:
7

---

Theorem 3,9 (R- R-SVGD). Suppose the R-SVGD system 13 with K₁ or K2 is targeting a centered
Gaussian distribution and initialized by (a, ) )im1 such that ph = Then we have x; (t) A, where As is the unique solution of the linear system
A -I- -))()) As, Ao=
(19)
and the sample covariance matrix C₁ given by
(00)
Again we observe that the particles at time / is time-changing linear transformation of the initializers.
Discrete-time Analysis for Finite Particles. Next we consider the algorithm in discrete time I. The SVGD updates according to the following equation:
(8) (t-1) = x; + N
N V ( (1) az)) m.
Σラ= K(x)" 20 ((o)))))
(21)
For simplicity, we only consider the case where both the target and initializers are centered, i.e., b 140 and show the convergence:
Theorem 3.10 (Discrete-time convergence). For a centered Gaussian target, suppose the particle
N system 21 with K1 or K₂ is initialized by (z: i=1 such that Po = O and CoQ = QCu. For 0 <<<< 0.5, we have pe, 0 and ||C, QII + 0 as long as all the eigenvalues of Q- lie in the interval (0,1+ + 1/c).
Furthermore, if we set the to be the smaller root of the equation f((u) = 1 - (it has distinct roots) where f((x) = (1 + e(1 =));,, then we have linear convergence, i.e.,
IC - Q1 (1 I) Q1I S e e Q11
as long as all the eigenvalues of Q- lie in the interval [u., 1/3 + 1/(3e)]
The above result illustrates that firstly the step sizes required are restricted by the largest eigenvalue of Q- In particular if Cp La then we need smaller step size if the smallest eigenvalue of Q is smaller, which corresponds to the 3--logsmmoothnes condition of the target distribution. Secondly we can potentially have faster convergence over iteration given larger step sizes. We helieve that the commutativity assumption in Theorem 10 can be relaxed and similar results can be obtained for general targets. Detailed examinations are left as future work.
Beyond Gaussian Targets
In this section we consider the Gaussian-SVGD with general target and have the following dynamics.
Theorem 4.1. Let 6² be the density of the target distribution with the potential function V(x) that satsifies Assumption C.1 Rad po be the density of N (Mpg Eo). The Gaussian-SVi SVGD with K₁ produces a Gaussian density Pt. with mean se, and covariance matrix En given by
L (1 (1,1 pe (1+ pu p.)m
(22)
Σ== =5,25, E, (E,Te + p.m.)) - (TeEa + m.1u.) Ex
VV2 where I' = LVV (x)) and m, = Ex (x).
Furthermore, suppose that 0* is the unique solution of the following optimization problem
min KL,(p II p²), where po is the Gaussian measure N E).
6=(,,))
Then we have Pt + Pg- N(u E²) as 00.
In particular, if the target is strongly log-concave, it gives rise to the following linear convergence:

---

Algorithm Density-based Gaussian-SVGI.
Algorithm 2 Particle-hased Gaussian-SVCC
for t in T do (1) Draw (w) )== from from.N (W Ex)
Draw (æe :=
(
from (N En)
for l in 0 do
ΣN
Update m, and I. using (24) or 25
14 s N
k=
2
Me+1 +- 14, +1 F(UD, Et, M4, Te)
Σ + 1
Ek=1
(t)
peli
N
M441 1 e2 G(u De, m4, l()
Update my and Te using
24 (55
Et+1 t- M4t1 M,
Update (æ; (1+1),
using 26
end for
end for
Theorem 4.2. Assume that the target o² is a-strongly log-concave and 8-log-smooth, i.e,, ad 1
(~) 3 BI. Then Pr. of Theorem 4. converges pa- at the following
l|,, --1=1| O0-
| Σ* 0 ((- Ve 0,
where 7/a the smallest eigenvalue of the matrix
La 60 Ex T 3 (2
f1" (E7)1/2 11 μ²) La
with lower bound y 3(L μ p²) +1
Typically the 3-log-smoothness condition is not required for continuous-time analyses. However, is required in the above statement, as our proof technique is based on comparing the decay of the energy function of the flow to that for WGF, following and we leave it as future work.
13 Relaxing this condition is interesting
Unifying Algorithms. For general targets, we propose two unifying algorithm frameworks where we can choose any bilinear kernel (c.g.. Kı-K4) to solve GVI with SVGD. The first framework is density-based where update 14, and E, according to the mean-field dynamics. It requires the closed-form of the ODE system
j4. F (Mrs De,
(23)
Σ, = Ex G (11 L,, m. Tτ) + G (14, E., m. m,,T,) TT) Σ,
where m, = Ex--, [V)) and I = L
(æ)]. and F and G are some closed-form func- tions. For example F (Je:: Σt, m, (,) = (1 -I-= r,E,) p4 (1 + pa p,)m, and G (J4, Σ,, m,, r,) = - E,T, - n.m, for K1 as shown in 22p. Note that m, and F, can be estimated from samples
using
m = - (V(a
N =
+ EN (2V(
(24)
or using the first-order estimator
Vi î I K=
(0) p,, TE.!
(25)
The second framework is particle-based and does not need the closed-form ODE of the mean and covariance, making it more flexible than Algorithm Here we intially draw particles and keep updating them over time using
(1+1) a
x
+ N
Vay K K(x (1) (0) (( - En K(n vV)
(26)
where VV is a (time-dependent) linear approximation of VV defined as (V = i (x m,) +m, Intuitively this is used instead of VV to ensure the Gaussianity of the particle system.
We now remark on the algorithms proposed in the literature that emerge as instances of the two proposed algorithm frameworks, thereby highlighting the unifying viewpoint offered by our analysis. The use of the kernel K2 for SVGD in variational inference dates back to [54]]. Algorithms and with the kernel K2 correspond precisely to the GF and GPF algorithms in (27) Moreover, K is chosen, Algorithm [ reproduces the BWGD algorithm in 43 (with N = 1) and shares some similarity with the FB-GVI algorithm in [19]. Detailed discussions on these variants, connections and convergence properties are deferred to Appendix D

---

5000
58.20
25.
Sw6D 06
54000 16
17.5
s0r 203
53.5
TGP 36.
15.9-
130
2009
12.5
12.9
16.
5.-
0.0
2.5
5.0
7.5
10.0
12.3
15.0
37.
23.0
252
300
750
1000
1290
2500
1757
2002
Lime
(a) F(po) over time.
(b) F(po) over iteration.
Figure 1: Convergence of Algorithms 1 and with bilinear kernels in Bayesian logistic regression.
5 Simulations
In this section, we conduct simulations to compare Gaussian-SVGD dynamics with different kernels and the performance of the algorithms mentioned in the previous section. We consider three settings, Bayesian logistic regression, and Gaussian and Gaussian mixture targets. Here we present the results for Bayesian logistic regression as involves non- non-Gausiian but unimodal target and is one of the typical setups such that GVI is preferred in practice. For sake of space we leave the simulations for the other two settings along with further discussions to Appendix
Bayesian Logistic Regression. The following generative model considered: Given a parameter
E Rdd, draw samples {(X;, Y.) (Rd 1))" such that X; lid. No Ia) and Y, X; ~ Bern(a((E,. X.))) where (.)) is the logistic function. Given the samples {(X;, ) and a uniform (improper) prior on E. the potential function of the posterior o² E is given by Ei=l (log(1 + exp((E, Po (O. Ia) to find the pa- that minimizes KL(pp || p²). Figure SBGD (Simple Bilinear Gradient Descent), GF (Gaussian Flow (27)). BWGD (Bures-Wassersee
VEE) = >
(x))), We run both Algorithms/ I and initialized at
Gradient Descent [43]), and RGF (Regularized Gaussian Flow) are density-based algorithms with the bilinear kernels K₁, K2, K3. and K4 (v 0.5) respectively; SBPF (Simple Bilinear Particle Flow), GPF (Gaussian Particle Flow (27)) BWPF (Bures Wasserstein Particle Flow), and RGPF (Regularized Gaussian Particle Flow) are particle-based algorithms with the bilinear kernels K1, K2. K3, and K₃ (v = 0.5) respectively. We use the sample estimator of F(po) = -J [[log Pa + V) dpg = KL(pe p²) + (C is some constant) to evaluate the learned parameters 0 (J,, Σ). For a fair comparison in Figure Tb we draw 2000 particles for particle-based algorithms and run all algorithms for 2000 iterations so that total of 2000 samples are drawn for the density-based algorithms. The largest safe step sizes are 0.02, 0.1, 2, 0.8, 0.02, 0.2, 4,
Figure shows the decay (pe) over time or iterations. For Figure la the same step size 0.01 is specified for all algorithms while for Figure [b] we choose the largest safe step size for each algorithm. In other words, Figure Ta provides the continuous-tima flow the dynamical system in each algorithm while Figure Ib emphasizes more on the discrete-time algorithmic behaviors. In Figure Ta there are roughly four distinct curves indicating four different kernels. K2 has the most rapid descent, followed by K2. K4, and K3.
From Figure Tb observe that BWPF and RGPF are the better choices for practical use. The difference in the largest step sizes shows that in terms of stability K3 is the best, K4 is almost as stable, but K2 and Kj are much worse. We also observe that particle-based algorithms are consistently counterparts (which are essentially stochastic gradient based). The
more stable than density- superiority of particle- based algorithms are even more evident in the Gaussian mixture experiment where the target is multi- modiil..d.......m We further remark that another recently proposed gensity-based algorithm, the FB-GVI [19] shows comparable performance to BWPF and RGPF with large step sizes. We conduct comparison of these three algorithms in Appendix[F]g] do not include it here for clarity. It would also really interesting to study the particle- bartecle-based analogue FB-GVI as future works.
10

---

Acknowledgements
We thank Lester Mackey and Jiaxin Shi for several clarifications about their work, 722. and for helpful discussions regarding the larger literature on SVGD. Promit Ghosal was supported in part by NSF grant DMS-2153661. Krishnakumar Balasubramanian was supported in part by NSF grant DMS-2053918, Natesh S. Pillai was supported by ONR grant N00014-21-1-264..
References
[1] Pierre Alquier and James Ridgway, Concentration of tempered posteriors and of their variational approximations, The Annals of Statistics 48 (2020), no. 3, 1475-1497.
[2] David Alvarez-Melis, Yair Schiff, and Youssef Mrouch, Optimizing functionals on the space of probabilities with input convex neural networks, Transactions on Machine Learning Research (2022).
[3] Michael Arbel, Anna Korba, Adil Salim, and Arthur Gretton, Maximum mean discrepancy gradient flow, Advances in Neural Information Processing Systems, vol. 32, 2019.
[4] Jimmy Ba, Murat A Erdogdu, Marzyeh Ghassemi, Shengyang Sun, Taiji Suzuki, Denny Wu, and Tianzong Zhang, Understanding the variance collapse of svgd in high dimensions, International Conference on Learning Representations, 2021.
[5] Julio Backhoff- Backhofft-Veragas,, Joaquin Fontbona, Gonzalo Rios, and Felipe Tobar, Stochastic gradi- descent in Wasserstein space, arXiv preprint arXiv:2201. 2idi:22201.0222222 (2022).
[6] Rajendra Bhatia, Tanvi Jain, and Yongdo Lim, On the Bures-Wassersie positive definite matrices, Expositiones Mathematicae 37 (2019), no. 2, 165-191.
distance between
[7] Rajendra Bhatia and Peter Rosenthal, How and why to solve the operator equation AX-XB =
Y, Bulletin of the London Mathematical Society 29 (1997), no. 1, 1-21.
[8] David M Blei, Alp Kucukelbir, and Jon D McAuliffe, Variational inference: A review for statisticians, Journal of the American statistical Association 1 112 (2017), no. 518, 859-877.
[9] Bobkov and M Ledoux, dimediimes empirical measures, order statistics. and Kan- torovich transport distances, Memoirs Amer. Math. Soc (2019).
[10] José Antonio Carrillo, Katy Craig, and Francesco S Patacchini. A blob method for diffusion, Calculus of Variations and Partial Differential Equations 58 (2019), 1-53.
[11] Anthony Caterini, Rob Cornish, Dino Sejdinovic, and Arnaud Doucet, Variational inference with continuously-indeeee normalizing flows, Uncertainty in Artificial Intelligence, PMLR. 2021, pp. 44-53.
[12] Peng Chen and Omar Ghattas, Projected Stein variational gradient descent, Advances in Neural Information Processing Systems 33 (2020), 1947-1958.
[13] Yifan Chen, Daniel Zhengyu Huang, Jiaoyang Huang. Sebastian Reich, and Andrew M Stuart, Gradient flows sampling: Mean-field models, Gaussian approximations and affine invariance, arXiv preprint arXiv:2302.1 1024 (2023).
[14] Badr-Eddine Chérief-Abdelle Pierre Alquier, and Mohammad Emtiyaz Khan, A generaliza- tion bound for online variational inference, Asian Conference on Machine Learning, PMLR,
2019, pp.
[15] Sinho Chewi, Thibaut Le Gouic, Chen Lu, Tyler Maunu, and Philippe Rigollet, SVGD as a kernelized Wasserstein gradient flow of the chi-squared divergence, Advances in Neural Information Processing Systems 33 (2020), 2098-2109
[16] Shui-Nee Chow, Wuchen Li, and Haomin Zhou, Wasserstein Hamiltonian flows, Journal of Differential Equations 268 (2020), no. 3, 1205-1219.
[17] Kacper Chwialkowski, Heiko Strathmann, and Arthur Gretton, A kernel test of goodness of fit, International conference on machine learning, PMLR, 2016, pp. 2606-2615.
[18] Gianluca Detommaso, Tiangang Cui, Youssef Marzouk, Alessio Spantini, and Robert Scheichl,
A Stein variational Newton method, Advances in Neural Information Processing Systems 31 (2018).
11

---

[19] Michael Diao, Krishnakumar Balasubramanian, Sinho Chewi, and Adil Salim, Forward- backward Gaussian variational inference via JKO in the Mures-Wasserstein space. arXiv preprint Siii:22304005399 (2023).
[20] Justin Domke, Provable smoothness guarantees for black-bor box variational inference, Interna- tional Conference on Machine Learning, PMLR, 2020, pp. 2587-2596.
[21] Justin Domke and Daniel R Sheldon, Importance weighting and variational inference, Advances in neural information processing systems 31 (2018).
[22] Andrew Duncan, Nikolas Nüsken, and Lukasz Szpruch, On the geometry of Stein variational gradient descent, Journal of Machine Learning Research 24 (2023), 1-39.
[23] Yihao Feng, Dilin Wang, and Qiang Liu, Learning to draw samples with amortized Stein variational gradient descent, Conference on Uncertainty in Artificial Intelligence, 2017.
[24] Nicolas Fournier. Convergence of the empirical measure in expected Wasserstein distance: Non asymptotic explicit bounds in IR arXiv preprint
Stt1:22000009222 (2022).
[25] Nicolas Fournier and Arnaud Guillin, On the rate of convergence in Wasserstein distance of the empirical measure, Probability theory and related fields 162 (2015), no. 3-4, 707-738.
[26] Oded Galor, Discrete dynamical systems, Springer Science & Business Media, 2007.
[27] Théo Galy-Fajou, Valerio Perrone, and Manfred Opper, Flexible and efficient inference with particles for the variational Gaussian approximation, Entropy 23 (2021), no. 8, 990.
[28] Chengyue Gong, Jian Peng, and Qiang Liu, Quantile Stein variational gradient descent for batch Bayesian optimization, International Conference on Machine Learning, PMLR, 2019, pp. 2347-2356.
[29] Jackson Gorham and Lester Mackey, Measuring sample quality with kernels, International Conference on Machine Learning, PMLR, 2017, pp. 1292-1301.
[30] Jackson Gorham, Anant Raj. and Lester Mackey, Stochastic Stein discrepancies, Advances in Neural Information Processing Systems 33 (2020), 17931-17942.
[31] Tuomas Haarnoja, Haoran Tang, Pieter Abbeel, and Sergey Levine, Reinforcement learning with deep energy- based policies, International conference on machine learning, PMLR, 2017,
pp. 1352-1361.
[32] Ye He, Krishnakumar Balasubramanian, Bharath K Sriperumbudur, and Jianfeng Lu, Regular- ized Stein variational gradient flow, arXiv preprint arXiv:221 1.07861 (2022).
[33] Richard Jordan, David Kinderlchrer, and Felix Otto, The variational formulation of the Fokker- Planck equation, SIAM journal on mathematical analysis 29 (1998), no. 1, 1-17.
[34] Mikolaj Kasprzak, Ryan Giordano, and Tamara Broderick, How good is your Gaussian approximation of the posterior? Fimite-sample computable error bounds for a variety of useful divergences, arXiv preprint arXiv:2209.114 14992 (2022).
[35] Anya Katsevich and Philippe Rigollet. On the approximation accuracy of Gaussian variational inference, arXiv preprint arXiv:2301.0211 (2023).
[36] Gabriel Khan and Jun Zhang, When optimal transport meets information geometry, Information Geometry 5 (2022), no. 47-78.
[37] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling, Improved variational inference with inverse autoregressive flow, Advances in Neural Information Processing Systems 29 (2016).
[38] Anna Korba, Pierre- Peere-Cyrril Aubin-Frankowski. Szymon Majewski, and Pierre Ablin, Kernel Stein discrepancy descent, International Conference on Machine Learning, PMLR, 2021, pp. 5719- 5730.
[39] Anna Korba and Adil Salim, Sampling as first-order optimization over a space of probubility measures, Tutorial at the International Conference of Machine Learning, 2022, https:// akorba. ghu... io/ wttp//oooreee/alaaiiiiiooi.eeeee Btt.//aaattittteeeeeeeeeeeeeeeee saiiimor__ull22222 Ic0.tttaclt.t...aallll... pdff
[40] Anna Korba, Adil Salim, Michael Arbel, Giulia Luise, and Arthur Gretton, A non-asymplotic analysis for Stein variational gradient descent, Advances Neural Information Processing Systems, vol. 33, 2020, pp. 4672-4682.
12

---

[41] Andreas Kriegl and Peter W Michor, The convenient setting of global analysis, vol. 53, American Mathematical Soc., 1997.
[42] John D Lafferty. The density manifold and configuration space quantizatiinn Transactions of the American Mathematical Society 305 (1988), no. 2, 699-741.
[43] Mare Lambert. Sinho Chewi, Francis Bach, Silvère Bonnabel, and Philippe Rigollet, Variational inference via Wasserstein gradient flows, Advances in Neural Information Processing Systems, vol. 35, 2022.
[44] Thomas Laurent, Local and global existence for an aggregation equation, Communications in Partial Differential Equations 32 (2007), no. 12, 1941-1964.
[45] Michel Ledoux, On optimal matching of Gaussian samples, Journal of Mathematical Sciences 238 (2019), 495-522.
[46] Michel Ledoux and Jie- die-Xiang Zhu, On optimal matching of Gaussian samples III, Probability and Mathematical Statistics 41 (2021).
[47] Jing Lei, Convergence and concentration of empirical measures under Wasserstein distance in unbounded functional spaces, Bernoulli 26 (2020), no. 1, 767-798.
[48] Wu Lin, Mohammad Emtiyaz Khan, and Mark Schmidt, Fast and simple ational inference with mixture of exponential-famie approximations, International Conference on Machine Learning, PMLR, 2019, pp. 3992-4002.
vari-
[49] Chang Liu and Jun Zhu, Riemannian Stein variational gradient descent for Bayesian inference, Proceedings of the AAAI Conference on Artificial Intelligence, vol. 32, 2018.
[50] Chang Liu, Jingwei Zhuo, Pengyu Cheng, Ruiyi Zhang, and Jun Zhu, Understanding and accel- erating particle-based variational inference, International Conference on Machine Learning, PMLR, 2019, pp. 4082-4092.
[51] Qiang Liu, Stein variational gradient descent as gradient flow, Advances in neural information processing systems, vol. 30, 2017.
[52] Qiang Liu, Jason Lee, and Michael Jordan, A kernelized Stein discrepancy for goodness-of-P tests, International conference on machine learning, PMLR, 2016, pp. 276-284.
[53] Qiang Liu and Dilin Wang, Stein variational gradient descent: A general purpose Bayesian inference algorithm, Advances in neural information processing systems, vol. 29, 2016.
[54]
Stein variational gradient descent as moment matching, Advances in Neural Informa- tion Processing Systems, vol. 31, 2018
[55] Xing Liu, Harrison Zhu, Jean-Francoif Ton, George Wynne, and Andrew Duncan, Grassmann Stein variational gradient descent, International Conference on Artificial Intelligence and Statistics, PMLR, 2022, pp. 2002-2021.
[56] Xingchao Liu, Xin Tong, and Qiang Liu, Sampling with trusthworthy constraints: A variational gradient framework, Advances in Neural Information Processing Systems 34 (2021), 23557- 23568.
[57] John Lott, Some geometric calculations on Wasserstein space, Communications in Mathematical Physics 277 (2008), no. 2, 423-437.
[58) Jianfeng Lu, Yulong Lu, and James Nolen, Scaling limit of the Stein variational gradient descent: The mean field regime, SIAM Journal on Mathematical Analysis 51 (2019), no. 648-671.
[59] Jan R. Magnus and Heinz Neudecker, Matrix differential calculus with applications in statistics and econometrics, Wiley Series in Probability and Statisties, John Wiley, 1988.
[60] Luigi Malagø, Luigi Montrucchio, and Giovanni Pistone, Wasserstein Riemannian geometry of Gaussian densities, Information Geometry 1 (2018), 137-179.
[61] Gael M Martin, David T Frazier, and Christian P Robert, Approximating Bayes in the 2 l century, Statistical Science 1 (2023), no. 1, 1-26.
[62] Boris Muzellec and Marco Cuturi, Generalizing point embeddings using the wasserstein space of elliptical distributions, Advances in Neural Information Processing Systems 31 (2018).
[63] Yurii Exgenévich Nesterov, A method of solving a convex programming problem with conver- gence rate ) Doklady Akademii Nauk 269 (1983), no. 3, 543-547.
13

---

[64]
Manfred Opper and Cédric Archambeau, The variational Gaussian approximation revisited, Neural Computation 21 (2009), no. 3. 786-792.
[65]
Felix Otto, The geomelry of dissipative evolution equations: the porous medium equation, Comm. Partial Differential Equations 26 (2001), 101-174.
[66]
Peter Petersen, Riemannian geometry, vol. 171, Springer, 2006.
[67] Matias Quiroz, David Nott, and Robert Kohn, Gaussian variational approximations for high-dis high--iimnsional state space models, Bayesian Analysis 1 (2022), no. 1, 1-28.
[68] Danilo Rezende and Shakir Mohamed, Variational inference with normalizing flows, Interna- tional Conference Machine Learning, PMLR, 2015, pp. 1530-1538.
[69] Adil Salim, Anna Korba, and Giulia Luise, The Wasserstein proximal gradient algorithm, Advances in Neural Information Processing Systems (H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, eds.), vol. 33, Curran Associates, Inc., 2020, pp. 12356-12366.
[70] Adil Salim, Lukang Sun, and Peter Richtarik, convergence theory for SVGD in the population limit under Talagrand's 's inequality TI, International Conference on Machine Learning, PMLR, 2022, pp. 19139-19191
[71] Jiaxin Shi, Chang Liu, and Lester Mackey, Sampling with mirrored Stein operators, International Conference on Learning Representations, 2022.
[72] Jiaxin Shi and Lester Mackey, Finite-Particle Convergence Rate for Stein Variational Gradient Descent, arXiv preprint arXiv:2211 1.09221 (2022).
[73] Vladimir Spokoiny, Dimension free non-asynponn asymppotii bounds on the accuracy of high dimensional Laplace approximation, arXiv preprint arXiv:2204.1 1038 (2022).
[74] Weijie Su, Stephen Boyd, and Emmanuel Candes, A differential equation for modeling Nes- terov's accelerated gradient method: theory and insights, Advances in neural information processing systems, vol. 27, 2014.
[75] Lukang Sun, Avetik Karagulyan, and Peter Richtarik, Convergence of Stein variational gra- dient descens under a weaker smoothness condition, International Conference on Artificial Intelligence and Statisties, PMLR, 2023, pp. 3693-3717.
[76] Amirhossein Taghvaei and Prashant Mehta, Accelerated flow for probability distributions, International Conference on Machine Learning, PMLR, 2019, pp. 6076-6085.
[77] Asuka Takatsu, Wasserstein geometry of Gaussian measures, Osaka Journal of Mathematics 48 (2011), no. 1005-1026.
[78] Linda SI. Tan and David J Nott, Gaussian variational upproximation with sparse precision matrices, Statistics and Computing 28 (2018), 259-275.
[79] Gerald Teschl, Ordinary differential equations and dynamical systems, Graduate Studies in Mathematics 140 (2000), 08855-8019,
[80] Marcin Tomczak, Siddharth Swaroop, and Richard Turner, Efficient low rank Gaussian varia- tional inference for neural networks, Advances in Neural Information Processing Systems 33
(2020), 4610-4622.
[81] Aad W van der Vaart, Asymptotic statistics, vol. 3, Cambridge University Press, 2000.
[82] Jesse van Oostrum, Stures-Wasserster geometry for positive-definite Hermitian matrices and their trace-one subset, Information Geometry 5 (2022), no. 405-425.
[83] Cedric Villani, Topics in optimal transportation, Graduate Studies in Mathematics 58, American Mathematical Society, 2003.
[84] Dilin Wang, Ziyang Tang, Chandrajit Bajaj, and Qiang Liu, Stein variational gradient descent with matrix-valued kernels, Advances in neural information processing systems 32 (2019).
[85] Yifei Wang, Peng Chen, and Wuchen Li, Projected Wasserstein gradient descent for high- dimensional Bayesian inference, SIAM/ASA Journal on Uncertainty Quantification 10 (2022), no. 4, 1513-1532.
[86] Yifei Wang and Wuchen Li, Accelerated information gradient flow, Journal of Scientific Computing 90 (2022), 1-47.
14

---

[87] Lantian Xu, Anna Korba, and Dejan Slepcev, Accurate quantization of measures via interacting particle-based optimization, International Conference on Machine Learning, PMLR, 2022,
pp. 24576-24595.
[88] Rentian Yao and Yun Yang, Mean field variational inference via Wasserstein gradient flow, .08074 (2022).
arXiv preprint
[89] Jingwei Zhuo, Chang Liu, Jiaxin Shi, Jun Zhu, Ning Chen, and Bo Zhang. Message passing Stein variational gradient descent, International Conference on Machine Learning, PMLR, 2018, pp. 2-18-607.
15

---

A Further Discussion
Other Related Works. There exist several other approaches for minimizing KL-divergence over the Wasserstein space (or over appropriately restricted subsets). For example, normalizing flows [68] the blob method [00,, variants of gradient descent in Wasserstein spaces [69
37
natural-
gradient based variational inference techniques [48], mean-fiele methods [88], neural networks based approaches [2] and MCMC methods. This is certainly not comprehensive list and summarizing the large literature on this topic is impossible in the limited space available. However, we emphasize that a majority of the above methods are nonparametric and do not come with strong theoretical guarantees. Our main focus in this work is on theoretically understanding (Gaussian)-SVCh and its connection to GVL.
Beyond Gaussian VI. While we present comprehensive study of Gaussian-SVGD, it is currently unclear how such relation between bilinear kernels and Gaussian family could be generalized. For example, it would be interesting to ask if there is a class (submanifold) C of probability densities such that for any initialization-tain -targtt pair Pu. p² E C, the dynamics of SVGD with a radius- kernel function (RBF-SVGD) would always remain in C. If such C is identified, we could similarly carry out the convergence analysis of RBF-SVGD and perform variational inference with respect to the class C. It remains open question whether the uniform in time propagation of chaos as shown in Theorem 3.7 would hold for general kernels. Moreover, there has been increasing interest on variational inference with respect other special density classes. Remarkably (43 considered Gaussian mixtures in additional to GVI, and [62] analyzed the family elliptical distributions, which has applications for performing variational inference under heavy-tails [21].
SVGD in High Dimensions. severe particle degeneracy in high dimensions and cannot guarantee good covariance estimation. However, such issue has not been observed for Craussian-SVGG in our simulations. It would be interesting future works to study whether this undesired phenomenon is tied to specific choices of kernels and to earry out theoretical analysis of (Gaussian-SSVGD in high dimensions.
well-known in literature 89 that RBF-SVGD suffers from
B Details on Gaussian-Steie Metrics
For any region Ω Rd. denote the set of probability densities on Q by
P(L) = {0 E F(!)) In dd = 1, 1, 0),
where F(2) is the set of - smooth functions on Q. As studied in [42). P(22) forms a Fréchet manifold called the density manifold. For any "point" p E P(2).. the tangent space is given by
T,P(O) fo E F(D):
dx = 0}.
And the cotangent space at p, P(?) consists of equivalent classes of F(!), each containing functions that differ by constant. A Riemannian metric (tensor) assigns each p P(?)) a positive definite inner product 9p (T)( T(P( + R. If we define the pairing between P(?) and TS?( by
(2, 0) = f f dæ,
where 4- P(?) and a E TTP()). Any Riemannian metric uniquely corresponds to isomorphism (called the canonical isomorphism) between the tangent and cotangent bundles [66]. i.., we have Gp T, P(?( > P(?( through
9o (2) (G,o1, (1)0 01.02 T2P(S))
Definition B.1 (Wasserstein metric). The Wasserstein metric is induced by the following canonical isomorphism Gwass T(P(S)) - TT such that
Wass - - (oVo), & PTT).
Note that the Wasserstein metric we define here corresponds exactly to the Wasserstein-2 distance studied in the literature of optimal transport 83 For details of such connection, please refer to (36) Also we need to clarify that the concepts we introduce follow from the convention in Riemannian geometry, where the manifold strictly speaking is required to be finite-dimensiona For
16

---

mathematically formal formulation of refer 42
calculus on the density manifold, please
65 57
The Wasserstein gradient flow can be seen as the natural gradient flow on the density manifold with respect to the Wasserstein metrie. In specifie, consider any energy functional E(p), e.g., the KL divergence from p to some fixed target p.. The Wasserstein gradient flow for E(p) is provided by
OE 5p:
L Öp:
(G
Wass
V
Pr
(27)
where 6E is the variational derivative of the functional with respect to Pt- If E(p) is the KL ap; divergence mentioned above, (27) gives the linear Fokker-Planck equation (33]
= (Vp: + VVV),
where ((x) - log p (z) is called the potential function. E(p) is the Wasserstein metric between and p², ( gives the geodesic flow the density manifold. If E(p) is the maximum mean discrepancy (MMD) or kernelized Stein discrepancy (KSD) between and p² it leads to the MMD descent and KSD descent algorithms 38)
Interestingly, the mean field flow of SVGD can also be seen as the gradient flows of the KI. divergence on the density manifold but with respect to the Stein metrie, where perform kemelization in the cotangent space before taking the divergence [51] (2).
Definition B.2 (Stein metric). The Stein metric is induced by the following canonical isomorphism Wass P(12) + P(2) such that
(a()) (())
PG S -1 -
K( ((.:))ol))( & E TT
In particular, the mean field PDE of the SVGD algorithm can be written as
(c(x)= -CGSscinn 8 KL P1)))
(x)
(28)
8p:
V (p((z)) / ((a,v)p(())
KT I1
(y) dy
5pt
(polac) / K(E (((e,v)p(()) V v(og pi (v) + V + +1) dy
V. (pa(ar)
( dy
where V(x) log p" (x)).
We set Ω = R and consider the multivariate Gaussian densities . (M, Σ) in Rd:
p(x, 0) /det(2=>))
exp (-
1 (x 11) T Σ -ー' pu))
where 0 (p, Σ) E θ and == Rdd x Sym (d, Here Sym + (d, IR) the set of (symmetric) positive definite d x d matrices, which is an open subset of the x symmetric matrix space (ym R). In this way. can identified as Riemannian submanifold of the density manifold P(Rd) with the induced Riemannian structure.
We first look into the Stein metric with the bilinear kernel K1 (1 y) = x y + 1, and derive the closed form of the induced metric tensor on O, which plays an essential role in characterizing the SVGD dynamics.
T
Theorem B.3 (Gaussian-Saee Stein metric with the simple bilinear kernel), Given (n, E) E O. let (,( denote the Gaussian-Steni metric tensor for the multivariate Gaussian model with the bilinear kernel K₁ (x, y) = x ITy + 1, and Go be the corresponding canonical isomorphism from T,O ~ Rdd x Sym(d, R) to T Rd x (d, R).
For any V E R", and Sym(d, R), the inverse map G is given by the following automorphism on Rdd Sym(d, RR);
(v, = (28Eμ + (1 2 4)v, (Z(2SS + pv + (2SE + Pu TE))
(29)
17

---

And for any E, E Tgo the Stein metric tensor can be written as
go(S. n) - tr(S,, + (b) S2 + b
(1+ 1 pbb, b2.
(30)
Here { = 11, Σ1 and 1 = (12 Ez) in which P.. E IR., E1 En E Sym(d, R). And b,, S; (i 1,2) are defined as
bi: 5.)) = Go( E.).
Then for the Gaussian-See Saassiann-teein metric with kernel (2 y) = (xx Ty - (1) + we have the following result:
Theorem B.4 (Gaussian-Stecd metric with the affine-invariant bilinear kernel). Given A = (n,, Σ) O let go(-.)) denote the affine-invarian Gaussian-Sted
metric tensor for the multivariate Gaus- sian model with the affine-inarrant bilinear kernel K2, and Ge be the corresponding canonical isomorphism from To I& x Sym(d, IR) to T θ = Idd x Sym(d, R).
For any V E R", and Sym(d, R), the inverse map Gg is given by the the following automorphism on R2d Iy):
G (v) (v. (E S+ SE²)).
(31)
And for any E,n E TyO the Stein metric tensor can be written as
go( tr(S + i M2.
(32)
Here = (م (p En and 1 (ïz E2). in which i B2 E Rd. E, E Sym(d,R) And St's (i 1,2) are defined as the symmetric solution of
+ -
Now if we further restrict the Gaussian family to have zero mean, it induces the submanifold O Sym + (d, IR). We have the following result for K₃ or K2 as a direct corollary Theorem B.3
or Theorem B.4]
Corollary B.5. Given Σ E On let gz(,,)) denote the CauusiaanSttein metric tensor for the centered Gaussian family with the bilinear kernel K (x, y) x= + 1, and G5 be the corresponding canonical isomorphism from TgQo = Sym(d, IR) O Sym(d, R).
For any S E Sym(d, IR), the inverse map G 1 is given by the following automorphism Sym IR):
( 2 + S
(33)
And for any Σr, Eg E TyOo the Stein metric tensor can be written as
g Σe) = tr(5)S,2²).
(34)
where = 1,2, St is the unique solution in Sym(d, R) that satisfies the Lyapunov equation
E S.y2.
Next we consider the Bures-Waie Wasserstein metric for Gaussian families, which is defined as the Wasserstein metric restricted to the Gaussian family. It has the following elegant expressions from
77 60]: :
Theorem B.6 (Bures Wasserstein metric). Given 0 EH,E) O, let go(-..) denote the Bures- Wasserstein metric tensor for the multivariate Gaussian model (or equivalently Gaussian-Stee metric with the kernel K3), and Go be the corresponding isomorphism from T6O Rdd x Sym(d, R) to Ty ~ Rdd x Sym(d, R.). For any v E TR",, and S E Sym(d, R), the inverse map Go given by the following automorphism on Bdd Sym(d, R):
(G (v, (v, 2(ES SΣ))
(35)
And for any E,n E ToQ the
metric tensor can be written as
go(E, m) = tr(S, + j P2.
(36)
18

---

Here ( and rj
(2n in which P1, f2 Rd, En, Eg e Sym(d, And S
(i 11 are defined us the symmetric solution of
Notably this metric coincides with Gaussian-Steie metric with the kernel Kz(æ,))
(x:
14) Σ 1D'E-'(0 (y- 11) +1
Note that the only difference between 32) and (36) the power of Σ.
Theorem B.7 (Regularized Gaussian-Stein metric with the affine-invariant bilinear kernel), Given 0 (H, 2) E θ, let go(:;)) .) denote the affine-inwer affine-invariant regularized Gaussian- -teinn metric tensor for Gaussian families, and Ga be the corresponding isomorphism from ToO = -2² x Sym(d, R) to Th θ = R2d Sym(d, IR). For any u E Ted and S E Sym(d, R), the inverse map G is given by the following automorphism on Idd x Sym(d,R) R):
(G (v,
DI - S((1 - DEΣ -(11)))= ( (v,
(37)
Notably this metric coincides with Gaussian-Sar Stein metric with the kernel K1(1, y) (x 44) - ((1 2) +2) (y- 11) +.
C Properties of SVGD Solutions with the Simple Bilinear Kernel
58 showed a few nice properties of the mean field PDE ( for SVGD with radius-based kernels that can be written K (x= y) which is symmetric and positive definite, meaning that
ΣΣ K(æ; - x,)SEK 2 O, Va, E Rdd, & E,R, m M.N.
Although their results covered large class of kernels commonly used in practice e.g., Gaussian kernels, they do not apply to the bilinear kernel (x, y) = x + 1. In this subsection we establish similar results for the bilinear kernel. In fact, some of their results for radius-base kernels still hold here while some do not.
Before showing the properties of the mean field PDE, we need the following assumption on the potential function V.
Assumption C.1.
The potential function V I IR satisfies the conditions below:
1. V E C (Id)), V 2 O, and V V(x) co as |æ||| + co.
2. For any a, 3 > O, there exists a constant Ca., >> O such that if [yll| S a||a + 8, then the following inequality always holds that
(1+ VVV + (1 I VV Ca. V(x)).
To make things precise although all vector norms in Rd and all matrix norms are equivalent, always choose the Euclidean norm for vectors and the spectral norm for matrices unless otherwise specified. Note that assumption here is little different from Assumption 2.1 in (58). We do not require their second formula but our second piece assumption is slightly stricter than their third one. It is straightforward check that any positive definite quadratic form satisfies all the assumptions above, corresponding to the case where the target is non-degeneaate Gaussian distribution.
We use Pv (TR:) and (P(RR) (p 1......) to denote the set of probability measure μ on Rd satisfying
|elle, = ( (1 Vxx)) du(x) 8
(38)
and
|||||p == f ||0|| dp(x) <
(39)
respectively.
19

---

Theorem C.2 Well-posednned and regularity of the mean field solution). Let V satisfy Assump-
tion
C.1
For any E Pv (R there is unique pt C( o)), , Pv (B.)) which is a weak solution
10 with initial condition Po v. Moreover, there exists C1 > 0 depending on V such that t t2
llpeleey ScClll
(40)
I E (R*) n (PR"), then for any t E o0), we have that P E PP((R n (V(RR) and there exists O2 > 0 depending on V such that
120. has a density po 2 0, then pt also has density. Furthermore, po (RR)) for some ki then have Pt E 7t (R*). Here
llellpe S [
(41)
(R= Wh
D (P((R) vlal Sk),
k21
denotes the Sobolev (Hilbert) space of order k. Theorem C.3 (Well-posedmng posedness of the finite-paris inite-particle solution), Let V satisfy Assumption C.1 Then
for any initial condition Xo = (x ,. E RdN, the system
13 has a unique solution
(t)
T N E c oo), (0),
(1)
X; = 1
and the measure 1
N 8 d!
is weak solution to the PDE
Finally, we show that if two initial probability measures are close to each other, the solutions of (2) up to time T are also close. We need to further impose the following assumption on V. Assumption C.4. There exists constant Cv and some index q such that [V(x)) Cv + V(z))
for every x E Rd
(42)
and that
V(x) + V(y) [[|=e|
sup V (0x + (1 (1- (1-00)) 1 Cv 1 ec(o..]
(43)
|
Remarkably Gaussian distribution satisfies both Assumptions observation is that (42) implies that there is Co such that
C.I and C.4
Secondly an important
V(w) < Co ||x) Vae Rd
(44)
where of = =q 1). Indeed, we note that
VV
1/9 I 1- 1- C./o.
0 1++1 V
[
q
t
Integrating from q/(g- -
| we get
44 which shows that Pp C Py for any p d'
Theorem C.5. Let V satisfy Assumptions C.I and C.4
Let R> 0. Assume that vı, V2 two initial probability measures in (RR.) satisfying |iller & (i 1,2). Let pra and p2, be the
associated weak solutions to (2) Then given any T > 0, there exists constant Cr >> 0 depending
on V, R and T such that
sup Wp (PI. p2,) S Cr (21,22).
tElorr
Remark.
1. Theorem C.5 implies the convergence of empirical measure to the mean field limit time E [0,T) T7. In fact, if we Vi,n be an empirical measure, which converges v2 as grows to infinity, then p1.n3, the solution of (2) at time with initializaiton Vi. will also converge to p2.t. any E [0, T}.
2. In general there is guarantee that the density Pt will converge the target density p* in ( with the bilinear kernel. Some counterexamples will be elaborated in the remarks following Theorem 3.6
3. We show in Theorem 3.1 that for Gaussian families Pr. always converges to the target density p² as + co and the convergence rate is linear in divergence. We also establish uniform in time convergence result (Theorem 3.7 for the empirical measure in this case.
20

---

D Details of the Gausiian-SGGD Algorithms
Different Ways of Estimating Ft. The first-order estimator of T' arises from
Te = po [+) po(x) dx =- - v pa Vpo(x)VV((x) dæ
(a(aa) (m-))Va\) dæ E- E - [1)(a))
Since I is symmetric we also have
Ep (VV()) -=
Note that using the first-order estimator also comes at cost as the inverse Σ is needed. However, for density-based Gaussian-SVGD this E-1 might cancel with Σ in computing.
Previous Algorithms Under the Proposed Frameworks. The use of K2 for SVGD variational inference dates back to [54]. Our Algorithm 2 slightly different from [54] in the sense that VV is replaced by a linear function to ensure Gaussianity, Moreover, Algorithms T and [ with the kernel K2 correspond precisely to the GF and GPF algorithms in 27 If K3 (Bures-Wane Wasserstein metric) chosen, Algorithm reproduces the BW-SGD algorithm in (with N 1). K3 BBures-Wasseree
19] also uses metric) but their energy function for gradient flow is different from others. Instead of directly performing the gradient descent to minimize KI. divergence, they separate the KL divergence into two parts and perform the proximal gradient descent.
Variants of Gradient Descent in Density-Based Gaussian-SVGI For the density-based SVGD, we draw new samples at each step. It is interesting to study the behaviorial difference between drawing one sample (stochastic) and efficiently many samples (almost deterministic). For example, in [43] only one sample is drawn at each time step and they study the stochastic properties arising from this design. [19 considers both settings. In general, they do not differ much in convergence rates but there could be huge gaps in the constants of the bounds and will actually impact practical performance. Another choice only draw N samples at time 0 and we use linear transformation of the same N points to serve new samples at time f, which becomes similar to the particle-based Gaussian-SVGG SVGD. Moreover, the vanilla gradient descent could also be replaced by accelerated ones or gradient descent with adaptive learning rate, e.g., AdaGrad, RMSProp.
Resampling Scheme and Particle-Level Convergence. As presented in the main text the Gaussian SVGD for general target is given by
N Σ
x x
(++1 (t)
Σ
K(x; (0)
(1)
10 K( "<aa"",)")vvvr (v) j=1
(45)
N
where (V(x) = fi M4) mc
The updating rules of Caussian-SVGR given above is totally deterministic, meaning that at each time step x²") is updated only using deterministic quantities and (a (!) k=1 without any external randomness. This is computationally more efficient but imposes difficulty in the analysis. On the other hand, we could also consider slightly modifying the updating rules by applying resampling scheme to get a better estimation of m, and Tt. In other words, we resample (yx) i.id. from (u,, Σe) and use the following estimators
M (V
M
1
Σ
m:
(46)
M
M
or first-order estimators
M Po
M
F
Σ
- p.) TT
1
Σ
(y)
-
(47)
M
M
k=
In this way, when M is large enough (V will be sufficiently close to l,(x m4) + m.
21

---

3
SBGD (Ki)
(K2)
BWGD (Ka)
RGF (Kg)
1-
(K)
GPF (Kg)
Kig lp-)
BWPF (Ks)
(Ka)
-2
-5
75
100
125
150
175
200
time
Figure 2: Convergence of Algorithms and
with bilinear kernels for a Gaussian target.
Now we replace m. and I by m. and r, and consider the continuous-time dynamics
N
(
K 1
)
(( (1) (0 (Te(()) 14) + mi)).
(48)
j=1 where be and Ng are sample mean and variance and
m, == En-N(n (VV(()) r I,
Theorem D.I (Equivalence of density-based and particle-based algorithms). The solution of the
finite-partich system 48 with K1 is given by a Ax mb) + ue where As is the unique solution of
- T,C, mn Ar,
Ao 11=1
where bh and En are the unique solution of the ODE system
f4: (T T.E+ M4 + p. p4 p.)m ン -2, (E,I' p.m..) (T-Ee + m.uu.) 2x
(49)
This can be proved using exactly the same technique as in the proof of Theorem 3.6 Here 449 is the
(22), and hence by Theorem 4.1 [4, and Eg converges to fa and Σ' that solves the GVI and the convergence rate is given in Theorem 4.2 when the target is strongly log-concave, We also conjecture that there is still uniform time convergence the mean-field limit for this particle system and leave it to future works.
E Details of the Simulations
Gaussian Targets. Following [19 we consider scenario where the target is Gaussian Mu, Σ)
where μ [ 1) and Σ-1 U diag ( >2,- Auo with U R 10x10 drawn from the Haar measure of the orthogonal matrices O(10) and λ₁ --, .-- - being geometric sequence such that 0.01 and Åιo 1. We run the eight different algorithms introduced in Section 5 and show the decay of log KL(pp || p²) over time in Figure 2 Clearly the algorithms with Ki show faster rate the same rate. This is actually confirmed from our theoretical analysis as in all the other dynamics except SBGD and SBPF, the mean converges at the rate of (ce-(AA) while the covariance converges at a faster rate, resulting in the fact that the KL divergence converges at O 2//A). But for K the rate is different and given in Theorem 3.1
time compared the others while the other algorithms eventually all converge
Gaussian Mixture Targets. Next we consider the 1-dimensional Gaussian mixture targets given by wnN(un + (1 w1 (H)2, oz). We run the aforementioned eight algorithms with initial
22

---

50G0
50.20
0.-
1GP 102
5.-
RC7 361
UC3
SBA (C)
-55
-20
-20
130
600
500
(a) F(po) over time.
(b) ((p)) over iteration.
Figure Convergence of Algorithms 1 and with bilinear kernels for a Gaussian mixture target.
a = 0 and = or particles drawn from (O, 1). Here again we plot the decay of F(pg) over time iteration shown in Figure In particular, the plots correspond to the specific setting of H1 = 5, P2 = 10, 01 5, d2 = 2 and p" (x) x 0.3 exp( - (x - 5)0/ + 0.7 exp(- (æ - 10))/8) These parameters are arbitrarily chosen. For the decay of F(po) over time, we fix the step size to be 0.1 and run 1000 iterations. For F(pe) over time, draw 500 particles for particle- baseee---bddddeee algorithms and run all algorithms for 500 iterations so that a total of 500 samples are drawn for the density- algorithms. The step sizes chosen to be the largest ones that still allow convergence. Specifically for these eight algorithms the step sizes are 0.02, 0.1, 1, 1, 0.2, 0.8, 8, 8. Consistent with the results of Bayesian logistic regression, the particle-based ones more stable and allow larger step sizes, with BWPF and RGPF clearly outperforming all the others. In fact, the constrast between particle-based and density-based algorithms particularly significant in this problem probably because of the non-log-concave target.
More the Bayesian Logistic Regression. We compare three so-far best performed algorithms. WWF,, RGPF. and FB-GVI [19] for the same problem with different step sizes. From Figure 4 that BWPF outperforms the other two. FB-GVT better than RGPF with larger learning rate but fluctuate a bit more when = 2. This is probably attributed to the stochastic gradients. Furthermore, it is interesting compare ordinary gradient descent (OGD) the variational parameters (mean and covariance) and SVGD with a radius based kernel function (RBF-SVGD)
(2, y) = exp -- - Firstly, OGD does not converge as fast as BWPF and it is not as stable. Secondly, RBF-SVGD is quite sensitive to the choice bandwidth and does not converge as fast as BWPF in general. We also notice that RBF-SVGD is significantly slower in computation compared to Gaussian SVGD. However, as a particle-based algorithm, RBF-SVGD does have the advantage of being stable even when the step size is large.
We show this comparison in Figure with the same step size 17
20
F Analogous Results for the Affine-Invarint Bilinear Kernels
For the Mures-Wassense metric (Gaussian-Stedl Stein metric with K3), the dynamics of natural gradient descent has already been studied in literature. See [60 86 13 for proofs for the following theorem.
Theorem F.1 (Wasserstein gradient flow for the Gaussian family). Let po ~ NN Eo) and p² ~ (b, Q) be two Gaussian measures. Then the solution of ([] converges to p* as - co. In particular, Pe is the density of N N(u,, Σx) where the mean lt and covariance matrix Σ, satisfies
P - (1r b)
(50)
Σ 21 -E Q-1
EoQ QEo, have ||m.: - b|l= Ole ((A) and |E Q11 = O(e-2t/A)) where X is the largest eigenvalue of Q.
Now we focus on the results for CoaussianSVGGI with the kernel K2. First remark that for K₂ the previous results on the well-posedness of mean-field PDE and finite-particle solutions in Appendix C
still hold and can be proved using similar techniques to Appendix H but for simplicity they
23

---

BWFF
20.0
BWPF In
BWPF In=
RGPF (n- 8)
17.5
(n-
15.0
FB-GVI (n
FB BB-GVI (n
12.5
FB-GWI (n 0.51
10.0
100
200
300
400
500
iteration
Figure 4: Performance of BWPF, RGPF, and FB-GVI with different step sizes 1).
BWPF
SVGD (h 1]
SVGD (b= 2]
SVGD (h 101
18
SVGD (b 201
14
100
300
400
500
Iteration
Figure 5: Performance of BWPF, OGD, and RBF-SVGD with bandwidth h.
omitted. The proofs of the following theorems will also be omitted unless different proof technique from the analogous results needs to be applied.
Theorem F.2 (Analogue Theorem 3.I.. For any 2 the solution Pt of SVGD 2 with the bilinear kernel remains Gaussian densiry with mean μ4 and covariance matrix E, given by
j4 (
(51)
ΣS =2221; -
- -2-10
which has unique global solution (0) given any pe R.d and Eo E Sym+ (y,+t(d,i (d, IR). And p converges weakly to p² as + oc. Ef = Q then we have the following rates
|| - bl = ( -t2). E2 Q1 = Ole Ole-2),
Ve > 0,
where X is the largest eigenvalue of Q.
The proof of the dynamics is similar to that of Theorem 3.1 The rate Ole (/e- is trivially true from the theory linear ODEs and the rate O(e 26 is given by Theorem 3.2
Theorem F3 (Analogue of Theorem 3.6). Suppose the initial particles satisfy that Co is non- There exists a unique solution of the finite particle system
13) given by
x
Ar(x; (0)
Mo) ps.
(52)
24

---

where A₁ is the unique (matrix) solution of the linear system
(- As, A = I,
(53)
and Me and C₁ are the unique solution of the ODE system
14. = - (14
(54)
C+ 14==2;; 2C, QiQ -0---Q-1- C
The proof is similar to that of Theorem 3.6
Theorem F,4 (Analogue of Theorem 3.7 Given the same setting as Theorem F.3
further suppose
the initial particles are drawn i.idd from N(uo Eo). Then there exists a constant C,.QUb,E such that for all I, for all N 2 2, with the empirical measure CN
Li=1 8 fo, the second moment of 2,
Wasserstein-2 distance between 2 and Pr converges:
N log log
if 1 if 2 if 2
W (5 p.) $ Ca.
N N- N)² N-2/d
(55)
The proof similar to that of Theorem 3.7 But for sake of completeness we provide more proof details in Appendix L
Theorem F.5 (Analogue of Theorem
Let p' be the density of a target distribution with the potential function V(x) = log p² (a) and Po be the density of N N( Eo). Then for any 2 0 the Gaussian-SVe SVGD produces Gaussian density Pt with mean μa, and covariance matrix Eg given by the following ODE system:
μi -Epppe [v))))
(56)
Σ = 2Σ, Ee (V=(()) (V
Furthermore, suppose that 0* the unique solution of the following oplimization problem
min KL(Po || p*), where Pe is the Gaussian measure N(u, E).
6=(,,E)
Then we have Pr + Pa- N(x"; Σ¹) as 1 cc.
The proof similar that of Theorem 4.1 particular, if the target strongly log-concave, it gives rise to the following linear convergence rate.
Theorem F.6 (Analogue of Theorem 4.2). Assume that the target p¹ is c--ttrongly log-concave and
3-log- smooth, Le., aT 3 V = 3T. Then Pt converges to po- at the following rate:
14 1 ((e at/
-at/ mix(
- =0 at/ max (8.1))
The proof is similar to that of Theorem 4.2
Theorem F.7 (Analogue of Theorem D.I.. The solution of the funit--pariiie system 48 with K2 is
given by x! (1) Ax (0) (o)) + nr where Ar is the unique solution of
A = (T- (T,CC Ao -
where ht and E, are the unique solution of the ODE system
j4 m+
8ラ 25, -
The proof is similar to that of Theorem 3.6 or Theorem D.I
25

---

G Proofs for Section B
2/2 (det(2)) (X)) (det(E))=112 dimensional normal random vector where is a positive definite matrix. Then for any real matrix A and d-dimensional real vector b, we have
Lemma G.1. Let p(z) (2x)-d/2
-1/2 exp (-1xTE-1z)
x) be the density of a d-
f Axp(x) dx = tr(A)))
f ('xxAp((z) do Alb.
Proof. Since E- is positive definite matrix, 5-1/22 and (x) (2m)- exp (- 2a (aTTar) æ). Then we have
find its positive definite root E- 1/2, Let
/ Aarp(ar) dar / /
521/2 / (det(E))) po(y) (det dy
12 A25 )on0
po(3) dy
12 13)
j=1
= =tr AΣ tr(AE).
The second equation is given by
f b" Amxp(()) dr = | |Aza' bp(ar) da
=A ( xx p(xx) dx b AEb.
Lemma G.2 (Lyapunov equation). The Lyapunov equation
PX - XP
has unique solution. If Q Sym(d, R), then the solution X Sym(d,
Proof. By Rylvester-Rosemblum theorem in control theory 77), PX XP Q has unique solution X. Note that if P, Q E (d, IR) and X a solution, then X is also a solution. Thus, we have X X which implies that E Sym IR).
Proof of Theorem B.3 Note that the tangent space each θ Rd Sym R) is ToO Rdd x Sym(d (d, R). If we define the inner product (pairing) on TpO by
(E. n) tr(E,E2) 1l P2
(57)
for any E, ToO, where (iz Σn) and n (iz, Σz), then the tangent bundle TO is trivial.
Since
P(Rd)
0
p(.,0)
provides an immersion from θ to P(Rd), we consider its pushforward dó given by
d
TeO
T,P(R*)
g H (Vop(:
26

---

On the other hand, for any d E (RR"), the inverse canonical isomorphism of Stein metric maps it
-v (((.))ol. T)P(R").
Thus, we obtain that
(Vop(z,0).E) ~ V...0 ((4)))))))))))
(58)
for any a R
Now we try to find suitable function & that satisfies the equation above. We compute that
(Vop(x, 0)
=tr(Vz 0) + vup(x,0) 1
- 1 tr ( (x - -1)- u)) - - (1) p(z,0) 0).
Letting (())= = f (xx, 4) 0) (3) dy, we we get
Ve p(x, / ((2 = Ve ()
= (V) E -x-(x - p1) V w(z)) p(x,
We choose (x)( (i 14) bj. where S E Sym R) and b R will be determined later. Note that S₁ needs to be symmetric because the gradient field is curl-free. We derive that
f (a y 1)p(y, (V)
(ar)
/ (a: y 1)p(), (S 11) + dy
Tyy / (x" (y + p1) 1) u. (Siy + b1) dy
La (x μ + 1)b1
=(S + bi )x + b
By comparison of the coefficients, we need
(5 - 6, 5 1) - (S ())++b)))))))) Σ- Ir(S + bi 1 Ar(E-1EE))
(59)
81Eμ + (1+ 1 pubi ih,
Note that the Lyapunov equation
PX + XP Q
where
P E 1 1 1 mn E.
Q Q= Σ 1 HTTT (Euii) + i E),
has a unique solution X = S E Sym(d, R). Together with
b =
14 (eg - S:2ip),
27

---

we find that (59) holds with the unique solution described above. Now from the calculations above straightforward check that (59) equivalent the following equation
G Go-(bi
1 51) (
Ei.:
where Go is the map defined that Go is isomorphism.
Thus, the existence and uniqueness of the solution indicates
Similarly, we let
( = 52(z - (4) + b2
where S2 Sym(d, R) the unique solution of the Lyapunov equation
PX+XP=O2
where
P E I
1 p u ui E,
Q= E2
+ μ (Epiia + T2z E.),
and
b = 1 + xe (it- - -S2Ep)
Now we compute the Riemannian tensor
go(E, n) 9p(5. = / (1,
/ (v)g(2)) p(x, 8) / (( yx 1) V (3) dy da
/ p( 0) (5) 11) + bi) ((52 + bap ) + b2) dæ
= tr (8) bz } + b (S + (1 μ b2)b))
-t((S, (b) S1) (1 f1 /1 [1)b] b2-
Finally, we show that Gy is indeed the canonical isomorphism corresponding to go(-.)). We check that
go(E. n) (S, bb4 )))) + 6T + (1+ p b2)
1
tr(S E2) + bj 0 (E)
Proof of Theorem B.7 Similar to the proof of Theorem B.3 we define the inner product on TyO by
(E,n) = tr(E + M1 P2
for any E E ToO where & - (74, Σn) and y = (72. Σ.) The map
+ + P(R")
p(-00,
where p(:, Σ) denotes the density of (p,, E), provides immersion from θ to P(IRO). We consider its pushforward dov given by
dΦo
ToO
T,P(R**
g
(oo(::0),E)
On the other hand, for any & E PRR(), the inverse canonical isomorphism of regularized Stein metric maps it to
-- TT + +v1 -Tx,p (V4) E T,P(RT)
28

---

Thus, we obtain that
(Vop(ar, 0) V (p(æ, ()(1 TTT,, - TK,, (Vo(x)))
(60)
for any R
Now we try to find a suitable function v that satisfies the equation above. We compute that (Vop(x, 0),E)
=tr(Vep(2, 0)Ei) + Vup(a, 0) T) 1
( Σ E (x-
(x-))²² --)²E= EE mi) +
+ E
(1) p(x, 8).
Letting V(x))= ((1 - u)Tk,p uT)-1Tkkp (V get that the RHS of (60) is equal to
V= pl., (x:) (v)z)) T (a: 11) v ((())) p(a,, E).
We choose (x) = S1 x - (1) bi where S1 E Sym(d, R) and b E Rd will be determined later. Note that S1 needs to be symmetric because the gradient field is curl-free. We derive that
1)p(y / ((x - 4) (y 1) ()V dy / ((x 11) (z- - p)) + 1)l(y, () μ) + b1) dy
T (Vbg(z))
=5 u)) +b.
Thus, we know (x) = 8,E (11 - U)Σ vl) -1 x + bj. need b1 ii, and
By comparison of the coefficients, we
Σ ( - +) "= -)) 330" tr
(61)
([1- (11
= -
Note that the first equation equivalent to the following Lyapunov equation
(1-= X ((1 - + - =
which has unique solution = S1 E Sym((, R), and the second equation is automatically satisfied once have the first one. Now from the calculations above it is straightforward to that
G (b,, = (b,, (((1 - =) + wz) - +
- (11= +
The proofs of Theorems B.4 and B.6 are similar to that of Theorems B.3 and B.7 In particular,
other proofs of Theorem B.6 also be found in literature, for example, 777 60
Finally,
Corollary B.5 is the direct corollary of Theorem B.3 or Theorem B.4
H Proofs for Section
Given probability measure μ and mocc--maasr map f, we denote by f the pushforward of the measure under the map f- Definition H.1 (Mean field characteristic flow). Given a probability measure v, we say that the map
xt) [0, co) x I Rd
a mean field characteristic flow associated the particle system or to the mean field PDE if X is C in time and solves the following problem
Xt, x, = (VK pr) (X v)) - (K VV)) (x(,
F4 X(t,,,v)+1
(62)
X (0, x, v) = x
29

---

Note that here (tt,,, : v) is the push-forward of y under the map H vV) and {X( v2) can be regarded as a family of maps from I to R., parameterized by t and v. In the lemma below we show that the mean field characteristic flow (62 is well-defined.
Lemma H.2 (Solution of the mean field characteristic flow). Assume the conditions of Assumption C.7 hold, and E (R.( For any I > O, there exists unique solution ( E C²¹ (10, [O, to the problem 62 where Y the function space given by
|2(x) x 1 1+ x
Y
u E C(TRR Rd): sup
#ER
Moreover, the measure 14+ X satisfies
1p I [e
for some constant that is independent v.
Proof. We prove the lemma in two steps. First we show local well well-posedness of the mean field characteristic flow. Second we extend the local solution to [0, oo).
Fix r > 0, and we define
|u(x)) x1 + |
Y = u u sup
S srj.
wERd
We show that there exists To such that (62) has unique solution X(1, 2r, v) E [0, Tol. and the solution the following function class
S. = C To), Yr),
which is a complete metric space equipped with the uniform metric
[u(x)) v(x))]
ds (u, v) == sup dy (u(t, (),, dy o)(
sup meRd
rE(U,2a]
1 |
Now we check the integral formulation of 62 given by
'/' V2K (X(,,,,,, u(dx))
X (1,, v) x
0
3d
1.' K (X((s,,, VV VV v(dx')
Rdd
S / VV (X((s,z', v) T X(s, v)u(da") ds. X(s,
(63)
Let us define the operator F u((,,)) F(u)(t,-) by
F(u)(t,)) = x 1' x() ds 1/'' / VV VV (u(s, u(ss u(s (dx') ds. JRd
We aim to show that F is a contraction map in S,, and thus, it has a unique fixed point. For this purpose we first prove that J maps S, into S,. It straightforward check that (t) Ft)) x). We now need to establish bound on (t. x If a S,, then for any [0, Tol and R
[|u(s, x) = ||0||| + [l)(s x ell s (r + 1) ||||
By Assumption C.I
have that
| (u(s, u(s, ( ||a(s
(u(s,
<(r + ((++1)(1+)) 1)(1 1|= VV (u(s,
<(r + 1)Cr+1, V(x)))
30

---

As consequence, we have
(t, =|| st((r + 1)|| + +r) + ( |)||=||| + (r + 1) f1+ (1+V(a)))(aa))
sŌ,t(1+ |x]]):
for some constant Čj, where we used the assumption that E Pv(RR). To I r/C, we get
Therefore, choosing
[F) x < Č,To
sup sup (E øERd
1+ |
which shows that F maps from S, to S, for sufficiently small To- Next, we prove that F is indeed contraction map. If u, u E S,, then for any t E [0, To) and a E Rdd
[(u)(t,x) - [t))
[ - v(s ds S.' J VVV (u(s, x's v(dx') [lu(s x))
1/"' { VV (u(s, x x2))| I)I 1|l((, - v(ss [lu(s x) ds
% VV (u(s, x))) VV (0((s,)))) |10(8, x')l [|[[, dsl I+ III IV.
Term can be upper-hounded by
1/(1 + x < /"
I ||u[s, x)) v(s,
ds S v).
(64)
1 |
Similarly we have
(1 + [|=|]) S tds(u, ur + 1)C,+1, / (1 L + (zd)))ddzz))
(65)
For the third term, we apply Assumption
C.T and get
|| V (u(s, x )) -|l(s,/
2(), || x0) v(ss x')
=(1 + || () VV (u(s,x')))
1 ||-|||
<(r + 1)C (1 (1+
Thus, we have
m/(1 + ||>||) <(r + 1)Cr+1, / ((r)))((dz))
J."
||0[s æ) ds
1 ||||||
<tds(u, cr + 1)2 / ((+
(66)
Finally applying Assumption C.I once again, we have
IV (u(s, - VV (v(s, x )) - (ee(( xx))
[[)])
=(1 + || VV (u(s. x²))
VV (o(s:
1 1 ||>|||
5(r+ + 1)(1 + || II)
max 1]
2 V(uu( x') + (1 A)o(s, x'))|| |lu(s, - v(s, x²)
lu(ss x') - v(s x')|
+(r 1)(1 +|
max AE
V²V(Au(s, x') + (1 ))u(s,
1 x<|
+
31

---

where in (*) we have used the fact that Au + (1 - A)e Sp, and thus, Au - (1 X)e also satisfies the inequality (3.3), which enables us apply Assumption C.I
Thus, Term IV is also bounded from above by (the same as the upper bound of Term III
IV/(1+ + | I]) S tdg(u, vr C) (1++((Z)))
(67)
Now combining (64 67 we conclude that F is a contraction on S, when To is small enough. By
the contraction mapping theorem, F has a unique fixed point (-...v) E Sy, which solves (63) After defining [0, Tol.
X (f, , v one sees that X (t, x, v) solves (62) in the small time interval
Now we proceed to the second step of extending the local solution. Define
T == sup (t E IR+ U [oo)
(62 has a (unique) solution on [0,
If oc, then have global solution. Otherwise suppose < co. After examining the bounds we have established in the previous step, we can see that supposing the local solution exists at some time To, it may be extended beyond To as long as the quantity
/ (1 2()
Me Hell
is finite at time To. We therefore establish an upper bound on this quantity.
/ (1 + ((X(((
X(s ((d(d)))u(da)) (8, X (V(X(s,,, (V
/ (1 ((
$C1, _
The last inequality follows from Assumption C.1 Therefore, by Grönwall's inequality we get [[elllyy
- (1+ ((X((, S Ci.ot [ (1 V(z))u(dx)) cec.o (1llpll
holds for all E T). Next we show an upper bound on [t, x v) We derive that
0, (1, 2)12 = 2X xx v) (t,
/
VV (X u) Xtt v) (Tlda')) x x(t, x)
= = 2x,X(t,
1+
Rdd
1 (1+ (3))
21,
1 Ci.0
xt,,x(,, (1 + Chooe
Again the first inequality follows from Assumption C.I Thus, by Grönwall's inequality we have [x x, w)ll I cxp (t + (cCr.at 1)Ilelpaa(e))) ||=|||
for all E (0.r). This also implies that
[| x)]l
< (1 + C1 Goot (100))) (x (X
< + Cloc C1.ot Grod (ull)))))) exp (t + Ce 1)(()(1Pa()))) ||~||,
32

---

for all E (0,r).
Next extend the solution to / E [0, For this purpose, first prove that for any sequence {4:) such that 0 < 4 < </2 <... < T and lim;- 4 = 7, the sequence of functions {X (4i, (,,;,)) is Cauchy in Y. Then by completeness of Y, there is limiting function for the sequence, and is straightforward to see that such function is unique (does not rely on the choice of the sequence (t))))
fact, we check that for any j
|IX x. v) - -X(t,,
=(1; - 14) 1. X(At, + +(1 (1-A)t;, x. 1. dA
<ty - t4) sup xtx
<(t - 44) (1+ oue ( exp (T + (e - (1)()))))(((4)) |x]||
Thus, for any e > 0 there exists N > O such that for any > i > N, we have
| (1,, a. v) X(1,, x,, 1/ll <6
(X), v), X (li, x, U)) sup
1 |<]|
In other words, (X(,,,,, is Cauchy sequence in Y.
Now know that 62 has unique solution [0, r]. Since (1 0O, can further find
unique solution of 62 on [T, 7 + Tol for some To small enough, which contradicts the definition of T. Therefore, we conclude that T = oc, and (62) has unique global solution. Finally, thanks to the integral formulation X is continuous on [0, oo) x Rdd The proof is complete.
Given v, let ((t, xr v) be the mean field characteristic flow defined in (2 and let P: = (t., vV Note that 2 can be rewritten as
Proof of Theorem C.2
Pe + (p.U(pt)) = O,
where U[p] is the vector field given by
VV(y) av(dy) - -f (())VV(y))' Tayyy aayy-
(p)((x)= x
Rd
Then Pr is a weak solution to ( in the sense that
sup ||p|| < <0, VT > 0,
tE[0,7]
and 1. J/a (Ó(1. az) + Vo(t, U((e(((()))) pe pe(dx)) di + I $10, ())xx))(dz)) =
holds for all C. (10, 00) Rd)) This either directly checked or follows immediately from Theorem 5.34 in [8
By Lemma
H.2 there exists some constant C1 such that
|| lp:llr ||l|ly.
Suppose that v E (RR) n Py (Rd), As shown in the proof of Lemma H.2 the map X(t, x, v) is an element of the space with dy (X (t, x) < Cit Therefore, since
x))l < 2| 2 |y (Xt,,
we have
| = [10llPr = | gd
(1, v) s c Cat|||u1ppr
for some constant C2 > and Pr (PRR) n Py (Rd).
33

---

We now explain that the uniqueness of the weak solution follows from the uniqueness of the mean field
characteristic flow. Suppose C (10, T Pv Pv (B. is another weak solution to ( By definition of
the weak solution, the vector field (t, x) + [9.] (x) is bounded over [0, T x R continuous in
and Lipschitz continuous in æ. Then we can define continuous family of maps (t,,, v) by
X = U (X) {0, 22
And the measure q =
(tt,,,))H1) -((,,, is weak solution to the transport equation
= 0
with initial condition go - -= 20- Uniqueness of the solution to this linear equation implies that = Qt- Thus, we have [(,, v) = Qt- In other words, (t,, æ. v) is the mean field characteristic flow for v. Uniqueness of the characteristic flow implies that X X, and hence qe - Ps- Thus, conclude that the weak solution is unique.
Lastly, show the regularity result: If has density po 0, then Ps also has density. Furthermore, if po E (IR") for some k, then we have Pt E H²(())
Note that we have already proven that Pr E (10, 27, Pv) and
letllp < |poolll t 820.
Noting that
U(p[(x) = x - / Vi(y)y a dply)
= + [ V(y))
/ (V(r) - du(
we have
U[A((()) S S [|lllll \VUp) s e ply
D& (g.()) 0 for
Thus, xt) (x) C TT.C (RR) where cb (Rd) is the space of continuous functions with bounded (k + 1) 1) order derivatives. Let (. X (1, xr denote the characteristic flow. Since + satisfies the ODE system
(x) U (pp
from the regularity theory of ODE systems (see Chapter 2 (79]) we know that both the map + 4; and its inverse & are C* Therefore, if pu has density, then pt also has a density and given by
(x(x) (9:) (o(0, exp (-f'
(V.VV--Uo..) Up,)) (z)( ds
Moreover, since p satisfies
p - (p.)
with the vector field Utt, x)) C (I0,T),C) ch (R.)). it follows from Lemma 2.8 of 44 that
E C (10, T1, H (R2))
for any T 2 0.
34

---

Proof of Theorem C.3 We show that the particle system (3) is well-poses and that the empirical measure is weak solution the mean field PDE. We introduce the function
N V
1 N
(x.
+1.
Since V is C hence locally Lipschitz, by Picard-Lindelöf theorem the problem (3) has a unique solution up to some time To > 0. Intuitively we only need to show that the solution does not blow up any finite time. We claim that for some constant C,
Hx(X,) & HN (Xp) .
(68)
To establish this, we first differentiate V (z:(1)) with respect to & and sum over
N
)
a,
Σ
N
1=1
VV
a
(1)
T
(1)
(0)
( VV
(o.")
VV (a)))
N
N2
i=1
1,j=1
N
<C1.0
Σ
V
+1
1=1
Note that here we have used Assumption C.I
By Grönwall's inequality, (68 holds.
Now to be rigorous, once again we define
= sup {t E R U {0o}
has a (unique) solution on [0, f)
T < 00, we define
(7)
(9) lim z
Then (13) has unique solution on [0, r]. Again by Picard-Lindelöf theorem, there exists some 0 such that (13) has unique solution on [r, T + c). which contradicts the definition of Thus, we conclude that oc, which means that there is global unique solution to (I3).
Having established the well-posedness of the finite particle system, it now follows from the definition of the characteristic flow X (t, æ, m that
x 20(0) X (t. mā
and
11 (dz) (x (1, m46))) )) H
Similar to the proof of Theorem (2)
C.2 conclude that nr is weak solution to the mean field PDE
Finally, we show Theorem
Proof of Theorem
C.3
Recall that o' 9/ 1). By assumption |и||er R and the
fact that Pp (Rd) Pv PVR.), we know that there exists C > 0 depending on R such that
||s||/y C < oc.
From the proof Theorem C.2 and Definition form
we know that the weak solutions μ4, take the
= (X(t,,, i= 1 1,
Now we bound (41,,, (2.)) using Wp 12). Let # be coupling between e1 and For 8 define (s(x)) == ([1==] + 8)2/2 be an approximation p | Given any two points
2T, D2 E Rdd we have that
35

---

0,0s ((t.2,,))) Vés (X(,, X (2,)
- x(2,
(Fem xx
/ Baa
= (X(t,x1,,1) X((,,x,
(1) Xx,,,x,,)))
W X X(t, 22, v2)) v2)) (da' dz/))
- v2)) (o X 222, 12) (dæ' dz'2)
330
(t, 22, M2 TT (dæ' da²))
R2a
1 + 12 + 13 + I1-
Below bound I, individually. First, noticing that
V cobs(0) (x))l = [||=]] + +8y0/2-1z0 s0| s |
(69)
we obtain
I S [|X(1, X(1,, 22, 12)
(70)
Next we bound Iz:
-(t,),, -x 1//Ba. (V( vs (dx')
1a <
& ( X Xt,,
(71)
(1
1/4 (F (1, (d)))))))
1/8
[x(t,a) (t,, 12) c/ /1 1/9 ||pp, - (1/4 | 11
|| | X| v1, - 22,
(72)
Here we have applied Hölder's inequality in a and Theorem C.2 Assumption C.4 and the definitions of the Po-morm and PP
The inequality is due to
36

---

Similarly for 13 we use Hölder's inequality again and get
13 S [X (1,11 - (2), || x22 12)]l
J Rai
[VV(X(;22,. (1x
x(1, ab2; 22)1 x (da' day))
- (1, -(1, ||- (1, u2)l
(L (VV( (n
1/9 (1 (1x X(t, ah, u2)l x (da/ dai))
1/8
( |x - (t, (1
2(,, ( - (dx,
1/P
(da) dx/))
(73)
Finally proceed bound 44. An application of the intermediate value theorem the difference of VV yields that
14 [t v1, X(t, 22, 12)l -I m2)1
Iew sup 0EjO.1]
12)
[x(t, - Xt, 12)1 || 2 12)l TT (dx' day2))
< (1((, - X(t, x2, 12)
|| (1, 22, v2)
1/9
(ee sup (6(0.11
12)1 x (dxx dx')) F'V (1- x22z (2)) [| 22,1111
(
1/p
m2) [(,,,,))
- X (t, 222,
(dzz dzzh)
C/ (11 -
((,, (t,
|| 22,
(
1/9
u2))l + || a/2, 22))l") x (day dab))
(fu 1x X(t, 222 21) 11 n² (dz' dzz;)
1/v
CV/ (111, + 11 || v1, xz,, (|, "Xtr
Veu (x((t, v1) x(t, 2by, 21) 1| n (dz) dzy))
1/P
S C2/4 (eit/qien + 01t//1112 1/ C21l|122 || pp X v1)) - X x2, x2) P-1
(Lre
1/p
| X(1, 2by, a1, 12) (da' dab))
(74)
Note that to get a we have applied (43 of Assumption C.4land b is implied by Theorem C.2 If we define
(
1/v
(s(
(X(s, - X(s 22)1 Ti da²))
300
combining 70
73
and 774 we obtain that
O,ps (X(t, v1) - (t,
/ (c/) gl plt #1,
((t, 12)
- (t, +
37

---

Now integrating the inequality above with respect to the coupling x (dx1, dzy) using the fact that
[ (s, (s, (x 22, 12) || (dx) dx²))
(J (1x(), X(s, 02,)) (<) dxz)) Rad
(0- 1)/r
(fenl |X(s, (2)))
1/v
Se (T))
Thus, we obtain that
0,os (X, x1, v1) - X (t, & 8C / (C1/99205)) D s Cv.Re CT (t)
Integrating t get
Φs (X(t, X((,,,,,)) S 0z(œ1 x2) + CV.Re [ (s) ds.
Finally letting 8 0 yields
(t()")(t) DB(*"(0) Cv.ne J (s)
By Grönwall's inequality, we obtain that
(t()")()) S (0) exp (Ov.Re
Now since TO [[V1, 12) and Hi,t (X - v1 1z)) 4/, v, the mapping Ed : (x1, IR2dd + (t, 2τ, v1), X(t, 22, v2)) R22d
satisties that (T)) T(41, As consequence we have that
inf / |||| xz x22 ( dx²))
xET(m),
inf (*())((t
pn)
< exp (Cv,ne CTT)
inf
(0)
exp (Cv.ne CTT) WM U2).
I Proofs for Section 3.1
Proof of Theorem
3.1 For any 8 (m,, E) O, define E(u, E) E(p). Then we have
=l p²) p²)= = (tr(Q- log det (Q- (Σ) d ( T)TTa Q- b)).
where p the density of N(,,, and p" is the density of (b. Q).
Now we consider the gradient flow on the submanifold θ
Ô - (Vo,
For clarity note that here
V D. E(O,) = (v ElH,, E,), V E(,,, Ex)
where E( Σ) denotes the standard matrix derivative (not the covariant derivative affine connection in the contexts of Riemannian geometry). We calculate that
v EnE( Σ) = Q-'(p ((u-b) -b), -)E( 1 (00-))1
38

---

Thus, the gradient flow on θ is equivalent to
(2Vx, Fiu E.) + (1 + p. V2)V Fllee 21t)
E1 - (22,, VE, Bl μb V ( 21r)
(2Vx, E(p E. w. E(H,, M1 El
ir = (1 - Q 12,) (1 + μ u4 1 (1:: - -b)
(75)
-0- Σg 2Σ, - Eg (E) + (M b) - 2E)(Ex (1:: b)p,, Σx
Note that it is trivial check that the functions the right-hand-side of [.1 are locally Lipschitz with respect to f4, and Σ, (continuously differentiable hence locally Lipschitz). By Picard-Lindelöt theorem, this ODE system given Mo RR² Eo Sym (d, IR) has a unique solution on & [0, e) for some 0. Let
: sup {t R U (oc}:
75 has (unique) solution on [0, s)}
For convenience we define the curve on Rdd x Sym + (d, IR) by
7 [0, 8) + Rdd x Sym(d, R)), 7(t) (1/: E,).
Next we consider the density flow given by
P: - SE(P))
(pat)) K(, K(. pr) dy
v
(76)
ôpt.
By Theorem C.2 we know that there is a unique solution pt in P(Rd) for E [0, oo). claim that
Claim, Pt is Gaussian density for e s).
In fact, if we let Pr ([=(()),, where @ is the immersion
P(R
0).
By uniqueness of the solution it suffices to prove that (76) holds for Pt Pt. This can be checked by direct calculation of course. But there is also a more elegant way to show it. We consider the following commutative diagram.
dôg > TRP((R))
Tge
Go
T
Vo
P(R()
Here dón is the pushforward of the immersion point 0 and gh is the inverse of the pullback map 6² TT TTR() + Tg θ restricted on Im G3 d O)). The diagram is commutative due to the fact that Gg is the canonical isomorphism the submanifold @ induced from P(Rd).
Now we show that áE E Im Gp. dog In the proof of Theorem B.3 we have shown that the maps Kp. (b, 5) to some & such that (2 = Six 41) + b, where μ is obtained from G (b, J.). Thus, Im contains all functions (or more precisely the equivalent classes of functions that differ by constant) with the form
(xx) / (x S(x 11) + bæ + C, S Ry,, b Rd, C' is any constant.
In other words, Im the contains all quadratic forms on Rdd Note that we can derive that
EE()))
log Pt + V
87:
39

---

is exactly a quadratic form. Thus,
E Lmmo, = Im tu. Ga, Im Gp. dóg.
8p:
Next since døo, maps the tangent vector
E
To,
op E Tp P(RR)) we have that
36
8F
Vo F.
6pt
Combining this with the fact that
SE
E
Im
Gpe døe, get
ops
8E
vol Vo, E).
dpt
Thus, we have
8E
G=1 8p:
G p:
vol (Ve. E) dós Go. ().
And conclude
Go. (Vo G 6E Pr 8p:
P. = doo 0. - -óoo Cdd
The claim is proven.
Now back to the original problem. Suppose co. Since we know that the mean field PDE 76 has a unique solution on [0, oo). in particular, it exists on [0, s). Note that the weak limit of Gaussian distributions is Gaussian and since Ps P(Rd) it does not degenerate. By definition of Φ we have that o (p.) E O. By letting (Mg, Again by Picard-Lindelöf theorem there exists a small neighborhood [s. + +e') such that 75) has unique solution. This together with the solution on [0. contradicts the definition of Therefore,
=s)) == Φ ( we obtain the solution of (75) 0, s]
conclude that oc. 775) has a unique global solution corresponding to the mean and covariance matrix of Pt.
Next, we prove that A converges weakly to p² as > co. We calculate the quantity A(n,, Ex). By Jacobi's formula in matrix caleulus (Theorem 8.1 in 59 we have
det M = det 2, tr Σt)
Thus, we derive that
É(H, Ex) = - ) + (H,- bb) T)TQQi it,
= tr (tQ=((1Q111))1 - 27 1222)) 2 ((Q-1 - -1 1)
(1 +1+ u4 (1) [)[Q"] 2 -b)
((1Q-1 - 1) Q- - -)) T ((Q-1-Xip , -M,
- -m4 - b) TQ Q02((1 - b) S.
Noticing that
0 -5-1fo /-' Blu E.) ds = Elpg Ex) E(u Ex) < oc,
we obtain that E( Ex) - 0 as + co, which is equivalent to hr - and Ex - Q by checking the expression above. Thus, we have shown that Pt converges weakly p" which is the density function of (b, Q).
Finally, we show the convergence rates of ne and Er Since we have already proven that Pt p², it implies that
He - b= o(1), Σ. Q o(1), Er¹ Q- o(1).
40

---

If we set 77, 7 b) and St = (Q-1 - Σ EE, then
-E(ME: Ent
tr ((Q)' - Ex 1 10-00'110 '11, 1111-b)) b))ut) T. (( + Q QQ'(0m b)y,))
+ + b) -M, (14-b)
tr ((S: + n,1)) T(S, + +1144 )) + :i Th
Ie
H4 2 la (1 (1+1 14 Pa
vv)((S))
(St) m []
m 8 la
b, @ la ba)
vec(S()
7i J]
+ + + I|md)))
be &
la
7h
On the other hand, E( Ex) can be written as
1
E( Σ,)
(tr(Q- E)
log (et((Q + (1 b) T) b))
2
1
(Ir(S.) log ded + S:) + ni Qn:)
2
tr(s, S.) + T Qn. + ((1(s.1112)
4
Le
vec(S:)
((s) nc' 7|
+ +((1
2Q]
7
Now we prove that Ve > 0 there exists T > 0 such that Flpp Ex) 2 4(7 F) E for 82 T.
It suffices to show that
Is
b, & la (1+b, Ta
2 1 201
b. 1s
2Q]
which is equivalent to
le
1 b 1/2
v2
b
Q-1
b1+b b (1+bb))Q
V2
This is true because by definition is the smallest eigenvalue of the matrix.
By Grönwall's inequality, we know El4, Σx) = ((ee0((eee-))))00 4(y-e)t). Thus, we conclude
l|+. bll O(e -()-())))) [, QI O(e (0(--(()-))) Ve > 0.
Finally, we provide a lower bound on Note that for any >> 0 b * 0 we have
1
La
b b
11/2
+ u
V2
> 0.
1 +
b
Q-1/2
b bQ
2
Thus,
L
1 2 1-2
> 1ta I
& 1-2
1/2
{1 + b b)Q
v2
(1- ub b)Q-1 (T-vubb))--1
= Sm.
Since Amax is the largest eigenvalue of Q, we know the smallest eigenvalue of S2 given by
ub Tb
min
where u > 0.
1 + u'
2Amax
41

---

We find u such that this quantity is maximized and get
1 ububl b
2
7 2 mix
4>0
1+ 2Amax
b.u b + 2Amax + bbb 21 )² 8Amax
1 b + 2Amax
If b = 0. then the smallest eigenvalue is given by
7 min
2Amax + 2Amax
Proof of Theorem 3.2 Equation is a direct corollary of Theorem 3.1
By Theorem there is unique global solution. Thus, we only need to check that the E, given by (6) and (7) satisfy the algebraic Riccati equation (5)
For
E 220 we take the derivative with respect to and get
+ (1 (1-6 20-
(77)
bー E
5-1
2e 21 (Q 1 1
(78)
Substituting 78 into 77
we get
2Σ, = 2Q- -2= Σ
Multiplying by 2g and using the fact that E, and Q commute, we can see that the algebraic Riccati equation (5) holds.
For
E, = - 17(1 E 2t
1 ne 21 U0
we apply the Sherman-Morrisn formula:
01-a 21)
711 2-2
[E, (E,)_1
1tne 21
vv =I
vv
+
m1 20)
1 +
1 1tge
Thus,
1 = 21 11(1 e 20
1 UV )
ne ne-21 2t
)(+
17 - c c-21
2 -Q-1
vo
-2 00
1+
+ 1Je
2ne 21
10 T 3 -3,(555'' (5,
7
Moreover,
770 -21
1 77 20
1Je
v0
00
+ mpe 20
+ re-
Thus,
2ne
+ 2 tr(vo" T) = n7c n--22 2t
2r 2tr( Q-))Ee)
tr(v v)
+ 1e--2
1 176 22
On the other hand,
det(Li,) = 1 17(1 e 22)
1+
v
1 -
+ 72-
Thus,
-2ne" 2
2ne -2
O log det((2)) = + me
ne-
1 + me -26
42

---

Therefore,
Pt (2x)
112 4/2 (2)))//((e((())) -1/2 exp
② _ x
1
with
D 1 + 1111 - 26
20 1 7/6 22
00
a solution with the initial condition En = Ia. The theorem follows by the uniqueness of the solution of the mean field PDE (Theorem C.2)
Proof of Theorem 3.4 Similar to the proof of Theorem 3.1 we have
14 =
E( Σx) = Q (14: - b) -V -
and
= -) ",,1111(""""3"3"""""000000""""""
+E, = 2((1 - D>Σ, + vz) "1=1 (11 vI) E?
Following the arguments similar to the proof of Theorem 3.1 we can show
has unique global solution,
Pt is the density of dfN(H, Σx) given by 8
((0)) s 0 and E(u, Ex) - 0 as + 00,
Pt converges weakly to p²
Finally suppose ΣgQ = Q.o. then E, also commutes with Q since 0 is a solution the ODE satisfied by Σ QE, and the solution is unique. Thus, diagonalize E, and Q simutaneously. Then there exists orthogonal matrix P such that
Σ PT diag (o
,0
(1) }.P. Q= pT diag P
And 9 reduces to
20 (1) (X: -0(1)
(1)
d (1 - v) + 11
Solving this ODE, we get
(a) 1,)
(,:0) - A,) (1- 21
(o, )0)
(0:0))-
Thus, have o0 x as and
0f Au = =
In particular, conclude | - -Q11 O (e -2/((12///+)) where is the largest eigenvalue of
Q.
Proof of Theorem 3.5
First, we define the Hamiltonian on the centered Gaussian submanifold by
H(E, (5. S.)) + E(E)).
By Corollary B.5 we have
G s = 2(E25 - SY²).
Thus, the Hamiltonian is reduced to
1
H(E, S.) 2r (
(tr(Q- log det((Q-E,)) d). 2
43

---

Therefore, we have
+ ST
Ig)) Vs,H((E, 2 2 + S,
and thus the Hamiltonian AIG flow the Gaussian submanifold is given by 12
Next we show E, is well-defined and remains positive definite. We check that He == H(E,,, St) is decreasing with respect to L.
dHd dt
= tr S. + Vz,
= Ir (vs, He (-a+S; Vs, H4) + V., Ho H.
= -4, tr 52) 5.
Let 03 be the smallest eigenvalue of Σg. Then
log det (0,0 = log det E, - log det Q 2 dlog ar log del Q.
Therefore, have
llg o. + 1) 2
log del Q s - (log del deeddel())) (E,Q-1)
<E(2:) S H4 < Ho.
which yields that
2 exp (log det Q/d - 2H d 1). This means that the smallest eigenvalue of E, has positive lower bound. Thus, E, Sym+ R) for any & 2 0.
Finally we show that the AIG flow on the centered Gaussian submanifold coincides with the one on the density manifold. Similar to the proof of Theorem 3.1 we have a commutative graph
døs
T T(P()))
G2
T
P(R()
Here ths S - f(x) TT + C, i.e., maps symmetrie matrix to the quadratic function ((x)) (or more precisely the equivalent classes of quadratic functions that differ by a constant). Now the only things we need show are that AH Im thy and that &H E Im dØs. 80,
opt
874 -G E Im døy is trivially true from the commutative graph. Now since SE 50:
log P V E
op:
Im th (Pr is centered Gaussian density and check the definition of thy), it suffices to prove
8
&G de dæ E Im 1m.
8p:
Note that we have
5p: f
8
Φ.V (p(t) / K( ((-B))())VV
N.G
d.
8p:
5p. / Vd. (ml.)) / ((
V / (((:))(())V
(Stz)) ((SE)
2 (SP + EnS?) x E Im
Thus, the proof is complete.
44

---

We remark that the fact the Stein AIG flow remains Gaussian is highly non-trivial. In fact it requires 5% to lie the cotangent space of the Gaussian submanifold. One sufficient condition is that the öp.
variational derivatives of both the kinetic and potential energies lie in this same space, and the former could be interpreted as: The Gaussian submanifold is totally geodesic under the given metric, meaning that any geodesic flow with an initial position and velocity chosen from the Gaussian submanifold remains Gaussian. Fortunately both the Wasserstein metric and the Stein metric satisfy this property.
J Proofs for Section 3.2
Proof of Theorem 3.6 From the proof of Theorem 3.1 we know that (16) has unique solution that
continuous in t and bounded for any t E [0, T" with T co. Thus, the linear system (15) also has unique solution.
Now check that the sample mean 14 and covariance matrix C. satisfy (16). We simplify the right-han-side of (3) using 14 and C₁-
N
RHS- = 2;
(4)
(a
T2.)) (0) + +1) (x,)" b) (t)
N
j=1
(t)
Q
(a)) (1) (b) (1)) xz +1)
(0)
j=1
2,
(
Q N
(1) ((x) ( (8) +1) + Q 2 1
Σ
((a))) xx +1)
(1)
(1)
N
j=1
(1 Q (Ce + mmphi + Q Pig )x, (0) Q pe + 7 b.
(79)
Let X, (m) uN ) Then we have that
XJ1
X,T X5 N mel. T.
14
C.
N
Then 79 can be written in the matrix form as
-14, X, xxx(1 (C2 + ) + n.b (011111 1 QTTQQ11 1-11-0000 16 16 Q-1
(80)
Multiplying by 1 //N on the left, we get
ir = u.. - u. (C2 + p. + ut u., 1o- ----Q----0 1-p- -)) -1
Thus, we have
ir QI- (1 + MT HTTHT.H ( b).
Note that Õ = =x, X, + X1 X.)NN uiù iu Substituting (80) we obtain
Ö =C, C. (C, + u4 (1 1 b) -Q- (C, + (M+ b) C..
Next show that 14 and 15
satisfies
13
RHS (1 Q- (C, + mm ) + Q-
)a (1) Q-1
u1 + 7Q-66
-I-Q- + m.H.) + 1- =b,,)
(t)
p1 +
(I-Q 10)) - (1 +1+11) m² 1.) - b)
-I-Q- + m.p..) + - " (0)
(b) + i..
Ax (0) -1)) i4, = (1())
By Theorem C.3
13
has a unique solution. Thus, the solution of 13 is given by
16
The proof of Theorem 3.7 is quite long and tedious defer it to Appendix D
45

---

Proof of Theorem 3.8 If Co is non-singular, this is a direct corollary of (heorem3.66amm Theorem 3.2
To also accommodate for the singular case, we provide a direct proof by calculation. Since Co QCo, we know that Co and Q are simultaneously diagonalizable. There exists some orthogonal matrix Po such that we have the spectral decompositions
Co = =P Doto = Qofo.
where Do = diag(A)
(0)
) and Qo = diag(g (11, ... qu) are diagonal matrices. Let D, :=
PoC,Po (18) can be rewritten as
a = Py diag(( ((e 2 + (1 - =))
) (0) /q) 1/2 , (e- 21 (1 - E 21 )
(0) /qa) p Pox!
Thus, by taking the derivative with respect to t, we obtain
1 1 10) /n1
1 - /// /
(100 1110= c diag
Poz! (0)
+ (c-21
- (1) (0) //1 (3)33/22
(c-2t 21 + (1 - 2) (1) /qa)3/2
-U(c 21 (1-- 20-11--111 -3/2 (0)
x;
where
U U=e-2p
2t PJ diag( (1 - (0)
/u1, ".
(0) / Po = 2I-
On the other hand, we check that
1 N
Σ (cx x (1) +1)00=af()
N
= =(e 21 1 + (1 - -20)0-116) 1/2 (0)
Σ
(aa)
(0) ) ( + - (1-- -- -1 (0)
"0²f" +1)-
x)
N
j=1
(- 21 (1 e "-QQ-11
-1/2
(0) xj
N
=(e 2 + (1 - 2010 (o) -1/2 (0)
Q
(e 2 + (1 - -------- -1/2 (0)
N j=1
1 N
Σ ) T (e 21 + (1 - c 21 1 C))
-1 (0)
x
x
N j=1
Q ( 1 (1 E -2 Co)
1/2 (0)
x;
=(e 2 + (1 - e 21 10-2Q-11100 Q1) -12 (0)
x
1 N
Σ Q- (e 22 + (1 E 21 -0- -1/2
(c-2t --=00)1)00000 (x)))
(1 - e 20 --)))--00) x ())
Ñ ==1
=(e 21 I +
2010 Co) -122 (0)
x;
( 21 + (1 20 CC) 1/2 Co 2 + (1 - 20- (c)) -1a,"0))
-1 (0)
=c 20 Q (e 21 + (1 - c 1 Co)
-3/2 (0) x;
21 (1 - Q Cp)
-3/2 (0)
Thus, we conclude that [88 is solution of By Theorem and hence the theorem follows.
13
C.3
the solution for
13 is unique
Proof of Theorem
3.9
For any particle system of the R-SVGF, we can derive that
X = (1 (1- v)
X,XT
+
11T
+ +VTa
Xx N (X.x, + +11
N
N
X,XT
T1T
-1
-11- v)
+vTd
+ vId
X₂(I -.Q
N
N
46

---

By Sherman- Morrison formula we have
11
-1
1
1 v 11
VIa (1-
Ia
N
V
N
By Woodbury matrix identity we derive
X
11
(1 v)
+
+ tuld
N
N
11
-1
11 T
-1
vla + (1 v)
vlaa + (1-v)
Xx:
N
N
11
11 -1
N 11 Ia 1 Xx
VIa +Ia+ (1-
X₁
Xi
(ula
(1la+(1-)) (1-1)
N
N
la 1 1-1 -
11
1
-1
X1
la + V
C
Xi.
N
Nu2
1
1
1
Substituting this into we have
X. H la 1- a 11
X1(1
-GQ-)
-2222 X
1
1 +o +/)))a)
x, x(x1(1--QQ7))
N
Nuē
1 X1 - 1, -
1
Ia+ 1
(,a-- -0-0))--
X1
Ci)
(
v
1-
1 X1 Ia - ( (1 Ia+ C.
1
-1
1
Ct
(I (I-CQ-))
-xT + -1-
(81)
Multiplying by x, on the left we get
x X
(vla + (1 (1))))-----))
Thus, the derivative of covariance matrix C. is given by
G = X, Xt X,
N +
N
=(vTa + (1- (101:000000000110000000000000000
=2(vIa (1-1),,,)- (vIa + +11-0011-00000-0000000+------1-1
Next we show that 19 and (0) satisfies 81
RHS A (wla + (1- XoAT
Similar to the proof of Theorem 3.4. could be shown that the R-SVGF also has a unique solution and the proof is complete.
Next we show Theorem 3.10
Lemma J.1. The covariance matrix C2 (t = 1,2,.... of the discrete-time finite particle system satisfies the following equation
C++1 ((+(((---"(()((((((((-000))))
Proof. Let X = =x11 ,NN) T Then (1 can be written as
X = x. -x
Thus, we have
x, X, C+11= N
= - -(+(([---ca)(((((((--000
47

---

46
y= f(x)
1+
Figure 6: Plot of the function fe(x) (1 c(1 x)))
Proof of Theorem 3.10 First since CoQ = QCo they simutancously diagonalizable. By Theo- rem 3.8 any Ct and Q are simutancously diagonalizable. Thus, without loss of generality we assume all Ct and Q are diagonal matrices. Then by Lemma J.I we know
Q- = (T + <( Q C-))
If we define fe(x) (1 + c(1 - =))²x then every entry (i.e., eigenvalue) of Q-1C, follows the fe fe-iteration trajectory of the corresponding entry of Q -Co.
The fixed points of fx) (1 + e(1 z))²x are {0. 1, 1/2 + 13. IF 1 then (f)( > 1. (/'(2/e + 1)] > 1, and ['(1)] <1. By Proposition 1.9 of [26]. 1 is an attracting fixed point while O and 2/c 1 are repelling fixed points. By definition there exists interval around such that for all initial points in that interval, the trajectory of any eigenvalue of Q- C2 converges to 1. We now quantify that result further.
Note that f((x) = (1 + - (er + E - 3Ex) is an upward-opening parabola whose zeros are 1/3 + 1/(3e) and 1 + Thus, f' is monotone decreasing on [0, 1/3 + 1/(3e)]. In particular, both equation (z) = 1 and (x) = 1-- e have two distinct roots, the smaller ones lying on (0. 1/3 1/(3e)). Define w. to be the smaller root of fz(z) 1 and ue, the smaller root of f((x) = 1-16. e. Since f. is monotone desereasing and f((1) = 2e, we have < w. < Ue < 1.
Thus, for any x E [w,, 1/3--1 1/(3c] we have 0 S fz < 1 and 0 S fe(x) S 11/3 + 1 1/(3)) < 1/3 + 1/(3)) (here we have used the condition c < 0.5). On the other hand, since (0) 2 1 and 0 and 1 are two fixed points, it holds that fc(r) 2 x for any E [0, 1]. which implies that for any x E [w:, 1/3 + 1/(3)) we have fe(x) 2 fe (w.) 2 w.. Hence we know that f is a contraction map and fe ([w,, 1/3 + 1/(3))] C [w,, 1/3 + 1/e)]]. which implies that there is a unique fixed point (i.e.. 1) such that the f--iteration trajectory converges to it.
Next we prove that for any æ E (0, 1+//e) the trajectory falls into the interval [w,, finite iterations. I E (1/3 + 1/(3e), 1 1/e) then after one iteration we get fz(z) E (0. (1/3 log w./r
after
1/(3c))) (0, 1/(3))) We claim that if (0, w:)) then after to = 2og + e(1 - w.))
we have f! (fo) (x) E [w,, 1). Firstly (x) <
for any t. It suffices to prove f ((t)) (x) 2 we Suppose f( (x) we Then for any 0 S S to we have f!() (xx < we By definition of to we know f((()()) > (1 + <(1 - w. 2 2 we.
This is contradiction. Thus, the claim holds and in conclusion for any x (0, 1 + 1/c) the fe-iteration trajectory converges to the fixed point 1.
Finally we consider the case when E [u,, 1/3 1/(3)). Since is monotone descreasing here have f((x) E [0, e-e. And by similar argument we know fe( 1/3+1/(3e)]) c [ue, 1/3+1/(3e)
Thus, we conclude
| (x) 1| s (1
- 1 s ... s - 1) -- - 1 S e e-c 1/-
48

---

K Proofs for Section 4
The following result was derived in 43 We provide the proof below for completeness.
Lemma K.1. Let p be probability measure and ps be a Gaussian measure with parameters
0 = (m, E). Then we have the following expressions:
1
Vi KL(po p²) = Le (V())),
Vy KLp( (1) (Bexrip [=][ (82)
Proof. We compute that
Vp KL ll p²) Vu / log po(xc) p² (x)
=V==va
po(x) dæ
V (upo(x))
po(az) dæ + / log Pe(x))
V (PPo(z) dæ
po(x)
p (xx
Ve I po po(x)) dæ - log Po V (x) dz
p- (x)
Ve log Po(ar) p=((x) po (2x) dar
/ V log (xx) (x) da
Vapo dx
E Pu (V(x))
where we have used the fact that pa is Gaussian density in Similarly, we have
pa(x) log
Vx. KL(po 11 r²)
V3
po( dæ
10 (xx
V VT,Po())
f log po(x)) (x dx (x(x)
Pal (x( dx
(x
log pa()
1 - (2 u))x 11) 1)1=-1
Σ-1
po po(x)) dx
p*(x)
2
2
1
log po(xr)
22 (2p0(1)
da
dæ
2
pa (x)
Pa (#(
2 log
po(2) da
P" x
< (-))
Here again in (a) and (b) we have used the closed-form expression of Gaussian densities. Thus,
1
V VuE(,, 2) Ex-po (v)),)
VzE(H, (Ba-no (Ba~no.V"V(x)) -
Proof of Theorem
4.1 By definition of Gaussian approximate gradient descent, we have that
0, - 6. (V,, E(O:)), where E(O.) KI[(p9, || p²) and (14, E).
Applying Theorem B.3 we obtain that 0, = - G (Vo, F((+)) given by
ju (2V (2V,, E(u, E,)E.4 + (1 + u1 V E(m,, Ex)
Σ (2E, Vs, F( Σx) 14, V E(H,, E.)) (2V Elu,, E.EE V F(P, E.),.) Ex
(83)
49

---

Substituting
82
into 83
we get
i4 = Ex-pe
V2 (x: Ex) 14, (1 14 14, Ba-pi
1 [V(x))
Σ, = 25, ES EEpEgopp (=x)) M.Bu VV
(Exopp (Vi((z)) Earp, A (zx)), Eg
Next we prove the convergence.
E(H, Σx) = tr(V E(H, I,) E.) V E(H,, I it, t1 ΣJE E, (2V:: E,, V.LBu V B1, 2) 5.
Noticing that
0 25-11 1/ E(u, E.) ds E(Ho. Eo) E(H, E,) < oc,
obtain that E( E2g - as 1 oc, which means that there exists Mo E such that Pe converges to Poo. Poo is the density of N(H Exe). (Since E',,, Ex) is given by the KL divergence between Pr. and p², by Lemma K.I will diverge te or 2, diverges.) In particular, it satisfies that
{ 22 [2Vx:B(u Exx) + V E( E.o ) =0
V EmMoo Ex0)= -0
which is equivalent
VI Eoo) = V Eoo) 0
and implies that Poo pa
Lemma K.2. Given otroongyy convex measure p", we define 8 as the unique minimizer of KL(pp p²) where pe denotes Gaussian measure with parameters 0. Then it holds that
E| (VV(2)] + tr ((E,a [<)] - (V))
2a (KI(Po ll p²) - KL(po- || p*))
This is proven in the Appendix D of [43] The proof idea is to consider the Gaussian approximate Wasserstein gradient flow from Pe with the target p²
Lemma K.3, Given stronngy convex measure p", we define 8" as the unique minimizer of (L p²) where pe denotes a Gaussian measure with parameters 0. Then the Wasserstein-2 distance between pa and pą- satisfies that
aW2 po.) I KL(pp ll p') KL p').
This is Lemma E.2 of
Proof of Theorem 4.2 From Lemma K.1 and the proof of Theorem 4.1 we know that
E(pe Ex) = tr ((2Vs Σ)Σ + E(,, Σ) TTV EE EnE(,,
V.JE(H, V
--rr ((Emol V2 + Epg [V(x))] (F,p VVV(2) 1+ Ep [V()))
- -| [V(=))]
For convenience we let 77 Epp [v())) and S, (Epu [2V(x)) - Σ-1))1/2 Then get
-E((4, Ex) = [vec (S.)) ni]
La X Σn
p. 2 1!/
vve((S))
me
(1 pº 1.)
7:
(ve((St) ni
=: [voeT (S.)
nol M,
50

---

Noting that M, + Mo. for any > 0 there exists T 0 such that the smallest eigenvalue /// of M, satisfies 2/ 2 7 e for any t T. where a/a is the smallest eigenvalue of M..
Moreover, by Lemma K.2 we have
vec(S:)
Ivec (S,)
mi]
2 2a<(KL(po Il p²) - KI,(po |l p²)).
7
Thus, for > the derivative of KL. divergence is controlled, i.e.,
(L p²) 2 2(7 - (KL((p |l p²) KL ll p²))
By Grönwall's inequality we know
KL p²) KI p²) =O(e 0(E-2(9--)6)
By Lemma
K.3 this implies
(p(0,, 00e 2((-0)())))).
Noting that
(p pa.) 11 14 11 + tr (221/2 - -127)
we conclude
114 14,-4"11 ==11==100- O ()--())))) | - Σ O(e O Ve x 0.
Finally, we provide lower bound on y. Note that for any
0 if μ * have
1
1 O o (E7)122
1+u
(1 + u)p Ia Ju @ (E
Thus,
La 1 D 1 T & (E
"¹ @ 122)
u Ia la X E
1tu
= Q2-
(1+ T "²)
(1 upe T
Since E" satisfies that
Epg
and (x) < Bld, we know the smallest eigenvalue of Σn is at least 1/8. Thus, the smallest eigenvalue of Ωr
B(1 + u) wr] where u 0, = TT n²"
We find such that this quantity is maximized and get
Y min
2
1
u>0
B(1 + u)
1 + ((1+r) + / (1 3(1 3(1 -4
3(1 + +1
If 4 = 0, we still have
? = -min ( {8/1/> > a
8 +1
L Proofs of Uniform in Time Convergence
To show Theorem 3.7 we need lemma on the convergence of empirical measures the i.i.d. setting. There are results for general measures on this given by 25 47 24 but for our purpose we always have a Gaussian distributions as the limit and there are tight results with faster convergence rates as shown in 45 46):
51

---

Lemma L.1 (Convergence of empirical measures for Gaussian distributions). Fix the dimension d 2 1. There exists constant Ca such that for all N Z I, with UN = N Ex= Ox, where {xx} is i..d. sequence drawn from μ ~ (O,, we have
N-1 log log N - N)² 2-2/d if 2
if 1 if
W (AN S Co x
Proof of Theorem 3.7 Suppose the sample mean and covariance time t is m. and C₁, and that the mean and covariance of the mean-field limit is fe and E,. We bound E[ (CN.P))) :c)] using Theorems 3.1 and and Lemma
L.] in six steps.
Step I. We prove that (z! has the same distribution (3/ where (5()) //2/2 CO 1/2 (z (0)
mo) + my, Note that according to Theorem 3.6 where we have a!!" Ar(æ; (0) mo) + m. Here A, is the unique (matrix) solution of the linear system
A = (T- (C, + (.H) + Q bp,,)) Ar. Ao =
(84)
and m, and C₁ are the unique solution of the ODE system
m = (1 --)) - -11 + m m m.)Q bb)
(85)
C = -C-20,-00 C (C, m.(m: - bb) Q - (Cr + (m: b)m,)c
Since C2 is the sample covariance, have
A.C.A C2.
which simplies that
1/2 (C, A 1. (C- A.C
Thus, P. = CT an orthogonal matrix, and we have
A. - C//2P.C5''''
Since the multivariate Gaussian distribution No Ia) is invariant under orthogonal transformation, the joint distribution of 1/2 (x2 -mo) 1=1 is also invariant under Pz. Thus, we have (ax m,))
N
has the same distribution as a (1) m.) and Step is proven.
Step II. Establish uniform decay rates for ||C, Q11 and Imm bl1. We begin by checking the energy function
0 < 1(m.:mm,,. C.) = (tx(Q- log det CC;) d (m, T) b)).
As shown in the proof of Theorem 3.1 we have
E(m, Ct) = --llc----- 1 m.(m, b) T)- -1l---1- -b)1 - 0.
Thus, E(m., C.) & E(mo; Co) for any 2 0. Furthermore, similar to the proof of Theorem 3.1 check that
É(mr, C.)
=[vecT (Q--((Q-),1) (m; -) b)TQ-11
Le
m. E la
vec (Q-+c,1 - !)
m.
la
(1 +1+ m. m
Q Qmmm b)
1
vec (Q-1C 1)
271 [vec" (Q-10- 1) (m, - b)TQ-1)
2Q
Q Q-I(m b)
(m:
tr (QQ-'C, - 1) T (Q C, - 1)) + (m;
b)
Q(m;: - -b)
22n (tr(Q- 1) log de (C, + (m: T) 1Q(m, b))
E4y, Ct).
52

---

where y is the smallest eigenvalue of
2 v2 m, @ Q-1/2
1 v2
m,
Q-1122
m, mi m
and as shown in the proof of Theorem 3.1 it has lower bound
m, mr + Qmax
where Imax is the largest eigenvalue of Q. Now since
1 (m - -b) b) (m - b) S E(m,, Ct) S E(mo, Co);
2
we know (m4- - 6)_TT (m - b) S 2
Emmo, Cu). Thus, [mm
upper bounded by some quantity
Fi mm. Hence "/ uniformly lower bounded by
1
Y inf Y 120
1+ 1+ F? + qmax
Thus, by Grönwall's inequality we have Elm,, Ct) S e- E-t'''E(mo, Co). There exists F2 F2,Q.b,Co, such that |m, bll s e F22
Now similarly
1 (tr(Q- 1) log E(m,, Ci) S " Co) 0 S 2
also renders an upper bound for ||C, (CC--Q1F with exponential decay noting that tr(A) log det +A) is quadratic in ||A|[ when ||A|]] small, i.e., there exists F& C. such that
||C - F S ---- (F.
27
By Theorem 3.1 we know that p4, E satisfy the same ODEs as m, C₁:
c4 (1 1--Q-12 Q-12),)) (1 ue 1
(p+ -b) b)
(86)
Σ = 2Σ, - Eg (Ee + M, - b) - 1-1---11-1 "Ie + (14: - bb) E.
Thus, similarly we have
|1 - bll S
c 22
[2 - 111 S e -2-
Step III. Show that ||C, Ex f and ||m, - 1el can controlled after sufficient time.
For any > we define
+ {(v, S) Id Sym (d, ||0 bll S e, and (1 =
and consider the relative energy function
1 (tr (c. - 1) - log -.) + (m: -.) 2m -.))
E(m;, Cty Me E.) 2
show that for small enough if (me, Ct), (149 Σx) E 0,, then
E(mi, Ct, p4: Σ,) < 0.
We derive that
i E(m, Ct, 14, Ex)
1
CT (E, Ci) (z - 1.)): + (m, TQQQ'mm - (.))
tr 2
53

---

Here note that
CF - Gi) (Σ, - - C,
CT ( 1el14 - -b) -C2-mi -m)))
+ + + Q- + (H,, 6)μl )Ee -
Q Q + (m: b)m,
2(G (m+ -((14-b)]] (0))Q=a 1-²
+ - Ex)
b)
20, -=)²0-1-
- -(m,(m; -b) -m((1 - (
2(C, εx)²Q-
(He(me :.) (m; - H))(m: )) ))7)²
On the other hand, we have
(m, - Q im))
-(m, - -Q- - 11) + (m: - QQ (Cmm
- -E(pe)
- -Im, m 14 T (m) 11) Q
2b+ +m -4)) 102(mmmm)))))) my - p p.)
mi
=:1 + 12 13 + 14.
Here we have
1 = (m: - TQ -m) (m, - p,) (- -
- mi)
2 (1 <m =T)TT
QQQQQ (m: 14) + (m, T Q- Ex)u
2 (1 c) m1) T Q (m++ - M.) 2
tr (1C( - 20)²Q-''771
(I, - m.) T)Q0-9)))))
2 (1 cm m1) TQ (--@0mmmp) (4)) tr ((C, - -1----------1
(1 + c) (1m m1) )))-))
and
Ig 14 (m, -1,) Q== T (H (m+ (1)) (mim (m )m: - bi).
Combining all these together we have
F(m,, Ct, M4, E+)
2 2m
(0)- m.) T4))(QQ= -22--()))mm -M:) - + e) tr ((H (1 TM) 2
1
2(1 -c) ir ((H((m: - MO)T + (m: m)) - -b) Q
1
+ (m: - m1) TQ (Miri (m: -11) + (m.m - unn ) - b))
c(2-c) d2-c) Q (mL m.) + 0(e))
-(mi-ui) QQ- -0---0Q-1 2(1 e)
Since 0 I I 3 Q 2 I, and mill and || are bounded by F1 shown in Step II, there exists 6o = mo. (Cg can be seen as a continuous function) such that for any c S co we have É(m,, Cr, M4: Ex) < 0 as long have (m, Ct) O. and (H, Ex) θ..
54

---

Given Q, b. Co. my, Eo. Mo suppose at time to we have (mi,. Cto) E eeo and (H1 Etoo E Dea Then for any >> to we know
mm,, Ct, Cr Ex) S E( Ctas Mep Ene)
S tr((Ct,
Etg CTT2) (mos 110
1QQTTQ--000 11o0)
1
$ 4(1 co) quiin
11Cro Itoll
|
14
qmin
Note that
tr(A) - log deg (1 + A)
lim M4
tr(AT A)
2
Fixing any 6 > 0 as long En is small enough we have
E(mi, Ct. M4: Σt)
tr((C; - -E))2 (m+ (1.) ( (m:
1m= 14)
4+8
1
(4+ 8) 2 11C, --el +
m+ 14,11
quax
Therefore, we conclude that for any given Q, b, Co. mo. Co and Mo there exists Eo as stated above and F4 Fa,q Co.. 0 such that long (mi:: Ctu) Oe, and (Mrp Eng) E Oc then for any t > to
|[ E11i s F4 (1(Cro -ir + || M4t 11)
m, 1ill Fa (11Cto
Eeg 112- + 1|mm 11)
Step IV. Uniformly bound ||C, Ee and ||m: - 1111 using ||0 - Eoll? + |mo - -oll..
Note that by definition of the Frobenius norm (or any matrix norm). given cn 0 there exists
E (0, cn) such that for any S E Sym(d, R) as long as the norm is small enough, i.e., |S
Ql( S 61, then have (1 <q = < (1 + en)Q. Then by Step II we know that if we
F5 = max F2, F3, F2, F3)) and to > 27 log Fo then the following bounds hold:
m4 bl <1. |[. - Q1 < c1
| - bll c1,
E - Q11 <..
Now it straight forward to check from 885 and (8) and the results in Step that there exists Fs = FoQ mo such that for any 0
d de
lm, 1,11 Follm, m. Hell,
[IC, E, s FollC, E.1
Thus, by Grönwall's inequality we have
Im 1211 < e Fot [|mo - moll? ||C, - E2 - Po - Eoll
Combining this with Step III, know for any 2 0, there exists F- Foto F4 (only depending Q, b, Co, my Eo, (b) such that
|[C, Eglli S Fr (I[Co Eoll + mmo 10
m: - 11 < Fr (I(Co Eo13 + m 100l3)
Step V. Let y: (0)
2 (x (0)
(a)
mo) + pr
Define 500
N EN
6. cn
and 7(0)
<N
N Dm 8 000 We show that
E
W (112
S
UN
(4)
log log N
55

---

Since o is trivial by the definition of the Wasserstein metric, we only need to check c. In fact,
N
1
N
N
3
Σ
(z) 2
Io 1/2 c/.32 (z. - Ho) |
N
/ 0 1/2 (mo - Hb)1| + +3 - mil?
+3
=3 | 2 1/2 c/ c./3112 +3 c Co 1/2 (mo- moll" +311 -mell
Is + 16 17.
Note that
1. S 3|1[cll 1C || mo - Holl"
I7 S 3F, (11C- Epll | + Imm M4 1)ol3)
By the central limit theorem VN 14) converges to (O, Eo) in distribution. Thus, N log No N (mo- Mo) converges to 0 in distribution and hence also in probability. (There could even be almost sure results using the law of iterated logarithm but converge in probability good enough.)
Similarly by CLT every entries of VN(Co 2i)) converges to a Gaussian distribution. Thus,
V log log (Co- Eg) converges zero matrix in probability. Therefore, have N Eg I. + in probability.
N
log log N ||C
+ F3. All these constants here (F3, F7, ||C (|e) can be or chosen as a continuous function of Q, b, Co, my. mo Mn and by continuous mapping theorem they converge to the values of the same function with Co Eo and mo Mo Thus, we conclude that N log N (4 + 17) + 0 in probability.
By Step II, have & S
Now we derive that
-1// 10- + 1/2 - /// 1|2
Is $
F <2112//- 11 || + || -
Now we show lemma: Suppose A, B E Sym (d, I) are two positive definite matrices. Then we have ||////2 - B1/2|| < 2/A I|A B|| where À is the smallest eigenvalue of A and B. Note that we are using the spectral norm here.
In fact, denote the largest eigenvector of A B1/2 by n, and let æ E IR be the corresponding eigenvector such that x = We have
||4 B|| (A B)x
x A / 12 + x T (A B222
= T (A1/2 + 11/2)a // 2 21/\A.
Thus, we have that ||A//2 13 S + |14 Bll. Moreover, the Frobenius norm is bounded by
Vd
Vd
1A 11/22BB/1115 S 4//2 11 5
A -BII
S IA Bl
2v
2V
Applying this lemma, we know
Vd
|
1-1/2
Dn CollE,
2 X
56

---

where À is the smallest eigenvalue of Co and Eo, which converges to the smallest eigenvalue of Σo as N goes infinity.
Next we need to show that the smallest eigenvalue of C. and En are uniformly (in time) lower bounded by some X > O (depending on Q, b, Co, ma, Eo, Mo). We revisit Step II, where we show that Emm, C+)= & mmm,, C.). Then
Lr(Q-'C,
log (e C,) S 2E(mo, Co)
leads to uniform lower bound of the smallest eigenvalue of C, since tr(Q tr(Q-'C, 'c, log dct(Q--(.)) + oc as the smallest eigenvalue of Ct goes down to zero. Thus, we have
1)
C - /
d IE, Cells 4x
S 4X'
dFr (1ICo
Eoll? + I m)ll3))
Following similar arguments conclude that in probability
N
la + 0 in probability N + Therefore,
log log
N
N
1
-10
N log log N
(1s + Lo + 17) 0,
log log IV
I
which implies that
E
W
SN
(20) -
log log
N
Step VI. Apply Lemma L.1 get the final result. Note that 2!! are i.i.d. from (O, and y, a linear function of at; (unlike m, and C. which are random, ,4, and E, are deterministic). By Lemma L.I and Step II (Mt, and Σ, are uniformly bounded) have
N- log log N N - N)² 2/d
if d if d = 2 if 3
E
W
SN pt)]
CQ.b. x
(87)
Thus, we derive that
WE(CH pe) (2) E W2(CH pe)
(*+) E
(MaLGH .500 + Wa(c( p.))
ww (5N) (70770 + 2E W3
EN < pr)
(1)
2E
N-1 log log N
d =1 if d =
(w(ss) S
CQ x
N (log N)² 2/d
if 2d23
Note that we have used Step I in (*), the triangle inequality in (*+), and Step V along with 87 in (* **).
Proof of Theorem F.4 The proof is roughly similar to that of Theorem 3.7 But Step III is little different (can be strengthened and simplified at the same time).
We show the global contraction of ||m, -1ell and also the contraction of ||C, Eell after sufficient time. Note here ||= the Frobenius norm.
Frist we check the derivative of squared Euclidean norm of m, p4:
1 d
mm - -Mel1 (m, - (1)) (in - j4:)
2 dt
(m: (1,) Q--Q-1 am 11) Imax
|| 11,112
57

---

Thus, ||m 14, S Mmaa
m=
Moll-
To bound ||Cr EellE F
define
S. = & (d, (1 (1++)Q)
This time we do not need the relative energy but can directly check the derivative of the squared Frobenius norm:
1 d |c, - 1.13 d 2dl
2 ir((C; - Ex) (C, - E,)(O, - En)
2C,,, Σ,}² 2(C, En)(C - =
22 E (C, E.)(C, [.)(C =2) (G. + =,)
2( Σ,)² - -C; E.) (C, + Ex) (10))-1-C0 - -CC, - -E, + E.)(C
ENO-"
(88)
where !f denotes equal in trace. Now if C, E, E S, then
(r((C - [., + E.)(C =) ir (Q 1/2 Ex)(C + E.)(C, - -1,)Q-1/2)
2 (( - c)---))222000 EnQ((G En - 2 2(1 - c) tr (CC, - -E,)))
(89)
and
tr((C, E-) + (G -11-)))))
tr ((C, - En)²)
=tr =tr((C, (C, - -()))=-)
1 tr ((C: Ex) (IG, E, - 2(1 - ()Q) + Q + E, 2(1 o)Q) (C. -x,) 2 0.
=2 2
(90)
Note that o is not trivially true. We show it as a lemma: Suppose A is a symmetric matrix and B is positive definite symmetric matrix. Then tr tr(ABBB-1) tr(A²).
In fact, we write B = PAPT where P is an orthogonal matrix and A = diag dd) is diagonal matrix. Then
tr tr(ABBB-))= = tr(APAPT APA-1PT) =
=|AA - APA APA-1/21
Denoting PT AP = (a:j) we get
2
VT
| APA 1P12-1111
1/2
Σ
dig
F
f,f=1
A
Σ
an ir(A)).
x
(,-=1
Thus, substituting 9 and (89)
into 888
we get
dे IG- -.- = 2tr (cC E (0,)) S -4(11 -26)]( - -2113 O.
d
Suppose at time to both Ct, and Eno lie in S., Then for any 2 to we have
||C, Eellr S se- 2(=) |1O, Meall 1r.
58


Reviews:
[{"review_id": "oZPkvVnIeF", "reviewer": "Reviewer_ZXYv", "summary": "This work studies the behavior of Gaussian-SVGD (SVGD projected to the space of Gaussian distributions, using bilinear kernels) and its variants.  The authors studied the behavior of the mean-field PDEs and established for Gaussian targets finite-particle convergence results which are significantly better than previous results for nonparametric kernels and targets. They have also derived density- and particle-based implementations of Gaussian-SVGD for general targets, which are shown to generalize recent works on Gaussian variational inference (GVI) with algorithmic guarantees.  One of the proposed algorithms outperform recent works on a Bayesian logistic regression experiment.", "questions": "For the experiments it would be helpful to compare with ordinary gradient descent on the variational parameters, to provide some ideas on the practical utility.", "limitations": "Yes.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The theoretical results appear interesting (to me as a non-expert, and note that I didn't check the proofs): the mean-field results sometimes improve over recent works on Gaussian variational inference with algorithmic guarantees, and the finite-particle results apply to standard SVGD, albeit with the other restrictions.\n\n- The discussions on connections between Gaussian SVGD and existing GVI approaches are also interesting, and one of the proposed methods appear promising empirically.", "weaknesses": "- While the restrictions to Gaussian families and affine kernels are similar to some of the recent works, it inevitably limits the scope of this work as much of the interest around SVGD arises from its flexibility.\n\n- While Gaussian variational inference in general has demonstrated competitive empirical performance, it is unclear if the proposed approach maintains this property, as it is not clear if the more complex algorithms interact nicely with common modifications such as minibatching.  This is also relevant because while some of the recent works (e.g., [17, 36]) also studied alternative GVI approaches, they come with full algorithmic guarantees, whereas for the present method there is no qualitative convergence guarantees for the discrete-time, finite-particle case.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the detailed review and incisive questions.\n\n1. **Restrictions to Gaussian families and affine kernels:**\n   - For a large class of Bayesian inference problems, thanks to the Bernstein von Mises theorem, the posterior (target) distribution is approximately Gaussian in the limit of large samples under appropriate regularity conditions. Notably a recent work [A] shows that in terms of mean and covariance estimation (of the posterior) using GVI has improved rates over other methods like Laplace approximation. There is an array of works providing **stochastic** algorithms (which is based on randomized discretizations of the BW Gradient flow) for solving GVI (see, e.g., [17, 36]), and [55, 68, 70, 57]  show the improved practical performance of GVI on various practical problems. We show in this work that there is a surprising **deterministic** discretization (BWPF) of the BW Gradient flow if we aim to solve GVI. Furthermore, we also have **quantitative rates** for mean and covariance estimation for **deterministic** SVGD algorithm in the finite-particle discrete-time setting (which we describe next), similar to the results for the **stochastic** algorithms mentioned above.\n\n[A] A Katsevich, P Rigollet. \"On the approximation accuracy of gaussian variational inference.\" \narXiv preprint arXiv:2301.02168, 2023.\n\n  - In terms of assumptions on kernel in prior works on SVGD, we also highlight that prior work by [49] requires that the kernel be radial which rules out the important class of bilinear kernels that we consider. The work of [27]  relaxed the radial kernel assumption. However, they required boundedness assumptions which we avoid in this work for the case of bilinear kernels. \n\n  **To sum up, understanding Gaussian-SVGD is an important step towards understanding the general situations with other kernels (note here that flexibilty offered by fully-nonparametric SVGD suffers from the issue of kernel choice). Surprisingly this special case has not been well studied in the literature and our work provides nearly-complete solution to this important special case.**\n\n2. **About discrete-time, finite-particle guarantees:** \n- We have an updated version of Theorem 3.9 (discrete-time, finite-particle) with the following **quantitative** convergence rates:\n\n     **Theorem:** For a centered Gaussian target, suppose the SVGD particle system with $K_{1}$ or $K_{2}$ is initialized by $\\bigl(\\boldsymbol{x}\\_{i}^{(0)}\\bigr)\\_{i=1}^{N}$ such that $\\boldsymbol{\\mu}\\_{0}=\\boldsymbol{0}$ and $C_{0}Q=QC\\_{0}$. For $0<\\epsilon< 0.5$, we have $\\boldsymbol{\\mu}\\_{t}=\\boldsymbol{0}$ and $\\lVert C_{t}- Q \\rVert\\to 0$ as long as all the eigenvalues of $Q^{-1}C\\_{0}$ lie in the interval $(0, 1+1/\\epsilon)$. Furthermore, if we set $u_{\\epsilon}$ to be the smaller root of the equation $f\\_{\\epsilon}'(u)=1-\\epsilon$ (it has $2$ distinct roots) where $f\\_{\\epsilon}(x) := (1 + \\epsilon(1-x))^2x$, then we have linear convergence, i.e.,  \n$$\n\\lVert C\\_{t}- Q \\rVert\\leq (1-\\epsilon)^{t}\\lVert C_{0}-Q \\rVert\\leq e^{-\\epsilon t}\\lVert C\\_{0}-Q \\rVert\n$$\n\nas long as all the eigenvalues of $Q^{-1}C\\_{0}$ lie in the interval $[u\\_{\\epsilon}, 1/3+1/(3\\epsilon)]$.\n\n\n  The above result is essentially direct refinement of the current result in Theorem 3.9. and is the first such result for finite-particle, discrete-time setting which matches the observed empirical performance. The only other comparable result for **deterministic** SVGD is by [62], where the results are not in alignment with practice. Obtaining similar results for general targets is extremely challenging and is one of the central open problems regarding SVGD. In view of this, our result provides a first concrete step towards developing a theory of SVGD that aligns with practice.\n\nWe also remark here that the fully algorithmic results for GVI in e.g., [17, 36] are for **stochastic** algorithms and are only proved in expectation. \n\n\n3. **About ordinary gradient descent:** We thank the reviewer for this great suggestion. We will be happy to add this result in our revision. \n"}, {"review_id": "Fb6X63mkYc", "reviewer": "Reviewer_1EPQ", "summary": "This work performs a theoretical investigation of Gaussian SVGD, a special case of SVGD restricted to the submanifold of Gaussian densities by means of a bilinear kernel. The authors characterize the mean-field dynamics of Gaussian SVGD, with a particular focus on Gaussian targets and obtain finite-particle guarantees in certain settings. Moving beyond Gaussian targets, the authors show that the dynamics of mean-field Gaussian SVGD converges to the best Gaussian approximation (In KL divergence) to the target density, and obtain convergence rates for the mean-field dynamics with the target density is log-smooth and log-concave", "questions": "1. Could you obtain discrete-time finite-particle rates for Gaussian SVGD under the settings considered in your results. If not, could you highlight what would be the technical challenges involved in obtaining such a result?\n\n\n2. Could you clarify why the mean-field continuous-time Gaussian SVGD dynamics is of interest, when (to the best of my knowledge) : 1) It does not satisfactorily explain the behavior of SVGD in practical scenarios, 2) It is outperformed by Bures-Wasserstein Gradient Flows (as well as the practically implementable Bures-JKO scheme) in most settings\n\n\n3. Could you comment on the technical novelty of these results? To the best of my understanding, the mean-field guarantees can be obtained by adapting the existing analysis of the Gradient Flow induced by SVGD [5,6] to the submanifold of Gaussian densities (a technique which is well-established for Wasserstein gradient flows in [3]). \n\n\n4. I would recommend reworking the overall presentation. Explicitly stating where each theorem is proven in the Appendix would be a good starting point. \n", "limitations": "I found the discussion on the limitations of this work to be somewhat unsatisfactory. The authors state that limitations are discussed on Section 6 (Other Related Works) and Section 7 (Conclusion). Section 6 discusses some of the prior literature and Section 7 highlights some future directions. Neither of these sections adequately discuss the limitations of this work (e.g. applicability is primarily restricted to centered Gaussian targets for most of the results, finite-particle discrete-time rates are absent in most settings)", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 1, "contribution": 2, "strengths": "The work provides a thorough treatment of the mean-field dynamics of Gaussian SVGD. The experimental evaluation is also satisfactory. ", "weaknesses": "My primary concerns are as follows :\n\n1. Most of the results presented in this work consider continuous-time mean-field dynamics under the assumption of a Bilinear kernel and Gaussian target density. In my opinion, these assumptions adversely impact the significance of these results, as they fail to explain one of the major appeals of SVGD, i.e., the fact that it provides a non-parametric approximation to a large class of target densities. (e.g. [1] shows discrete-time mean-field guarantees for subgaussian densities and [2] extends it to the class of densities satisfying a generalized Talagrand’s inequality). While the exponentially fast convergence guarantees for the mean-field dynamics are seemingly appealing, they fall short of providing a satisfactory explanation of the behavior of SVGD in practically relevant settings. \n\n\n2. Despite the simplifying assumption of Gaussian targets and bilinear kernels, the work does not provide any quantitative discrete-time finite-particle convergence rates. To the best of my understanding, Theorem 3.9 is the only result that considers finite-particle discrete-time dynamics, but it does not provide any quantitative rates.\n\n\n3. It is not clear why one should prefer Gaussian SVGD over the Bures-Wasserstein Gradient Flow [3] since the only setting where it seems to outperform Bures-Wasserstein Gradient Flow is that of centered Gaussians, which is a highly restrictive condition. \n\n\n4. The results in Section 4 leave much to be desired. Considering the fact that the best known results for Bures-Wasserstein gradient flow [3] cover both log-concave and strongly log-concave densities (see also [4] for a computable JKO discretization of the Bures-Wasserstein flow with convergence guarantees for log-concave and log-strongly concave densities), it is not at all clear what the benefits of Gaussian SVGD are when the target is not Gaussian. \n\n\n5. The overall presentation of the results requires a lot of work. For instance, it is not clearly stated where exactly in the Appendix each Theorem is proved. \n\n\n6. In addition to the limited applicability of these results in explaining the behavior of SVGD (as the most practically relevant case is that of non-logconcave targets) and the absence of satisfactory finite-particle discrete-time guarantees, I am somewhat concerned about the technical novelty of these results. To the best of my understanding, It seems like the mean-field continuous-time guarantees of this work can be easily obtained from the well-established results on SVGD [5,6] by restricting to the (finite-dimensional) Gaussian submanifold. While limited theoretical contribution in itself is not a major weakness, I find it difficult to recommend acceptance given the unsatisfactory practical usefulness of these results. \n\n\n[1] Salim et. al., “A Convergence Theory for SVGD in the Population Limit under Talagrand's Inequality T1”\n\n[2] Sun et. al., “Convergence of Stein Variational Gradient Descent under a Weaker Smoothness Condition”\n\n[3] Lambert et. al., “Variational inference via Wasserstein gradient flows”\n\n[4] Diao et. al., “Forward-backward Gaussian variational inference via JKO in the Bures-Wasserstein Space”\n\n[5] Duncan et. al., “On the geometry of Stein variational gradient descent”\n\n[6] Liu, “Stein Variational Gradient Descent as Gradient Flow”\n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "1.***Could you obtain discrete-time finite-particle rates***\n\nYes. We have the following **quantitative** result:\n\n**Theorem:** For a centered Gaussian target, suppose the SVGD particle system with $K_{1}$ or $K_{2}$ is initialized by $\\bigl(\\boldsymbol{x}\\_{i}^{(0)}\\bigr)\\_{i=1}^{N}$ such that $\\boldsymbol{\\mu}\\_{0}=\\boldsymbol{0}$ and $C_{0}Q=QC\\_{0}$. For $0<\\epsilon< 0.5$, we have $\\boldsymbol{\\mu}\\_{t}=\\boldsymbol{0}$ and $\\lVert C_{t}- Q \\rVert\\to 0$ as long as all the eigenvalues of $Q^{-1}C\\_{0}$ lie in the interval $(0, 1+1/\\epsilon)$. Furthermore, if we set $u_{\\epsilon}$ to be the smaller root of the equation $f\\_{\\epsilon}'(u)=1-\\epsilon$ (it has $2$ distinct roots) where $f\\_{\\epsilon}(x) := (1 + \\epsilon(1-x))^2x$, then we have linear convergence, i.e.,  \n$$\n\\lVert C\\_{t}- Q \\rVert\\leq (1-\\epsilon)^{t}\\lVert C_{0}-Q \\rVert\\leq e^{-\\epsilon t}\\lVert C\\_{0}-Q \\rVert\n$$\n\nas long as all the eigenvalues of $Q^{-1}C\\_{0}$ lie in the interval $[u\\_{\\epsilon}, 1/3+1/(3\\epsilon)]$.\n\nThis is a direct refinement of Theorem 3.9 and is the 1st such result for finite-particle, discrete-time setting matching practice. The only other rates (for deterministic SVGD) by [62] do not align with practice. Obtaining similar results for general targets is extremely challenging and is one of the central open problems in SVGD. \n\n2.**Bilinear kernel and Gaussian family assumptions**.\n\nNone of the existing results in the literature cover this simple case. [49] requires that radial kernels, and bilinear kernels are not. [27]  relaxed this. But they need boundedness, which is also not satisfied by bilinear kernels. \n\nFor Gaussian targets, we have shown uniform in time propagation of chaos in Theorem 3.7. Prior works only derived bounds for the empirical measure that **grow with time $t$**. This is unrealistic because SVGD only exactly recovers the target distribution or solves the GVI problem when it converges at $t=\\infty$, resulting in $\\infty$ in their bounds. \n\n **Understanding SVGD on Gaussian families is an important step towards understanding the general situation. This special case has not been well studied in the literature and our work provides nearly-complete solution to this important special case.**\n\n\n3.**Relation with BWGF**. There is a little bit misunderstanding here we guess and we apologize for lack of clarity previously.\n\n- Gaussian-SVGD is not 1 algorithm or 1 flow. **It is a family of flows with different bilinear kernels.** BWF is one of them (the one with kernel $K_3$, See Page 3). WGF and $K_3$-SVGD agree but only on the Gaussian submanifold and they induce BW flow. \n\n- We show each flow in this family induces two algorithms, one updates the density parameters ( so **stochastic**), the other updates the particles (so **deterministic**). The particle-based one is more stable in practice. **BWGD is the stochastic algorithm from BW flow appeared in [3] whereas BWPF is the particle-based one proposed by us.** It is hard to derive the deterministic particle-based counterpart of BWGD from the perspective of WGF. But from the perspective of $K_3$-SVGD, it is natural.\n\n4. **With all due respects, we simply disagree with the claim on lack of theory novelty.**\n\n- We prove the **first** result on **uniform in time** propagation of chaos in SVGD. Previously such results were not available even for the setting we consider. Our proof ideas are generalizable to other setting although admittedly more challenging.\n  \n- We explicitly solved the finite-particle dynamic system for Gaussian targets (See Theorem 3.6 and Theorem 3.8). Although there is no novelty in checking a given solution, **finding the solution is not trivial** in the first place.\n\n- We point out that our results with bilinear kernels **cannot be obtained as special cases of previous works on non-parametric SVGD**. [1] and [2] did not cover bilinear kernels. They require their kernels to be bounded and have bounded derivatives, which **is essentially required by their proof** and not satisfied by bilinear kernels. **Interestingly this is not an issue of technique but rather an issue of bilinear kernels.** The fact that Gaussian-SVGD using bilinear kernels is capable of solving GVI for general targets is **non-trivial** and that it provides deterministic algorithms for GVI is even more striking. No such **deterministic** particle method was available for GVI prior to our work. \n\n- For the mean field guarantee, [5] (Thm 22) or [6] is not applicable to our setting. Their condition cannot be easily checked in general (in fact not hold) for the Gaussian-SVGD setting we consider. Deriving tight results for special settings **rather than imposing unrealistic conditions for general settings** is a small but more solid step towards solving the general problem. \n\n- Directly applying the technique of [3] to SVGD would make the calculation too complicated to proceed. Therefore, we choose a different approach by calculating the explicit form of the restricted Riemannian metric tensor (e.g. See Theorem A.3). This approach is standard in Riemannian geometry but has added some novelty to our context. \n\n5. The results in **[3] and [4]** provide rates for discrete-time finite-particle algorithms in the log-concave and non-log-concave settings. However, their algorithm is stochastic and their results only hold in expectation. The fluctuations of their **stochastic algorithms** are not quantified.\n\nOur finite-time discrete-time result is an exact guarantee for the **deterministic** algorithm! Compared to [3] and [4], it holds for Gaussian targets. So our work complements the results in the works of [3] and [4]. Extending the results in [3] and [4] to quantifying the fluctuations and extending our results for non-Gaussian targets are both interesting future directions to make a fair comparison between these different algorithms. \n\n6.**Thanks for your suggestion. We will explicitly state where each theorem is proven in the main section of the paper.**"}, {"review_id": "S5Nbhv9Qfh", "reviewer": "Reviewer_Bt7r", "summary": "This paper studies the Stein variational gradient descent and its variants with a bilinear kernel on the space of Gaussian measures. The authors prove the rate of convergence of the dynamics, proposed finite particle algorithms and proved a uniform-in-time propagation of chaos,  and finally prove the convergence rate of the finite particle algorithms.", "questions": "1. Intuitive explanation of the convergence rate: as proved in the paper, the rate of convergence of SVGD for centered Gaussian is $O(e^{- t})$, which do not depend on $\\lambda$. This seems to be interesting, but I'm wondering if there's an intuitive explanation for this result? The dependence of $O(e^{- t/\\lambda})$ in the Wasserstein gradient flow case makes more sense to me intuitively, since larger $\\lambda$ will give a more \"flat\" potential.\n\n2. Regarding the non-Gaussian target:\n    - It is shown in Theorem 4.1 that when $\\rho_{\\theta^*}$ be the unique Gaussian measure that minimize $D_{KL}(\\rho_{\\theta^*}||\\rho^*)$, then the Gaussian-SVGD will converge to $\\rho_{\\theta^*}$. But when there are multiple $\\rho_{\\theta^*}$ that achieve the minimization, will the dynamic converge to one of them?\n\n    - Does the particle dynamics Algorithm 2 yield a similar uniform-in-time propagation of chaos results as in Theorem 3.7 with the same rate for Gaussian-SVGD?\n  \n\n\nMinor points and typos:\n\n1. Some ambiguity in definitions： \n    - The definition of the kernels on page 3 (lines 111-118) is not entirely clear to me, in particular, what is $\\mu$ and $\\Sigma$ here? If I understand the latter part correctly, I think $\\mu$ and $\\Sigma$ will be the mean and covariance of $\\rho_t$, which is actually time variant. I suggest the authors clarify this point in the definition.\n    - The same problem appears in equation (26) when defining $\\widehat{\\nabla V}$, since $\\widehat{\\nabla V}$ depends on $t$, which is a time-dependent vector field.\n\n\n2. Definition 2, line 101: $G_{\\rho}^{Wass}$ should be $G_{\\rho}^{Stein}$", "limitations": "No limitations and potential negative societal impact.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The study of the dynamics properties of  SVGD and its variants on Gaussian space seems to be interesting, which is a natural analog to the Bures-Wasserstein gradient flow. This paper gives a comprehensive analysis of such algorithms, including the well-posedness of the dynamics and the rate of convergence. \n\n2. The proof and the results seem to be correct. I didn't check all the details in the proof, but both the proof and the results make sense to me.  \n\n\n3. The paper is in general well-written, the theorems are stated properly, and I think it is not hard for readers to understand the main idea.", "weaknesses": "I do not observe an obvious weakness of this paper, but I do have a few questions, please refer to the \"Questions\" part.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the detailed review and incisive questions.\n\n1.**Intuitive explanation of the convergence rate:** \n\nThe intuition could be obtained from the looking at corresponding ODE arising in the analysis. The mean-field dynamics of SVGD with bilinear kernel (for $K_1, K_2$ or $K_3$) takes form of linear ODEs $\\dot{\\Sigma}_t \\approx - \\Sigma_t\\Sigma^{-1} $ where $\\Sigma$ is the scaling matrix of the kernel and $\\Sigma_t$ is the covaraince matrix of the solution of the mean-field dynamics of SVGD started from Gaussian initial data. The above linear ODEs warrant the rate of convergence of $\\Sigma_t$ to be $O(e^{-t})$ for $K_1,K_2$ and $O(e^{-\\frac{t}{\\lambda}})$ for $K_3$ since the scaling matrix $\\Sigma$ is equal to the identity matrix in the first two cases and equal to the covaraince matrix of the target distribution for the case of $K_3$.\n\nFurthermore, this also explains the fact that by choosing different regimes of $\\lambda$, regularized SVGD can interpolate between the two extremes, WGF and SVGD (See Table 1). \n\n2.**Regarding the non-Gaussian target:** \n\nIn the context of Theorem 4.1, when there are multiple $\\rho_{\\theta^*}$ that achieve the minimization, the dynamic will converge to one of them, but we can't distinguish where it converges as it depends on the initialization. We will clarify this further in the revision.\n\n3.**Uniform-in-time propagation of chaos for particle dynamics** \n\nCurrently, we don't have results of uniform-in-time propagation of chaos for Algorithm 2 with a general target, but we believe this should hold for the class of log-concave targets. Rigorously proving this is an interesting and challenging direction for future work.\n\n4.**Typos:** \n\nThanks very much for pointing out the typos. We will fix them, along with the ambiguities mentioned, in the revision.\n\n"}, {"review_id": "jWggddTdFR", "reviewer": "Reviewer_EWuB", "summary": "This paper delves into an examination of Gaussian-SVGD through the analysis of mean-field PDE and discrete particle systems. It offers evidence for finite particle convergence and supports these findings with empirical validations.\n", "questions": "Can you elaborate some advantages of Gaussian-SVGD compared with RBF-SVGD?\n", "limitations": "The setup of the article is somewhat artificial. It would be more meaningful if the significance of this linear system could be demonstrated.", "rating": 5, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "Through clear modeling, the article provides a comprehensive analysis for Gaussian-SVGD.\n\nThe paper gives new interpretations for Wasserstein Gradient Flow and SVGD under linear systems, which are very interesting.\n", "weaknesses": "The setup of the paper is relatively artificial, because in actual scenarios, linear kernels are rarely used to implement SVGD.\n\nThe paper does not put much effort into justifying the significance of Gaussian variational inference, such as making comparisons with other distribution classes. If there were clear conclusions, it would have a significant impact on the significance of this paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "1.**Linear kernels are rarely used to implement SVGD. It would be more meaningful if the significance of this linear system could be demonstrated. Signficance of GVI**\n\nFirst, we should emphasize again that Gaussian-SVGD is **NOT** SVGD with a linear kernel when the target is not Gaussian. We should use Gaussian-SVGD in order to perform GVI and cannot directly use RBF-SVGD for the same purpose. Moreover, there is no simple modification to make RBF-SVGD work for the GVI task. (The linear approximation in Equation (26) does not guarantee Gaussian dynamics for RBF kernel.) The linear kernel is used in our paper due to the fact that it leads to a **deterministic** algorithm that can be used to solve GVI. \n\n\nSecond, we elaborate on why GVI (or more generally variational inference) is preferred to exact sampling in many situations. The key argument is that **if one is only interested in summary statistics of the target such as the mean and covariance (as is often the case in Bayesian inference), generating samples may not be the most suitable way to achieve this goal**. In contrast to exact sampling, variational inference (VI), aims to find, among all measures in a certain parameterized family P, the closest measure to the target distribution. In general, VI algorithms usually show better performance both computationally and statistically than sampling algorithms. In specific, classical MCMC methods like MALA can be very expensive in terms of time cost, and it is notoriously difficult to identify clear-cut stopping criteria for the algorithm; The particle-based RBF-SVGD suffers from severe particle degeneracy especially in high dimensions [C] and cannot guarantee good covariance estimation.\n\nFor a large class of Bayesian inference problems, thanks to the Bernstein von Mises theorem, the posterior distribution is approximately Gaussian in the limit of large samples under appropriate regularity conditions. Thus, Gaussian VI is considered an important and basic task in this field. Notably a recent work [A] shows that in terms of mean and covariance estimation (of the posterior) using GVI has improved rates over other methods like Laplace approximation. In fact, there is an array of works providing **stochastic** algorithms (which is based on randomized discretizations of the BW Gradient flow) for solving GVI (see, e.g., [17, 36]), and [55, 68, 70, 57]  show the improved practical performance of GVI on various practical problems. We, however, present **deterministic** algorithms to do the same task, which are even more practically useful.\n\nFinally, the flexibliy offered by non-parametric nature has hindered researchers from obtaining theoretical results that align with practice. Given this, our contributions in this work as Reviewer BMUk pointed out provides \"a nice balance between establishing important theoretical results and providing a usable algorithm for users to implement in practice\". **Understanding SVGD on Gaussian families is an important step towards understanding the general situation. Surprisingly this special case has not been well studied in the literature and our work provides nearly-complete solution to this important special case.**\n\n \n[A] A Katsevich, P Rigollet. \"On the approximation accuracy of gaussian variational inference.\" arXiv preprint arXiv:2301.02168, 2023.\n\n[C] Zhuo, Jingwei, et al. \"Message passing Stein variational gradient descent.\" ICML, 2018.\n\nTo conclude, the **advantages of Gaussian-SVGD compared with RBF-SVGD** come in the following two perspective:\n\n- Gaussian-SVGD performs GVI, which is a more preferred approach for Bayesian inference.\n\n- Gaussian-SVGD provides a theoretically principled framework for understanding SVGD that **aligns with practice**. Gaussian-SVGD provides a **deterministic algorithm** to implement GVI. Uniform-in-time propogation type result (i.e., Theorem 3.7) are available for Gaussian-SVGD. Such results are not available currently for RBF-SVGD. \n\nIf a practitioner asks \"When should I use Gaussian-SVGD\", our paper provides a nearly-complete answer to this question. However, if a practitioner asks the question, \"When should I use RBF-SVGD\", currently no convincing answer could be given. The best we can think of is something along the lines of try it out and see if it works well, which is extremely ad-hoc. \n\n"}, {"review_id": "ccPmLEbjqo", "reviewer": "Reviewer_BMUk", "summary": "This paper focuses on the theoretical understanding and algorithmic contributions of Gaussian-Stein Variational Gradient Descent (Gaussian-SVGD) for Gaussian Variational Inference (GVI). The paper discusses the dynamics of Gaussian-SVGD and provides convergence rates in both mean-field and finite-particle settings when the target is Gaussian. It also shows that Gaussian-SVGD converges to the best Gaussian approximation to the target in KL-divergence for non-Gaussian targets. The authors propose two algorithmic frameworks for density-based and particle-based implementations of Gaussian-SVGD, which encompass various GVI algorithms proposed in the literature as special cases. The Gaussian-SVGD, and other related algorithms, are compared against one another on three examples: Gaussian target, Bayesian logistic regression, and a Gaussian mixture model. The authors highlight future directions for research, including establishing convergence results for particle-based Gaussian-SVGD algorithms under log-concave and more general targets, and exploring acceleration techniques. ", "questions": "* This sort of relates to my previous question approach motivation. It is well-known that SVGD does not scale well to high-dimensional target distributions. Is it possible to show, either theoretically (i.e. results like Theorem 3.7) or empirically, that Gaussian-SVGD is better than SVGD when $d$ is large?\n* Similarly, what is known about how many particles are needed compared to SVGD?\n* Particle-based methods are preferred over standard variational methods which fit a Gaussian distribution to the target because of their nonparametric nature. In the experiments section, it seems like the Gaussian mixture model would have been the more interesting example to include in the paper, particularly in the context of how good the approximation to the target is compared against SVGD. Could this be included? ", "limitations": "The paper could discuss more the theoretical limitation of this work. However, from an ethical perspective, I do not believe that this paper touches on issues that would have a negative societal impact. ", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper is very well written and the mathematical content is clearly articulated. The authors have done a good job of relating this work to other similar works from the literature. The paper is well-structured and the authors' particular contributions to this line of research are clearly outlined. \n\nThe technical detail, as far as I can see, is correct and the contribution of this work to the ML community, in particular the new convergence rates, are of significant interest. The paper is closely linked to other related works but is still highly original. \n\nThe paper provides a nice balance between establishing important theoretical results and providing a usable algorithm for users to implement in practice. ", "weaknesses": "This is a nice piece of work and there are no areas where I can see significant weaknesses. However there are a few things that the authors may wish to consider to improve their paper:\n\n* Starting from a very simple question, \"Why should I use Gaussian-SVGD?\" I think the authors have skipped some of the motivation for this work. If you know that the target is Gaussian, then of course you would just fit a Gaussian to it and would not require a particle-based variational inference approach. If the target is non-Gaussian, then why would you use Gaussian-SVGD over standard SVGD? \n* It is not clear why regularized SVGD is introduced. What is the motivation here? Is it simply to show that you can derive a Gaussian-SVGD algorithm with the regularized Stein metric?\n* The kernel $K(\\cdot,y)$ is introduced in Def 2.2 without being defined. Additionally, I think in eq. (1) you need to replace $\\rho$ with $\\rho_t$.\n* The simulation section is unfortunately a bit weak compared to other similar papers. The authors only consider three simple models and the details are largely glossed over, (e.g. what is $d$? How many particles? etc.). The authors seem to have run out of space and only provide a short simulation study for the logistic regression example. It would be better to include a more detailed simulation study in the main paper and to consider more challenging models.\n* Additionally, on the point of experiments, it is a bit surprising that SVGD is not also included in this list.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the detailed review and incisive questions.\n\n1.**About why we should use Gaussian-SVGD; \"If the target is non-Gaussian, why would you use Gaussian-SVGD over standard SVGD?\"**:\n\nWe interpret the question more as **when** we should use Gaussian-SVGD rather than **why**. It is important to note that **standard SVGD (e.g. using a bilinear or RBF kernel) does not solve GVI**. On one hand, if our goal is to find the best Gaussian approximation, we should no doubt use Gaussian-SVGD instead of SVGD. Furthermore, if the goal is to estimate the mean and variance from a target distribution, Gaussian-SVGD is still preferred because standard SVGD can subject to significant error in estimating the covariance especially when the dimension is high [B].\n\nHowever, as shown in [A] for a class of log-concave targets, samples generated via GVI (i.e., doing a Gaussian targets) would provide accurate estimators of the mean vector and the covariance matrix, and Gaussian-SVGD is essentially a deterministic algorithm to efficiently and accurately perform GVI. \n\nFor a large class of Bayesian inference problems, thanks to the Bernstein von Mises theorem, the posterior distribution is approximately Gaussian in the limit of large samples under appropriate regularity conditions. In fact, there is an array of works providing **stochastic** algorithms (which is based on randomized discretizations of the BW Gradient flow) for solving GVI (see, e.g., [17, 36]), and [55, 68, 70, 57]  show the improved practical performance of GVI on various practical problems. We show in this work that there is a surprising **deterministic** discretization (BWPF) of the BW Gradient flow if we aim to solve GVI. \n\nFurthermore, when one refers to standard SVGD, one also needs to specify the choice of kernels. In a way, our work suggests a principled way of selecting the kernel (which could be one of the bi-linear kernels) in standard SVGD for provably estimating the posterior mean and covariance for log-concave targets. \n  \n[A] A Katsevich, P Rigollet. \"On the approximation accuracy of gaussian variational inference.\" arXiv preprint arXiv:2301.02168, 2023.\n\n[B] Ba, Jimmy, et al. \"Understanding the variance collapse of SVGD in high dimensions.\" ICLR. 2021.\n\n\n2.**About why regularized SVGD is introduced**:  What we want to show is that the regularized Stein metric with kernel $K_2$ coincides with the Stein metric with a different bilinear kernel, which we can call as $K_4$. This phenomenon is interesting and also makes the algorithm fall under our general framework of Gaussian-SVGD. We will modify the paragraph on \"Three Bilinear Kernels\" on Page 3 to include $K_4$. Moreover, $K_4$ is important because it interpolates between $K_2$ and $K_3$. Note that $K_2$ and $K_3$ both have advantages and disadvantages, and $K_4$ can strike a balance between them by choosing $\\lambda$ (See Table 1 and Figure 1).\n\n3.**Is Gaussian-SVGD is better than SVGD when d is large? and What is known about how many particles are needed compared to SVGD?**\n\nThese two are very deep question, actually (from a theoretical perspective). For the first question, we do observe empirically that Gaussian-SVGD is better than SVGD (say with Gaussian or Matern kernels) when d is large, from our preliminary experiments in terms of estimating posterior covariances of log-concave targets. The intuitive reasoning is the parametric nature of SVGD incurs significantly less variance compared to nonparametric SVGD. \n\nIn an ongoing work, we are attempting to have a high-dimensional analysis (i.e., double asymptotic analysis) of Gaussian SVGD. Comparing to general SVGD theoretically in this regime, and providing a theoretical answer to the second question, are harder because we currently don't have tools to analyze SVGD well, even in fixed or low-dimensional setup. Nevertheless, a positive answer to these question would be great, and we hope to get there someday!\n\n4.**About Gaussian mixture model, simulation section and SVGD**:\n\nFirst, we would like to point out that we have a Gaussian mixture model experiment in Section D of the appendix.  Taking your suggestion, we plan to move all experiments to the main paper if we are allowed an additional page. If not, we will have the Gaussian mixture model in the main draft and other experiments in the appendix.\n\nOur primary goal is to make progress towards understanding SVGD theoretically, given its widespread usage in practice. Nevertheless, we will be happy to add more numerical experiments, with different dimensions for the current models and also some new models.\n\nFor the experiments, when one refers to SVGD it is important to specify which kernel to use. Given that you simply to refer to SVGD, we take it that you refer to some non-bi-linear kernels? That is a great suggestion, and we will be happy to add a comparison to SVGD with a Gaussian or Matern kernel in our experiments. \n\n\n\n\n"}]

Author rebuttal:
General Clarifications
=================

1.**Quantitative results for finite particle, discrete-time setting** 

We have an updated version of Theorem 3.9 (discrete-time, finite-particle) with the following **quantitative** convergence rates:

 **Theorem:** For a centered Gaussian target, suppose the SVGD particle system with $K_{1}$ or $K_{2}$ is initialized by $\bigl(\boldsymbol{x}\_{i}^{(0)}\bigr)\_{i=1}^{N}$ such that $\boldsymbol{\mu}\_{0}=\boldsymbol{0}$ and $C_{0}Q=QC\_{0}$. For $0<\epsilon< 0.5$, we have $\boldsymbol{\mu}\_{t}=\boldsymbol{0}$ and $\lVert C_{t}- Q \rVert\to 0$ as long as all the eigenvalues of $Q^{-1}C\_{0}$ lie in the interval $(0, 1+1/\epsilon)$. Furthermore, if we set $u_{\epsilon}$ to be the smaller root of the equation $f\_{\epsilon}'(u)=1-\epsilon$ (it has $2$ distinct roots) where $f\_{\epsilon}(x) := (1 + \epsilon(1-x))^2x$, then we have linear convergence, i.e.,  
$$
\lVert C\_{t}- Q \rVert\leq (1-\epsilon)^{t}\lVert C_{0}-Q \rVert\leq e^{-\epsilon t}\lVert C\_{0}-Q \rVert
$$

as long as all the eigenvalues of $Q^{-1}C\_{0}$ lie in the interval $[u\_{\epsilon}, 1/3+1/(3\epsilon)]$.


This is a direct refinement of Theorem 3.9 and is the 1st such result for finite-particle, discrete-time setting matching practice. The only other rates (for deterministic SVGD) by [62] do not align with practice. Obtaining similar results for general targets is extremely challenging and is one of the central open problems in SVGD. 


2.**Gaussian-SVGD versus bilinear SVGD**

We like to clarify that our algorithm (i.e, Gaussian-SVGD) is **different** from SVGD with a bilinear kernel, when the target is not Gaussian.  Standard (bilinear or RBF) SVGD can sample from the target but does not solve Gaussian variational inference (GVI) whereas Gaussian-SVGD can solve GVI. To achieve this, in Equation (26) we use a linear approximation $\widehat{\nabla V}$ instead of the actual $\nabla V$. (Note that this linear approximation does not guarantee Gaussian dynamics for RBF kernel.)

- If one wants to perform GVI, one essentially need to use Gaussian-SVGD instead of bilinear SVGD or RBF-SVGD.
  
- To estimate the posterior mean or covariance of a **not necessarily Gaussian target**, one can use either Gaussian-SVGD or standard (bilinear or RBF) SVGD. The particle-implementation of Gaussian-SVGD provides a **deterministic** algorithm for GVI, and is shown to perform well in our experiments. Furthermore, [A] shows the GVI provably does better than many other existing methods like Laplacian approximation. In contrast, there are no such theoretical guarantees for bilinear SVGD, RBF-SVGD or general SVGD in this settings. 


[A] A Katsevich, P Rigollet. "On the approximation accuracy of gaussian variational inference." arXiv preprint arXiv:2301.02168, 2023.

3.**Further clarifications from prior works**  

We now explain further the advantages of studying the Gaussian-SVGD and the different types of results available in the literature on SVGD. First note that the flexibility offered by the \emph{nonparametric} aspect of SVGD also leads to unintended consequences. Indeed, on one hand, from a practical perspective, the question of how to pick the right kernel for implementing the SVGD algorithm is unclear. Existing approaches are mostly ad-hoc and do not provide clear instructions on the selection of kernels. On the other hand, developing a deeper theoretical understanding of SVGD dynamics is challenging due to its nonparametric formulation. [49] derived the continuous-time PDE for the evolving density that emerges as the mean-field limit of the finite-particle SVGD systems, and shows the well-posedness of the PDE solutions. Furthermore, the following different types of convergences could be examined regarding SVGD, some of which has been analyzed previously in the literature:
  
- (A) Unified convergence of the empirical measure for $N$ finite particles to the continuous target as time $t$ and $N$ jointly grow to infinity;
  
- (B) Convergence of mean-field SVGD to the target distribution over time;
- (C) Convergence of the empirical measure for finite particles to the mean-field distribution at any finite given time $t\in [0,\infty)$;
- (D) Convergence of finite-particle SVGD to the equilibrium over time;
- (E) Convergence of the empirical measure for finite particles to the continuous target at time $t=\infty$.

From a practical point of view (A) is the ideal type of result that fully characterizes the algorithmic behavior of SVGD, which could be obtained by combining either (B) and (C) or (D) and (E). Regarding (B), [44] showed the convergence of mean-field SVGD in kernel Stein discrepancy, which is known to imply weak convergence under appropriate assumptions. The works of [19],[34],[60] and [65] sharpened the results with weaker conditions or explicit rates. The work [27]  extended the above result to the stronger Fisher information metric and Kullback-Leibler divergence based on a regularization technique. The works [44] and [34] obtained time-dependent mean-field convergence (C) of $N$ particles under various assumptions using techniques from the literature of propagation of chaos. The work of [62] obtained even stronger results for (C) and combined (B) to get the first unified convergence (A) in terms of KSD. However, they have a rather slow rate $1/\sqrt{\log\log N}$, resulting from the fact that their bounds for (C) still depends on the time $t$ (sum of step sizes) double-exponentially. Moreover, there has not been any work that studies the convergence (D) and (E) for SVGD, which illustrate a new way to characterize the unified convergence (A).

4.**Prior work** by [49] requires that the kernel be radial which rules out the important class of bilinear kernels that we consider. The work of [27]  relaxed the radial kernel assumption. However, they required boundedness assumptions which we avoid in this work for the case of bilinear kernels. 

Meta-review: There has been thorough discussion of this paper in the reviewing process. The main object being studied is a sampling algorithm called Gaussian SVGD, i.e., SVGD with bilinear kernels and Gaussian initialization. There are two main theoretical contributions:

* A **continuous-time** analysis of **mean-field** Gaussian-SVGD for general targets (linear convergence for both gaussian and strongly log-concave targets)
* Both continuous and discrete-time analysis of **finite-particle** Gaussian-SVGD for **Gaussian targets** (quantitative rate provided in rebuttal)

While concerns were raised regarding the finite-particle analysis being specific to Gaussian targets, the solid theoretical contribution of the first part analyzing general targets satisfied the reviewers. Additionally, it was noted that the uniform-in-time propagation of chaos results, though for Gaussian targets, improve upon prior work. Based on these, I recommend acceptance.
Decision: Accept (poster)