4. Continuous time
The reverse SDE and the probability-flow ODE Draft
The forward equation, three times over
As DDPM's steps shrink, its forward chain approaches the Ornstein–Uhlenbeck process. Working in continuous time makes the reverse process easier to see. Anderson's theorem gives a reverse diffusion whose drift uses the score. A related ordinary differential equation follows the same distributions without adding fresh noise.
Both constructions follow from the Fokker–Planck calculation in Section 1.2. The deterministic flow also gives us a way to evaluate likelihoods by tracking how volumes change along a path.
The forward equation is the one of Proposition 1.1.5, \(dx = f_t\, x\, dt + g_t\, dW\) with \(f_t = \dot\alpha_t/\alpha_t\) and \(g_t^2 = -\sigma_t^2\dot\lambda_t\). The table gives its coefficients on the schedules of the series; the names VP and VE are those of Song et al.
| schedule | \(\alpha_t\) | \(\sigma_t\) | \(f_t\) | \(g_t^2\) | the forward SDE |
|---|---|---|---|---|---|
| variance preserving (VP) | \(e^{-B(t)/2}\) | \(\sqrt{1 - e^{-B(t)}}\) | \(-\tfrac12\beta_t\) | \(\beta_t\) | \(dx = -\tfrac12\beta_t x\, dt + \sqrt{\beta_t}\, dW\) |
| variance exploding (VE) | \(1\) | \(\sigma_t\) | \(0\) | \(\tfrac{d}{dt}\sigma_t^2\) | \(dx = \sqrt{\tfrac{d}{dt}\sigma_t^2}\, dW\) |
| the EDM clock | \(1\) | \(t\) | \(0\) | \(2t\) | \(dx = \sqrt{2t}\, dW\) |
| the linear path | \(1 - t\) | \(t\) | \(-\tfrac{1}{1 - t}\) | \(\tfrac{2t}{1 - t}\) | \(dx = -\tfrac{x}{1 - t}\, dt + \sqrt{\tfrac{2t}{1 - t}}\, dW\) |
Remark 4.1.1 (takeaway: one process, many magnifications). In the scaled coordinate \(\bar x = x/\alpha_t\) every row of the table is the same process: \(d\bar x = (g_t/\alpha_t)\, dW\), a Brownian motion run on the clock \(\rho_t = \sigma_t/\alpha_t\), since \((g_t/\alpha_t)^2 = \frac{d}{dt}\rho_t^2\). The schedule chooses a magnification \(\alpha_t\) and a clock, nothing else; the data law seen through the noise at ratio \(\rho\) is \(p_0 * \mathcal N(0, \rho^2 I)\) on every schedule. Karras et al. made this the organizing principle of their design space, and Section 5.2 uses it to show that DDIM is Euler's method.
Anderson's theorem and the family of reverse processes
Reversing a diffusion requires more than reversing its drift. We must also account for how the added noise has spread the density. The score supplies that correction. To keep the signs clear, we first write the reversed process using a new clock that increases from zero.
Reversal is an inference problem. A particle observed at a particular location could have arrived from several neighbouring locations, and the crowded neighbours supplied more candidates. Bayes' rule therefore biases the reverse motion toward where probability was concentrated. The score expresses that bias locally. Reversing only the sign of the forward drift would miss it.
left the law of the origin y one step h = 0.02 before t = 0.1 for a particle seen at x = 0.7 on R2:
the candidates' law p_{t-h} (grey), the likelihood of x from y as a dashed bell at x/a = 0.713
(sd 0.196) and the posterior over origins (green), mean 0.616 and sd 0.208, weights 0.000,
0.818, 0.182 on the three modes; on the axis the particle 0.700, the drift-undone point 0.713,
the Bayes origin 0.616 and Anderson's Euler point 0.607, the pull (b^2/a) s_t(x) = −0.097
right the pull per unit time against h, −4.86 at h = 0.02 on its way to g^2 s_t(x) = −5.38, and the
spread of the origin per sqrt(h), 1.47 on its way to g = 1.45
(the live figure needs JavaScript)
Theorem 4.2.1 (Anderson). Let \(x_t\) solve \(dx = f(x, t)\, dt + g(t)\, dW\) on \([0, 1]\) with marginals \(p_t\) positive and \(C^2\), and \(f\) of class \(C^1\). Then the process \(y_s = x_{1 - s}\), \(s \in [0, 1]\), is a weak solution of
\[\begin{aligned} dy_s &= \Bigl[-f(y_s, 1 - s) + g(1 - s)^2\, \nabla \log p_{1 - s}(y_s)\Bigr] ds \\ &\quad + g(1 - s)\, d\bar W_s, \qquad y_0 \sim p_1, \end{aligned} \tag{4.1}\]for a Brownian motion \(\bar W\). Written in the original time, with \(dt < 0\): \(dx = [f(x, t) - g(t)^2 \nabla \log p_t(x)]\, dt + g(t)\, d\bar W\).
Proof of the marginal statement. What a sampler needs is that the solution of (4.1) started from \(p_1\) has law \(p_{1 - s}\) at time \(s\), and this follows from the Fokker–Planck equation alone. Let \(q_s = p_{1 - s}\). By Theorem 1.2.2 for the forward process, \(\partial_s q_s = -(\partial_t p)\big\rvert_{t = 1 - s} = \nabla \cdot (f p) - \tfrac12 g^2 \Delta p\), everything evaluated at \(t = 1 - s\). The Fokker–Planck equation of (4.1), whose drift is \(b = -f + g^2 \nabla \log q\), applied to \(q_s\), has right side
\[\begin{aligned} -\nabla \cdot (b\, q) + \tfrac12 g^2 \Delta q &= \nabla \cdot (f q) - \nabla \cdot \bigl(g^2\, q\, \nabla \log q\bigr) + \tfrac12 g^2 \Delta q \\ &= \nabla \cdot (f q) - g^2 \Delta q + \tfrac12 g^2 \Delta q \\ &= \nabla \cdot (f q) - \tfrac12 g^2 \Delta q, \end{aligned}\]using \(q \nabla \log q = \nabla q\). This is \(\partial_s q_s\). So \(q_s\) solves the Fokker–Planck equation of (4.1) with the initial law \(q_0 = p_1\), and when that equation has a unique solution, which holds for the Gaussian mixtures of the series and under the regularity assumed, the marginals of (4.1) are \(q_s\). ∎
The calculation proves that the proposed reverse equation has the right distribution at every time. Anderson's full theorem goes further: the time-reversed process itself solves that equation. Haussmann and Pardoux give general conditions for this stronger statement. For generation, the marginal result above supplies what we need.B. D. O. Anderson, "Reverse-time diffusion equation models", Stochastic Processes and their Applications 12 (1982), 313–326; U. G. Haussmann and E. Pardoux, "Time reversal of diffusions", Annals of Probability 14 (1986), 1188–1205. The use for generative modelling is Song et al. (2021), Section 3.2.
In the increasing reverse clock, the drift has two parts: the negative of the forward drift, and \(g^2\) times the score. For our contracting forward processes, the first undoes the contraction and the second counters the spread caused by diffusion. Anderson's formula fixes these coefficients for the time-reversed process. If we ask only for the same marginals, however, there is more freedom.
The coefficient of the score also pays for the noise added during reversal. A deterministic flow needs enough inward current to undo the forward spreading. A stochastic reverse process adds new spreading of its own, which requires additional inward current. This is why changing the noise amplitude changes the drift, and why its contribution involves the square of the amplitude.
Theorem 4.2.2 (the family of reverse processes). Under the assumptions of Theorem 4.2.1, for every \(\lambda \ge 0\) the process
\[\begin{aligned} dx &= \Bigl[f(x, t) - \frac{1 + \lambda^2}{2}\, g(t)^2\, \nabla \log p_t(x)\Bigr] dt \\ &\quad + \lambda\, g(t)\, d\bar W, \qquad dt < 0, \quad x_1 \sim p_1, \end{aligned} \tag{4.2}\]has marginals \(p_t\). The case \(\lambda = 1\) is Anderson's reverse SDE and the case \(\lambda = 0\) is the ordinary differential equation \(dx/dt = f(x, t) - \tfrac12 g(t)^2 \nabla \log p_t(x)\).
Proof. The same computation with \(b = -f + \tfrac{1 + \lambda^2}{2} g^2 \nabla \log q\) and diffusion coefficient \(\lambda g\) in the reversed time: \(-\nabla \cdot (b q) + \tfrac{\lambda^2}{2} g^2 \Delta q = \nabla \cdot (f q) - \tfrac{1 + \lambda^2}{2} g^2 \Delta q + \tfrac{\lambda^2}{2} g^2 \Delta q = \nabla \cdot (f q) - \tfrac12 g^2 \Delta q\), which is \(\partial_s q_s\) for every \(\lambda\). ∎
So the one-time laws of the forward process are shared by a one-parameter family of reverse processes, from the deterministic flow to processes noisier than Anderson's; this is Remark 1.3.6 made precise. Karras et al. write the same family with a time-dependent rate, \(dx = -\dot\sigma\sigma \nabla \log p\, dt \pm \beta_t \sigma^2 \nabla \log p\, dt + \sqrt{2\beta_t}\, \sigma\, dW\) on the clock \(\sigma\), the \(\beta_t\) term being a Langevin diffusion that leaves \(p_t\) invariant at a fixed time; in the notation of (4.2) their \(\beta_t\sigma^2\) is \(\tfrac{\lambda^2}{2} g^2 = \lambda^2 \sigma \dot\sigma\).Karras, Aittala, Aila and Laine (2022), Section 4, Equation 6, and their Algorithm 2, whose parameter \(S_{\mathrm{churn}}\) sets how much of this noise is injected at each step. Stochastic interpolants, Section 6.5, reach the same family from the side of interpolation: M. S. Albergo, N. M. Boffi and E. Vanden-Eijnden, "Stochastic interpolants: a unifying framework for flows and diffusions", Journal of Machine Learning Research 26 (2025), 1–80 (arXiv 2303.08797), "a family of forward and backward Fokker-Planck equations with tunable diffusion coefficient".
left 40 of 2,000 paths of the reverse process (4.2) on R2 for the chosen lambda, the same Brownian
increments for every lambda, over the exact density p_t, which does not depend on lambda
right the exact law of each discretized process's endpoints against p_0, with their W_1 and KL,
or KL(p_t || q_t) along the run for the four lambdas
(the live figure needs JavaScript)
Example 4.2.3 (the family on R2). In the figure's default run the exact laws of the four discretized processes' endpoints are within \(0.0025\), \(0.0016\), \(0.0022\) and \(0.0056\) of R2 in Wasserstein-1 distance for \(\lambda = 0, 0.5, 1, 2\) (relative entropies \(0.0001\), \(0.0001\), \(0.0002\) and \(0.0007\)), the discretization error of 300 Euler–Maruyama steps, and the shading behind the paths, the exact \(p_t\), is the same picture four times. The paths are not: at \(\lambda = 0\) they are the smooth lines of the flow and never cross; at \(\lambda = 1\) their quadratic variation over the run is \(9.87\) and they forget where they started; at \(\lambda = 2\) the score term is \(2.5\) times the flow's and the noise twice Anderson's, and the paths are rougher still.
The second view of the right panel follows the relative entropy of each discretized process's marginal from the exact \(p_t\) along the run: for every member, the flow included, it rises from about \(10^{-6}\) at the start as the drift steepens, peaks near \(t = 0.05\) and falls over the last steps to the four values above, the member with \(\lambda = 0.5\) ending below the flow and the two noisier members above it; with the exact score the Langevin term of Remark 4.2.5 has nothing to correct. The rank correlation between a path's start and its end, one for the flow, is \(-0.007\) at \(\lambda = 1\) and \(0.012\) at \(\lambda = 2\). The shared marginals are a theorem; the distinct paths are a choice, and Section 5 shows that the choice is not free of consequences once the time is discretized. ◁
Algorithm 4.2.4 REVERSE-EULER-MARUYAMA
(the family (4.2); lambda = 0 is Euler for the flow)
Picture: one step backward in time, from t to s < t.
The drift moves the point against the forward drift
and up the score; the noise, if any, is scaled by
lambda and by the step's length.
x_t --[ (f x - (1+lam^2)/2 g^2 s_theta)(s - t) ]--> .
|
+ lam g sqrt(t - s) z
v
x_s
Input a score s_theta(x, t); the coefficients f_t, g_t;
a grid 1 - eps = t_0 > t_1 > ... > t_N = eps;
the noise multiplier lambda >= 0.
Output a sample.
1. x ~ N(0, I)
2. for k = 0, 1, ..., N-1:
3. h <- t_{k+1} - t_k (negative)
4. z ~ N(0, I)
5. x <- x + h (f_{t_k} x - (1 + lam^2)/2 g_{t_k}^2 s_theta(x, t_k))
+ lam g_{t_k} sqrt(|h|) z
6. return D_theta(x, t_N) (the final denoising, Tweedie)
The program below runs the algorithm on R2 with its own random numbers, so that its Wasserstein distances differ from the figure's by the sampling noise of 2,000 draws and not otherwise.
# The family of reverse processes on R2, variance preserving:
# dx = [ f x - (1 + lam^2)/2 g^2 score ] dt + lam g dW, dt < 0,
# by Euler-Maruyama from 2,000 draws of N(0, 1), for lam = 0, 1, 2, with the
# Wasserstein-1 distance of the endpoints to R2 (exact CDF on a fine grid).
#
# Picture: three processes, one family of densities.
#
# lam = 0 smooth flow lines --------\___
# lam = 1 Anderson's reverse SDE ~~\/~\_/~\__
# lam = 2 twice the noise /\/\/\/\/\__
# all three end in p_0, and all three pass
# through p_t at every t
import math, random
W, M, S = (0.45, 0.35, 0.20), (-1.4, 0.2, 1.6), (0.25, 0.15, 0.30)
def B(t): return 0.1 * t + 9.95 * t * t
def beta(t): return 0.1 + 19.9 * t
def alpha(t): return math.exp(-0.5 * B(t))
def sigma(t): return math.sqrt(1 - math.exp(-B(t)))
def score(x, a, s):
logits, grads = [], []
for w, m, sd in zip(W, M, S):
v = a * a * sd * sd + s * s
logits.append(math.log(w) - 0.5 * (x - a * m) ** 2 / v - 0.5 * math.log(2 * math.pi * v)); grads.append((a * m - x) / v)
top = max(logits); z = sum(math.exp(l - top) for l in logits)
return sum(math.exp(l - top) / z * g for l, g in zip(logits, grads))
def cdf0(x):
return sum(w * 0.5 * (1 + math.erf((x - m) / (sd * math.sqrt(2)))) for w, m, sd in zip(W, M, S))
def w1(points):
pts = sorted(points); n = len(pts); grid = [-4 + 8 * i / 1600 for i in range(1601)]
acc, j = 0.0, 0
for x in grid:
while j < n and pts[j] <= x: j += 1
acc += abs(cdf0(x) - j / n) * (8 / 1600)
return acc
steps, n = 300, 2000
times = [1 - 1e-3 - (1 - 2e-3) * i / steps for i in range(steps + 1)]
for lam in (0.0, 1.0, 2.0):
rng = random.Random(5)
x = [rng.gauss(0, 1) for _ in range(n)]
for t0, t1 in zip(times, times[1:]):
h, a, s, g2, f = t1 - t0, alpha(t0), sigma(t0), beta(t0), -0.5 * beta(t0)
noise = lam * math.sqrt(g2 * abs(h))
x = [xi + h * (f * xi - 0.5 * (1 + lam * lam) * g2 * score(xi, a, s)) + noise * rng.gauss(0, 1) for xi in x]
print(f"lambda = {lam}: W_1(endpoints, R2) = {w1(x):.3f}")
lambda = 0.0: W_1(endpoints, R2) = 0.036
lambda = 1.0: W_1(endpoints, R2) = 0.031
lambda = 2.0: W_1(endpoints, R2) = 0.037
Remark 4.2.5 (takeaway: why one would add noise at all). With an exact score and exact integration, every member of (4.2) has the right marginals. At a finite step count they differ: Example 4.2.3 gives errors between \(0.002\) and \(0.006\) in Wasserstein-1 distance. A learned score adds another reason for the choice to matter.
Freeze the noise level at \(t\) and consider an auxiliary Langevin process, with its own increasing clock \(h\),
\[dx=\tfrac{\lambda^2}{2}g_t^2s_t(x)\,dh+\lambda g_t\,dW_h.\]Its exact stationary law is \(p_t\): the score drift balances its spreading. Mixing under this process can reduce a distributional error left by earlier steps. The rate depends on the geometry of \(p_t\), and with a learned score, \(p_t\) need no longer be stationary, so fresh noise does not guarantee an improvement.
Applying the same invertible deterministic flow to two initial laws preserves their relative entropy, even though distances between individual trajectories can grow or shrink. The stochastic correction offers a different way to alter that discrepancy. Song et al. built this into their predictor–corrector samplers, and Karras et al. found that on ImageNet at \(64 \times 64\) a stochastic sampler with 511 network evaluations reached an FID of \(1.36\) where their best deterministic one reached \(1.55\) at far fewer evaluations; the figure of Section 5.3 shows the same two effects on R3, the deterministic samplers converging faster and the stochastic ones recovering from their coarse steps.Song et al. (2021), Section 4.2 and Algorithm 1 for predictor–corrector sampling; Karras et al. (2022), Section 4 and Figure 5: stochasticity "corrects errors made in earlier sampling steps", with ImageNet-64 "FID 1.36, NFE 511", the "new SOTA of 1.36" of the abstract, against \(1.55\) for the deterministic sampler on the pretrained model.
left p_t of R2 at t = 0.1 under the variance-preserving schedule, the start q_0, which is p_t shifted
by 0.2, and the law q_h that the frozen-time process dx = D s_t(x) dh + sqrt(2 D) dW_h, with
D = g^2/2 = 1.045, has carried it to by h = 0.1
right KL(q_h || p_t) against h on a log scale for four starts (none, the shift, a blur of 0.15 and the
weights (0.30, 0.50, 0.20)): the shift 0.0739 at h = 0, 0.0139 at h = 0.1 and 0.0027 at h = 1,
halving by h = 0.032; the blur 0.0034 to 0.0003 by h = 0.1, halving by 0.023; the misweighting
0.0517 to 0.0403 by h = 0.1 and 0.0090 by h = 1, halving by 0.32; p_t itself below 1e-4; under
a score whose weights are off by 0.05, p_t drifts to 0.0024 by h = 1 toward a stationary 0.0057
(the live figure needs JavaScript)
The probability-flow ODE
Set the fresh-noise term to zero in the family above. What remains is a deterministic velocity field that carries the same distributions from noise to data. We can solve it with an ordinary ODE solver and, under suitable regularity, run it in either direction.
Definition 4.3.1 (probability-flow ODE). The member \(\lambda = 0\) of (4.2) is the probability-flow ODE,
\[\frac{dx}{dt} = u_t(x) := f_t\, x - \tfrac12\, g_t^2\, s_t(x) . \tag{4.3}\]Its solution map from time \(s\) to time \(t\) is written \(X_{s \to t}\); generation integrates from \(t = 1\) to \(t = 0\), so a sample is \(X_{1 \to 0}(z)\) with \(z \sim \mathcal N(0, I)\).
The velocity (4.3) is the one of Proposition 1.3.2, and Theorem 1.3.3 already said that its flow carries \(p_1\) to \(p_0\). What is new is a second form of the velocity, which does not mention the score at all.
For the next calculation, hold one pair \((x_0,z)\) fixed and let \(\alpha_t x_0+\sigma_t z\) trace a smooth curve. Its derivative is available from the schedule. This curve has the same one-time distributions as the forward SDE, but it is a different way of coupling times; Brownian sample paths do not have ordinary velocities. The proposition averages the velocities of these smooth curves.
Proposition 4.3.2 (the velocity is a conditional expectation). For the Gaussian kernel of Section 1,
\[\begin{aligned} u_t(x) &= \mathbb E\bigl[\dot x_t \,\big|\, x_t = x\bigr] = \dot\alpha_t\, D_t(x) + \dot\sigma_t\, \hat\varepsilon_t(x), \\ \dot x_t &:= \dot\alpha_t x_0 + \dot\sigma_t z . \end{aligned}\]Proof. By Tweedie's formula (2.2),
\[\begin{aligned} \dot\alpha_t D_t + \dot\sigma_t \hat\varepsilon_t &= \dot\alpha_t\, \frac{x + \sigma_t^2 s_t}{\alpha_t} - \dot\sigma_t \sigma_t s_t \\ &= \frac{\dot\alpha_t}{\alpha_t}\, x + \Bigl(\frac{\dot\alpha_t\sigma_t^2}{\alpha_t} - \dot\sigma_t\sigma_t\Bigr) s_t = f_t x - \tfrac12 g_t^2 s_t, \end{aligned}\]the last by \(g_t^2 = 2\sigma_t\dot\sigma_t - 2\sigma_t^2\dot\alpha_t/\alpha_t\) of Proposition 1.1.5. The first equality is the linearity of conditional expectation applied to \(\dot x_t = \dot\alpha_t x_0 + \dot\sigma_t z\). ∎
At a given point and time, the flow velocity is the conditional average of the velocities of the kernel paths \(t \mapsto \alpha_t x_0 + \sigma_t z\) passing through it. Section 6 takes this observation as the starting point for flow matching; here it follows from Tweedie's formula.
On the linear path, \(\dot\alpha = -1\) and \(\dot\sigma = 1\), so \(u_t = \hat\varepsilon_t - D_t = \mathbb E[z - x_0 \mid x_t]\). Each conditional path points from a data sample toward a noise sample. The flow follows the average direction consistent with its current position.
space-time of R2 under the variance-preserving schedule, noise at the left and data at the right, p_t shaded; through the probe (x, t) = (0.4, 0.4) the fan of fifteen kernel curves alpha_t x_0 + sigma_t z, one per mode and quantile of the posterior over x_0, whose weights are 0.314, 0.437 and 0.248 on the modes at −1.366, 0.204 and 1.585; the curves from those centres pass the probe with slopes 3.417, −0.060 and −3.120, whose posterior-weighted average 0.274 is u_t(x), the slope of the flow line through the probe; a Brownian path through the same point, whose chords over dt = 0.1, 0.01 and 0.001 are typically about 9.0, 28 and 90 in size (the live figure needs JavaScript)
Theorem 4.3.3 (properties of the flow). Suppose \(p_0\) has finite second moment and the velocity (4.3) is locally Lipschitz in \(x\), uniformly for \(t\) in compact subsets of \((0, 1]\); this holds for every Gaussian mixture. Then: (a) the ODE has a unique solution from every \((x, t)\), and \(X_{s \to t}\) is a diffeomorphism of \(\mathbb R^d\) with inverse \(X_{t \to s}\); (b) two solutions never occupy the same point at the same time; (c) \(X_{1 \to t}\) pushes \(p_1\) forward to \(p_t\) for every \(t \in (0, 1]\), and \(X_{t \to 0}\) pushes \(p_t\) to \(p_0\) when the velocity extends continuously to \(t = 0\).
Proof. (a) and (b) are the Picard–Lindelöf theorem and the uniqueness it gives; the inverse of the flow map is the flow map in the other direction. (c) is Theorem 1.3.3. ∎
Invertibility gives the deterministic flow a useful interpretation. Run (4.3) forward from a data point \(x\) to obtain \(z = X_{0 \to 1}(x)\). In the exact flow, running backward from that same noise point recovers \(x\). The noise serves as a latent code.
This is a continuous normalizing flow learned without solving an ODE during training. Song et al. use the map for image interpolation and likelihood evaluation. DDIM inversion, discussed in Section 5.1, approximates it with discrete steps.
left the density of 2,048 flow lines of the probability-flow ODE of R3 on the linear path, from
draws of N(0, I), forty-eight of them as lines, the particles at the chosen t over the exact
p_t, and a probe's line integrated to t = 1 with the image of a small disc along it
strip the particles at t = 1, 0.75, 0.5, 0.25 and 0 over the contours of p_t
right log p_t along the probe's line, exact and by the divergence integral of Theorem 4.4.1, with
the divergence itself
(the live figure needs JavaScript)
Example 4.3.4 (the flow of R3). The figure's two thousand flow lines run from a disc of noise to the two moons, drawn as a density with forty-eight of them as lines. Around \(t = 0.5\) the particles have gathered into a bent cloud that is \(p_{0.5}\); by \(t = 0.2\) the two moons are separate; at \(t = 0\) the particles lie on them. No two lines cross, by Theorem 4.3.3(b), although they come close where the flow splits between the moons: the lines that pass near the gap are the ones whose posterior, in the language of Section 2, was long undecided. ◁
Remark 4.3.5 (confusion: "deterministic" sampling is random). The flow is deterministic as a map; the sample \(X_{1 \to 0}(z)\) is random because \(z\) is. All its diversity comes from the initial noise. A stochastic sampler injects additional randomness, so its endpoint also depends on the later draws. How strongly it remembers its starting point depends on the process and the amount of noise. In the last figure of the introduction the reverse SDE's endpoints have a rank correlation of \(-0.009\) with their starts; the flow's is \(1.000\), which is one up to the solver's error.
Likelihood by the change of variables
An invertible flow changes density by expanding or compressing volume. Its divergence measures that local change. Integrating the divergence along a trajectory lets us recover a data point's log density from the density at the noise end.
Follow a tiny parcel of probability along the flow. If its volume doubles while its probability stays fixed, its density halves. Divergence measures the instantaneous relative expansion of that parcel: positive divergence lowers density along the path, and negative divergence raises it. Integrating this local accounting gives the change-of-variables formula below.
Theorem 4.4.1 (the instantaneous change of variables). Let \(x(t)\) solve (4.3) and let \(p_t\) be the marginals it transports. Then, along the curve and after integration,
\[\begin{aligned} \frac{d}{dt}\log p_t\bigl(x(t)\bigr) &= -\nabla \cdot u_t\bigl(x(t)\bigr), \\ \log p_0(x) &= \log p_1\bigl(X_{0 \to 1}(x)\bigr) + \int_0^1 \nabla \cdot u_t\bigl(X_{0 \to t}(x)\bigr)\, dt . \end{aligned}\]Proof. Along the curve, \(\frac{d}{dt}\log p_t(x(t)) = \frac{\partial_t p_t + u_t \cdot \nabla p_t}{p_t}\), and the continuity equation gives \(\partial_t p_t = -\nabla \cdot (p_t u_t) = -p_t \nabla \cdot u_t - u_t \cdot \nabla p_t\); the gradient terms cancel. Integrate from \(0\) to \(1\). ∎R. T. Q. Chen, Y. Rubanova, J. Bettencourt and D. Duvenaud, "Neural ordinary differential equations", NeurIPS 2018 (arXiv 1806.07366), Theorem 1; the trace estimator that makes the divergence affordable in high dimension is from W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever and D. Duvenaud, "FFJORD: free-form continuous dynamics for scalable reversible generative models", ICLR 2019 (arXiv 1810.01367), after M. F. Hutchinson (1990). Song et al. (2021), Section 4.3 and Appendix D.2, apply both to the probability-flow ODE and report \(2.99\) bits per dimension on CIFAR-10.
Computing the divergence exactly requires the trace of a \(d \times d\) Jacobian. Hutchinson's estimator replaces that trace by \(\mathbb E_v[v^\top \nabla u_t\, v]\), where the random vector \(v\) has identity covariance. Each estimate needs one vector–Jacobian product instead of \(d\).
The calculation therefore has two sources of numerical uncertainty: the ODE solver's error and the trace estimator's variance. It provides a likelihood through the flow's change of variables, complementing the variational bound in Theorem 3.4.1.
Example 4.4.2 (the probe of the figure). For the probe at \((0.20, 0.55)\) in the figure's default view, the flow carries the point to \(z = (0.65, 1.36)\), where \(\log \mathcal N(z; 0, I) = -2.971\); the integral of the divergence along the line is \(1.962\); their sum, \(-1.008\), is the estimate of \(\log p_0(x)\), against the exact value \(-1.007\) from the mixture formula. The same integral is the log determinant of the flow's Jacobian along the line, \(1.963\) by the variational equation, so the disc drawn along the probe's line grows in area by the factor \(e^{1.962} = 7.1\) on its way to noise. The right panel runs the identity at every intermediate time: the exact \(\log p_t(x(t))\) and its reconstruction from the integral lie on one curve. The divergence is positive along most of the line because, read from data to noise, the flow spreads a moon out into a disc. ◁
Remark 4.4.3 (warning: likelihood and sample quality disagree). The weights of Section 2.4 that give the best samples are not the ones that give the best likelihoods. Song, Durkan, Murray and Ermon showed that the denoising objective weighted by \(g_t^2\) is an upper bound on the negative log-likelihood of the probability flow, and that training with it improves the bits per dimension while it does not improve, and can worsen, the FID; Kingma and Gao's theorem of Section 2.4 explains why the two can disagree, since every weighted objective is a likelihood for a differently augmented law. Dupuis, Farghly, Haddouche, Durmus and Simsekli have since bounded the gap between the score-matching loss and the sample quality in three metrics and found that it is the low noise levels at which the score's accuracy matters most for closing it.Song, Durkan, Murray and Ermon (2021): "for a specific weighting scheme, the objective upper bounds the negative log-likelihood", with \(2.83\) and \(3.76\) bits per dimension on CIFAR-10 and ImageNet \(32 \times 32\). B. Dupuis, T. Farghly, M. Haddouche, A. Durmus and U. Simsekli, "Tightening the score matching gap for diffusion models", arXiv 2607.04442 (5 July 2026, a preprint when read on 24 October 2026): "the quality of the score approximation has more impact on closing the score matching gap for low noise scales".
Remark 4.4.4 (warning: a learned score need not be a score). The family of exact scores \(s_t\) satisfies a partial differential equation of its own, obtained by taking the gradient of the logarithm of the Fokker–Planck equation, and the velocity (4.3) built from it is a gradient field at each time, by Proposition 2.2.3. A network trained by denoising score matching satisfies neither property exactly; Lai et al. measured the violation and showed that penalizing it improves the likelihood and the conservativity of the learned field. The samplers of Section 5 and the theory of Section 10 assume a score that is close to a true score in \(L^2(p_t)\), and nothing more.C.-H. Lai, Y. Takida, N. Murata, T. Uesaka, Y. Mitsufuji and S. Ermon, "FP-Diffusion: improving score-based diffusion models by enforcing the underlying score Fokker-Planck equation", ICML 2023 (arXiv 2210.04296): "scores learned through denoising score matching (DSM) fail to fulfill the underlying score FPE", and satisfying it "improves the likelihood and the degree of conservativity".
Example 4.4.5 (the Jacobian of the learned denoiser). We compare the \(\varepsilon\)-model of Section 11, widths (64, 64, 128), trained on FashionMNIST at \(28 \times 28\) on the 150-step chain, \(\beta_i\) linear from \(10^{-4}\) to \(0.02\) (Example 3.1.3), and the same architecture trained on DDPM's chain of 1,000 steps. Four held-out images are noised with fixed draws to six indices of each chain. At every such point, the Jacobian of the learned denoiser \(D_\theta = (x_t - \sigma_i \varepsilon_\theta)/\alpha_i\), computed by differentiating the network, is set against \((\alpha_i/\sigma_i^2)\, \mathrm{Cov}(x_0 \mid x_t)\) under the empirical law of the 60,000 training images, the value the Jacobian would have if the network had learned the training set's score exactly.
The figure compares these two operators, with the six leading eigen-images of both operators side by side and the tile \(D_\theta - D_{\mathrm{emp}}\), where the network's answer leaves the empirical posterior mean. The learned denoiser's Jacobian is not symmetric: the relative asymmetry \(\|J - J^\top\|/\|J\|\) is \(0.147\) at the clean end of the 150-step chain and \(0.697\) at its noise end, so the learned field is not a gradient and Proposition 2.2.3's identity cannot hold for it exactly.
Where the empirical law's posterior is owned by one training image, at low noise, its exact covariance is zero and the learned Jacobian's top eigenvalue is \(1.397\): the network is not the lookup table of Section 2.5, which would be flat there. Where the posterior straddles a handful of training images, at the noise end of the 150-step chain with \(6.5\) of them in effect, the exact covariance's top eigenvalue is \(4.657\) and the learned Jacobian's \(2.358\), the network the smoother of the two.
At the noise end of DDPM's chain of 1,000 steps the exact posterior averages \(59{,}722\) training images and its covariance has top eigenvalue \(0.515\), while the learned Jacobian's is \(21.97\): at the noise levels where the whole training set is visible the network is more expansive than the law it was trained on, which is where the samplers of Section 5 start. The exact counterpart on R3 is in Example 2.2.4: at its probe the Jacobian is symmetric to the last digit and its larger eigenvalue is \(1.386\). ◁
for a held-out image noised to a chain index: the learned denoiser's output, the exact posterior mean over the training set and their difference D_theta - D_emp as a change map; the six leading eigen-images and eigenvalues of the symmetrized Jacobian of the learned denoiser and of the exact posterior covariance (alpha/sigma^2) Cov(x_0 | x_t) of the empirical law, one eigen-image followed across the index; the Jacobian's asymmetry, from 0.1-0.2 at low noise to 0.6-0.8 at high noise, and the effective number of training points in the exact posterior, from 1 to tens of thousands; the default view, image 1 at index 119 of the 150-step chain: top eigenvalues 1.764 learned and 6.162 exact, asymmetry 0.564, 11.4 training images in effect, D_theta - D_emp of root mean square 0.111 (the live figure needs JavaScript and loads about 0.86 MB of data)
Discretization, and the shape of what is proved
The continuous equations assume an exact score and an exact solution. A working sampler has neither. It also usually starts from a Gaussian approximation to the terminal law. Keeping these three errors separate helps us see what better training can fix and what requires a better sampler. Algorithm 4.2.4 gives the simplest discretization; Section 5 improves it, and Section 10.1 states the relevant guarantees and their assumptions.
These distinctions suggest three diagnostic experiments. Start from the true terminal law to isolate initialization, replace the learned field by the exact one in a toy example to isolate learning, and refine the time grid to expose numerical error. A change that helps only one experiment identifies a particular cause. Increasing the step count alone cannot establish that the learned field or the starting law is correct.
Remark 4.5.1 (takeaway: three errors). Consider a sampler built from (4.2) with a learned score \(s_\theta\), initialized at \(\gamma = \mathcal N(0, I)\) and run for \(N\) steps. Its error has three sources.
The first is the starting distribution. The prior error, \(\mathrm{KL}(p_1 \| \gamma)\), is bounded by Theorem 1.4.3 and enters the final bound through data processing. The examples also measure the conditional prior term: on the 150-step chain it is \(70.05\) nats per image (Example 3.1.5). Example 5.4.2 shows the effect of the starting mismatch on samples.
The second is the learned score, measured by \(\int g_t^2\, \mathbb E_{p_t}\|s_\theta - s_t\|^2\, dt\). Section 2.3 relates this quantity to denoising loss above its irreducible floor. The third is discretization: the error from following the field on a finite grid, which depends on the solver and the field's regularity.
For the reverse SDE, Chen, Chewi, Li, Li, Salim and Zhang proved that these add up to a polynomial bound with no assumption on \(p_0\) beyond moments, Benton, De Bortoli, Doucet and Deligiannidis sharpened the step count to \(\tilde O(d\log^2(1/\delta)/\varepsilon^2)\) for an accuracy \(\varepsilon^2\) in relative entropy on a law smoothed by noise of variance \(\delta\), and Li, Wei, Chen and Chi showed that the probability-flow ODE converges at rate \(1/N\) where the stochastic sampler converges at \(1/\sqrt N\).S. Chen, S. Chewi, J. Li, Y. Li, A. Salim and A. R. Zhang, "Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions", ICLR 2023 (arXiv 2209.11215); J. Benton, V. De Bortoli, A. Doucet and G. Deligiannidis, "Nearly \(d\)-linear convergence bounds for diffusion models via stochastic localization", ICLR 2024 (arXiv 2308.03686); G. Li, Y. Wei, Y. Chen and Y. Chi, "Towards non-asymptotic convergence for diffusion-based generative models", ICLR 2024 (arXiv 2306.09251): "a convergence rate proportional to \(1/T\)" for the ODE sampler and "\(1/\sqrt T\)" for the DDPM-type sampler, "no smoothness assumption is imposed". Section 10.1 tabulates these with the 2024 to 2026 refinements. The decomposition separates the learning problem from the sampling problem. Once the weighted \(L^2\) score error is controlled, the remaining questions concern initialization and numerical approximation.