Reverse SDEs, the probability-flow ODE and samplers

Before this lesson: lesson 06 (the score, Langevin dynamics, denoising score matching), lesson 07 (the forward process, ancestral sampling, $\beta_t$ vs $\tilde\beta_t$) and lesson 08 (conversions between $\varepsilon$, $x_0$ and the score; Tweedie's formula; the EDM form $\sigma = \sigma_t/a_t$). For the order of a numerical method, see time stepping in the ML for science track. The lesson gives the minimum about SDEs that you need; the ML for science track treats Itô calculus and the Fokker–Planck equation in depth. Time: about 115 minutes, plus the notebook.

---

Why this matters #

Training fixes the network. The sampler decides how often you call it, and part of the sample quality. The original DDPM used 1,000 network evaluations per sample. With the same trained network and no retraining, DDIM, Heun or DPM-Solver reach similar quality with 10–50 evaluations. Every image, video and audio diffusion library (k-diffusion, diffusers) offers these samplers by name.

In science the cost multiplies. GenCast, a diffusion model for 15-day global weather ensembles, calls its denoiser 39 times per 12-hour forecast step (a second-order solver on 20 noise levels), for 30 steps and many ensemble members. For a diffusion or flow emulator of physical fields, the sampler is the main control of inference cost. "How do you make sampling 20 times faster without retraining?" is a standard interview question. This lesson is the answer, and it is also the basis of flow matching (lesson 10), guidance (lesson 11) and few-step generation (lesson 12).

1. Predict first #

2. The picture #

Sampling solves a differential equation backward in noise level #

At a noisy point $x$ with noise level $\sigma$, the trained network gives the denoiser $D(x; \sigma) \approx \mathbb E[x_0 \mid x]$ (lesson 08), and so the score $\nabla_x \log p_\sigma(x) = (D(x;\sigma) - x)/\sigma^2$ by Tweedie's formula. A sampler uses these values to move a point from high noise to low noise. It makes three choices:

• the solver: how to move from one noise level to the next with the help of $D$ (Euler, Heun, ...);
• the noise: deterministic steps (an ODE) or fresh random noise in every step (an SDE);
• the step schedule: at which noise levels to call the network.

The cost of a sampler is its number of network function evaluations (NFE). The network is by far the most expensive part, so we always compare samplers at equal NFE, not at equal numbers of steps.

Same distributions, different paths #

The figure below shows the central fact of this lesson in one dimension. The data are two narrow Gaussians at $-1$ and $+1$. The blue background is the density $p_t$ of the noised data at each time $t$ (time runs from data at $t = 0$ to noise at $t = 1$). Both samplers start from the same noise points at $t = 1$ and run backward to $t = 0$. The score is computed exactly, so no network is involved.

Interactive figure: loading…

So the two samplers agree on everything a sampler must get right (the distribution at each time) and disagree on the paths. Section 3 proves this with the Fokker–Planck equation.

Three samplers in one sentence each #

• Ancestral sampling (DDPM): take one small step down with the posterior mean, then add posterior noise. This is a discretization of the reverse SDE.
• DDIM: predict $\hat x_0$, then put it back at the next, lower noise level, mixing the predicted noise and some fresh noise. With no fresh noise the step is deterministic.
• ODE solvers (Euler, Heun, DPM-Solver): integrate the probability-flow ODE with a standard numerical method. DDIM without fresh noise is one of them: Euler's method.

The next figure runs these samplers on the 2D "8 Gaussians" data with the exact score. All differences you see come from the sampler, because the score has no error.

Interactive figure: loading…

3. Definitions and main results #

The forward process as an SDE #

Lesson 07 adds noise in $T$ discrete steps: $x_t = \sqrt{1 - \beta_t}\,x_{t-1} + \sqrt{\beta_t}\,z_t$. Let the number of steps grow and each step shrink. Write $\beta_t = \beta(t)\,\Delta t$ for a smooth function $\beta(t)$ and a small time step $\Delta t$. Then $\sqrt{1 - \beta(t)\Delta t} \approx 1 - \tfrac12\beta(t)\Delta t$, so

A Gaussian increment with variance $\Delta t$ is an increment of Brownian motion $w$: $dw = \sqrt{dt}\,z$. In the limit we get a stochastic differential equation (SDE):

The drift $f$ is the deterministic part of the motion; the diffusion coefficient $g$ scales the random part. This is the variance-preserving (VP) SDE (Song et al. 2021). Its marginals are the ones you know: $x_t = a_t x_0 + \sigma_t\varepsilon$ with $a_t = \exp\big(-\tfrac12\int_0^t \beta(s)\,ds\big)$ and $\sigma_t^2 = 1 - a_t^2$.

The variance-exploding (VE) SDE has no drift, only noise: $x_t = x_0 + \sigma(t)\varepsilon$, which is the SDE

(Variances of independent increments add, so the noise added between $t$ and $t + dt$ must have variance $d\sigma^2$.) EDM uses the VE SDE with $\sigma(t) = t$, so $g(t) = \sqrt{2t}$: the time is the noise level.

Euler–Maruyama is the simplest way to simulate an SDE. On a grid with step $h$:

It is Euler's method for the drift plus a Gaussian kick with the variance $g^2 h$ that the noise adds in time $h$. Exercise 1 of the notebook checks that it reproduces the closed-form marginals of the VP SDE.

Theorem 1 (reverse-time SDE, Anderson 1982) #

Intuition first. Run the forward SDE backward. The noise must now remove the spread it added, which no random kick can do on average. The fix is a drift toward high density: the score $\nabla\log p_t$ points from a point toward where the density of $p_t$ is larger. A drift of size $g^2\nabla\log p_t$ exactly undoes the spreading that noise of strength $g$ causes, while a fresh noise term keeps the spread right.

Theorem. Let $x$ solve the forward SDE $dx = f(x, t)\,dt + g(t)\,dw$ on $[0, T]$ with $x(0) \sim p_0$ and marginals $p_t$ (smooth and positive). Then the process that runs from $t = T$ to $t = 0$,

has marginal $p_t$ at every time. Here $dt < 0$ (time runs backward) and $\bar w$ is a Brownian motion in reversed time.

In practice, one Euler–Maruyama step from $t$ down to $t - h$ (with $h > 0$) is

The score enters with a plus sign: $+g^2 s\,h$ pushes the point toward high density. For the VP SDE the step is $x_t + \big[\tfrac12\beta(t)\,x_t + \beta(t)\,s(x_t, t)\big]h + \sqrt{\beta(t)h}\;z$.

Proof (with the Fokker–Planck equation)

We prove that the marginals agree, which is all a sampler needs. (Anderson proved more: the whole reversed path has the law of the forward path run backward.)

1. The Fokker–Planck equation. The density $p_t$ of the solution of $dx = f\,dt + g\,dw$ satisfies

$$\partial_t p = -\nabla\cdot(f\,p) + \tfrac12 g^2 \Delta p .$$

The first term moves mass along the drift; the second spreads it like heat. (A derivation is in the ML for science track; here we use it as a known fact.) 2. Reverse the clock. Let $\tau = T - t$ and $q_\tau = p_{T - \tau}$. Then $\partial_\tau q = -\partial_t p = \nabla\cdot(f\,q) - \tfrac12 g^2 \Delta q$. 3. A candidate. An SDE $dy = b\,d\tau + g\,dw_\tau$ in the clock $\tau$ has the Fokker–Planck equation $\partial_\tau q = -\nabla\cdot(b\,q) + \tfrac12 g^2\Delta q$. 4. Match the two. We need $-\nabla\cdot(b\,q) + \tfrac12 g^2 \Delta q = \nabla\cdot(f\,q) - \tfrac12 g^2\Delta q$, that is $-\nabla\cdot(b\,q) = \nabla\cdot(f\,q) - g^2\Delta q$. 5. The key identity. $\Delta q = \nabla\cdot(\nabla q) = \nabla\cdot(q\,\nabla\log q)$, because $\nabla \log q = \nabla q / q$. So the choice $b = -f + g^2\nabla\log q$ gives $-\nabla\cdot(b\,q) = \nabla\cdot(f q) - g^2\nabla\cdot(q\nabla\log q) = \nabla\cdot(fq) - g^2\Delta q$. ✓ 6. Uniqueness. Both processes start from $p_T$ and their densities solve the same equation, so (under the usual uniqueness of solutions) they are equal at every $\tau$. 7. Written in the original clock ($d\tau = -dt$), the drift $b = -f + g^2\nabla\log p_t$ becomes $f - g^2\nabla\log p_t$ multiplying a negative $dt$: the formula of the theorem. $\blacksquare$

Theorem 2 (the probability-flow ODE, Song et al. 2021) #

Theorem. The ordinary differential equation

has the same marginals $p_t$ as the forward SDE: if $x(T) \sim p_T$ and you solve it backward, then $x(t) \sim p_t$ for every $t$ (and the same holds forward from $x(0) \sim p_0$).

Proof.

1. Write the diffusion term of the Fokker–Planck equation as a flux, with the identity of step 5 above: $\tfrac12 g^2\Delta p = \nabla\cdot\big(\tfrac12 g^2\,p\,\nabla\log p\big)$.
2. Then the Fokker–Planck equation reads $\partial_t p = -\nabla\cdot\Big(\big(f - \tfrac12 g^2\nabla\log p\big)\,p\Big)$.
3. This is a continuity equation $\partial_t p + \nabla\cdot(p\,u) = 0$ with the velocity $u = f - \tfrac12 g^2\nabla\log p$, and no diffusion term. A continuity equation is exactly the equation for the density of points that move by the ODE $dx/dt = u(x, t)$ (lesson 04, continuous flows).
4. Same equation and same initial density, so the same marginals. $\blacksquare$

The ODE moves mass like a fluid: it never adds randomness, it only pushes points along a velocity field. The factor $\tfrac12$ is not optional: without it the ODE removes the spread twice as fast as the noise added it. (The notebook shows this for Gaussian data: from $\sigma = 80$, the right ODE ends with standard deviation 1, the ODE without the $\tfrac12$ with standard deviation 0.0125.)

A whole family of samplers. The same flux trick shows that for every $\gamma(t) \ge 0$ the backward process

has the marginals $p_t$. With $\gamma = g$ you get the reverse SDE; with $\gamma = 0$ the probability-flow ODE. (Check: its Fokker–Planck equation in reversed time has the drift term $\nabla\cdot(fq) - \tfrac12(g^2 + \gamma^2)\Delta q$ and the diffusion term $+\tfrac12\gamma^2\Delta q$, which add up to $\nabla\cdot(fq) - \tfrac12 g^2\Delta q$, as in step 2 of the proof of Theorem 1.) The amount of noise is a free design choice. One way to see the extra noise: it is Langevin dynamics (lesson 06) at a fixed noise level, added on top of the ODE. Langevin dynamics keeps $p_t$ unchanged, so it can be added in any amount.

The probability-flow ODE in EDM coordinates. Take the VE SDE with $\sigma(t) = t$: $f = 0$ and $g^2 = 2t$. The ODE is $dx/dt = -t\,\nabla\log p_t(x)$. Tweedie's formula for $x = x_0 + \sigma\varepsilon$ is $D(x; \sigma) = x + \sigma^2\nabla\log p_\sigma(x)$, so $\nabla\log p_\sigma = (D - x)/\sigma^2$ and

The slope points away from the denoiser output. Read backward (decreasing $\sigma$), the point moves toward $D$. This form is the one used by EDM and by most modern samplers.

Ancestral sampling is a discretization of the reverse SDE #

The DDPM step of lesson 07 is $x_{t-1} = \frac{1}{\sqrt{1 - \beta_t}}\big(x_t - \frac{\beta_t}{\sqrt{1 - \bar\alpha_t}}\hat\varepsilon\big) + \sqrt{\beta_t}\,z$. The score is $s = -\hat\varepsilon/\sqrt{1 - \bar\alpha_t}$, so the step is $x_{t-1} = (x_t + \beta_t\,s)/\sqrt{1 - \beta_t} + \sqrt{\beta_t}\,z$. With $1/\sqrt{1 - \beta_t} = 1 + \tfrac12\beta_t + O(\beta_t^2)$:

With $\beta_t = \beta(t)\Delta t$ this is the Euler–Maruyama step of the reverse VP SDE above. So DDPM sampling is a first-order discretization of Theorem 1. It is accurate only when every $\beta_t$ is small, which is why DDPM needs many steps.

DDIM: a family of samplers for one trained model #

Song, Meng and Ermon (2021) asked which samplers can use a network trained with the DDPM loss. The loss $\mathbb E\|\varepsilon - \varepsilon_\theta(x_t, t)\|^2$ only sees pairs $(x_0, x_t)$ at single times, through $q(x_t \mid x_0) = \mathcal N(a_t x_0, \sigma_t^2 I)$. Any process with these same marginals has the same optimal network $\mathbb E[\varepsilon \mid x_t]$. They found a family of such processes with one free noise scale $\nu_t$ per step. Sampling from it gives the DDIM update (Song et al. 2021, eq. 12). From level $t$ to any lower level $t'$ (written $t - 1$):

with $z \sim \mathcal N(0, I)$ and, with a parameter $\eta \in [0, 1]$ (eq. 16 of the paper),

(Song et al. call $\nu_t$ "$\sigma_t$"; we keep $\sigma_t$ for the noise scale.)

• $\eta = 0$: no fresh noise. The step is deterministic, and so is the whole sample once $x_T$ is fixed.
• $\eta = 1$: $\nu_t^2$ equals the DDPM posterior variance $\tilde\beta$ of this step, and the mean equals the DDPM posterior mean. So $\eta = 1$ is ancestral DDPM sampling with variance $\tilde\beta$. (Worked example (a) checks this with numbers.)
• Any subsequence of levels. Nothing in the derivation needs consecutive steps, so the same formula jumps between any two levels, for example 50 levels out of 1,000.

Why every value of η keeps the DDPM marginals (the non-Markovian forward process)

Song et al. define a forward process that runs down in $t$ and may look at $x_0$:

$$q_\nu(x_{t-1} \mid x_t, x_0) = \mathcal N\Big(a_{t-1}x_0 + \sqrt{\sigma_{t-1}^2 - \nu_t^2}\;\frac{x_t - a_t x_0}{\sigma_t},\ \nu_t^2 I\Big),$$

started from $q(x_T \mid x_0) = \mathcal N(a_T x_0, \sigma_T^2 I)$.

Claim. $q_\nu(x_t \mid x_0) = \mathcal N(a_t x_0, \sigma_t^2 I)$ for every $t$ and every choice of $\nu$.

1. Assume it holds at level $t$. Write $x_t = a_t x_0 + \sigma_t\varepsilon$; then $(x_t - a_t x_0)/\sigma_t = \varepsilon$.
2. A draw from the kernel is $x_{t-1} = a_{t-1}x_0 + \sqrt{\sigma_{t-1}^2 - \nu_t^2}\;\varepsilon + \nu_t z$, with $z$ independent of $\varepsilon$.
3. Given $x_0$, this is a sum of independent Gaussians: mean $a_{t-1}x_0$, variance $(\sigma_{t-1}^2 - \nu_t^2) + \nu_t^2 = \sigma_{t-1}^2$.
4. The claim holds at $t = T$ by definition. Induction downward gives every $t$. $\blacksquare$

Consequences. The process is non-Markovian ($x_{t-1}$ depends on $x_0$, not only on $x_t$), but the training loss never sees this: it uses only $q(x_t \mid x_0)$. Song et al. (Theorem 1) also show that the ELBO of every member is a weighted $L_{\text{simple}}$ plus a constant. Sampling replaces the unknown $x_0$ by $\hat x_0$ and $\varepsilon$ by $\hat\varepsilon$, which gives the update above.

$\eta = 1$ is DDPM. Let $\beta' = 1 - a_t^2/a_{t-1}^2$ be the $\beta$ of the step. Then $\nu_t^2 = \frac{\sigma_{t-1}^2}{\sigma_t^2}\beta' = \tilde\beta'$, the posterior variance of lesson 07. The coefficient of $x_t$ in the mean becomes $\sqrt{\sigma_{t-1}^2 - \tilde\beta'}/\sigma_t$, and its square is $\frac{\sigma_{t-1}^2}{\sigma_t^2}\big(1 - \frac{\beta'}{\sigma_t^2}\big) = \frac{\sigma_{t-1}^2}{\sigma_t^4}(\sigma_t^2 - \beta') = \frac{\sigma_{t-1}^4 (1 - \beta')}{\sigma_t^4}$, using $\sigma_t^2 - \beta' = 1 - a_t^2 - 1 + a_t^2/a_{t-1}^2 = (1 - \beta')\,\sigma_{t-1}^2$. So the coefficient is $\sqrt{1 - \beta'}\,\sigma_{t-1}^2/\sigma_t^2$: exactly the DDPM coefficient of $x_t$ in $\tilde\mu$.

Theorem 3 (DDIM with $\eta = 0$ is Euler's method for the probability-flow ODE) #

In the EDM coordinates $\tilde x = x/a$ and $\tilde\sigma = \sigma/a$, the $\eta = 0$ DDIM step is

one Euler step of the probability-flow ODE $\dfrac{d\tilde x}{d\tilde\sigma} = \dfrac{\tilde x - D(\tilde x; \tilde\sigma)}{\tilde\sigma}$, with the network's $\hat x_0$ in place of $D$.

Proof.

1. Set $\nu_t = 0$ in the DDIM update and divide by $a_{t-1}$: $\tilde x_{t-1} = \hat x_0 + \frac{\sigma_{t-1}}{a_{t-1}}\hat\varepsilon = \hat x_0 + \tilde\sigma_{t-1}\hat\varepsilon$.
2. Rearrange the definition of $\hat x_0$: $x_t = a_t\hat x_0 + \sigma_t\hat\varepsilon$. Divide by $a_t$: $\tilde x_t = \hat x_0 + \tilde\sigma_t\hat\varepsilon$ and $\hat\varepsilon = (\tilde x_t - \hat x_0)/\tilde\sigma_t$.
3. Put this into step 1: $\tilde x_{t-1} = \hat x_0 + \tilde\sigma_{t-1}(\tilde x_t - \hat x_0)/\tilde\sigma_t$. Add and subtract $\tilde x_t$: $\tilde x_{t-1} = \tilde x_t + (\tilde\sigma_{t-1} - \tilde\sigma_t)(\tilde x_t - \hat x_0)/\tilde\sigma_t$.
4. In these coordinates $\tilde x = x_0 + \tilde\sigma\varepsilon$ is the VE form, whose probability-flow ODE is $d\tilde x/d\tilde\sigma = (\tilde x - D)/\tilde\sigma$ (the boxed equation). Step 3 is one Euler step of it from $\tilde\sigma_t$ to $\tilde\sigma_{t-1}$, with the slope evaluated at the start. $\blacksquare$

Interactive figure: loading…

What Theorem 3 tells you.

• DDIM is a first-order method: its global error is $O(h)$ in the step size $h$. Second-order methods on the same ODE need fewer NFE for the same error.
• In these coordinates the slope points at the denoiser output: one Euler step to $\tilde\sigma = 0$ returns $D(\tilde x; \tilde\sigma)$. A step is exact when $D$ does not change along it. If the data are a single point $x^\ast$, then $D \equiv x^\ast$, every path is a straight line, and one DDIM step is exact.
• The coordinates matter. Euler steps in $t$ on the VP form of the ODE give a different, worse sampler. In $(\tilde\sigma, \tilde x)$ the paths are nearly straight at high noise, which is why DDIM works with few steps. Lesson 10 builds on the same idea: flow matching chooses a path that is as straight as possible.

The Gaussian toy: exact paths #

Take 1D data $x_0 \sim \mathcal N(0, \sigma_{\text{data}}^2)$ in EDM coordinates. The denoiser is linear, $D(x; \sigma) = \sigma_{\text{data}}^2 x/(\sigma_{\text{data}}^2 + \sigma^2)$ (lesson 08), so the ODE is

(Separate the variables: $dx/x = \sigma\,d\sigma/(\sigma_{\text{data}}^2 + \sigma^2)$, so $\ln x = \tfrac12\ln(\sigma_{\text{data}}^2 + \sigma^2) + C$.) Each start value gives exactly one data value: the ODE is a deterministic, invertible map from noise to data. For $\sigma \gg \sigma_{\text{data}}$, $x(\sigma)$ is almost proportional to $\sigma$: a straight line that one large step follows almost exactly. Near $\sigma \approx \sigma_{\text{data}}$ the path bends. Easy at high noise, hard at low noise: this is the reason for the step schedules below.

ODE or SDE? #

• ODE ($\eta = 0$; Euler, Heun, DPM-Solver). Deterministic and invertible: DDIM inversion maps an image to its noise (used for editing and interpolation). It gives exact likelihoods with the instantaneous change of variables of lesson 04. It is usually the best choice at low NFE. But errors (network error plus discretization error) are carried along, never corrected.
• SDE (reverse SDE, ancestral, DDIM with $\eta > 0$). The fresh noise and the extra score drift act like Langevin steps at a fixed noise level. They pull samples back toward $p_t$ and so partly correct earlier errors. The price: every step must remove the noise it adds, so it needs more steps.
• Numbers (Song et al. 2021, Table 1, CIFAR-10 FID, lower is better): with 10 steps, $\eta = 0$ gives 13.36 and $\eta = 1$ gives 41.07; with 100 steps, 4.16 and 5.78; with 1,000 steps, 4.04 and 4.73. Ancestral sampling with the larger variance $\beta_t$ is the best at 1,000 steps (3.17) and useless at 10 steps (367).
• EDM's stochastic sampler (Karras et al. 2022, Algorithm 2) adds a controlled amount of noise ("churn") in a band of noise levels, then takes a deterministic Heun step. Too much churn loses detail.

In the 8-Gaussians figure the score is exact, so there is no network error for the noise to correct. There, at 10 NFE on the Karras schedule, Euler puts about 84% of the samples on a mode and the reverse SDE about 62%; at 50 NFE both are above 95%.

Solvers: accuracy per NFE #

From here on we work in EDM coordinates and write $x$ and $\sigma$ (no tildes). The ODE is $dx/d\sigma = d(x, \sigma)$ with the slope $d(x, \sigma) = (x - D(x;\sigma))/\sigma$, solved on a decreasing grid $\sigma_0 > \sigma_1 > \dots > \sigma_N = 0$ with step $h_i = \sigma_{i+1} - \sigma_i < 0$.

• Euler (= DDIM): $x_{i+1} = x_i + h_i\,d(x_i, \sigma_i)$. One NFE per step. Global error $O(h)$: halving all steps halves the error.
• Heun (EDM, Algorithm 1): an Euler predictor $x' = x_i + h_i d_i$, then the trapezoidal corrector

Two NFE per step. Global error $O(h^2)$: halving the steps divides the error by about 4.

• The last step. The slope $(x - D)/\sigma$ is undefined at $\sigma = 0$, so the last step into $\sigma = 0$ is a plain Euler step (it returns $D(x; \sigma_{N-1})$). So $N$ Heun steps that end at 0 cost $2N - 1$ NFE: 35 for EDM's 18 levels, which reach FID 1.79 on class-conditional CIFAR-10.

The order of a method is the power $p$ in "global error $\approx C h^p$". On log–log axes, error against NFE is a line with slope $-p$.

Interactive figure: loading…

DPM-Solver: solve the linear part exactly

Lu et al. (2022) change the variable to $u = \log\sigma$. Then $dx/du = \sigma\,dx/d\sigma = x - D(x; e^u)$. The term $x$ is linear and can be integrated exactly (variation of constants):

$$x(u_{i+1}) = \frac{\sigma_{i+1}}{\sigma_i}\,x(u_i) - \int_{u_i}^{u_{i+1}} e^{\,u_{i+1} - u}\,D\big(x(u); e^{u}\big)\,du .$$

Only $D$ is approximated, by a polynomial in $u$. A constant $D = \hat x_0$ gives

$$x_{i+1} = \frac{\sigma_{i+1}}{\sigma_i}\,x_i + \Big(1 - \frac{\sigma_{i+1}}{\sigma_i}\Big)\hat x_0 ,$$

which is the DDIM step again (DPM-Solver-1 = DDIM, Lu et al., section 4.1). Check: the Euler step $x_i + (\sigma_{i+1} - \sigma_i)(x_i - \hat x_0)/\sigma_i$ is the same expression. DPM-Solver-2 and -3 use 2 and 3 NFE per step for orders 2 and 3. DPM-Solver++(2M) makes $D$ linear in $u$ from the current and the previous denoiser outputs: order 2 for 1 NFE per step. Lu et al. write the time as $\lambda = -\log\sigma$ (half the log-SNR) and report good samples with 10–20 NFE.

Step schedules #

Karras et al. (2022, eq. 5) choose $N$ levels between $\sigma_{\max}$ and $\sigma_{\min}$ with a parameter $\rho$:

• $\rho = 1$: uniform in $\sigma$.
• $\rho \to \infty$: uniform in $\log\sigma$ (a geometric grid). Proof: with $u = i/(N-1)$ and $c^{1/\rho} = 1 + (\ln c)/\rho + O(\rho^{-2})$, the bracket is $1 + [(1-u)\ln\sigma_{\max} + u\ln\sigma_{\min}]/\rho + O(\rho^{-2})$, and its $\rho$-th power tends to $\sigma_{\max}^{1-u}\sigma_{\min}^{u}$.
• EDM default: $\sigma_{\min} = 0.002$, $\sigma_{\max} = 80$, $\rho = 7$, for data scaled to $\sigma_{\text{data}} = 0.5$. EDM found that $\rho = 3$ nearly equalizes the error per step, but $\rho$ between 5 and 10 gives better images: errors near $\sigma_{\min}$ have a large effect on what you see.

Interactive figure: loading…

Why the steps crowd at low noise.

1. The local error of a step of size $h$ is about $\tfrac12 h^2 |x''|$ for Euler (and $\propto h^3$ for Heun): it is large where the path bends.
2. In the Gaussian toy, $x(\sigma) = c\sqrt{\sigma_{\text{data}}^2 + \sigma^2}$ gives $x''(\sigma) = c\,\sigma_{\text{data}}^2/(\sigma_{\text{data}}^2 + \sigma^2)^{3/2}$. It is largest at $\sigma = 0$ and decays like $\sigma^{-3}$. For $\sigma_{\text{data}} = 0.5$ the curvature at $\sigma = 0$ is $(100.25/0.25)^{3/2} \approx 8{,}000$ times the curvature at $\sigma = 10$. For the same Euler error, a step at $\sigma = 10$ may be $\sqrt{8000} \approx 90$ times longer.
3. Real data have structure at many scales. The path bends whenever $\sigma$ passes a length scale of the data: the distance between modes, the width of a mode, fine detail. At high noise the denoiser output stays near the data mean, so the path is nearly straight.

An adaptive ODE solver finds this spacing by itself (section 5): with no schedule given, it chooses an almost geometric grid.

Where the schedule starts. The first level must be the highest level the model was trained at ($\sigma_{\max}$, or $t = T$). With $T = 1000$ and 10 steps, "leading" spacing uses $t = 901, 801, \dots, 1$ and never visits $t = T$; "trailing" spacing uses $t = 1000, 900, \dots, 100$ (Lin et al. 2024). Use trailing spacing; see lesson 08 for the related zero-terminal-SNR problem.

Answers to section 1. (1) No retraining. Every $\eta$ keeps the DDPM marginals $q(x_t \mid x_0)$, and the loss sees only these, so the same network is optimal. With $\eta = 0$ the output is a deterministic function of $x_T$: an ODE solve (Theorem 3). (2) Not the same point: the SDE path is random. But the same distribution: both have the marginals $p_t$ at every time (Theorems 1 and 2). (3) Not uniformly. Most levels go to low noise, where the paths bend. With $\rho = 7$, 7 of the 18 levels lie below $\sigma = 1$; a uniform grid puts only its last level there.

4. Worked example #

(a) One DDIM step by hand. 1D, $\bar\alpha_t = 0.5$, $\bar\alpha_{t-1} = 0.8$, $x_t = 1.0$, network output $\hat\varepsilon = 0.6$. So $a_t = \sigma_t = \sqrt{0.5} = 0.70711$, $a_{t-1} = \sqrt{0.8} = 0.89443$, $\sigma_{t-1} = \sqrt{0.2} = 0.44721$.

1. $\hat x_0 = (x_t - \sigma_t\hat\varepsilon)/a_t = (1 - 0.70711 \cdot 0.6)/0.70711 = 0.57574/0.70711 = 0.81421$.
2. $\eta = 0$: $x_{t-1} = a_{t-1}\hat x_0 + \sigma_{t-1}\hat\varepsilon = 0.89443 \cdot 0.81421 + 0.44721 \cdot 0.6 = 0.72825 + 0.26833 = 0.99658$.
3. Check as Euler. $\tilde\sigma_t = 0.70711/0.70711 = 1$, $\tilde\sigma_{t-1} = 0.44721/0.89443 = 0.5$, $\tilde x_t = 1/0.70711 = 1.41421$. The slope is $(\tilde x_t - \hat x_0)/\tilde\sigma_t = (1.41421 - 0.81421)/1 = 0.6 = \hat\varepsilon$. Euler: $\tilde x_{t-1} = 1.41421 + (0.5 - 1) \cdot 0.6 = 1.11421$. Back to VP: $0.89443 \cdot 1.11421 = 0.99658$ ✓.
4. $\eta = 1$. The step has $\beta' = 1 - 0.5/0.8 = 0.375$. $\nu = \frac{0.44721}{0.70711}\sqrt{0.375} = 0.63246 \cdot 0.61237 = 0.38730$, so $\nu^2 = 0.15$. The DDPM posterior variance is $\tilde\beta' = \frac{1 - 0.8}{1 - 0.5}\cdot 0.375 = 0.4 \cdot 0.375 = 0.15$ ✓.
5. The $\eta = 1$ mean: $0.72825 + \sqrt{0.2 - 0.15}\cdot 0.6 = 0.72825 + 0.22361 \cdot 0.6 = 0.86242$. The DDPM posterior mean (lesson 07) is $\frac{a_{t-1}\beta'}{1 - \bar\alpha_t}\hat x_0 + \frac{\sqrt{1 - \beta'}\,(1 - \bar\alpha_{t-1})}{1 - \bar\alpha_t}x_t = 0.67082 \cdot 0.81421 + 0.31623 \cdot 1 = 0.86242$ ✓.

(b) Euler vs Heun for one step, with the exact answer. Data $\mathcal N(0, 1)$ ($\sigma_{\text{data}} = 1$). Go from $\sigma = 2$, $x = 1$ to $\sigma = 1$, so $h = -1$. The slope is $d(x, \sigma) = \sigma x/(1 + \sigma^2)$.

• Exact: $x(1) = 1 \cdot \sqrt{(1 + 1)/(1 + 4)} = \sqrt{0.4} = 0.63246$.
• Euler (1 NFE): $d_0 = 2 \cdot 1/5 = 0.4$, so $x_E = 1 + (-1)(0.4) = 0.6$. Error $-0.0325$.
• Heun (2 NFE): $d_0' = d(0.6, 1) = 0.6/2 = 0.3$, so $x_H = 1 - (0.4 + 0.3)/2 = 0.65$. Error $+0.0175$.

(c) Karras grids for $N = 5$ ($\sigma_{\max} = 80$, $\sigma_{\min} = 0.002$). For $\rho = 7$: $80^{1/7} = 1.8701$ and $0.002^{1/7} = 0.4116$; the middle level is $\big((1.8701 + 0.4116)/2\big)^7 = 1.14085^7 = 2.515$. Then $\sigma_5 = 0$:

$\rho$	$\sigma_0$	$\sigma_1$	$\sigma_2$	$\sigma_3$	$\sigma_4$
1 (uniform)	80	60.0005	40.001	20.0015	0.002
3	80	34.75	10.90	1.609	0.002
7 (EDM)	80	17.53	2.515	0.1698	0.002
$\infty$ (geometric)	80	5.657	0.4	0.02828	0.002

The geometric middle level is $\sqrt{80 \cdot 0.002} = 0.4$, the geometric mean of the end points.

(d) Schedule and solver together. Data $\tfrac12\mathcal N(-1, 0.1^2) + \tfrac12\mathcal N(1, 0.1^2)$ with the exact denoiser; 500 paths from $\sigma = 80$ to $\sigma = 0.002$ on 18 levels (17 steps, no step to 0). The table gives the largest error against an accurate reference solve:

grid	Euler (17 NFE)	Heun (34 NFE)
$\rho = 1$	0.86	8.5
$\rho = 3$	0.46	0.14
$\rho = 7$	0.26	0.048
geometric	0.25	0.039

• The schedule comes first. On the uniform grid the last step, from $\sigma = 4.7$ to $0.002$, is far too long, and Heun is even worse than Euler. At $\sigma = 4.7$ the denoiser output lies between the modes, so the Euler predictor stops short of the nearest mode; the slope at the predicted point then points the other way, and the average of the two slopes ends far from the data.
• On a good grid Heun wins at equal cost. With $\rho = 7$, Heun on 10 levels (18 NFE) has error 0.18, against 0.26 for Euler on 18 levels (17 NFE); Heun on 18 levels (34 NFE) has 0.048, against 0.11 for Euler on 40 levels (39 NFE).

5. In code #

NumPy checks the worked examples (a)–(c). The notebook has these cells and the exercises.

import numpy as np
np.set_printoptions(suppress=True)

# (a) One DDIM step (VP), and the same step as Euler in the EDM coordinates
ab_t, ab_p, xt, eps_hat = 0.5, 0.8, 1.0, 0.6          # alpha_bar at t and t-1, x_t, network output
a_t, s_t, a_p, s_p = np.sqrt(ab_t), np.sqrt(1 - ab_t), np.sqrt(ab_p), np.sqrt(1 - ab_p)
x0_hat = (xt - s_t * eps_hat) / a_t
x_prev = a_p * x0_hat + s_p * eps_hat                  # eta = 0
st_t, st_p = s_t / a_t, s_p / a_p                      # sigma tilde = sigma_t / a_t
euler = xt / a_t + (st_p - st_t) * (xt / a_t - x0_hat) / st_t
print(round(x0_hat, 5), round(x_prev, 5), round(euler * a_p, 5))   # 0.81421 0.99658 0.99658

# eta = 1: the DDIM noise and mean equal the DDPM posterior of this (respaced) step
beta = 1 - ab_t / ab_p
nu = 1.0 * (s_p / s_t) * np.sqrt(1 - ab_t / ab_p)
mean_ddim = a_p * x0_hat + np.sqrt(s_p**2 - nu**2) * eps_hat
mean_ddpm = a_p * beta / (1 - ab_t) * x0_hat + np.sqrt(1 - beta) * (1 - ab_p) / (1 - ab_t) * xt
print(round(nu**2, 5), round((1 - ab_p) / (1 - ab_t) * beta, 5))   # 0.15 0.15
print(round(mean_ddim, 5), round(mean_ddpm, 5))                    # 0.86242 0.86242

# (b) Euler and Heun for one step of dx/dsigma = sigma x / (1 + sigma^2) (data N(0, 1))
d = lambda x, s: s * x / (1 + s**2)
x_e = 1 + (1 - 2) * d(1, 2)
x_h = 1 + (1 - 2) * (d(1, 2) + d(x_e, 1)) / 2
print(x_e, x_h, round(np.sqrt(2 / 5), 5))                          # 0.6 0.65 0.63246

# (c) Karras grids
def karras(n, rho, smax=80.0, smin=0.002):
    u = np.arange(n) / (n - 1)
    if np.isinf(rho):
        return smax * (smin / smax) ** u
    return (smax ** (1 / rho) + u * (smin ** (1 / rho) - smax ** (1 / rho))) ** rho

for rho in (1, 3, 7, np.inf):
    print(rho, np.round(karras(5, rho), 4))
# 1 [80.     60.0005 40.001  20.0015  0.002 ]
# 3 [80.     34.7465 10.9031  1.6087  0.002 ]
# 7 [80.     17.5278  2.5152  0.1698  0.002 ]
# inf [80.      5.6569  0.4     0.0283  0.002 ]
print((karras(18, 7) < 1).sum(), (karras(18, 1) < 1).sum())        # 7 1

Second, give SciPy's adaptive Runge–Kutta solver (RK45) the probability-flow ODE of the two-mode mixture of example (d), with no schedule at all, and look at where it puts its steps. The denoiser is exact (Tweedie for a Gaussian mixture, with responsibilities in log space):

import numpy as np
from scipy.integrate import solve_ivp

mu, sd = np.array([-1.0, 1.0]), 0.1         # data: 0.5 N(-1, 0.1^2) + 0.5 N(1, 0.1^2)

def denoiser(x, sigma):
    """D(x; sigma) = E[x0 | x0 + sigma * eps = x] for the mixture."""
    v = sd**2 + sigma**2
    logw = -(x[:, None] - mu) ** 2 / (2 * v)
    w = np.exp(logw - logw.max(1, keepdims=True))
    w /= w.sum(1, keepdims=True)
    return (w * (mu + sd**2 * (x[:, None] - mu) / v)).sum(1)

def pf_ode(sigma, x):                       # dx/dsigma = (x - D(x; sigma)) / sigma
    return (x - denoiser(x, sigma)) / sigma

x_start = 80.0 * np.random.default_rng(0).standard_normal(500)
sol = solve_ivp(pf_ode, (80.0, 0.002), x_start, method="RK45", rtol=1e-5, atol=1e-5)
levels = sol.t[:-1]                         # sigma at the start of each accepted step
print("steps:", len(levels), " NFE:", sol.nfev)                       # steps: 25  NFE: 224
print(np.sum(levels >= 1), np.sum((levels >= 0.1) & (levels < 1)), np.sum(levels < 0.1))   # 18 5 2
print(round(np.median(sol.t[1:] / sol.t[:-1]), 2))                    # 0.75
print(np.mean(np.abs(np.abs(sol.y[:, -1]) - 1) < 0.4))               # 1.0: every path ends near a mode

The solver takes 25 steps. Their median ratio $\sigma_{\text{next}}/\sigma$ is 0.75, so it chooses an almost geometric grid by itself: 7 of its 25 steps lie below $\sigma = 1$, where a uniform grid would put none. The last two steps are long again: below $\sigma \approx 0.1$, the width of each mode, the denoiser barely changes $x$ and the path is straight. But 224 NFE is a lot (RK45 uses 6 evaluations per step, plus rejected steps). EDM tried RK45 and found that its cost outweighs its benefits; a fixed Heun grid is cheaper.

6. Where it shows up #

• Image and video models. Libraries offer these samplers by name: k-diffusion has sample_euler, sample_heun, sample_euler_ancestral, sample_dpmpp_2m and get_sigmas_karras; diffusers has DDIM, Euler, Heun and DPM-Solver++ schedulers. 20–50 steps is a common default.
• Weather. GenCast samples with a second-order DPM-Solver++ variant on a Karras grid ($\rho = 7$) with 20 noise levels (39 NFE) plus EDM-style churn. The ensemble cost is members × forecast steps × NFE.
• Flow matching and rectified flow (lesson 10). Sampling integrates $dx/dt = v_\theta(x, t)$ with Euler or Heun; the NFE accounting is the same. Straighter paths need fewer steps, for the reason of Theorem 3: a step is exact when the path is straight.
• Inversion, editing and likelihood. DDIM inversion runs the ODE from data to noise. Likelihoods (bits per dimension, lesson 14) need an accurate, usually adaptive ODE solve with a Hutchinson trace estimate (lesson 04).
• Guidance (lesson 11) changes the denoiser inside the same samplers; classifier-free guidance doubles the NFE of each step.
• Below 10–20 NFE, retrain: progressive distillation and consistency models (lesson 12).

7. Common mistakes #

• Comparing samplers at equal steps instead of equal NFE. A Heun step costs 2 NFE; classifier-free guidance doubles every count.
• Skipping steps in ancestral sampling but keeping the original $\beta_t$. On a subsequence, recompute $\beta' = 1 - \bar\alpha_t/\bar\alpha_{t_{\text{prev}}}$ for each step. DDIM's formula does this automatically.
• Forgetting the $\tfrac12$ in the probability-flow ODE. $dx/dt = f - g^2\nabla\log p_t$ without noise shrinks the spread twice as fast as it should. The reverse SDE has $g^2$; the ODE has $\tfrac12 g^2$.
• Heun's corrector on the last step into $\sigma = 0$. It divides by zero; use Euler there.
• Mixing coordinates. DDIM is Euler in $(\tilde\sigma, \tilde x)$. Euler on the VP ODE in $t$ is a different sampler with more curved paths.
• Starting below the top level ("leading" spacing) or with a nonzero terminal SNR (lesson 08).
• An SDE sampler at very low NFE. At 10 steps, $\eta = 1$ gives FID 41 where $\eta = 0$ gives 13 (CIFAR-10).
• A uniform $\sigma$ grid. Most steps are then spent where the path is straight (example (d)).
• Clipping $\hat x_0$ on unbounded data. Image samplers often clip $\hat x_0$ to $[-1, 1]$ (in diffusers, clip_sample=True is the default of the DDPM and DDIM schedulers). On physical fields this silently cuts the extremes; turn it off.

8. Check yourself #

1. One DDIM step with $\eta = 0$: $\bar\alpha_t = 0.36$, $\bar\alpha_{t-1} = 0.64$, $x_t = 0.8$, $\hat\varepsilon = 0.5$. Compute $x_{t-1}$ and check it as an Euler step.
   Answer

   $a_t = 0.6$, $\sigma_t = 0.8$, $a_{t-1} = 0.8$, $\sigma_{t-1} = 0.6$. $\hat x_0 = (0.8 - 0.8 \cdot 0.5)/0.6 = 0.6667$; $x_{t-1} = 0.8 \cdot 0.6667 + 0.6 \cdot 0.5 = 0.8333$. EDM coordinates: $\tilde\sigma_t = 4/3$, $\tilde\sigma_{t-1} = 3/4$, $\tilde x_t = 1.3333$, slope $(\tilde x_t - \hat x_0)/\tilde\sigma_t = 0.5$; $\tilde x_{t-1} = 1.3333 + (0.75 - 1.3333) \cdot 0.5 = 1.0417$, and $0.8 \cdot 1.0417 = 0.8333$ ✓.
2. Write the reverse-time SDE and the probability-flow ODE for the VP SDE with $\beta(t)$.
   Answer

   Reverse SDE: $dx = \big[-\tfrac12\beta(t)x - \beta(t)\nabla\log p_t(x)\big]dt + \sqrt{\beta(t)}\,d\bar w$ (with $dt < 0$). PF-ODE: $dx/dt = -\tfrac12\beta(t)\big[x + \nabla\log p_t(x)\big]$.
3. Why may the DDIM forward process be non-Markovian?
   Answer

   The training loss uses only the pairs $(x_0, x_t)$ at single times, through $q(x_t \mid x_0)$. Every member of the DDIM family has the DDPM marginals $q(x_t \mid x_0)$, so the optimal $\varepsilon_\theta(x_t, t) = \mathbb E[\varepsilon \mid x_t]$ is the same, whatever the joint law across different times.
4. For $N = 3$ levels between $\sigma_{\max} = 80$ and $\sigma_{\min} = 0.002$, what is the middle level for $\rho = 1$, $\rho = 7$ and $\rho \to \infty$?
   Answer

   $\rho = 1$: $40.001$. $\rho = 7$: $\big((1.8701 + 0.4116)/2\big)^7 = 2.515$. $\rho \to \infty$: $\sqrt{80 \cdot 0.002} = 0.4$, the geometric mean.
5. The data are a single point $x^\ast$. Show that DDIM with $\eta = 0$ reaches $x^\ast$ exactly in one step from any $\sigma$.
   Answer

   $D(x; \sigma) = x^\ast$ for every input, so the Euler step to $\sigma = 0$ gives $x + (0 - \sigma)(x - x^\ast)/\sigma = x^\ast$. The true path $x(\sigma) = x^\ast + \sigma\varepsilon$ is a straight line, so Euler has no discretization error.
6. In the reverse SDE, the score term is $g^2\nabla\log p_t$; in the probability-flow ODE it is $\tfrac12 g^2\nabla\log p_t$. Explain the factor in one or two sentences.
   Answer

   The forward Fokker–Planck equation has the spreading term $\tfrac12 g^2\Delta p$. The ODE has no noise, so its drift must undo exactly this term: $\tfrac12 g^2\nabla\log p$. The reverse SDE adds noise of strength $g$ again (another $\tfrac12 g^2\Delta p$), so its drift must undo twice as much: $g^2\nabla\log p$.

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. Write the reverse-time SDE for the forward SDE $dx = f\,dt + g\,dw$.

   Show answer

2. Write the probability-flow ODE for the forward SDE $dx = f\,dt + g\,dw$.

   Show answer

3. Give the VP SDE and say what discrete process it is the limit of.

   Show answer

4. Why do the reverse SDE and the PF-ODE have the same marginals?

   Show answer

5. Write the DDIM update with $\eta = 0$.

   Show answer

6. In which coordinates is DDIM ($\eta = 0$) exactly Euler's method, and for which ODE?

   Show answer

7. What are the order and the NFE count of Heun's method with $N$ steps ending at $\sigma = 0$?

   Show answer

8. Give the Karras step schedule.

   Show answer

9. ODE or SDE sampler: what does each do better?

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Make sampling 20 times faster without retraining." — DDIM, Heun or DPM-Solver++(2M) on a Karras grid: 10–50 NFE instead of 1,000. If retraining is allowed: distillation or consistency models.
• "What is the probability-flow ODE, and why does it give the right distribution?" — $dx/dt = f - \tfrac12 g^2\nabla\log p_t$. The Fokker–Planck equation of the forward SDE can be written as a continuity equation with this velocity.
• "Is DDIM deterministic? Why is that useful?" — With $\eta = 0$, yes: inversion, interpolation in noise space, reproducible samples, likelihoods.
• "What is DDIM in ODE terms?" — Euler on the PF-ODE in $(\sigma_t/a_t, x_t/a_t)$, so first order; it is also the first-order exponential integrator (DPM-Solver-1).
• "ODE or SDE sampler?" — ODE at low NFE. Noise corrects earlier errors when you can afford more steps; tune the amount empirically.
• "Heun with 18 noise levels: how many network calls?" — 35; the last step is Euler.
• "Why are the noise levels not uniform?" — The paths bend at low noise; Karras $\rho = 7$ puts short steps there.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. Simulate the VP SDE with Euler–Maruyama and compare the mean and variance of $x_t$ with the closed form.
2. Sample the 8 Gaussians with the reverse SDE and with the probability-flow ODE (exact score); compare the counts per mode and the spread.
3. Write the DDIM step with $\eta$; check that $\eta = 0$ equals an Euler step in $(\tilde\sigma, \tilde x)$ to rounding error, and that $\eta = 1$ has the DDPM posterior variance.
4. Euler vs Heun on the Gaussian toy: error against NFE; fit the slopes of the log–log lines.
5. Karras $\rho$ vs a uniform grid: sample quality at 10 NFE with the exact 8-Gaussians denoiser.

11. Watch and read #

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

Watch how flows and diffusion models are both defined by an ODE or SDE from noise to data, and how they are simulated with Euler and Euler–Maruyama.

DDIM Explained — Fast and Deterministic Diffusion ModelsDeepBean · 28 min

Follow the non-Markovian forward process and the role of the noise parameter; compare with the deeper box of section 3.

Diffusion Models From Scratch | Score-Based Generative Models Explained | Math ExplainedOutlier · 38 min

A step-by-step derivation of score-based models, the SDE view and the probability-flow ODE.

• Song, Sohl-Dickstein, Kingma, Kumar, Ermon, Poole, Score-Based Generative Modeling through Stochastic Differential Equations (ICLR 2021): section 3 (VE and VP SDEs, the reverse SDE), section 4.3 (the probability-flow ODE), Appendix D (the family of samplers).
• Song, Meng, Ermon, Denoising Diffusion Implicit Models (ICLR 2021): eq. 7 (the non-Markovian forward process), Theorem 1, eq. 12 (sampling), eq. 16 ($\eta$), section 4.3 (the ODE view) and Table 1.
• Karras, Aittala, Aila, Laine, Elucidating the Design Space of Diffusion-Based Generative Models (EDM, 2022): section 3 (eq. 5 time steps, Algorithm 1 Heun, Figure 3 curvature), section 4 (Algorithm 2, the stochastic sampler).
• Lu, Zhou, Bao, Chen, Li, Zhu, DPM-Solver (NeurIPS 2022): section 3 (the exact solution of the linear part), section 4.1 (DPM-Solver-1 = DDIM).
• Anderson, Reverse-time diffusion equation models, Stochastic Processes and their Applications 12 (1982): the original time-reversal result.
• Holderrieth and Erives, An Introduction to Flow Matching and Diffusion Models (MIT 6.S184 lecture notes, 2025): chapters on SDEs, the Fokker–Planck equation and sampling.
• k-diffusion sampling.py: short reference implementations of get_sigmas_karras, sample_euler, sample_heun, sample_euler_ancestral and sample_dpmpp_2m.