Module 4 · Reading 4.3 · read before S25: Diffusion Models

The PDF posted in Canvas is the assigned version of this reading. It is built from this page, so the text is identical; anything here marked as deep reading is optional and is not in the PDF.

The last two sections left us with a dilemma. VAEs train stably but blur; GANs sharpen but train on a knife-edge. Diffusion models resolve the dilemma with a move that sounds almost too crude to work: destroy the data with noise, gradually, over many steps; then train a network to undo the destruction one step at a time. Generation becomes a chain of small denoising problems, and each denoising problem is ordinary supervised regression: stable, well-understood, no adversary anywhere.

This recipe, crystallized in the Denoising Diffusion Probabilistic Models paper (Ho, Jain & Abbeel, NeurIPS 2020), is the engine of essentially every flagship image generator since: Stable Diffusion, DALL-E 2, Imagen, and their video successors. The mathematics can be dressed up in variational bounds or stochastic differential equations, but the working core fits in two short scripts, and both run on your laptop in seconds.

How to read this

The forward process: destroying data on a schedule

Start with what is not learned. The forward process takes a data point \(x_0\) and produces progressively noisier versions \(x_1, \ldots, x_T\) by mixing in Gaussian noise, a little per step, controlled by a noise schedule \(\beta_1, \ldots, \beta_T\) (small numbers, increasing). One step scales the signal down by \(\sqrt{1 - \beta_t}\) and adds noise with variance \(\beta_t\). Because each step is Gaussian, the whole chain collapses into a closed form; you can jump straight from \(x_0\) to any \(x_t\) without simulating the steps in between:

\[x_t = \sqrt{\bar{\alpha}_t}\, x_0 + \sqrt{1 - \bar{\alpha}_t}\, \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, I)\]

where \(\bar{\alpha}_t\) is the running product of \((1 - \beta_s)\). As \(\bar{\alpha}_T\) approaches zero, the endpoint approaches \(\mathcal{N}(0, I)\) regardless of the original data scale; standardizing to a convenient unit scale is not a mathematical requirement for that limit, but it makes the variance-preserving interpretation and optimization well conditioned.

“Driven to nearly zero” depends entirely on the schedule, and this is worth seeing in numbers rather than trusting. The DDPM paper uses T = 1000 steps with \(\beta\) rising linearly from 1e-4 to 0.02. The script below tracks three statistics of \(x_t\) on a two-moons dataset: the signal coefficient \(\sqrt{\bar{\alpha}_t}\), the overall standard deviation, and the correlation between \(x_t\) and the original \(x_0\). It then repeats the measurement for a schedule cut to T = 200: first naively, then rescaled.

import torch

torch.manual_seed(0)

# Two-moons data, generated in-code and standardized
def two_moons(n):
    t = torch.rand(n) * torch.pi
    outer = torch.stack([torch.cos(t[: n // 2]), torch.sin(t[: n // 2])], dim=1)
    inner = torch.stack([1 - torch.cos(t[n // 2 :]), 0.5 - torch.sin(t[n // 2 :])], dim=1)
    x = torch.cat([outer, inner]) + 0.05 * torch.randn(n, 2)
    return (x - x.mean(0)) / x.std(0)

x0 = two_moons(4000)

def report(T, beta_max, ts):
    betas = torch.linspace(1e-4, beta_max, T)
    alpha_bar = torch.cumprod(1 - betas, dim=0)
    print(f"T = {T}, betas 1e-4 .. {beta_max}")
    print(f"{'t':>4} {'sqrt(a_bar)':>12} {'std(x_t)':>9} {'corr(x_t, x0)':>14}")
    for t in ts:
        eps = torch.randn_like(x0)
        xt = alpha_bar[t].sqrt() * x0 + (1 - alpha_bar[t]).sqrt() * eps
        corr = torch.corrcoef(torch.stack([xt[:, 0], x0[:, 0]]))[0, 1]
        print(f"{t:>4} {alpha_bar[t].sqrt().item():>12.4f} "
              f"{xt.std().item():>9.4f} {corr.item():>14.4f}")

# The DDPM paper's schedule: T=1000, betas from 1e-4 to 0.02
report(1000, 0.02, [0, 249, 499, 749, 999])
print()
# Same beta range but only 200 steps: noising never finishes
report(200, 0.02, [0, 99, 199])
print()
# 200 steps with beta_max rescaled: noising completes again
report(200, 0.05, [0, 99, 199])

Running this prints:

T = 1000, betas 1e-4 .. 0.02
   t  sqrt(a_bar)  std(x_t)  corr(x_t, x0)
   0       0.9999    0.9998         0.9999
 249       0.7239    0.9988         0.7192
 499       0.2803    0.9973         0.2891
 749       0.0579    0.9998         0.0432
 999       0.0064    0.9902        -0.0082

T = 200, betas 1e-4 .. 0.02
   t  sqrt(a_bar)  std(x_t)  corr(x_t, x0)
   0       0.9999    0.9998         1.0000
  99       0.7762    1.0014         0.7770
 199       0.3636    1.0120         0.3533

T = 200, betas 1e-4 .. 0.05
   t  sqrt(a_bar)  std(x_t)  corr(x_t, x0)
   0       0.9999    1.0000         1.0000
  99       0.5321    1.0119         0.5261
 199       0.0782    1.0049         0.0720

The first table shows the schedule driving the signal coefficient near zero. Because this toy data was standardized to unit variance, std(x_t) stays near 1.0; for non-unit-scale data the marginal variance transitions toward 1 instead of remaining exactly constant.

The second table is the cautionary one. Keep the same \(\beta\) range but run only 200 steps, and the endpoint still carries a 0.35 correlation with the data: the “pure noise” you would start sampling from is not pure at all, and samples inherit a systematic bias toward the training data’s large-scale layout. The third table shows the repair: raise beta_max to 0.05 so that the same total destruction fits in 200 steps (endpoint correlation 0.07). The lesson generalizes: T and the \(\beta\) range are one design decision, not two, and any change to the step count must be checked against the statistics of \(x_T\).

Historical note. The core idea predates its fame by five years. Sohl-Dickstein and colleagues published “Deep Unsupervised Learning using Nonequilibrium Thermodynamics” at ICML 2015 (diffusion models in recognizable form, inspired by physics), and the paper drew modest attention while GANs absorbed the field’s energy. A separate line by Song and Ermon (NeurIPS 2019) reached a nearly equivalent method through score matching — learning the gradient of the data density. DDPM in 2020 connected the threads and demonstrated image quality competitive with GANs, and within two years diffusion had taken over. It is a standing reminder that “didn’t beat the benchmark in year one” and “wrong idea” are different statements.

The reverse process: regression on noise

Generation runs the corruption backwards: start from \(x_T \sim \mathcal{N}(0, I)\) and repeatedly apply a learned denoising step. The network’s job at each step is to look at a noisy \(x_t\) (and the step index \(t\)) and estimate what noise was mixed in. DDPM’s central simplification is the parameterization: rather than predicting the clean \(x_0\) or the mean of \(x_{t-1}\), the network \(\varepsilon_\theta(x_t, t)\) predicts the noise vector \(\varepsilon\) itself, trained with plain mean squared error:

\[L = \mathbb{E}\bigl[\, \lVert \varepsilon - \varepsilon_\theta(x_t, t) \rVert^2 \,\bigr]\]

Why predict the noise rather than the clean point? The parameterizations are algebraically convertible, but their loss weighting and numerical behavior differ across timesteps. DDPM found epsilon prediction effective in its setup; modern systems also use \(x_0\) and velocity (\(v\)) parameterizations successfully. Treat target choice and weighting as part of the recipe, not an inherent ranking that noise prediction is always easier.

Training is startlingly mundane. Sample a data point, sample a random step \(t\), sample noise, form \(x_t\) with the closed-form mix, and regress the network’s output onto the noise you injected. Every minibatch is a fresh supervised problem with known ground truth — no adversary whose improvement degrades your gradient signal, no game to balance, none of the fragility from the GAN section. One subtlety deserves a flag: the loss will not go to zero, and should not. At small \(t\), the injected noise is tiny next to the data’s own variation and cannot be separated from it, even in principle, so the optimal predictor has substantial residual error there; at large \(t\), \(x_t\) is almost pure noise and the noise is easy to predict. A plateauing diffusion loss usually means the model has converged, not stalled: judge it by its samples.

At sampling time, this toy uses the common DDPM “fixed large” variance choice:

\[x_{t-1} = \frac{1}{\sqrt{\alpha_t}} \left( x_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}}\, \hat{\varepsilon} \right) + \sqrt{\beta_t}\, z\]

with fresh noise \(z \sim \mathcal{N}(0, I)\) except at the last step. Other valid DDPM variance choices use the posterior variance, and learned variance is possible. Noise is required by this ancestral update; deterministic DDIM uses a different update and is not obtained by simply deleting \(z\) from this formula.

A tiny diffusion model, trained and sampled

Here is the whole pipeline on two moons: a small MLP denoiser conditioned on \(t\) through a handful of sinusoidal features, trained for 4,000 steps, then run backwards from pure noise. Because the data is 2D, we can compare the sampled distribution to the real one with plain statistics instead of eyeballing images: per-dimension means and standard deviations, the correlation between the two coordinates (the moons’ interlocking shape makes it strongly negative), and the average distance from each sample to its nearest data point.

import torch
import torch.nn as nn

torch.manual_seed(0)

def two_moons(n):
    t = torch.rand(n) * torch.pi
    outer = torch.stack([torch.cos(t[: n // 2]), torch.sin(t[: n // 2])], dim=1)
    inner = torch.stack([1 - torch.cos(t[n // 2 :]), 0.5 - torch.sin(t[n // 2 :])], dim=1)
    x = torch.cat([outer, inner]) + 0.05 * torch.randn(n, 2)
    return (x - x.mean(0)) / x.std(0)

data = two_moons(4000)

T = 200
betas = torch.linspace(1e-4, 0.05, T)   # beta_max rescaled so noising completes in 200 steps
alphas = 1 - betas
alpha_bar = torch.cumprod(alphas, dim=0)

def t_embed(t):  # t: (n,) integer timesteps -> (n, 5) features
    ts = t.float() / T
    return torch.stack([ts, torch.sin(2 * torch.pi * ts), torch.cos(2 * torch.pi * ts),
                        torch.sin(8 * torch.pi * ts), torch.cos(8 * torch.pi * ts)], dim=1)

net = nn.Sequential(nn.Linear(7, 128), nn.ReLU(), nn.Linear(128, 128), nn.ReLU(),
                    nn.Linear(128, 2))          # predicts the noise eps
opt = torch.optim.Adam(net.parameters(), lr=1e-3)

for step in range(4000):
    x0 = data[torch.randint(0, len(data), (256,))]
    t = torch.randint(0, T, (256,))
    eps = torch.randn_like(x0)
    xt = alpha_bar[t, None].sqrt() * x0 + (1 - alpha_bar[t, None]).sqrt() * eps
    pred = net(torch.cat([xt, t_embed(t)], dim=1))
    loss = ((pred - eps) ** 2).mean()
    opt.zero_grad(); loss.backward(); opt.step()
    if step % 1000 == 0 or step == 3999:
        print(f"step {step:4d}  denoising loss {loss.item():.4f}")

# Ancestral sampling: start from pure noise, denoise step by step
with torch.no_grad():
    x = torch.randn(2000, 2)
    for t in reversed(range(T)):
        tt = torch.full((2000,), t)
        eps_hat = net(torch.cat([x, t_embed(tt)], dim=1))
        x = (x - betas[t] / (1 - alpha_bar[t]).sqrt() * eps_hat) / alphas[t].sqrt()
        if t > 0:
            x = x + betas[t].sqrt() * torch.randn_like(x)

def stats(name, s):
    corr = torch.corrcoef(s.T)[0, 1]
    print(f"{name:>10}  mean ({s[:, 0].mean():+.4f}, {s[:, 1].mean():+.4f})  "
          f"std ({s[:, 0].std():.4f}, {s[:, 1].std():.4f})  corr(x1,x2) {corr:+.4f}")

stats("data", data)
stats("samples", x)
stats("N(0, I)", torch.randn(2000, 2))

# Average distance from each sample to its nearest data point
holdout = two_moons(2000)
for name, s in [("held-out data", holdout), ("model samples", x),
                ("N(0, I) noise", torch.randn(2000, 2))]:
    d = torch.cdist(s, data).min(dim=1).values.mean()
    print(f"mean nearest-data distance, {name:>14}: {d.item():.4f}")

Running this prints:

step    0  denoising loss 0.9931
step 1000  denoising loss 0.3516
step 2000  denoising loss 0.3581
step 3000  denoising loss 0.3481
step 3999  denoising loss 0.3213
      data  mean (+0.0000, +0.0000)  std (1.0000, 1.0000)  corr(x1,x2) -0.4531
   samples  mean (+0.0429, -0.0370)  std (0.9788, 1.0005)  corr(x1,x2) -0.4410
   N(0, I)  mean (+0.0095, +0.0371)  std (0.9841, 0.9846)  corr(x1,x2) -0.0200
mean nearest-data distance,  held-out data: 0.0162
mean nearest-data distance,  model samples: 0.0480
mean nearest-data distance,  N(0, I) noise: 0.2544

The training trace behaves exactly as promised: the loss falls from 0.99 (the variance of unit noise, which is what an untrained predictor scores) to about 0.33 within the first thousand steps and then plateaus, because the remainder is irreducible. The interesting rows are the comparisons. Started from noise with zero correlation between coordinates, the sampler produces points whose coordinate correlation is -0.44 against the data’s -0.45, with means and standard deviations close to the data’s. The nearest-neighbor metric sharpens the picture: model samples land an average of 0.048 from the data manifold — five times closer than Gaussian noise (0.254), but three times farther than genuinely held-out data (0.016). Sized to what this toy can show: the model has clearly learned the shape of the distribution, and equally clearly generates a somewhat thicker, fuzzier version of it than nature does. Both facts came from a table, not a scatter plot (a habit worth keeping when the data is too high-dimensional to plot).

People. Two graduate students anchor this story. Jonathan Ho was a Berkeley PhD student (advised by Pieter Abbeel) when he led the DDPM paper in 2020, and stayed central to the line — he later co-authored classifier-free guidance and Google’s Imagen and video-diffusion work. Yang Song, then a Stanford PhD student with Stefano Ermon, built the score-based route: his 2019 generative score-matching paper and the 2021 SDE framework unified diffusion and score models into one continuous-time picture, work recognized with an ICLR 2021 outstanding paper award.

Sampling in fewer steps: DDIM

The toy sampler above calls the network 200 times to produce one batch of samples; the DDPM paper’s models call it 1,000 times. That cost is diffusion’s structural weakness (a GAN spends one forward pass per sample), and the first big attack on it came from a change of perspective rather than a change of training. Denoising Diffusion Implicit Models (DDIM; Song, Meng & Ermon, ICLR 2021) observed that the training objective you just used never actually commits you to the ancestral sampler. The same trained \(\varepsilon_\theta\) supports a whole family of reverse processes, including a deterministic one: at each step, use the noise estimate to form a best guess of the clean data, \(\hat{x}_0 = (x_t - \sqrt{1 - \bar{\alpha}_t}\, \hat{\varepsilon}) / \sqrt{\bar{\alpha}_t}\), then re-noise that guess analytically down to the next timestep you plan to visit — which no longer needs to be the adjacent one. No fresh noise is injected (the stochastic term gets weight \(\eta = 0\)), so the sampler becomes a deterministic map from \(x_T\) to a sample, and the timestep sequence becomes a free choice: 50 steps, 20, 10.

The script below trains the identical denoiser as before, then compares the ancestral sampler against DDIM at decreasing step budgets. Quality is measured as in the last section: mean nearest-neighbor distance to the training data (held-out data scores about 0.015, the floor) and the correlation between the two coordinates (the moons make it strongly negative).

import torch
import torch.nn as nn

torch.manual_seed(0)

def two_moons(n):
    t = torch.rand(n) * torch.pi
    outer = torch.stack([torch.cos(t[: n // 2]), torch.sin(t[: n // 2])], dim=1)
    inner = torch.stack([1 - torch.cos(t[n // 2 :]), 0.5 - torch.sin(t[n // 2 :])], dim=1)
    x = torch.cat([outer, inner]) + 0.05 * torch.randn(n, 2)
    return (x - x.mean(0)) / x.std(0)

data = two_moons(4000)

T = 200
betas = torch.linspace(1e-4, 0.05, T)
alphas = 1 - betas
alpha_bar = torch.cumprod(alphas, dim=0)

def t_embed(t):
    ts = t.float() / T
    return torch.stack([ts, torch.sin(2 * torch.pi * ts), torch.cos(2 * torch.pi * ts),
                        torch.sin(8 * torch.pi * ts), torch.cos(8 * torch.pi * ts)], dim=1)

net = nn.Sequential(nn.Linear(7, 128), nn.ReLU(), nn.Linear(128, 128), nn.ReLU(),
                    nn.Linear(128, 2))
opt = torch.optim.Adam(net.parameters(), lr=1e-3)

for step in range(4000):
    x0 = data[torch.randint(0, len(data), (256,))]
    t = torch.randint(0, T, (256,))
    eps = torch.randn_like(x0)
    xt = alpha_bar[t, None].sqrt() * x0 + (1 - alpha_bar[t, None]).sqrt() * eps
    pred = net(torch.cat([xt, t_embed(t)], dim=1))
    loss = ((pred - eps) ** 2).mean()
    opt.zero_grad(); loss.backward(); opt.step()

def ddpm_sample(n):
    x = torch.randn(n, 2)
    for t in reversed(range(T)):
        tt = torch.full((n,), t)
        eps_hat = net(torch.cat([x, t_embed(tt)], dim=1))
        x = (x - betas[t] / (1 - alpha_bar[t]).sqrt() * eps_hat) / alphas[t].sqrt()
        if t > 0:
            x = x + betas[t].sqrt() * torch.randn_like(x)
    return x

def ddim_sample(n, n_steps):
    ts = torch.linspace(T - 1, 0, n_steps).long()      # strided timestep subsequence
    x = torch.randn(n, 2)
    for i, t in enumerate(ts):
        tt = torch.full((n,), t.item())
        eps_hat = net(torch.cat([x, t_embed(tt)], dim=1))
        x0_hat = (x - (1 - alpha_bar[t]).sqrt() * eps_hat) / alpha_bar[t].sqrt()
        ab_prev = alpha_bar[ts[i + 1]] if i + 1 < len(ts) else torch.tensor(1.0)
        x = ab_prev.sqrt() * x0_hat + (1 - ab_prev).sqrt() * eps_hat   # eta = 0
    return x

def quality(s):
    d = torch.cdist(s, data).min(dim=1).values.mean().item()
    corr = torch.corrcoef(s.T)[0, 1].item()
    return d, corr

with torch.no_grad():
    holdout = two_moons(2000)
    d, corr = quality(holdout)
    print(f"{'sampler':>14} {'NN dist':>8} {'corr(x1,x2)':>12}  (denoiser calls)")
    print(f"{'held-out data':>14} {d:>8.4f} {corr:>+12.4f}")
    d, corr = quality(ddpm_sample(2000))
    print(f"{'DDPM 200':>14} {d:>8.4f} {corr:>+12.4f}  (200)")
    for k in [50, 20, 10, 5, 2]:
        d, corr = quality(ddim_sample(2000, k))
        print(f"{f'DDIM {k}':>14} {d:>8.4f} {corr:>+12.4f}  ({k})")

Running this prints:

       sampler  NN dist  corr(x1,x2)  (denoiser calls)
 held-out data   0.0152      -0.4166
      DDPM 200   0.0485      -0.4417  (200)
       DDIM 50   0.0482      -0.4485  (50)
       DDIM 20   0.0521      -0.4755  (20)
       DDIM 10   0.0574      -0.4899  (10)
        DDIM 5   0.0900      -0.5854  (5)
        DDIM 2   0.2489      +0.0386  (2)

The table is the whole economics of diffusion sampling in seven rows. DDIM at 50 steps matches the 200-step ancestral sampler (0.0482 versus 0.0485) — a 4x speedup for free. Twenty steps costs a few percent; ten steps is visibly worse but usable; and below five the process falls apart, with 2-step samples farther from the data than they are structured (the correlation even flips sign). The degradation is graceful until it is not, because each DDIM step is effectively an integration step of an underlying continuous-time process: larger strides accumulate integration error slowly, then catastrophically. Production systems live exactly in this trade space: 20–50 steps with refined solvers, or a further distillation stage that trains a new network to jump the whole trajectory in 1–4 steps (progressive distillation and consistency models are the marquee versions, and the adversarial variant appeared at the end of the GAN section).

Deterministic DDIM enables an approximate inversion procedure, but determinism alone does not make a discrete map bijective. Reversing the numerical trajectory with an imperfect denoiser accumulates error, and classifier-free guidance can enlarge it. Editing pipelines use DDIM inversion and stronger optimization/null-text variants precisely because reconstruction is not automatically exact.

Steering the denoiser: classifier-free guidance

Everything so far generates unconditionally — you get a sample, not the sample you asked for. The mechanism that turns a denoiser into a promptable system is conditioning, and the specific trick that makes conditioning strong is classifier-free guidance (CFG), now as much a part of the standard recipe as the noise schedule. The idea arrived in two stages. Dhariwal and Nichol (NeurIPS 2021) steered a trained denoiser with a separately trained classifier (nudge \(x_t\) in the direction that makes the classifier more confident of the target class), which worked but required an extra classifier robust to noisy inputs. Ho and Salimans (“Classifier-Free Diffusion Guidance,” 2022) eliminated the classifier with a training-time trick: randomly drop the conditioning signal (10–20% of batches see a null token), so one network learns both the conditional predictor \(\varepsilon(x_t, c)\) and the unconditional \(\varepsilon(x_t)\). At sampling time, run both and extrapolate past the conditional prediction:

\[\varepsilon_{\text{guided}} = \varepsilon_{\text{uncond}} + s\,(\varepsilon_{\text{cond}} - \varepsilon_{\text{uncond}})\]

At \(s = 1\) this is plain conditioning; at \(s \gt 1\) it amplifies exactly the component of the prediction that the condition is responsible for. The name “classifier-free” is precise: by Bayes’ rule, the difference between the conditional and unconditional scores is the gradient of an implicit classifier’s log-probability, \(\nabla_x \log p(c \mid x_t)\). That is the same quantity the earlier classifier-guidance method computed with a separate network, now obtained by subtracting two outputs of one network. Scaling that difference by \(s\) sharpens the implicit classifier’s influence, which is why guidance behaves like turning up a “match the condition” force and why its side effects are classifier-shaped: samples migrate toward whatever the model finds most unambiguously conditioned. The script conditions our two-moons denoiser on which moon a sample belongs to, drops the label 15% of the time in training, and then sweeps the guidance scale, asking for the outer moon. Three statistics per scale: purity (fraction of samples whose nearest data point is on the requested moon), nearest-neighbor distance to that moon, and spread (average per-coordinate standard deviation, a diversity proxy; the outer moon’s own spread is 0.73).

import torch
import torch.nn as nn

torch.manual_seed(0)

def two_moons_labeled(n):
    t = torch.rand(n) * torch.pi
    outer = torch.stack([torch.cos(t[: n // 2]), torch.sin(t[: n // 2])], dim=1)
    inner = torch.stack([1 - torch.cos(t[n // 2 :]), 0.5 - torch.sin(t[n // 2 :])], dim=1)
    x = torch.cat([outer, inner]) + 0.05 * torch.randn(n, 2)
    y = torch.cat([torch.zeros(n // 2), torch.ones(n // 2)])   # 0 = outer, 1 = inner
    mean, std = x.mean(0), x.std(0)
    return (x - mean) / std, y

data, labels = two_moons_labeled(4000)

T = 200
betas = torch.linspace(1e-4, 0.05, T)
alphas = 1 - betas
alpha_bar = torch.cumprod(alphas, dim=0)

def t_embed(t):
    ts = t.float() / T
    return torch.stack([ts, torch.sin(2 * torch.pi * ts), torch.cos(2 * torch.pi * ts),
                        torch.sin(8 * torch.pi * ts), torch.cos(8 * torch.pi * ts)], dim=1)

# Conditioning channel: +1 = outer moon, -1 = inner moon, 0 = "no condition" (dropped)
net = nn.Sequential(nn.Linear(8, 128), nn.ReLU(), nn.Linear(128, 128), nn.ReLU(),
                    nn.Linear(128, 2))
opt = torch.optim.Adam(net.parameters(), lr=1e-3)

for step in range(4000):
    idx = torch.randint(0, len(data), (256,))
    x0, y = data[idx], labels[idx]
    c = torch.where(y == 0, 1.0, -1.0)
    c = c * (torch.rand(256) > 0.15).float()      # drop the label 15% of the time
    t = torch.randint(0, T, (256,))
    eps = torch.randn_like(x0)
    xt = alpha_bar[t, None].sqrt() * x0 + (1 - alpha_bar[t, None]).sqrt() * eps
    pred = net(torch.cat([xt, t_embed(t), c[:, None]], dim=1))
    loss = ((pred - eps) ** 2).mean()
    opt.zero_grad(); loss.backward(); opt.step()

def sample_guided(n, scale, target=1.0):
    """DDPM sampling with classifier-free guidance toward c = target."""
    x = torch.randn(n, 2)
    for t in reversed(range(T)):
        tt = torch.full((n,), t)
        feats = t_embed(tt)
        eps_c = net(torch.cat([x, feats, torch.full((n, 1), target)], dim=1))
        eps_u = net(torch.cat([x, feats, torch.zeros(n, 1)], dim=1))
        eps_hat = eps_u + scale * (eps_c - eps_u)
        x = (x - betas[t] / (1 - alpha_bar[t]).sqrt() * eps_hat) / alphas[t].sqrt()
        if t > 0:
            x = x + betas[t].sqrt() * torch.randn_like(x)
    return x

outer_pts = data[labels == 0]
with torch.no_grad():
    print(f"{'scale':>6} {'purity':>7} {'NN dist':>8} {'spread':>7}")
    for scale in [0.0, 1.0, 2.0, 4.0, 8.0]:
        s = sample_guided(2000, scale, target=1.0)   # condition: outer moon
        d_all = torch.cdist(s, data)
        nearest = d_all.argmin(dim=1)
        purity = (labels[nearest] == 0).float().mean().item()
        nn_dist = torch.cdist(s, outer_pts).min(dim=1).values.mean().item()
        spread = s.std(dim=0).mean().item()
        print(f"{scale:>6.1f} {purity:>7.3f} {nn_dist:>8.4f} {spread:>7.4f}")
    ref = outer_pts.std(dim=0).mean().item()
    print(f"outer-moon data spread: {ref:.4f}")

Running this prints:

 scale  purity  NN dist  spread
   0.0   0.468   0.4389  1.0041
   1.0   0.997   0.0232  0.7140
   2.0   1.000   0.0161  0.6390
   4.0   1.000   0.0140  0.5776
   8.0   1.000   0.0137  0.5155
outer-moon data spread: 0.7294

Scale 0 is the unconditional model: samples split between the moons (purity 0.468, spread 1.00, the whole dataset). Scale 1, plain conditioning, is already nearly perfect here: purity 0.997, and the spread (0.714) matches the real outer moon’s (0.729); the model reproduces the conditional distribution faithfully. The interesting rows are above 1. Guidance pushes purity to 1.000 and pulls samples tighter onto the moon (nearest-neighbor distance falls from 0.023 to 0.014) — but the spread keeps shrinking past the data’s own value, down to 0.52 at scale 8, 30% below reality. That is the CFG trade-off in miniature and it is not a toy artifact: guided samples are more typical of the condition than real data is, concentrating in the distribution’s safest regions. In image models this appears as higher fidelity and prompt adherence at the price of diversity, and eventually as oversaturated, “airbrushed” outputs at extreme scales. Production text-to-image systems default to scales around 5–10, chosen by exactly the kind of sweep you just ran.

Pitfall. Guidance changes the effective distribution you sample from, so any evaluation of a guided model is an evaluation at one particular scale — comparing two models sampled at different guidance scales confounds the model with the knob. When you report generative metrics (or compare against a paper), fix the sampler, step count, and guidance scale, and state them; a “better” FID achieved by tuning guidance is a different claim than a better model.

From two moons to Stable Diffusion

Scaling this toy to photographs changes the components but not the logic. The MLP becomes a U-Net or transformer; the five-feature time embedding becomes a proper sinusoidal embedding; the schedule gets refinements like the cosine schedule of Nichol and Dhariwal (ICML 2021), which fixed the linear schedule’s habit of destroying information too early. Three additions carry most of the distance to a product you have used.

Two of the three you have now built yourself: DDIM-style fast sampling and classifier-free guidance are toy versions of production components. The third, latent diffusion (Rombach et al., CVPR 2022), runs denoising in a compressed autoencoder latent, greatly reducing spatial compute. This explains Stable Diffusion’s computational efficiency and architectural lineage; “Stable” is a product/model name from Stability AI, not a technical abbreviation for the latent space.

The trade-off table you should carry out of these three sections: VAEs give stable training and a usable latent space but blur; GANs give sharp one-shot generation but unstable training and mode dropping; diffusion gives stable training and excellent coverage and fidelity, but pays with slow, iterative sampling. Modern systems mix the ingredients (a VAE for the latent space, diffusion for generation, sometimes adversarial distillation for speed), which is why you needed all three sections, not just the newest one.

Pitfall. The most common way to break a working diffusion implementation is changing T (usually downward, to make sampling faster) while keeping the \(\beta\) schedule fixed. Nothing errors: training still converges, sampling still runs. But as the middle table in the first example shows, the forward process no longer reaches noise (corr(x_T, x0) = 0.35 in our truncated run), while the sampler still starts from \(\mathcal{N}(0, I)\), a distribution the model was never trained to denoise from. The resulting bias is quiet and systematic. Whenever you touch T or the schedule, reprint the statistics of \(x_T\); it is a three-line check.

Going deeper

One process, three languages. The discrete noising chain you trained is one description of diffusion; the field now fluently uses two others, and reading current papers requires recognizing all three. The score language says the network learns \(\nabla_x \log p(x_t)\), the direction of increasing data density at each noise level; the epsilon-predictor is that score in disguise (they differ by a known scaling). The differential equation language, unified by Song and colleagues (ICLR 2021), takes the step count to infinity: the forward process becomes a stochastic differential equation, and generation becomes solving either a reverse SDE (the DDPM-like stochastic sampler) or a deterministic probability-flow ODE (the DDIM-like sampler) — which is why fewer sampling steps means larger integration error, exactly the degradation your DDIM table measured. Once sampling is an ODE solve, the numerical-methods toolbox opens: higher-order solvers (Heun, DPM-Solver-class methods) buy accuracy per step, which is how modern systems make 20 steps look like 200.

Flow matching: a major current formulation. A newer language asks: if what we ultimately train is a network defining an ODE from noise to data, why take the detour through noising chains at all? Flow matching (Lipman et al., ICLR 2023) trains the vector field directly along a chosen probability path. Rectified-flow variants often use straight interpolants, which can make low-step numerical integration attractive. Stable Diffusion 3 is one prominent example; diffusion, score, and flow-matching formulations all remain active, so do not read “current” as “the only modern objective.”

Beyond images. The recipe’s indifference to what \(x\) is has made it a general-purpose engine. Audio diffusion arrived early. Video models treat time as one more dimension of the tensor. Protein-structure generation adapted the framework to 3D coordinates with rotational symmetries baked in (RFdiffusion, published in Nature in 2023, designs novel proteins with a diffusion model over backbone geometry). Robotics uses diffusion to model action distributions (“diffusion policies”), where covering multiple valid action modes (the very property your GAN section showed adversarial training struggles with) is the selling point. And discrete diffusion over text tokens is an active research line challenging autoregressive generation on its home turf. The through-line: wherever “average of two good answers is a bad answer” describes the output space, mode-covering iterative refinement is a natural fit.

What the simple story leaves out. Two honest gaps between this section and practice. First, our MLP scales to two dimensions; real systems spend most of their engineering on the denoiser architecture: U-Nets with attention at multiple resolutions, and increasingly pure transformers over latent patches (the DiT line of Peebles & Xie, ICCV 2023, whose scaling behavior echoes the scaling laws story). Second, the loss the toy uses weights all timesteps equally, but which timesteps matter most for perceptual quality is a real design axis — schedule and loss-weighting choices (as in the SD3 paper’s timestep-sampling analysis) measurably move image quality in ways the toy cannot show.

Common pitfalls

  • The training loss plateaus at a large value and you keep tuning to push it down. Cause: the epsilon-prediction objective has a large irreducible floor — much of the injected noise is unrecoverable in principle. Fix: evaluate by sampling and comparing distribution statistics, not by the absolute loss value.
  • Samples are biased toward the data’s global layout even from “pure noise” starts. Cause: schedule/T mismatch — the forward process never fully destroys the data, but sampling starts from \(\mathcal{N}(0, I)\) anyway. Fix: check sqrt(a_bar_T) and corr(x_T, x0); rescale beta_max when you change T.
  • An ancestral DDPM implementation gives over-smoothed samples. Cause: its stochastic term may have been omitted or the wrong variance used. Fix: implement the selected ancestral posterior/variance exactly; do not mix its coefficients with deterministic DDIM.
  • Training diverges or starts at an enormous loss. Cause: data scale may be poorly matched to the schedule/model parameterization. Fix: normalize to a documented convenient range or unit scale and invert that preprocessing after sampling; it is an optimization convention, not an endpoint theorem.
  • A fast sampler built by striding the ancestral update produces degraded samples even at generous step counts. Cause: the DDPM update’s coefficients are only correct for adjacent timesteps; skipping steps inside it is integrating with the wrong step size. Fix: use the DDIM update (or a proper ODE solver), which is derived for arbitrary timestep subsequences; our 50-step DDIM matches the 200-step ancestral sampler.
  • Conditional samples score beautifully on adherence but users complain everything looks the same. Cause: guidance scale set high — the sampler concentrates on the condition’s most typical region, shrinking diversity below the data’s own spread (measured directly in our sweep). Fix: sweep the scale and report both an adherence metric and a diversity metric; pick the knee, not the max.

Try it yourself

  1. Rerun the sampler with the noise-injection line (x = x + betas[t].sqrt() * ...) deleted, and recompute the comparison table. Which statistics degrade, and in which direction?
  2. Implement the strided sampler: denoise using only every 4th timestep (50 steps instead of 200), reusing the trained network. How much does the nearest-data distance degrade? This is the intuition behind DDIM-style fast samplers.
  3. Replace the linear schedule with a cosine one (\(\bar{\alpha}_t = \cos^2\bigl((t/T + 0.008)/1.008 \cdot \pi/2\bigr)\), from Nichol & Dhariwal) and reprint the forward-process table. How does the destruction profile differ at small \(t\)?
  4. Compute the nearest-data distance separately for samples whose nearest data point is on the upper moon vs. the lower moon. Is the model equally good at both — the diffusion analogue of the mode-coverage check from the GAN section?
  5. Add the stochastic term back to DDIM (eta > 0): after computing x0_hat, mix in fresh noise scaled by eta * sqrt(1 - ab_prev) and reduce the deterministic noise share to match. Rebuild the step-count table at eta = 0.5 and 1.0. Which step budgets are hurt most by stochasticity?
  6. In the guidance experiment, sample with scale = -1 (guidance away from the outer moon) and compute purity with respect to each moon. Relate what you find to negative prompting in image generators.

Further reading

  • Ho, Jain & Abbeel, “Denoising Diffusion Probabilistic Models,” NeurIPS 2020 — the paper this section’s toy implements: epsilon-prediction, the training loop, and the ancestral sampler. arxiv.org/abs/2006.11239
  • Sohl-Dickstein et al., “Deep Unsupervised Learning using Nonequilibrium Thermodynamics,” ICML 2015 — the physics-inspired original, five years ahead of its fame. arxiv.org/abs/1503.03585
  • Song & Ermon, “Generative Modeling by Estimating Gradients of the Data Distribution,” NeurIPS 2019 — the score-matching route to the same method. arxiv.org/abs/1907.05600
  • Song et al., “Score-Based Generative Modeling through Stochastic Differential Equations,” ICLR 2021 — the continuous-time unification: reverse SDE, probability-flow ODE, and why samplers are solvers. arxiv.org/abs/2011.13456
  • Song, Meng & Ermon, “Denoising Diffusion Implicit Models” (DDIM), ICLR 2021 — the deterministic, stride-able sampler you built above. arxiv.org/abs/2010.02502
  • Ho & Salimans, “Classifier-Free Diffusion Guidance,” 2022 — the two-in-one conditioning trick behind every modern text-to-image sampler. arxiv.org/abs/2207.12598
  • Rombach et al., “High-Resolution Image Synthesis with Latent Diffusion Models,” CVPR 2022 — diffusion moved into an autoencoder’s latent space; the Stable Diffusion architecture. arxiv.org/abs/2112.10752
  • Lipman et al., “Flow Matching for Generative Modeling,” ICLR 2023 — the simulation-free vector-field objective that is displacing the noising-chain formulation. arxiv.org/abs/2210.02747
  • Esser et al., “Scaling Rectified Flow Transformers for High-Resolution Image Synthesis,” ICML 2024 — rectified flow plus a transformer backbone at production scale (Stable Diffusion 3). arxiv.org/abs/2403.03206
  • Prince, Understanding Deep Learning, Chapter 18 (Diffusion models) — a free textbook chapter with the variational derivation this section deliberately skipped. udlbook.github.io/udlbook/

Before class you should be able to

  • Write the closed-form forward step \(x_t = \sqrt{\bar{\alpha}_t}\, x_0 + \sqrt{1 - \bar{\alpha}_t}\, \varepsilon\) and explain what each factor does.
  • Explain why changing T without rescaling the beta schedule silently biases samples, and name the three-line check that catches it.
  • State the DDPM training loss and explain why it is stable where GAN training is not.
  • Explain why a plateauing denoising loss usually means convergence, and what to judge the model by instead.
  • Describe what DDIM changes (deterministic sampler, free choice of timestep subsequence) and roughly how quality degrades as steps shrink.
  • Write the classifier-free guidance combination \(\varepsilon_{\text{uncond}} + s\,(\varepsilon_{\text{cond}} - \varepsilon_{\text{uncond}})\) and describe the fidelity-versus-diversity trade-off as \(s\) grows.
  • Recite the three-way trade-off table: what VAEs, GANs, and diffusion each buy and pay, and where latent diffusion fits.

Connects to

  • This session is covered by Quiz 4 (Module 4: S23–S26), given at S27.
  • Fall Recess is next week — no class; S26 resumes the following Monday.
  • If your semester project involves generation, the comparison-statistics pattern here (held-out reference, noise baseline, nearest-neighbor distance) is directly reusable for your report (due the week of S27).

← S24 · GANs · All readings · S26 · Efficiency and Scale →