Stein Variational Gradient Descent: Particles That Repel Their Way to the Posterior

Bayesian inference is simple to state and hard to compute. Given data, we want the posterior distribution over parameters. In almost every interesting model, that posterior has no closed form, its normalising constant is intractable, and direct sampling is expensive. Two classical families of tools attack the problem from opposite directions. Markov chain Monte Carlo (MCMC) gives asymptotically exact samples but can be slow to mix. Variational inference (VI) turns inference into optimisation, but usually forces the posterior into a chosen parametric family, such as a Gaussian, and then hopes that family is flexible enough.

Stein variational gradient descent (SVGD), introduced by Liu and Wang, offers a third way. It keeps a set of particles, moves them with a deterministic gradient-like rule, and lets the particles themselves represent the posterior. There is no parametric family to fit and no accept-or-reject step. The trick that makes it work is a repulsive force between particles. Without it, every particle would collapse to the same mode. With it, the particles spread out and trace the shape of the posterior.

The problem, stated precisely

Let $\theta$ be a parameter vector and let $p(\theta)$ be the target posterior, known only up to a constant. That is, we can evaluate $\log p(\theta)$ and its gradient, often called the score function, but we cannot evaluate the normalising constant and we cannot sample from $p$ directly.

Suppose we have $n$ particles $\theta_1, \dots, \theta_n$ drawn from some initial distribution $q$ . Our goal is to move these particles, step by step, so that their empirical distribution gets closer to $p$ . Closeness here is measured by Kullback-Leibler (KL) divergence:

$$\text{KL}(q \parallel p) = \mathbb{E}_{\theta \sim q}[\log q(\theta) - \log p(\theta)]$$

If we apply a smooth transformation $T(\theta) = \theta + \epsilon \phi(\theta)$ to every particle, with a small step size $\epsilon$ and a vector field $\phi$ to be chosen, the distribution of the transformed particles changes. The question is: which $\phi$ reduces the KL divergence fastest?

The Stein identity and the best direction

The answer comes from Stein's identity. For a smooth test function $\phi$ that vanishes suitably at infinity, the following expectation under $p$ is zero:

$$\mathbb{E}{\theta \sim p}[\nabla{\theta} \log p(\theta) \cdot \phi(\theta) + \nabla_{\theta} \cdot \phi(\theta)] = 0$$

This identity is the reason SVGD needs only the score function $\nabla_{\theta} \log p(\theta)$ , not the normalising constant. The constant disappears when we take the log gradient.

Stein variational gradient descent restricts $\phi$ to lie in a reproducing kernel Hilbert space (RKHS) defined by a positive definite kernel $k$ . In that space, the direction of steepest descent of the KL divergence has a closed form. It is the expectation, over the current particle distribution, of a kernel-weighted score plus a kernel gradient term:

$$\phi^*(\theta) = \mathbb{E}{\theta' \sim q}[k(\theta', \theta) \nabla{\theta'} \log p(\theta') + \nabla_{\theta'} k(\theta', \theta)]$$

That single equation contains the whole algorithm. Read its two terms separately, because they do opposite jobs.

The first term, $k(\theta', \theta) \nabla_{\theta'} \log p(\theta')$ , is an attraction term. It is a kernel-weighted average of the score at the other particles. It pulls a particle uphill, toward regions where the posterior density is high, with nearby particles contributing more weight than distant ones.

The second term, $\nabla_{\theta'} k(\theta', \theta)$ , is a repulsion term. For the common radial basis function (RBF) kernel, this term points away from neighbouring particles. It stops particles from piling up on top of each other at a single mode. Without it, SVGD would be little more than parallel gradient ascent, and every particle would converge to the nearest mode. With it, the particles behave like a system of mutually repelling charges that still feels the pull of the posterior.

In practice we replace the expectation by the average over the finite particle set. The update for particle $i$ at iteration $t$ is:

$$\theta_i \leftarrow \theta_i + \epsilon \, \hat{\phi}^(\theta_i), \quad \hat{\phi}^(\theta_i) = \frac{1}{n} \sum_{j=1}^{n} [k(\theta_j, \theta_i) \nabla_{\theta_j} \log p(\theta_j) + \nabla_{\theta_j} k(\theta_j, \theta_i)]$$

That is the entire algorithm. Initialise particles, repeatedly apply this update, and use the final particles as approximate posterior samples.

The kernel matters more than it looks

The usual choice is the RBF kernel:

$$k(\theta, \theta') = \exp\left(-\frac{|\theta - \theta'|^2}{2 h^2}\right)$$

The bandwidth $h$ controls the range of interaction. If $h$ is too small, each particle only sees itself, the repulsion term vanishes, and the particles collapse. If $h$ is too large, the kernel is nearly constant, the repulsion is weak and indiscriminate, and the particles move almost in lockstep.

The standard practical fix is the median heuristic: set $h$ to the median pairwise distance between particles, recomputed at each iteration. This adapts the interaction range to the current spread of the particles, which is exactly what you want as the particles move from a diffuse initialisation to a concentrated posterior. It is a heuristic, not a theorem, but it works well in a wide range of problems and it removes one fragile tuning knob.

For the RBF kernel, the gradient term in the update has an explicit form:

$$\nabla_{\theta_j} k(\theta_j, \theta_i) = -\frac{(\theta_j - \theta_i)}{h^2} k(\theta_j, \theta_i)$$

Notice the sign. When $\theta_j$ is close to $\theta_i$ , this term pushes $\theta_i$ away from $\theta_j$ , with a strength that decays with distance. That is the repulsion, written out in full.

A concrete example in code

Consider a simple target: a mixture of two Gaussians in two dimensions, with modes that are well separated. A single Gaussian variational approximation would be forced to pick one mode or sit uselessly between them. SVGD particles, by contrast, should split across both modes. The code below implements the update in NumPy for this toy target. It is short on purpose: the algorithm really is this short once the score function is available.

import numpy as np

def log_target_and_score(x):
    # Mixture of two 2D Gaussians, equal weights.
    # Returns log density (up to a constant) and its gradient.
    mu1 = np.array([2.0, 2.0])
    mu2 = np.array([-2.0, -2.0])
    d1 = x - mu1
    d2 = x - mu2
    log_p1 = -0.5 * np.sum(d1 ** 2, axis=-1)
    log_p2 = -0.5 * np.sum(d2 ** 2, axis=-1)
    m = np.maximum(log_p1, log_p2)
    log_p = m + np.log(0.5 * np.exp(log_p1 - m) + 0.5 * np.exp(log_p2 - m))
    w1 = np.exp(log_p1 - log_p)
    w2 = np.exp(log_p2 - log_p)
    score = -(w1[:, None] * d1 + w2[:, None] * d2)
    return log_p, score

def svgd(x, n_iter=500, eps=0.1):
    n = x.shape[0]
    for _ in range(n_iter):
        log_p, score = log_target_and_score(x)
        # Pairwise squared distances and the median heuristic.
        diff = x[:, None, :] - x[None, :, :]
        dist2 = np.sum(diff ** 2, axis=-1)
        h = np.median(np.sqrt(dist2 + 1e-8))
        h2 = h ** 2 + 1e-8
        K = np.exp(-dist2 / (2.0 * h2))
        # Kernel gradient term: -K * diff / h^2, summed over j.
        grad_K = -K[:, :, None] * diff / h2
        phi = (K @ score + np.sum(grad_K, axis=1)) / n
        x = x + eps * phi
    return x

rng = np.random.default_rng(0)
particles = rng.normal(size=(100, 2))
particles = svgd(particles)

Run this and plot the particles and you will see the intended behaviour: particles start as a single Gaussian blob, then split and settle around both modes, with their density roughly tracking the target. No particle is wasted sitting in the valley between the modes, because the score pulls it toward one mode or the other, and repulsion keeps the two groups from collapsing into points.

For a real Bayesian model, the only part that changes is log_target_and_score. For Bayesian logistic regression, for example, the score is the gradient of the log prior plus the log likelihood, which is available in closed form or by automatic differentiation. SVGD then gives you a particle approximation to the posterior over weights, from which you can compute predictive means, predictive variances, and credible intervals by averaging over particles.

Why SVGD is not just parallel MCMC or plain VI

It helps to place SVGD precisely relative to its neighbours.

Compared with MCMC, SVGD is deterministic given the initial particles. There is no random proposal and no rejection, so it can be much more sample efficient in the early phase, and it parallelises trivially across particles. The price is the loss of asymptotic exactness guarantees in the finite particle regime: with a fixed number of particles, SVGD converges to an approximation whose quality depends on the kernel, the bandwidth, and the particle count. MCMC, given infinite time, is exact. SVGD, given infinite particles, is also well understood in the mean-field limit, but you never have infinite particles.

Compared with parametric VI, SVGD is nonparametric in the sense that matters. The particle set can represent multiple modes, skewness, and other shapes that a Gaussian family cannot. The price is computational cost: each iteration costs $O(n^2 d)$ for $n$ particles in $d$ dimensions, because every particle interacts with every other particle through the kernel. For a few hundred particles this is fine. For tens of thousands, it is not, and you need approximations such as mini-batch subsampling of particles or a structured kernel.

There is also a deeper connection worth knowing. In the limit of many particles, SVGD can be viewed as a gradient flow in the space of distributions, with the kernel defining the geometry of that flow. This is the same style of thinking that underlies Wasserstein gradient flows, except that the kernel geometry replaces the Wasserstein geometry. That viewpoint explains why the choice of kernel is not a minor implementation detail: it defines what "steepest descent" means in distribution space.

Limitations and common pitfalls

First, SVGD degrades in high dimensions. The RBF kernel suffers from the curse of dimensionality: in high dimensions, pairwise distances concentrate, the median heuristic becomes less informative, and the repulsion term weakens. Practical remedies include using a kernel on a lower dimensional representation, using an anisotropic kernel with a learned length scale per dimension, or combining SVGD with a dimension reduction step. Do not expect out-of-the-box SVGD with a single isotropic RBF kernel to work well on a thousand dimensional posterior without such adjustments.

Second, too few particles is a silent failure mode. With, say, 20 particles on a multimodal target, SVGD may still look plausible: the particles spread out and find high density regions. But the particle distribution is a coarse approximation, and estimates of posterior variance or tail probabilities will be unreliable. The diagnostic is simple: increase the particle count and check whether your quantity of interest changes. If it moves materially, you did not have enough particles.

Third, the step size needs care. A fixed step size that is too large causes oscillation or divergence, especially early on when gradients are large. In practice, practitioners often use an adaptive optimiser such as AdaGrad or Adam on the SVGD update direction, treating $\hat{\phi}^*$ as a gradient. This is a heuristic layered on top of the theory, but it stabilises training considerably. Monitor the movement of particles, not just the objective, because SVGD has no simple scalar loss that tells you it has converged.

Fourth, SVGD needs a differentiable log density. If your model has discrete parameters, or if the score is unavailable or noisy, the basic algorithm does not apply directly. Variants exist for discrete and gradient-free settings, but they change the theory, so do not force the standard update onto a problem where the score does not exist.

Fifth, do not confuse SVGD particles with independent posterior samples. They are a deterministic, interacting approximation. For predictive averaging this distinction rarely matters much, but for estimating Monte Carlo error bars it does. If you need honest uncertainty about your Monte Carlo estimate itself, SVGD particles alone do not give it to you.

When to reach for SVGD

SVGD earns its place when three conditions hold: you can compute the score function, the dimension is moderate, and the posterior may be multimodal or otherwise non-Gaussian in a way that matters for your decision. Bayesian neural networks with a modest number of parameters, Bayesian logistic regression with structured priors, and probabilistic models where predictive uncertainty drives an action are all natural fits.

It is less compelling when you need asymptotically exact samples and can afford a long MCMC run, when the model is high dimensional and a well tuned parametric VI approximation is already good enough, or when the score function is unavailable. As with most inference methods, the right question is not whether SVGD is universally better, but whether its particular trade-off, flexibility for cost, matches your problem.

The mental image to keep is the one from the update equation: each particle is pulled toward high posterior density and pushed away from its neighbours. Attraction alone gives you optimisation. Repulsion alone gives you diversity. Together, and only together, they give you inference.

---

Connect with me on LinkedIn: linkedin.com/in/umeshtharukaofficial