Diffusion Model

A diffusion model generates images by learning to undo a gradual corruption. A fixed forward process adds a little Gaussian noise to an image at each of $T$ steps until nothing of the image is left; a neural network is trained to take one step back, from a noisier image to a slightly cleaner one. To generate, the model starts from pure noise and applies the learned reverse process $T$ times. Each reverse step is a form of image denoising, and the network that performs it is usually a U-Net told which step it is at. In 2021 diffusion models overtook GANs in ImageNet sample quality as measured by the Fréchet Inception Distance (FID) (Dhariwal & Nichol, 2021), and latent diffusion became the basis of Stable Diffusion (Rombach et al., 2022).

The generative model article places diffusion among the other families and compares their likelihoods, sample quality and speed. This article covers the ideas shared by every diffusion method and fixes the notation of the diffusion articles that go deeper: DDPM, score matching, diffusion sampling, classifier-free guidance and latent diffusion.

Core Idea

Destroying structure is easy: add enough independent Gaussian noise to any image and it becomes indistinguishable from a sample of $\mathcal{N}(0, I)$. Sohl-Dickstein et al. (2015) noted, citing Feller (1949), that when each noising step is small, the reverse of the step is also approximately Gaussian, so a network only has to predict the mean (and possibly the variance) of a Gaussian at every step. A long chain of such easy problems models a very complicated distribution.

Three facts make this practical. First, the noisy image at any step can be drawn in one shot from the clean image, so training picks a random step for each image. Second, the mean of each reverse step depends on the unknown clean image only through an estimate of it, so training reduces to a regression: predict the noise that was added (Ho, Jain & Abbeel, 2020). Third, that noise prediction is, up to a scale factor, the score, the gradient of the log-density of noisy images, which links diffusion to score-based models (Song & Ermon, 2019) and to differential equations (Song et al., 2021).

Mathematical Formulation

Notation

The diffusion articles on this site follow Ho et al. (2020). Images are scaled to $[-1, 1]$ and treated as vectors in $\mathbb{R}^D$, $D = H \times W \times C$.

Symbol	Meaning
$x_0$	A clean image, $x_0 \sim q(x_0)$, the data distribution
$x_t$, $t = 1, \dots, T$	The image after $t$ noising steps; $x_T$ is close to pure noise
$\beta_t \in (0, 1)$	Variance of the noise added at step $t$ (the noise schedule)
$\alpha_t = 1 - \beta_t$	Fraction of the previous step’s variance kept
$\bar\alpha_t = \prod_{s=1}^{t} \alpha_s$	Fraction of the clean image’s variance left at step $t$
$\epsilon \sim \mathcal{N}(0, I)$	Standard Gaussian noise
$\epsilon_\theta(x_t, t)$	Network predicting the noise in $x_t$
$\hat x_0$	Prediction of the clean image
$p_t(x)$	Density of $x_t$ under the forward process (the noisy-image marginal)
$s_\theta(x_t, t) \approx \nabla_{x_t} \log p_t(x_t)$	Score: gradient of the log-density of noisy images at step $t$
$\mathrm{SNR}(t) = \bar\alpha_t / (1 - \bar\alpha_t)$	Signal-to-noise ratio of $x_t$

Forward process

The forward process is a Markov chain with fixed Gaussian transitions (Ho et al., 2020, Eq. 2):

$$q(x_t \mid x_{t-1}) = \mathcal{N}\big(x_t;\ \sqrt{1 - \beta_t}\, x_{t-1},\ \beta_t I\big), \qquad q(x_{1:T} \mid x_0) = \prod_{t=1}^{T} q(x_t \mid x_{t-1}).$$

The factor $\sqrt{1 - \beta_t}$ shrinks the image as noise is added, so that data with unit variance keeps unit variance. Because each step is a linear map plus independent Gaussian noise, $x_t$ given $x_0$ is Gaussian, with the closed form (Eq. 4)

$$q(x_t \mid x_0) = \mathcal{N}\big(x_t;\ \sqrt{\bar\alpha_t}\, x_0,\ (1 - \bar\alpha_t) I\big), \qquad x_t = \sqrt{\bar\alpha_t}\, x_0 + \sqrt{1 - \bar\alpha_t}\, \epsilon.$$

The mean follows by multiplying the factors $\sqrt{\alpha_s}$. For the variance, if $x_{t-1}$ has variance $1 - \bar\alpha_{t-1}$ around its mean, one more step gives $\alpha_t (1 - \bar\alpha_{t-1}) + \beta_t = 1 - \alpha_t \bar\alpha_{t-1} = 1 - \bar\alpha_t$. The schedule is chosen so that $\bar\alpha_T \approx 0$, making $x_T$ essentially $\mathcal{N}(0, I)$ whatever the image was; the first example checks the closed form numerically.

Reverse process

The generative model runs the chain backward from $p(x_T) = \mathcal{N}(0, I)$ with learned Gaussian transitions (Eq. 1):

$$p_\theta(x_{0:T}) = p(x_T) \prod_{t=1}^{T} p_\theta(x_{t-1} \mid x_t), \qquad p_\theta(x_{t-1} \mid x_t) = \mathcal{N}\big(x_{t-1};\ \mu_\theta(x_t, t),\ \Sigma_\theta(x_t, t)\big).$$

The noisy images $x_1, \dots, x_T$ are latent variables with the same dimension as the image, and the forward process plays the role of a fixed encoder with no parameters. Training therefore maximizes an evidence lower bound on $\log p_\theta(x_0)$, which decomposes into KL divergences between Gaussians (the DDPM article derives it). Its targets are the forward-process posteriors $q(x_{t-1} \mid x_t, x_0) = \mathcal{N}(\tilde\mu_t(x_t, x_0), \tilde\beta_t I)$, with $\tilde\beta_t = \frac{1 - \bar\alpha_{t-1}}{1 - \bar\alpha_t} \beta_t$ (Eqs. 6–7). Ho et al. fixed the variance to $\Sigma_\theta = \sigma_t^2 I$ with $\sigma_t^2 = \beta_t$ or $\tilde\beta_t$, which gave similar results. Nichol and Dhariwal (2021) later learned it, which let them sample with an order of magnitude fewer steps at little cost in quality.

Three parameterizations

The network could output $\mu_\theta$ directly, but Ho et al. found it better to predict the noise. Given $x_t$, any one of the following determines the others through the closed form, so they are three readings of the same denoiser:

• Noise, $\epsilon_\theta(x_t, t)$, trained with the simple loss $\mathbb{E}_{x_0, t, \epsilon} \lVert \epsilon - \epsilon_\theta(\sqrt{\bar\alpha_t}\, x_0 + \sqrt{1 - \bar\alpha_t}\, \epsilon,\ t) \rVert^2$ (Ho et al., Eq. 14). The reverse mean becomes $\mu_\theta = \frac{1}{\sqrt{\alpha_t}}\big(x_t - \frac{\beta_t}{\sqrt{1 - \bar\alpha_t}} \epsilon_\theta(x_t, t)\big)$ (Eq. 11).
• Clean image, $\hat x_0 = \big(x_t - \sqrt{1 - \bar\alpha_t}\, \epsilon_\theta(x_t, t)\big) / \sqrt{\bar\alpha_t}$. Ho et al. report that predicting $x_0$ directly gave worse samples in their early experiments.
• Score, $s_\theta(x_t, t) = -\epsilon_\theta(x_t, t) / \sqrt{1 - \bar\alpha_t}$.

The score reading follows from the closed form: $\nabla_{x_t} \log q(x_t \mid x_0) = -(x_t - \sqrt{\bar\alpha_t}\, x_0)/(1 - \bar\alpha_t) = -\epsilon / \sqrt{1 - \bar\alpha_t}$. Vincent (2011) showed that regressing onto this conditional score, denoising score matching, is equivalent to matching the score of the noisy marginal $p_t$ itself. In our reasoning, the minimizer of the squared error is the conditional mean, so the optimal $\epsilon_\theta$ is $\mathbb{E}[\epsilon \mid x_t]$, the optimal $\hat x_0$ is $\mathbb{E}[x_0 \mid x_t]$, and the optimal score is the true $\nabla \log p_t(x_t)$; Salimans and Ho (2022, Eq. 12, arXiv v2) use the same relation, $\nabla \log p_t(x_t) \approx (\sqrt{\bar\alpha_t}\, \hat x_0 - x_t)/(1 - \bar\alpha_t)$ in our notation. The second code example verifies these identities numerically.

The parameterizations differ in how the loss weights the noise levels. Since $\epsilon - \epsilon_\theta = \sqrt{\mathrm{SNR}(t)}\,(\hat x_0 - x_0)$,

$$\lVert \epsilon - \epsilon_\theta \rVert^2 = \mathrm{SNR}(t)\, \lVert x_0 - \hat x_0 \rVert^2,$$

so the noise loss is an image-reconstruction loss weighted by the signal-to-noise ratio, a weighting that, as Salimans and Ho (2022) point out, gives zero weight at zero SNR. Compared with the true variational bound, Ho et al. describe their simple loss as down-weighting the terms for small noise levels, so that the network concentrates on the harder, noisier steps; it gave their best samples. Salimans and Ho also introduced v-prediction, a target that mixes the two, $v = \sqrt{\bar\alpha_t}\, \epsilon - \sqrt{1 - \bar\alpha_t}\, x_0$, from which $\hat x_0 = \sqrt{\bar\alpha_t}\, x_t - \sqrt{1 - \bar\alpha_t}\, \hat v$. It was one of the parameterizations they found to work well with their distillation method, which samples in very few steps.

The score view and SDEs

Song and Ermon (2019) reached the same models from the other side: they trained one network to estimate the score at a geometric sequence of noise scales by denoising score matching and sampled with annealed Langevin dynamics, from the largest scale to the smallest. Noise makes the score well defined, since images lie near a low-dimensional set off which the score of clean data is undefined. Ho et al. (2020) noted that their loss resembles denoising score matching over many noise levels, and their sampler Langevin dynamics.

Song et al. (2021) unified the two in continuous time. The forward process becomes a stochastic differential equation (SDE), $dx = f(x, t)\, dt + g(t)\, dw$, whose time reversal depends on the data only through the score:

$$dx = \big[f(x, t) - g(t)^2\, \nabla_x \log p_t(x)\big]\, dt + g(t)\, d\bar w.$$

The DDPM chain becomes the variance-preserving SDE, $f = -\tfrac{1}{2}\beta(t)\, x$ and $g = \sqrt{\beta(t)}$, and Song and Ermon’s noise scales the variance-exploding SDE; DDPM’s ancestral sampler is one discretization of the reverse SDE. A deterministic probability flow ODE, $dx = [f - \tfrac{1}{2} g^2 \nabla_x \log p_t(x)]\, dt$, has the same marginals and gives exact likelihoods. The score matching and diffusion sampling articles develop these views.

Noise schedules and signal-to-noise ratio

The schedule $\beta_1, \dots, \beta_T$ decides how quickly information is destroyed. Ho et al. used $T = 1000$ and a linear schedule from $\beta_1 = 10^{-4}$ to $\beta_T = 0.02$, small relative to data in $[-1, 1]$. Nichol and Dhariwal (2021) found it suboptimal for $64 \times 64$ and $32 \times 32$ images: the end of the forward process is so noisy that skipping up to 20% of the reverse steps barely hurt FID. Their cosine schedule defines $\bar\alpha_t$ directly,

$$\bar\alpha_t = \frac{f(t)}{f(0)}, \qquad f(t) = \cos^2\!\left(\frac{t/T + s}{1 + s} \cdot \frac{\pi}{2}\right), \qquad \beta_t = \min\!\left(1 - \frac{\bar\alpha_t}{\bar\alpha_{t-1}},\ 0.999\right),$$

with offset $s = 0.008$, so that $\bar\alpha_t$ falls linearly in the middle and changes slowly at both ends. Kingma et al. (2021) showed that in continuous time the variational bound depends on the schedule only through the signal-to-noise ratio at its two endpoints, and that the same holds for losses weighted by a function of SNR, which makes $\mathrm{SNR}(t)$ the natural coordinate for comparing schedules. The schedule still affects the variance of the Monte Carlo estimate of the loss, which they reduced by learning it, and, in our reading, it sets the effective weighting of a loss such as Ho et al.’s, which samples $t$ uniformly rather than SNR.

Building Blocks

The denoiser network. The network maps a noisy image and a step index to an output of the same shape, which makes an encoder–decoder with skip connections the natural choice. Ho et al. (2020) used a PixelCNN++-style U-Net with two residual blocks per resolution, four resolutions for $32 \times 32$ images and six for $256 \times 256$, and group normalization. Their CIFAR-10 model had 35.7 million parameters, and their LSUN and CelebA-HQ models 114 million (Appendix B, arXiv v2).

Time embedding. One network serves every step, so it must be told $t$. Ho et al. add a Transformer-style sinusoidal embedding of $t$ (see positional encoding) into each residual block. Dhariwal and Nichol (2021) replaced the addition with adaptive group normalization, which scales and shifts the normalized activations by amounts predicted from the time and class embeddings, and found it improved FID in their ablation (Table 3, arXiv v4).

Attention. A denoiser at high noise must decide the global layout of the image, which convolutions alone see poorly. Ho et al. placed self-attention blocks at the $16 \times 16$ resolution; Dhariwal and Nichol added attention at $32 \times 32$ and $8 \times 8$ and more heads, which improved FID on ImageNet $128 \times 128$. Text-conditioned models add cross-attention (see Variants).

Transformer denoisers. Peebles and Xie (2023) replaced the U-Net with a plain transformer over patches, in the manner of the Vision Transformer, applied to the latent of a pretrained autoencoder. Their Diffusion Transformer (DiT) conditions each block on time and class through adaptive layer normalization initialized so that every block starts as the identity (adaLN-Zero). FID fell steadily as forward-pass compute grew, and their largest model, DiT-XL/2 with classifier-free guidance at scale 1.5, reached an FID of 2.27 on class-conditional ImageNet $256 \times 256$ (Table 2, arXiv v2).

How Data Flows

Training step. Draw an image $x_0$, a step $t$ uniformly from $\{1, \dots, T\}$ and noise $\epsilon$; form $x_t = \sqrt{\bar\alpha_t}\, x_0 + \sqrt{1 - \bar\alpha_t}\, \epsilon$ in one operation; pass $(x_t, t)$ through the denoiser; and take a gradient step on $\lVert \epsilon - \epsilon_\theta(x_t, t) \rVert^2$ (Ho et al., Algorithm 1). No chain is simulated during training.

Sampling. Draw $x_T \sim \mathcal{N}(0, I)$. For $t = T, \dots, 1$, predict the noise, compute the reverse mean and add fresh noise $\sigma_t z$ (except at the last step):

$$x_{t-1} = \frac{1}{\sqrt{\alpha_t}}\left(x_t - \frac{\beta_t}{\sqrt{1 - \bar\alpha_t}}\, \epsilon_\theta(x_t, t)\right) + \sigma_t z, \qquad z \sim \mathcal{N}(0, I)$$

(Algorithm 2). Generating one image costs $T$ network evaluations. Ho et al. observed that large-scale features appear first and details last: with $x_t$ frozen, repeated sampling preserves only large-scale features when $t$ is large and all but fine details when $t$ is small (§4.3, Figure 7). The DDPM article derives this sampler, and the diffusion sampling article covers faster deterministic and higher-order alternatives.

Why It Works

Each reverse step is a regression with a fixed target, so training minimizes a stationary squared error rather than chasing a GAN’s moving target; Dhariwal and Nichol (2021) list distribution coverage, a stationary training objective and easy scalability as the properties that made diffusion attractive. In our reading, the objective is a reweighted likelihood bound, and like maximum likelihood it penalizes missing parts of the data distribution, which a GAN’s discriminator need not do. One shared network acts as a layout generator at high noise and a detail-preserving denoiser at low noise. Finally, many small steps keep each reverse conditional close to Gaussian (Sohl-Dickstein et al., 2015); with large steps the true denoising distribution becomes multimodal, because several clean images are plausible for one noisy image (Xiao, Kreis & Vahdat, 2022).

Training Considerations

• Scaling. Pixels are mapped to $[-1, 1]$ so that the data and the prior $\mathcal{N}(0, I)$ have comparable scales; the schedule is set relative to this range (Ho et al., 2020, §3.3 and §4).
• Number of steps. $T = 1000$ is common in pixel space. Training cost does not grow with $T$; sampling cost does.
• Schedule and weighting. Both decide how much training effort goes to each signal-to-noise ratio; learning the variances also needs a variational term in the loss (Nichol & Dhariwal, 2021).
• Parameterization. DDPM and the latent diffusion objective (Rombach et al., 2022, Eq. 3) predict $\epsilon$; for distillation to very few steps, Salimans and Ho (2022) switched to $\hat x_0$, mixed or $v$ targets, whose implied $\hat x_0$ stays stable at low SNR.
• Compute. Rombach et al. (2022) cite 150 to 1,000 V100 days to train the strongest pixel-space models of Dhariwal and Nichol, and about five days on one A100 GPU to draw 50,000 samples, which motivated latent diffusion.

Variants

Conditioning. A conditional model $\epsilon_\theta(x_t, t, y)$ takes the condition as an extra input. Class labels are embedded like the time step and injected through adaptive normalization (Dhariwal & Nichol, 2021). Text features enter the U-Net through cross-attention, with queries from the image and keys and values from the text; spatially aligned conditions, such as semantic maps or low-resolution images, are instead concatenated to the input (Rombach et al., 2022).

Guidance. Sampling can be pushed toward the condition at the cost of diversity. Classifier guidance adds the gradient of a classifier trained on noisy images to the score (Dhariwal & Nichol, 2021). Classifier-free guidance trains one network with and without the condition and extrapolates, $\tilde\epsilon = (1 + w)\, \epsilon_\theta(x_t, t, y) - w\, \epsilon_\theta(x_t, t)$ (Ho & Salimans, 2022, Eq. 6, in our notation); the classifier-free guidance article covers it.

Latent diffusion. Rombach et al. (2022) ran diffusion on the latent of a pretrained autoencoder, several times smaller than the image in each spatial direction, which cuts the cost of every step; the latent diffusion article covers it, and the variational autoencoder article the kind of autoencoder.

Fast sampling. DDIM (Song, Meng & Ermon, 2021) samples deterministically in tens of steps from a network trained with the same objective; progressive distillation repeatedly halves the steps of a deterministic sampler and reached an FID of 3.00 on CIFAR-10 with four steps (Salimans & Ho, 2022, Table 2, arXiv v2).

Variance-preserving and variance-exploding models. DDPM shrinks the image as it adds noise; Song and Ermon’s models do not. Kingma et al. (2021) showed that the two are equivalent in continuous time up to a rescaling of $x_t$.

Use in Computer Vision

Diffusion models generate images unconditionally, from class labels and from text. Sohl-Dickstein et al. (2015) already denoised CIFAR-10 images and inpainted a bark texture; Song et al. (2021) performed inpainting and colorization with one unconditional score model by conditioning the reverse SDE; and latent diffusion models reported strong results in inpainting, class-conditional synthesis, text-to-image synthesis and super-resolution (Rombach et al., 2022). As of October 2026, the CompVis repository describes Stable Diffusion as a latent text-to-image diffusion model that builds on latent diffusion, with an 860M-parameter U-Net conditioned on a frozen CLIP ViT-L/14 text encoder. Text-to-image generation, Stable Diffusion, image inpainting and image-to-image translation have their own articles.

Limitations

• Slow sampling. Generation needs many sequential network evaluations; Dhariwal and Nichol (2021) list slower sampling than GANs as a limitation, and the comparison table in the generative model article puts numbers on it.
• Likelihood and sample quality diverge. The loss that gives the best samples is not the one that gives the best likelihood bound (Ho et al., 2020, Table 1); the generative model article discusses the numbers.
• Few-step sampling needs care. The Gaussian reverse step is accurate only for small steps, so cutting steps without a better solver or distillation degrades samples.
• Guidance trades diversity for fidelity. Both classifier and classifier-free guidance improve sample fidelity at the cost of mode coverage (Dhariwal & Nichol, 2021; Ho & Salimans, 2022).

Important Models

Diffusion probabilistic models (2015). Sohl-Dickstein et al. introduced forward and learned reverse diffusion, Gaussian and binomial, with thousands of steps trained on a variational bound.

Noise-conditional score networks (2019). Song and Ermon reported a CIFAR-10 Inception score of 8.87 for their score model.

DDPM (2020). Ho et al. added noise prediction, a U-Net with self-attention and a linear schedule, and obtained an FID of 3.17 and an Inception score of 9.46 on unconditional CIFAR-10, and LSUN $256 \times 256$ samples similar in quality to ProgressiveGAN.

Score SDE (2021). Song et al. unified both lines, added predictor–corrector samplers and the probability-flow ODE, and reported an FID of 2.20 on CIFAR-10 and $1024 \times 1024$ samples.

Improved DDPM and ADM (2021). Nichol and Dhariwal introduced the cosine schedule and learned variances; their improved U-Net (ADM) with classifier guidance surpassed BigGAN-deep on ImageNet, a result the generative adversarial network article discusses.

Latent diffusion (2022) and DiT (2023). Latent diffusion with cross-attention became the basis of Stable Diffusion; DiT replaced the U-Net with a transformer.

Example

The closed form on a shape image, and two schedules. A $16 \times 16$ disk image, $+1$ inside and $-1$ outside, is noised by iterating $q(x_t \mid x_{t-1})$ in 1,000 independent copies. At selected steps the code compares the mean coefficient and variance of the copies with $\sqrt{\bar\alpha_t}$ and $1 - \bar\alpha_t$, then tabulates $\mathrm{SNR}(t)$ for the linear and cosine schedules.

import numpy as np

T = 1000
t_idx = np.arange(1, T + 1)

# Linear schedule of Ho et al. (2020): beta_1 = 1e-4 ... beta_T = 0.02.
beta_lin = np.linspace(1e-4, 0.02, T)

# Cosine schedule of Nichol & Dhariwal (2021): alpha_bar_t = f(t) / f(0), s = 0.008,
# converted to betas and clipped at 0.999.
s = 0.008
f = np.cos((np.arange(T + 1) / T + s) / (1 + s) * np.pi / 2) ** 2
beta_cos = np.minimum(1 - f[1:] / f[:-1], 0.999)

def alpha_bar(beta):
    return np.cumprod(1 - beta)  # alpha_bar[t-1] is the value at step t

ab_lin, ab_cos = alpha_bar(beta_lin), alpha_bar(beta_cos)

# A 16x16 shape image: a disk at +1 on a background of -1 (pixels scaled to [-1, 1]).
yy, xx = np.mgrid[:16, :16]
x0 = np.where((yy - 7.5) ** 2 + (xx - 7.5) ** 2 < 5.5 ** 2, 1.0, -1.0).ravel()

# Run the chain q(x_t | x_{t-1}) step by step for 1,000 independent copies of x0
# and compare the moments of x_t with the closed form N(sqrt(ab_t) x0, (1 - ab_t) I).
rng = np.random.default_rng(0)
N = 1000
x = np.tile(x0, (N, 1))
check = [1, 10, 100, 250, 500, 750, 1000]
print("linear schedule: iterated chain vs closed form q(x_t | x_0)")
print(f"{'t':>5s}{'sqrt(ab_t)':>12s}{'mean coef':>11s}{'1 - ab_t':>10s}{'variance':>10s}")
for t in t_idx:
    b = beta_lin[t - 1]
    x = np.sqrt(1 - b) * x + np.sqrt(b) * rng.standard_normal(x.shape)
    if t in check:
        m = x.mean(0)                       # per-pixel mean over the copies
        coef = m @ x0 / (x0 @ x0)           # least-squares fit of m = c * x0
        var = x.var(0).mean()               # per-pixel variance, averaged
        print(f"{t:5d}{np.sqrt(ab_lin[t-1]):12.4f}{coef:11.4f}{1-ab_lin[t-1]:10.4f}{var:10.4f}")

def snr(ab):
    return ab / (1 - ab)

print("\nsignal-to-noise ratio SNR(t) = ab_t / (1 - ab_t)")
print(f"{'t':>5s}{'linear':>12s}{'cosine':>12s}")
for t in check:
    print(f"{t:5d}{snr(ab_lin[t-1]):12.3g}{snr(ab_cos[t-1]):12.3g}")
for name, ab in [("linear", ab_lin), ("cosine", ab_cos)]:
    t_half = t_idx[np.argmax(snr(ab) < 1)]
    frac = np.mean(snr(ab) < 1e-3)
    print(f"{name}: SNR falls below 1 at t = {t_half}; "
          f"{frac:.0%} of steps have SNR < 0.001; ab_T = {ab[-1]:.2e}")

Output:

linear schedule: iterated chain vs closed form q(x_t | x_0)
    t  sqrt(ab_t)  mean coef  1 - ab_t  variance
    1      0.9999     0.9999    0.0001    0.0001
   10      0.9991     0.9990    0.0019    0.0019
  100      0.9471     0.9472    0.1030    0.1028
  250      0.7239     0.7228    0.4759    0.4762
  500      0.2803     0.2790    0.9214    0.9233
  750      0.0579     0.0552    0.9966    0.9950
 1000      0.0064     0.0086    1.0000    0.9983

signal-to-noise ratio SNR(t) = ab_t / (1 - ab_t)
    t      linear      cosine
    1       1e+04    2.42e+04
   10         527    1.58e+03
  100        8.71        34.8
  250         1.1        5.54
  500      0.0853       0.976
  750     0.00336       0.169
 1000    4.04e-05    2.43e-09
linear: SNR falls below 1 at t = 260; 17% of steps have SNR < 0.001; ab_T = 4.04e-05
cosine: SNR falls below 1 at t = 497; 2% of steps have SNR < 0.001; ab_T = 2.43e-09

The iterated chain matches the closed form to within sampling error, about 0.003, so one draw of $\epsilon$ can replace $t$ steps during training. At $t = 1000$ only $0.0064$ of the disk’s amplitude remains under unit-variance noise. The SNR table shows what Nichol and Dhariwal objected to: the linear schedule crosses $\mathrm{SNR} = 1$ at step 260 and spends 17% of its steps below $\mathrm{SNR} = 0.001$, where in our reading almost nothing of a small image is left to denoise; the cosine schedule crosses it at step 497 and spends 2% of its steps below 0.001. The script runs in about 1 second.

An exact denoiser recovers the modes. For a toy distribution the optimal noise predictor $\epsilon^\ast$ can be written down, which isolates sampling from learning. The data are four $16 \times 16$ shape images (disk, square, plus sign, diagonal bar), each a mode of a Gaussian mixture, $x_0 = \mu_k + 0.1\,\xi$ with $\xi \sim \mathcal{N}(0, I)$, chosen with weights 0.4, 0.3, 0.2 and 0.1. Every noisy marginal $q(x_t)$ is then a mixture of $\mathcal{N}(\sqrt{\bar\alpha_t}\, \mu_k,\ (0.01\, \bar\alpha_t + 1 - \bar\alpha_t) I)$, whose score is exact. The code checks that $\hat x_0$ computed from $\epsilon^\ast$ equals the posterior mean computed from the mixture, then runs Algorithm 2 with $\sigma_t^2 = \beta_t$ for 1,000 samples under each of seeds 0 to 4 and reports how often each shape was the nearest and the root-mean-square (RMS) distance to it.

import numpy as np

T = 1000
beta = np.linspace(1e-4, 0.02, T)          # linear schedule of Ho et al. (2020)
alpha = 1 - beta
ab = np.cumprod(alpha)

# Toy data: four 16x16 shape images (disk, square, plus sign, diagonal bar) in [-1, 1].
# Each is a mode of a Gaussian mixture: x0 = mu_k + 0.1 * noise, chosen with weight w_k.
yy, xx = np.mgrid[:16, :16]
shapes = [(yy - 7.5) ** 2 + (xx - 7.5) ** 2 < 5.5 ** 2,          # disk
          (abs(yy - 7.5) < 5) & (abs(xx - 7.5) < 5),              # square
          (abs(yy - 7.5) < 2) | (abs(xx - 7.5) < 2),              # plus sign
          abs(yy - xx) < 2.5]                                     # diagonal bar
mu = np.array([np.where(m, 1.0, -1.0).ravel() for m in shapes])   # (4, 256)
w = np.array([0.4, 0.3, 0.2, 0.1])
s2 = 0.1 ** 2                                                    # within-mode variance

def posterior(x, t):
    """Responsibilities r_k(x_t) and the noisy-marginal variance v_t, for step t (1-based)."""
    a = ab[t - 1]
    v = a * s2 + 1 - a                     # p_t is a mixture of N(sqrt(a) mu_k, v I)
    d2 = ((x[:, None, :] - np.sqrt(a) * mu[None]) ** 2).sum(-1)
    logits = np.log(w) - d2 / (2 * v)
    r = np.exp(logits - logits.max(1, keepdims=True))
    return r / r.sum(1, keepdims=True), v

def score(x, t):                           # exact grad_x log p_t(x)
    r, v = posterior(x, t)
    return -(x - np.sqrt(ab[t - 1]) * (r @ mu)) / v

def eps_star(x, t):                        # optimal noise predictor
    return -np.sqrt(1 - ab[t - 1]) * score(x, t)

# The three parameterizations agree: x0 recovered from eps equals the posterior
# mean E[x0 | x_t] computed directly from the mixture.
rng = np.random.default_rng(0)
t = 400
x0 = mu[rng.choice(4, 500, p=w)] + 0.1 * rng.standard_normal((500, 256))
xt = np.sqrt(ab[t - 1]) * x0 + np.sqrt(1 - ab[t - 1]) * rng.standard_normal(x0.shape)
x0_from_eps = (xt - np.sqrt(1 - ab[t - 1]) * eps_star(xt, t)) / np.sqrt(ab[t - 1])
r, v = posterior(xt, t)
comp_mean = mu[None] + (s2 * np.sqrt(ab[t - 1]) / v) * (xt[:, None, :] - np.sqrt(ab[t - 1]) * mu[None])
x0_direct = (r[:, :, None] * comp_mean).sum(1)
print(f"t={t}: max |x0 from eps - E[x0|x_t]| = {np.abs(x0_from_eps - x0_direct).max():.1e}")

# Ancestral sampling (Ho et al., Algorithm 2) with sigma_t^2 = beta_t and the exact eps.
def sample(n, rng):
    x = rng.standard_normal((n, 256))
    for t in range(T, 0, -1):
        z = rng.standard_normal(x.shape) if t > 1 else 0.0
        x = (x - beta[t - 1] / np.sqrt(1 - ab[t - 1]) * eps_star(x, t)) / np.sqrt(alpha[t - 1]) \
            + np.sqrt(beta[t - 1]) * z
    return x

print("seed  disk  square  plus   bar   RMS to nearest mode")
for seed in range(5):
    x = sample(1000, np.random.default_rng(seed))
    d2 = ((x[:, None, :] - mu[None]) ** 2).mean(-1)
    k = d2.argmin(1)
    frac = np.bincount(k, minlength=4) / len(k)
    rms = np.sqrt(d2.min(1).mean())
    print(f"{seed:4d}" + "".join(f"{f:6.3f}" for f in frac) + f"   {rms:.4f}")
print("weights" + "".join(f"{f:6.3f}" for f in w) + "   0.1000 (data)")

Output:

t=400: max |x0 from eps - E[x0|x_t]| = 2.0e-15
seed  disk  square  plus   bar   RMS to nearest mode
   0 0.377 0.304 0.214 0.105   0.1007
   1 0.381 0.314 0.206 0.099   0.1008
   2 0.434 0.306 0.181 0.079   0.1008
   3 0.418 0.288 0.195 0.099   0.1006
   4 0.373 0.291 0.227 0.109   0.1007
weights 0.400 0.300 0.200 0.100   0.1000 (data)

The clean image recovered from the noise prediction agrees with the posterior mean to rounding error, so the noise, clean-image and score readings are one function. Started from pure noise, the sampler generates every shape in proportions that track the mixture weights: over the five seeds the disk’s share ranged from 0.373 to 0.434 against a weight of 0.4 and the bar’s from 0.079 to 0.109 against 0.1, deviations within about two standard errors of a share estimated from 1,000 draws (0.015 for the disk and 0.009 for the bar, in our calculation). The RMS distance to the nearest shape, 0.1006 to 0.1008, matches the data’s within-mode spread of 0.1: the samples are new draws around the templates, not copies of them. A trained network only approximates $\epsilon^\ast$; the gap between the two is what training reduces. The script runs in about 14 seconds.