August 19, 2026

The Evolution of Generative Models

From counting to energy, to latents, to adversaries, to autoregression, to score and flow — one problem, seven answers.

Every generative model in this post is an answer to the same question: given samples from an unknown distribution, how do you build a machine that produces new samples from it? What changes across seventy years is what object the model represents — a table of counts, an energy, a latent-variable decoder, a sampler, a conditional, a score, or a velocity field — and each choice determines the training objective and the sampling procedure that follow from it.

My aim here is to write each family in a single consistent notation so the lineage is visible: EM's ELBO is literally the VAE's ELBO; the Langevin sampler that made Boltzmann machines intractable is the sampler that makes diffusion work; diffusion is a hierarchical VAE with the encoder frozen; flow matching is a continuous normalizing flow with the simulation removed from training. These are not analogies. They are the same equations, reused.

1. The problem and the notation

Fix a data space $\mathcal{X}$ — $\mathbb{R}^D$ for images and audio, $\mathcal{V}^T$ (sequences of tokens from a finite vocabulary $\mathcal{V}$) for text. There is a true but unknown distribution $p_{\text{data}}$ on $\mathcal{X}$, and we observe a finite dataset

A generative model is a parametric family $\{p_\theta : \theta \in \Theta\}$ together with a procedure for drawing $x \sim p_\theta$. Two capabilities are worth distinguishing, because most of the history is about trading one for the other:

• Density evaluation: given $x$, compute (or lower-bound) $\log p_\theta(x)$.
• Sampling: draw fresh $x \sim p_\theta$.

A model that can do both cheaply and exactly is rare. Autoregressive models get exact density but slow sampling; GANs get fast sampling but no density at all; energy-based models get a flexible unnormalized density but neither exact likelihood nor fast sampling.

Notation used throughout

One convention worth stating explicitly because it causes endless confusion: in this post "score" always means $\nabla_x \log p(x)$, the gradient with respect to the data, not the statistician's score $\nabla_\theta \log p_\theta(x)$. The whole diffusion literature rests on the first meaning.

2. One objective and three obstacles

Nearly every model before 2014 — and most after — is trained by maximum likelihood:

This is not an arbitrary choice. As $N \to \infty$ the empirical average converges to $\mathbb{E}_{p_{\text{data}}}[\log p_\theta(x)]$, and

so maximizing likelihood is exactly minimizing the forward KL divergence $D_{\mathrm{KL}}(p_{\text{data}} \Vert p_\theta)$. The direction matters. Forward KL pays an infinite penalty wherever $p_{\text{data}} > 0$ but $p_\theta \to 0$, so maximum-likelihood models are mode-covering: they would rather smear probability over empty regions than drop a mode. This single fact explains why VAEs produce blurry images, why language models assign non-trivial probability to nonsense, and why GANs — which minimize a different divergence — behave in the opposite way (sharp, but mode-dropping).

If maximum likelihood is the goal, why are there so many families? Because $\log p_\theta(x)$ is hard to write down. Exactly three obstacles arise, and each family is defined by which one it accepts and how it dodges the others.

1. The normalization obstacle. If you parameterize an unnormalized score $\tilde p_\theta(x) = e^{-E_\theta(x)}$, then $\log p_\theta(x) = -E_\theta(x) - \log Z(\theta)$ and $Z(\theta)$ is a $D$-dimensional integral you cannot compute. Energy-based models live here.
2. The marginalization obstacle. If you introduce latents, $p_\theta(x) = \int p_\theta(x \mid z)p(z)\,dz$ is again an intractable integral, and the posterior $p_\theta(z \mid x)$ needed to reason about it is intractable too. Mixture models, VAEs, and diffusion live here.
3. The expressiveness obstacle. You can sidestep both by restricting the architecture so that $\log p_\theta(x)$ is closed-form — a factorized chain rule (autoregressive) or an invertible map (normalizing flows) — but you pay in architectural freedom or sampling speed.

And there is a fourth path, taken in 2014: give up on likelihood entirely. Define the model implicitly as a sampler $x = G_\theta(z)$ and train it by comparing samples to data. That is the GAN, and it is why it needs an entirely different theory.

The lineage at a glance

3. Era I — Classical probabilistic models (1950s–2000s)

The classical era establishes the three templates that everything else inherits: the chain rule, the ELBO, and the Gibbs distribution. There are no deep networks, but almost every modern objective is already present in embryo.

3.1 Count-based models: the chain rule with a lookup table

Setup. Data are sequences $x = (x_1,\dots,x_T)$ with $x_t \in \mathcal{V}$. The chain rule of probability is exact and assumption-free:

The problem is that $x_{\lt t}$ has $|\mathcal{V}|^{t-1}$ possible values. The $n$-gram assumption truncates the context to the last $n-1$ symbols:

with parameters $\theta = \{\theta_{w \mid c}\}$ constrained to the simplex, $\theta_{w \mid c} \ge 0$ and $\sum_{w \in \mathcal{V}} \theta_{w \mid c} = 1$ for each context $c$.

Training objective

Maximum likelihood over the corpus, subject to the simplex constraints. Writing $C(c,w)$ for the number of times symbol $w$ follows context $c$ in $\mathcal{D}$, the log-likelihood collapses into a sum over contexts:

Introducing a Lagrange multiplier $\lambda_c$ per context and setting the derivative to zero gives $C(c,w)/\theta_{w\mid c} = \lambda_c$, and the constraint fixes $\lambda_c = C(c) := \sum_w C(c,w)$. The maximum-likelihood solution is a closed-form ratio of counts:

Train:
  1. C <- empty hash map                    # sparse count table
  2. for each sequence x in D:
  3.     for t = 1..T:
  4.         C[(x[t-n+1:t-1], x[t])] += 1
  5. theta[w|c] <- (C[c,w] + alpha) / (C[c] + alpha*|V|)   # add-alpha smoothing

Sample:
  1. c <- start context (BOS symbols)
  2. for t = 1..T:
  3.     x[t] ~ Categorical(theta[ . | c])
  4.     c    <- shift(c, x[t])
  5. return x

Inference / generation

Ancestral sampling: draw $x_1 \sim \theta_{\cdot \mid c_1}$, append it to the context, draw $x_2$, and so on. Every draw is an exact sample from $p_\theta$ — a property autoregressive models keep forever and that no other family (except flows) enjoys. Density evaluation is a single pass of table lookups.

Why it breaks

The table has $O(|\mathcal{V}|^n)$ entries, so it grows exponentially in context length while the data does not: most contexts are seen zero times, and MLE assigns them probability zero, making the test log-likelihood $-\infty$. The entire smoothing literature (add-$\alpha$, Good–Turing, Kneser–Ney, backoff and interpolation) exists to patch this. But smoothing is a symptom, not a cure. The real disease is that the model has no notion of similarity between contexts: "the cat sat on the" and "a cat sat upon a" are unrelated table keys.

3.2 Mixture models and EM: where the ELBO comes from

Setup. Introduce a discrete latent $z \in \{1,\dots,K\}$ indicating which of $K$ components generated $x$. The Gaussian mixture model is

The log-likelihood $\sum_i \log \sum_k \pi_k \mathcal{N}(x^{(i)}; \mu_k,\Sigma_k)$ has a logarithm outside a sum, so it has no closed-form maximizer. The fix — and this is the most consequential derivation in the whole post — is to introduce an arbitrary distribution $q(z)$ and apply Jensen's inequality:

$\mathcal{L}(q,\theta)$ is the evidence lower bound (ELBO). The gap is exactly a KL divergence:

Expectation–Maximization is coordinate ascent on $\mathcal{L}$. The E-step maximizes over $q$ (the gap is zero when $q = p_\theta(z\mid x)$, which for a mixture is a tractable $K$-way softmax); the M-step maximizes over $\theta$ with $q$ held fixed. Because each step increases a lower bound that touches the true objective after the E-step, the likelihood is monotonically non-decreasing.

Init: pi, mu, Sigma  (e.g. k-means)
repeat until the log-likelihood stops improving:

  E-step   (exact posterior; responsibilities)
      r[i,k] <- pi[k] * N(x_i; mu[k], Sigma[k])
                / sum_j pi[j] * N(x_i; mu[j], Sigma[j])

  M-step   (closed-form weighted moments)
      N_k      <- sum_i r[i,k]
      pi[k]    <- N_k / N
      mu[k]    <- (1/N_k) * sum_i r[i,k] * x_i
      Sigma[k] <- (1/N_k) * sum_i r[i,k] * (x_i - mu[k])(x_i - mu[k])^T

Inference / generation

Ancestral again, now through the latent: sample $k \sim \mathrm{Categorical}(\pi)$, then $x \sim \mathcal{N}(\mu_k, \Sigma_k)$. Density evaluation is an exact $K$-term sum.

Relatives worth naming

• Hidden Markov models put the latent in a chain, $p_\theta(x,z) = p(z_1)\prod_t p(z_t \mid z_{t-1}) p(x_t \mid z_t)$. The E-step becomes the forward–backward algorithm; dynamic programming replaces the $K$-way softmax. This is the ancestor of linear state-space sequence models.
• Factor analysis and probabilistic PCA take $z \in \mathbb{R}^d$ continuous with a linear decoder, $x = Wz + \mu + \varepsilon$. This is exactly a VAE with a linear decoder — and because everything is Gaussian, the posterior is available in closed form, so no encoder network is needed.
• Latent Dirichlet allocation made variational inference (rather than exact EM) standard practice for models whose posterior is not tractable.

3.3 Markov random fields: probability as energy

The third classical template comes from statistical physics. Any strictly positive distribution can be written as a Gibbs distribution

where $E_\theta$ is an energy: low energy means high probability. The Ising model takes $x \in \{-1,+1\}^D$ on a lattice with $E_\theta(x) = -\sum_{(i,j)\in \mathcal{E}} W_{ij}x_i x_j - \sum_i b_i x_i$, and the Hammersley–Clifford theorem says the Gibbs form and conditional-independence structure are equivalent.

This parameterization is maximally convenient (any neural network can be an energy, with no architectural constraints) and maximally inconvenient (the sum defining $Z(\theta)$ has $2^D$ terms). That tension is obstacle #1 from §2, and resolving it drives the next thirty years.

4. Era II — Early neural generative models (1980s–2012)

4.1 Boltzmann machines: learning an energy

Setup. A Boltzmann machine is an Ising model with learned weights and hidden units. Split the state into visible $v \in \{0,1\}^{D}$ and hidden $h \in \{0,1\}^{M}$, define a joint energy $E_\theta(v,h)$, and marginalize:

where $F_\theta$ is the free energy. Note that this model faces obstacles #1 and #2 simultaneously.

Training objective and its gradient

Maximum likelihood, as always. Differentiating $\log p_\theta(v) = -F_\theta(v) - \log Z(\theta)$ and using $\nabla_\theta \log Z(\theta) = -\mathbb{E}_{\tilde v \sim p_\theta}[\nabla_\theta F_\theta(\tilde v)]$ gives the defining equation of energy-based learning:

Read it as a contrast: push energy down on real data, push energy up on the model's own samples. Training halts when the two phases balance, i.e. when the model's samples are statistically indistinguishable from the data under the energy's sufficient statistics. The negative phase requires sampling from $p_\theta$ — which requires MCMC — and that is the bottleneck.

4.2 Restricted Boltzmann machines and contrastive divergence

Removing the within-layer connections yields the RBM, a bipartite energy

whose conditionals factorize completely:

with $\sigma$ the logistic function. Factorized conditionals mean an entire layer can be resampled in one vectorized operation — block Gibbs sampling. The free energy is also closed-form, $F_\theta(v) = -b^\top v - \sum_j \mathrm{softplus}(c_j + W_{:,j}^\top v)$, and the weight gradient reduces to a difference of correlations, $\nabla_W \log p_\theta(v) = \mathbb{E}_{\text{data}}[v h^\top] - \mathbb{E}_{\text{model}}[v h^\top]$.

Contrastive divergence (CD-$k$) approximates the intractable model expectation by running the Gibbs chain for only $k$ steps — often $k=1$ — initialized at the data rather than at random. It is a biased gradient estimator (it does not follow the gradient of any objective exactly), but it works well enough in practice and made deep unsupervised learning viable in 2006.

repeat:
  1. sample minibatch v0 ~ D
  2. h0 ~ p(h | v0)                          # positive phase, one vectorized pass

  3. v_k, h_k <- v0, h0                      # negative phase: k Gibbs sweeps
     for s = 1..k:
         v_k ~ p(v | h_k)
         h_k ~ p(h | v_k)

  4. grad_W <- v0 @ h0^T  -  v_k @ h_k^T     # <vh> under data  -  under model
     grad_b <- v0 - v_k
     grad_c <- h0 - h_k
  5. theta <- theta + eta * grad             # ascent: we are maximizing log-likelihood

Inference / generation

There is no ancestral shortcut. To sample you must run the Gibbs chain from a random state until it mixes, which for a multimodal energy can take arbitrarily long — the chain gets trapped in a basin and rarely crosses between modes. This mixing problem is the fundamental limitation of energy-based generation, and it is precisely what diffusion solves fifteen years later.

Stacking RBMs greedily gives deep belief networks, which in 2006 were the first practical way to train deep networks; layer-wise generative pretraining followed by discriminative fine-tuning was the standard recipe until better initialization, ReLUs, and more data made it unnecessary.

4.3 Score matching: the trick that removes $Z$

In 2005 Hyvärinen asked a sharp question: what if we never need $Z(\theta)$ at all? Observe that the partition function is a constant in $x$, so it vanishes under the data gradient:

So if we are content to learn the score $s_\theta(x) \approx \nabla_x \log p_{\text{data}}(x)$ rather than the density, normalization simply does not arise. Score matching minimizes the Fisher divergence

which looks unusable because $\nabla_x\log p_{\text{data}}$ is unknown — but integration by parts eliminates it, leaving an objective computable from data alone:

The remaining cost is the trace of the Jacobian, which is $O(D)$ backward passes — prohibitive for images. Two workarounds were found: sliced score matching (Hutchinson's estimator) and, far more importantly, denoising score matching (Vincent, 2011). Perturb the data with Gaussian noise, $\tilde x = x + \sigma \varepsilon$ with $\varepsilon \sim \mathcal{N}(0,I)$, and match the score of the noisy distribution. Since $\nabla_{\tilde x}\log q_\sigma(\tilde x \mid x) = -(\tilde x - x)/\sigma^2$ is available in closed form, the objective becomes a plain regression:

4.4 Neural autoregressive models and the Helmholtz machine

Two other threads from this era matter. First, Bengio's neural probabilistic language model (2003) replaced the $n$-gram lookup table with an embedding plus MLP: $p_\theta(x_t \mid x_{\lt t}) = \mathrm{softmax}\big(g_\theta(\mathbf{e}_{x_{t-n+1}},\dots,\mathbf{e}_{x_{t-1}})\big)$. Contexts now live in $\mathbb{R}^d$ and share statistical strength; the exponential table is gone. NADE (2011) applied the same idea to arbitrary binary vectors with weight sharing across conditionals. The factorization is unchanged from §3.1 — only the function class changed.

Second, the Helmholtz machine (1995) and its wake–sleep algorithm introduced a separate recognition network $q_\phi(z \mid x)$ trained alongside the generative network $p_\theta(x \mid z)$: the wake phase updates $\theta$ using latents inferred by $q_\phi$; the sleep phase updates $\phi$ on the generative model's own dreams. This is amortized inference, eighteen years before the VAE. What was missing was a single objective that both networks optimize and a low-variance way to backpropagate through sampling.

5. Era III — The deep generative era (2013–2019)

5.1 Variational autoencoders: amortized EM with a reparameterization

Setup. A continuous latent $z \in \mathbb{R}^d$ with a fixed prior $p(z) = \mathcal{N}(0,I_d)$, and a neural decoder producing the parameters of the observation model, e.g. $p_\theta(x \mid z) = \mathcal{N}\big(x; \mu_\theta(z), \sigma^2 I\big)$ for images or a product of Bernoullis. Marginalizing $z$ is intractable (obstacle #2), and unlike a mixture the posterior $p_\theta(z\mid x)$ has no closed form because $\mu_\theta$ is nonlinear.

Training objective

Take the ELBO from §3.2 and make $q$ the output of an encoder network: $q_\phi(z \mid x) = \mathcal{N}\big(z; \mu_\phi(x), \operatorname{diag}\sigma^2_\phi(x)\big)$. Rearranging the ELBO into its standard two-term form,

with the bound gap equal to $D_{\mathrm{KL}}(q_\phi(z\mid x)\Vert p_\theta(z\mid x))$. For diagonal Gaussians the KL term is analytic, $\tfrac12\sum_{j=1}^d\big(\mu_{\phi,j}^2 + \sigma_{\phi,j}^2 - \log\sigma_{\phi,j}^2 - 1\big)$, so only the reconstruction term needs sampling.

The reparameterization trick

The obstacle to gradient descent is that $z$ is sampled from a distribution that depends on $\phi$, so $\nabla_\phi \mathbb{E}_{q_\phi}[\,\cdot\,]$ cannot pass through the sampler. Rewrite the sample as a deterministic function of a parameter-free noise variable:

Now the randomness is in $\varepsilon$, which does not depend on $\phi$, and the gradient flows straight through $\mu_\phi$ and $\sigma_\phi$. The resulting estimator has dramatically lower variance than the score-function (REINFORCE) alternative, and that variance reduction — not the ELBO, which was thirty years old — is what made the VAE work.

Train:
  repeat:
    1. x        ~ minibatch from D
    2. mu, logvar <- Encoder_phi(x)
    3. eps      ~ N(0, I)
    4. z        <- mu + exp(0.5*logvar) * eps            # reparameterization
    5. L_rec    <- log p_theta(x | z)                    # e.g. -||x - Dec(z)||^2 or -BCE
    6. L_kl     <- 0.5 * sum(mu^2 + exp(logvar) - logvar - 1)
    7. loss     <- -(L_rec - beta * L_kl)                # beta = 1 is the exact ELBO
    8. (theta, phi) <- Adam step on loss

Sample:
  1. z ~ N(0, I)
  2. x <- Decoder_theta(z)        # or sample from p_theta(x|z); one forward pass

Inference / generation

Sampling is a single decoder pass from a prior draw — fast and simple. Density is only available as the lower bound $\mathcal{L}$ (tightened, if needed, by importance weighting as in IWAE). The encoder also gives you a useful representation for free, which is why VAEs remain ubiquitous as components even where they lost as standalone image generators.

Failure modes, and what they teach

• Blurriness. With a Gaussian likelihood the reconstruction term is a squared error, whose optimum is a conditional mean. Averaging over plausible completions produces a blur. Combined with the mode-covering bias of forward KL from §2, VAE samples are smooth rather than sharp.
• Posterior collapse. With a powerful autoregressive decoder, the model can drive $q_\phi(z\mid x) \to p(z)$, zeroing the KL term and ignoring $z$ entirely. The latent carries no information.

5.2 Generative adversarial networks: abandoning likelihood

Setup. Define the model implicitly. A generator $G_\theta : \mathcal{Z}\to\mathcal{X}$ pushes a simple prior forward, $x = G_\theta(z)$ with $z\sim p(z)$, inducing a distribution $p_\theta$ that you can sample from but whose density you cannot evaluate — there is no $\log p_\theta(x)$ anywhere in the method. Since likelihood is off the table, we need a different way to compare $p_\theta$ to $p_{\text{data}}$: train a critic to tell them apart.

Training objective

With a discriminator $D_\psi : \mathcal{X}\to(0,1)$, the two-player minimax game is

For fixed $G_\theta$, the inner maximization has a closed-form solution obtained pointwise:

and substituting it back shows what the generator is really minimizing:

So a GAN performs distribution matching under Jensen–Shannon divergence rather than forward KL. JS is symmetric and bounded, and in practice the generator behaves in a mode-seeking way: it is not infinitely punished for ignoring a mode, only for producing samples the critic can flag as fake. That is the sharpness/diversity trade-off in one line, and it is the exact mirror image of the VAE's.

Two practical amendments matter. The non-saturating loss replaces the generator's $\log(1-D)$ with $-\log D$, because early in training $D$ rejects everything and the original loss has vanishing gradient. And Wasserstein GANs replace JS with the earth-mover distance via the Kantorovich–Rubinstein dual, $W_1(p_{\text{data}},p_\theta) = \sup_{\Vert f\Vert_L\le 1}\mathbb{E}_{p_{\text{data}}}[f] - \mathbb{E}_{p_\theta}[f]$, enforced with weight clipping or a gradient penalty. Wasserstein distance gives useful gradients even when the two distributions have disjoint support — the common case for data on a low-dimensional manifold.

repeat:
  # --- critic: n_critic steps (n_critic = 1 for vanilla GAN, ~5 for WGAN-GP)
  for i = 1..n_critic:
      x  ~ minibatch from D
      z  ~ N(0, I)
      L_D <- -( log D_psi(x) + log(1 - D_psi(G_theta(z))) )
      psi <- Adam step on L_D

  # --- generator: one step
  z  ~ N(0, I)
  L_G <- -log D_psi(G_theta(z))          # non-saturating; NOT +log(1 - D)
  theta <- Adam step on L_G

Inference / generation

One forward pass: $z\sim\mathcal{N}(0,I)$, $x = G_\theta(z)$. This is the fastest sampler of any family, and by 2019 StyleGAN had made it the best-looking one too. The price is that there is no likelihood, no encoder (unless you add one), and training is a non-stationary game rather than an optimization — hence mode collapse (the generator maps many $z$ to one $x$), oscillation, and sensitivity to almost every hyperparameter.

5.3 Normalizing flows: exact likelihood by construction

Setup. Let $f_\theta : \mathcal{Z}\to\mathcal{X}$ be a diffeomorphism — invertible and differentiable, with $\dim\mathcal{Z} = \dim\mathcal{X} = D$. Setting $x = f_\theta(z)$ with $z\sim p_Z = \mathcal{N}(0,I)$, the change-of-variables formula gives the density exactly:

No bound, no partition function, no latent integral. The catch is the $\log\det$ of a $D\times D$ Jacobian, which naively costs $O(D^3)$; flows are therefore an exercise in designing expressive maps whose Jacobian is triangular. The affine coupling layer of NICE/RealNVP splits $x = (x_a, x_b)$ and transforms only half:

which is trivially invertible and has $\log|\det J| = \sum_j s_\theta(x_a)_j$ regardless of how complicated $s_\theta$ and $t_\theta$ are. Alternate which half is transformed, stack many layers, and add $1\times1$ convolutions (Glow) to mix channels.

Training and sampling

Training is direct maximum likelihood: minimize $-\log p_\theta(x)$ averaged over a minibatch, with both terms above computed in a single forward pass. Sampling is $z\sim p_Z$ followed by $x = f_\theta(z)$. Flows are the only family that gives exact density and exact one-pass sampling. What they give up is architecture: invertibility and dimension preservation are severe constraints, and in practice flows need many more parameters than a VAE or diffusion model of comparable quality.

Continuous normalizing flows

Take the number of layers to infinity and the residual step size to zero, and the stack of coupling layers becomes an ODE:

with the density evolving according to the instantaneous change-of-variables (continuity) equation

Now $v_\theta$ can be any network: invertibility comes free from the ODE, and the $\log\det$ becomes a divergence (estimable with Hutchinson's trick). Beautiful — but training by maximum likelihood requires integrating the ODE and backpropagating through the solver on every step, which is slow and stiff. FFJORD works; it just costs too much.

5.4 Deep autoregressive models and the Transformer

Setup. Back to the exact chain rule of §3.1, now with a deep conditional:

The design question is only how to compute the context representation $h_t$. RNNs and LSTMs carry a recurrent state; PixelRNN/PixelCNN impose a raster-scan order on images with masked convolutions; WaveNet uses dilated causal convolutions for audio; the Transformer uses causal self-attention, in which $h_t$ attends over all of $x_{\lt t}$ with a learned, content-based weighting.

Training objective

Plain maximum likelihood, which for discrete $x_t$ is token-level cross-entropy:

Two properties make this the most scalable objective in machine learning. First, the loss is exact — no bound, no adversary, no MCMC, so the training curve means something and improvements are monotone in a real quantity (perplexity, equivalently bits per token, equivalently compression). Second, with causal masking all $T$ conditionals are computed in a single parallel forward pass under teacher forcing: position $t$ predicts $x_t$ while conditioning on the ground-truth prefix. That is the property that turned likelihood training into a hardware-utilization problem, and it is the real reason Transformers scaled.

Train (one parallel pass over the whole sequence):
  1. x = (x_1..x_T)  ~ minibatch
  2. H     <- Transformer_theta(embed(x), causal_mask)    # H[t] sees only x_<t
  3. logits<- H @ W_o^T
  4. loss  <- mean_t CrossEntropy(logits[t], x[t])        # teacher forcing
  5. theta <- AdamW step on loss

Sample (sequential; T passes, KV-cached):
  1. context <- prompt ; cache <- {}
  2. for t = 1..T:
  3.     logits, cache <- Transformer_theta(last_token, cache)
  4.     logits <- logits / tau                           # temperature
  5.     logits <- top_k / top_p filter(logits)
  6.     x[t] ~ Categorical(softmax(logits))
  7.     if x[t] == EOS: break
  8. return x

Inference / generation

Sampling is inherently sequential: $T$ forward passes, each depending on the last. Caching keys and values makes each step cheap, but latency remains $O(T)$ — the exact mirror of training's parallelism. Temperature $\tau$, top-$k$, and nucleus (top-$p$) truncation are post-hoc corrections to the fact that a mode-covering, forward-KL-trained model puts real probability mass on the long tail of bad continuations.

The other structural weakness is the fixed factorization order. A raster-scan image model has no reason to believe pixel $(0,0)$ should be decided first, and errors made early are never revised. This is precisely the weakness diffusion does not have: diffusion decides everything at once, coarse-to-fine, and revisits its own output at every step. It is also why autoregression stayed dominant for text (which really is sequential and discrete) while losing images to diffusion.

6. Era IV — The diffusion era (2020–2023)

The idea is disarmingly simple. Learning to map noise to data in one jump is hard; learning to remove a little noise is easy. So define a fixed process that gradually destroys the data, learn to reverse one small step of it, and generate by chaining many small reversals. Diffusion converts one intractable problem into a thousand tractable regressions.

6.1 The forward process

Setup. Let $x_0 \sim p_{\text{data}}$ and define a Markov chain of increasing corruption over $t = 1,\dots,T$ with a fixed variance schedule $\beta_t \in (0,1)$:

Nothing here is learned. Because Gaussians compose, the $t$-step marginal is available in closed form — the property the whole method depends on. With $\alpha_t = 1-\beta_t$ and $\bar\alpha_t = \prod_{s\le t}\alpha_s$,

So we can jump to any noise level in one line of code, with no simulation. The schedule is chosen so $\bar\alpha_T \approx 0$, making $q(x_T)\approx\mathcal{N}(0,I)$: the process destroys all information and lands on a distribution we can sample from trivially.

6.2 Training objective: three derivations, one loss

The reverse of a Gaussian diffusion with small steps is itself approximately Gaussian, so we parameterize $p_\theta(x_{t-1}\mid x_t) = \mathcal{N}\big(x_{t-1};\mu_\theta(x_t,t),\sigma_t^2 I\big)$ and fit it. Three different starting points converge on the same objective, which is why the method feels inevitable.

• As a hierarchical VAE. Treat $x_{1:T}$ as latents with a fixed encoder $q$. The ELBO of §3.2 becomes $\mathbb{E}_q\big[\log p(x_T) + \sum_t \log \frac{p_\theta(x_{t-1}\mid x_t)}{q(x_t\mid x_{t-1})}\big]$, which rearranges into a sum of per-step KL divergences $D_{\mathrm{KL}}\big(q(x_{t-1}\mid x_t,x_0)\,\Vert\,p_\theta(x_{t-1}\mid x_t)\big)$. Both sides are Gaussian, so each KL is a squared difference of means.
• As denoising score matching. Since $q(x_t\mid x_0)$ is Gaussian, its conditional score is $\nabla_{x_t}\log q(x_t\mid x_0) = -\varepsilon/\sqrt{1-\bar\alpha_t}$, and §4.3's DSM objective applies directly at every noise level simultaneously.
• As Tweedie's formula. For $x_t = \sqrt{\bar\alpha_t}x_0 + \sqrt{1-\bar\alpha_t}\varepsilon$, the posterior mean satisfies $\mathbb{E}[x_0\mid x_t] = \big(x_t + (1-\bar\alpha_t)\nabla_{x_t}\log p_t(x_t)\big)/\sqrt{\bar\alpha_t}$. Predicting the noise, predicting the clean data, and predicting the score are three reparameterizations of one quantity.

Ho et al. showed that dropping the per-step weights from the ELBO gives a simpler and empirically better objective — the loss that launched the era:

Read what this is: a mean-squared-error regression. One network, one loss, no adversary, no partition function, no encoder to learn, no ODE to integrate. The score is recovered whenever you need it via $s_\theta(x_t,t) = -\varepsilon_\theta(x_t,t)/\sqrt{1-\bar\alpha_t}$, and the random choice of $t$ each step is what lets a single network cover all noise levels.

Train:
  repeat:
    1. x0  ~ D
    2. t   ~ Uniform{1..T}
    3. eps ~ N(0, I)
    4. xt  <- sqrt(abar[t]) * x0 + sqrt(1 - abar[t]) * eps      # closed-form jump
    5. loss<- || eps - eps_theta(xt, t) ||^2
    6. theta <- Adam step on loss    (keep an EMA copy of theta for sampling)

Sample (ancestral, T steps):
  1. x_T ~ N(0, I)
  2. for t = T..1:
  3.     z <- N(0, I) if t > 1 else 0
  4.     x_{t-1} <- (1/sqrt(alpha[t])) * ( x_t - (beta[t]/sqrt(1-abar[t])) * eps_theta(x_t, t) )
                    + sigma[t] * z
  5. return x_0

6.3 The SDE view and the probability-flow ODE

Song et al. showed that the discrete chain is a discretization of a stochastic differential equation, $dx = f(x,t)\,dt + g(t)\,dw$, and that this SDE has an exactly reversible counterpart driven by the score:

More useful still, there is a deterministic ODE with the same time-marginals $p_t$ — the probability-flow ODE:

Three consequences follow immediately, and they close several loops at once.

1. A trained diffusion model is a continuous normalizing flow (§5.3): integrate the PF-ODE and you get exact likelihoods via the divergence formula. Diffusion trains a CNF without ever simulating one.
2. Sampling is now an ODE-solving problem, so decades of numerical integration apply. DDIM, DPM-Solver, and Heun-style solvers cut sampling from 1000 steps to 10–50 with no retraining. This is where nearly all early inference-speed progress came from.
3. Deterministic sampling makes the map $x_T \mapsto x_0$ invertible, giving latent interpolation and real-image editing — the properties GANs had and VAEs promised.

6.4 Guidance and latent diffusion

Conditioning. To sample from $p(x\mid c)$ for a prompt $c$, Bayes' rule on scores gives $\nabla_x \log p(x\mid c) = \nabla_x\log p(x) + \nabla_x\log p(c\mid x)$, which classifier guidance implements with a separately trained classifier. Classifier-free guidance removes the classifier by training one network with the condition randomly dropped (say 10% of the time, replaced by a null token $\varnothing$) and then extrapolating at sampling time:

This sharpens the conditional distribution — effectively sampling from something like $p(x)p(c\mid x)^w$ — and trades diversity for prompt fidelity. It is the single most important practical trick in text-to-image generation, and it is worth noticing that it is a deliberate departure from maximum likelihood: we correct the mode-covering bias of §2 at inference time rather than at training time.

Latent diffusion. Running diffusion on $512\times512\times3$ pixels is wasteful, because most of that capacity models imperceptible high-frequency detail. Latent diffusion first trains an autoencoder (with perceptual and adversarial losses — §5.2 returning as a component) to compress images by roughly $8\times$ per side, then runs the entire diffusion process in that latent space. A one- to two-order-of-magnitude compute reduction is what made Stable Diffusion runnable on a consumer GPU, and the pattern — learn a compact space, then model the distribution inside it — is now standard for video and audio too.

7. Era V — Flow matching and hybrids (2022–)

7.1 Flow matching: regress the velocity, skip the simulation

Setup. Return to the CNF of §5.3, whose model is a velocity field $v_\theta(x,t)$. We want the ODE $dx_t/dt = v_\theta(x_t,t)$ to transport a simple $p_0$ into $p_1 = p_{\text{data}}$. The ideal objective would be flow matching,

where $u_t$ is a vector field generating the desired path $p_t$ — but $u_t$ is unknown. The key theorem is that we may replace it with a conditional field defined per data point. Choose an interpolation between a noise sample $x_0\sim p_0$ and a data sample $x_1\sim p_{\text{data}}$,

whose conditional velocity is immediate by differentiation, $u_t(x_t \mid x_0,x_1) = \dot\alpha_t x_1 + \dot\sigma_t x_0$. Then the conditional flow matching loss

has the same gradients in $\theta$ as $\mathcal{L}_{\text{FM}}$. The reason is the standard conditional-expectation argument: the minimizer of a squared error against a random target is the conditional mean, and $\mathbb{E}[u_t(x_t\mid x_0,x_1)\mid x_t] = u_t(x_t)$. The network never sees the intractable marginal field, but converges to it anyway.

Rectified flow is the simplest instance: take the straight line $\alpha_t = t$, $\sigma_t = 1-t$, so $x_t = (1-t)x_0 + t x_1$ and the regression target is just the constant displacement $x_1 - x_0$:

Train:
  repeat:
    1. x1 ~ D                              # data
    2. x0 ~ N(0, I)                        # noise (or any tractable source dist.)
    3. t  ~ Uniform[0,1]                   # (logit-normal schedules work better in practice)
    4. xt <- (1 - t) * x0 + t * x1
    5. loss <- || v_theta(xt, t) - (x1 - x0) ||^2
    6. theta <- Adam step on loss

Sample (Euler, N steps; any ODE solver works):
  1. x <- N(0, I) ; h <- 1/N
  2. for i = 0..N-1:
  3.     t <- i * h
  4.     x <- x + h * v_theta(x, t)
  5. return x

Set that beside Algorithm 7 and the punchline is clear: flow matching and diffusion are the same algorithm with a different interpolation schedule and a different regression target parameterization. Given a Gaussian source, one can convert between $\varepsilon$-prediction, $x_0$-prediction, $v$-prediction, and score. What flow matching adds is not a new capability but a cleaner formulation, and three real benefits: the path is a design choice rather than an inherited SDE, straight paths are cheaper to integrate (fewer steps for the same quality), and $p_0$ need not be Gaussian — you can flow between two arbitrary distributions, which is what makes image-to-image and cross-modal transport natural. The current generation of large image and video models is trained this way.

7.2 Few-step generation: distillation

The remaining weakness of transport-based models is that sampling needs many network evaluations. Three families of fixes, all of which recycle earlier ideas:

• Consistency models train $f_\theta(x_t,t)$ to map any point on a PF-ODE trajectory to the same endpoint, enforcing self-consistency $f_\theta(x_t,t) = f_\theta(x_{t'},t')$ along a trajectory with the boundary condition $f_\theta(x,0) = x$. One evaluation then jumps from noise to data.
• Reflow retrains a rectified flow on its own (noise, sample) pairs. Because the coupling is now deterministic rather than independent, the induced paths straighten, and a straight path is exactly integrated by a single Euler step.
• Adversarial and distribution-matching distillation add a discriminator or a reverse-KL score-difference term to a few-step student. This is §5.2 returning as the finishing tool: the GAN objective is excellent at making a fixed, small number of steps produce sharp samples, and the diffusion teacher supplies the mode coverage the GAN alone would lose.

7.3 Discrete data: diffusion for text

Gaussian noise is meaningless on tokens, so discrete corruption processes replace it. The approach that works is absorbing-state (masked) diffusion: the forward process replaces each token independently with [MASK] with probability increasing in $t$, and the model predicts the original tokens at masked positions. The objective is a weighted cross-entropy over masked positions,

which is a principled variational bound and, notably, is BERT-style masked language modeling generalized over all masking rates. Generation unmasks tokens progressively in any order, trading the strict left-to-right constraint of §5.4 for parallel decoding and the ability to revise. Whether this displaces autoregression for language is, as of now, unsettled: the AR objective still extracts more signal per token, but masked diffusion decodes many tokens per pass.

7.4 Hybrids: what people actually deploy

No production system is a pure member of one family. The dominant patterns compose them:

• Autoencoder + diffusion (latent diffusion): a VAE/VQGAN compresses, trained with adversarial and perceptual losses; a flow-matching or diffusion model handles the distribution in latent space. Three families in one system.
• Tokenizer + autoregression (VQ-VAE → VQGAN → discrete AR): a vector-quantized autoencoder turns images into token sequences, then a Transformer models them with cross-entropy. This is how the AR lineage competes on images, and how unified multimodal models put pixels and text in one sequence.
• Autoregressive backbone + diffusion head: a Transformer produces context vectors and a small diffusion or flow head models the continuous distribution of the next element, avoiding quantization loss while keeping sequential structure.
• Diffusion teacher + adversarial student: as in §7.2, for real-time generation.

8. The unifying picture

Strip away the architectures and every family is a choice of what the network outputs and how that output is scored against data.

The connections, stated plainly

1. Counting $\to$ Transformers. Identical factorization $\prod_t p(x_t\mid x_{\lt t})$, identical cross-entropy objective, identical ancestral sampler. Only the estimator of the conditional changed — from a table to a function that generalizes across contexts.
2. EM $\to$ VAE. The same ELBO. EM's E-step computes the exact posterior per datapoint; the VAE amortizes an approximate posterior into a network and makes it differentiable via reparameterization.
3. Energy $\to$ score $\to$ diffusion. Taking $\nabla_x$ of $\log p_\theta = -E_\theta - \log Z$ annihilates $Z$, converting an intractable normalization problem into a regression. Adding a noise schedule turns single-level denoising score matching into diffusion, and turns slow-mixing MCMC into a fixed-length reverse process that reaches all modes because it starts from pure noise.
4. VAE $\to$ diffusion. A diffusion model is a hierarchical VAE with $T$ latents, a fixed non-learned encoder, and latents of the same dimension as the data. Its loss is that ELBO with the per-step weights dropped.
5. CNF $\to$ flow matching. Same model object (a velocity field), same sampler (an ODE solve). Flow matching removes ODE simulation from training by supplying a closed-form conditional regression target.
6. Diffusion $\leftrightarrow$ flow matching. Two parameterizations of the same transport. The probability-flow ODE of a diffusion model is a CNF; a Gaussian-source flow-matching model is a diffusion model with a different schedule. Choosing straight paths is what buys fewer sampling steps.
7. GAN $\to$ everything's last mile. The adversarial objective survives inside the autoencoders that latent diffusion depends on, and inside the distillation that makes few-step samplers sharp.

Three things worth remembering

First, the objective determines the failure mode. Forward KL is mode-covering and gives blur and a heavy tail of implausible samples; JS/Wasserstein under an adversary is mode-seeking and gives sharpness with dropped modes; the denoising regression that underlies diffusion and flow matching is likelihood-flavored but reweighted to emphasize the noise levels that matter perceptually, which is why it gets both. Guidance and truncated sampling are inference-time corrections applied on top — deliberate, useful departures from the distribution the model was fit to.

Second, progress came from replacing intractable quantities with regressions. The partition function, the posterior, the marginal velocity field: each was at some point the blocker, and each was removed by finding an equivalent least-squares problem whose target is computable from a single data point plus a single noise draw. That is the shape of nearly every breakthrough above.

Third, the split that remains is data type, not method. Discrete and genuinely sequential data (language, code) still belongs to autoregression, where the exact chain rule and fully parallel teacher-forced training are hard to beat. Continuous, high-dimensional, simultaneously-structured data (images, audio, video, molecules, protein structure) belongs to transport — diffusion and flow matching — where coarse-to-fine refinement beats committing to a scan order. The interesting frontier is where the two meet: discrete diffusion pulling parallel decoding into language, and autoregressive backbones with continuous heads pulling sequential structure into images and video.

Seventy years, one question, and a short list of tricks used over and over: factorize with the chain rule, bound with Jensen, differentiate away the normalizer, reparameterize the sampler, and regress against a conditional target you can compute in closed form.