Score matching, the reverse-time SDE and the probability-flow ODE

Before this lesson: SDEs, Euler–Maruyama and the Fokker–Planck equation (lesson 12); the multivariate Gaussian and Gaussian conditioning (probability 05). For the machine-learning side of diffusion models (networks, training in practice), see generative models 06, 07 and 09; this lesson focuses on the mathematics and does not repeat them. Time: about 110 minutes, plus the notebook.

---

Why this matters #

• Ensemble weather forecasts. GenCast (Price et al., Nature 2025) produces each ensemble member of a 15-day global forecast as one sample of a conditional diffusion model, conditioned on the two most recent weather states. The cost of an ensemble is members × sampler steps × network cost, so the choice of sampler matters.
• Downscaling and super-resolution. Conditional diffusion models generate fine-scale fields (precipitation, wind) that are consistent with a coarse simulation, with realistic small-scale variability instead of a blurred mean.
• Priors for inverse problems and data assimilation. By Bayes' rule, the score of a posterior is the score of the prior plus the score of the likelihood. A learned score is a prior for the state of the atmosphere or the subsurface (score-based data assimilation, Rozet and Louppe 2023).
• Seismic waveforms. Denoising diffusion models generate synthetic earthquake records conditioned on magnitude, distance and site (Palgunadi et al. 2024; Jung et al. 2025).
• Every practical choice in these systems (SDE or ODE, which solver, how many steps, how to weight the loss over noise levels) follows from the three theorems of section 3.

1. Predict first #

2. The picture #

The score points toward high density #

The score of a density $p$ on $\R^d$ is the vector field

A Gaussian example. For $p = \mathcal N(\mu, v)$ in 1D, $\log p(x) = -(x - \mu)^2/(2v) + \text{const}$, so $s(x) = -(x - \mu)/v$. The score points toward the mean, and it is stronger when the variance is smaller. The constant disappears under the derivative.

This is true in general. If $p = \tilde p/Z$ with an unknown normalising constant $Z$, then $\log p = \log\tilde p - \log Z$, and $\log Z$ does not depend on $x$. So $\nabla_x\log p = \nabla_x\log\tilde p$. A score model never needs $Z$. This removes the main difficulty of energy-based models (generative models 06).

Noise smooths the data #

Add Gaussian noise of level $\sigma$: $x_\sigma = x_0 + \sigma\varepsilon$, with $x_0 \sim p_0$ (the data) and $\varepsilon \sim \mathcal N(0, I)$. The density of $x_\sigma$ is the convolution $p_\sigma = p_0 * \mathcal N(0, \sigma^2I)$: every data point is smeared into a small Gaussian cloud.

• For small $\sigma$, $p_\sigma$ is close to the data distribution.
• For large $\sigma$, $p_\sigma$ is close to a wide Gaussian $\mathcal N(0, \sigma^2 I)$, which is easy to sample.

A diffusion model walks from large $\sigma$ down to small $\sigma$ and uses the score of $p_\sigma$ at every level. For a Gaussian mixture $p_0 = \sum_k w_k\mathcal N(m_k, C_k)$ everything is exact: noising keeps a mixture, $p_\sigma = \sum_k w_k\,\mathcal N(m_k, C_k + \sigma^2 I)$, and its score is a weighted average of Gaussian scores,

The weights $r_k(x)$ are the responsibilities: the probability that $x$ came from component $k$. All figures and the notebook of this lesson use such closed forms instead of a trained network.

Interactive figure: loading…

The score is a denoiser #

Section 3 proves Tweedie's formula:

The left side is the best guess (in mean squared error) of the clean sample, given the noisy one. So the score at a noisy point, times $\sigma^2$, is the step from that point to the average clean point it came from. Learning to denoise and learning the score are the same task.

Interactive figure: loading…

Sampling with the score: Langevin dynamics #

Lesson 12 showed that the SDE $dx = \nabla\log\pi(x)\,dt + \sqrt2\,dW$ has $\pi$ as its stationary density. Its Euler–Maruyama discretisation is the unadjusted Langevin algorithm (ULA):

So a score alone is enough to sample. But Langevin dynamics with the score of the clean data has two problems. Between separated modes the density is almost zero, the score carries no information there, and a chain rarely crosses from one mode to another. And the chain needs a very long time to reach the right weights of the modes. Smoothing the target with noise fixes both problems, at the price of sampling a blurred distribution. The solution of diffusion models is to use all noise levels in sequence, from large to small.

Interactive figure: loading…

3. Definitions and main results #

Score matching #

We want a model $s_\theta(x)$ of the score of $p_\sigma$. The natural loss is explicit score matching (ESM):

We cannot compute it, because the target $\nabla\log p_\sigma$ is unknown. Hyvärinen (2005) removed the target by integration by parts; the result needs the trace of the Jacobian of $s_\theta$, which costs $d$ backward passes in dimension $d$. Denoising score matching avoids both problems.

Theorem 1 (denoising score matching; Vincent 2011) #

Let $q(x \mid x_0) = \mathcal N(x; x_0, \sigma^2I)$ be the noising kernel. Its score is $\nabla_x\log q(x \mid x_0) = -(x - x_0)/\sigma^2 = -\varepsilon/\sigma$. Define the denoising score matching loss

If $E\|s_\theta(x)\|^2 < \infty$ for $x \sim p_\sigma$, then $J_{\text{DSM}}(\theta) = J_{\text{ESM}}(\theta) + C$, where $C$ does not depend on $\theta$.

Why this is true (intuition). A noisy point $x$ could have come from many clean points $x_0$. Each one gives a different target $-\varepsilon/\sigma$. The model cannot know which $x_0$ was used, so the best it can do is predict the average target at $x$. That average is exactly the score of $p_\sigma$ (step 3 below).

Proof

1. Expand both squares: $\|s - y\|^2 = \|s\|^2 - 2\langle s, y\rangle + \|y\|^2$, with the target $y$.
2. The terms $E\|s_\theta(x)\|^2$ are equal in the two losses, because in both $x$ has the law $p_\sigma$.
3. Key identity. Differentiate $p_\sigma(x) = \int p_0(x_0)\,q(x \mid x_0)\,dx_0$ under the integral and use $\nabla q = q\,\nabla\log q$: $\nabla p_\sigma(x) = \int p_0(x_0)\,q(x \mid x_0)\,\nabla_x\log q(x \mid x_0)\,dx_0$. Divide by $p_\sigma(x)$; the factor $p_0(x_0)q(x \mid x_0)/p_\sigma(x)$ is the posterior $p(x_0 \mid x)$ by Bayes' rule:

$$\nabla\log p_\sigma(x) = E\big[\nabla_x\log q(x \mid x_0)\,\big|\,x\big].$$

The marginal score is the posterior mean of the conditional score. 4. The cross terms are equal. By step 3 and the tower property $E[E[\,\cdot \mid x]] = E[\,\cdot\,]$: $E\langle s_\theta(x), \nabla\log p_\sigma(x)\rangle = E\langle s_\theta(x), E[\nabla_x\log q \mid x]\rangle = E\langle s_\theta(x), \nabla_x\log q(x \mid x_0)\rangle$. 5. The remaining terms do not contain $\theta$, so $C = E\|\varepsilon/\sigma\|^2 - E\|\nabla\log p_\sigma(x)\|^2 = d/\sigma^2 - E\|\nabla\log p_\sigma\|^2$. $\blacksquare$

Consequences.

• Both losses have the same gradient in $\theta$, so DSM trains toward the true score of $p_\sigma$. Over all functions, the minimiser is the conditional mean $s^*(x) = E[-\varepsilon/\sigma \mid x] = \nabla\log p_\sigma(x)$.
• The loss does not go to zero. Its minimum is $C$, the conditional variance of the target given $x$: $C = E\,\|\varepsilon/\sigma - E[\varepsilon/\sigma \mid x]\|^2 \ge 0$. No model can predict the part of the noise that $x$ does not reveal.
• Tweedie's formula. From $x_0 = x - \sigma\varepsilon$ and step 3, $E[x_0 \mid x] = x - \sigma E[\varepsilon \mid x]$, and this equals $x + \sigma^2\nabla\log p_\sigma(x)$. One derivative more gives the posterior covariance, $\operatorname{Cov}[x_0 \mid x] = \sigma^2\big(I + \sigma^2\nabla^2\log p_\sigma(x)\big)$.
• All noise levels at once. Train one network $s_\theta(x, \sigma)$, draw $\sigma$ at random (for example log-uniform), and weight the loss by $\lambda(\sigma) = \sigma^2$: $\sigma^2\|s_\theta + \varepsilon/\sigma\|^2 = \|\sigma s_\theta + \varepsilon\|^2$. This is $\varepsilon$-prediction: a network $\hat\varepsilon_\theta$ predicts the noise and $s_\theta = -\hat\varepsilon_\theta/\sigma$. The DDPM loss is the same loss for the kernel $x = \sqrt{\bar\alpha}\,x_0 + \sqrt{1 - \bar\alpha}\,\varepsilon$ (generative models 08 compares the parameterisations).

Noising SDEs and their marginals #

A diffusion model adds noise continuously with a forward SDE

with a drift $f$, a scalar diffusion coefficient $g(t)$, and marginal densities $p_t$. Two choices are standard (Song et al. 2021):

• Variance exploding (VE): $f = 0$. Then $x_t = x_0 + \sigma(t)\varepsilon$ with $\sigma(t)^2 = \int_0^t g(s)^2\,ds$. In this lesson, as in the EDM convention of Karras et al. (2022), $g(t) = \sqrt{2t}$, so $\sigma(t)^2 = t^2$ and the noise level equals the time: $\sigma = t$.
• Variance preserving (VP): $f = -\tfrac12\beta(t)x$, $g = \sqrt{\beta(t)}$. This is an OU process with a time-dependent rate (lesson 12): $x_t = \alpha_t x_0 + \sqrt{1 - \alpha_t^2}\,\varepsilon$ with $\alpha_t = \exp\big(-\tfrac12\int_0^t\beta(s)\,ds\big)$.

In both cases $p_t$ is the data distribution scaled and convolved with a Gaussian, so a training sample $x_t$ is drawn in one step and the DSM loss applies at every $t$.

Theorem 2 (reverse-time SDE; Anderson 1982) #

Let $x$ solve the forward SDE on $[0, T]$ with marginals $p_t$. Assume $f$, $g$ and $p_t$ are smooth, $p_t > 0$, and the Fokker–Planck equation has a unique solution for a given initial density. In the reversed clock $\tau = T - t$, the process

has $y(\tau) \sim p_{T-\tau}$ for every $\tau$. Written in the original clock (so $dt < 0$), this is the familiar form

Why the extra term. Running the forward SDE backward would undo the drift but would add more noise, not remove it. The term $g^2\nabla\log p_t$ pushes every particle up the density, toward the data, strongly enough to undo the spreading of the noise and to compensate the fresh noise of the reverse process.

Proof (through the Fokker–Planck equation)

1. In $d$ dimensions, with noise $g(t)\,dW$, the forward marginals solve $\partial_t p = -\nabla\cdot(f p) + \tfrac12 g^2\Delta p$ (lesson 12, the $d$-dimensional form).
2. Let $q_\tau = p_{T-\tau}$. Then $\partial_\tau q = -\partial_t p = \nabla\cdot(f q) - \tfrac12 g^2\Delta q$.
3. A candidate SDE $dy = b\,d\tau + g\,dW_\tau$ has the Fokker–Planck equation $\partial_\tau q = -\nabla\cdot(b q) + \tfrac12 g^2\Delta q$.
4. Make the two right-hand sides equal: we need $-\nabla\cdot(bq) = \nabla\cdot(fq) - g^2\Delta q$. Use $\Delta q = \nabla\cdot(\nabla q) = \nabla\cdot(q\,\nabla\log q)$. Then the right side is $\nabla\cdot\big((f - g^2\nabla\log q)\,q\big)$, and $b = -f + g^2\nabla\log q$ works.
5. Both processes start from $p_T$ and solve the same Fokker–Planck equation, so by uniqueness they have the same marginals. $\blacksquare$

This proof shows that the marginals agree, which is all a sampler needs. Anderson proved more: the law of the whole reversed path equals the law of the forward path run backward.

Euler–Maruyama for the reverse SDE. For the VE SDE of this lesson ($f = 0$, $g^2 = 2t$), on a decreasing grid $t_0 > t_1 > \dots > t_N$ with $\Delta_n = t_n - t_{n+1} > 0$ and $z_n \sim \mathcal N(0, I)$:

Interactive figure: loading…

Theorem 3 (probability-flow ODE; Song et al. 2021) #

The ODE

has the same marginals $p_t$ as the forward SDE, in either direction of time.

Proof

1. Write the diffusion term as a flux: $\tfrac12 g^2\Delta p = \nabla\cdot\big(\tfrac12 g^2\,p\,\nabla\log p\big)$, because $p\nabla\log p = \nabla p$.
2. The Fokker–Planck equation becomes $\partial_t p = -\nabla\cdot\big((f - \tfrac12 g^2\nabla\log p)\,p\big)$. This is a continuity equation with the velocity $u = f - \tfrac12 g^2\nabla\log p$, and no diffusion term.
3. A continuity equation $\partial_t p + \nabla\cdot(p\,u) = 0$ is exactly the density equation of the ODE $dx/dt = u(x, t)$ (lesson 12 with $b = 0$). Same equation, same initial density, so the same marginals. $\blacksquare$

For the VE SDE of this lesson the PF-ODE is $dx/dt = -t\,s(x, t)$. It has no noise, so a higher-order solver works well. Heun's method (two score evaluations per step) has global error $O(h^2)$: with $F(x, t) = -t\,s(x, t)$ and a step $h_n = t_{n+1} - t_n$ (negative when sampling),

Interactive figure: loading…

A family of samplers. The same flux trick shows that, for every $\tilde g(t) \ge 0$, the backward process $dx = \big[f - \tfrac12\big(g^2 + \tilde g^2\big)\nabla\log p_t\big]\,dt + \tilde g\,d\bar W$ has the marginals $p_t$. With $\tilde g = g$ you get the reverse SDE; with $\tilde g = 0$ you get the PF-ODE. The amount of noise is a free design choice.

ODE or SDE?

• The ODE is a deterministic bijection between noise and data. It gives encodings, and exact log-likelihoods through the instantaneous change of variables $\frac{d}{dt}\log p_t(x(t)) = -\nabla\cdot u$ (as for neural ODEs, lesson 03). With Heun and a good time grid it needs few steps (tens).
• The SDE adds fresh noise. The noise and the score drift pull samples back toward $p_t$, so earlier errors (score errors, a wrong start distribution) are partly forgotten. It usually needs more steps; plain Euler–Maruyama is the usual solver.
• The time grid is geometric in $\sigma$ (a constant ratio $t_{n+1}/t_n$), because the score changes on the scale of $\sigma$ itself.
• The start. You cannot sample $p_T$ exactly. You start from $\mathcal N(0, T^2I)$, which is close to $p_T$ when $T$ is much larger than the data scale (the prior mismatch error).

4. Worked example #

(a) A Gaussian. Data $p_0 = \mathcal N(2, 1)$ in 1D, noise level $\sigma = 1$.

1. Noised density: variances add, so $p_\sigma = \mathcal N(2, 1 + 1) = \mathcal N(2, 2)$.
2. Score at $x = 4$: $s(4) = -(4 - 2)/2 = -1$.
3. Tweedie: $E[x_0 \mid x = 4] = 4 + \sigma^2 s(4)$, that is $4 + 1\cdot(-1) = 3$.
4. Check with Gaussian conditioning: prior variance 1, noise variance 1, so the posterior mean is $(1\cdot 4 + 1\cdot 2)/(1 + 1) = 3$. ✓
5. Posterior variance from Tweedie's second-order formula: $s'(x) = -\tfrac12$, so $\sigma^2(1 + \sigma^2 s') = 1\cdot(1 - \tfrac12) = \tfrac12$. Gaussian conditioning: $1\cdot 1/(1 + 1) = \tfrac12$. ✓
6. DSM floor: $C = d/\sigma^2 - E[s(x)^2]$. Here $E[s(x)^2] = E[(x - 2)^2]/4 = 2/4$, so $C = 1 - \tfrac12 = \tfrac12$. Even the exact score has $J_{\text{DSM}} = 0.5$.
7. PF-ODE with $\sigma = t$: $p_t = \mathcal N(2, 1 + t^2)$ and $dx/dt = -t\,s = t(x - 2)/(1 + t^2)$. Separate variables: $\frac{d(x - 2)}{x - 2} = \frac{t\,dt}{1 + t^2}$, so $\log|x - 2| = \tfrac12\log(1 + t^2) + \text{const}$ and $x(t) - 2 = (x(T) - 2)\sqrt{(1 + t^2)/(1 + T^2)}$.
8. From $T = 10$ and $x(T) = 12$: $x(0) = 2 + 10/\sqrt{101} = 2.9950$. One noise value gives one data value: this is the deterministic decoder.

(b) A mixture. $p_0 = \tfrac12\mathcal N(-2, 0.25) + \tfrac12\mathcal N(2, 0.25)$ and $\sigma = 0.5$.

1. Noising adds $\sigma^2 = 0.25$ to each variance: $p_\sigma = \tfrac12\mathcal N(-2, 0.5) + \tfrac12\mathcal N(2, 0.5)$.
2. At $x = 1$, the two Gaussian scores are $-(1 - 2)/0.5 = 2$ and $-(1 + 2)/0.5 = -6$.
3. The responsibility of the right component: the log-densities differ by $\big[(1 + 2)^2 - (1 - 2)^2\big]/(2\cdot 0.5)$, that is $(9 - 1)/1 = 8$. So $r_+ = 1/(1 + e^{-8}) = 0.999665$ and $r_- = 0.000335$.
4. Score: $0.999665\cdot 2 + 0.000335\cdot(-6) = 1.99732$.
5. Tweedie: $E[x_0 \mid x = 1] = 1 + 0.25\cdot 1.99732 = 1.4993$. This is almost the posterior mean inside the right component alone, $(0.25\cdot 1 + 0.25\cdot 2)/0.5 = 1.5$.

Far from both components, both Gaussian densities underflow to 0 in floating point (at $x = 1000$ each is about $e^{-10^6}$), and $r_\pm$ becomes $0/0$. Compute the responsibilities from the log-densities with a log-sum-exp (a softmax in log space).

5. In code #

The mixture score in closed form against a finite difference of $\log p_\sigma$, Tweedie's formula, Theorem 1 by Monte Carlo, the PF-ODE map of the worked example, and a check of the factor $\tfrac12$:

import numpy as np
from scipy.special import logsumexp
from scipy.integrate import solve_ivp

w, m, v0 = np.array([0.5, 0.5]), np.array([-2.0, 2.0]), 0.25   # p_0 = 0.5 N(-2, 0.25) + 0.5 N(2, 0.25)

def log_p(x, sigma):
    """log p_sigma(x): each component variance grows by sigma^2."""
    V = v0 + sigma**2
    return logsumexp(np.log(w) - 0.5 * np.log(2 * np.pi * V) - (x - m)**2 / (2 * V))

def score(x, sigma):
    """Closed form: responsibility-weighted Gaussian scores, responsibilities in log space."""
    V = v0 + sigma**2
    lw = np.log(w) - (x - m)**2 / (2 * V)
    r = np.exp(lw - logsumexp(lw))
    return np.sum(r * (-(x - m) / V))

sigma, x = 0.5, 1.0
fd = (log_p(x + 1e-5, sigma) - log_p(x - 1e-5, sigma)) / 2e-5
print(f"{score(x, sigma):.5f} {fd:.5f}")          # 1.99732 1.99732
print(f"{x + sigma**2 * score(x, sigma):.4f}")    # Tweedie: E[x0 | x = 1] = 1.4993

# Theorem 1 by Monte Carlo: p_0 = N(2, 1), sigma = 1. For every linear model s(x) = a x + b,
# J_DSM - J_ESM is the same constant 1/sigma^2 - 1/(1 + sigma^2) = 0.5.
rng = np.random.default_rng(0)
x0 = 2 + rng.standard_normal(2_000_000)
eps = rng.standard_normal(2_000_000)
xs = x0 + eps                                          # samples of p_sigma = N(2, 2)
for a, b in [(0.0, 0.0), (-0.5, 1.0), (1.0, -3.0)]:
    s = a * xs + b
    J_dsm = np.mean((s + eps / 1.0)**2)                # target -eps / sigma
    J_esm = np.mean((s + (xs - 2) / 2)**2)             # target: the true score of N(2, 2)
    print(f"a = {a:4.1f}, b = {b:4.1f}: J_DSM - J_ESM = {J_dsm - J_esm:.3f}")
# a =  0.0, b =  0.0: J_DSM - J_ESM = 0.500
# a = -0.5, b =  1.0: J_DSM - J_ESM = 0.500
# a =  1.0, b = -3.0: J_DSM - J_ESM = 0.502

# PF-ODE dx/dt = -t * score for p_0 = N(2, 1) under dx = sqrt(2t) dW, so p_t = N(2, 1 + t^2)
sol = solve_ivp(lambda t, y: t * (y - 2) / (1 + t**2), (10.0, 0.0), [12.0], rtol=1e-10, atol=1e-12)
print(f"{sol.y[0, -1]:.4f} {2 + 10 / np.sqrt(101):.4f}")   # 2.9950 2.9950

The third Monte Carlo value differs from 0.5 only by sampling error. Now the factor $\tfrac12$. In the reversed clock $\tau = T - t$, a sampler $dy = \kappa(t)\,(y - 2)\,d\tau + c\,dW$ changes the variance by $dV/d\tau = 2\kappa V + c^2$. The right variance is that of $p_t$, $V = 1 + t^2$, so $dV/d\tau = -2t$:

t = np.linspace(0.1, 5, 50)
V, dV_dtau = 1 + t**2, -2 * t                     # tau = T - t runs backward
g2 = 2 * t
kappa = -g2 / V                                   # g^2 * score = kappa * (y - 2)
for name, kap, c2 in [("reverse SDE", kappa, g2), ("PF-ODE", kappa / 2, 0 * t), ("PF-ODE without 1/2", kappa, 0 * t)]:
    print(f"{name:20s} max residual {np.max(np.abs(2 * kap * V + c2 - dV_dtau)):.1e}")
# reverse SDE          max residual 1.8e-15
# PF-ODE               max residual 8.9e-16
# PF-ODE without 1/2   max residual 1.0e+01

Both correct samplers keep the variance of $p_t$ exactly. Without the $\tfrac12$ the residual is $-2t$: the variance shrinks too fast, and the samples come out too narrow.

6. Where it shows up #

• Diffusion models. DDPM, score SDEs and EDM are all DSM plus a sampler for the reverse SDE or the PF-ODE. Fast samplers such as DDIM and DPM-Solver are ODE solvers for the PF-ODE; DPM-Solver is an exponential integrator that treats the linear part exactly (the idea from lesson 02; details in generative models 09).
• Weather and climate. GenCast (Price et al., Nature 2025) samples each ensemble member with a conditional diffusion model on the sphere, conditioned on the two most recent atmospheric states, and rolls forward in 12-hour steps. As of 2026, diffusion-based ensembles are one of the main families of ML weather models (lesson 16).
• Seismology. Denoising diffusion models generate three-component earthquake waveforms conditioned on magnitude, distance and site parameters (Palgunadi et al. 2024; Jung et al. 2025). The loss is the DSM loss of this lesson with the conditions as extra inputs: $s_\theta(x, \sigma \mid \text{magnitude}, \text{distance}, \dots)$.
• Data assimilation and inverse problems. By Bayes' rule,

the normaliser $p(y)$ drops out. A learned score is a prior for the state; the likelihood score comes from the observation model. Score-based data assimilation (Rozet and Louppe 2023) and diffusion posterior sampling (Chung et al. 2023) use Tweedie's formula to approximate $\nabla\log p_t(y \mid x_t)$ cheaply. Compare with the classical filters of lesson 06 and the inverse problems of lesson 04.

• Physics. For a Boltzmann density $p \propto e^{-U/(k_BT)}$, the score is the force divided by $k_BT$. Langevin dynamics $dx = \nabla\log p\,dt + \sqrt2\,dW$ has $p$ as stationary density (lesson 12). It is the sampler family above with $f = 0$, $g = 0$ (so $p_t = p$ does not change) and $\tilde g = \sqrt2$.

7. Common mistakes #

• Wrong sign or wrong clock in the reverse SDE. In the original clock the drift is $f - g^2\nabla\log p_t$ and $dt < 0$; in the reversed clock it is $-f + g^2\nabla\log p$ with $d\tau > 0$. Mixing the two gives a process that adds noise instead of removing it.
• Dropping the $\tfrac12$ in the PF-ODE. With drift $f - g^2\nabla\log p$ and no noise, the samples contract too fast and come out too narrow (section 5).
• Using the loss value as a quality measure. The DSM loss has a floor $C(\sigma)$ that depends on $\sigma$ and on the data. Compare against a closed-form score when you can.
• One noise level, or $\sigma$ uniform on a linear scale. Most samples then sit at large $\sigma$. Use a log-uniform (or log-normal) $\sigma$ and a geometric sampling grid.
• Dividing by zero at $\sigma \to 0$. The score scales like $1/\sigma$. Stop the sampler at a small $\sigma_{\min} > 0$ and let the network predict $\varepsilon$ (unit scale) instead of the score.
• Mixture responsibilities in linear space. They underflow far from the components. Use log-sum-exp.
• Judging samplers with small samples. A distance between two independent samples of the same distribution (for example the energy distance) is not zero. Compare against that noise floor.

8. Check yourself #

1. Why does a score model not need the normalising constant of $p$?
   Answer

   If $p = \tilde p/Z$, then $\log p = \log\tilde p - \log Z$, and $\log Z$ does not depend on $x$. So $\nabla_x\log p = \nabla_x\log\tilde p$.
2. Data $p_0 = \mathcal N(0, 4)$ in 1D, noise level $\sigma = 3$, weight 1. What is the minimum value of $J_{\text{DSM}}$?
   Answer

   $C = 1/\sigma^2 - E[s^2] = 1/\sigma^2 - 1/(s_0^2 + \sigma^2) = 1/9 - 1/13 = 4/117 \approx 0.0342$, where $s_0^2 = 4$ is the data variance. It is reached by the exact score of $p_\sigma = \mathcal N(0, 13)$.
3. Write the probability-flow ODE of the variance-preserving SDE $dx = -\tfrac12\beta(t)\,x\,dt + \sqrt{\beta(t)}\,dW$.
   Answer

   $f = -\tfrac12\beta x$ and $g^2 = \beta$, so $dx/dt = -\tfrac12\beta(t)\big(x + \nabla\log p_t(x)\big)$.
4. A network $\hat\varepsilon_\theta(x, \sigma)$ predicts the noise. Give the score and the denoiser $D(x) = E[x_0 \mid x]$ in terms of $\hat\varepsilon_\theta$.
   Answer

   $s_\theta = -\hat\varepsilon_\theta/\sigma$, and by Tweedie $D(x) = x + \sigma^2 s_\theta = x - \sigma\,\hat\varepsilon_\theta(x, \sigma)$.
5. In step 4 of the proof of Theorem 2, why does $b = -f + g^2\nabla\log q$ follow?
   Answer

   We need $-\nabla\cdot(bq) = \nabla\cdot(fq) - g^2\Delta q$. Since $\Delta q = \nabla\cdot(q\nabla\log q)$, the right side is $\nabla\cdot\big((f - g^2\nabla\log q)\,q\big)$, so $bq = -(f - g^2\nabla\log q)\,q$ works.
6. Why can the PF-ODE give exact log-likelihoods while the reverse SDE cannot do so directly?
   Answer

   The ODE flow is a bijection, so the instantaneous change of variables applies: $\log p_0(x(0)) = \log p_T(x(T)) + \int_0^T\nabla\cdot u(x(t), t)\,dt$. The SDE map is random, not a bijection, so there is no such formula.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Define the score of a density $p$.

   Show answer

2. Why does a score model not need the normalising constant $Z$ of $p = \tilde p/Z$?

   Show answer

3. Write the denoising score matching loss at noise level $\sigma$.

   Show answer

4. State Tweedie's formula.

   Show answer

5. Why does DSM learn the marginal score although its target is the conditional score $-\varepsilon/\sigma$?

   Show answer

6. Write the reverse-time SDE of $dx = f\,dt + g\,dW$.

   Show answer

7. What is the drift of the probability-flow ODE?

   Show answer

8. VE and VP noising SDEs: give $f$ and $g$.

   Show answer

9. When would you sample with the PF-ODE instead of the reverse SDE?

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "What does a diffusion model actually learn?" The score of the noised data at every noise level, by denoising. Equivalently the noise $\varepsilon$, or the denoiser $E[x_0 \mid x_\sigma]$; these are one function seen through linear maps.
• "Why add noise at all instead of doing score matching on clean data?" Data lie near a low-dimensional set, where the clean score is undefined or badly estimated away from the data. Noise fills the space and connects the modes, so sampling can move between them.
• "You have a budget of about 20 network evaluations. Which sampler?" The PF-ODE with a second-order solver (Heun) on a geometric (EDM) time grid; the SDE needs more steps.
• "Why does DSM work without the true score?" DSM = ESM + a constant (Theorem 1): the cross terms match because the marginal score is the posterior mean of the conditional score.
• "What is Tweedie's formula and where is it used?" $E[x_0 \mid x] = x + \sigma^2\nabla\log p_\sigma(x)$; in denoisers, and in likelihood approximations for posterior sampling in inverse problems and data assimilation.

10. Practice #

The notebook for this lesson runs in Colab or in the browser. It has five exercises, each with a check cell and a folded solution:

1. The closed-form score of a noised 2D Gaussian mixture (responsibilities with log-sum-exp), checked against finite differences of $\log p_\sigma$.
2. Tweedie's denoiser; check it on the worked example and against Gaussian conditioning.
3. Theorem 1 by Monte Carlo: the DSM and ESM losses differ by the same constant for several models.
4. Two samplers with the exact score: the reverse SDE (Euler–Maruyama) and the PF-ODE (Heun); compare their samples with the data by the energy distance and its noise floor.
5. Learn the score from samples: fit a model that is linear in its parameters (random Fourier features) to the DSM targets by least squares, then sample with it.

11. Watch and read #

Diffusion and Score-Based Generative ModelsMITCBMM · 92 min

Yang Song explains score matching, the noise levels, the reverse SDE and the probability-flow ODE.

Score-based Diffusion Models | Generative AI AnimatedDeepia · 19 min

An animated overview of score-based diffusion models and the SDE view.

Score Matching Explained - The Key Idea Behind Diffusion ModelsDeepBean · 28 min

Why score matching works without the true score, the idea behind Theorem 1.

MIT 6.S184: Flow Matching and Diffusion Models - Lecture 03A - Score Functions (2026)Peter Holderrieth · 46 min

Conditional and marginal scores, and why the marginal score is a posterior mean.

MIT 6.S184: Flow Matching and Diffusion Models - Lecture 01 - Flow and Diffusion Models (2026)Peter Holderrieth · 78 min

The ODE and SDE view of generative models that this lesson builds on.

• Song, Sohl-Dickstein, Kingma, Kumar, Ermon and Poole, Score-Based Generative Modeling through Stochastic Differential Equations (ICLR 2021): §3.2 (reverse SDE), §3.3 (the training loss), §3.4 (VE and VP SDEs), §4.3 and Appendix D (probability-flow ODE).
• Yang Song, Generative Modeling by Estimating Gradients of the Data Distribution (blog post): the sections on perturbing data with an SDE, reversing the SDE and the probability-flow ODE.
• Holderrieth and Erives, An Introduction to Flow Matching and Diffusion Models (MIT 6.S184 notes): conditional and marginal scores, the SDE extension trick, and score matching = denoising score matching + a constant.
• Vincent, A Connection Between Score Matching and Denoising Autoencoders (2011): the original proof of Theorem 1. Hyvärinen, Estimation of Non-Normalized Statistical Models by Score Matching (JMLR 2005): explicit score matching by integration by parts.
• Karras, Aittala, Aila and Laine, Elucidating the Design Space of Diffusion-Based Generative Models (EDM, 2022): the $\sigma = t$ convention, Heun sampling and the time grid.
• Applications: Price et al., Probabilistic weather forecasting with machine learning (GenCast, Nature 2025); Rozet and Louppe, Score-based Data Assimilation (2023); Chung et al., Diffusion Posterior Sampling for General Noisy Inverse Problems (2023); Palgunadi et al., High Resolution Seismic Waveform Generation using Denoising Diffusion (2024); Jung et al., Broadband Ground Motion Synthesis by Diffusion Model with Minimal Condition (ICML 2025).