Diffusion Models as Generative Particle Dynamics

Diffusion models power most text-to-image systems today (Stable Diffusion¹, Imagen²). At a high level the recipe is simple: start from noise, apply a learned update many times, and you get an image. Under the hood, this is a transport process of “measure”/”mass”.

In an earlier post I covered the SDE formalism. Here I want to build geometric intuition instead: denoising defines a vector field, and generation describes particles following that field.

Brief history. Sohl-Dickstein et al.³ proposed learning the reverse of a noising process. Ho et al.⁴ (DDPM) simplified the loss to noise prediction and matched GANs on image quality. Song & Ermon⁵ and Song et al.⁶ connected the whole thing to score functions and continuous-time SDEs. That unified view opened ODE-based samplers (e.g., Flow Matching⁷⁸⁹) and is the backbone of current systems.

In this post I want to answer a simple question: why does repeatedly applying a small denoising correction produce good samples?

Start from denoising

Let’s start from the denoising problem itself. Take a clean sample $x$ and corrupt it:

A denoiser trained on this task learns more than a point estimate of $x$. Tweedie’s identity¹⁰ says:

Rearranging the identity gives:

The denoising residual is the score (up to $\sigma^2$). So if a denoiser works across noise levels, it is implicitly parameterizing the score field $\nabla_x \log p_\sigma(x)$ over the space¹¹. Milanfar & Delbracio¹² give a useful operator-theoretic interpretation: denoisers behave like gradient fields of smoothed energies.

Now think in particles

Once you view denoising as a score field, the sampling step becomes very concrete. You initialize $N$ particles $x_1, \dots, x_N$ at noise and run Langevin dynamics:

where $s_\theta(x,t) \approx \nabla_x \log p_t(x)$ is the learned score at noise level $t$. The noise term keeps exploration alive and helps avoiding mode and point collapse.

Sampling code (click to expand):

import numpy as np

def langevin_step_euler(x, score_fn, eta):
    """One Euler-Maruyama Langevin step.

    Args:
        x: Particle locations with shape (N, D).
        score_fn: Callable that maps x -> score with shape (N, D).
        eta: Step size.
    """
    grad = score_fn(x)
    noise = np.random.randn(*x.shape)
    return x + eta * grad + np.sqrt(2.0 * eta) * noise

def annealed_langevin(initial_x, score_fn, sigmas, steps_per_sigma=100, eta_factor=0.5, seed=42):
    """Run Langevin updates across a noise schedule.

    Uses eta = eta_factor * sigma**2 at each noise level.
    """
    np.random.seed(seed)
    x = initial_x.copy()
    for sigma in sigmas:
        eta = eta_factor * (sigma ** 2)
        for _ in range(steps_per_sigma):
            x = langevin_step_euler(x, lambda y: score_fn(y, sigma), eta)
    return x

def dummy_score(x, sigma):
    return -x / (1.0 + sigma ** 2)

np.random.seed(42)
particles = np.random.randn(500, 2) * 10.0
sigmas = np.geomspace(10.0, 0.05, num=10)
samples = annealed_langevin(particles, dummy_score, sigmas, steps_per_sigma=100)

Each step has a clear role: the score term moves particles toward higher density, and the noise term keeps exploration alive so the cloud does not collapse too early. Over many steps, the particles travels from the initial noise distribution toward $p_{\mathrm{data}}$.

Franceschi et al.¹³ call this viewpoint “generative particle dynamics.” In that language, the model defines a velocity field on sample space, and generation is just integrating that field. Whether you use Langevin noise, an ODE solver, or something in between is a discretization choice. The underlying object is the same: transport under a learned field.

A toy example: analytic score on a mixture

Before training a score network, it helps to isolate the transport story with an exact score. A Gaussian mixture is perfect for that. Take a $K$-component mixture of Gaussians (MoG) $p(x)=\sum_k \pi_k\,\mathcal{N}(x;\mu_k,\Sigma_k)$.

Adding noise $\mathcal{N}(0,\sigma^2 I)$ just adds up on each component’s covariance¹⁴:

Differentiating $\log p_\sigma$ gives the score in closed form:

We can read it as: the score at $x$ is a weighted average of pulls toward each mode mean $\mu_k$, weighted by how likely component $k$ explains $x$ (the posterior responsibility $w_k$), and scaled by that component’s precision $S_k^{-1}$.

At large $\sigma$ the pulls are weak and roughly uniform ($S_k\approx\sigma^2 I$, so precision is tiny). As $\sigma$ shrinks the pulls sharpen into separate basins around modes. Particles follow this field, and injected noise helps them jump across low-density barriers between basins.

In this MoG example, removing the noise term makes particles settle near mode centers and under-sample the lower-density regions between modes.

I use a 6-component mixture: six modes at irregular positions in $[-10, 10]^2$, non-uniform weights ($\pi_0\approx 0.55$ dominates), and anisotropic variances. No symmetry to exploit.

Because the score is available in closed form, there is no network training in this toy setup. That is intentional: we can study transport dynamics directly, without approximation error from a learned model.

The score field at three noise levels:

With that setup in place, we can run annealed Langevin dynamics: a geometric schedule from $\sigma=10$ down to $\sigma=0.05$, with 100 Langevin steps per level and step size $\eta = 0.5 \cdot \sigma^2$.

As a sanity check, compare the KDE of final Langevin samples to true samples drawn directly from the mixture:

Adding batch interactions

So far, each particle evolves independently. Each $x_i$ only sees the external score $\nabla_x \log p_\sigma(x_i)$, and the rest of the batch does not enter its update.

A natural next step is to let nearby particles interact through a kernel. This is exactly the SVGD idea¹⁵. The deterministic update is

and then moves particles by

You can read this update as a balance between attraction and diversity. The first term is a kernel-smoothed score that pulls particles toward high-density regions. The second term is the kernel-gradient interaction. Pairwise, its direction depends on which variable you differentiate with respect to, so it is better to think of it as a batch-level spreading effect rather than a literal point-to-point repulsive force. If we also add Langevin noise, we get a stochastic variant:

With a radial basis function (RBF) kernel $k_h$ at bandwidth $h>0$, nearby pairs interact strongly and far-away pairs barely interact. In aggregate over the batch, this term counteracts particle collapse and improves coverage. Without this interaction you get something closer to mean shift. With it, you get SVGD plus optional Langevin noise.

For the kernel used in the code,

the gradient with respect to the first argument is

The batch effect is easier to see by summing over neighbors:

So the interaction pulls each particle away from its kernel-weighted local mean $\bar{x}_{k,i}$ when viewed through the full batch update.

So the interaction is just the original kernel value scaled by a displacement and $1/h^2$. Larger bandwidth weakens this spreading effect everywhere, and the exponential still ignores interactions for distant pairs.

Interacting sampling code (click to expand):

import numpy as np

def rbf_kernel(x, bandwidth):
  """Compute an RBF kernel matrix and pairwise offsets.

  Args:
    x: Particle locations with shape (N, D).
    bandwidth: RBF bandwidth h.

  Returns:
    kernel: Shape (N, N), where kernel[i, j] = k(x_i, x_j).
    diff: Shape (N, N, D), where diff[i, j] = x_i - x_j.
  """
  diff = x[:, None, :] - x[None, :, :]
  sqdist = np.sum(diff ** 2, axis=-1)
  kernel = np.exp(-sqdist / (bandwidth ** 2))
  return kernel, diff


def stein_interacting_step(x, score_fn, sigma, eta, bandwidth, repulsion=1.0):
  """One noisy SVGD-style update.

  Args:
    x: Particle locations with shape (N, D).
    score_fn: Callable score_fn(x, sigma) -> array with shape (N, D).
    sigma: Current noise level.
    eta: Step size.
    bandwidth: RBF bandwidth h.
    repulsion: Weight on the kernel-gradient interaction term.
  """
  score = score_fn(x, sigma)
  kernel, diff = rbf_kernel(x, bandwidth)
  n_particles = x.shape[0]

  pooled_score = (kernel @ score) / n_particles
  repulsive_term = (2.0 / (bandwidth ** 2)) * np.sum(kernel[:, :, None] * diff, axis=1)
  repulsive_term /= n_particles

  drift = pooled_score + repulsion * repulsive_term
  noise = np.random.randn(*x.shape)
  return x + eta * drift + np.sqrt(2.0 * eta) * noise


def annealed_interacting_sampler(
  initial_x,
  score_fn,
  sigmas,
  steps_per_sigma=100,
  eta_factor=0.5,
  bandwidth_factor=1.0,
  repulsion=1.0,
  seed=42,
):
  """Run annealed interacting sampling with an SVGD-style drift.

  Uses h = bandwidth_factor * sigma and eta = eta_factor * sigma**2.
  """
  np.random.seed(seed)
  x = initial_x.copy()
  for sigma in sigmas:
    eta = eta_factor * (sigma ** 2)
    bandwidth = bandwidth_factor * sigma
    for _ in range(steps_per_sigma):
      x = stein_interacting_step(
        x,
        score_fn=score_fn,
        sigma=sigma,
        eta=eta,
        bandwidth=bandwidth,
        repulsion=repulsion,
      )
  return x

Are we double-counting with score and kernel gradients?

At first glance, it can feel like double-counting: both the MoG score and the RBF kernel gradient look like an exponential weight times a displacement vector. But they play different roles.

The score term is attraction to high-density regions of the target:

It points toward where probability mass lives.

The kernel-gradient term is interaction between particles:

For fixed $(i,j)$, this vector points from $x_j$ toward $x_i$ because the derivative is taken with respect to $x_j$. But in the summed SVGD update for each particle, the net effect is to discourage particle collapse and improve coverage across modes.

You can think of it this way: the score asks, “where should this particle move if it were alone?” The kernel term asks, “given where everyone else already is, how should this finite set rearrange to cover the target better?”

So the two terms are not redundant. One encodes the target density, and the other controls how a finite cloud of particles covers that density. If you remove the score term, particles only spread and never lock onto the right regions. If you remove the interaction term, particles follow the score independently and can pile into the dominant basin.

This also clarifies the Langevin connection. Classical Langevin is

Vanilla SVGD is deterministic (no injected Gaussian noise): it replaces the single-particle drift with a kernel-coupled batch drift. So SVGD is not “Langevin where noise prevents collapse”. In SVGD, anti-collapse comes primarily from the deterministic kernel-gradient interaction. If you add Gaussian noise to SVGD, you get a stochastic SVGD/Langevin-SVGD hybrid.

These distinctions become concrete when we run all three settings side by side on the same imbalanced MoG:

(A) collapses almost entirely into the dominant basin — score-only drift has no mechanism to maintain diversity. (B) gives good coverage out of the box, reaching TV $= 0.0037$ to the target mode weights. (C) can also reach TV $= 0.037$, but only with a broader kernel and stronger repulsion. With weak repulsion it still collapses, because the kernel interaction is local and cannot move particles back across basins once mass has pooled. Noise and repulsion both encourage exploration, but differently: noise is stochastic and global, while the kernel interaction is deterministic and local to nearby particles.

Why this framing helps

The particle picture strips diffusion models down to mechanics: the score is a vector field, sampling is integration, and the noise schedule controls how far particles can jump between basins. You do not need the full SDE formalism to build good intuition for what happens during generation.

The practical recipe is still short: corrupt data at many noise levels, train a denoiser (which is a score estimator in disguise), initialize particles from noise, and iterate the learned update. Independent Langevin is the simplest version. If you want particles to coordinate within a batch, kernel-coupled updates (SVGD-style¹⁵) give a natural extension, at the cost of an extra bandwidth parameter.

References