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.