Score Matching

The normalizing constant is the enemy

Write an unnormalized model — an energy function E_θ(x) computed by any network you like. The density it induces is

p_θ(x) = exp(−E_θ(x)) / Z(θ)      Z(θ) = ∫ exp(−E_θ(x)) dx

Maximum likelihood needs log p_θ(x) = −E_θ(x) − log Z(θ), and Z(θ) is an integral over all of R^d. For an image at d = 3072 that integral has no closed form and no cheap estimate; contrastive divergence approximates its gradient with an inner MCMC loop that is slow and unstable. This is the structural reason energy-based models were hard for decades. Score matching sidesteps it entirely. Because Z(θ) does not depend on x, differentiating with respect to x annihilates it: ∇_x log p_θ(x) = −∇_x E_θ(x). The obstacle vanishes by construction, not by approximation.

What the score function actually is

The score is a vector field, s(x) = ∇_x log p(x), with the same shape as the data: if x: [d], then s(x): [d]. At every point it says which way is uphill in log-density, and how steeply. For p = N(μ, σ^2 I) it is available in closed form:

log p(x) = −||x − μ||^2 / (2σ^2) + const     ⇒   s(x) = −(x − μ) / σ^2

A vector pointing back at the mean, scaled by distance — a restoring force. Two clarifications. First, this is not the statistician’s Fisher score ∇_θ log p_θ(x), which differentiates with respect to parameters; the generative-modelling score differentiates with respect to the data. Second, the score is strictly less information than the density: it determines p only up to a constant, which is precisely the constant we were trying to escape.

Fisher divergence and the integration-by-parts trick

The natural objective compares vector fields directly — the Fisher divergence:

J(θ) = ½ E_{p(x)} || s_θ(x) − ∇_x log p(x) ||^2

Useless as written: it needs the true score we do not have. Hyvärinen (2005) showed the cross term can be rewritten. Expand the square; the only problematic piece is −E_p[ s_θ(x) · ∇_x log p(x) ]. Substituting ∇log p = ∇p / p cancels the density in the expectation, leaving −∫ s_θ(x) · ∇p(x) dx, and integrating by parts flips the derivative onto s_θ:

J(θ) = E_{p(x)} [ tr( ∇_x s_θ(x) ) + ½ || s_θ(x) ||^2 ] + const

Everything left is computable from samples. The boundary term must vanish, which requires p(x) s_θ(x) → 0 as ||x|| → ∞ — a mild but real condition.

The trace is the bottleneck

That identity is exact and it is still impractical. ∇_x s_θ(x) is the [d, d] Jacobian of the network’s output with respect to its input, and reverse-mode autodiff gives you one row per backward pass. Extracting all d diagonal entries therefore costs d backward passes per sample: for a 32×32 RGB image, 3072 backward passes for one training example. That is a non-starter on a GPU cluster and unthinkable on a CPU budget.

Sliced score matching (Song et al., 2019) attacks this with Hutchinson’s trace estimator: for random v with E[vv^T] = I, tr(A) = E_v[ v^T A v ]. The quantity v^T ∇_x s_θ v needs a single Jacobian-vector product, so cost drops from O(d) passes to O(1) at the price of estimator variance. Better — but the practical answer turned out to be a different objective entirely.

Denoising score matching: Vincent’s identity

Vincent (2011) changed the target. Perturb each data point with a Gaussian kernel, x˜ = x + σε with ε ∼ N(0, I), and match the score of the conditional instead. That conditional is a Gaussian centred at x, so its score is known exactly:

∇_x˜ log q(x˜ | x) = −(x˜ − x) / σ^2 = −ε / σ

L(θ) = ½ E_{x, ε} || s_θ(x˜) + ε/σ ||^2

The identity says this equals Fisher divergence against the smoothed density q_σ, up to a constant. No Jacobian, no trace, no boundary condition — a plain regression whose target is a known vector. The catch is stated in the name of the density: you learn ∇ log q_σ, not ∇ log p, and the bias disappears only as σ → 0, exactly where the target −ε/σ blows up. That tension drives everything below.

A worked example, and why it is the diffusion loss

Take p = N(0, 1), σ = 0.5. Draw x = 1.2, ε = −0.8, so x˜ = 0.8. The regression target is −ε/σ = 1.6. But the true smoothed density is q_σ = N(0, 1.25), whose score at 0.8 is −0.8/1.25 = −0.64. The per-sample target is off by 2.24 — and that is fine, because averaging over every (x, ε) pair that lands on x˜ = 0.8 recovers it: the posterior mean of x is 0.64, giving −(0.8 − 0.64)/0.25 = −0.64 exactly. High variance, zero bias.

Now compare DDPM, where x_t = √α-bar_t · x_0 + √(1−α-bar_t) · ε:

s_θ(x_t, t) = − ε_θ(x_t, t) / √(1 − α-bar_t)

DDPM’s E||ε − ε_θ||^2 is denoising score matching with weight λ(t) = 1 − α-bar_t. Noise prediction and score prediction are the same network, rescaled.

Sampling: Langevin dynamics

A score is not obviously a sampler, but it is enough. Langevin dynamics walks uphill in log-density while injecting noise:

x_{k+1} = x_k + (η/2) · s_θ(x_k) + √η · z_k ,   z_k ∼ N(0, I)

As η → 0 and K → ∞ the chain’s law converges to p. Drop the noise and you get gradient ascent, which collapses onto modes; the √η term is what makes it sample rather than optimize. Note what never appears: Z(θ). The sampler consumes precisely the object score matching can learn.

Two honest caveats. Finite η leaves discretization bias unless you add a Metropolis accept/reject step. And mixing between modes separated by low-density valleys is exponentially slow — worse, if two modes are fully separated the score inside each is identical regardless of their relative weights, so the chain cannot learn how often to visit them.

Why one noise level is never enough

Naive score matching on clean data fails for two linked reasons. The manifold hypothesis: real images occupy a thin low-dimensional sheet in R^d, so p is zero almost everywhere and log p — hence its gradient — is undefined off the sheet. And coverage: the loss is an expectation under p, so regions of low density receive almost no gradient signal — yet those are exactly the regions a Langevin chain initialized from noise starts in.

Song and Ermon (2019) fixed both with a noise-conditional score network s_θ(x, σ), trained on a geometric ladder σ_1 > σ_2 > … > σ_L — say 50 down to 0.01. Large σ smears mass everywhere, curing both problems; small σ preserves detail. Weight each level by λ(σ) = σ^2 so the terms have comparable magnitude, since ||s|| ∼ √d/σ. Then run annealed Langevin: sample at σ_1, use the result to initialize σ_2, and walk down the ladder.

The bridge to flow matching

Score-based models and flow-based models are two coordinate systems on one object. Every diffusion SDE has a deterministic probability-flow ODE with identical marginals, whose drift is built from the score:

dx/dt = f(x, t) − ½ g(t)^2 · s_θ(x, t)

Once you are in ODE form you are doing flow matching — regressing a velocity field — and the velocity and the score are related by an affine change of variables fixed by the noise schedule. The companion article on rectified flow takes that road: linear interpolants, straight paths, and reflow to cut the step count. The practical difference is where the effort goes. Score matching gives you an SDE sampler that trades many steps for stochastic correction; flow parameterizations aim to make the trajectory straight enough that a handful of Euler steps — the deciding factor on a CPU, where each step is a full forward pass — suffice.