Diffusion and large vision models, from forward process to few-step samplers

Why this subject matters now

Five years ago a practitioner could treat image generation as a niche. GANs produced sharp faces, VAEs produced blurry reconstructions, and neither mattered much outside graphics research. Diffusion models ended that. Between Ho et al.'s DDPM in 2020 and the latent-diffusion systems of 2022, likelihood-based generation caught and then passed GANs on fidelity while keeping the things GANs never had, stable training with a plain MSE loss, no discriminator collapse, mode coverage you can measure, and a controllable trade between compute and quality at inference time. Today the same mathematical object, in its flow-matching reformulation, is the generative backbone of essentially every frontier image and video system, among them Stable Diffusion 3 and Flux, Imagen and Veo at Google DeepMind, Sora at OpenAI, Movie Gen at Meta, Wan at Alibaba, Kling at Kuaishou. It also escaped media generation. Diffusion policies in robotics, protein structure generation, audio synthesis, and discrete-diffusion language models all reuse the machinery on this page.

What a practitioner is expected to know has shifted accordingly. In 2020 it was enough to recite the DDPM training loop. The current interview bar is structural. It asks why the epsilon objective is a reweighted variational bound and what the reweighting does, why DDIM can jump fifty steps at once when ancestral sampling cannot, why classifier-free guidance sharpens and where it breaks, why the field moved from U-Nets to transformers and from discrete diffusion to rectified flow, and how a 12-billion-parameter model produces a megapixel image in under a second on one GPU. Each of those answers is a derivation, not a fact, and the derivations share one spine. Everything is a way of estimating and then integrating the score of a noised data distribution. The page follows that spine from the forward process to the one-step distilled samplers.

Core theory

The forward process, a Gaussian Markov chain

Fix a data distribution \( q(x_0) \) on \( \R^d \) and a variance schedule \( \beta_1, \dots, \beta_T \in (0, 1) \). The forward process is the Markov chain that repeatedly shrinks the signal and adds Gaussian noise,

The shrink factor \( \sqrt{1-\beta_t} \) is chosen so that second moments are controlled. If \( \E[\|x_{t-1}\|^2] = d \) then \( \E[\|x_t\|^2] = (1-\beta_t) d + \beta_t d = d \). The chain neither inflates nor collapses the scale of its input, which is why this parameterization is called variance preserving. Writing \( \alpha_t = 1 - \beta_t \) and \( \bar\alpha_t = \prod_{s=1}^{t} \alpha_s \), the marginal of \( x_t \) given only \( x_0 \) has a closed form, and this single identity is what makes diffusion training cheap, since any timestep can be sampled directly without simulating the chain.

Claim. \( q(x_t \mid x_0) = \N\!\big(x_t;\ \sqrt{\bar\alpha_t}\, x_0,\ (1-\bar\alpha_t) I\big) \).

Proof by induction. The base case \( t = 1 \) is the definition, since \( \bar\alpha_1 = \alpha_1 \) and \( 1 - \bar\alpha_1 = \beta_1 \). For the inductive step, assume the claim at \( t - 1 \), so there is an \( \varepsilon \sim \N(0, I) \) with \( x_{t-1} = \sqrt{\bar\alpha_{t-1}}\, x_0 + \sqrt{1 - \bar\alpha_{t-1}}\, \varepsilon \). One more forward step with fresh noise \( \varepsilon' \sim \N(0, I) \), independent of \( \varepsilon \), gives

The mean coefficient is \( \sqrt{\alpha_t \bar\alpha_{t-1}} = \sqrt{\bar\alpha_t} \) by definition of the running product. The noise term is a sum of two independent zero-mean Gaussians, hence Gaussian with variance equal to the sum of variances,

So \( x_t = \sqrt{\bar\alpha_t}\, x_0 + \sqrt{1-\bar\alpha_t}\, \bar\varepsilon \) with \( \bar\varepsilon \sim \N(0, I) \), which is the claim. The collapse of the cross term is the only place the specific \( \sqrt{1-\beta_t} \) scaling is needed. An arbitrary shrink factor would still give a Gaussian marginal but would not telescope into a single product. If \( \beta_t \gt 0 \) for all \( t \) then \( \bar\alpha_T \to 0 \) as the schedule accumulates, and \( q(x_T \mid x_0) \to \N(0, I) \) regardless of \( x_0 \). The chain forgets its input, which is what lets sampling start from pure noise.

I verified this numerically rather than trusting the algebra. One million independent chains were simulated step by step from the fixed point \( x_0 = (1, -0.5) \) under the linear schedule \( \beta_t \) from \( 10^{-4} \) to \( 0.02 \) with \( T = 1000 \) (the DDPM defaults), on the H100 this site's benchmarks run on. At \( t = 500 \) the closed form predicts mean \( (0.28033, -0.14017) \) and per-coordinate standard deviation \( 0.9599 \) (from \( \bar\alpha_{500} = 0.078587 \)). The simulated chains gave mean \( (0.28008, -0.14186) \) and standard deviations \( (0.96033, 0.96024) \). At \( t = 1000 \), where \( \bar\alpha_T = 4.0 \times 10^{-5} \), the predicted mean is \( (0.00635, -0.00318) \), the simulated mean \( (0.00576, -0.00294) \), and both standard deviations agree with 1.0 to three decimals. The chain and the closed form are the same distribution to Monte Carlo precision.

Variance preserving versus variance exploding

The chain above is one of two standard conventions. The variance preserving (VP) family shrinks the signal while adding noise, so for unit-variance data the total variance stays near 1: \( \Var(x_t) = \bar\alpha_t \Var(x_0) + (1-\bar\alpha_t) I \). The variance exploding (VE) family, introduced with the noise conditional score networks of Song and Ermon (2019) and adopted by Karras et al.'s EDM (2022), never shrinks the signal,

In VE the terminal distribution is \( \N(x_0, \sigma_{\max}^2 I) \), which approximates a tractable prior only because \( \sigma_{\max} \) is chosen very large relative to the data scale. EDM uses \( \sigma_{\max} = 80 \) for images scaled to \( [-1, 1] \). The two are related by a change of variables. Dividing the VP state by \( \sqrt{\bar\alpha_t} \) gives \( x_t / \sqrt{\bar\alpha_t} = x_0 + \sqrt{(1-\bar\alpha_t)/\bar\alpha_t}\,\varepsilon \), a VE process with \( \sigma_t = \sqrt{(1-\bar\alpha_t)/\bar\alpha_t} \). Nothing statistical distinguishes them. What differs is numerics (VE keeps the clean image at its native scale, which simplifies preconditioning) and the induced discretization when a solver takes finite steps. A useful habit is to reduce every parameterization to the pair \( (\alpha_t, \sigma_t) \) with \( x_t = \alpha_t x_0 + \sigma_t \varepsilon \). VP has \( \alpha_t^2 + \sigma_t^2 = 1 \), VE has \( \alpha_t = 1 \), and the flow-matching linear path later on the page has \( \alpha_t + \sigma_t = 1 \). The signal-to-noise ratio \( \mathrm{SNR}(t) = \alpha_t^2 / \sigma_t^2 \) is invariant under all such rescalings, which is why Kingma et al.'s variational diffusion models (2021) showed the continuous-time bound depends on the schedule only through the SNR endpoints, a fact used repeatedly below.

The reverse process and the variational bound

Generation requires the reverse conditionals \( q(x_{t-1} \mid x_t) \), which are intractable because they depend on the unknown data distribution through Bayes' rule. The model therefore posits a learned reverse Markov chain

Choosing a Gaussian for each reverse step is not an arbitrary convenience. A classical result going back to the non-equilibrium-thermodynamics framing of Sohl-Dickstein et al. (2015), and originally to Feller's work on diffusion processes, is that the true reversal of a diffusion with sufficiently small steps is itself Gaussian to first order in the step size. For small \( \beta_t \), \( q(x_{t-1} \mid x_t) \) is approximately normal with a mean shifted from \( x_t \) by a term involving the score of the marginal. Small steps are what purchase the Gaussian approximation. Large learned jumps require the heavier machinery of distillation covered later.

The training objective is a variational bound on the negative log-likelihood. Since \( p_\theta(x_0) = \int p_\theta(x_{0:T})\, dx_{1:T} \), multiply and divide by the forward chain and apply Jensen's inequality to the concave logarithm,

Expanding both chains,

The forward factors point the wrong way. They condition on the past while the model conditions on the future. The fix is the observation that conditioning additionally on \( x_0 \) costs nothing (the chain is Markov, so \( q(x_t \mid x_{t-1}) = q(x_t \mid x_{t-1}, x_0) \)) and lets Bayes' rule flip the direction, for \( t \ge 2 \),

Substituting this for every \( t \ge 2 \), the marginal ratios \( q(x_t \mid x_0) / q(x_{t-1} \mid x_0) \) telescope. Every numerator cancels the next denominator, leaving only \( q(x_T \mid x_0) / q(x_1 \mid x_0) \), and the surviving \( \log q(x_1 \mid x_0) \) from the \( t = 1 \) term cancels the denominator. Grouping what remains by which model factor it faces,

\( \L_T \) has no parameters. It measures how close the forward chain gets to the prior, and with \( \bar\alpha_T = 4\times 10^{-5} \) it is a fraction of a nat. \( \L_0 \) is a one-step decoder term, handled in DDPM by a discretized Gaussian likelihood over the 256 pixel levels. The heart of the objective is the \( T - 1 \) middle terms, each a KL divergence between two Gaussians. One is the model's reverse step, the other a distribution \( q(x_{t-1} \mid x_t, x_0) \) that is fully tractable because it conditions on the clean image. The entire training signal is "match the reverse step you would take if you already knew the answer, averaged over what the answer might be."

The tractable posterior, by completing the square

Everything reduces to computing \( q(x_{t-1} \mid x_t, x_0) \). By Bayes' rule it is proportional (in \( x_{t-1} \)) to \( q(x_t \mid x_{t-1})\, q(x_{t-1} \mid x_0) \), a product of two Gaussians whose exponents are both quadratics in \( x_{t-1} \), so the product is Gaussian and its parameters fall out of completing the square. Write only the exponent, dropping every factor that does not involve \( x_{t-1} \),

where the squares are applied coordinatewise. Each coordinate separates because all covariances are isotropic. Expand and collect powers of \( x_{t-1} \). The quadratic coefficient is

using the same identity \( \alpha_t(1-\bar\alpha_{t-1}) + \beta_t = 1 - \bar\alpha_t \) that closed the induction. The linear coefficient (of \( -2 x_{t-1} \cdot \) inside the square) is

A quadratic \( -\frac{1}{2}(A x^2 - 2 B x) \) is the exponent of \( \N(B/A,\ 1/A) \), so the posterior is Gaussian, \( q(x_{t-1} \mid x_t, x_0) = \N\big(\tilde\mu_t(x_t, x_0),\ \tilde\beta_t I\big) \), with

Both facts deserve a reading. The posterior variance \( \tilde\beta_t \) is \( \beta_t \) shrunk by the factor \( (1-\bar\alpha_{t-1})/(1-\bar\alpha_t) \lt 1 \). Knowing the clean image removes some but not all uncertainty about the intermediate state, and near \( t = T \) where \( \bar\alpha \approx 0 \) the shrink factor approaches 1, since at high noise even the clean image barely constrains one step. The posterior mean is a convex-looking blend of the current noisy state and the clean image, with weights that shift from almost entirely \( x_t \) at small \( t \) toward a heavier \( x_0 \) contribution as \( \beta_t \) grows.

Training needs the mean in terms of what the network can predict. Solving the closed-form marginal for \( x_0 = (x_t - \sqrt{1-\bar\alpha_t}\,\varepsilon)/\sqrt{\bar\alpha_t} \) and substituting into \( \tilde\mu_t \), the \( x_0 \) coefficient contributes \( \frac{\sqrt{\bar\alpha_{t-1}}\,\beta_t}{(1-\bar\alpha_t)\sqrt{\bar\alpha_t}} = \frac{\beta_t}{(1-\bar\alpha_t)\sqrt{\alpha_t}} \) on \( x_t \) (using \( \sqrt{\bar\alpha_{t-1}/\bar\alpha_t} = 1/\sqrt{\alpha_t} \)), so the total \( x_t \) coefficient is

once more by the induction identity, while the \( \varepsilon \) term carries coefficient \( -\frac{\beta_t \sqrt{1-\bar\alpha_t}}{(1-\bar\alpha_t)\sqrt{\alpha_t}} = -\frac{\beta_t}{\sqrt{\alpha_t}\sqrt{1-\bar\alpha_t}} \). The posterior mean collapses to

This is the central formula of the subject. The ideal reverse step rescales the current state and subtracts a fixed multiple of the noise that was mixed in. Parameterize the model the same way, \( \mu_\theta(x_t, t) = \frac{1}{\sqrt{\alpha_t}}\big(x_t - \frac{\beta_t}{\sqrt{1-\bar\alpha_t}}\, \varepsilon_\theta(x_t, t)\big) \), so that the network's only job is predicting the noise.

From KL terms to the epsilon objective, and what gets dropped

The KL divergence between two isotropic Gaussians with the same variance \( \sigma_t^2 \) is \( \|\mu_1 - \mu_2\|^2 / (2\sigma_t^2) \) (the log-determinant and trace terms cancel exactly when covariances match). With both means in the parameterization above, their difference is \( \frac{\beta_t}{\sqrt{\alpha_t}\sqrt{1-\bar\alpha_t}} (\varepsilon - \varepsilon_\theta) \), so each middle term of the bound is

Ho et al.'s empirical finding, the one that made diffusion practical, is that discarding the time-dependent weight and minimizing

produces better samples than the exact bound. It is worth being precise about which weighting was dropped. With the common choice \( \sigma_t^2 = \tilde\beta_t \), the weight simplifies to \( \frac{\beta_t}{2 \alpha_t (1-\bar\alpha_{t-1})} \), which is large at small \( t \) (as \( t \to 1 \), \( 1 - \bar\alpha_{t-1} \to \beta \)-scale and the ratio blows up) and small at large \( t \). The exact bound therefore spends most of its gradient on the final, nearly-noiseless steps, which dominate log-likelihood but contribute almost nothing to perceptual quality. \( \L_{\text{simple}} \) reweights toward mid and high noise levels, where the model learns global structure. In SNR language, using the identity \( \|\varepsilon - \varepsilon_\theta\|^2 = \mathrm{SNR}(t)\, \|x_0 - \hat{x}_0\|^2 \) derived in the parameterization section, uniform epsilon weighting is an SNR-weighted reconstruction loss. The choice of weighting is not cosmetic bookkeeping. It is the main lever the later literature (VDM, EDM, min-SNR, and the flow-matching timestep distributions) keeps adjusting, and every one of those papers can be read as a different answer to "which noise levels deserve gradient."

One more piece completes classic DDPM, the reverse variance. Ho et al. fixed \( \sigma_t^2 \) to either \( \beta_t \) or \( \tilde\beta_t \) (the two extremes for a data distribution that is respectively unit-Gaussian or deterministic) and reported similar results. Nichol and Dhariwal's improved DDPM (2021) learned an interpolation \( \sigma_t^2 = \exp\big(v_\theta \log\beta_t + (1-v_\theta)\log\tilde\beta_t\big) \) trained by the exact bound on the variance only, which cut the step count needed for good likelihood and mattered again once fewer than 50 steps were taken.

The score-based view, denoising score matching

A parallel line of work, Song and Ermon's noise conditional score networks (2019) building on Hyvärinen's score matching (2005) and Vincent's denoising reformulation (2011), arrives at the same algorithm from a different starting point. The score of a density is \( \nabla_x \log p(x) \), a vector field pointing toward higher probability. If the score of the noised marginals \( q_t(x_t) = \int q(x_t \mid x_0)\, q(x_0)\, dx_0 \) were known, samples could be drawn by Langevin dynamics,

whose stationary distribution is \( q_t \) as \( \eta \to 0 \). Learning the score naively by regressing \( \| s_\theta(x_t) - \nabla \log q_t(x_t) \|^2 \) is impossible because \( q_t \) is exactly the unknown. Denoising score matching resolves this with an identity. Regressing onto the score of the conditional \( q(x_t \mid x_0) \), which is a known Gaussian, has the same minimizer.

Derivation. Expand the intractable objective and keep only the terms that depend on \( \theta \),

The cross term is where the trick lives,

using \( \nabla q_t = q_t \nabla \log q_t \) in both directions and exchanging gradient and integral. Adding back \( \E \| s_\theta \|^2 \) and completing the square in the joint expectation shows

with \( C \) independent of \( \theta \), so the two objectives have identical gradients. The conditional score is elementary. With \( q(x_t \mid x_0) = \N(\sqrt{\bar\alpha_t}\, x_0, (1-\bar\alpha_t) I) \) and \( x_t = \sqrt{\bar\alpha_t}\, x_0 + \sigma_t \varepsilon \) where \( \sigma_t = \sqrt{1-\bar\alpha_t} \),

So the DSM regression target is \( -\varepsilon/\sigma_t \), and a network trained to predict \( \varepsilon \) already is a score model under the identification \( s_\theta(x_t, t) = -\varepsilon_\theta(x_t, t)/\sigma_t \). The two objectives differ only by the constant factor \( \sigma_t^{-2} \) per noise level, which is a choice of loss weighting, not of minimizer at any single \( t \). The DDPM training loop and the NCSN training loop are the same computation written in different units. The minimizer itself has a name. Since the regression target given \( x_t \) is \( \E[-\varepsilon/\sigma_t \mid x_t] \), the optimal \( \varepsilon^*(x_t, t) = \E[\varepsilon \mid x_t] \), and by Tweedie's formula (worked problem 2) this is equivalent to computing the posterior mean \( \E[x_0 \mid x_t] \). Everything a diffusion model learns is one conditional expectation.

This identity is also checkable, and I checked it. An epsilon model was trained on data drawn from \( \N((0.7, -0.3),\ 0.36\, I) \), for which the noised marginal \( q_t \) is Gaussian with known parameters, so the true score is available in closed form. Comparing \( -\varepsilon_\theta/\sigma_t \) against the analytic score on 50,000 fresh samples per noise level gave median relative error 2.1% at \( t = 50 \), 0.9% at \( t = 250 \), 0.5% at \( t = 500 \), 0.25% at \( t = 900 \), with mean cosine similarity 0.9999 or better at every level. A plain MSE noise predictor really does converge to the score field, most accurately at high noise where the marginal is smoothest.

The SDE formulation and the probability-flow ODE

Song et al. (2021) unified both discrete constructions as discretizations of stochastic differential equations. The forward corruption becomes an Itô SDE \( dx = f(x, t)\, dt + g(t)\, dw \), and the VP chain corresponds to

The correspondence is exact in the small-step limit. Discretize with step \( \Delta \) and per-step variance \( \beta_t = \beta(t)\Delta \), giving \( x_{t+\Delta} = (1 - \tfrac{1}{2}\beta_t)\, x_t + \sqrt{\beta_t}\, z \), and compare with the DDPM step \( \sqrt{1-\beta_t}\, x_t + \sqrt{\beta_t}\, z \). Since \( \sqrt{1-\beta_t} = 1 - \tfrac{1}{2}\beta_t + O(\beta_t^2) \), the two agree to first order. The closed-form marginal also survives the limit, as \( \bar\alpha_t = \prod_s (1 - \beta_s) = \exp\big(\sum_s \log(1-\beta_s)\big) \to \exp\big(-\int_0^t \beta(s)\, ds\big) \), which is exactly the variance decay of the linear SDE.

Two reverse-time processes share the forward marginals, and both can be derived in one computation from the Fokker-Planck equation, which any density \( p_t \) evolved by the forward SDE satisfies,

The diffusion term can be rewritten as a transport term, because \( \tfrac{1}{2} g^2 \Delta p_t = \tfrac{1}{2} g^2\, \nabla \cdot (\nabla p_t) = \nabla \cdot \big( \tfrac{1}{2} g^2\, p_t\, \nabla \log p_t \big) \), using \( \nabla p_t = p_t \nabla \log p_t \) once more. Substituting,

which is a continuity equation. The deterministic probability-flow ODE

transports \( p_0 \) into exactly the marginals \( p_t \) of the stochastic forward process, with no randomness anywhere. More generally, for any \( \lambda \ge 0 \), splitting the diffusion term as \( \tfrac12 g^2 = \tfrac12(1+\lambda^2) g^2 - \tfrac12 \lambda^2 g^2 \) and converting only the first part into transport shows that

run backward in time has the same marginals for every \( \lambda \). The extra noise injected by the \( \lambda g\, d\bar w \) term is exactly compensated by the stronger score drift. \( \lambda = 0 \) is the probability-flow ODE, and \( \lambda = 1 \) is Anderson's reverse-time SDE (1982),

which is what ancestral DDPM sampling discretizes. Substituting the VP drift and the learned score \( \nabla \log p_t \approx -\varepsilon_\theta / \sigma_t \) into the reverse SDE and taking an Euler-Maruyama step of size \( \beta_t \) reproduces, to first order, the DDPM update \( x_{t-1} = \frac{1}{\sqrt{\alpha_t}}\big( x_t - \frac{\beta_t}{\sqrt{1-\bar\alpha_t}} \varepsilon_\theta \big) + \sqrt{\beta_t} z \). Expand \( 1/\sqrt{\alpha_t} = 1 + \tfrac12 \beta_t + O(\beta_t^2) \) and the drift terms match term by term. This is the rigorous content of "DDPM, NCSN, and score SDEs are the same model". The variational chain is a first-order discretization of the reverse VP-SDE, its training loss is denoising score matching in disguise, and the family of samplers, from fully stochastic ancestral steps to the deterministic ODE, is a one-parameter family sharing the same learned vector field and the same marginals. Which member you integrate, and how accurately, is the entire subject of fast sampling below.

The stochastic-deterministic choice is not a tie in practice. The ODE gives a bijection between noise and images (the basis of inversion and editing) and admits high-order solvers. The SDE's injected noise contracts accumulated error toward the correct marginals, which is why EDM-style stochastic churn helps when the model is imperfect and many steps are affordable, and pure ODE solving wins when steps are scarce.

Noise schedules and the log-SNR axis

A schedule is a choice of \( (\alpha_t, \sigma_t) \), best summarized by the signal-to-noise ratio \( \mathrm{SNR}(t) = \alpha_t^2/\sigma_t^2 \) or its logarithm \( \lambda_t = \log \mathrm{SNR}(t) \). The DDPM linear schedule (\( \beta_t \) linear from \( 10^{-4} \) to \( 0.02 \), \( T = 1000 \)) spans \( \lambda \) from about \( +9.2 \) at \( t = 1 \) (where \( \bar\alpha \approx 0.9999 \)) down to \( -10.1 \) at \( t = T \) (using the measured \( \bar\alpha_T = 4\times10^{-5} \)). Its defect is visible in the measured curve. \( \bar\alpha \) falls to 0.524 by \( t = 250 \) and 0.079 by \( t = 500 \), so the second half of the chain is spent at nearly-destroyed signal where there is little left to learn. Nichol and Dhariwal's cosine schedule defines the cumulative product directly,

which makes \( \lambda_t \) fall roughly linearly in \( t \) through the midrange. Equal spacing in time is close to equal spacing in log-SNR, so training effort and sampler steps spread evenly across perceptually distinct noise levels. The small offset \( s \) keeps \( \beta_1 \) from being vanishingly small, and implementations clip \( \beta_t \le 0.999 \) at the noisy end. Kingma et al.'s VDM result makes the log-SNR axis canonical. Reparameterizing time leaves the continuous variational bound unchanged except through the endpoints \( \lambda_{\min}, \lambda_{\max} \), so "schedule" during training really means "importance distribution over \( \lambda \)", and schedules should be compared as densities on the log-SNR axis, not as curves in \( t \).

Schedule choice interacts with resolution, and this bit of folk wisdom has a clean argument, made explicit in Hoogeboom et al.'s simple diffusion (2023) and adopted as the shift parameter in Stable Diffusion 3 and Flux. Consider downsampling a noised image by averaging \( k \times k \) pixel blocks. Averaging preserves the signal (natural images are locally correlated, so the block mean is close to each pixel's value) but averages \( k^2 \) independent noise draws, shrinking the noise standard deviation by \( k \). The downsampled image therefore has SNR larger by \( k^2 \). A per-pixel noise level that thoroughly destroys a \( 64\times64 \) image still leaves the global structure of a \( 1024\times1024 \) image legible in its low frequencies. A schedule tuned at low resolution, transplanted to high resolution, spends almost no time in the regime where coarse structure is actually uncertain, and models trained that way produce incoherent layouts. The fix is to shift the schedule so the same information destruction happens at the same \( t \), moving log-SNR down by \( 2 \log(k) \) when resolution grows by \( k \) per side. SD3 expresses the same correction as a timestep shift \( t' = \frac{c\, t}{1 + (c - 1)\, t} \) with \( c \approx 3 \) for \( 1024^2 \) training, and video models shift further still because temporal redundancy compounds the spatial kind.

Parameterizations, from epsilon to x-zero to v

At any fixed \( (x_t, t) \), the quantities \( x_0 \), \( \varepsilon \), and the score are affinely related through \( x_t = \alpha_t x_0 + \sigma_t \varepsilon \), so a network may predict any of them. The choice changes conditioning of the regression and the implicit loss weighting, not the expressible model class. Predicting \( \varepsilon \) is well scaled everywhere (the target is always unit variance) but degrades as \( t \to T \), since reconstructing \( \hat{x}_0 = (x_t - \sigma_t \varepsilon_\theta)/\alpha_t \) divides by \( \alpha_t \to 0 \), so tiny epsilon errors become huge image errors precisely where samplers take their first, largest steps. Predicting \( x_0 \) has the mirrored failure at \( t \to 0 \), where the target becomes trivially recoverable and the loss stops teaching denoising.

Salimans and Ho's v-prediction (2022) removes both failure modes. For a VP schedule, \( \alpha_t^2 + \sigma_t^2 = 1 \) invites the angle substitution \( \alpha_t = \cos\phi_t,\ \sigma_t = \sin\phi_t \), so \( x_t = \cos\phi\, x_0 + \sin\phi\, \varepsilon \) is a rotation in the plane spanned by the signal and its noise. Define the velocity as the derivative along that rotation,

\( (x_t, v) \) is an orthogonal rotation of \( (x_0, \varepsilon) \), so both directions invert cleanly,

At the clean end (\( \sigma \to 0 \)), \( v \to \varepsilon \), so the target is the noise, which is the informative quantity there. At the noisy end (\( \alpha \to 0 \)), \( v \to -x_0 \), so the target is the image, again the informative quantity. Neither reconstruction ever divides by a vanishing coefficient, the target has unit variance at every \( t \), and a single-step prediction from pure noise is meaningful, which is why progressive distillation requires v-parameterization. The loss weightings also line up neatly. Fixing \( x_t \), errors in the three targets are proportional (\( \varepsilon \)-error \( = \alpha_t \cdot v \)-error, \( x_0 \)-error \( = \sigma_t \cdot v \)-error), so measured in \( x_0 \)-MSE units, the epsilon loss carries weight \( \mathrm{SNR}(t) \), the v loss carries \( \mathrm{SNR}(t) + 1 \) (since \( 1/\sigma_t^2 = (\alpha_t^2+\sigma_t^2)/\sigma_t^2 \)), and the plain \( x_0 \) loss carries weight 1. Hang et al.'s min-SNR-\( \gamma \) weighting (2023) truncates the first, using weight \( \min(\mathrm{SNR}(t), \gamma) \) in \( x_0 \) units with \( \gamma = 5 \), so the near-clean steps (SNR in the thousands) stop dominating the epsilon objective. They report roughly 3× faster convergence on ImageNet DiT training. EDM reaches a similar endpoint differently, choosing its \( c_{\text{skip}}, c_{\text{out}}, c_{\text{in}} \) preconditioning so the effective target has unit variance and the loss weight is flat where the data manifold is actually ambiguous.

Sampling I, ancestral DDPM and the DDIM family

Ancestral sampling runs the learned chain. From \( x_T \sim \N(0, I) \), repeat

with \( \sigma_t^2 = \tilde\beta_t \) and no noise at \( t = 1 \). It is faithful to the bound and costs \( T \) network evaluations, one thousand in the original recipe. The route to fewer steps begins with Song et al.'s DDIM (2021) observation that the training objective never used the forward chain's Markov structure. \( \L_{\text{simple}} \) depends only on the marginals \( q(x_t \mid x_0) \). Any family of joint distributions with those marginals trains the identical network, so one may pick the family whose reverse process is cheapest to sample. DDIM chooses, for a free parameter \( \sigma \ge 0 \),

The marginal condition checks directly. If \( x_t \mid x_0 \sim \N(\sqrt{\bar\alpha_t} x_0, (1-\bar\alpha_t) I) \), the normalized residual \( (x_t - \sqrt{\bar\alpha_t} x_0)/\sqrt{1-\bar\alpha_t} \) is a standard Gaussian, so \( x_{t-1} \mid x_0 \) has mean \( \sqrt{\bar\alpha_{t-1}}\, x_0 \) and variance \( (1 - \bar\alpha_{t-1} - \sigma^2) + \sigma^2 = 1 - \bar\alpha_{t-1} \), which is exactly \( q(x_{t-1} \mid x_0) \), for every \( \sigma \). This is a non-Markovian forward family (the implied \( q_\sigma(x_t \mid x_{t-1}, x_0) \) depends on \( x_0 \)) indexed by how much fresh randomness each reverse step injects. The sampler replaces \( x_0 \) with the model's estimate \( \hat{x}_0 = (x_t - \sqrt{1-\bar\alpha_t}\, \varepsilon_\theta)/\sqrt{\bar\alpha_t} \),

The conventional knob is \( \sigma = \eta \sqrt{(1-\bar\alpha_{t-1})/(1-\bar\alpha_t)} \sqrt{1 - \bar\alpha_t/\bar\alpha_{t-1}} \). Setting \( \eta = 1 \) recovers ancestral DDPM (the expression equals \( \sqrt{\tilde\beta_t} \)), and \( \eta = 0 \) is the deterministic limit, where the update is "estimate the clean image, then re-noise it analytically to the previous level using the same predicted noise." Because nothing random is added, consecutive steps can be composed across large gaps. The update is well defined between any pair \( \bar\alpha_t \to \bar\alpha_{t'} \), so a 1000-step training chain can be sampled on any sub-sequence of 50 or 20 timesteps. The same determinism gives DDIM inversion, the encoder used by image-editing pipelines. Running the update with time reversed finds the latent that reproduces a given image.

Sampling II: the ODE view and higher-order solvers

Why does skipping steps degrade gracefully instead of catastrophically? Because deterministic DDIM is a first-order integrator of the probability-flow ODE, and its error is a discretization error, controlled by step size and curvature rather than by any per-step validity condition. The sharp way to see it, due to Lu et al.'s DPM-Solver (2022), is that the PF-ODE is semi-linear. In \( (\alpha, \sigma) \) notation with \( \lambda = \log(\alpha/\sigma) \) (half the log-SNR), the VP probability-flow ODE \( \dot x = -\tfrac12 \beta x + \tfrac{\beta}{2\sigma} \varepsilon_\theta \) has an exact variation-of-constants solution

in which the linear drift has been integrated exactly and all remaining error lives in the integral of the network output against \( e^{-\lambda} \). Approximating \( \hat\varepsilon_\theta \) as constant over the interval gives \( x_t = \frac{\alpha_t}{\alpha_s} x_s - \alpha_t \big( \frac{\sigma_s}{\alpha_s} - \frac{\sigma_t}{\alpha_t} \big) \varepsilon_\theta = \frac{\alpha_t}{\alpha_s} x_s + \big( \sigma_t - \frac{\alpha_t \sigma_s}{\alpha_s} \big) \varepsilon_\theta \), which is algebraically identical to the deterministic DDIM update. DDIM is the first-order exponential integrator. Taylor-expanding \( \hat\varepsilon_\theta \) in \( \lambda \) to first or second order, with derivatives estimated from previous evaluations (multistep) or midpoints (single-step), gives the second- and third-order DPM-Solver rules whose local error is \( O(h^3) \) or \( O(h^4) \) in the step size \( h = \Delta\lambda \). The practical consequence is large. Error that shrinks like \( h^2 \) or \( h^3 \) instead of \( h \) means 10 to 20 well-placed steps do the work of hundreds of Euler steps.

The refinements that ship in production samplers are small but consequential. DPM-Solver++ rewrites the expansion around the \( \hat{x}_0 \) prediction instead of \( \varepsilon \), which keeps iterates near the data range and composes correctly with the thresholding used under strong guidance. Its 2M (second order, multistep) variant is the workhorse default in Stable Diffusion tooling. Zhao et al.'s UniPC (2023) adds a corrector step of matching order that reuses the model call, buying roughly one order of accuracy at 5 to 10 steps. Karras et al.'s EDM analysis showed the discretization grid matters as much as the rule. Their \( \sigma \)-spacing (\( \sigma_i^{1/\rho} \) linear, \( \rho = 7 \)) concentrates steps where the ODE trajectory curves, and their Heun predictor-corrector (a second-order trapezoidal rule) reached then-state-of-the-art ImageNet FID at 35 network evaluations and matched hundreds-of-steps baselines. Zhang and Chen's DEIS at Georgia Tech reached similar conclusions with exponential integrators independently of the Tsinghua line.

My measured sweep on the two-moons model makes the tradeoff concrete, using sliced 2-Wasserstein distance to 20,000 held-out points (metric noise floor 0.0041, measured by comparing two independent data draws). Ancestral sampling with all 1000 steps scores 0.0151. Deterministic DDIM scores 0.0136 at 1000 steps, 0.0155 at 200, 0.0200 at 100, 0.0254 at 50, 0.0569 at 20, 0.1345 at 10, 0.3344 at 5, 1.3168 at 2. The \( \eta = 1 \) column is consistently worse below 200 steps (0.0809 versus 0.0569 at 20 steps), measuring exactly the "stochasticity hurts when steps are scarce" effect predicted by the ODE analysis. The knee between 20 and 50 steps, with degradation graceful above it and steep below, is the same shape production image models exhibit, and it is why the few-step methods at the end of this page exist. Below about 10 steps, better integration cannot save a curved trajectory, and the model itself must change.

Guidance: classifier and classifier-free

Conditional generation wants samples from \( p(x \mid c) \) for a label or prompt \( c \). At every noise level, Bayes' rule splits the conditional score,

since \( \log p_t(x_t \mid c) = \log p_t(x_t) + \log p_t(c \mid x_t) - \log p(c) \) and the last term has no \( x_t \) dependence. Dhariwal and Nichol's classifier guidance (2021) implements the second term with a separate classifier trained on noised images and, crucially, scales it by \( s \gt 1 \), giving \( \tilde\varepsilon = \varepsilon_\theta - s\, \sigma_t \nabla_{x_t} \log p_\phi(c \mid x_t) \), using the score-epsilon dictionary from above. Scaling by \( s \) means sampling from a distribution proportional to \( p_t(x_t)\, p_t(c \mid x_t)^s \). The classifier term is exponentiated, concentrating mass where the class is unambiguous. This traded diversity for fidelity well enough to beat GANs on ImageNet FID, at the cost of training and backpropagating through a noise-robust classifier at every sampling step.

Ho and Salimans' classifier-free guidance (2022) removes the classifier by estimating the same quantity from the generative model itself. Train one network with the conditioning randomly dropped (10 to 20 percent of examples see a null token \( \varnothing \)), so it learns both \( \varepsilon_\theta(x_t, c) \) and \( \varepsilon_\theta(x_t, \varnothing) \). Since \( \nabla \log p_t(c \mid x_t) = \nabla \log p_t(x_t \mid c) - \nabla \log p_t(x_t) \), the implicit classifier's gradient is the difference of the two heads, and guided sampling uses

In score form this is \( (1+w) \nabla \log p_t(x_t \mid c) - w \nabla \log p_t(x_t) = \nabla \log \big[ p_t(x_t \mid c)^{1+w} / p_t(x_t)^{w} \big] \). At each noise level the sampler follows the score of an over-sharpened conditional, the conditional density tilted by an extra \( w \) powers of the implicit classifier. Two caveats make this an approximation rather than a theorem about the final distribution. The tilted densities at different \( t \) are not the diffused versions of any single tilted data distribution (powering and convolving with a Gaussian do not commute), so CFG sampling does not exactly target \( p(x \mid c)^{1+w} p(x)^{-w} \) at \( t = 0 \), and the network pair is only an estimate of the two scores. In practice the qualitative account is accurate. Guidance scales of 5 to 15 (in the common convention where the scale is \( 1 + w \)) markedly improve prompt adherence and per-image coherence while visibly shrinking diversity and pushing colors and contrast toward saturation, because the mean of the tilted distribution overshoots the conditional mean and its variance contracts. Worked problem 4 computes both effects exactly for Gaussians, and my measured class-conditional sweep shows the same signature. At guidance weight \( w = 3 \) the per-axis sample standard deviation contracts from \( (0.81, 0.62) \) to \( (0.45, 0.46) \) against a data value of \( (0.82, 0.63) \), and at \( w = 7 \) the sampler destabilizes outright, flinging mass off the manifold (std inflates to \( (0.72, 1.87) \) while distributional distance grows 8×).

The known fixes target specific mechanisms. Saturation comes from iterates leaving the trained data range, so Imagen's dynamic thresholding rescales \( \hat{x}_0 \) per sample so a chosen percentile (99.5) fits in \( [-1, 1] \). It requires an \( x_0 \)-space sampler, one reason DPM-Solver++ predicts \( x_0 \). Lin et al.'s guidance rescaling matches the standard deviation of the guided prediction back to the conditional one. Kynkäänniemi et al. (2024) showed most of guidance's benefit comes from a middle interval of noise levels, and that disabling it at the high-noise end restores much of the lost diversity at equal fidelity. Guidance-interval scheduling is now common in production samplers. Karras et al.'s autoguidance (2024) replaces the unconditional branch with a smaller or undertrained version of the same model, decoupling "sharpen toward better samples" from "sharpen toward the conditional", and won further FID at high guidance. Because every CFG step costs two forward passes, distilled models (below) usually bake guidance into the student, either by distilling a guided teacher or by conditioning the student on \( w \) directly.

Flow matching and rectified flow

Flow matching, introduced by Lipman et al. (2023) with parallel formulations by Liu et al. (rectified flow, 2022) and Albergo and Vanden-Eijnden (stochastic interpolants, 2023), reframes generation as learning the velocity field of a chosen deterministic transport. Fix a probability path, a family \( p_t \) interpolating from a tractable \( p_0 = \N(0, I) \) at \( t = 0 \) to the data distribution at \( t = 1 \), generated by an unknown marginal velocity field \( u_t \) through the continuity equation \( \partial_t p_t = -\nabla \cdot (u_t p_t) \). Regressing on \( u_t \) directly is intractable for the same reason the score was. The conditional flow matching objective conditions on the data endpoint. Choose per-example Gaussian paths \( p_t(x \mid x_1) = \N(\mu_t(x_1), \sigma_t^2 I) \) whose mixture over \( x_1 \) is \( p_t \), with conditional velocity \( u_t(x \mid x_1) \) known in closed form, and minimize

The gradient equivalence. The marginal field is the posterior average of conditional fields. Dividing the mixture continuity equation by \( p_t \) shows \( u_t(x) = \E_{x_1 \mid x}\big[ u_t(x \mid x_1) \big] \) with weights \( p_t(x \mid x_1) q(x_1) / p_t(x) \). Expanding both objectives, the \( \|v_\theta\|^2 \) terms agree because both expectations are over the same marginal \( p_t \), and the cross terms agree by the tower property,

The remaining terms do not involve \( \theta \), so \( \nabla_\theta \L_{\text{CFM}} = \nabla_\theta \L_{\text{FM}} \). The tractable per-example regression finds the true marginal field. This is the same maneuver as denoising score matching, with "condition on \( x_0 \) to make the target computable" replaced by "condition on \( x_1 \)". The minimizer is again a conditional expectation, \( v^*(x, t) = \E[u_t(x \mid x_1) \mid x] \).

For a Gaussian path \( x = \alpha_t x_1 + \sigma_t \varepsilon \), differentiating the sample map gives the conditional velocity \( u_t(x \mid x_1) = \frac{\dot\sigma_t}{\sigma_t}(x - \alpha_t x_1) + \dot\alpha_t x_1 \). The linear (optimal-transport) path takes \( \alpha_t = t,\ \sigma_t = 1 - t \), i.e. \( x_t = (1-t)\, x_0 + t\, x_1 \) with \( x_0 \) the noise draw, and the velocity collapses to a constant along each segment,

The training loop becomes simple. Draw a noise-data pair, pick a uniform \( t \), interpolate linearly, regress the network on the difference vector. Sampling integrates \( \dot x = v_\theta(x, t) \) from noise at \( t = 0 \) to data at \( t = 1 \) with any ODE solver.

Diffusion is a point in this family, not a rival. Choose instead the VP path \( \alpha_t = \sqrt{\bar\alpha_t},\ \sigma_t = \sqrt{1-\bar\alpha_t} \). Substituting Tweedie's identity \( \E[x_1 \mid x] = (x + \sigma_t^2 \nabla \log p_t(x))/\alpha_t \) into the posterior-averaged conditional velocity gives

and for the VP schedule, where \( \dot\alpha/\alpha = -\tfrac12\beta \) and \( \dot\sigma\sigma = -\alpha\dot\alpha = \tfrac12 \beta \alpha^2 \), the score coefficient is \( -\tfrac12\beta\sigma^2 - \tfrac12\beta\alpha^2 = -\tfrac12 \beta \), giving \( u_t(x) = -\tfrac12 \beta x - \tfrac12 \beta \nabla \log p_t(x) \), exactly the probability-flow ODE drift. Flow matching with the VP Gaussian path is diffusion's deterministic sampler, learned directly as a velocity instead of assembled from an epsilon model. The linear path is a different, flatter route through the same space of probability paths, one whose velocity target is better conditioned at both endpoints (it is \( v \)-prediction up to a rotation and sign) and whose trajectories have lower curvature, so fixed-step ODE solvers lose less. This, plus the cleaner \( t \in [0,1] \) formulation, is why Stable Diffusion 3, Flux, Meta's Movie Gen, Alibaba's Wan, and most current video systems train rectified flow objectives. SD3's ablation across 61 formulation/schedule variants found rectified flow with logit-normal timestep sampling the most reliable configuration.

Why rectified flow straightens. Liu et al.'s deeper observation concerns the coupling. Training pairs each noise draw with a random data point, so conditional segments cross, and the marginal field, their posterior average, is curved even though every segment is straight. But the trained flow itself defines a better coupling. Transport each \( x_0 \) along the learned ODE to its endpoint \( \hat{x}_1 \), and retrain on the coupled pairs \( (x_0, \hat{x}_1) \). This "reflow" operation provably does not increase any convex transport cost and strictly reduces the crossing that causes curvature. Iterated, the coupling approaches one whose paths are non-crossing and nearly straight, and a straight path is integrated exactly by a single Euler step. My measurements make this concrete. Straightness, measured as \( \E\|v_\theta(x_t, t) - (\hat x_1 - x_0)\|^2 \) along 100-step trajectories, fell from 0.914 to 0.0004 after one reflow round, a drop of three orders of magnitude, effectively straight. In sample quality (sliced Wasserstein, floor 0.0041), the base flow degraded from 0.0123 at 100 Euler steps to 0.0990 at 8 steps and 0.9816 at 1 step, useless. The reflowed model scored 0.0158 at 8 steps, 0.0192 at 2, and 0.0191 at one step, matching the base model's hundred-step quality with a single network evaluation. The cost is a small drift in the target distribution (the student learns the teacher's samples, compounding the teacher's error, visible in the reflowed 100-step score of 0.0163 versus 0.0123), which is why production systems reflow once or twice, not to convergence.

Worked problems

Implementation

Everything below was run on this machine (H100 80GB HBM3, PyTorch 2.7 / CUDA 12.8, JAX 0.6 on the same GPU) against the two-moons distribution, 200,000 standardized points, a 166k parameter MLP with a 128-dimensional sinusoidal time embedding standing in for the U-Net or DiT, batch 8192. The full measurement set is stored in classes/data/diffusion.json. First, the complete DDPM training loop. Precompute the schedule, sample a random timestep per example, form \( x_t \) with the closed-form marginal, regress on the noise. This loop trained to \( \varepsilon \)-MSE 0.218 in 9.9 s (805 steps/s) in PyTorch and to 0.217 in the JAX version. The two produced statistically indistinguishable samples (sliced Wasserstein 0.0254 versus 0.0235 with a 50-step DDIM sampler).

import math, torch
import torch.nn as nn

T = 1000
beta = torch.linspace(1e-4, 0.02, T, device="cuda")   # linear schedule
alpha = 1.0 - beta
abar = torch.cumprod(alpha, dim=0)                    # closed-form marginal coeffs

class TimeMLP(nn.Module):
    """eps_theta(x, t): 2-D data + 128-dim sinusoidal time embedding."""
    def __init__(self, h=256, emb=128):
        super().__init__()
        half = emb // 2
        self.register_buffer("freqs", torch.exp(
            -math.log(10000.0) * torch.arange(half) / (half - 1)))
        self.net = nn.Sequential(
            nn.Linear(2 + emb, h), nn.SiLU(),
            nn.Linear(h, h), nn.SiLU(),
            nn.Linear(h, h), nn.SiLU(), nn.Linear(h, 2))
    def forward(self, x, t):                          # x: (B,2), t: (B,) float
        ang = t[:, None] * self.freqs[None, :]        # (B, 64)
        temb = torch.cat([ang.sin(), ang.cos()], 1)   # (B, 128)
        return self.net(torch.cat([x, temb], 1))      # (B, 2)

model = TimeMLP().cuda()
opt = torch.optim.Adam(model.parameters(), lr=1e-3)
sched = torch.optim.lr_scheduler.CosineAnnealingLR(opt, 8000)

for step in range(8000):                              # 9.9 s on the H100
    x0 = data[torch.randint(0, len(data), (8192,), device="cuda")]
    t = torch.randint(0, T, (8192,), device="cuda")
    eps = torch.randn_like(x0)
    at = abar[t][:, None]                             # (B, 1)
    xt = at.sqrt() * x0 + (1 - at).sqrt() * eps       # q(x_t | x_0) directly
    loss = ((model(xt, t.float()) - eps) ** 2).mean() # L_simple
    opt.zero_grad(); loss.backward(); opt.step(); sched.step()

import math, jax, jax.numpy as jnp

T = 1000
beta = jnp.linspace(1e-4, 0.02, T)
alpha = 1.0 - beta
abar = jnp.cumprod(alpha)

def time_embed(t, dim=128):                     # t: (B,) float
    half = dim // 2
    freqs = jnp.exp(-math.log(10000.0) * jnp.arange(half) / (half - 1))
    ang = t[:, None] * freqs[None, :]
    return jnp.concatenate([jnp.sin(ang), jnp.cos(ang)], axis=1)

def init_mlp(key, sizes):                       # [130, 256, 256, 256, 2]
    params = []
    for kin, kout in zip(sizes[:-1], sizes[1:]):
        key, k = jax.random.split(key)
        params.append((jax.random.normal(k, (kin, kout)) / jnp.sqrt(kin),
                       jnp.zeros(kout)))
    return params

def eps_model(params, x, t):                    # (B,2),(B,) -> (B,2)
    h = jnp.concatenate([x, time_embed(t)], axis=1)
    for w, b in params[:-1]:
        h = jax.nn.silu(h @ w + b)
    w, b = params[-1]
    return h @ w + b

def loss_fn(params, key, x0):
    kt, ke = jax.random.split(key)
    t = jax.random.randint(kt, (x0.shape[0],), 0, T)
    eps = jax.random.normal(ke, x0.shape)
    at = abar[t][:, None]
    xt = jnp.sqrt(at) * x0 + jnp.sqrt(1.0 - at) * eps
    return jnp.mean((eps_model(params, xt, t.astype(jnp.float32)) - eps) ** 2)

@jax.jit
def adam_step(params, m, v, i, key, x0, lr, b1=0.9, b2=0.999):
    loss, g = jax.value_and_grad(loss_fn)(params, key, x0)
    m = jax.tree.map(lambda a, b: b1 * a + (1 - b1) * b, m, g)
    v = jax.tree.map(lambda a, b: b2 * a + (1 - b2) * b * b, v, g)
    def upd(p, m_, v_):
        mh = m_ / (1 - b1 ** (i + 1)); vh = v_ / (1 - b2 ** (i + 1))
        return p - lr * mh / (jnp.sqrt(vh) + 1e-8)
    return jax.tree.map(upd, params, m, v), m, v, loss

key = jax.random.PRNGKey(0)
params = init_mlp(key, [2 + 128, 256, 256, 256, 2])
m = jax.tree.map(jnp.zeros_like, params)
v = jax.tree.map(jnp.zeros_like, params)
steps = 8000
for i in range(steps):
    key, kb, ks = jax.random.split(key, 3)
    idx = jax.random.randint(kb, (8192,), 0, data.shape[0])
    lr = 0.5e-3 * (1 + math.cos(math.pi * i / steps))   # cosine decay
    params, m, v, loss = adam_step(params, m, v, i, ks, data[idx], lr)

Next the samplers. The ancestral loop is the trained chain run backward. The DDIM function subsumes it (\( \eta = 1 \) on the full grid reproduces ancestral statistics, \( \eta = 0 \) is deterministic and step-skipping), and the classifier-free guidance sampler wraps DDIM around two model evaluations per step with the null-token trick. These passed the following correctness checks. Ancestral 1000-step samples scored sliced-Wasserstein 0.0151 against held-out data (floor 0.0041), DDIM \( \eta = 0 \) at 1000 steps scored 0.0136, and DDIM \( \eta = 1 \) tracked ancestral quality on the full grid (0.0186), degrading faster as steps shrink, exactly as the ODE analysis predicts.

@torch.no_grad()
def sample_ddpm(model, n):
    """Ancestral sampling: T stochastic reverse steps."""
    x = torch.randn(n, 2, device="cuda")
    for t in range(T - 1, -1, -1):
        eps = model(x, torch.full((n,), float(t), device="cuda"))
        ab, ab_prev = abar[t], abar[t - 1] if t > 0 else torch.tensor(1.0)
        mean = (x - beta[t] / (1 - ab).sqrt() * eps) / alpha[t].sqrt()
        if t > 0:                                   # tilde-beta posterior variance
            var = beta[t] * (1 - ab_prev) / (1 - ab)
            x = mean + var.sqrt() * torch.randn_like(x)
        else:
            x = mean
    return x

@torch.no_grad()
def sample_ddim(model, n_steps, n, eta=0.0):
    """DDIM on a sub-grid; eta=0 deterministic, eta=1 ancestral-like."""
    ts = torch.linspace(T - 1, 0, n_steps).long().tolist()
    x = torch.randn(n, 2, device="cuda")
    for i, t in enumerate(ts):
        eps = model(x, torch.full((n,), float(t), device="cuda"))
        ab_t = abar[t]
        ab_p = abar[ts[i + 1]] if i + 1 < len(ts) else torch.tensor(1.0)
        x0_hat = (x - (1 - ab_t).sqrt() * eps) / ab_t.sqrt()   # Tweedie estimate
        sig = eta * ((1 - ab_p) / (1 - ab_t)).sqrt() * (1 - ab_t / ab_p).sqrt()
        x = ab_p.sqrt() * x0_hat + (1 - ab_p - sig**2).clamp_min(0).sqrt() * eps
        if i + 1 < len(ts):
            x = x + sig * torch.randn_like(x)
    return x

@torch.no_grad()
def sample_cfg(model, y, w, n_steps, n):
    """Classifier-free guidance: model(x, t, y) trained with 10% of labels
    replaced by a null token; eps_cfg = (1+w) eps_cond - w eps_uncond."""
    ts = torch.linspace(T - 1, 0, n_steps).long().tolist()
    x = torch.randn(n, 2, device="cuda")
    yv = torch.full((n,), y, device="cuda", dtype=torch.long)
    null = torch.full((n,), NULL_CLASS, device="cuda", dtype=torch.long)
    for i, t in enumerate(ts):
        tv = torch.full((n,), float(t), device="cuda")
        eps = (1 + w) * model(x, tv, yv) - w * model(x, tv, null)
        ab_t = abar[t]
        ab_p = abar[ts[i + 1]] if i + 1 < len(ts) else torch.tensor(1.0)
        x0_hat = (x - (1 - ab_t).sqrt() * eps) / ab_t.sqrt()
        x = ab_p.sqrt() * x0_hat + (1 - ab_p).clamp_min(0).sqrt() * eps
    return x

def sample_ddim(params, key, n, n_steps, eta=0.0):
    ts = jnp.linspace(T - 1, 0, n_steps).astype(jnp.int32)
    key, k0 = jax.random.split(key)
    x = jax.random.normal(k0, (n, 2))
    for i in range(n_steps):
        t = ts[i]
        ab_t = abar[t]
        ab_p = abar[ts[i + 1]] if i + 1 < n_steps else jnp.array(1.0)
        eps = eps_model(params, x, jnp.full((n,), t, dtype=jnp.float32))
        x0_hat = (x - jnp.sqrt(1 - ab_t) * eps) / jnp.sqrt(ab_t)
        sig = (eta * jnp.sqrt((1 - ab_p) / (1 - ab_t))
                   * jnp.sqrt(1 - ab_t / ab_p))
        x = (jnp.sqrt(ab_p) * x0_hat
             + jnp.sqrt(jnp.clip(1 - ab_p - sig**2, 0.0)) * eps)
        if i + 1 < n_steps:
            key, kn = jax.random.split(key)
            x = x + sig * jax.random.normal(kn, x.shape)
    return x

def sample_cfg(params, key, y, w, n, n_steps):
    """eps_model_c(params, x, t, y) with a null label index for dropout."""
    ts = jnp.linspace(T - 1, 0, n_steps).astype(jnp.int32)
    key, k0 = jax.random.split(key)
    x = jax.random.normal(k0, (n, 2))
    yv = jnp.full((n,), y, dtype=jnp.int32)
    null = jnp.full((n,), NULL_CLASS, dtype=jnp.int32)
    for i in range(n_steps):
        t = jnp.full((n,), ts[i], dtype=jnp.float32)
        eps = ((1 + w) * eps_model_c(params, x, t, yv)
               - w * eps_model_c(params, x, t, null))
        ab_t = abar[ts[i]]
        ab_p = abar[ts[i + 1]] if i + 1 < n_steps else jnp.array(1.0)
        x0_hat = (x - jnp.sqrt(1 - ab_t) * eps) / jnp.sqrt(ab_t)
        x = jnp.sqrt(ab_p) * x0_hat + jnp.sqrt(jnp.clip(1 - ab_p, 0.0)) * eps
    return x

The flow-matching objective replaces all schedule bookkeeping with a linear interpolation, and its sampler is a plain Euler loop. The reflow variant differs from base training in one line, the source of pairs. The velocity network converged to MSE 1.373 in 9.7 s. The irreducible part of that number is the variance of \( x_1 - x_0 \) given \( x_t \), which is what reflow removes (post-reflow, trajectories are straight and the conditional target is nearly deterministic).

def fm_loss(model, x1, x0=None):
    """Conditional flow matching, linear (rectified-flow) path.
    x1: data batch (B,2); x0: paired noise for reflow, else fresh noise."""
    if x0 is None:
        x0 = torch.randn_like(x1)              # independent coupling
    t = torch.rand(x1.size(0), 1, device=x1.device)
    xt = (1 - t) * x0 + t * x1                 # straight interpolant
    v_target = x1 - x0                         # constant along each segment
    v = model(xt, t.squeeze(1) * 1000.0)       # reuse sinusoidal embedding
    return ((v - v_target) ** 2).mean()

@torch.no_grad()
def sample_euler(model, n_steps, n, z0=None):
    x = torch.randn(n, 2, device="cuda") if z0 is None else z0.clone()
    dt = 1.0 / n_steps
    for i in range(n_steps):
        t = torch.full((x.size(0),), i * dt, device="cuda")
        x = x + dt * model(x, t * 1000.0)      # dx/dt = v_theta(x, t)
    return x

# Reflow: regenerate the coupling from the trained flow, retrain on pairs.
z_pool = torch.randn(200_000, 2, device="cuda")
x_pool = sample_euler(flow1, 100, None, z0=z_pool)   # (noise, endpoint) pairs
# ... second model trained with fm_loss(model2, x_pool[idx], z_pool[idx])

def fm_loss(params, key, x1, x0=None):
    k0, kt = jax.random.split(key)
    if x0 is None:
        x0 = jax.random.normal(k0, x1.shape)   # independent coupling
    t = jax.random.uniform(kt, (x1.shape[0], 1))
    xt = (1 - t) * x0 + t * x1                 # straight interpolant
    v = eps_model(params, xt, t[:, 0] * 1000.0)  # same MLP, now predicts v
    return jnp.mean((v - (x1 - x0)) ** 2)

def sample_euler(params, key, n, n_steps, z0=None):
    x = jax.random.normal(key, (n, 2)) if z0 is None else z0
    dt = 1.0 / n_steps
    for i in range(n_steps):
        t = jnp.full((x.shape[0],), i * dt)
        x = x + dt * eps_model(params, x, t * 1000.0)
    return x

# Reflow: pair each noise draw with its own generated endpoint.
z_pool = jax.random.normal(jax.random.PRNGKey(7), (200_000, 2))
x_pool = sample_euler(params1, None, 200_000, 100, z0=z_pool)
# ... retrain with fm_loss(params2, key, x_pool[idx], x0=z_pool[idx])

Finally the DiT block with adaLN-zero, the conditioning mechanism examined in the architecture section below. The timestep and class/text-pooled embeddings are summed into one conditioning vector \( c \). A per-block MLP maps \( c \) to six vectors, shift, scale, and gate for each of the two sublayers. The final projection is initialized to zero, so every gate starts at zero and the block is exactly the identity at initialization. I verified both properties on this machine. Output equals input at initialization to float precision, and shapes are \( (B, N, d) \to (B, N, d) \) with \( c \in \R^{B \times d} \).

class DiTBlock(nn.Module):
    """DiT block with adaLN-zero conditioning (Peebles & Xie 2023)."""
    def __init__(self, d, n_head, mlp_ratio=4):
        super().__init__()
        self.norm1 = nn.LayerNorm(d, elementwise_affine=False, eps=1e-6)
        self.attn = nn.MultiheadAttention(d, n_head, batch_first=True)
        self.norm2 = nn.LayerNorm(d, elementwise_affine=False, eps=1e-6)
        self.mlp = nn.Sequential(
            nn.Linear(d, mlp_ratio * d), nn.GELU(approximate="tanh"),
            nn.Linear(mlp_ratio * d, d))
        self.adaLN = nn.Sequential(nn.SiLU(), nn.Linear(d, 6 * d))
        nn.init.zeros_(self.adaLN[-1].weight)   # the "-zero": gates start at 0,
        nn.init.zeros_(self.adaLN[-1].bias)     # so the block is the identity

    def forward(self, x, c):                    # x: (B,N,d) tokens, c: (B,d)
        sh1, sc1, g1, sh2, sc2, g2 = self.adaLN(c).chunk(6, dim=1)
        h = self.norm1(x) * (1 + sc1[:, None]) + sh1[:, None]
        x = x + g1[:, None] * self.attn(h, h, h, need_weights=False)[0]
        h = self.norm2(x) * (1 + sc2[:, None]) + sh2[:, None]
        x = x + g2[:, None] * self.mlp(h)
        return x

blk = DiTBlock(384, 6).cuda()
x = torch.randn(4, 256, 384, device="cuda")     # 256 latent patch tokens
c = torch.randn(4, 384, device="cuda")          # timestep + label embedding
out = blk(x, c)
assert out.shape == (4, 256, 384)
assert torch.allclose(out, x)                   # identity at init: verified

def init_dit_block(key, d, mlp_ratio=4):
    k1, k2, k3 = jax.random.split(key, 3)
    s = 1.0 / jnp.sqrt(d)
    return {
        "wqkv": jax.random.normal(k1, (d, 3 * d)) * s,
        "wo":   jax.random.normal(k2, (d, d)) * s,
        "w1":   jax.random.normal(k3, (d, mlp_ratio * d)) * s,
        "b1":   jnp.zeros(mlp_ratio * d),
        "w2":   jnp.zeros((mlp_ratio * d, d)), "b2": jnp.zeros(d),
        "ada_w": jnp.zeros((d, 6 * d)),        # adaLN-zero: zero-init
        "ada_b": jnp.zeros(6 * d),
    }

def layer_norm(x, eps=1e-6):                   # no learned affine in DiT
    mu = x.mean(-1, keepdims=True)
    var = x.var(-1, keepdims=True)
    return (x - mu) / jnp.sqrt(var + eps)

def dit_block(p, x, c, n_head):                # x: (B,N,d), c: (B,d)
    ada = jax.nn.silu(c) @ p["ada_w"] + p["ada_b"]
    sh1, sc1, g1, sh2, sc2, g2 = jnp.split(ada, 6, axis=-1)
    h = layer_norm(x) * (1 + sc1[:, None]) + sh1[:, None]
    B, N, d = h.shape
    qkv = (h @ p["wqkv"]).reshape(B, N, 3, n_head, d // n_head)
    q, k, v = jnp.moveaxis(qkv, 2, 0)          # each (B, N, H, dh)
    att = jax.nn.dot_product_attention(q, k, v)
    x = x + g1[:, None] * (att.reshape(B, N, d) @ p["wo"])
    h = layer_norm(x) * (1 + sc2[:, None]) + sh2[:, None]
    h = jax.nn.gelu(h @ p["w1"] + p["b1"]) @ p["w2"] + p["b2"]
    return x + g2[:, None] * h                 # identity at init: gates are 0

Measured results

The complete step-count sweeps, all with 20,000 generated points scored by sliced 2-Wasserstein distance against 20,000 held-out points. Two independent draws of real data score 0.0041 against each other. That is the floor any sampler is chasing. Lower is better throughout.

Diffusion step counts (1000/50/20/10/5/2) apply to the DDIM columns. Flow columns use the Euler counts (100/32/16/8/4/2/1). Three structural facts are visible. Deterministic beats stochastic at every reduced step count. The un-reflowed flow model degrades with the same shape as DDIM, because both are first-order solvers on curved trajectories, though the flow's flatter path buys roughly a 2× step advantage in the mid-range. And the reflowed model is essentially flat from 100 steps down to one. At 1 NFE it scores 0.0191, better than DDIM achieves at 20 steps, because after straightening (curvature metric 0.914 → 0.0004) a single Euler step integrates the ODE almost exactly. This 60-second experiment reproduces, in miniature, the entire trajectory of the few-step-generation literature.

How it is done in practice

Architectures: the U-Net and its inductive biases

The denoiser in DDPM through Stable Diffusion XL is a U-Net, an encoder that halves resolution while widening channels, a bottleneck, and a decoder that mirrors the encoder with skip connections carrying each resolution's features straight across. The residual blocks receive the timestep embedding through per-block scale-shift (FiLM-style) modulation, and at the lower resolutions self-attention layers let distant regions coordinate, with cross-attention to the text encoder interleaved at the same depths.

x_t (64x64xC)
  │ ResBlock+Attn ──────────────────────────────► skip ─┐
  ▼ down                                                │
 32x32x2C ── ResBlock+Attn+CrossAttn ─────────► skip ─┐ │
  ▼ down                                              │ │
 16x16x4C ── ResBlock+Attn+CrossAttn ───────► skip ─┐ │ │
  ▼ down                                            │ │ │
  8x8x4C ─── bottleneck: ResBlock + Attn            │ │ │
  ▲ up   ◄──────────────── concat ◄─────────────────┘ │ │
 16x16 ...   (t-embedding modulates every ResBlock)   │ │
  ▲ up   ◄──────────────── concat ◄───────────────────┘ │
 32x32 ...                                              │
  ▲ up   ◄──────────────── concat ◄─────────────────────┘
eps_theta (64x64xC)

Its inductive biases are convolutional locality (edges and textures come nearly free), multi-scale processing matched to the coarse-to-fine order in which diffusion resolves an image, and skip connections that let the high-SNR steps act as a near-identity with small corrections. Those biases buy sample efficiency at small scale and become a straitjacket at large scale. Channel widths and attention placements are hand-tuned per resolution, global context only exists at the bottleneck, and the architecture has no clean scaling knob comparable to "add layers, widen, train longer."

Diffusion transformers and adaLN-zero

Peebles and Xie's DiT (2023) discards all of it. Patchify the latent into tokens (patch 2 over a 32×32×4 latent gives 256 tokens), run a standard ViT-style transformer, and unpatchify to predict noise. The subtle part is conditioning. They compared in-context conditioning (append timestep and class as tokens), cross-attention, and adaptive LayerNorm, and the winner, adaLN-zero, is precise in a way worth spelling out. Each block's LayerNorms carry no learned affine parameters. Instead a small MLP maps the conditioning vector \( c \) (timestep embedding plus class or pooled-text embedding) to six per-channel vectors \( (\gamma_1, \beta_1, g_1, \gamma_2, \beta_2, g_2) \). The block computes \( x \mathrel{+}= g_1 \odot \mathrm{Attn}\big( (1+\gamma_1) \odot \mathrm{LN}(x) + \beta_1 \big) \) and likewise for the MLP sublayer. The "-zero" is the initialization. The modulation MLP's output layer starts at exactly zero, so every gate \( g \) is zero and every block is the identity at initialization, the whole network computing zero residual regardless of depth. This mirrors the zero-init residual-branch trick from large ResNet and GPT training and measurably beats plain adaLN and cross-attention conditioning at equal FLOPs. My implementation above verifies the identity property exactly. The scaling evidence is what moved the field. Across DiT-S through DiT-XL, FID tracks model GFLOPs almost monotonically regardless of how compute is arranged (depth, width, or token count), giving diffusion the same "spend more, predictably get more" scaling law language models enjoy, with none of the U-Net's per-resolution hand-tuning. Every current frontier system, SD3's MMDiT, Flux, and the video transformers, is a descendant. The MMDiT variant runs separate parameter sets for text and image tokens joined by shared attention, which Esser et al. found clearly better than a single shared stream at equal parameters.

Latent diffusion: the autoencoder stage

Rombach et al.'s latent diffusion (2022) split generation into a perceptual compression stage and a semantic generation stage, on the observation that most of an image's bits encode imperceptible detail that a likelihood-weighted diffusion model wastes capacity on. The autoencoder (f = 8 downsampling, 4 or more recently 16 channels) is trained once, with three losses, pixel reconstruction, a perceptual (LPIPS) loss computed as a distance in the feature space of a pretrained network, and a patch adversarial loss. The combination is load-bearing. A plain L2 autoencoder at 48× compression (\( 1024^2 \times 3 \to 128^2 \times 4 \)) averages over the fine detail it cannot store and reconstructs blur. The perceptual loss makes it preserve what feature detectors respond to, and the adversarial critic forces the decoder to resynthesize plausible texture rather than the conditional mean. A tiny KL penalty (weight around \( 10^{-6} \)) keeps the latent distribution loosely Gaussian without imposing a real VAE bottleneck. The diffusion model then never sees a pixel.

train:  x ──E──► z (128x128x4) ──[diffusion trains here]
sample: N(0,I) ──sampler──► z_hat ──D──► x_hat (1024x1024x3)
        text ──T5/CLIP──► tokens ──cross-attn──► every block

The arithmetic in worked problem 5 quantifies the payoff, roughly 1500× less transformer compute per sampling pass at \( 1024^2 \) for the dominant attention-bearing stage, which is the entire difference between interactive generation and batch rendering. The trend is toward heavier latents. SD3 and Flux moved from 4 to 16 latent channels because 4-channel latents measurably lose text glyphs and fine structure the decoder cannot re-hallucinate, and video models add 4× to 8× temporal compression in a causal 3D VAE. The autoencoder's reconstruction ceiling silently upper-bounds every downstream benchmark number, a fact worth remembering when two systems with different VAEs are compared on FID.

Conditioning mechanisms

Text conditioning enters through cross-attention. Image tokens form queries, the text encoder's output sequence forms keys and values, so each spatial location retrieves the prompt content it needs. The choice of encoder shows up directly in capability. CLIP encoders (contrastively trained, used alone in SD1/2) give strong style and object-category signal but a bag-of-words-like representation that mangles attribute binding and word order. T5-class encoders (span-corruption language models, adopted by Imagen and standard since SD3, which concatenates two CLIPs with T5-XXL) carry compositional and spelling information, and ablations in both Imagen and SD3 attribute most of the gain in text rendering and complex binding to the language-model encoder, with SD3 showing T5 can even be dropped at inference for a modest typography-heavy penalty. Current systems increasingly condition on decoder-only LLM embeddings or full VLM rewrites of the prompt (DALL-E 3's recaptioning result showed dense synthetic captions alone lift prompt following substantially).

Structural control composes with the same backbone. Zhang et al.'s ControlNet (2023) clones the U-Net encoder into a trainable copy that ingests a spatial map (edges, depth, pose), injecting its features into the frozen base through zero-initialized 1×1 convolutions, the same start-as-identity trick as adaLN-zero, so training cannot initially damage the base model. T2I-Adapters achieve similar control with a far smaller side network. IP-Adapter conditions on a reference image by adding a parallel cross-attention branch over CLIP image embeddings (about 22M parameters), giving subject or style transfer without fine-tuning. Image-to-image needs no training at all. Noise the input to an intermediate \( t \) (strength 0.3 to 0.8) and denoise from there, trading faithfulness against edit magnitude along the same SNR axis derived above. Inpainting either resamples with the known region clamped to its noised ground truth at every step (with a resampling correction, as in RePaint from Lugmayr et al. at ETH Zurich) or, in production, fine-tunes with the mask and masked image as extra input channels. LoRA fine-tuning of the attention projections is the dominant personalization mechanism, typically rank 4 to 64 on the cross- and self-attention weights.

Evaluation, and why the numbers are shaky

The standard metric is Fréchet Inception Distance (Heusel et al., 2017). Embed 50k real and 50k generated images with an Inception-V3 trained on ImageNet, fit Gaussians to both clouds, and report the closed-form Fréchet distance \( \|\mu_r - \mu_g\|^2 + \tr\big( \Sigma_r + \Sigma_g - 2(\Sigma_r \Sigma_g)^{1/2} \big) \). Its pathologies are well documented. The embedding is an ImageNet classifier, so FID is most sensitive to the categories and textures that network cares about. Kynkäänniemi et al. (2023) showed FID can be reduced substantially by nudging generations toward ImageNet class statistics with no human-visible quality change, and Stein et al. (2023) found DINOv2 embeddings rank models closer to human judgment. The estimator is also biased with strong sample-size dependence (Chong and Forsyth, 2020). FID at 10k images is systematically higher than at 50k, so cross-paper comparisons at unstated sample counts are meaningless, and resize/JPEG preprocessing differences alone shift scores by more than many claimed improvements (Parmar et al.'s clean-FID work). A single distance also conflates two failures, so Kynkäänniemi et al.'s precision/recall (2019) splits them with k-NN manifold estimates. Precision is the fraction of generated samples inside the real manifold (fidelity), recall the fraction of real samples inside the generated manifold (coverage). Guidance visibly trades recall for precision in exactly the way the Gaussian computation in problem 4 predicts. Text-image alignment is scored by CLIP similarity, which saturates and can be gamed, so the field increasingly leans on human-preference models (HPSv2, PickScore, ImageReward), structured prompt suites (PartiPrompts, GenEval, T2I-CompBench), and human arenas. Within one codebase, FID deltas are informative. Across papers, treat every table as directional at best, and trust blinded human comparisons over all of it.

Efficiency and deployment

Take the deployment arithmetic for a Flux-class model. In bf16, 12B parameters is 24 GB of weights, plus the T5 encoder (4.7B, 9.4 GB) and VAE, so the model fits an 80 GB H100 with abundant activation room but needs quantization for consumer 24 GB cards. Unlike LLM decoding, diffusion sampling is compute-bound. Every step reprocesses all tokens (4096 image tokens plus text at \( 1024^2 \), no KV-cache reuse across steps), so the levers are different from language serving. The measured attention numbers from this machine's benchmark set show why fused attention is non-negotiable at video scale. At 8192 tokens, naive attention takes 98.9 ms and 35 GB of peak memory against FlashAttention's 3.6 ms and 0.57 GB, and at 16k tokens (a modest video clip's spatiotemporal token count) the naive path cannot run at all while the fused kernel takes 13.7 ms at 640 TFLOPS.

The remaining levers, in rough order of leverage. Fewer steps via distillation (next section), the largest single win, 25 to 50×. Step caching relies on the fact that adjacent steps produce nearly identical intermediate features, so DeepCache (U-Nets) and TeaCache/FORA (DiTs) reuse deep features across steps and skip most of the network on a schedule, for 2 to 4× at mild quality cost. In quantization, weight-only int8 is essentially free. Activation quantization is harder than in LLMs because activation statistics drift across timesteps, which Q-Diffusion handles with timestep-aware calibration, and MIT-Han-lab's SVDQuant (2024) reaches 4-bit weights and activations on Flux by absorbing outliers into a low-rank branch, reporting about 3.5× memory and 3× latency improvement on consumer GPUs. Token merging exploits spatial redundancy, and compilation adds tens of percent. Stacked, a 1024×1024 generation that cost 5 to 10 s at SDXL launch runs well under a second today, and the distilled tier ships real-time preview-as-you-type products.

The current research frontier

Few-step generation: distilling the trajectory

All few-step methods start from the same fact, demonstrated twice above (problem 3 analytically, the measured sweep empirically). The bottleneck is not the model's knowledge but the curvature of the path a first-order solver must follow. Salimans and Ho's progressive distillation (2022) halves the path repeatedly. A student is trained so one of its DDIM steps reproduces two of the teacher's, then the student becomes the teacher, so \( N \) rounds cut \( 2^N \) steps to one. With v-parameterization (required, since an epsilon student cannot represent a large jump from pure noise) they reached 4-step samplers matching 8192-step teachers on CIFAR-scale FID. Its drawbacks are the serial rounds and error compounding across them.

Song et al.'s consistency models (2023) collapse the recursion into a single self-consistency condition. Learn \( f_\theta(x_t, t) \) mapping any point on a PF-ODE trajectory directly to the trajectory's origin, so \( f_\theta(x_t, t) = f_\theta(x_{t'}, t') \) for any two points on the same trajectory, with the boundary condition \( f_\theta(x, t_{\min}) = x \) enforced architecturally by the parameterization \( f_\theta = c_{\text{skip}}(t)\, x + c_{\text{out}}(t)\, F_\theta(x, t) \) with \( c_{\text{skip}}(t_{\min}) = 1, c_{\text{out}}(t_{\min}) = 0 \). Consistency distillation trains by taking adjacent points. From \( x_{t_{n+1}} \), one teacher ODE step estimates \( \hat x_{t_n} \), and the loss is \( d\big( f_\theta(x_{t_{n+1}}, t_{n+1}),\ f_{\theta^-}(\hat x_{t_n}, t_n) \big) \) with \( \theta^- \) an EMA copy as target, chaining local consistency into global. Consistency training removes the teacher entirely by substituting the unbiased trajectory estimate available from the data. Since \( x_t = x_0 + t z \) (VE form), the pair \( (x_0 + t_{n+1} z,\ x_0 + t_n z) \) with shared \( z \) lies approximately on a common trajectory, and the same loss applies, making a standalone one-step generative model trained from scratch. One step gives usable samples, two steps (generate, re-noise, regenerate) recovers much of the remaining gap. Latent consistency models (Luo et al., 2023) applied distillation in Stable Diffusion's latent space with guidance baked in, giving the 4-step SD workflows now ubiquitous, and LoRA-only LCM distillation made it a downloadable plugin. The successor line (Kim et al.'s CTM, improved/easy consistency tuning, and MeanFlow from Geng et al. 2025, a CMU and MIT collaboration that regresses the average velocity over an interval via an identity relating it to the instantaneous field) now posts one-step ImageNet FIDs competitive with multi-step teachers. Shortcut models (Frans et al., 2024) condition a rectified-flow network on the intended step size so one model serves every step budget.

The adversarial branch trades likelihood grounding for sharper single steps. Sauer et al.'s adversarial diffusion distillation (SDXL-Turbo, 2023) combines score distillation from the teacher with a feature-space discriminator, producing 1-step samples and 4-step samples preferred over the 50-step teacher in human studies. Latent adversarial distillation (LADD, behind SD3-Turbo) and Flux's schnell variant industrialized the recipe. What is lost at 1 to 4 steps is consistent across methods and visible in my toy sweep's residual gap. The losses are fine texture and high-frequency detail (adversarial variants mask this best), sample diversity (distilled students mode-seek, and baked-in guidance compounds it), the compute-quality dial itself, and, for adversarial students, the clean correspondence between the model and any likelihood or score interpretation, which matters the moment you want inversion or editing.

Video and 3D

Video diffusion began by inflating image models. Blattmann et al.'s video LDM (2023) froze a text-to-image U-Net and interleaved temporal layers (temporal convolutions and attention over the frame axis, with attention factorized as spatial-then-temporal to keep cost near \( O(F \cdot N^2 + N \cdot F^2) \) rather than \( O((NF)^2) \)), plus a temporal VAE decoder to remove flicker. Stable Video Diffusion scaled the recipe with a three-stage data curriculum. The Sora technical report (Brooks et al., 2024) marked the architectural break. Encode video into spacetime patches with a causal 3D VAE and run a single diffusion transformer with full 3D attention over all spatiotemporal tokens, trained across durations, resolutions, and aspect ratios. Open implementations of the recipe (CogVideoX at Zhipu, HunyuanVideo at Tencent, Wan at Alibaba, Mochi at Genmo) converged on the same shape, a 3D-VAE with roughly 8×8×4 spatial-temporal compression, a rectified-flow DiT of 2B to 14B parameters, full or windowed 3D attention, and T5/LLM text conditioning. The token counts explain the engineering pressure. Five seconds of \( 720p \) at 16 latent fps is tens of thousands of tokens, which is why the measured 16k-token attention number above (13.7 ms fused, impossible naive) is the relevant regime, and why sparse/windowed attention and step caching matter more for video than for images.

3D generation via 2D priors is built on score distillation sampling from Poole et al.'s DreamFusion (2022). A differentiable renderer \( g(\theta, \pi) \) (NeRF then, 3D Gaussian splatting now) produces an image from parameters \( \theta \) at camera \( \pi \). The diffusion model judges noised renders. Formally, take the diffusion training loss with \( x = g(\theta) \), \( \L(\theta) = \E_{t, \varepsilon}\big[ w(t) \| \varepsilon_\phi(\alpha_t g(\theta) + \sigma_t \varepsilon;\, y, t) - \varepsilon \|^2 \big] \), and differentiate through the chain,

The U-Net Jacobian is dropped for two stated reasons. It costs a backward pass through the diffusion network per render, and empirically it is poorly conditioned at the noise levels that matter. The move is principled rather than ad hoc. Poole et al. show the Jacobian-free gradient is exactly the gradient of the KL divergence \( \KL\big( q(x_t \mid g(\theta)) \,\|\, p_\phi(x_t \mid y) \big) \) (the score of the Gaussian \( q \) contributes the \( -\varepsilon \) term, and the problematic Jacobian appears only in a term whose expectation vanishes), i.e. probability density distillation, mode-seeking by construction. That explains the characteristic SDS failures, over-saturated colors and low diversity, worsened by the very high guidance (\( w \approx 100 \)) needed to make the mode sharp, and the Janus problem of a face on every side, which is a prior mismatch across cameras rather than an optimization bug. ProlificDreamer's variational score distillation replaces the fixed noise target with a learned score of the render distribution, restoring diversity and normal guidance scales. The renderer side moved to 3D Gaussian splatting (DreamGaussian and successors) for 10 to 100× faster iteration, and the frontier is bypassing per-scene optimization entirely. Large reconstruction models and multi-view diffusion (Zero123, MVDream, and their industrial descendants) amortize 3D generation into a feed-forward pass, with SDS surviving as a refinement stage.

Autoregressive image generation and unified models

The competing lineage tokenizes images and applies next-token prediction, VQ-VAE/VQGAN discretization, then Parti (20B, 2022) and LlamaGen showing plain AR transformers scale on image tokens. Two 2024 results made the comparison serious again. Tian et al.'s VAR reframed AR generation as next-scale prediction (predict a coarse token grid, then successively finer ones), restoring a coarse-to-fine order that raster-scan AR lacks and beating DiT baselines on ImageNet at lower inference cost. Li et al.'s MAR dropped quantization, using masked autoregression over continuous tokens with a small per-token diffusion head for the output distribution, outperforming both discrete AR and plain diffusion at equal compute and clarifying that "AR versus diffusion" is really two independent axes, factorization order versus per-step output model. The practical trade today is that diffusion/flow still holds photorealistic texture and resolution economics (attention over 4k latent tokens beats sequential generation of 4k tokens for wall-clock), while AR composes naturally with LLM training infrastructure, interleaved multimodal context, and instruction following. Unified multimodal models bet on that composability. Chameleon (Meta) trains one early-fusion token stream for text and images. Transfusion (Meta) mixes next-token loss on text with diffusion loss on image latents in one transformer, the recipe behind several production systems. DeepSeek's Janus-Pro decouples the understanding and generation vision encoders. ByteDance's BAGEL and Alibaba's Qwen-Image push open-weight unified models with strong text rendering. And GPT-4o and Gemini's native image generation demonstrated the in-context editing and world-knowledge benefits of generating inside the language model. The open question is whether generation should live inside the reasoning model or behind an interface. Current evidence suggests hybrid designs, an LLM planning for a flow-based renderer, will persist for a while.

Open source to read

huggingface/diffusers is the reference map of the whole ecosystem, with every scheduler, guidance trick, and adapter in one API. Open src/diffusers/schedulers/scheduling_ddpm.py first and match each line against the ancestral update on this page, then read scheduling_dpmsolver_multistep.py with the exponential-integrator section in hand.

lucidrains/denoising-diffusion-pytorch is the cleanest end-to-end DDPM to study whole. The single file denoising_diffusion_pytorch/denoising_diffusion_pytorch.py contains schedule, U-Net, losses, and both samplers. NVlabs/edm is the Karras et al. reformulation. Start at training/loss.py to see the \( \sigma \)-space preconditioning and weighting made concrete, then generate.py for the Heun sampler. crowsonkb/k-diffusion is where many production samplers were first implemented. k_diffusion/sampling.py is a one-file catalog of Euler, Heun, DPM-Solver++ and their SDE variants.

CompVis/latent-diffusion is the historical source of Stable Diffusion. Read ldm/models/diffusion/ddpm.py for the training wrapper and the autoencoder losses discussed above. Its successor Stability-AI/generative-models (SDXL, SVD) modularizes the same design. Begin with sgm/modules/diffusionmodules/sampling.py. facebookresearch/DiT holds the diffusion transformer exactly as published. models.py contains the adaLN-zero block this page's implementation mirrors. black-forest-labs/flux shows the current rectified-flow frontier in readable form. src/flux/model.py is the MMDiT-style backbone with double- and single-stream blocks.

ashawkey/stable-dreamfusion reimplements score distillation against Stable Diffusion. The SDS gradient with its dropped Jacobian is a dozen lines in guidance/sd_utils.py. huggingface/peft is where diffusion LoRA lives in practice. Read src/peft/tuners/lora/layer.py to see exactly what a rank-16 adapter does to a cross-attention projection. And comfyanonymous/ComfyUI is the de facto production sampler graph. comfy/samplers.py shows how guidance, samplers, and schedules compose when users chain them arbitrarily.

Common misconceptions

"The simple loss is the ELBO." It is the ELBO with its per-timestep weights deliberately discarded, which changes the optimum under finite capacity. Training with the exact bound improves log-likelihood and worsens FID. The reweighting shifts capacity from the near-clean steps that dominate likelihood to the mid-noise steps that determine perceptual structure. Weighting is a modeling decision, and VDM, EDM, min-SNR, and SD3's timestep sampling are all different stances on it.

"The network removes the noise at each step." The network computes a conditional expectation, \( \E[\varepsilon \mid x_t] \), equivalently the score. Its one-shot clean-image estimate is a posterior mean, provably blurry (problem 3 exhibits the collapse exactly, with a single step retaining 0.6% of the data variance). Detail emerges only from iterating estimate, step, re-estimate, which is why "just predict \( x_0 \) once" is not a sampler and why few-step models must be trained, not merely stepped, into few steps.

"DDIM skipping steps is an approximation of the trained chain." Backwards. The training objective only constrains marginals, DDIM is an exact reverse process for a different (non-Markovian) forward family with the same marginals, and its deterministic limit is a first-order ODE solver. Step skipping is not a hack on the chain. The chain was one arbitrary member of the family.

"More steps always help." Above the point where discretization error falls below model error, extra steps do nothing. My measured sweep is flat from 200 to 1000 steps (0.0155 to 0.0136 against a floor of 0.0041), with all remaining error coming from the network, not the solver. Meanwhile added stochasticity helps only when steps are plentiful. At 20 steps the \( \eta = 1 \) sampler is 42% worse than \( \eta = 0 \).

"Guidance scale is a free quality knob." The guided score targets a tilted distribution whose variance contracts and whose mean overshoots the conditional (computed exactly in problem 4, measured in the CFG sweep). Moderate scales buy adherence. High scales buy saturation, diversity collapse, and eventually instability, and every production fix (thresholding, rescaling, guidance intervals, autoguidance) is a way of spending the sharpening where it helps.

"Flow matching is a different theory that replaced diffusion." Flow matching with the VP Gaussian path reproduces the probability-flow ODE of diffusion exactly, as derived above. The linear path is a different member of the same family with lower curvature and better-conditioned targets. The frontier models switched paths and objectives, not mathematical frameworks, and score, epsilon, v, and velocity remain affine translations of one another.

Self-check

References

1. Bishop and Bishop, Deep Learning: Foundations and Concepts, Springer, 2024. Chapter 20 is a clean textbook treatment of diffusion. bishopbook.com
2. Sohl-Dickstein, Weiss, Maheswaranathan, Ganguli. Deep Unsupervised Learning using Nonequilibrium Thermodynamics. 2015. arXiv:1503.03585
3. Ho, Jain, Abbeel. Denoising Diffusion Probabilistic Models. 2020. arXiv:2006.11239
4. Nichol, Dhariwal. Improved Denoising Diffusion Probabilistic Models. 2021. arXiv:2102.09672
5. Song, Ermon. Generative Modeling by Estimating Gradients of the Data Distribution. 2019. arXiv:1907.05600
6. Song, Sohl-Dickstein, Kingma, Kumar, Ermon, Poole. Score-Based Generative Modeling through Stochastic Differential Equations. 2021. arXiv:2011.13456
7. Song, Meng, Ermon. Denoising Diffusion Implicit Models. 2021. arXiv:2010.02502
8. Dhariwal, Nichol. Diffusion Models Beat GANs on Image Synthesis. 2021. arXiv:2105.05233
9. Ho, Salimans. Classifier-Free Diffusion Guidance. 2022. arXiv:2207.12598
10. Kingma, Salimans, Poole, Ho. Variational Diffusion Models. 2021. arXiv:2107.00630
11. Vincent. A Connection Between Score Matching and Denoising Autoencoders. Neural Computation, 2011. doi:10.1162/NECO_a_00142
12. Karras, Aittala, Aila, Laine. Elucidating the Design Space of Diffusion-Based Generative Models. 2022. arXiv:2206.00364
13. Rombach, Blattmann, Lorenz, Esser, Ommer. High-Resolution Image Synthesis with Latent Diffusion Models. 2022. arXiv:2112.10752
14. Peebles, Xie. Scalable Diffusion Models with Transformers. 2023. arXiv:2212.09748
15. Esser, Kulal, Blattmann, et al. Scaling Rectified Flow Transformers for High-Resolution Image Synthesis (Stable Diffusion 3). 2024. arXiv:2403.03206
16. Lipman, Chen, Ben-Hamu, Nickel, Le. Flow Matching for Generative Modeling. 2023. arXiv:2210.02747
17. Liu, Gong, Liu. Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow. 2022. arXiv:2209.03003
18. Albergo, Vanden-Eijnden. Building Normalizing Flows with Stochastic Interpolants. 2023. arXiv:2209.15571
19. Salimans, Ho. Progressive Distillation for Fast Sampling of Diffusion Models. 2022. arXiv:2202.00512
20. Song, Dhariwal, Chen, Sutskever. Consistency Models. 2023. arXiv:2303.01469
21. Lu, Zhou, Bao, Chen, Li, Zhu. DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps. 2022. arXiv:2206.00927
22. Poole, Jain, Barron, Mildenhall. DreamFusion: Text-to-3D using 2D Diffusion. 2022. arXiv:2209.14988
23. Zhang, Rao, Agrawala. Adding Conditional Control to Text-to-Image Diffusion Models (ControlNet). 2023. arXiv:2302.05543
24. Heusel, Ramsauer, Unterthiner, Nessler, Hochreiter. GANs Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium (FID). 2017. arXiv:1706.08500
25. Kynkäänniemi, Karras, Laine, Lehtinen, Aila. Improved Precision and Recall Metric for Assessing Generative Models. 2019. arXiv:1904.06991