If you understand the score — the gradient of log p(x) — you understand the engine behind DDPM, Stable Diffusion, FLUX, Sora, and every score-SDE descendant. Yang Song's 2021 SDE paper is arguably the single most-cited foundation of modern generative AI. Score matching is the math that powers a $30B/year image-generation industry.
Learning Objectives
After this lesson, you will be able to:
Define the score function ∇_x log p(x) and explain why knowing the score is enough to sample from the distribution — even when you can't compute p(x) directly
Use Langevin dynamics to walk noisy gradient steps toward high-density regions, and recognize why a single noise level fails on real-world data
Train a noise-conditioned score network with denoising score matching, then sample with annealed Langevin from high noise down to low noise
Connect score-based models to diffusion (DDPM) and flow matching through the SDE framework — they are different parameterizations of the same underlying idea
Build this --> Implement Langevin sampling on a 2D mixture of Gaussians, then train a tiny score network with denoising score matching and watch your samples converge toward the modes; you'll see annealing in action when you compare single-σ vs multi-σ training
Don't worry if "score function" sounds intimidating — it is just the gradient of the log-density. Once you see Langevin dynamics turn a noisy gradient into samples, the whole framework clicks.
Drag a point around the density and watch the score arrows always point toward where the data lives.
Loading visualization...
The score function is defined as:
sθ(x)=∇xlogp(x)
Why the gradient and not the density? Computing p(x) directly requires the partition function Z = ∫ exp(-E(x)) dx — the normalizing constant that makes p(x) integrate to 1. In high dimensions, Z is intractable. But when you take the gradient of log p(x), Z disappears (it is a constant in x). The score is partition-function-free.
If you know the score, how do you turn it into samples? The answer comes from physics: Langevin dynamics, originally a model of molecular motion in fluids. The discrete update rule is:
xt+1=xt+2η∇xlogp(xt)+ηzt,zt∼N(0,I)
This is beautiful in theory, but on real data it fails for one critical reason: the score is unreliable in low-density regions. Most of high-dimensional space is empty (the manifold hypothesis). In those empty regions, the training data gives the network almost no signal about where to go. Langevin starting from random noise wanders aimlessly because the score it sees is junk.
The fix: add noise to the data and train the score network on noisy versions across many noise levels.
Vincent (2011) showed an elegant trick: if you take a clean data point x, perturb it with Gaussian noise to get x̃ = x + σz, then the score of the perturbed distribution has a closed form:
∇x~logpσ(x~∣x)=−σ2x~−x
The denoising score matching loss is then:
LDSM(θ)=Ex,σ,ϵ[λ(σ)sθ(x+σϵ,σ)−σ−ϵ2]
The crucial detail: train on multiple noise levels σ, from large (σ_max ≈ 50) down to tiny (σ_min ≈ 0.01). This is the noise-conditional score network (NCSN) of Song & Ermon 2019. The network learns the score of progressively noisier versions of the data, giving it dense supervision across the entire input space — even in regions far from the data manifold, where high noise reaches.
Now to sample. Start from pure noise. Run Langevin dynamics conditioned on σ_max for a few steps. Then drop to a smaller σ. Run more Langevin steps. Keep dropping σ until you reach σ_min. By the end, the samples have moved from "noise" through "blurry data" to "sharp data."
What Do You Think?
You train a score network with denoising score matching but only on a single noise level σ=0.1. You then try to sample with annealed Langevin from σ=10 down to σ=0.01. What goes wrong?
The right answer is the second option. The whole point of NCSN is that the score network needs to know what to do at every noise level the sampler will encounter. If you trained only at σ=0.1, the σ=10 forward pass is wildly out-of-distribution; the score outputs are garbage and Langevin walks off into nowhere. You must train across the same noise schedule the sampler uses.
Annealed Langevin: for σi∈[σmax,…,σmin]:ηi=ϵ⋅(σi/σmin)2
Yang Song's 2021 paper made the leap: instead of discrete noise levels, treat the noise schedule as a continuous time variable t ∈ [0, T] and model the corruption as a stochastic differential equation (SDE) — see the math-track lesson on stochastic calculus and SDEsStochastic CalculusStochastic calculus extends derivatives to random processes — Brownian motion, Itô integrals, and stochastic differential equations. The math diffusion models, flow matching, and score-based generative models all live in.Learn more → for the Brownian-motion and Itô machinery this rests on:
This is the punchline: DDPM, NCSN, and continuous-time diffusion are all the same algorithm in different parameterizations. DDPM predicts the noise ε; the score view predicts ∇log p; they differ only by a scaling factor σ. Once you grasp this, every diffusion paper becomes readable.
The SDE view also unlocks the probability flow ODE — a deterministic counterpart to the reverse SDE that gives the same marginal distributions but is easier to integrate with adaptive solvers. This is what enables fast samplers like DDIM and DPM-Solver that can produce high-quality samples in 10-50 steps instead of 1000.
The previous section presented the forward SDE, the reverse SDE, and the probability-flow ODE as a fait accompli. That is not a derivation — it is an assertion. The actual mathematics is one of the most beautiful results in stochastic calculus, and the rest of score-based generative modeling is just engineering on top of it. The forward process is defined (we choose how to destroy the data); the reverse process is derived (it falls out of the Fokker–Planck equation and a 1982 paper by Brian Anderson). The math-track lesson on stochastic calculus and SDEsStochastic CalculusStochastic calculus extends derivatives to random processes — Brownian motion, Itô integrals, and stochastic differential equations. The math diffusion models, flow matching, and score-based generative models all live in.Learn more → covers Brownian motion and the Itô formula at the level we need; this section assumes that background and builds on it.
A stochastic differential equation in our setting has the general form:
dx=f(x,t)dt+g(t)dw,t∈[0,T]
Score-based literature uses two canonical choices for f and g:
Variance-Exploding (VE) SDE. Set f(x, t) = 0 and g(t) = σ_max · (σ_min / σ_max)^t · √(2 log(σ_max / σ_min)). There is no drift, so the SDE is pure noise injection. The marginal variance grows unboundedly with t — concretely the marginal at time t is N(x_0, σ²(t) I) with σ²(t) = σ_min² · (σ_max / σ_min)^(2t). This is the continuous-time limit of the NCSN noise schedule.
Variance-Preserving (VP) SDE / DDPM. Set f(x, t) = -½ β(t) x and g(t) = √β(t) for some noise schedule β(t) > 0. The drift pulls x toward the origin while the diffusion injects noise; the two balance so that the marginal variance stays bounded. The Euler–Maruyama discretization of this SDE is exactly DDPM (Ho et al. 2020). If you sample t = 1, 2, …, N instead of continuous t, you recover the discrete forward chain x_n = √(1 - β_n) x_{n-1} + √β_n ε.
What Do You Think?
You train two diffusion models — one on the VE SDE (no drift, exploding variance), one on the VP SDE (drift pulls toward zero, bounded variance). You inspect the marginal p_T at the final time T. What do you see?
Every SDE induces a partial differential equation that governs how its marginal density p_t(x) evolves over time. For the forward SDE above, the marginal satisfies the Fokker–Planck equation (also called the Kolmogorov forward equation):
∂t∂pt(x)=−∇x⋅[f(x,t)pt(x)]+21g(t)2∇x2pt(x)
The math-track stochastic-calculus lesson derives Fokker–Planck from Itô's formula applied to a test function φ(x): take expectation of dφ(x_t), integrate by parts, and identify the result with ∫φ · ∂p/∂t · dx. We state the result here and move on. The crucial property: the Fokker–Planck equation is time-symmetric in form — its structure does not privilege the forward direction of time. That symmetry is the door Anderson walked through.
In 1982, Brian D. O. Anderson ("Reverse-time diffusion equation models", Stochastic Processes and their Applications 12) proved a remarkable theorem: any diffusion process with smooth marginals has a time-reversed counterpart that is also a diffusion process, and the reversed drift is determined by the score of the marginal density. Concretely, if x_t solves the forward SDE above with marginal p_t, then the time-reversed process satisfies:
dx=[f(x,t)−g(t)2∇xlogpt(x)]dt+g(t)dwˉ
The derivation sketch: start from Fokker–Planck and rewrite the diffusion term as a divergence using the identity g² ∇²p_t = ∇·(g² ∇p_t) = ∇·(g² p_t ∇log p_t) (using ∇log p_t = ∇p_t / p_t). Then:
∂p_t/∂t = -∇·(f p_t) + ½ ∇·(g² p_t ∇log p_t) = -∇·([f - ½ g² ∇log p_t] p_t) + ½ ∇·(g² p_t ∇log p_t)
Re-grouping, the marginal evolution can be expressed as a different drift v(x,t) = f(x,t) - g(t)² ∇log p_t(x) plus the same diffusion coefficient g. Reversing time τ = T - t flips the sign of the time derivative and gives a reverse-time SDE with drift f - g² ∇log p_t and a reverse-time Brownian motion \bar{w}. The full proof requires care with the reverse Wiener process — see Anderson 1982 or Haussmann & Pardoux (1986) for the rigorous treatment — but the structure of the answer falls out of the Fokker–Planck rewriting above. Yang Song et al. 2021 ("Score-Based Generative Modeling through SDEs", ICLR) rediscovered this for the generative-modeling community and showed that the score ∇log p_t is the only missing ingredient.
The reverse SDE is stochastic — every reverse trajectory is different. There is a beautiful deterministic alternative that shares the same marginal at every t. Start again from Fokker–Planck and write the diffusion term as a divergence in a different way:
dtdx=f(x,t)−21g(t)2∇xlogpt(x)
The derivation is one line of algebra. Take Fokker–Planck and add-then-subtract ½ ∇·(g² p_t ∇log p_t):
∂p_t/∂t = -∇·(f p_t) + ½ ∇·(g² p_t ∇log p_t) − ½ ∇·(g² p_t ∇log p_t) + ½ g² ∇² p_t
The last two terms cancel (using the same identity as before, in the other direction). What remains is:
∂p_t/∂t = -∇·([f − ½ g² ∇log p_t] p_t)
This is the continuity equation for a deterministic flow with velocity v(x,t) = f(x,t) − ½ g(t)² ∇log p_t(x). So the ODE dx/dt = v(x,t) produces the same marginals p_t as the SDE — but with no stochasticity. Same marginals, different sample paths.
Deterministic generation. Fix the initial noise; you always get the same image. Useful for editing, interpolation between samples, and reproducibility.
Exact log-likelihoods. Because it is an ODE, you can apply the instantaneous change-of-variables formula (Chen et al. 2018, Neural ODEs) to compute log p_0(x_0) from log p_T(x_T) and the divergence of the velocity field integrated along the trajectory. SDE-based diffusion only gives a variational lower bound; the probability-flow ODE gives the exact number.
Faster sampling. Off-the-shelf high-order ODE solvers (Heun, RK45, DPM-Solver, DEIS, UniPC) can hit comparable sample quality in 10–50 steps that the SDE needs 1000 to match. DPM-Solver (Lu et al. 2022) and DEIS (Zhang & Chen 2022) are specifically designed exponential-integrator solvers for this ODE.
Everything in the reverse SDE and the probability-flow ODE is known except one quantity: the time-conditioned score ∇_x log p_t(x). Train a network s_θ(x, t) to approximate it. The objective is denoising score matching (Vincent 2011, building on Hyvärinen 2005):
The trick is that the conditional score ∇log q(x_t | x_0) has a closed form whenever q(x_t | x_0) is Gaussian — which it is for both VE and VP. So the regression target is computable in one line, even though the marginal score ∇log p_t(x_t) is intractable. Vincent's theorem says the minimizer of this loss is exactly the marginal score we want.
That's the whole framework. Forward SDE = choice of corruption process. Fokker–Planck = how the marginal evolves. Anderson 1982 = reverse-time SDE exists, and the score is its only unknown. Probability-flow ODE = deterministic sibling, same marginals, exact likelihoods. Denoising score matching = how to train the network without ever touching p_t directly.
Tests · The four target modes should each receive a roughly equal share of samples (around 100 of 400). If one mode hogs everything, your sampler step size is too small. If samples are scattered uniformly, the score network has not learned the modes.
The score is the gradient of the log-density. A vector field that points uphill on the data distribution. Knowing the score is enough to sample from p(x) without ever computing p(x) directly, because the partition function vanishes after taking the gradient.
Langevin dynamics turns the score into samples. Noisy gradient ascent: walk in the direction of the score, add Gaussian noise to avoid getting stuck, and after many steps the trajectory follows p(x). Single-noise-level Langevin fails on real data because the score is unreliable in low-density regions.
Denoising score matching gives you a tractable training target. The conditional score of x̃ = x + σε has the closed form -(x̃-x)/σ². Train a noise-conditional score network on this regression target across many σ, and you get a network that knows the score of the noisy data at every level.
Annealed Langevin samples from coarse to fine. Start from pure noise. Run Langevin at σ_max for a few steps. Drop σ. Repeat. By the end, samples have moved from noise through blurry to sharp data. Step size scales with σ².
Score-based ≡ diffusion. DDPM (predict the noise) and NCSN (predict the score) differ by a scaling factor and are mathematically equivalent under the SDE framework. The reverse SDE is denoising; the only thing it needs is the score; that is what every modern diffusion model is.
Why does the score function avoid the partition function Z = ∫ p(x) dx?
The score gives you everything you need: a vector field, a sampler, and a unified theory of diffusion. Next up: Diffusion Models — DDPM, classifier-free guidance, and the practical sampling tricks (DDIM, DPM-Solver) that make Stable Diffusion run in 20 steps instead of 1000.