Year: 2021

Abstract:
Stein variational gradient descent (SVGD) is a deterministic inference algorithm that evolves a set of particles to fit a target distribution. Despite its computational efficiency, SVGD often underestimates the variance of the target distribution in high dimensions. In this work we attempt to explain the variance collapse in SVGD. On the qualitative side, we compare the SVGD update with gradient descent on the maximum mean discrepancy (MMD) objective; we find that the variance collapse phenomenon relates to the bias from deterministic updates present in the “driving force” of SVGD, and empirically verify that removal of such bias leads to more accurate variance estimation. On the quantitative side, we demonstrate that the variance collapse of SVGD can be accurately predicted in the proportional asymptotic limit, i.e., when the number of particles n and dimensions d diverge at the same rate. In particular, for learning high-dimensional isotropic Gaussians, we derive the exact equilibrium variance for both SVGD and MMD-descent, under certain empirically verified near-orthogonality condition on the converged particles, and confirm that SVGD suffers from the “curse of dimensionality”.

heading: 1 INTRODUCTION
text: A typical challenge in Bayesian learning is to efficiently and accurately learn a complex posterior distribution. Markov Chain Monte Carlo (Robert and Casella, 2013) provides asymptotically accurate samples, but simulating the chain till convergence can be time-consuming. On the other hand, variational inference (Wainwright et al., 2008) approximates the intractable target distribution with parametric families by minimizing the KL divergence, which can be more computationally efficient. Among a plethora of variational approximations, particle-based inference algorithms are particularly compelling as they bypass the inherent distributional assumptions on the variational family.
As a noticeable example, Stein Variational Gradient Descent (SVGD) (Liu and Wang, 2016) is a deterministic particle-based inference algorithm that iteratively transports the particles by the functional gradient of KL divergence in the reproducing kernel Hilbert space (RKHS). The functional gradient takes the form of a kernelized Stein’s operator and only requires access to the unnormalized target density. Despite the empirical successes of SVGD (Liu, 2017; Haarnoja et al., 2017; Kim et al., 2018), very few convergence guarantees have been established except in the mean-field limit, i.e., the number of particles n → ∞ under fixed dimensionality d, for which it has been shown that the distribution of particles converges to the true invariant solution (Liu, 2017; Lu et al., 2019). Moreover, it has been observed that as the problem dimensionality d becomes larger, the variance estimated by SVGD can be much smaller than variance of the target distribution (Zhuo et al., 2017). This observation, which we refer to as the variance collapse phenomenon (see Figure 1), is highly undesirable for practitioners due to two reasons: (i) underestimating the variance leads to failure in explaining the uncertainty of model predictions, which is a key benefit of being Bayesian; (ii) modern Bayesian inference problems are usually high-dimensional; for instance, training Bayesian neural networks (BNNs) (MacKay, 1992) requires inferring the posterior distribution of network weights, which could be more than millions of dimensions in real-world problems (Krizhevsky et al., 2012).
Our Contribution. We study the algorithmic bias of SVGD that leads to variance collapse in high dimensions. We focus on the commonly-used Euclidean distance kernel (e.g., Gaussian RBF kernel), and provide: (i) qualitative understanding on the pitfall of SVGD by comparing it to another interacting particle algorithm termed MMD-descent; (ii) quantitative understanding on the variance underestimation in learning simple target distributions. Our findings can be summarized as follows.
• In Section 3 we connect SVGD with MMD-descent, a kernel-based particle inference algorithm that performs maximum mean discrepancy (MMD) minimization, and empirically show that despite their similar updates, MMD-descent does not collapse the variance in high dimensions.
• In Section 4 we identify the log derivative driving force as the problematic term in SVGD, and design experiments to illustrate that the combination of this driving force and the deterministic bias, i.e. the absence of particle resampling, prevents accurate variance estimation in high dimensions.
• In Section 5 we argue that the proportional asymptotic limit: n, d → ∞, d/n → γ ∈ (0,∞), is more relevant to understanding the variance collapse phenomenon. As an example, we derive the exact dimension-averaged variance of SVGD and MMD-descent in learning Gaussian distribution in the proportional limit, under certain concentration assumption which we empirically verify. Our analysis confirms that variance estimated by SVGD scales inversely with the dimensionality.

heading: 2 BACKGROUND

heading: 2.1 INTEGRAL PROBABILITY METRIC
text: To measure the “closeness” between distributions on X ⊆ Rd, one may consider the maximum discrepancy between the target p and sample distribution q over test functions F : DF (p, q) := suph∈F Eq[h(x)]− Ep[h(y)], known as the integral probability metric (IPM) (Müller, 1997). If F is a unit ball in the reproducing kernel Hilbert space (RKHS)H, the resulting DF is termed the maximum mean discrepancy (MMD) (Gretton et al., 2012), and its squared value can be evaluated as: MMD2(p, q) = Ex,x′ [k(x,x′)] + Ey,y′k(y,y′)− 2Ex,y[k(x,y)] (1) where x,x′ ∼ p, y,y′ ∼ q, and the kernel k : Rd × Rd → R satisfies E √ k(x,x) <∞.
Stein’s Discrepancy. When integration under the target p is intractable, Stein’s method (Stein et al., 1972) can be used to construct zero-mean test functions w.r.t. p. For differentiable h in the Stein Class of p, i.e., ∫ x ∇x(h(x)p(x))dx = 0, the Stein’s discrepancy (Gorham and Mackey, 2015) is given as DSteinF (p, q) = suph∈F Eq[h(x)>∇x log p(x) +∇>xh(x)]. Note that the Stein’s discrepancy only involves the score of p and thus the normalization constant is not required. When h is restricted in the product RKHS Hd with inner product 〈f, g〉Hd =∑d i=1〈fi, gi〉H, the corresponding maximum discrepancy, known as kernel Stein discrepancy (KSD), can be estimated efficiently from samples (Liu et al., 2016; Chwialkowski et al., 2016).

heading: 2.2 (DETERMINISTIC) PARTICLE INFERENCE ALGORITHMS
text: We now consider the approximation of p using a set of particles constructed as follows: starting from initial particles X = {xi}ni=1, we iteratively optimize their positions using the deterministic update: xi = xi + η∆(xi), where η is the step size, and ∆(·) : Rd → Rd represents the update direction. Motivated by the SVGD algorithm, we focus on the setting where the update ∆ is deterministic, in which case particles may converge to a deterministic fixed point. Note that there are other particlebased inference algorithms that do not fit into this description, such as sequential Monte Carlo (Doucet et al., 2001), and stochastic variants of particle gradient method (Chen et al., 2018a; Gallego and Insua, 2018), which our current analysis does not cover (see Appendix A.4 for more discussion).
Stein Variational Gradient Descent. SVGD constructs the update direction ∆ as the optimal perturbation in the RKHS that decreases Kullback-Leibler divergence. In particular, constrain ∆ in RKHS unit ball, and take q = 1n ∑n i=1 δxi , then the update for each particle x can be computed as:
∆SVGD(x) = Ex′∼q[k(x′,x)∇x′ log p(x′)︸ ︷︷ ︸ driving force +∇x′k(x′,x)︸ ︷︷ ︸ repulsive force ] := 1 n n∑ i=1 [S1(xi,x) + S2(xi,x)]. (2)
Intuitively, the log derivative term in the update rule S1(xi,x) := k(xi,x)∇xi log p(xi) corresponds to a driving force that guides particles towards high likelihood regions, whereas the kernel derivative term S2(xi,x) := ∇xik(xi,x) provides a repulsive force to prevent the particles from collapsing into the mode. Typically-used kernels in SVGD include the Gaussian RBF kernel k(x,x′) = exp(−‖x− x′‖22 /2σ2) (Liu and Wang, 2016; Zhuo et al., 2017), the linear kernel k(x,x′) = x>x′ + 1 (Liu and Wang, 2018), and the inverse multi-quadratic (IMQ) kernel k(x,x′) = 1/ √ 1 + ‖x− x′‖22/(2σ2) (Gorham and Mackey, 2017).

heading: 3 CONNECTING SVGD WITH MMD MINIMIZATION
text: In this section, we introduce another particle inference algorithm termed MMD-descent, whose update rule closely relates to that of SVGD. Despite the similarity, we empirically observe MMD-descent accurately estimates the target variance independent of the dimensionality. By analyzing the similarity of the update rules and contrast in the algorithmic performances, we identify factors that lead to the variance collapse of SVGD. In this work we focus on the following class of kernel functions. Definition 1 (Euclidean Distance Kernel). A positive semi-definite kernel is referred to as a Euclidean distance kernel if it can be represented as: k(x,y) = f ( ‖x− y‖22/σ2 ) .
For example, the Gaussian RBF kernel and IMQ kernel both satisfy this definition. Note that σ is the tunable bandwidth that usually scales with the distance between particles.
We now construct a kernel-based particle inference algorithm termed MMD-descent, which draws inspiration from the kernel herding algorithm (Welling, 2009) and particle gradient descent (Chizat and Bach, 2018). Instead of greedily reducing the MMD by adding one particle a time, we consider an update rule similar to SVGD that transport all particles together to approximate the target distribution, by minimizing the MMD between the particles and the target via gradient descent:
∆MMD(x) = Ey∼p[∇xk(x,y)︸ ︷︷ ︸ driving force ] + Ex′∼q[−∇xk(x,x′)︸ ︷︷ ︸ repulsive force ].
This update can be seen as the finite-particle discretization of the MMD gradient flow (Arbel et al., 2019). For Euclidean distance kernels we have −∇xk(x,x′) = ∇x′k(x,x′) = S2(x′,x); hence, MMD-descent and SVGD share the same repulsive force. Furthermore, when k is in the Stein class of p, the driving force of MMD-descent can be connected SVGD via a simple integration by parts:
Ey∼p[∇xk(x,y)] = −Ey∼p[∇yk(x,y)] = Ey∼p[k(x,y)∇y log p(y)] = Ey∼p[S1(y,x)]. (3) Therefore, the MMD-descent update for a set of particles {xi}ni=1 can be equivalently written as:
∆MMD(x) = Ey∼p[S1(y,x)] + 1
n n∑ i=1 S2(xi,x). (4)
We remark that MMD-descent is not a practical algorithm due to the required integration under p, which is typically intractable. Instead, the purpose of introducing MMD-descent is to compare the update with SVGD and understand the cause of the variance collapse phenomenon.
SVGD vs. MMD-descent. By comparing SVGD (2) and MMD-descent (4), we observe that:
• SVGD and MMD-descent have identical repulsive force.
• In MMD-descent, the driving force is integrated under the target distribution p, whereas in SVGD the expectation is under the current particle distribution q.
At the infinite particle limit, p = q is a fixed point for both update rules. In fact, under certain conditions, uniqueness of the fixed point has been established for both updates in this mean-field limit (Lu et al., 2019; Arbel et al., 2019). However, in the practical setting where the particle size is not significantly greater than the dimensionality, it is not clear if the two algorithms: (i) provide a reasonable “approximation” to the target distribution; (ii) converge to similar fixed points. Given the analogous updates, one natural question to ask is: do SVGD and MMD-descent reliably estimate the target variance in high dimensions, and do they give similar solutions?
The empirical answer is in the negative: SVGD and MMD-descent converge to completely different solutions: SVGD is known to underestimate the marginal variance in high dimensions, even for simple Gaussian targets (Zhuo et al., 2017). In contrast, MMD-descent does not exhibit the same issue empirically. Figure 1 illustrates this discrepancy in a Bayesian neural network (BNN) experiment. We consider a two-hidden-layer BNN with 100 hidden units each, so that inference in the weight space is a high-dimensional problem. We obtain (approximate) target samples via Hamiltonian Monte Carlo (HMC) (Neal et al., 2011), and then initialize 100 particles from the target to be optimized by SVGD or MMD-descent with the Gaussian RBF kernel (for detailed setup see Appendix D). Observe that although we used the same number of particles for SVGD and MMD-descent, SVGD severely underestimates the variance, whereas MMD-descent generates diverse predictions similar to HMC.

heading: 4 UNDERSTANDING THE PITFALL OF SVGD
text: We now qualitatively characterize the algorithmic bias of SVGD. Motivated by the observation that SVGD and MMD-descent only differs in the driving force which involves S1, we first compare S1 and S2 in SVGD and argue that the former term is more likely problematic in high dimensions due to larger fluctuation. We then design controlled experiments by modifying the updates of SVGD and MMD-descent to show that when the log derivative S1 is coupled with the deterministic bias (i.e., absence of particle resampling), the algorithm does not reliably estimate the target variance.

heading: 4.1 HIGH VARIANCE FROM INTEGRATION BY PARTS
text: Recall that the update rule of SVGD involves a driving force (log derivative) term 1n ∑n i=1 S1(xi,x)
and a repulsive force (kernel derivative) term 1n ∑n i=1 S2(xi,x). From integration by parts we have the following equality: Ey∼p[S1(y,x)] = −Ey∼p[S2(y,x)]. Thus the target distribution p is a fixed point of SVGD when the number of particles tends to infinity. However, in the finite sample setting, the fluctuations in S1 and S2 might have impact on the converged solution.
In many settings, we expect the log derivative driving force S1 to have higher variance. For instance, in the case of unit Gaussian target p and Gaussian RBF kernel k(x,x′) = exp(−‖x−x′‖22/2σ2), we know that S1(y,x) = −yk(y,x) and S2(y,x) = σ−2(x−y)k(y,x). Roughly speaking, when the norm of x,y and the bandwidth σ scale with √ d, ‖S1‖22 increases with the dimensionality, whereas ‖S2‖22 remains bounded. The large magnitude of S1 results in large fluctuation in the driving force term in SVGD. In Figure 2(a), we visualize the distributions of S1 and S2 when the particles x are drawn i.i.d. from Gaussian distributions, and indeed observe the higher variance of S1. The following proposition quantifies this discrepancy in an idealized setting.
Proposition 2. Define the mean squared error as: MSEp[f(y)] = Ep(y) ‖f(y)− Ep[f(y)]‖ 2 2. Then for y ∼ N (a, Id) where ‖a‖22 = O(d), Gaussian RBF kernel k(x,y)=exp(−‖x−y‖ 2 2 /2σ
2) with bandwidth σ = Θ( √ d), and x ∈ Sd−1( √ d), we have
MSEp[S2(y,x)] = Θ(d −1); MSEp[S1(y,x)] = Θ(d).
More generally, for x ∈ Sd−1( √ d), strongly log-concave p, and Euclidean distance kernel with Lipschitz f and lower-bounded by a scaled Gaussian kernel, MSEp[S1(y,x)] = Ω(d) ·MSEp[S2(y,x)].
The proposition indicates that when particles are i.i.d. sampled from some target p, then S1 would have much larger fluctuation than S21. Intuitively speaking, the higher variance in S1 suggests that more samples are required to estimate the term accurately. Therefore, when the dimensionality is large compare to the particle size, poor estimation of S1 may relate to the variance collapse. The following subsection empirically demonstrate that this is indeed the case; in particular, we show that in the presence of deterministic bias, the high-variance S1 leads to “biased” particles in SVGD.

heading: 4.2 BIAS FROM DETERMINISTIC UPDATE
text: In this subsection we isolate a cause of variance collapse that relates to the deterministic nature of the update rule. Observe that in the derivation of SVGD, particles {xi}ni=1 are assumed to be sampled
1Note that particles optimized by SVGD cannot be considered as i.i.d. samples from some distribution, and thus Proposition 2 does not rigorously apply to particles obtained by the actual algorithm.
from some underlying distribution q. But due to the deterministic update, q is entirely represented by the same set of particles, and redrawing i.i.d. samples is not feasible during optimization. We now demonstrate that this failure to resample, which we refer to as the deterministic bias, when combined with the log derivative driving force S1, may cause the algorithm to converge to the wrong target.
Deterministic Bias in MMD-descent. We start with an illustration of deterministic bias in MMDdescent. Given samples {yi}mi=1 ∼ p drawn from the target, we have two distinct MMD-gradient updates (connected via integration by parts (3)) that differs in the driving force.
∆MMD1 (x) = 1
m m∑ i=1 S1(yi,x)+ 1 n n∑ i=1 S2(xi,x). ∆ MMD 2 (x) = − 1 m m∑ i=1 S2(yi,x)+ 1 n n∑ i=1 S2(xi,x).
We refer to ∆MMD1 and ∆ MMD 2 as the log derivative update and the kernel derivative update, respectively. Section 4.1 suggests that ∆MMD1 tend to have larger fluctuation than ∆ MMD 2 . For both update rules, we consider the two variations: (i) fixed particles, where we draw {yi}mi=1 in the beginning and use this same set of target samples throughout optimization; (ii) resampled particles, where we redraw new i.i.d. samples {yi}mi=1 from p at each iteration. In this construction, the fixed-particle update (i) emulates the deterministic bias in S1, due to the absence of resampling from p.
We simulate the two variants of MMD-descent in Figure 2(b), in which the target is a 50-dimensional unit Gaussian, n = 50 and we vary m from 50 to 1000. First observe that the log derivative update with fixed particles (solid red), which combines the high variance S1 and the deterministic bias, results in diverged particles. In contrast, the log derivative update with resampled particles (dashed red), in which S1 is present but not the deterministic bias, converges to the desired variance. Moreover, if the driving force is estimated by S2, then the variance is accurate even in the presence of deterministic bias, as shown in the kernel derivative update with fixed particles (solid blue). This indicates that the combination of deterministic bias and the log derivative driving force S1 leads to biased convergence2.
Algorithm 1 SVGD with Particle Resampling Input: Initial density q0. Number of Steps S. for s = 1 to S do
Sample {xsi}ni=0 ∼ q0(x). for t = 1 to s− 1 do
either Compute S1 from resampled particles: xsi ← xsi+ ηn ∑n j=1 S1(x t j ,x s i )+S2(x s j ,x s i ). or Compute S2 from resampled particles: xsi ← xsi+ ηn ∑n j=1 S1(x s j ,x s i )+S2(x t j ,x s i ).
end for end for Output {xSi }ni=1. Particle Resampling in SVGD. The previous experiment on MMD-descent suggests that the deterministic bias arises from the algorithm not being able to “redraw” i.i.d. samples. We further validate this observation in SVGD by constructing a variant that achieves particle resampling. The modified update exhibits a doubleloop structure (see Algorithm 1): at iteration i of the outer loop, we obtain a new set of particles indexed as qi after (i−1) inner loop steps. In the inner loop, we first draw n i.i.d. particles q̂i from the initial distribution q0, and then update q̂i via (i−1) SVGD steps to obtain qi; at each step j of the inner loop, either the driving force S1 or the repulsive force S2 is computed between q̂i and qj , as opposed to between particles within q̂i as in the original SVGD.
2Note that due to the different driving force than SVGD, the log derivative update of MMD-descent with fixed particles results in divergence instead of variance collapse (see Appendix A.3 for more discussion).
We remark that this resampled update resembles transport-based particle algorithms (e.g. Nitanda and Suzuki (2017)), where old particles define a transport map that the new particles (initialized from q0) follow. By construction, at each outer loop iteration, a new set of independent particles are generated and updated; importantly, one of S1 and S2 is not evaluated on the same set of particles themselves, and hence deterministic bias is not present. However, the complexity of such resampled updates scales quadratically with the iterations, which renders the algorithm computationally prohibitive3.
Figure 2(c) compares the modified SVGD updates (with Gaussian RBF kernel) in learning a Gaussian target. Observe that when the driving force S1 is computed on resampled particles (blue), SVGD accurately estimates the target variance even with small particle size (n = 20) at each step. In contrast, if only the repulsive force S2 is evaluated on particles resampled from q0, then the variance is still underestimated (green), due to deterministic bias in S1. This confirms our argument that combination of log derivative S1 and deterministic bias leads to the algorithmic bias of SVGD.

heading: 5 SVGD IN THE PROPORTIONAL LIMIT
text: While the previous discussion provides qualitative understanding of the algorithmic bias of SVGD, no quantitative characterization is provided. In this section, we make use of the deterministic nature of the SVGD update, and directly analyze the particle fixed point; this allows us to precisely compute the variance of SVGD particles and confirm the variance collapse in simple settings.

heading: 5.1 SCALING LIMIT AND BASIC ASSUMPTIONS
text: Previous works have considered the infinite-particle limit of SVGD under fixed dimensionality, which is known as the mean-field limit (Lu et al., 2019; Duncan et al., 2019). In this regime, the particle dynamics of SVGD is described by a nonlinear PDE: ∂tρ = ∇ · (ρ(K ∗ (∇ρ + ∇V ρ))), which converges weakly to the unique target measure as t→∞ under certain non-degeneracy assumptions. However, since the invariant solution is the desired target distribution, the mean-field limit does not capture the variance collapse phenomenon. Moreover, similar guarantees can also be derived for MMD-descent under similar assumptions (Arbel et al., 2019). Therefore, the mean-field limit does not explain the observed discrepancy between SVGD and MMD-descent (Figure 1 in Section 3).
Instead, we propose to analyze the algorithm in the proportional asymptotic limit:
• (A1) Proportional Limit. n, d→∞ and d/n→ γ; we also assume γ > 1. Note that larger γ is analogous to higher dimensionality or smaller particle size; we focus on the γ > 1 regime, where the number of particles is less than the problem dimensionality — this is a common feature of modern Bayesian inference problems. In this asymptotic limit, prior works studied spectral properties of the kernel matrix (e.g., El Karoui et al. (2010)), which the SVGD update crucially depends upon. However, these random matrix results typically require the particles to be i.i.d., which is not satisfied by SVGD due to the interacting update. To overcome this technical difficulty, we make an additional assumption on the fixed point:
• (A2) Near-orthogonality. The optimized particles {xi}ni=1 satisfy |x>i xi−dυ| < υ d, |x>i xj | < υ d for all i 6= j, some υ > 0, and d−1/2 d → 0 as d→∞ with probability 1.
The above assumption simply ensures that the particles tend to equilibrate at a fixed point where they “repel” one another as much as possible (due to the repulsive force), as shown in Figure 3 for the case of Gaussian target (for more empirical verification see Appendix A.2). We emphasize that this is a condition only on the fixed point of the SVGD algorithm, not on the particles produced by the SVGD iterations (unless initialized at the fixed point).
Under assumptions (A1) and (A2), we can approximate the Gram matrix via Taylor expansion of the kernel matrix, which allows us to simplify the fixed point equation. We restrict ourselves to the following class of kernels:
3Resampling can also be achieved via kernel density estimation (KDE) at each iteration, as in Dai et al. (2016). However, KDE is also known to be less effective in high dimensions.
• (A3) Kernel Function. k(x,y) = f ( ‖x−y‖22
2σ2
) , where f : R≥0 → R+ is bounded, monotonically
decreasing, and differentiable4 on R≥0. Remark. (A3) covers common choices of Euclidean distance kernels with decaying tail, such as the Gaussian RBF kernel, the IMQ kernel and the log-inverse kernel. Our analysis also allows for different choices of bandwidth σ, which we specify in the sequel.

heading: 5.2 LEARNING HIGH-DIMENSIONAL GAUSSIANS
text: We now consider a simple example in which SVGD does not reliably estimate the variance – learning an isotropic Gaussian target (empirically observed in Zhuo et al. (2017)). Note that our goal is to demonstrate a negative result: we do not claim that SVGD underestimates the variance in all settings; instead, we show that variance collapse is present even in learning the simplest target distribution.
• (A4) Gaussian Target. p(x) ∝ exp ( − 12x >x ) .
Remark. We assume (A4) for simple and concise presentation, and we believe that the result can be extended to more general sub-Gaussian targets. In addition, since Gaussian prior is widely used in Bayesian inference, it is natural to expect that target potentials in many real-world problems are “Gaussian-like” (i.e., exhibit quadratic growth) outside of some radius (Cheng et al., 2018a). Indeed, our empirical observations on BNN align with our theoretical predictions for Gaussian target.
We aim to compute the dimension-averaged marginal variance: v = 1d ∑d j=1 Varj({xi}ni=1), where Varj(·) is the particle variance at the j-th coordinate. In the following subsections we derive the equilibrium variance of SVGD and MMD-descent under two different choices of kernel bandwidth.
Kernels with Adaptive Bandwidth. We first consider the case where the bandwidth σ is adaptively tuned based on the optimized particles. In particular, we analyze the median heuristic (Scholkopf and Smola, 2001): σ = √
Med{‖xi − xj‖22}/2, which is the most common choice made in practice. Under the previous assumptions, we have the following result for SVGD with median bandwidth.
Proposition 3. Given (A1-4) and the choice of bandwidth σ = √
Med{‖xi − xj‖22}/2, particles driven by SVGD (Equation (2)) equilibriate at the dimension-averaged variance
vSVGD = γ−1 · f ′(1) [f(1)− f(0)]−1 . Remark. The proposition precisely captures our observation that SVGD underestimates the dimension-averaged variance inversely proportional to the problem dimensionality (γ = d/n). In contrast, previous analysis (Zhang et al., 2020) predicts that SVGD would collapse to zero variance as t→∞, which does not align with the empirical observations.
For Gaussian kernel, we also provide a precise contrast between SVGD and MMD-descent. Corollary 4. Given (A1-4), for the Gaussian RBF kernel f(x) = exp(−x) with the median bandwidth, SVGD converges to the following dimension-averaged variance: vSVGD = (e−1)−1 ·γ−1 < 1, whereas MMD-descent (Equation (4)) leads to vMMD = 1.
In other words, as the problem dimensionality γ increases, more particles are required for SVGD to reliably estimate the dimension-averaged variance (specifically, the number of particles should grow linearly with d). In contrast, the variance estimated by MMD-descent remains accurate and is independent to γ5; this aligns with the empirical observations in Section 3. Remarkably, Figure 4(a) demonstrates that Proposition 3 is accurate even for reasonably small n, d: observe that once γ > 1 (i.e., d > n), the prediction (black) becomes well-aligned with the empirical value.
Kernels with Fixed Bandwidth. One possibility remains that the culprit of variance collapse is the median bandwidth. Specifically in kernel Stein discrepancy (KSD), it has been argued that the IMQ kernel can outperform the Gaussian RBF kernel (Gorham and Mackey, 2017), and the IMQ kernel is often employed without an adaptive bandwidth. Due to the connection between KSD and SVGD, one may speculate that the IMQ kernel with fixed bandwidth can alleviate the variance collapse problem.
4For the median bandwidth, we only require local differentiability around 1, similar to El Karoui et al. (2010). 5In Appendix A.4 we demonstrate that Langevin Monte Carlo also does not suffer from variance collapse.
We show that this is unfortunately not the case. We first provide a general characterization: Proposition 5. Given (A1-4) and fixed bandwidth6 σ = √ d, the SVGD variance satisfies
f ′(vSVGD) = γ · [ f(vSVGD)− f(0) ] .
In addition, if f ′ is also monotone on R≥0, then vSVGD decreases as γ > 1 increases.
Note that the monotonicity assumption is again satisfied by many Euclidean distance kernels of interest. While the equation may not provide an explicit expression of the stationary variance, we can verify that the variance decreases as the problem becomes more high-dimensional (i.e., larger γ). Our next proposition specifically handles the IMQ kernel considered in Gorham and Mackey (2017). Corollary 6. Given (A1-4), for the IMQ kernel f(x) = (1 + x)−1/2 with fixed bandwidth, we have the following stationary variance of SVGD under two different scalings: • When σ = √ d, vSVGD < 1 and is decreasing as γ > 1 increases.
• When σ = 1, vSVGD → 0 as n, d→∞ at a rate of d−1/3.
This corollary suggests that the IMQ kernel with fixed bandwidth is not a remedy to the variance collapse problem. We remark that the second setting (σ = 1) is considered in Gorham and Mackey (2017) for KSD. In both cases (and including the median bandwidth covered by Proposition 3), SVGDIMQ underestimates the target variance as the dimensionality increases (large γ). The agreement between the theoretical predictions and the empirical simulations (using finite particles) is illustrated in Figure 4(b)(c). Finally, in Appendix A we include a more general (but less precise) characterization of variance collapse beyond Gaussian target, as well as additional empirical evidence.

heading: 5.3 A MODIFICATION OF SVGD
text: We have thus far shown that in the simple setting of learning high-dimensional Gaussian, SVGD underestimates the dimension-averaged variance unless the number of particles is larger than the dimensionality. Now we further validate our theoretical findings by introducing a heuristic modification of SVGD that corrects for this variance collapse in the overparameterized regime7.
The starting observation is that the variance collapse indicates that the deterministic bias causes the driving force term to dominate. Because during each update every particle xi is most “correlated” with itself, one should expect S1(xi,xi) to contribute significantly to this bias. We thus consider a modification of SVGD which simply shrinks (damps) the term S1(xi,xi) by λ = min{1, (f(1)− γ−1f ′(1))/f(0)}, where λ is chosen such that when d > n, the equilibrium variance matches the target variance in the setup of Proposition 3. We refer to this update as damped SVGD:
∆Damp(xi) = ∑ j 6=i [S1(xj ,xi) + S2(xj ,xi)] + λS1(xi,xi). (5)
Remark. At high level, (5) resembles the annealed SVGD algorithm (D’Angelo and Fortuin, 2021), in that they both ”weaken” the driving force S1. However, annealed SVGD uses a heuristic learning rate schedule (target-independent) to modify the strength of S1, whereas we derive our λ to ensure that the dimension-averaged variance is correct under (A1-4) in the d > n regime.
6Note that the fixed bandwidth can be arbitrary, as constants can be absorbed into the function f . 7The modified update in Section 4 also alleviate the variance collapse, but is computationally intractable.
Simulating Toy Densities. We demonstrate the usefulness of this bias-correction term in learning simple target distributions. We stack 250 independent 2- dimensional distributions to create a high-dimensional target distribution, and visualize the first two dimensions in Figure 5. We optimize 50 particles using SVGD or SVGD (damped) with the Gaussian RBF kernel and median bandwidth.
As shown in Figure 5 (top), SVGD collapses its particles to the mode, whereas damped SVGD generates more dispersed particles with variance exactly matching that of the Gaussian target. It is also possible that the heuristic modification leads to more accurate samples for other “Gaussian-like” target distributions8, as shown in the second row of Figure 5. We provide additional experimental results in Appendix A.3.

heading: 6 RELATED WORKS
text: Guarantees and Applications of SVGD. Liu and Wang (2018); Lu et al. (2019) characterized SVGD in the mean-field limit and showed the weak convergence to the target distribution. Under further assumptions on the potential and kernel, quantitative convergence rate can be derived (Duncan et al., 2019; Korba et al., 2020). In addition, Liu and Wang (2018) shows SVGD using kernels with finite-dimensional feature maps, exactly estimates the expectations for some set of functions, casting SVGD as a moment matching method. In the high-dimensional setting, Zhuo et al. (2017); Wang et al. (2018) observed in experiments that particles driven by SVGD tend to underestimate the marginal variance, but did not provide any quantitative understanding of the phenomenon.
On the application side, Haarnoja et al. (2017); Liu (2017) adopt SVGD to learn a stochastic sampling network to approximate the policy in Q-learning (Sutton et al., 1998), whereas in Gangwani et al. (2018), SVGD encourages diverse policies for exploration. SVGD can be used in meta-learning to quickly obtain parameter samples from training sets (Yoon et al., 2018). Recent works also applied SVGD in Batch Bayesian optimization (Gong et al., 2019) and the learning of mixture models (Wang and Liu, 2019). Leveraging the Markov blanket structure, SVGD also achieves strong performance in learning graphical models (Zhuo et al., 2017; Wang et al., 2018).
Stein’s Method. Stein’s method provides powerful tools in approximating probability distributions and specifying convergence rates (Erdogdu, 2016; Gorham et al., 2016). Liu and Wang (2016) utilizes the connection between Stein’s operator and the gradient of KL divergence to construct particle inference algorithm. Stein’s lemma is also useful in implicit variational inference (Huszár, 2017) to estimate the score using samples from an implicit distribution (Li and Turner, 2017; Shi et al., 2018). Related to our analysis in Section 4, Erdogdu et al. (2016) observed that algorithms that are equivalent in expectation via Stein’s lemma might have different convergence properties. In addition, the “curse of dimensionality” of kernel Stein estimators has also been studied in Oates et al. (2016).

heading: 7 CONCLUSION
text: We analyzed the variance collapse of SVGD in high dimensions based on a connection between SVGD and a proposed MMD-descent algorithm. We qualitatively identified factors that lead to this phenomenon, and also quantitatively characterized the equilibrium variance in the proportional limit for simple models. Looking forward, we believe that understanding interacting particle systems (SVGD, two-layer neural nets, etc.) in the proportional limit instead of the mean-field limit is an important direction. In addition, while our analysis confirms the variance collapse of the original algorithm, it remains possible that certain dimension reduction schemes could alleviate the issue. Chen and Ghattas (2020); Gorham et al. (2020); Gong et al. (2020) recently proposed dimension-reduced versions of SVGD, and it would be interesting to theoretically justify these approaches.
8We however note that due to the distribution-specific derivation of λ, we should not expect the proposed modification to be a general solution to the variance collapse problem of SVGD.

heading: 1 Introduction 1

heading: 2 Background 2
text: 2.1 Integral Probability Metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2 (Deterministic) Particle Inference Algorithms . . . . . . . . . . . . . . . . . . . . 2

heading: 3 Connecting SVGD with MMD Minimization 3

heading: 4 Understanding the Pitfall of SVGD 4
text: 4.1 High Variance from Integration by Parts . . . . . . . . . . . . . . . . . . . . . . . 4
4.2 Bias from Deterministic Update . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

heading: 5 SVGD in the Proportional Limit 6
text: 5.1 Scaling Limit and Basic Assumptions . . . . . . . . . . . . . . . . . . . . . . . . 6
5.2 Learning High-dimensional Gaussians . . . . . . . . . . . . . . . . . . . . . . . . 7
5.3 A Modification of SVGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

heading: 6 Related Works 9

heading: 7 Conclusion 9

heading: A Additional Results 14
text: A.1 Additional Results for Section 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.2 Empirical Verification of (A2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
A.3 Additional Results for Modified SVGD . . . . . . . . . . . . . . . . . . . . . . . 16
A.4 Discussion on Other Sampling Algorithms . . . . . . . . . . . . . . . . . . . . . . 17
A.5 Additional Figures for Section 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

heading: B Proof of Propositions in Section 4 19

heading: C Proof of Propositions in Section 5 20
text: C.1 Gaussian Kernel with the Median Bandwidth . . . . . . . . . . . . . . . . . . . . 21
C.2 General Kernel with Adaptive Bandwidth (Median Heuristic) . . . . . . . . . . . . 23
C.3 General Kernel with Fixed Bandwidth . . . . . . . . . . . . . . . . . . . . . . . . 24
C.4 “Almost-Gaussian” Target . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
D Experiment Setup 25

heading: A ADDITIONAL RESULTS
text: A.1 ADDITIONAL RESULTS FOR SECTION 5
We include additional figures to illustrate our analytic results in Section 5. In Figure 6 we plot the stationary variance of SVGD and MMD-descent (with the Gaussian RBF kernel and median bandwidth heuristic) under fixed dimensionality d against varying number of particles n. As shown in Figure 6(a), the particle variance of SVGD scales linearly with 1/d when γ > 1 (i.e., d > n), as predicted by Proposition 3 and Corollary 4; this again confirms the variance collapse phenomenon. When n > d, our analysis cannot predict the stationary variance accurately, but the value empirically approaches the target variance from below as γ decreases.
Non-asymptotic Correction for MMD-descent. As noted in Section 5, there is a small gap between the empirical particle variance of MMD-descent and the prediction of Proposition 3, due to finite-particle error in the empirical simulation. To characterize this discrepancy, we now consider a regime where n is fixed and d is large. In this case, we show that SVGD collapses all the particles to 0, whereas MMD-descent with the Gaussian RBF kernel can still estimate the marginal variance up to O(1/n) error (see Appendix C for derivation). Proposition 7. Consider the setting of fixed n and d→∞, then given (A2)(A3), particles driven by SVGD (2) collapses to zero variance: vSVGD → 0, whereas MMD-descent (4) with the Gaussian RBF kernel leads to vMMD → (n− 1)/(n+ 1).
Figure 6(b) demonstrates the agreement between the proposition and finite-particle behavior of MMDdescent. Observe that the finite-particle error decays at 1/n and is independent to the dimensionality. In other words, our analysis suggests that to achieve non-trivial estimation of the target variance, SVGD requires Θ(d) particles, whereas for MMD-descent Od(1) particles suffice.
Based on this observation, we speculate that to match the accuracy of MMD-descent with n particles (in terms of marginal variance), SVGD requires N = nd particle. This conjecture is empirically confirmed in Figure 6(c): when we increase the number of particles in SVGD by an additional factor of d, then the equilibrium variance on longer depends on the dimensionality and becomes similar to that of MMD-descent with n particles. This again suggests that SVGD may not be reliable unless the number of particles scale with the problem dimensionality.
Variance Collapse in Almost-Gaussian Target. In Section 5 we focus on Gaussian target for precise characterization of the stationary variance. For more general target distributions, if the potential behaves “Gaussian-like” outside of a certain radius (smaller than the true target variance), then under the same assumption on the fixed point particles (A3), we intuitively expect SVGD to underestimate the marginal variance when γ = d/n is large. Proposition 8. Let p(x) ∝ exp(−f(x)) be an isotropic density, i.e., Ep[x] = 0,E[xx>] = Id, and assume f exhibits Gaussian tail growth outside a Euclidean ball with radius 1, then given (A1-3), there exists γ∗ = d/n > 1 such that SVGD underestimates the dimension-averaged variance for γ > γ∗.
While the proposition above only considers quadratic growth of the potential (similar to Mattingly et al. (2002)), we empirically verify that SVGD underestimates the target variance under anisotropy as well as different tail growth condition. In Figure 7(a) we construct the target distribution to be a mixture of two Gaussians with mean µ1, µ2 such that µ1 + µ2 = 0, ‖µ1 − µ2‖2 = 1, whereas in Figure 7(b) we consider a factorized distribution with a cubic-growth potential: log p(x) ∝ − ∏d i=1(x
(i))3/3. In both cases we observe similar variance collapse phenomenon of SVGD across different choices of kernel and bandwidth.
Log-inverse Kernel. We also illustrate the variance collapse phenomenon of SVGD with the loginverse kernel proposed in (Chen et al., 2018b). This heavy-tailed kernel is defined as k(x,x′) =( α+ log(1 + ‖x− x′‖22 /σ2) )−1 . One can easily verify that it is a Euclidean distance kernel satisfying assumption (A2), and thus we also expect variance collapse to occur in learning the unit Gaussian target (A3). This is indeed confirmed in Figure 7(c)(d), where we see that the log-inverse kernel with both the median heuristic and fixed bandwidth (σ = 1) underestimates the marginal variance, and in the case of median bandwidth, we can analytically predict the equilibrium variance for γ > 1 via Proposition 3.
A.2 EMPIRICAL VERIFICATION OF (A2)
We empirically validate the near-orthogonality assumption in Section 5 in learning Gaussian targets. In particular, we plot the following quantities for particles driven by SVGD and MMD-descent till convergence (from random initialization):
(i) d−1/2maxi{‖xi‖22 − υd}; (ii) d −1/2maxi,j{x>i xj}.
Note that for (A2) to hold, the above quantity should decay with the dimensionality d.
In Figure 8 we validate this hypothesis for both SVGD and MMD-descent on two different kernels: the Gaussian RBF kernel and the IMQ kernel, and two target distributions: unit Gaussian and mixture of Gaussians (MoG). For the Gaussian RBF kernel we use the median bandwidth heuristic, and for the IMQ kernel we set the bandwidth σ = √ d. The MoG model consists of two Gaussians with mean µ1 +µ2 = 0, ‖µ1 − µ2‖2 = 1. For both algorithms we initialize the particles fromN (0, 2Id). Observe that in all cases the quantities of interest indeed decrease with d. We leave the rigorous justification of this observation as future work.
A.3 ADDITIONAL RESULTS FOR MODIFIED SVGD
In Section 5 we proposed a modification of SVGD based on derivations of learning a unit Gaussian target distribution in high dimensions (and kernels with the median bandwidth),
∆Damp(xi) = ∑ j 6=i [S1(xj ,xi) + S2(xj ,xi)] + λS1(xi,xi),
where λ = max{1, (f(1)− γ−1f ′(1))/f(0)}. Note that λ = 1 recovers the original SVGD update, whereas under extreme overparameterization (p n) we have λ = f(1)/f(0) < 1 due to (A3). In addition, larger λ implies smaller stationary variance of SVGD, vice versa. In the figures below we refer to the SVGD update with λmin = min{1, (f(1)− γ−1f ′(1))/f(0)} as “fully-damped”, and the update with λ in between λmin and 1 as “intermetidate”.
Bayesian Neural Network. We now apply this modified update to the BNN setting discussed in Section 3 (for detailed setup see Section D). We initialize 100 particles from a standard normal distribution and then optimize them via (modified) SVGD with the Gaussian RBF kernel and median bandwidth for 50000 steps with learning rate η = 5e− 3. As shown in Figure 9, the original SVGD update (Figure (b), original) significantly underestimates the target variance (estimated using HMC), and as we decrease λ, the variance of the predictions gradually increases (Figure (c), intermediate). Observe that at our proposed damping factor (Figure (d), fully-damped), the modified update produces diverse samples similar to HMC.
In real-world problems, it is not likely that our proposed modification (derived for Gaussian target) always results in optimal performance. This being said, Figure 9 suggests that practitioners could tune the damping term λ in the range [f(1)/f(0), 1] and adjust the desired level of uncertainty. Theoretical justification of our modified update is left as future work.
Bayesian Logistic Regression. We also conduct additional experiment of the modified (damped) SVGD update in a Bayesian logistic regression problem. Given input Z = {zi}mi=1 ∈ Rm×d, labels y = {yi}mi=1 ∈ Rm, and parameter θ ∈ Rd, we model the Bernoulli conditional distribution with probability Pr(yi=1|zi) = 1/(1 + exp(−θ>zi)). We place an isotropic Gaussian prior on θ; the posterior density is given as
p(θ) ∝ exp ( y>Zθ −
m∑ i=1 log(1 + exp(−θ>zi))− α 2 ‖θ‖22
) .
Following Dalalyan (2014), we sample the coordinates of zi from a Rademacher distribution and then normalize the vector by its Euclidean norm; the labels are generated from a Bernoulli distribution with true parameters θ∗ = 1d; we set m = 500, d = 100, and the regularization parameter α = 1.
We report the dimension-averaged marginal variance of the particles (in the n < d regime) optimized by variants of SVGD in Figure 10(b). We use the Gaussian RBF kernel with the median bandwidth, and for the ”intermediate” update we set λmin = min{1, (f(1) − 2γ−1f ′(1))/f(0)}. The ground truth variance is estimated by Langevin Monte Carlo (see Appendix A.4 for more discussion). Figure 10(b) illustrates that for this problem, the original SVGD update (red) underestimates the target variance, similar to the Gaussian case. However, our heuristic modification (fully damped,
blue) does not remove this bias completely – this is not surprising because the correction term is derived only for the isotropic Gaussian target. In particular, the damped update might overestimate the dimension-averaged variance in the small-particle regime. Nevertheless, we observe that the proposed modification generates more diverse samples (with variance closer to the true value).
To further quantify the convergence of the particles, in Figure 10(a) we report the MMD between the SVGD particles and (approximate) ground-truth samples obtained from LMC; we employ the V-statistics MMD with the IMQ kernel (fixed bandwidth). Observe that the modified update (blue, green) results in smaller MMD compared to the original algorithm (red).
Complete Removal of S1(xi,xi). Recall that the SVGD update can be decomposed into three terms: 1) driving force from each particle itself S1(xi,xi); 2) driving force from other particles ∑ j 6=i S1(xj ,xi); 3) repulsive force ∑ j S2(xj ,xi) (note that S2(xi,xi) = 0). Specifically, S1(xi,xi) has a large impact on the update and the variance collapse: in the modification above, we are able to generate diverse samples just by decreasing the strength of this single term.
One might speculate that instead of using the specific λ derived for Gaussian target, an easier alternative would be to completely remove S1(xi,xi), i.e., setting λ = 0. This is however not the case. Following the same derivation as in Appendix C, one can see that when λ = 0, Equation 13 does not equilibriate at any finite variance for γ > 1. We empirically validate this finding in Figure 10(c), where we fix d = 50 and vary the particle size n from 25 to 500, and optimize the particle using (modified) SVGD with the Gaussian RBF kernel and the median bandwidth. Observe that while the original SVGD (red) underestimates the target variance, setting λ = 0 (blue) leads to significant overestimation or even divergence when d > n.
One minor remark is that heuristically speaking, SVGD with λ = 0 roughly resembles MMD-descent with fixed log derivative S1 (see Figure 2(b) in Section 4), due to the presence of deterministic bias, and that the driving force in both updates involves the log derivative S1, but does not include the “self” term S1(xi,xi). As a result, we observe an overestimation of marginal variance (or diverging particles) in both modified updates.
A.4 DISCUSSION ON OTHER SAMPLING ALGORITHMS
In addition to SVGD and MMD-descent, there are many other approaches to draw (approximate) samples from a given target distribution. We briefly discuss two popular methods: Langevin Monte Carlo (LMC) and Herding. In particular, we show that LMC does not suffer from variance collapse in high dimensions; in fact, under some conditions, one particle suffices to estimate the dimensionaveraged variance of the target.
Langevin algorithm. Langevin Monte Carlo (LMC) is a time discretization of the Langevin diffusion, which can be interpreted as the (Wasserstein) gradient flow of relative entropy (Jordan et al., 1998). Convergence rate of the discrete-time algorithm has been analyzed under various metrics (Dalalyan, 2014; Durmus and Moulines, 2017; Erdogdu et al., 2021) and assumptions on the target distribution (Vempala and Wibisono, 2019; Erdogdu and Hosseinzadeh, 2020).
Here we provide a short sketch that existing quantitative convergence guarantee of LMC implies that the algorithm does not underestimate the dimension-averaged variance. For strongly convex and smooth potentials (which includes the Gaussian setting in Section 5), Erdogdu et al. (2021, Theorem 4) showed that when the strong convexity parameter is dimension-free, then under appropriate choice of step size, running LMC for T = Õ(d) iterations leads to χ2(qT ||p) = O(1). We denote the covariance Ep[xx>] = Σ, and assume Ep[x] = 0,Tr(Σ) = d WLOG. By properties of the χ2-divergence and strongly log-concave concentration, given one sample xq drawn from qT ,
Pr (∣∣ ‖xq‖2 − Ep ‖x‖2 ∣∣ > t) . exp (−t2) .
This concentration directly implies that in the high-dimensional limit d→∞ (Assumption 1), we have ‖xq‖22 /Tr(Σ)
p→ 1. In other words, the dimension-averaged variance Tr(Σ)/d can be reliably estimated using one single particle driven by LMC for Õ(d) steps. We remark that similar results are expected to hold true under more general conditions such as the log-Sobolev inequality on the target distribution, as well as other variants of the algorithm (Roberts and Tweedie, 1996; Cheng et al., 2018b). We therefore speculate that variance collapse is a unique property of certain interacting particle algorithms such as SVGD.
Kernel herding. Consider the approximation of an intractable distribution p(x) with a set of particles X = {xi}ni=1. To generate these particles, the kernel herding algorithm was introduced by Welling (2009) to minimize the MMD between the particles and the target distribution. The algorithm proceeds in a greedy manner: given the current set of selected particles {x1, · · · ,xn−1}, the next particle is chosen based on the following:
xn ← argmin x
MMD2 ( p, 1
n ( n−1∑ i=1 δxi + δx )) = argmax x Ey∼p[k(x,y)]− 1 n n−1∑ i=1 k(x,xi).
Intuitively, the first term encourages sampling in high density areas for the target, whereas the second term discourages drawing samples close to existing ones. It has been shown that the kernel herding algorithm reduces the MMD at a rate O( 1N ) for finite-dimensional Hilbert spacesH (Welling, 2009; Bach et al., 2012; Huszár and Duvenaud, 2012).
Note that the same procedure can also be used to greedily minimize the kernel Stein discrepancy (Chen et al., 2018b; 2019). In addition, the non-greedy gradient descent update on KSD (analogous to MMD-descent) has been recently analyzed (Korba et al., 2021). We leave the high-dimensional characterization of these algorithms as future work.
A.5 ADDITIONAL FIGURES FOR SECTION 4
In Section 4 we conducted experiments that qualitatively identified the cause of variance collapse, primarily for the Gaussian RBF kernel. Here we reproduce these findings for the IMQ kernel: k(x,y) = 1/ √
1 + ‖x− y‖22 /2σ2), which is also commonly-used in SVGD. As shown in Figure 11, results are qualitatively similar to that of Gaussian RBF kernel.
Finally, we provide additional empirical evidence that SVGD with proper resampling procedure (see algorithm 1) reliably estimates the target variance. We consider the cubic-growth potential: p(x) ∝ exp ( − ∏d i=1(x (i))3/3 )
, and optimize the particles using SVGD with the Gaussian RBF kernel (fixed bandwidth σ = √ d). As shown in Figure 12, the original SVGD update (red) underestimates the dimension-averaged marginal variance inversely proportional to the problem dimensionality; in contrast, when S1 is resampled (blue), then the estimated variance is stable as d increases.

heading: B PROOF OF PROPOSITIONS IN SECTION 4
text: Proof of Proposition 2. For simplicity we first assume the target is zero-centered: y ∼ N (0, Id). In this case when the kernel is Gaussian RBF, the expectation of S1, which we denote as µx, has the following closed-form: µx = − σ d
(1+σ2)d/2+1 exp
( − ‖x‖ 2 2
2+2σ2
) x. We compute the mean squared error
of interest as follow: MSEp[S1(y,x)] = ∫ y ‖S1(y,x)− µx‖22 p(y)dy
= ∫ y ‖ − yk(x,y)− µx‖22 p(y)dy
= ∫ y k2(x,y)yTy p(y)dy + 2µTx ∫ y k(x,y)y p(y)dy + µTxµx ∫ y p(y)dy
= e − ‖x‖ 2 2 2+σ2 σd
(2 + σ2)d/2+1 (2‖x‖22 + dσ2)−
2e − ‖x‖ 2 2 1+σ2 σ2d (1 + σ2)d+2 ‖x‖22 + e − ‖x‖ 2 2 1+σ2 σ2d (1 + σ2)d+2 ‖x‖22
=‖x‖22 [2e− ‖x‖222+σ2 σ2 + 2 ( σ2 2 + σ2 )d/2 − e− ‖x‖221+σ2 (σ2 + 1)2 ( σ2 1 + σ2 )d] + de − ‖x‖ 2 2 2+σ2 ( σ2 2 + σ2 )d/2+1 .
Similarly, for the kernel derivative S2 we have, MSEp[S2(y,x)] = ∫ y ‖S2(y,x)− (−µx)‖22 p(y)dy
= ∫ y ‖x− y σ2 k(x,y) + µx‖22 p(y)dy
= 1
σ4 ∫ y k2(x,y)yTy p(y)dy − 2x σ4 ∫ y k2(x,y)y p(y)dy − 2µx σ2 ∫ y k(x,y)y p(y)dy
+ xTx
σ4 ∫ y k2(x,y) p(y)dy + 2xTµx σ2 ∫ y k(x,y) p(y)dy + µTxµx ∫ y p(y)dy
= e − ‖x‖ 2 2 2+σ2 σd−4
(2 + σ2)d/2+1 (2‖x‖22 + dσ2)−
4e − ‖x‖ 2 2 2+σ2 σd−4 (2 + σ2)d/2+1 ‖x‖22 + 2e − ‖x‖ 2 2 1+σ2 σ2d−2 (1 + σ2)d+2 ‖x‖22
+ e − ‖x‖ 2 2 2+σ2 σd−4
(2 + σ2)d/2+1 (2 + σ2)‖x‖22 −
2e − ‖x‖ 2 2 1+σ2 σ2d−2
(1 + σ2)d+1 ‖x‖22 +
e − ‖x‖ 2 2 1+σ2 σ2d
(1 + σ2)d+2 ‖x‖22
= e − ‖x‖ 2 2 2+2σ2 σ4 ( σ2 2 + σ2 )d/2+1 (d+ ‖x‖22) + 3e − ‖x‖ 2 2 1+σ2 (1 + σ2)2 ( σ2 1 + σ2 )d‖x‖22.
The simplification above largely follows from Ex∼N (µ,Σ)[‖x‖22] = µTµ+ Tr(Σ). For non-centered y ∼ N (a, Id), note that replacing ‖x‖22 with ‖x− a‖ 2 2 does not affect the order dependence on d as long as ‖a‖22 = O(d), and therefore the order that we aim to estimate remains the same. Given the bandwidth heuristic σ = Θ( √ d) and ‖x‖22 = d, one can easily verify that:
MSEp[S2(y,x)] ∈ Θ(d−1), MSEp[S1(y,x)] ∈ Θ(d).
To extend the result to general Euclidean distance kernel k(x,y) = f ( ‖x−y‖22
2σ2
) and strong log-
concave distributions satisfying ∥∥Ep(y)[y]∥∥22 = O(d) and Ep(y)[‖y‖22] = O(d), first note that
Ey∼p[‖µx‖2] = Ey∼p[‖∇yk(y,x)‖2] = C1Ey∼p [∥∥∥∥x− yσ2 ∇zf(z) ∥∥∥∥ 2 ] , z = ‖x− y‖22 2σ2
(i) ≤ C2O(d−1) √ Ey∼p [ ‖x− y‖2 ] = O(d−1/2),
where (i) is by Cauchy-Schwarz and the Lipschitzity of f . Similarly, following the expansion of MSEp[S2] and the assumptions that k is upper-bounded, one can verify that MSEp[S2(y,x)] = O(1). Finally, we lower-bound the variance of S1 via the following calculation,
Ey∼p [ ‖S1(y,x)− µx‖22 ] = Ey∼p [ ‖∇y log p(y)‖22 k 2(x,y) ] − ‖µx‖22
≥ C1Ey∼p [ ‖y‖22 k 2(x,y) ] − ‖µx‖22
(i) ≥ C2Ey∼N (0,I) [ ‖y‖22 e − ‖x−y‖ 2 2 σ2 ] − ‖µx‖22 (ii) = Ω(d),
where (i) is by the strongly log-concavity of p and the fact that the kernel being lower-bounded by a scaled Gaussian RBF kernel, and (ii) directly follows from result in the Gaussian case. Combining the calculations yields the desired statement.

heading: C PROOF OF PROPOSITIONS IN SECTION 5
text: In this section we aim to calculate the equilibrium variance of SVGD and MMD-Descent under the proportional asymptotics. We first restate and comment on the Assumptions in Section 5.
• (A2) Near-orthogonality. Particles at fixed point of SVGD (or MMD-descent) {xi}ni=1 satisfy |x>i xi − dυ| < υ d, |x>i xj | < υ d for all i 6= j, some υ > 0, and d−1/2 d → 0 as d→∞ with probability 1.
• (A4) Gaussian Target. p(x) ∝ exp ( − 12x >x ) .
Under (A4) the goal of SVGD or MMD-descent is to draw samples from a unit Gaussian distribution, and we aim to calculate the dimension-averaged variance v of the particles evolved by the two updates. Note that since both SVGD and MMD-descent form an interacting particle system, one can no longer treat the converged particles as i.i.d Gaussian samples, i.e. xi ∼ N (0, vI). This significantly complicates the analysis in the proportional limit.
However, empirical results (see Figure 8) demonstrate that starting from non-degenerate random initialization, the equilibrium solution of both algorithms admits thin-shell concentration around the radius, i.e., the norm of xi concentrates around √ dv, and particles are almost orthogonal to one another. This observation is captured by (A2). We conjecture that such property holds for the updates we consider due to the repulsive force pushing the particles away from one another, and thus leading to this near-orthogonal configuration.
Under assumptions (A2)(A4), we are able to compute the stationary variance of both SVGD and MMD-descent in the asymptotic limit. We first provide a detailed derivation for the case of Gaussian RBF kernel, and then extend the analysis to more general settings.
C.1 GAUSSIAN KERNEL WITH THE MEDIAN BANDWIDTH
Stationary Variance of SVGD in Proportional Limit. We aim to solve the stationary point of SVGD update (2). Under (A2) and the median heuristic σ = √ Med{‖xi − xj‖22}/2, the fixed point condition for the Gaussian RBF kernel is given as
∆(xk) = 1
n n∑ i=1 [ − k(xi,xk)xi + 1 dv k(xi,xk)(xk − xi) ] = 0, (6)
for all k, or equivalently n∑ i=1 k(xi,xk)xi = 1 dv + 1 n∑ i=1 k(xi,xk) · xk. (7)
Note that for the left hand side of (7), we have the following equivalence,
LHS = n∑ i=1 k(xi,xk)xi = Xkk,
where X = [x1, · · · ,xn] ∈ Rd×n is the data matrix, kk = [k(x1,xk), · · · , k(xn,xk)]> ∈ Rn is the k-th column of kernel Gram matrix K ∈ Rn×n. As for the RHS of Equation (7), note that for i 6= k, (A2) allows us to take the following Taylor expansion on entries of the kernel matrix around its concentrated value, i.e.,
k(xi,xk) = exp
( − ‖xi − xk‖22
2dv
) = e−1 +O( ).
Where · d1/2 → 0. Similarly for i = k we have k(xk,xk) ≈ 1. Combining the equations above,
RHS = 1
dv + 1 n∑ i=1 k(xi,xk) · xk = ( n+ e− 1 (dv + 1)e +O( ) ) xk.
Equating the RHS and LHS of Eq (7) in matrix form (over all k) we have
X ·K = n+ e− 1 (dv + 1)e X +X · diag( ).
where diag( ) is a squared matrix where the i-th diagonal element is the error from Taylor expansion i = O( ). Define
m = n+ e− 1 dv + 1 .
The fixed point of SVGD thus simplifies to
X · (K −mIn − diag( )) = 0. (8)
Denote A = K − mIn − diag( ), Recall that the K is an Euclidean kernel matrix with Kij = k(xi,xk) = exp ( −(2dv)−1 ‖xi − xk‖22 ) . From Theorem 4 in Bordenave et al. (2013), it follows
that the empirical spectrum of A, which we write as µ(A) = n−1 ∑n i=1 δλi(A), converges weakly to the following quantity,
µ(A)→ (
1− 2 e −m
) + 1
e µ
( 1
dv X>X
) + µ(diag( )),
In addition, define S = n/(n− 1)In − 1/(n− 1)1n1>n , then by the Hoffman-Wielandt inequality (see Bordenave et al. (2013, Lemma 6)) we have
W2
( µ ( 1
dv X>X
) , µ(S) ) ≤ √ 1
n tr
( 1
dv X>X − S
)2 = √ n−1 · (nO( ))2 → 0,
where we used (A2), and W2(·, ·) is the 2-Wasserstein distance. Hence, we know that µ(A)→ (
1− 2 e −m
) + 1
e µ
( 1
dv X>X
) + µ(diag( ))
→ (
1− 2 e −m
) + 1
e µ (S) + µ(diag( )). (9)
When γ > 1 (i.e. d > n), Equation (8) requires µ(A)→ 0. Consequently,
1− 2 e + n e(n− 1) −m = 0 ⇔ m→ 1− e−1.
Therefore, from the definition of m we have the desired result.
vSVGD → n d(e− 1) = 1 e− 1 1 γ . (10)
Stationary Variance of MMD-descent in Proportional Limit. First note that for the Gaussian RBF kernel, the driving force of MMD-descent (4) admits the following closed-form,
Ey∼p[S2(y,x)] = − σd
(1 + σ2)d/2+1 exp
( − ‖x‖ 2 2
2 + 2σ2
) x.
This can be verified via simple numerical calculation:∫ p(y)k(x,y)∇y log p(y)dy
= ∫ e− ‖x−y‖22
2σ2 (−y) 1√ (2π)d e− y>y 2 dy
=− 1√ (2π)d ∫ y ye− ( 1√ 1+σ2 x− √ 1+σ2y )2 2σ2 e − x>x 2+2σ2 dy
=− σ d (1 + σ2)d/2+1 e − ‖x‖ 2 2 2+2σ2 x.
Therefore, the stationary point of MMD-descent satisfies
∆xk = − σd (1 + σ2)d/2+1 e −‖ xk‖22 2+2σ2 xk + 1 nσ2 ∑ i 6=k k(xk,xi)(xk − xi) = 0,
∀k, or equivalently, n∑ i=1 k(xk,xi)xk − ( dv 1 + dv )d/2+1 e− ‖xk‖22 2+2dv nxk = n∑ i=1 k(xk,xi)xi. Under assumptions (A2)(A3), similar to the SVGD case, we have the matrix form of the fixed point,[ 1 + e−1(n− 1)− ( dv
1 + dv
)d/2+1 e− dv 2+2dv n ] X +Xdiag( ) = XK.
with = o(1). Following an anlogous calculation, as n, d→∞ with d/n = γ ∈ (1,∞) we obtain, 1 + e−1(n− 1)− ( dv
1 + dv
)d/2+1 e− dv 2+2dv n → 1− 1
e . (11)
Moreover, observe that limd→∞ (dv/(1 + dv)) d/2+1 = e−1/(2v). We arrive at the desired result,
vMMD → 1.
Following the exact same reasoning, one can show that for fixed bandwidth σ = c √ d for some constant c ∈ Θ(1), MMD-descent also estimates the variance correctly, i.e., vMMD = 1.
SVGD and MMD-descent under finite n and large d. From (A2) we have the following decomposition for i 6= j,
‖xi − xj‖22 = ‖xi‖ 2 2 + ‖xj‖ 2 2 − 2x > i xj = 2dv + 2(n− 1)−1dv +O( ) =
2n
n− 1 dv +O(d ).
For finite n, this indicates that the Euclidean distance between two particles concentrates around 2(n − 1)−1ndv (instead of 2dv); this is to say, the Gaussian RBF kernel admits the the following Taylor expansion for i 6= j,
k(xi,xj) = exp
( − ‖xi − xk‖22
2dv
) = e− n n−1 +O( ).
Therefore, for the MMD-descent algorithm, (11) reduces to
1 + e− n n−1 (n− 1)−
( dv
1 + dv
)d/2+1 e− dv 2+2dv n→ 1− e− n n−1 .
Hence the converged particles of MMD-descent with finite n and infinite d has variance
vMMD → n− 1 n+ 1 .
And for SVGD, one can use the exact same argument to obtain that
vSVGD → 0. Note that this result agrees with our characterization in the proportional asymptotic limit by taking γ →∞ (i.e., d n).
C.2 GENERAL KERNEL WITH ADAPTIVE BANDWIDTH (MEDIAN HEURISTIC) Recall the definition of the (bandwidth-adjusted) Euclidean distance kernel: k(x,y) = f ( ‖x−y‖22
2σ2
) .
Under (A2) and the median heuristic σ = √
Med{‖xi − xj‖22}/2, we know that the kernel can be equivalently written as k(x,y) = f ( ‖x−y‖22
2dv
) . Following the same procedure, we have the
following fixed point condition,
∆(xk) = 1
n n∑ i=1
[ −f ( ‖xi − xk‖22
2dv
) xi − 1
dv f ′
( ‖xi − xk‖22
2dv
) (xk − xi) ] = 0, (12)
for all k, or equivalently n∑ i=1 [ f ( ‖xi − xk‖22 2dv ) − 1 dv f ′ ( ‖xi − xk‖22 2dv )] xi = − 1 dv n∑ i=1 f ′ ( ‖xi − xk‖22 2dv ) xk.
Similar to the previous calculation on the Gaussian RBF kernel, applying (A2) and the differentiability of f around 1, we Taylor-expand the kernel and obtain the following equivalence:
X · (K − 1 dv K ′) = −f(0) + (n− 1)f ′(1) dv ·X +X · diag( ).
where K is the Gram matrix of the kernel k, K ′ is the Gram matrix of its derivative with K ′ij = f ′ ( ‖xi−xj‖22
2dv
) . Apply Theorem 4 of Bordenave et al. (2013) to K and K ′, we have the following
relation under (A1)(A2),
f(0)− f(1)− 1 dv
(f ′(0)− f ′(1)) = −f(0) + (n− 1)f ′(1)
dv , (13)
which gives the expression of the equilibrium variance,
vSVGD → f ′(1)
f(1)− f(0) · 1 γ .
Heuristic derivation of the modified (damped) update. We provide a brief sketch of the modified algorithm in Section 5 for the median bandwidth. Recall that our proposed update introduces a damping term λf(0)xk in the driving force S1. To derive the optimal λ for unit Gaussian target, we incorporate this damping term into the fixed point Equation (12), which leads to the following modification of Equation (13):
f(0)− f(1)− 1 dv
(f ′(0)− f ′(1)) = (1− λ)f(0)− f(0) + (n− 1)f ′(1)
dv .
Recall the true target variance v = 1; this gives following choice of λ: λ = f(0)−1 · ( f(1)− γ−1f ′(1) ) . (14)
Since our goal is to “weaken” the driving force S1(xk,xk), we take the minimum between the derived value in (14) and the default λ = 1.
C.3 GENERAL KERNEL WITH FIXED BANDWIDTH
We now consider Euclidean distance kernel with invariant (fixed) bandwidth that scales with the dimensionality d, which we write as k(x,y) = f(‖x− y‖22 /2d). Following (12), we have the following equilibrium condition for the particles,
n∑ i=1 [ f ( ‖xi − xk‖22 /2d ) − 1 d f ′ ( ‖xi − xk‖22 /2d )] xi = − 1 d n∑ i=1 f ′ ( ‖xi − xk‖22 /2d ) xk.
Under Taylor expansion around v, the fixed point condition entails that the stationary variance satisfies,
f(0)− f(v) + 1 d (f ′(0)− f ′(v)) = −1 d (f(0) + (n− 1)f ′(v)).
⇐⇒ f ′(v)
f(v)− f(0) = γ. (15)
Note that given (A3) and monotone f ′, the numerator of the RHS of (15) takes a negative value that increases with v, whereas the denominator is also negative but decreases with v. This implies that when γ becomes larger, v needs to decay towards 0 in order to satisfy the equation. For the Gaussian RBF kernel with fixed bandwidth σ = √ d, the equation can be easily solved as,
vRBF = log ( 1 + 1
γ
) < 1,
which is a decreasing function of γ > 1. On the other hand, for the IMQ kernel with fixed bandwidth σ = √ d, standard calculation yields,
vIMQ = 1
6
( 4γ(γ + 3)(
3γ4 √ 3(8γ + 27) + 8γ6 + 36γ5 + 27γ4 )1/3
+
( 3γ4 √ 3(8γ + 27) + 8γ6 + 36γ5 + 27γ4 )1/3
γ2
) − 2
3 .
One can numerically verify that the value is less than 1 for γ > 1 and also non-increasing.
IMQ kernel with dimension-independent bandwidth. Finally, we note that in the context of kernel Stein discrepancy (KSD), the IMQ kernel is often employed without the dimensiondependent bandwidth (Gorham and Mackey, 2017). In this case we write the kernel as k(x,y) = f(‖x− y‖22 /2), which gives the following stationary condition,
n∑ i=1 [ f ( ‖xi − xk‖22 /2 ) − f ′ ( ‖xi − xk‖22 /2 )] xi = − n∑ i=1 f ′ ( ‖xi − xk‖22 /2 ) xk.
Taking Taylor expansion around dv gives f(0)− f(vd)− f ′(0) + f ′(vd) = −f(0)− (n− 1)f ′(vd).
Note that the LHS is Θ(1) by (A2); in order for the inequality to hold asymptotically, we need to have f ′(vd) = Θ(n−1). On the other hand, for the fixed-bandwidth IMQ kernel,
f(a) = (1 + a)−1/2; 2f ′(a) = −(1 + a)−3/2. This implies that as d increases, the stationary variance decays to 0 at a rate of v = Θ(d−1/3). Following the exact same procedure, one can show that the log-inverse kernel with fixed dimensionindependent bandwidth also asymptotically collapses the variance of SVGD particles to 0 when γ > 1; we omit the derivation.
C.4 “ALMOST-GAUSSIAN” TARGET
For general p(x) ∝ exp(−f(x)), the fixed point equation of xk is given as,
∆(xk) = 1
n n∑ i=1
[ −f ( ‖xi − xk‖22
2σ2
) ∇f(xi)− 1
σ2 f ′
( ‖xi − xk‖22
2σ2
) (xk − xi) ] = 0.
Assume that SVGD does not underestimate that marginal variance, then by (A2) and the quadratic growth, we know that there exists some α > 0 such that for every particle xk and coordinate m ∈ [d],[ n∑ i=1 αf ( ‖xi − xk‖22 2σ2 ) xi − 1 σ2 f ′ ( ‖xi − xk‖22 2σ2 ) xi ] m = − 1 σ2 [ n∑ i=1 f ′ ( ‖xi − xk‖22 2σ2 ) xk ] m . For kernels with the median bandwidth, followings the same simplification as in (13), we get
X · ( αK − 1 dv K ′ ) = −f(0) + (n− 1)f ′(1) dv ·X +X · diag( ).
Solving the inequality yields,
v = α γ · f
′(1)
f(1)− f(0) .
Therefore, given any dimension-independent growth of the potential α, there exists a large enough γ such that v < 1, i.e., SVGD underestimates the marginal variance. The case for fixed bandwidth follows from the same line of reasoning, the details of which we omit.

heading: D EXPERIMENT SETUP
text: Bayesian Neural Network. We consider a BNN with two hidden layers of 100 units. In each layer, the preactivations st+1 are computed via st+1 = (W tat + bt)/ √ ht + 1, where at ∈ Rht are the input activations. The target function f is a BNN with the same architecture, whose weights and biases are randomly generated from standard normal distributions. For the training set, we sampled 10 input locations uniformly from [−2.5,−1.5] and [1.5, 2.5], respectively. We add random noises to the observations, y = f(x) + , ∼ N (0, 0.01). We first adopt the Hamiltonian Monte Carlo (HMC) (Neal et al., 2011) to generate asymptotic posterior samples. HMC used 50 independent chains and each chain selected 50 particles with the frequency of 100 iterations after 5k burn-in iterations. HMC generates diverse particles as shown in Figure 1(a). We then simulate SVGD and MMD-descent dynamics on this problem. For both updates we use the Gaussian RBF kernel with the median bandwidth heuristic. For SVGD, we evolve 100 particles for 50k iterations using learning rate η = 5e − 3, whereas for MMD-descent, we evolve 10 particles using HMC particles as approximate target samples to compute the driving force term. The particles are either initialized from the (approximate) target distribution (as in Figure 1), or from standard normal distribution (as in Figure 9). Note that the plots do not include observation variance.
Hyperparameter Setting in Section 5. For all experiments in Section 5, we initialize the particles from N (0, 0.8Id), and run SVGD (or MMD-descent) with learning rate η = 10−1 for 20k iterations. For experiments in the proportional limit (Figure 4(a)(b)), we fix n = 50 and vary d. For the median heuristic (following Garreau et al. (2017)), we compute the median of Euclidean distance between all particles at each iteration.

Review summaries:
In this paper, the authors analyze the underestimation issue of stein variational gradient descent (SVGD), and propose the maximum mean discrepancy (MMD) descent. From the perspective of the decomposition of the gradient term (driving force and repulsive force), this paper suggested to use another driving force term instead of the original one in SVGD. But the new driving force makes MMD-descent impractical since it depends on an intractable integral of the desired distribution $p$. In addition, the paper identify the log derivative driving force as the problematic term in SVGD, and propose a modified SVGD with particle resampling. They also argue that the proportional asymptotic limit is more relevant to understanding the variance collapse phenomenon. The theoretical dimensional analysis of SVGD on Gaussian also suggested another modified (damped) SVGD.

This paper provides an understanding of the variance collapse phenomenon of SVGD. The paper first (1) introduces the reader to the most important concepts and phenomena, then (2) gives an explanation for why this problem occurs, thanks to a comparison with an accurate (yet computationally intensive) algorithm they call MMD-descent. Finally (3) the paper shows how to fix SVGD with damping. The paper provides experiments and theory, nicely combined.

This work studies the variance collapse phenomenon of SVGD. By comparing to MMD-descend, the authors argue that the driving force of SVGD suffers from a bias caused by reusing data, and thus tends to underestimate the variance of the target distribution. Theory are developed in the setting of estimating standard Gaussian with a proposal limit (i.e., $d/n \to \gamma$), and explains the understanding in the overparameterized/high-dim setting (i.e., $\gamma > 1$). Experiments are also conducted to verify the understanding. Finally, motivated by the understanding, new algorithm is proposed to fix the issue of SVGD by damping the driving force term in SVGD.

This paper analyzed the curse-of-dimensionality problem of the vanilla SVGD with Euclidean distance kernel in a qualitative and quantitative way. Specifically, the author first built a connection of SVGD to MMD-descent, where they share identical repulsive forces with different driving forces (if Euclidean distance kernel is adopted). Then, the author argued that the variance collapse problem is rooted in (1) high variance and (2) the deterministic bias of the driving force, which were confirmed by sampling from the isotropic Gaussian. Quantitatively, the author analyzed the stationary variance of MMD-descent and SVGD with isotropic Gaussian under the proportional limit, which confirms the curse-of-dimensionality problem of SVGD.