Antonio Franca
Two Tracks · 03

The Denoising Perspective

Score matching, diffusion models, and the probability flow ODE

Every physical system evolves toward equilibrium. Drop a drop of ink into water and it spreads, slowly and irreversibly, until the mixture is uniform. The direction of time is the direction of increasing entropy — of increasing disorder. This is not a rule that can be escaped; it is encoded in the statistical mechanics of the world.

But the direction of time can be reversed in the imagination, if not in fact. Given the final, disordered state of the ink and the water, can you recover the initial configuration — the single drop, before it spread? In general, no. But if the process is diffusion — which has a precise stochastic structure — then the reverse process also has a precise stochastic structure, one that can be characterized exactly in terms of how the probability density of configurations evolves.

This is the insight from non-equilibrium thermodynamics that Sohl-Dickstein, Weiss, Maheswaranathan, and Ganguli brought to machine learning in 2015. Data is the signal. Noise is the equilibrium. The generative process is the reversal of diffusion. Their paper — “Deep Unsupervised Learning using Nonequilibrium Thermodynamics” — was largely ignored for five years, during which the transport track was busy building CNFs and Glow. Then, in 2019 and 2020, the idea came back with a force that would permanently reshape the field.

Understanding why requires following a parallel mathematical thread that begins not with transport maps or Jacobians, but with a much older question in statistical estimation: can you estimate a probability density without ever computing its normalization constant?

3.1 The Score Function: Density Without Normalization

Many probability models in machine learning arise as p(x)p~(x)p(x) \propto \tilde{p}(x), where p~(x)\tilde{p}(x) is an unnormalized energy function that can be evaluated pointwise, but the normalization constant Z=p~(x)dxZ = \int \tilde{p}(x)\, dx is intractable. Energy-based models, Boltzmann machines, Markov random fields — all share this structure. Maximum likelihood training requires logp(x)=logp~(x)logZ\log p(x) = \log \tilde{p}(x) - \log Z, and estimating ZZ by Monte Carlo is itself a hard problem.

The escape route is the score function. Define:

where the last equality holds because xlogZ=0\nabla_x \log Z = 0 — the normalization constant does not depend on xx, so its gradient vanishes. The score is the gradient of the log-density, and it is computable from the unnormalized model without knowing ZZ.

The score has a direct geometric interpretation: at any point xx, s(x)s(x) points in the direction of steepest ascent of the log-density — toward the nearest mode, along the path of highest probability. In thermodynamics, this gradient is the force that drives a particle toward equilibrium. In statistics, it encodes the local shape of the distribution more finely than any scalar summary.

[!definition] 3.1 — Score Function and Score Model The score function of a distribution p(x)p(x) is s(x)=xlogp(x)Rds(x) = \nabla_x \log p(x) \in \mathbb{R}^d. A score model sθ(x)s_\theta(x) is a neural network trained to approximate the true score xlogpdata(x)\nabla_x \log p_{\text{data}}(x). A generative model built around a score model is called a score-based generative model.

The immediate challenge: xlogpdata(x)\nabla_x \log p_{\text{data}}(x) is unknown. We have samples from pdatap_{\text{data}}, not an analytic expression for its density. Training sθs_\theta by direct regression against the true score creates a circular dependency — to fit the score, you need the score. Breaking this circle is the central problem of score estimation, and its resolution became the mathematical foundation for the entire denoising track.

3.2 Score Matching: Learning the Gradient (Hyvärinen, 2005)

Aapo Hyvärinen resolved the circular dependency in 2005 with a trick that is simple to state and not obvious to discover. The natural training objective for a score model is the Fisher divergence — the expected squared distance between model and truth:

J(θ)=12Epdata ⁣[sθ(x)xlogpdata(x)2].J(\theta) = \frac{1}{2}\,\mathbb{E}_{p_{\text{data}}}\!\left[\|s_\theta(x) - \nabla_x \log p_{\text{data}}(x)\|^2\right].

This cannot be minimized directly, since xlogpdata(x)\nabla_x \log p_{\text{data}}(x) is unknown. But Hyvärinen showed that, after expanding the square and applying integration by parts, J(θ)J(\theta) is equivalent — up to a constant independent of θ\theta — to:

This implicit score matching objective involves only sθs_\theta and samples from pdatap_{\text{data}} — no knowledge of the true score is required. The derivation proceeds by writing the cross term Ep[sθxlogp]=Ep[sθ(xp)/p]\mathbb{E}_p[s_\theta \cdot \nabla_x \log p] = \mathbb{E}_p[s_\theta \cdot (\nabla_x p)/p] and applying integration by parts, yielding Ep[tr(xsθ)]-\mathbb{E}_p[\mathrm{tr}(\nabla_x s_\theta)], which cancels the unknown score and leaves the computable form (3.2).

The first term, tr(xsθ(x))\mathrm{tr}(\nabla_x s_\theta(x)), is the trace of the Jacobian of sθs_\theta at a data point. Computing it exactly requires dd backward passes through the network — one per output dimension — at a cost of O(d)O(d) times that of a single forward pass. For moderate dimensions this is manageable. For images, where d105d \sim 10^5, it is not. Score matching found immediate applications in density estimation for tabular data, but the Jacobian trace was a hard wall against image-scale problems.

[!remark] 3.1 The integration by parts in Hyvärinen’s derivation requires pdatap_{\text{data}} to vanish at the boundaries of Rd\mathbb{R}^d — or equivalently, that s(x)0s(x) \to 0 sufficiently fast as x\|x\| \to \infty. This is satisfied for all practical continuous data distributions with finite variance. An alternative proof via Stein’s identity makes the same assumption explicit: Ep[g(x)xlogp(x)+xg(x)]=0\mathbb{E}_p[g(x) \nabla_x \log p(x) + \nabla_x g(x)] = 0 for suitable test functions gg. The identity (3.2) is a special case with g=sθg = s_\theta.

3.3 Denoising as Score Matching (Vincent, 2011)

Pascal Vincent’s 2011 insight converted the intractable Jacobian-trace objective into a simple regression, and in doing so revealed the deep connection between density estimation and denoising.

Define the noise-perturbed distribution at noise level σ>0\sigma > 0 as the convolution of pdatap_{\text{data}} with a Gaussian kernel:

pσ(x~)=pdata(x)N(x~;x,σ2I)dx.p_\sigma(\tilde{x}) = \int p_{\text{data}}(x)\, \mathcal{N}(\tilde{x};\, x,\, \sigma^2 I)\, dx.

This is the marginal distribution of x~=x+σϵ\tilde{x} = x + \sigma\epsilon when xpdatax \sim p_{\text{data}} and ϵN(0,I)\epsilon \sim \mathcal{N}(0, I). The score of pσp_\sigma at a noisy point x~\tilde{x} has a closed-form expression in terms of conditional expectations, known as the Tweedie formula:

x~logpσ(x~)=Epdata(xx~) ⁣[x]x~σ2.\nabla_{\tilde{x}} \log p_\sigma(\tilde{x}) = \frac{\mathbb{E}_{p_{\text{data}}(x|\tilde{x})}\![x] - \tilde{x}}{\sigma^2}.

The score at a noisy point x~\tilde{x} points toward the conditional mean of the underlying clean data, scaled by 1/σ21/\sigma^2. Substituting this into the Fisher divergence for pσp_\sigma and expanding gives the denoising score matching objective:

No trace, no Jacobians, no extra backward passes. Each training step requires one forward pass of sθs_\theta to evaluate the squared norm and one backward pass for the gradient with respect to θ\theta. The cost is exactly that of supervised regression — as cheap as training a denoising autoencoder.

[!theorem] 3.1 — Equivalence of Score Matching and Denoising For any σ>0\sigma > 0, the denoising score matching objective JDSM(θ;σ)J_{\text{DSM}}(\theta;\, \sigma) and the implicit score matching objective (3.2) applied to pσp_\sigma have the same minimizer θ\theta^*. In other words, minimizing the denoising loss (3.3) is equivalent to fitting the score of the noise-perturbed distribution, x~logpσ(x~)\nabla_{\tilde{x}} \log p_\sigma(\tilde{x}).

The theorem converts a hard density-estimation problem into an easy regression problem. Train a network to predict the noise ϵ\epsilon given the noisy observation x+σϵx + \sigma\epsilon — and you are, exactly and without approximation, learning the score function of the perturbed distribution. This equivalence is not heuristic; it is algebraic. And it removes the only computational barrier to training score models at scale.

3.4 Non-Equilibrium Thermodynamics: The 2015 Foundation

In 2015 — the same year NICE appeared — Sohl-Dickstein, Weiss, Maheswaranathan, and Ganguli published their thermodynamics paper. Their framing was different from score matching: no Fisher divergence, no Tweedie formula, no explicit score functions. Instead, they described the generative process using the language of Markov chains and the theory of non-equilibrium statistical mechanics.

The forward process. Define a Markov chain x0,x1,,xTx_0, x_1, \ldots, x_T that gradually corrupts a data point x0x_0 into pure noise:

where {βt}t=1T(0,1)\{\beta_t\}_{t=1}^T \subset (0,1) is a noise schedule chosen so that q(xT)N(0,I)q(x_T) \approx \mathcal{N}(0,I). The factor 1βt\sqrt{1-\beta_t} contracts the signal toward zero, while βt\beta_t injects Gaussian noise. The two together preserve the variance: if xt1x_{t-1} has unit variance, so does xtx_t.

The key algebraic property of this specific construction is that the marginal q(xtx0)q(x_t | x_0) is Gaussian for all tt simultaneously:

q(xtx0)=N ⁣(xt;  αˉtx0,  (1αˉt)I),αˉt=s=1t(1βs).q(x_t \mid x_0) = \mathcal{N}\!\left(x_t;\; \sqrt{\bar{\alpha}_t}\, x_0,\; (1-\bar{\alpha}_t)\, I\right), \qquad \bar{\alpha}_t = \prod_{s=1}^t (1 - \beta_s).

Equivalently, any noisy sample at time tt can be written as xt=αˉtx0+1αˉtϵx_t = \sqrt{\bar{\alpha}_t}\, x_0 + \sqrt{1-\bar{\alpha}_t}\, \epsilon for ϵN(0,I)\epsilon \sim \mathcal{N}(0,I). This closed-form marginal — a direct consequence of the Gaussian structure — is the algebraic miracle that makes diffusion models tractable. Without it, sampling from q(xtx0)q(x_t|x_0) would require running the Markov chain step by step.

The reverse process. The time-reversed chain q(xt1xt,x0)q(x_{t-1}|x_t, x_0) has a known closed form when x0x_0 is given:

q(xt1xt,x0)=N ⁣(μ~t(xt,x0),  β~tI),q(x_{t-1} \mid x_t, x_0) = \mathcal{N}\!\left(\tilde{\mu}_t(x_t, x_0),\; \tilde{\beta}_t I\right),

where μ~t\tilde{\mu}_t and β~t\tilde{\beta}_t are functions of βt\beta_t and αˉt\bar{\alpha}_t given by explicit formulas. The problem: x0x_0 is not available at generation time. A learned model pθ(xt1xt)p_\theta(x_{t-1} | x_t) must approximate this reverse without access to x0x_0.

Training maximizes the evidence lower bound on logpθ(x0)\log p_\theta(x_0), which decomposes as a sum of KL divergences between true and learned reverse conditionals at each step:

logpθ(x0)t=1TDKL ⁣(q(xt1xt,x0)    pθ(xt1xt))+C.\log p_\theta(x_0) \geq -\sum_{t=1}^T D_{\mathrm{KL}}\!\left(q(x_{t-1}\mid x_t, x_0)\;\|\; p_\theta(x_{t-1}\mid x_t)\right) + C.

Sohl-Dickstein et al. showed this works in principle. But the results — blurry samples from low-resolution datasets — could not compete with GANs or VAEs of the era. The paper went largely uncited for five years.

[!remark] 3.2 The 2015 paper lay dormant for two reasons. First, the training procedure required evaluating KL divergences at all TT noise levels, and TT had to be large (hundreds to thousands) for the Gaussian approximation of the reverse to be accurate at each step — expensive for networks of the time. Second, the network architecture — fully connected layers on small images — was insufficiently expressive to model the complex reverse transitions. Both problems would be resolved by 2020 through a combination of deeper networks, U-Net architectures, and an objective simplification that removed the explicit KL computation entirely.

3.5 Noise-Conditional Score Networks (Song & Ermon, 2019)

The dormancy ended with Yang Song and Stefano Ermon’s “Generative Modeling by Estimating Gradients of the Data Distribution” (NeurIPS 2019). Their approach dispensed with the Markov chain framing and returned to score matching — but with a key innovation that made it competitive for the first time.

The problem with training a score model at a single noise level σ\sigma is that the score is useful only near that noise level. For generation, you need to start from pure noise (high σ\sigma) and gradually refine toward data (low σ\sigma). A score model trained at high σ\sigma learns nothing about the fine-grained structure of pdatap_{\text{data}}; a model trained at low σ\sigma gives vanishing gradients in the high-noise regime where sampling must begin. Neither alone is sufficient.

Song and Ermon’s solution was to train a single neural network sθ(x,σ)s_\theta(x, \sigma) to estimate xlogpσ(x)\nabla_x \log p_\sigma(x) simultaneously across a geometric sequence of noise levels σ1>σ2>>σL\sigma_1 > \sigma_2 > \cdots > \sigma_L, with σ1\sigma_1 large enough to cover the full data space and σL\sigma_L small enough to be close to pdatap_{\text{data}}. The combined training objective is:

L(θ)=i=1Lλ(σi)JDSM(θ;σi),λ(σi)σi2.\mathcal{L}(\theta) = \sum_{i=1}^L \lambda(\sigma_i)\, J_{\text{DSM}}(\theta;\, \sigma_i), \qquad \lambda(\sigma_i) \propto \sigma_i^2.

The weight λ(σi)σi2\lambda(\sigma_i) \propto \sigma_i^2 normalizes the loss across scales: without it, small-σ\sigma terms would dominate and the network would learn a poor score at high noise levels, where the coarse structure of the distribution lives.

Annealed Langevin dynamics. Sampling with a trained sθs_\theta uses Langevin MCMC — a stochastic gradient ascent on the log-density:

xk+1=xk+ϵi2sθ(xk,σi)+ϵizk,zkN(0,I).x_{k+1} = x_k + \frac{\epsilon_i}{2}\, s_\theta(x_k, \sigma_i) + \sqrt{\epsilon_i}\, z_k, \qquad z_k \sim \mathcal{N}(0, I).

Run at noise level σi\sigma_i with step size ϵi\epsilon_i, this Markov chain converges (in the limit of small ϵi\epsilon_i and many steps) to pσip_{\sigma_i}. The full annealed Langevin procedure runs a short Langevin chain at σ1\sigma_1, then uses the terminal samples as initial conditions for σ2\sigma_2, and so on down to σL0\sigma_L \approx 0. The progressive coarsening from global structure to fine detail mirrors the multi-scale architecture of RealNVP, but here it emerges from the noise schedule rather than from architectural design.

The Noise-Conditional Score Network (NCSN), trained with L=10L = 10 noise levels and a U-Net architecture, produced CIFAR-10 samples competitive with contemporary GANs — the first score-based model to do so. NCSNv2 (Song & Ermon, NeurIPS 2020) refined the noise schedule and architecture further, establishing score-based models as a serious alternative to the adversarial paradigm.

3.6 DDPM: The Simplification That Changed Everything (Ho et al., 2020)

Jonathan Ho, Ajay Jain, and Pieter Abbeel’s “Denoising Diffusion Probabilistic Models” (NeurIPS 2020) returned to the Sohl-Dickstein Markov chain framework — but with a simplification so effective that it transformed diffusion models from a research curiosity into the dominant paradigm for generative modeling.

The key observation: the optimal reverse posterior mean μ~t(xt,x0)\tilde{\mu}_t(x_t, x_0) can be rewritten entirely in terms of the noise ϵ\epsilon that was added at step tt. Because xt=αˉtx0+1αˉtϵx_t = \sqrt{\bar{\alpha}_t}\, x_0 + \sqrt{1-\bar{\alpha}_t}\, \epsilon, we have x0=(xt1αˉtϵ)/αˉtx_0 = (x_t - \sqrt{1-\bar{\alpha}_t}\, \epsilon)/\sqrt{\bar{\alpha}_t}, and substituting:

μ~t(xt,x0)=11βt ⁣(xtβt1αˉtϵ).\tilde{\mu}_t(x_t, x_0) = \frac{1}{\sqrt{1-\beta_t}}\!\left(x_t - \frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}}\, \epsilon\right).

Learning the reverse process is therefore equivalent to learning to predict ϵ\epsilon from xtx_t. Parametrizing the learned reverse as pθ(xt1xt)=N(μθ(xt,t),βtI)p_\theta(x_{t-1}|x_t) = \mathcal{N}(\mu_\theta(x_t, t), \beta_t I) with μθ\mu_\theta in the noise-prediction form, and substituting into the ELBO, Ho et al. showed that the dominant term simplifies to:

This training objective is extraordinarily simple: sample a data point x0x_0, sample a timestep tt uniformly from {1,,T}\{1, \ldots, T\}, corrupt x0x_0 with the appropriate Gaussian noise to obtain xtx_t, and train the network to predict which noise was added. No explicit KL divergences, no Jacobians, no ODE solvers, no Langevin chains during training. Just noise prediction by regression.

The connection to score matching is exact: since xt=αˉtx0+1αˉtϵx_t = \sqrt{\bar{\alpha}_t}\, x_0 + \sqrt{1-\bar{\alpha}_t}\,\epsilon, the optimal noise predictor satisfies

ϵθ(xt,t)1αˉt  xtlogpt(xt),\epsilon_\theta(x_t, t) \approx -\sqrt{1 - \bar{\alpha}_t}\;\nabla_{x_t} \log p_t(x_t),

making the DDPM training objective a disguised form of denoising score matching — the same objective Vincent derived in 2011 from the Fisher divergence, applied now at discrete timesteps with the specific Gaussian forward process (3.4). NCSN and DDPM, developed from entirely different motivations, are training the same function.

The results were immediately recognized as a step change. FID scores on CIFAR-10 and LSUN surpassed contemporaneous GANs. Samples from 256×256256 \times 256 face generation were visually indistinguishable from photographs. And — unlike GANs — training was stable, reproducible, and required no adversarial game between competing networks. Within months, DALL-E 2, Stable Diffusion, and Imagen had all been built on DDPM foundations.

[!remark] 3.3 The DDPM sampling procedure generates samples by running the reverse chain pθ(xt1xt)p_\theta(x_{t-1}|x_t) backward from xTN(0,I)x_T \sim \mathcal{N}(0,I) to x0x_0. This requires T=1000T = 1000 sequential neural network evaluations — one per denoising step — making generation slow compared to a single normalizing flow forward pass. The problem of reducing the number of function evaluations (NFE) without degrading sample quality became a central research direction. Song et al. (DDIM, 2020) showed that TT could be reduced to 50 with minimal quality loss by using a deterministic non-Markovian sampler; the systematic theory of fast sampling is the subject of Chapter 8.

3.7 The Score SDE: A Continuous-Time Unification (Song et al., 2021)

By early 2021, the field had two distinct score-based approaches — NCSN (continuous noise levels, annealed Langevin sampling) and DDPM (discrete Markov chain, reverse diffusion sampling) — that were clearly related but whose precise connection had never been articulated. Yang Song, Jascha Sohl-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole’s “Score-Based Generative Modeling through Stochastic Differential Equations” (ICLR 2021, outstanding paper) provided the unification by passing to the continuous-time limit.

As the number of discrete noise steps TT \to \infty and the step size simultaneously goes to zero, the forward Markov chain (3.4) converges to a stochastic differential equation (SDE). In general, the forward process is:

where WW is a standard Brownian motion, f:Rd×[0,T]Rdf : \mathbb{R}^d \times [0,T] \to \mathbb{R}^d is the drift, and g:[0,T]Rg : [0,T] \to \mathbb{R} is the diffusion coefficient. The SDE (3.6) defines a continuous-time stochastic process {xt}0tT\{x_t\}_{0 \leq t \leq T} whose marginal distribution at time TT approaches N(0,I)\mathcal{N}(0,I) for appropriate choices of ff and gg.

Song et al. identified two canonical cases corresponding to the two established approaches:

VP-SDE (Variance Preserving, corresponding to DDPM): f(x,t)=12β(t)xf(x,t) = -\tfrac{1}{2}\beta(t)\, x and g(t)=β(t)g(t) = \sqrt{\beta(t)}, where β(t)\beta(t) is the continuous-time noise schedule interpolating the discrete {βt}\{\beta_t\}. The signal is contracted linearly while noise is injected proportionally — a variance-preserving process because the variance of xtx_t remains bounded.

VE-SDE (Variance Exploding, corresponding to NCSN): f(x,t)=0f(x,t) = 0 and g(t)=d[σ2(t)]/dtg(t) = \sqrt{d[\sigma^2(t)]/dt}, where σ(t)\sigma(t) is the noise level schedule. The drift is zero — the signal is not contracted — and noise is injected until the variance explodes. The endpoint distribution is not Gaussian in the strict sense but has variance σT21\sigma_T^2 \gg 1, which dominates any fixed-variance data distribution.

Both SDEs transform x0pdatax_0 \sim p_{\text{data}} into xTN(0,I)x_T \approx \mathcal{N}(0, I) as TT \to \infty. The key quantity connecting forward and reverse is the score of the marginal distribution at each time: xlogpt(x)\nabla_x \log p_t(x).

The reverse-time SDE. Brian Anderson (1982) showed that every forward Itô SDE has a corresponding reverse-time SDE that runs from t=Tt = T back to t=0t = 0 with the same marginal distributions. Applied to (3.6), the reverse process satisfies:

where Wˉ\bar{W} is a Brownian motion running backward in time. The score function xlogpt(x)\nabla_x \log p_t(x) appears explicitly as the correction term that reverses the forward diffusion. Approximating this score with a learned network sθ(x,t)xlogpt(x)s_\theta(x, t) \approx \nabla_x \log p_t(x) — trained by denoising score matching across all tt — gives a complete generative model.

[!definition] 3.2 — Score-Based SDE Model A score-based SDE generative model consists of a forward SDE (3.6), a learned score network sθ(x,t)xlogpt(x)s_\theta(x, t) \approx \nabla_x \log p_t(x) trained by denoising score matching, and the reverse SDE (3.7) with xlogpt\nabla_x \log p_t replaced by sθs_\theta. Sampling proceeds by numerically integrating the reverse SDE from xTN(0,I)x_T \sim \mathcal{N}(0, I) to x0x_0.

The continuous-time training objective is a natural generalization of (3.3): for a weight function λ(t)>0\lambda(t) > 0,

L(θ)=EtU[0,T],  x0pdata,  xtq(xtx0) ⁣[λ(t)sθ(xt,t)xtlogq(xtx0)2].\mathcal{L}(\theta) = \mathbb{E}_{t \sim \mathcal{U}[0,T],\; x_0 \sim p_{\text{data}},\; x_t \sim q(x_t|x_0)}\!\left[\lambda(t)\,\|s_\theta(x_t, t) - \nabla_{x_t} \log q(x_t \mid x_0)\|^2\right].

Because q(xtx0)q(x_t | x_0) is Gaussian for both the VP- and VE-SDE, xtlogq(xtx0)\nabla_{x_t} \log q(x_t | x_0) is a simple linear function of xtx_t and x0x_0 — and the objective reduces to denoising score matching. No Jacobians. No ODE simulation during training. The DDPM objective (3.5) is the VP-SDE version of this, discretized to T=1000T = 1000 steps.

3.8 The Probability Flow ODE

The most striking result of Song et al. (2021) was not the unification of NCSN and DDPM — that was expected — but the discovery that every diffusion model has an equivalent deterministic counterpart.

For any forward SDE (3.6) with marginals {pt}\{p_t\}, there exists an ordinary differential equation — no Brownian motion, no randomness — whose solution trajectories {xt}\{x_t\} share the same marginal distribution ptp_t at every time tt:

This is the probability flow ODE. It is a deterministic ordinary differential equation whose right-hand side involves the score function. Compared to the reverse SDE (3.7), the noise term g(t)dWˉg(t)\, d\bar{W} has been removed and the score drift has been halved — the factor of 12\tfrac{1}{2} is not an approximation but an exact result from the Fokker-Planck equation governing the evolution of ptp_t.

[!theorem] 3.2 — Probability Flow ODE Let {pt}0tT\{p_t\}_{0 \leq t \leq T} be the marginal densities of the forward SDE (3.6). The probability flow ODE (3.8) defines a deterministic flow ϕt:RdRd\phi_t : \mathbb{R}^d \to \mathbb{R}^d such that ϕt(x0)pt\phi_t(x_0) \sim p_t whenever x0p0x_0 \sim p_0. Equivalently, integrating the ODE backward from xTN(0,I)x_T \sim \mathcal{N}(0,I) produces a sample distributed as p0=pdatap_0 = p_{\text{data}}.

The implications of this theorem are profound, and they become clearer when stated as a connection between the two tracks. The probability flow ODE (3.8) is exactly a continuous normalizing flow in the sense of Chapter 2. Its vector field is

ut(x)=f(x,t)12g(t)2xlogpt(x).u_t(x) = f(x,t) - \tfrac{1}{2}\,g(t)^2\, \nabla_x \log p_t(x).

If we know the exact score xlogpt\nabla_x \log p_t at every time and position, we can integrate this ODE forward or backward to generate samples or compute exact likelihoods — exactly as FFJORD does, with the same NFE cost. The simulation bottleneck of Chapter 2 reappears: generating one sample requires numerically integrating the probability flow ODE, which costs as many network evaluations as a CNF.

But the training is free of that bottleneck. DDPM trains by denoising regression. The probability flow ODE is a consequence of that training — the deterministic dynamics that emerge from a trained score model, without any ODE solve during the training loop.

[!remark] 3.4 The probability flow ODE enables deterministic sampling: fixing the initial condition xTx_T fixes the generated output completely. This gives diffusion models a latent space — the space of initial Gaussian noise vectors — with meaningful geometry. Two data points x0(1)x^{(1)}_0 and x0(2)x^{(2)}_0 correspond to initial conditions xT(1)x^{(1)}_T and xT(2)x^{(2)}_T; interpolating linearly between these in the latent space and integrating the ODE produces a smooth perceptual morph between the two data points. This latent interpolation, unavailable from the stochastic reverse SDE sampler, made the probability flow ODE an immediate tool for image editing applications (DALL-E 2, InstructPix2Pix, and their successors all use related ideas).

The probability flow ODE also enables exact likelihood computation: running the ODE forward from x0x_0 to xTx_T and applying the instantaneous change-of-variables formula (Theorem 2.1) gives logp0(x0)\log p_0(x_0) exactly, as the sum of logN(xT;0,I)\log \mathcal{N}(x_T; 0, I) and the integrated divergence of utu_t. Diffusion models — which were thought to offer no tractable likelihood — turn out to have exactly computable likelihoods via their probability flow ODE, using the same Hutchinson trace estimator as FFJORD.

3.9 What the Denoising Track Got Right

The denoising track arrived at the same destination as the transport track but via different mathematics, and its tools had properties that the transport track’s lacked. Understanding the contrast clarifies why diffusion models displaced CNFs in practice — and why the question of Chapter 4 is so natural.

Scalability by construction. The DDPM training objective (3.5) requires exactly one forward pass of ϵθ\epsilon_\theta per training step, regardless of data dimension, distribution complexity, or model depth. There is no ODE to simulate, no Jacobian to estimate, no adjoint to integrate backward. The cost of a training step for DDPM equals the cost of supervised regression — and this cost does not grow as the target distribution becomes more complex. The simulation bottleneck of Chapter 2 is entirely absent from training.

The score as the universal object. The score function xlogpt(x)\nabla_x \log p_t(x) is not merely a computational device — it is the sufficient statistic for the reverse dynamics. Everything needed to generate from any distribution reachable by diffusion is encoded in the score, at every noise level. This is both mathematically clean and practically powerful: it means the architecture and training procedure are the same regardless of the target, and theoretical guarantees can be stated in terms of score estimation accuracy alone.

The SDE framework and its connections. By casting diffusion models as SDEs, Song et al. connected the denoising track to a rich body of mathematical theory: Fokker-Planck equations, Langevin dynamics, and — via the probability flow ODE — continuous normalizing flows. The framework is extensible: new noise processes, new reverse samplers (predictor-corrector, DDIM, DPM-Solver), and new training objectives all fit within the SDE language.

The implicit multi-scale decomposition. Training a score model across all noise levels t[0,T]t \in [0,T] implicitly teaches the model to decompose the data distribution at all scales simultaneously. At high tt, the score model learns global structure — which region of the space a sample belongs to. At low tt, it learns fine-grained local detail. This coarse-to-fine decomposition happens automatically from the training objective, without architectural engineering — and gives diffusion models their characteristic ability to generate globally coherent samples with faithful local texture.

3.10 Two Tracks, One Step From Convergence

Chapter 1 left us with normalizing flows: exact likelihoods, tractable training, but restricted architectures that impose artificial coordinate orderings on the data. Chapter 2 removed the architectural restriction via continuous-time dynamics, but introduced the simulation bottleneck: training a CNF requires ODE simulation, and this cost grows with distribution complexity. This chapter introduced the denoising track, which removes the simulation bottleneck from training entirely — at the cost of slow, multi-step sampling.

The situation in late 2021 was this:

  • CNFs: free architectures, exact likelihoods, simulation bottleneck at training.
  • Diffusion models: cheap training (one forward pass per step), many function evaluations at sampling.

Both bottlenecks trace to the same root: both methods define a complex, curved path from noise to data — the Riemannian path induced by the SDE, or the constrained path imposed by the Jacobian structure — and that path requires many steps to traverse accurately. Is there a simpler path?

The answer is yes, and it is almost embarrassingly direct: define the transportation path yourself, as a straight line between noise and data, and derive the vector field from that path. No diffusion, no score functions, no ODE simulation during training. Just regression.

This is the flow matching insight of Lipman et al. (2022) and Albergo & Vanden-Eijnden (2022). The probability flow ODE (3.8) is a particular choice of continuous normalizing flow — the one induced by a specific diffusion SDE. Flow matching says: why not choose a different flow, one that is simpler to train and cheaper to sample? The score function is not fundamental to the vector field; it is one way to derive a vector field that transports p0p_0 to p1p_1. There are others.

Understanding precisely what those others are, and how they unify the entire landscape of generative transport, is the subject of Chapter 4.

3.11 Historical Notes

The mathematical foundation for score matching is Hyvärinen (2005) (“Estimation of Non-Normalized Statistical Models by Score Matching”, JMLR), who proved the implicit score matching identity (3.2) and demonstrated its use for estimating energy-based models. Vincent (2011) (“A Connection Between Score Matching and Denoising Autoencoders”, Neural Computation) established Theorem 3.1 — that denoising and score matching are equivalent — converting the intractable Jacobian-trace objective into a simple regression.

The Tweedie formula connecting the posterior mean to the score appears in Robbins (1956) and was connected to diffusion models explicitly by Efron (2011) (“Tweedie’s Formula and Selection Bias”, JASA) and formalized in the diffusion context by Luo (2022) (“Understanding Diffusion Models: A Unified Perspective”).

The non-equilibrium thermodynamics framing is due to Sohl-Dickstein, Weiss, Maheswaranathan, and Ganguli (ICML 2015) (“Deep Unsupervised Learning using Nonequilibrium Thermodynamics”). The physics background — fluctuation theorems, Jarzynski’s equality, the thermodynamic arrow of time — is surveyed in Seifert (2012) (“Stochastic thermodynamics, fluctuation theorems and molecular machines”, Reports on Progress in Physics).

The revival of diffusion models began with Song and Ermon (NeurIPS 2019) (“Generative Modeling by Estimating Gradients of the Data Distribution”) and NCSNv2 (NeurIPS 2020) (“Improved Techniques for Training Score-Based Generative Models”). The DDPM simplification is from Ho, Jain, and Abbeel (NeurIPS 2020). The v-prediction reparametrization — predicting the velocity v=αˉtϵ1αˉtx0v = \sqrt{\bar\alpha_t}\epsilon - \sqrt{1-\bar\alpha_t}x_0 — appears in Salimans and Ho (2022) (“Progressive Distillation for Fast Sampling of Diffusion Models”).

The continuous-time SDE unification is Song, Sohl-Dickstein, Kingma, Kumar, Ermon, and Poole (ICLR 2021, outstanding paper) (“Score-Based Generative Modeling through Stochastic Differential Equations”). The reverse-time SDE identity (3.7) is Anderson’s result (Anderson, 1982, “Reverse-time diffusion equation models”, Stochastic Processes and their Applications). The mathematical background for SDEs and Fokker-Planck equations is in Øksendal (2003) (Stochastic Differential Equations) and Risken (1989) (The Fokker-Planck Equation).

Theoretical convergence rates for score-based models — relating total variation distance between model and pdatap_{\text{data}} to score estimation error and discretization step size — are developed in De Bortoli (2022), Chen, Chewi, Li, Li, Salim, and Zhang (2022), and Lee, Lu, and Tan (2022). The sharpest current results give polynomial dependence on dimension under log-Sobolev or Poincaré inequalities on pdatap_{\text{data}}.

The large-scale applications that made diffusion models dominant — DALL-E 2 (Ramesh et al., 2022), Stable Diffusion (Rombach et al., CVPR 2022), Imagen (Saharia et al., NeurIPS 2022) — all build on the DDPM backbone with classifier-free guidance (Ho & Salimans, 2021) and latent-space diffusion (Rombach et al.). The diffusion transformer architecture (DiT, Peebles & Xie, 2023), which replaced U-Nets with transformers as the score network backbone, is now the standard architecture for large-scale image and video generation.