Deep generative models: autoregressive, latent variable, flow, and score based

Why this subject matters now

Ten years ago generative modeling was a subfield with a credibility problem. Samples were blurry thumbnails and likelihoods were numbers nobody outside the community could interpret. Today the generative model is the product. Language models are autoregressive generative models of token sequences. Image, video, and audio systems are diffusion or flow-matching models, which this page will derive as the modern resolution of a sampling problem energy-based models posed decades earlier. Multimodal systems compress images into discrete codes with a VQ-VAE-style quantizer and then model the codes autoregressively. The practical stack of 2026 is built almost entirely out of the four or five ideas this page covers.

What a practitioner is expected to know has changed accordingly. It is no longer enough to know one family. The families are trade-offs against each other along a small set of axes, and interviews and design reviews probe exactly those axes, asking why an autoregressive model gives you an exact likelihood but slow sampling, why a VAE gives you fast sampling but only a bound on the likelihood, why a flow gives you both but pays in architectural constraints, why a GAN gives you neither but for years gave the sharpest images, and why diffusion took over continuous media by splitting one hard generation problem into many easy denoising problems. The families have also visibly converged. Flow matching, diffusion, consistency models, and normalizing flows are now understood as different ways to learn a transport map from noise to data, and the discrete side (autoregressive models over learned token vocabularies) is the backbone of every frontier lab's unified multimodal effort. A practitioner who understands the unifying structure can read any new paper in the area as a point in a known design space rather than a new invention.

This page is the map. Two neighboring pages go deep on the two families that earned their own treatment, diffusion and large vision models and adversarial generative modeling. Here they appear as members of the taxonomy, derived far enough to see exactly where they sit and why, with pointers forward for the rest.

The unifying question

Fix a dataset \( \D = \{x^{(1)}, \dots, x^{(n)}\} \) drawn i.i.d. from an unknown distribution \( p_{\text{data}} \) over a space \( \mathcal{X} \), which might be \( \R^{3072} \) for CIFAR images or a discrete space of token sequences. A generative model is a family \( \{p_\theta\} \) of distributions together with a training procedure that picks \( \theta \) so that \( p_\theta \approx p_{\text{data}} \). Everything else is a design decision about what you can afford to compute, and the entire field is organized by four axes.

Explicit versus implicit density. An explicit model writes down a formula for \( p_\theta(x) \), or at least for something provably related to it. An implicit model only defines a sampler. Draw \( z \sim p(z) \), push it through a network, and call the output distribution \( p_\theta \). GANs are the canonical implicit model, so you can draw from \( p_\theta \) all day but cannot evaluate it at a point.

Tractable versus approximate likelihood. Among explicit models, autoregressive models and normalizing flows give you \( \log p_\theta(x) \) exactly, in one or a few network passes. VAEs give a lower bound (the ELBO). Energy-based models give an unnormalized density, \( \log p_\theta(x) \) up to a constant \( \log Z_\theta \) that is itself intractable. Diffusion models give a bound (or an exact value through an ODE solve, at real cost). This axis decides whether you can do model comparison, compression, and out-of-distribution scoring with the model, not just sampling.

Sampling cost. A flow or a GAN samples in one forward pass. A VAE samples in one decoder pass. An autoregressive model needs one pass per dimension, which is why language model inference is an engineering field of its own. A diffusion model needs tens to thousands of passes, and an energy-based model needs an MCMC chain of unknown length. The measured gap is not subtle. The MADE model trained for this page evaluates the exact likelihood of a batch of images in a 0.43 ms forward pass but needs 172 ms of sequential passes, 404 times longer, to sample the same batch.

Latent structure. Some models introduce a latent variable \( z \) meant to capture the underlying factors of the data, giving you a representation, an interpolation space, and a knob for controllable generation. Others (autoregressive models, EBMs) model \( x \) directly and offer no such handle. Whether the latent is continuous (VAE, flow) or discrete (VQ-VAE) turns out to matter enormously for what you can build on top.

The rest of this page fills in this taxonomy.

One more observation frames everything that follows. The reason this is hard at all is that \( \mathcal{X} \) is enormous and \( p_{\text{data}} \) is concentrated on a vanishingly small, complicated subset of it, the manifold hypothesis. A 28×28 binary image space has \( 2^{784} \approx 10^{236} \) points, and the set that looks like a digit is a very thin sliver. Every family below is a different strategy for spending network capacity on that sliver. Slice the joint into conditionals (autoregressive), squeeze it through a bottleneck (VAE), warp a simple density onto it (flow), carve it out with an energy landscape (EBM), or learn the vector field that flows toward it (score-based).

Maximum likelihood and its geometry

Maximizing likelihood is minimizing forward KL

Most of the families train by maximum likelihood, so it pays to know exactly what that objective optimizes at the population level. The average log-likelihood over the dataset is a Monte Carlo estimate of \( \E_{x \sim p_{\text{data}}}[\log p_\theta(x)] \). Expand the forward KL divergence from data to model,

The first term is the negative entropy of the data distribution and does not depend on \( \theta \). Therefore

and maximum likelihood is exactly forward KL minimization, up to the finite-sample approximation of the expectation. This identity also calibrates what a likelihood number means. The best achievable average log-likelihood is \( -H(p_{\text{data}}) \), and the shortfall from it is precisely the KL in nats. When the VAE trained below reports a test ELBO of −97.08 nats per image, that number is (a bound on) the model's average code length for a binarized MNIST digit. The entropy of the data source is the unreachable floor.

Forward versus reverse KL, precisely

The KL divergence is asymmetric, and the asymmetry is not a technicality. It is the most reused piece of intuition on this page. Compare the two directions as functionals of the model \( q \) fit to a fixed target \( p \),

Forward KL is zero-avoiding and mode-covering. In \( \KL(p \| q) \) the integrand is weighted by \( p \). At any region where \( p(x) > 0 \) but \( q(x) \to 0 \), the ratio \( p/q \to \infty \) and the divergence blows up. The optimal \( q \) in a limited family must therefore put mass everywhere the data does, even at the price of putting mass where the data has none, in the valleys between modes. Concretely, for an exponential family \( q \), minimizing forward KL is moment matching. Setting the gradient with respect to the natural parameters to zero forces the expected sufficient statistics under \( q \) to equal those under \( p \). A single Gaussian fit by maximum likelihood to a bimodal density lands its mean between the modes and inflates its variance to cover both, generating many samples that look like neither mode. Worked Problem 1 makes this quantitative. The maximum-likelihood Gaussian fit to a well-separated two-mode mixture puts about 700 times more density at the empty midpoint than the data does.

Reverse KL is zero-forcing and mode-seeking. In \( \KL(q \| p) \) the integrand is weighted by \( q \). The model is only charged where it puts mass. If \( q(x) > 0 \) where \( p(x) \approx 0 \), the ratio \( q/p \) explodes, so the optimizer forces \( q \to 0 \) wherever \( p \) is small. But regions where \( p > 0 \) and \( q = 0 \) cost nothing. The optimal constrained \( q \) therefore locks onto one mode and ignores the others entirely. A reverse-KL Gaussian fit to the same mixture sits tightly on one component, matching its local shape and pretending the other component does not exist.

This dichotomy recurs through the whole page. Maximum-likelihood models (autoregressive, flows, VAEs, diffusion) are mode-covering. They never miss a mode, and their failure mode is over-dispersion, mass smeared into implausible regions, which the eye reads as blur or incoherence. The variational posterior in a VAE is trained with a reverse KL, \( \KL(q_\phi \| p) \), which is why it systematically underestimates posterior variance. GAN training approximates a symmetric Jensen-Shannon-type divergence but in practice behaves mode-seekingly. The generator's failure mode is dropping modes, which the eye cannot detect in any single sample, and which is exactly why GAN samples looked better than their coverage was. No divergence direction is right. They fail differently, and evaluation (a later section) has to measure both failure directions separately.

Autoregressive models

The chain rule factorization

The chain rule of probability holds for any distribution and any ordering of the \( D \) dimensions,

This is not an approximation but an identity. The modeling assumption enters only in how each conditional is parameterized, typically by a single network that maps a prefix to the parameters of a simple one-dimensional (or one-token) distribution, such as a softmax over a vocabulary, a Bernoulli probability per pixel, or a mixture of logistics per subpixel. The payoff is the cleanest likelihood story in the field. The exact log-likelihood is a sum of \( D \) one-dimensional log-probabilities, all computable in a single forward pass if the architecture is arranged so that the network's output at position \( i \) depends only on inputs before \( i \). Training is then plain maximum likelihood with none of the bounds, estimators, or auxiliary networks the other families need.

The cost is the ordering and the sampling loop. The ordering is arbitrary for images (raster order is a convention, not a truth) yet baked into the architecture, and generation must proceed sequentially, sampling \( x_1 \), feeding it back, sampling \( x_2 \), and so on, \( D \) network evaluations that cannot be parallelized across positions. Everything architectural in this family is about enforcing the dependency constraint during parallel training while living with the sequential sampling at generation time.

MADE, masking a fully connected autoencoder

Germain et al. (2015) turn an ordinary MLP autoencoder into an autoregressive density estimator purely with binary masks on the weight matrices. Assign each input dimension its position \( 1, \dots, D \) in the ordering. Assign each hidden unit \( k \) a degree \( m(k) \in \{1, \dots, D-1\} \), interpreted as "this unit may depend on inputs with position \( \le m(k) \)". Then a weight from unit \( j \) (or input \( j \)) to unit \( k \) is allowed iff \( m(k) \ge m(j) \), and the output unit that parameterizes \( p(x_i \mid x_{\lt i}) \) may connect to hidden unit \( k \) iff \( i > m(k) \), the strict inequality being what removes the self-dependence. Composing the masked layers, the output for dimension \( i \) is a function of inputs with positions strictly less than \( i \) only, so the product of the output conditionals is a valid autoregressive model, evaluated for all \( i \) in one pass. This is checkable, not just plausible. For the MADE trained below, the \( 784 \times 784 \) Jacobian of output logits with respect to inputs has maximum absolute value exactly 0.0 on and above the diagonal. The mask argument is a structural zero, not an approximate one.

PixelCNN and the blind spot

For images, van den Oord et al. (2016) enforce the same constraint with masked convolutions. The kernel is zeroed at the center pixel's raster-order successors, so each output position sees only pixels above it and to its left in the same row. Stacking such layers grows the receptive field, but with a defect. The combination of "same row, strictly left" and "any earlier row within the kernel height" composes into a receptive field that permanently excludes a triangular region above and to the right of the current pixel, no matter how many layers are stacked. This is the blind spot, and it means the model silently ignores conditioning information it is entitled to use. The fix in the gated PixelCNN is to split the network into two stacks, a vertical stack that sees all rows strictly above (a convolution masked to past rows, with the growing receptive field of an unmasked column), and a horizontal stack that sees the current row's left prefix and reads from the vertical stack. Each output then depends on the full visible history and nothing else. The same paper replaced ReLUs with a gated unit, \( \tanh(W_f * x) \odot \sigma(W_g * x) \), which measurably improved likelihoods.

Masked 3x3 stack (blind spot)          Two-stack fix (gated PixelCNN)

  rows above:  ███████░░  <- blind      vertical stack:  █████████  all rows above
               █████░░░░     region                      █████████
  current row: ██▓ . . .                horizontal stack: ██▓ . . .  left prefix
                 ^ pixel being predicted                    ^ + reads vertical stack
  ░ = never visible however deep the stack
  █ = visible history   ▓ = current pixel's left neighbor

WaveNet, dilation instead of masking

Audio at 16 kHz needs receptive fields spanning tens of thousands of steps. WaveNet (van den Oord et al., 2016) keeps the causal-convolution idea but dilates. Layer \( \ell \) uses a convolution with holes, reading inputs at offsets \( 0, 2^\ell \) with dilation doubling per layer, \( 1, 2, 4, \dots, 512 \), then repeating the block. The receptive field grows exponentially with depth, \( 1024 \) samples per 10-layer block for kernel size 2, at linear parameter cost. The lesson generalizes beyond audio. The autoregressive constraint is about the dependency graph, and any architecture whose dependency graph is causal works, which is why the family so smoothly absorbed the transformer. A decoder-only transformer with a causal attention mask is exactly this factorization with each conditional computed by attention over the full prefix, removing the receptive-field problem entirely. The mechanics are on the attention page.

The likelihood head, discretized logistic mixtures

Pixels are discrete (256 levels), and the natural head is a 256-way softmax per subpixel, which is what the original PixelCNN used. It works but wastes capacity. The softmax has no notion that level 127 is close to level 128, so nearby intensities are learned as unrelated classes. PixelCNN++ (Salimans et al., 2017) replaces it with a discretized mixture of logistics. The network outputs mixture weights \( \pi_k \), means \( \mu_k \), and log-scales for \( K \approx 10 \) logistic components, and the probability of the discrete level \( x \in \{0, \dots, 255\} \) (rescaled to \( [-1, 1] \)) is the probability mass the continuous mixture assigns to the bin around \( x \),

with the edge bins extended to \( -\infty \) and \( +\infty \) so the 256 masses sum to one, and \( \sigma \) the logistic CDF. The head is ordinal, smooth in its parameters, needs about 100 output channels instead of 768 per pixel, and improved CIFAR-10 density estimation from 3.14 to 2.92 bits per dimension in the paper's reported figures. The same construction, a continuous density integrated over quantization bins, is the honest way to put any continuous likelihood on discrete data, and it returns in the evaluation section.

Teacher forcing and exposure bias

Training maximizes \( \sum_i \log p_\theta(x_i \mid x_{\lt i}) \) with the true prefix \( x_{\lt i} \) from the data, which is teacher forcing. Generation conditions on the model's own sampled prefix. These are different input distributions, and the mismatch is called exposure bias. An early sampling error drags the prefix off the data manifold, where the conditionals were never trained, and errors can compound. Scheduled sampling (Bengio et al., 2015) mixes model samples into training prefixes, though doing so naively no longer optimizes a consistent likelihood. Sequence-level RL objectives are the heavyweight fix. Honest accounting requires saying that for large language models trained on enough data the practical severity of exposure bias is disputed. The models generalize off-manifold better than the small-RNN-era intuition predicted. But the train/sample asymmetry itself is structural to the family, and it is the right first thing to check when an autoregressive model's samples degrade over length.

The sequential cost, by contrast, is not disputed. It is arithmetic. Sampling \( D \) dimensions costs \( D \) forward passes. For the 784-dimensional MADE trained for this page, one forward pass on a batch of 64 takes 0.43 ms on an H100 80GB and full sampling takes 172.3 ms, a 404× gap, and that is with the entire model resident in cache-friendly conditions. For language models the same gap is why KV caching, speculative decoding, and batched serving exist as disciplines.

Latent variable models and the VAE

The marginal likelihood and why it is intractable

A latent variable model posits an unobserved \( z \) behind each \( x \). Draw \( z \sim p(z) \) from a simple prior, usually \( \N(0, I) \), then draw \( x \sim p_\theta(x \mid z) \) from a decoder network that maps \( z \) to the parameters of a simple likelihood. The model of the data is the marginal

The integrand is easy but the integral is not. For a nonlinear decoder there is no closed form, and naive Monte Carlo, \( \frac{1}{K}\sum_k p_\theta(x \mid z_k) \) with \( z_k \sim p(z) \), is hopeless in practice. For a given \( x \), almost all prior samples \( z_k \) decode to something unlike \( x \), so \( p_\theta(x \mid z_k) \) is astronomically small except on a tiny region of \( z \)-space that random sampling essentially never hits. What the estimator needs is samples from the posterior \( p_\theta(z \mid x) = p_\theta(x \mid z) p(z) / p_\theta(x) \), which requires the very quantity being computed. Every training method for this family is a way around this circularity, and the VAE's answer (Kingma & Welling, 2014, and Rezende, Mohamed & Wierstra, 2014, developed concurrently) is to learn the approximate posterior with a second network.

The ELBO, derived twice

Introduce any distribution \( q_\phi(z \mid x) \) with the same support as the posterior. First derivation, Jensen's inequality. Multiply and divide inside the marginal, then use concavity of \( \log \),

This is the evidence lower bound. Splitting \( p_\theta(x, z) = p_\theta(x \mid z) p(z) \) gives its two working parts,

Second derivation, the exact decomposition. Jensen tells you there is a gap but not what it is. Start instead from the definition of the posterior and take the expectation of \( \log p_\theta(x) \) (a constant in \( z \)) under \( q_\phi \),

which is exactly

The gap between the true log-likelihood and the ELBO is the KL from the approximate posterior to the true one, term by term, per data point. Three consequences follow immediately. First, maximizing the ELBO over \( \phi \) with \( \theta \) fixed tightens the bound by driving \( q_\phi \) toward the true posterior, so inference is itself an optimization. Second, maximizing over \( \theta \) both raises the likelihood and can slacken or tighten the bound, so the ELBO is a moving target that pushes the model toward configurations whose posteriors the encoder family can actually represent. Third, the gap term is a reverse KL in the sense of the earlier section. Here \( q \) is the model being optimized, the true posterior is the target, so the learned posterior is mode-seeking and systematically underestimates posterior spread. The importance weighted bound (IWAE) tightens all of this by averaging \( K \) weighted samples inside the log. For the model trained below, the ELBO of −97.08 nats improves to an IWAE-100 estimate of −93.07 nats, directly measuring about 4 nats of slack attributable to the factorized Gaussian posterior.

The reparameterization trick, and why it beats REINFORCE

The reconstruction term is an expectation under \( q_\phi \), and \( \phi \) sits inside the sampling distribution, so the gradient cannot move inside the expectation directly. The general-purpose answer is the score-function (REINFORCE) estimator, derived by differentiating the integral and using \( \nabla_\phi q_\phi = q_\phi \nabla_\phi \log q_\phi \),

It is unbiased and needs only the ability to evaluate \( \log q_\phi \), but its variance is set by the raw magnitude of \( f \). The estimator learns from correlation between \( f(z) \) and the score, jittering \( \phi \) and watching which samples happened to score well. It uses no information about how \( f \) changes with \( z \), even when \( f \) is a differentiable neural network. The reparameterization trick uses exactly that information. Write the sample as a deterministic, differentiable function of \( \phi \) and a parameter-free noise variable,

so the expectation is over \( \varepsilon \) and the gradient passes through,

The pathwise estimator's variance scales with the variance of \( \nabla_z f \) across the noise, typically orders of magnitude smaller than the variance of \( f \cdot \nabla \log q \), because a well-behaved decoder has similar gradients at nearby \( z \) even when its values vary. Problem 6 computes both variances exactly in a small discrete case where the same comparison can be done by hand. The price of reparameterization is its requirement. A differentiable sampling path exists for location-scale families and a few others, but not for discrete \( z \), which is the entire reason the Gumbel-softmax and straight-through machinery later on this page exists.

The Gaussian KL term in closed form

With \( q = \N(\mu, \diag(\sigma^2)) \) and \( p = \N(0, I) \), the KL term needs no samples at all. Derive it in one dimension, then independence sums it across dimensions. The KL is \( \E_q[\log q(z) - \log p(z)] \) with

Under \( q \), \( \E_q[(z-\mu)^2] = \sigma^2 \) makes the first quadratic term \( -\tfrac12 \), and \( \E_q[z^2] = \mu^2 + \sigma^2 \) makes the second \( \tfrac12(\mu^2 + \sigma^2) \). Collecting,

Each piece penalizes a way of diverging from the prior, means away from zero and variances away from one (from either side, since \( \sigma^2 - \log\sigma^2 \) is minimized at \( \sigma^2 = 1 \)). Closed forms deserve verification, and this one gets it in the implementation section. On a trained encoder's output for a test batch, the analytic KL is 21.0147 nats and a 200,000 sample Monte Carlo estimate of \( \E_q[\log q - \log p] \) is 21.0144 nats, agreeing to 0.0002 nats, within the estimator's standard error of 0.0004.

Amortized inference

Classical variational inference optimizes a separate \( q^{(i)} \) per data point, an inner loop per example. The VAE amortizes. One encoder network maps any \( x \) directly to its posterior parameters \( (\mu_\phi(x), \sigma_\phi(x)) \), so inference on a new point is a forward pass. The saving is enormous and the cost has a name, the amortization gap. The encoder's answer for a given \( x \) is generally worse than what per-example optimization of the same variational family would find, on top of the approximation gap from the family itself. Semi-amortized methods refine the encoder output with a few gradient steps at the cost of the inner loop returning. The deeper point is architectural. The VAE is two coupled networks trained on one objective, a recognition model and a generator, and the ELBO is precisely the objective that makes their cooperation honest, since by the exact decomposition neither can improve it by cheating the other.

Posterior collapse

The ELBO's KL term is minimized, at zero, by \( q_\phi(z \mid x) = p(z) \) for every \( x \), an encoder that ignores its input. If the decoder is powerful enough to model the data unconditionally, autoregressive decoders being the classic case, then \( \theta \) can achieve a good reconstruction term without using \( z \) at all, and the optimizer happily takes the free KL reduction. The latent becomes noise, the "latent variable model" degenerates into the decoder alone, and every downstream use of \( z \) (representation, interpolation, control) silently breaks. This is posterior collapse, and it is a property of the objective's optima, not a bug in any implementation. Early in training the reconstruction gradient through a still-random decoder is weak, the KL gradient is strong and immediate, and collapse is a stable point since a decoder that ignores \( z \) provides no gradient to resurrect it.

The standard fixes each attack one part of that story. KL annealing scales the KL term from 0 to 1 over early training, letting reconstruction establish an informative encoder before regularization bites. Free bits changes the objective to \( \sum_j \max(\lambda, \KL_j) \) per latent dimension (or group), so dimensions carrying less than \( \lambda \) nats are not penalized further and there is no gradient pressure to squeeze the last bits out. Weakening the decoder, PixelCNN decoders with restricted receptive fields, or dropout on the autoregressive context, removes the unconditional escape route. Architectural fixes make the latent hard to ignore, with skip connections from \( z \) into every decoder layer, or the hierarchical top-down designs below. And the \( \delta \)-VAE constrains \( q \) so its KL to the prior is bounded below by construction.

Beta-VAE and the rate-distortion view

Higgins et al. (2017) scale the KL term, \( \L_\beta = \E_q[\log p_\theta(x \mid z)] - \beta \KL(q \| p) \). For \( \beta > 1 \) this is no longer a bound on the likelihood. It is a different trade-off, and the clean way to see what it trades is rate-distortion (Alemi et al., 2018). Define the rate \( R = \E_{x}\!\left[\KL(q_\phi(z\mid x) \| p(z))\right] \), the average number of nats the encoder transmits about \( x \) through the latent channel, and the distortion \( D = -\E_x \E_{q_\phi}[\log p_\theta(x \mid z)] \), the reconstruction cost. The negative ELBO is \( D + R \), and \( \beta \) is the Lagrange multiplier tracing the achievable \( (R, D) \) frontier. Large \( \beta \) buys low rate (aggressive compression, few used latent dimensions) at high distortion, small \( \beta \) the reverse. The same total ELBO can be achieved at many points on the frontier, including the degenerate \( R = 0 \) collapse corner when the decoder family is strong, which restates posterior collapse as a rate-allocation fact rather than a mystery. The trained VAE below sits at \( R = 21.56 \) nats and \( D = 75.52 \) nats on binarized MNIST. Higgins et al.'s empirical claim was that the low-rate regime pressures dimensions to align with independent generative factors, disentanglement. The disentanglement section returns to what of that survives scrutiny.

VQ-VAE, discrete latents by vector quantization

van den Oord et al. (2017) replace the Gaussian latent with a learned discrete codebook \( \{e_1, \dots, e_K\} \subset \R^d \). The encoder produces continuous vectors \( z_e(x) \), a spatial grid of them for images, and each is snapped to its nearest code,

and the decoder reconstructs from \( z_q \). The argmin has zero gradient almost everywhere, so the encoder would receive no training signal. The straight-through estimator supplies one by copying the decoder's gradient at \( z_q \) unchanged to \( z_e \), implemented as \( z_q = z_e + \operatorname{sg}[z_q - z_e] \) with \( \operatorname{sg} \) the stop-gradient. The estimator is biased, it pretends the quantizer is the identity, but the bias is small exactly when quantization is fine, when \( z_e \) sits close to its code. The loss makes that self-fulfilling with three terms,

(as a maximization, with \( \beta \approx 0.25 \)). The codebook term moves each code toward the mean of the encoder outputs assigned to it, k-means with gradients, and touches only the codes. The commitment term moves the encoder toward its assigned code, preventing the encoder outputs from drifting or growing without bound while the codebook chases them. Without it, training visibly oscillates. Many implementations replace the codebook term with an exponential moving average update of the codes, which is the same fixed point with better conditioning. During training the prior over the code grid is uniform, so the KL term of the ELBO is the constant \( \log K \) per position and drops out of the gradient. Generation quality then comes from stage two, fitting an autoregressive model over the discrete code indices (a PixelCNN in the original, a transformer ever since), which turns a \( 32\times32\times3\times 8 \)-bit image modeling problem into a short sequence of categorical tokens. In the measured compression of the implementation below, MNIST digits at 6272 bits of pixels pass through a \( 7 \times 7 \) grid of 6-bit codes, 294 bits, a 21× reduction at a per-pixel reconstruction MSE of 0.0035, with all 64 codes in use at a usage perplexity of 30.2.

VQ-GAN (Esser, Rombach & Ommer, 2021) is the same skeleton pushed to compression rates where MSE reconstruction fails. A pixel-space L2 loss at 16× spatial downsampling produces mush, so the reconstruction term is replaced by a perceptual loss plus a patch discriminator, an adversarial term used purely as a learned distortion metric. The payoff is a codebook whose 16×16 grids are faithful enough that a transformer over 256 tokens generates megapixel-coherent images. This encoder-quantizer-transformer sandwich is the direct ancestor of the image tokenizers inside today's unified multimodal models, and the adversarial ingredient is covered properly on the GAN page.

Hierarchical VAEs

A single-layer VAE with a factorized Gaussian posterior is a weak density model. The fix that kept VAEs competitive is depth in the latent, not just the networks. A hierarchical VAE stacks latents \( z_1, \dots, z_L \) with a top-down generative path \( p(z_L) \prod_\ell p(z_{\ell} \mid z_{>\ell}) \) and, crucially, a top-down inference path sharing those parameters, so the posterior at each level conditions on both the data and the levels above (Ladder VAE lineage). NVAE (Vahdat & Kautz, 2020) engineered this to 30+ latent groups with depthwise separable convolutions and spectral regularization to keep the KL terms stable. VDVAE (Child, 2021) showed a 70+-layer hierarchy trains stably and, in reported figures, reached 2.87 bits per dimension on CIFAR-10, better than PixelCNN++'s 2.92, the first VAE to pass the autoregressive baseline it had trailed for years, while sampling in a single top-down pass. Two readings of that result both matter. The ELBO framework was never the bottleneck, capacity allocation was. And a deep enough chain of conditional Gaussians is itself a kind of learned, coarse-to-fine transport, which in retrospect looks a great deal like a few-step diffusion model with learned rather than fixed noising, one of several places where the families converge ahead of the final section.

Normalizing flows

The change of variables formula, derived

A flow gets exact likelihood and one-pass sampling at once by making the model a deterministic, invertible warp of a simple base density. Let \( z \sim p_z \) on \( \R^D \) and \( x = g(z) \) with \( g \) a diffeomorphism, inverse \( f = g^{-1} \). Probability mass is conserved, so for any region \( A \), \( \P(x \in A) = \P(z \in f(A)) \). In one dimension with \( f \) increasing, differentiating \( F_x(x) = F_z(f(x)) \) gives \( p_x(x) = p_z(f(x))\, f'(x) \), and a decreasing \( f \) flips the sign, so in general \( p_x(x) = p_z(f(x)) |f'(x)| \). In \( D \) dimensions the derivative of the substitution \( z = f(x) \) is the Jacobian matrix, and the volume of an infinitesimal cube maps through a linear map \( J \) to volume \( |\det J| \) times as large, which is the defining property of the determinant. Hence

The second term is the accounting. Where \( f \) compresses \( x \)-space into \( z \)-space (\( |\det J| > 1 \)), the density is inflated, and where it expands, deflated, so the total mass stays one. Composing flows \( f = f_L \circ \dots \circ f_1 \) multiplies Jacobians (chain rule), so log-dets add,

with \( u_0 = x \) and \( u_\ell = f_\ell(u_{\ell-1}) \). Training is direct maximum likelihood on this exact quantity. Two requirements shape every flow architecture. Each layer must be invertible (and the inverse cheap if you want fast sampling), and its Jacobian determinant must be computable in much less than the \( O(D^3) \) a general determinant costs. The whole subfield is a catalog of ways to buy expressive warps under those constraints.

Planar and radial flows

The first neural flows (Rezende & Mohamed, 2015) used layers whose determinant is cheap by the matrix determinant lemma. A planar flow is \( f(z) = z + u\, h(w\T z + b) \), a rank-one perturbation of the identity that bends density across the hyperplane \( w\T z + b = 0 \). Its Jacobian is \( I + u\, h'(w\T z + b)\, w\T \), and \( \det(I + u v\T) = 1 + v\T u \) gives

an \( O(D) \) computation. Radial flows do the same around a point. Both were designed to enrich variational posteriors, and for that they work, but each layer moves density along essentially one direction, so dozens are needed for even 2-D targets, and they have no closed-form inverse. Modern flows won by restructuring the layer, not stacking more of these.

Coupling layers, NICE and RealNVP

The coupling layer (Dinh et al., NICE 2015 and RealNVP 2017) is the decisive trick. Split \( x \) into two blocks \( (x_a, x_b) \). Pass \( x_a \) through unchanged, and transform \( x_b \) elementwise by parameters computed from \( x_a \),

where \( s, t : \R^{d} \to \R^{D-d} \) are arbitrary neural networks, as deep and nonlinear as desired. Invertibility is by construction and never requires inverting the networks. Given \( y \), recover \( x_a = y_a \), recompute the very same \( s(x_a), t(x_a) \), and undo the elementwise map, \( x_b = (y_b - t(x_a)) \odot \exp(-s(x_a)) \). The Jacobian, ordering the blocks \( (a, b) \), is block triangular,

The messy off-diagonal block, the full network Jacobian of \( s \) and \( t \), multiplies into the determinant not at all. The determinant of a block-triangular matrix is the product of the diagonal blocks' determinants, here \( 1 \cdot \prod_j e^{s_j} \). All the nonlinearity is free. NICE is the special case \( s \equiv 0 \), additive coupling, with \( |\det J| = 1 \) exactly (volume preserving, hence the final learned scaling layer in that paper). Since half the coordinates pass through untouched, layers alternate which block is transformed, or permute coordinates between layers. RealNVP's checkerboard and channel masks are this idea mapped onto image tensors, with multi-scale factoring-out of dimensions to keep cost down. The RealNVP trained for this page, 8 coupling layers alternating a 2-D mask, reaches a test NLL of 1.183 nats on standardized two-moons data (a fitted single Gaussian reaches 2.838 nats), reconstructs \( f^{-1}(f(x)) \) to a maximum error of \( 8.6 \times 10^{-6} \) in float32, and its analytic \( \sum_j s_j \) matches the log-determinant of the autograd Jacobian to \( 4.6 \times 10^{-5} \) across 200 test points.

Autoregressive flows, MAF and IAF

A coupling layer conditions half the coordinates on the other half. The fully autoregressive version conditions each coordinate on all of its predecessors. The masked autoregressive flow (Papamakarios, Pavlakou & Murray, 2017) defines the data-to-base direction coordinatewise as

with \( \mu_i, \sigma_i \) computed by one MADE-style masked network. The Jacobian \( \partial z / \partial x \) is triangular (each \( z_i \) depends only on \( x_{\le i} \)), so \( \log|\det| = -\sum_i \log \sigma_i(x_{\lt i}) \), and, because all the conditioners read the observed \( x \), density evaluation is one parallel pass. Sampling inverts the recursion, so \( x_i = \mu_i(x_{\lt i}) + \sigma_i(x_{\lt i}) z_i \) must be computed sequentially, \( D \) passes. In fact a MAF is exactly a Gaussian autoregressive model reread as a flow, which makes the connection between the two families an identity rather than an analogy. The inverse autoregressive flow (Kingma et al., 2016) transposes the asymmetry by making the conditioners read \( z \) instead, \( x_i = \mu_i(z_{\lt i}) + \sigma_i(z_{\lt i}) z_i \). Now sampling is one pass (all of \( z \) is known upfront) and density evaluation of an arbitrary given \( x \) is sequential. The rule of thumb falls out directly. Use MAF for density estimation, IAF for sampling, and for a variational posterior, where the flow only ever evaluates densities of its own samples, IAF is free in both directions, which is precisely why it was introduced to enrich VAE posteriors. Coupling layers are the compromise that is one pass both ways, paying with less expressive conditioning per layer. Parallel WaveNet ran this playbook at industrial scale, distilling a slow-to-sample WaveNet (MAF-like teacher, fast scoring) into an IAF student (fast sampling) by minimizing the KL between them, using each model in the direction it is fast.

Glow and invertible 1×1 convolutions

Alternating fixed masks is a rigid way to mix coordinates between couplings. Glow (Kingma & Dhariwal, 2018) replaces the permutation with a learned invertible \( 1 \times 1 \) convolution, one \( c \times c \) matrix \( W \) applied at every spatial position of a \( h \times w \times c \) tensor, a generalized, learnable channel shuffle. Its log-determinant is \( h \cdot w \cdot \log|\det W| \), and parameterizing \( W = P L U \) with fixed permutation \( P \), unit-lower-triangular \( L \), and upper-triangular \( U \) makes the determinant the product of \( U \)'s diagonal, \( O(c) \) to read off, while keeping inversion \( O(c^2) \) by triangular solves. With activation normalization and affine couplings this trained at the time to the best flow likelihoods on images and famously smooth latent interpolations, and it marks the practical ceiling of the discrete-layer flow program on images. Competitive likelihoods required very deep stacks, and sample fidelity still trailed GANs and, soon after, diffusion.

Continuous normalizing flows

Instead of composing discrete invertible layers, let a vector field flow the points, \( \frac{dz(t)}{dt} = f_\theta(z(t), t) \), integrating from \( t_0 \) to \( t_1 \) (Chen et al., 2018, neural ODEs). Invertibility is automatic, integrate backwards, as long as \( f_\theta \) is Lipschitz (by Picard's theorem trajectories cannot cross). The density evolves by the instantaneous change of variables. As a derivation sketch, over a small step \( \varepsilon \), \( z(t + \varepsilon) = z(t) + \varepsilon f(z(t), t) + O(\varepsilon^2) \) is a near-identity map with Jacobian \( I + \varepsilon \frac{ \partial f}{\partial z} \), and \( \log|\det(I + \varepsilon A)| = \varepsilon \tr A + O(\varepsilon^2) \) (expand the determinant, where off-diagonal products contribute only at order \( \varepsilon^2 \)). Dividing by \( \varepsilon \) and taking the limit,

so the log-likelihood is the base log-density at \( z(t_1) \) plus an integral of a trace along the trajectory. The hard determinant became a trace, computable for any architecture, no triangularity required. A trace of a \( D \times D \) Jacobian still costs \( D \) vector-Jacobian products done exactly. FFJORD (Grathwohl et al., 2019) estimates it unbiasedly with a single one via Hutchinson's estimator,

which holds because \( \E[\varepsilon\T A \varepsilon] = \E[\tr(A \varepsilon \varepsilon\T)] = \tr(A\, \E[\varepsilon\varepsilon\T]) = \tr(A) \), and \( \varepsilon\T A \varepsilon = \varepsilon\T (\nabla_z (f\T \varepsilon)) \) is one backward pass. The result is a flow with unrestricted architecture, at the price of ODE solver calls whose count is set by the stiffness of the learned field, frequently hundreds. That trade sat dormant until the diffusion literature met it from the other side. The probability flow ODE of a diffusion model is a continuous normalizing flow, and flow matching (Lipman et al., 2023) finally made CNF training scale by regressing the vector field directly against analytically known conditional targets instead of differentiating through the solver. That story continues on the diffusion page.

The expressiveness ledger

Flows pay for exactness with structure, and the costs are worth naming. Dimension is preserved. The latent space is \( \R^D \), the same size as the data, so there is no compression and no low-dimensional representation. If the data concentrates near a lower-dimensional manifold, a smooth bijection must stretch enormous volumes to place mass there, which shows up as numerically extreme log-scales. Individual layers are weak by design, so depth substitutes for per-layer capacity. And bits-per-dimension leaderboards reward covering every mode of pixel-level noise, which is not the same as sample fidelity, an evaluation theme formalized below. Where flows are unambiguously the right tool is where exact, differentiable densities are the product, such as variational posteriors, simulation-based inference, and the scientific sampler applications in the closing section.

Energy-based models

The Boltzmann form and the partition function

The most permissive explicit model drops every structural constraint. Let a network assign an unnormalized score to every configuration,

Any density can be written this way (take \( E = -\log p \)), the energy \( E_\theta \) can be any architecture with a scalar output, and composition is trivial, since sums of energies are products of experts. The entire difficulty is \( Z_\theta \), a \( D \)-dimensional integral over the whole space, intractable for any interesting \( E_\theta \). The family's history, Boltzmann machines through modern deep EBMs, is the history of training and sampling without ever computing it.

The two-phase maximum likelihood gradient

The log-likelihood \( \log p_\theta(x) = -E_\theta(x) - \log Z_\theta \) has a gradient with a remarkable structure. Start with the second term,

Substituting back, the gradient ascent direction on the data expectation is

The positive phase pushes energy down on data. The negative phase pushes energy up wherever the model currently puts its mass. At the optimum the two expectations cancel, and the model's own samples are statistically indistinguishable from data as far as \( \nabla_\theta E \) can tell. The formula is exact, elegant, and contains a trap. The negative phase needs samples from \( p_\theta \), which is the intractable sampling problem again, now in the inner loop of training. MCMC is the only general tool, and long chains per gradient step are unaffordable.

Contrastive divergence and its persistent variant

Hinton's contrastive divergence (2002) replaces converged chains with \( k \) MCMC steps started at the data, CD-\( k \) with \( k \) as small as 1. The bet is that early chain motion already points away from where the model overallocates mass near data, so the truncated negative phase is a useful, if biased, gradient. Formally CD approximates the gradient of a difference of KLs and is not the gradient of any function, but it trained RBMs well enough to matter historically. Persistent CD (Tieleman, 2008) instead maintains a set of chains across gradient steps. Each update advances the stored particles a few steps under the current energy and uses them for the negative phase, exploiting the fact that \( \theta \) moves slowly so the chains stay approximately equilibrated. Modern deep EBMs (Du & Mordatch, 2019) are recognizably PCD with Langevin dynamics as the sampler and a replay buffer as the persistent state. The lasting lesson is that EBM training quality is gated by sampler quality, which is exactly the observation that makes the score-based resolution in the next section feel inevitable.

Score matching, making the normalizer irrelevant

Hyvärinen (2005) sidesteps sampling entirely with one observation. The score \( s(x) = \nabla_x \log p(x) \), the gradient with respect to the input rather than the parameters, never sees the normalizer, because \( \nabla_x \log p_\theta(x) = -\nabla_x E_\theta(x) - \nabla_x \log Z_\theta \) and \( Z_\theta \) is constant in \( x \), so its gradient term is zero. Match scores instead of densities,

This still contains the unknown data score, but integration by parts removes it. Expand the square. The term \( \tfrac12 \E\|\nabla \log p_{\text{data}}\|^2 \) is a constant in \( \theta \) and drops. The cross term, written per coordinate \( i \) with \( p = p_{\text{data}} \), is

where the last step is integration by parts in \( x_i \), valid when \( p(x) s_\theta(x) \to 0 \) at infinity (true for any decaying data density and polynomially growing model score). The unknown \( \nabla \log p \) is gone, and only expectations of model quantities under the data remain. Summing over coordinates,

An estimable objective from data samples alone, no \( Z \), no MCMC. The intuition of the two terms is to make the score small on data (data should sit near critical points of the log-density) and to make its divergence negative (those critical points should be local maxima, density curving down in every direction). The catch is the trace of a Jacobian, \( D \) backward passes per example. Two descendants fix the cost. Denoising score matching (Vincent, 2011) corrupts the data with \( q_\sigma(\tilde x \mid x) = \N(\tilde x; x, \sigma^2 I) \) and matches the score of the corrupted marginal. A short computation shows this equals, up to a constant, the regression

whose target is just "point back toward the clean data", trace-free and trivially minibatched. The optimum is the score of the \( \sigma \)-smoothed data density, exact only as \( \sigma \to 0 \) where the estimator's variance explodes (the \( 1/\sigma^2 \) in the target), a tension resolved shortly by using many \( \sigma \) at once. Sliced score matching (Song et al., 2019) keeps the clean-data objective and Hutchinson-izes the trace with random projections \( v\T \nabla_x (v\T s_\theta) \), one backward pass per projection. On the 2-D mixture trained below, DSM at \( \sigma = 0.05 \) recovers the analytic mixture score to a relative squared error of 1.6 percent.

Noise-contrastive estimation

Gutmann & Hyvärinen (2010) attack the normalizer from a different angle, turning density estimation into classification. Draw noise from a known tractable \( p_n \), label data 1 and noise 0, and fit a logistic classifier whose logit is constrained to be \( \log p_\theta(x) - \log p_n(x) \), with \( \log Z \) treated as one more learnable parameter \( c \). The Bayes-optimal classifier for this task has logit exactly \( \log p_{\text{data}}(x) - \log p_n(x) \), so at the optimum (with enough capacity and data) \( p_\theta = p_{\text{data}} \) and \( c \) converges to the true log normalizer. Self-normalization falls out of the classification constraint rather than being imposed. The estimator's efficiency degrades as \( p_n \) diverges from \( p_{\text{data}} \), noise must be hard enough to force the model to learn, which in high dimensions is its practical limit. NCE's legacy is broad. It trained the original word embedding models under softmax normalizers too big to compute, and its ratio-estimation trick is the direct intellectual ancestor of the GAN discriminator below.

Score-based models and diffusion, briefly

The EBM section left a precise unfinished problem. The two-phase gradient and the very idea of sampling from an energy model need MCMC, and the natural chain in continuous spaces is Langevin dynamics,

which converges to \( p \) as \( \eta \to 0 \) and needs only the score. Score matching provides exactly that quantity without ever normalizing. So the pieces snap together. Learn \( s_\theta \approx \nabla \log p_{\text{data}} \) by (denoising) score matching, then sample by Langevin. Run naively this fails for an instructive reason. Where the data has no mass, the training objective, an expectation under the data, never constrains the score, so a chain initialized far from the manifold follows arbitrary directions, and separated modes mix impossibly slowly. Song & Ermon (2019) fixed both with one move, learning the scores of the data smoothed at many noise levels \( \sigma_1 > \dots > \sigma_L \), then sampling with annealed Langevin, coarse levels guiding the chain from anywhere in space onto the manifold, fine levels sharpening it. Heavily smoothed densities have well-defined scores everywhere and connected modes, while lightly smoothed ones carry the detail.

Diffusion models are this idea with a fixed forward noising process and a variational derivation, and the two meet in an identity. The denoising network of a diffusion model trained to predict the noise \( \varepsilon \) from \( x_t = \alpha_t x_0 + \sigma_t \varepsilon \) is, by the DSM argument above, estimating the score of the noised marginal, \( s_\theta(x_t) = -\hat\varepsilon_\theta(x_t) / \sigma_t \), so the "diffusion objective" and multi-scale denoising score matching are the same objective with a particular weighting over noise levels. Taking the noise schedule continuous turns the forward process into an SDE, the generative direction into a reverse-time SDE driven by the learned score, and its deterministic counterpart, the probability flow ODE, into a continuous normalizing flow with exact likelihoods, which closes the loop with the flow section. The full treatment, DDPM and its ELBO, samplers, guidance, latent diffusion, distillation to few steps, lives on the diffusion page. The trained toy system here is on this page because it is the smallest complete instance of the resolution. A DSM-trained score on a 2-D mixture, sampled with Langevin from pure noise, lands 49.79 percent of 20,000 chains in the right-hand mode of an equal mixture with marginal standard deviations within 2.5 percent of the truth.

Implicit models, the GAN in one derivation

Every family so far represents a density or its score. The GAN (Goodfellow et al., 2014) discards representation entirely. The map \( x = G_\theta(z) \), \( z \sim \N(0, I) \), defines a sampler whose density is never written, and training replaces likelihood with a two-sample test. A discriminator \( D \) is trained to distinguish data from generator samples while the generator is trained to fool it,

The inner maximization has a closed-form solution. Pointwise, the integrand is \( p_{\text{data}}(x) \log D + p_g(x) \log(1 - D) \), of the form \( a \log D + b \log(1-D) \), maximized (set the derivative \( a/D - b/(1-D) \) to zero) at

Substituting \( D^\ast \) back into the objective and adding and subtracting \( \log 2 \) in each term rearranges it to

the Jensen-Shannon divergence. At discriminator optimality, generator training minimizes a symmetric, bounded divergence estimated purely from samples. That is the framework-level content of GANs, likelihood-free training against a learned test, and it explains both their strength and their failure mode. The strength is that nothing in the objective rewards covering data modes the discriminator cannot detect are missing from finite batches, so all capacity goes to per-sample fidelity, which is why GAN images were sharp years before likelihood models caught up. The failure is that, for the same reason, mode collapse, the generator concentrating on a subset of the data distribution, is not directly penalized and must be fought with auxiliary machinery. In the language of the KL section, practical GAN training sits firmly on the mode-seeking side of the ledger. The dynamics of the minimax game, non-saturating losses, Wasserstein critics, spectral normalization, R1 penalties, and where GANs still win in 2026, are the subject of the adversarial modeling page. Here the essential point is that the discriminator is a density-ratio estimator (compare NCE above, where the same trick is used with a fixed noise distribution and yields a consistent density estimate, while the GAN makes the noise adversarial and gives up the density).

Evaluation, what the numbers do and do not mean

Likelihood, nats, and bits per dimension

For explicit models the natural score is held-out log-likelihood, reported as bits per dimension to normalize across resolutions, \( \text{bpd} = \text{NLL}_{\text{nats}} / (D \ln 2) \). Comparability is narrower than the shared units suggest. A discrete model (softmax or discretized-logistic pixels) reports a genuine probability mass. A continuous model (a flow) reports a density, which can be made arbitrarily large by concentrating near grid points, so continuous models must be evaluated on dequantized data, \( x + u \) with \( u \sim \mathrm{Uniform}[0,1)^D \), and Theis et al. (2016) showed the resulting continuous likelihood lower-bounds the discrete one, making the comparison legitimate in exactly one direction. Flow++ (Ho et al., 2019, from the Berkeley line of flow work) tightened the same bound by learning the dequantization noise variationally, itself worth several hundredths of a bit per dimension. Even then, numbers are only comparable within one data treatment. Dynamically binarized MNIST, fixed binarization, and 8-bit grayscale are three different tasks with three incompatible leaderboards, and an ELBO (a bound) is not exchange-comparable with an exact autoregressive likelihood at the same value. When the bound gap itself is the question, annealed importance sampling gives a practical estimate of the true marginal likelihood of decoder-based models (Wu et al., 2017, a Toronto and CMU collaboration). The trained models on this page illustrate the discipline. The MADE's 97.36 nats and the VAE's bound of 97.08 nats on the same dynamically binarized MNIST are legitimately comparable (and the IWAE-100 estimate of 93.07 nats shows the VAE is actually the better density model here), while the RealNVP's 1.183 nats on continuous 2-D data lives in a different universe entirely.

Likelihood does not predict sample quality

Theis, van den Oord & Bethge (2016) made precise what practitioners had observed. Log-likelihood and sample quality can be decoupled in both directions, and in high dimensions the decoupling is numerically severe. For the key construction, let \( q \) be an excellent model and consider the contaminated mixture \( \tilde q = 0.01\, q + 0.99\, \nu \) with \( \nu \) pure noise. Samples from \( \tilde q \) are 99 percent garbage. Its log-likelihood barely moves,

For CIFAR-10 images, where totals run to thousands of nats (3.0 bpd \( \times \) 3072 dimensions \( \times \ln 2 \approx 6390 \) nats), a 4.6-nat penalty is a rounding error, so a model can be near-optimal in likelihood while almost never producing a good sample. The reverse direction is just as easy. A model that memorizes the training set produces perfect-looking samples and assigns \( -\infty \) (or catastrophically low) log-likelihood to any held-out point. High-dimensional likelihood is dominated by how well the model covers the bulk of imperceptible pixel-level variability, forward-KL mode-covering again, while perception is a measure concentrated on structure. Neither number is lying. They measure different projections of the same object, and any evaluation that reports one as if it certified the other is wrong.

IS and FID

Sample-based metrics route both model samples and (for FID) real samples through a pretrained classifier. The Inception Score, \( \exp( \E_x \KL(p(y \mid x) \| p(y)) ) \), rewards samples that are individually classifiable (sharp \( p(y \mid x) \)) and collectively diverse (flat marginal \( p(y) \)). Its flaws are structural. It never touches the data distribution, so a model emitting one perfect prototype per class scores near the maximum, and it inherits the classifier's training distribution, making it meaningless off ImageNet-like data. The Fréchet Inception Distance fits Gaussians to real and generated feature activations and computes the closed-form Wasserstein-2 distance between them,

FID does compare against data and correlates with human judgment well enough to have become the field's default, but its assumptions deserve suspicion in proportion to its ubiquity. Features are not Gaussian, so models differing only in higher moments can tie. The estimator is strongly biased at small sample counts (standard practice fixes \( n = 50{,}000 \) precisely because values are not comparable across \( n \)). It inherits Inception's texture and class biases, and as a single scalar it cannot say whether a bad score means low fidelity or low coverage. The two-number fix, precision and recall for generative models (Sajjadi et al., 2018, refined into a nearest-neighbor manifold estimate by Kynkäänniemi et al., 2019), estimates manifold overlap in both directions. Precision, the fraction of generated samples landing on the (estimated) real manifold, reads fidelity. Recall, the fraction of real samples covered by the generated manifold, reads mode coverage. On these axes the historical stereotype becomes measurable. GANs sit high precision low recall, likelihood models the reverse, and a report of either alone can be gamed by trading the other.

Two-sample tests

The statistically principled frame for all sample-based evaluation is the two-sample test. Given \( \{x_i\} \sim p_{\text{data}} \) and \( \{\tilde x_j\} \sim p_\theta \), test \( H_0 : p_{\text{data}} = p_\theta \). Kernel maximum mean discrepancy gives a closed-form statistic with an honest null distribution. Classifier two-sample tests train a discriminator on held-out splits and read its accuracy (near 50 percent under \( H_0 \)), which is exactly the GAN discriminator recycled as an evaluator. The uncomfortable general finding is that for high-dimensional data these tests are either weak (simple kernels) or become model-comparison problems themselves (learned critics). Evaluation of generative models is genuinely unsolved. The practical discipline is to report metrics that fail differently, a likelihood or bound, an FID-class metric, precision and recall, and to state the sample sizes and data treatment, because every one of these numbers can be moved by protocol alone.

Discrete data and gradient estimation

The Gumbel-max trick and its softmax relaxation

Reparameterization needs a differentiable path from parameters to samples, and discrete variables do not have one. This single obstruction shapes VQ-VAE training, discrete latent models, and RL-adjacent objectives alike. The starting point is the Gumbel-max trick, an exact reparameterization of the categorical whose only flaw is an argmax. Let \( g_i = -\log(-\log u_i) \) with \( u_i \sim \mathrm{Uniform}(0,1) \) be i.i.d. Gumbel noise, and then

proved as Problem 6 below. Jang, Gu & Poole and Maddison, Mnih & Teh (both 2017) relax the argmax to a temperature-controlled softmax, the Gumbel-softmax or Concrete distribution,

a point on the simplex rather than a vertex, with a fully differentiable path from \( \pi \) to \( y \). As \( \tau \to 0 \), \( y \) converges to the one-hot argmax sample, and the estimator's gradient variance diverges. As \( \tau \) grows, \( y \) smears toward uniform and the surrogate's bias against the true discrete objective grows. Temperature is a bias-variance dial, typically annealed downward during training, and when the downstream computation genuinely needs a hard sample, the straight-through Gumbel-softmax uses the argmax forward and the soft \( y \)'s gradient backward, stacking the straight-through bias on top of the relaxation bias in exchange for exact discreteness on the forward pass.

REINFORCE and baselines

The score-function estimator from the VAE section is the fully general alternative, requiring nothing but the ability to score samples, \( \nabla_\phi \E[f] = \E[ f(z) \nabla_\phi \log q_\phi(z) ] \). Any baseline \( b \) independent of the sample can be subtracted without bias, because

and a well-chosen \( b \) (a running mean of \( f \), a learned value network, or the multi-sample leave-one-out baselines of VIMCO and RLOO) can collapse the variance by orders of magnitude. It can also raise it. Problem 7 works a case where the folk choice \( b = \E[f] \) makes the estimator 16 times noisier while the optimal baseline drives the variance to exactly zero. The modern division of labor is clean. Use reparameterization wherever a continuous path exists, relaxations when a discrete node must pass useful gradients and bias is tolerable, REINFORCE with strong baselines when the sample must be exactly discrete and unbiasedness matters, and straight-through when a biased but cheap and stable signal empirically wins, as it does in VQ-VAE, whose quantizer is the most widely deployed straight-through estimator in existence.

Controllable and conditional generation

Conditioning mechanisms

Everything above extends to \( p_\theta(x \mid c) \) for conditioning information \( c \), and the architectural options are shared across families. Concatenate an embedding of \( c \) to the input or latent (the original conditional VAE and GAN recipe), modulate normalization statistics with \( c \)-dependent scale and shift (FiLM, adaptive group norm, the class-conditional batch norm of BigGAN, the adaLN blocks of diffusion transformers), or cross-attend from the generation stream to a sequence encoding of \( c \), the mechanism carrying text prompts into every modern image model, covered mechanically on the cross-attention page. In autoregressive models, conditioning degenerates gracefully into prefixing, where \( c \) is just earlier tokens, which is why one decoder-only architecture absorbed instruction following, translation, and multimodal inputs without structural change.

Guidance

Score-based models admit a form of post-hoc conditioning that has become its own subfield. Bayes' rule on scores reads

so a sampler following the unconditional score can be steered by adding the input-gradient of any classifier evaluated on the partially generated \( x \), which is classifier guidance. Training a classifier on noised inputs is a nuisance, so classifier-free guidance trains one network with \( c \) randomly dropped, obtaining both conditional and unconditional scores, and samples along the extrapolated direction

with \( w > 1 \) overshooting the conditional score. Formally \( w \) exponentiates the implied classifier term \( p(c \mid x)^w \), sharpening the conditional at the cost of diversity and calibration, a knowing, controllable retreat from the learned distribution toward its most \( c \)-typical regions. The practice and its artifacts are detailed on the diffusion page. The framework-level point is that guidance is possible because score-based models expose the gradient of the log-density. Families that expose a likelihood (autoregressive) steer instead by re-weighting token distributions, and implicit models must bake all control into training.

Disentanglement and its impossibility result

The recurring hope for latent-variable models is that individual latent coordinates align with human-meaningful factors of variation, and beta-VAE-era papers reported exactly that. Locatello et al. (2019) tested the hope at scale, 12,000 trained models across methods, datasets, and seeds, and paired it with a theorem, that unsupervised disentanglement is impossible without inductive biases. The argument's core is a construction. For any model with latent \( z \sim p(z) \) factorized, there exists a smooth invertible map \( h \) with \( p(h(z)) = p(z) \) marginally (for an isotropic Gaussian, any rotation works) that entangles every coordinate. Composing the decoder with \( h^{-1} \) yields a second model with identical marginal \( p(x) \) and identical objective value whose latents are maximally entangled. No purely observational objective can distinguish the two, so any observed disentanglement is contributed by architecture, regularizer, data regularities, or luck of the seed, and their experiments found the seed's share embarrassingly large. The constructive program that survives is weak supervision. Pairs differing in few factors, temporal structure, or interventional data provably restore identifiability, which connects to independent-component-analysis theory, and, further out, to the causal-representation-learning agenda.

Latent arithmetic

Weaker than disentanglement and far more robust is the empirical smoothness of learned latent spaces. Interpolations between latents decode to semantic interpolations, and direction vectors computed as differences of group means, the average latent of smiling faces minus that of neutral faces, transfer across individuals when added, the word-analogy arithmetic of embedding spaces reproduced in generative latents (systematically explored in DCGAN-era work and StyleGAN's style space). Two cautions keep it honest. For a Gaussian prior, linear interpolation between two typical latents passes through points of atypically small norm, so spherical interpolation is the correct default. And a direction that edits one attribute cleanly for one region of latent space routinely entangles others elsewhere, which is the local, practical shadow of the Locatello result.

The current picture

As of 2026 the division of labor across families is fairly stable, and worth stating plainly because it is the map everything else on this page justifies. Autoregressive transformers own language outright and have re-entered images and interleaved multimodal generation through learned discrete tokenizers. A VQ-GAN-descended encoder turns pixels into tokens, one decoder-only transformer models text and image tokens in a single stream, and generation is next-token prediction throughout, the architecture of the unified multimodal models every frontier lab now ships. Diffusion and its flow-matching reformulation own continuous media, such as images, video, audio waveforms, molecular conformations, and robot action trajectories, with the center of gravity moving from pixel-space U-Nets to latent-space diffusion transformers, and from thousand-step samplers to few-step consistency and distilled models. Pure VAEs survive less as generators than as infrastructure. The latent space that latent diffusion runs in is a (KL- or VQ-regularized) autoencoder, and the ELBO survives as the training objective underneath diffusion's derivation. Normalizing flows hold the scientific niche where their exact, cheap likelihood is the entire point, with Boltzmann generators proposing equilibrium molecular configurations, trivializing maps in lattice field theory, and neural posterior estimation in simulation-based inference for astrophysics and neuroscience. GANs persist where single-forward-pass sampling and high per-sample fidelity beat distributional coverage, in super-resolution, codec enhancement, and as adversarial loss terms inside other systems, VQ-GAN being the canonical case.

The deeper trend is convergence. Diffusion's probability flow ODE is a continuous normalizing flow. Flow matching trains that flow by regression, and its optimal-transport variants make the flows straighter and the samplers shorter. Consistency models distill the flow into a one- or two-step map, which is functionally the GAN's contract, a single pass from noise to data, reached from the likelihood side. Autoregressive sampling composes conditional inverse-CDF maps, a triangular transport, which is what MAF makes explicit. So the 2026 synthesis reads all the families as one statement, a generative model is a learned transport from a simple distribution to the data distribution, and the families differ in how the transport is structured (triangular, bijective, stochastic, iterative) and which functional of it the training loss constrains (its likelihood, a bound, its score field, or a critic's verdict). That single sentence is the compression of this page, and it is the frame in which new methods are now proposed and evaluated.

Worked problems

Implementation

Everything in this section was trained on the H100 80GB used throughout the site (PyTorch 2.7.0, with JAX versions as direct translations of the same models). All quoted numbers are from those runs and are stored in classes/data/generative.json. Each snippet pairs the model with the numerical check that catches its characteristic bug, the VAE's closed-form KL against Monte Carlo, the flow's log-det against autograd, the MADE mask against the Jacobian, and the score model against the analytic score of a known density.

A VAE with the exact ELBO, and a Monte Carlo check on the KL

The model is deliberately the classic one, a 784-400-400 encoder to a 20-dimensional Gaussian latent, mirrored decoder to Bernoulli logits, dynamically binarized MNIST. Trained 30 epochs (20 seconds of wall clock), it reaches a test ELBO of −97.08 nats per image, split as 75.52 nats of reconstruction and 21.56 nats of KL. An IWAE-100 evaluation of the same model gives \( \log p_\theta(x) \ge -93.07 \) nats. The check at the bottom verifies the closed-form KL derivation against a 200,000-sample Monte Carlo estimate of \( \E_q[\log q - \log p] \), analytic 21.0147 nats against Monte Carlo 21.0144 ± 0.0004, agreeing to 0.0002 nats. If those two numbers disagree beyond the standard error, the KL formula, a sign, a factor of two, or a \( \log \sigma^2 \) vs \( \log \sigma \) convention, is wrong. This half-line check catches all four classic bugs.

import math, torch, torch.nn as nn, torch.nn.functional as F

class VAE(nn.Module):
    def __init__(self, d_lat=20, d_h=400):
        super().__init__()
        self.enc = nn.Sequential(nn.Linear(784, d_h), nn.ReLU(),
                                 nn.Linear(d_h, d_h), nn.ReLU())
        self.mu, self.logvar = nn.Linear(d_h, d_lat), nn.Linear(d_h, d_lat)
        self.dec = nn.Sequential(nn.Linear(d_lat, d_h), nn.ReLU(),
                                 nn.Linear(d_h, d_h), nn.ReLU(),
                                 nn.Linear(d_h, 784))
    def forward(self, x):                      # x: (B, 784) in {0,1}
        h = self.enc(x)
        mu, logvar = self.mu(h), self.logvar(h)          # (B, 20) each
        z = mu + torch.exp(0.5 * logvar) * torch.randn_like(mu)  # reparam
        return self.dec(z), mu, logvar         # logits: (B, 784)

def neg_elbo(logits, x, mu, logvar):
    # reconstruction: -log p(x|z), summed over pixels  -> (B,)
    bce = F.binary_cross_entropy_with_logits(logits, x,
                                             reduction="none").sum(-1)
    # KL( N(mu, sigma^2) || N(0, I) ), closed form, summed over latents
    kl = 0.5 * (mu.pow(2) + logvar.exp() - 1.0 - logvar).sum(-1)
    return bce + kl                            # (B,) negative ELBO in nats

vae = VAE().cuda()
opt = torch.optim.Adam(vae.parameters(), lr=1e-3)
for epoch in range(30):
    for x in loader:                           # x: (B, 784) in [0,1]
        xb = torch.bernoulli(x)                # dynamic binarization
        logits, mu, logvar = vae(xb)
        loss = neg_elbo(logits, xb, mu, logvar).mean()
        opt.zero_grad(); loss.backward(); opt.step()
# test ELBO: -97.08 nats  (recon 75.52 + KL 21.56); IWAE-100: -93.07 nats

# --- verify the closed-form KL against Monte Carlo -----------------
with torch.no_grad():
    h = vae.enc(xb); mu, logvar = vae.mu(h), vae.logvar(h)
    kl_closed = (0.5 * (mu.pow(2) + logvar.exp() - 1 - logvar)
                 .sum(-1)).mean()
    eps = torch.randn(200_000, *mu.shape, device=mu.device)  # (K, B, 20)
    z = mu + torch.exp(0.5 * logvar) * eps
    log_q = -0.5 * (eps.pow(2) + math.log(2*math.pi) + logvar).sum(-1)
    log_p = -0.5 * (z.pow(2) + math.log(2*math.pi)).sum(-1)
    kl_mc = (log_q - log_p).mean()
    # measured: closed 21.0147, MC 21.0144 +/- 0.0004  -> formula verified

import math, jax, jax.numpy as jnp, optax

def init_vae(key, d_h=400, d_lat=20):
    ks = jax.random.split(key, 7)
    g = lambda k, i, o: jax.random.normal(k, (i, o)) * jnp.sqrt(2.0 / i)
    return {"e1": (g(ks[0], 784, d_h), jnp.zeros(d_h)),
            "e2": (g(ks[1], d_h, d_h), jnp.zeros(d_h)),
            "mu": (g(ks[2], d_h, d_lat), jnp.zeros(d_lat)),
            "lv": (g(ks[3], d_h, d_lat), jnp.zeros(d_lat)),
            "d1": (g(ks[4], d_lat, d_h), jnp.zeros(d_h)),
            "d2": (g(ks[5], d_h, d_h), jnp.zeros(d_h)),
            "d3": (g(ks[6], d_h, 784), jnp.zeros(784))}

def lin(p, x): W, b = p; return x @ W + b

def encode(p, x):                              # x: (B, 784)
    h = jax.nn.relu(lin(p["e2"], jax.nn.relu(lin(p["e1"], x))))
    return lin(p["mu"], h), lin(p["lv"], h)    # (B, 20), (B, 20)

def decode(p, z):
    h = jax.nn.relu(lin(p["d2"], jax.nn.relu(lin(p["d1"], z))))
    return lin(p["d3"], h)                     # logits (B, 784)

def neg_elbo(p, x, key):
    mu, logvar = encode(p, x)
    eps = jax.random.normal(key, mu.shape)
    z = mu + jnp.exp(0.5 * logvar) * eps       # reparameterization
    logits = decode(p, z)
    bce = jnp.sum(jnp.maximum(logits, 0) - logits * x
                  + jnp.log1p(jnp.exp(-jnp.abs(logits))), axis=-1)
    kl = 0.5 * jnp.sum(mu**2 + jnp.exp(logvar) - 1.0 - logvar, axis=-1)
    return jnp.mean(bce + kl)

opt = optax.adam(1e-3)

@jax.jit
def step(p, opt_state, x, key):
    loss, grads = jax.value_and_grad(neg_elbo)(p, x, key)
    updates, opt_state = opt.update(grads, opt_state)
    return optax.apply_updates(p, updates), opt_state, loss

# --- closed-form KL vs Monte Carlo ---------------------------------
def kl_check(p, x, key, K=200_000):
    mu, logvar = encode(p, x)
    closed = jnp.mean(0.5 * jnp.sum(mu**2 + jnp.exp(logvar) - 1 - logvar,
                                    axis=-1))
    eps = jax.random.normal(key, (K, *mu.shape))       # (K, B, 20)
    z = mu + jnp.exp(0.5 * logvar) * eps
    log_q = -0.5 * jnp.sum(eps**2 + math.log(2*math.pi) + logvar, -1)
    log_p = -0.5 * jnp.sum(z**2 + math.log(2*math.pi), -1)
    return closed, jnp.mean(log_q - log_p)     # should agree to ~1e-3

A RealNVP coupling flow, checked against autograd

Eight affine coupling layers on 2-D two-moons data, alternating which coordinate passes through, with the log-scale bounded by a \( 2\tanh \) for stability and the final linear layer zero-initialized so the flow starts at the identity. Trained 6000 steps (about a minute), test NLL is 1.183 nats per point against 2.838 for the best single Gaussian. The two checks are the ones that matter for any flow. Composing forward and inverse returns the input to \( 8.6 \times 10^{-6} \) (float32 roundoff through 8 layers), and the network's cheap \( \sum_j s_j \) log-det matches \( \log|\det| \) of the brute-force autograd Jacobian to \( 4.6 \times 10^{-5} \) over 200 points. The second check is cheap only in low dimension, which is precisely why toy dimensions are the right place to certify the layer before scaling it.

import math, torch, torch.nn as nn

class Coupling(nn.Module):
    """y = mask*x + (1-mask)*(x*exp(s(x_a)) + t(x_a)),  x_a = mask*x."""
    def __init__(self, mask, d_h=128):
        super().__init__()
        self.register_buffer("mask", mask)     # (2,), 1 = pass through
        self.net = nn.Sequential(nn.Linear(2, d_h), nn.Tanh(),
                                 nn.Linear(d_h, d_h), nn.Tanh(),
                                 nn.Linear(d_h, 4))
        self.net[-1].weight.data.zero_()       # start at identity
        self.net[-1].bias.data.zero_()
    def _st(self, xa):
        s, t = self.net(xa).chunk(2, dim=-1)   # (B,2), (B,2)
        s = 2.0 * torch.tanh(s)                # bounded log-scale
        return s * (1 - self.mask), t * (1 - self.mask)
    def forward(self, x):                      # data -> base
        xa = x * self.mask
        s, t = self._st(xa)
        y = xa + (1 - self.mask) * (x * torch.exp(s) + t)
        return y, s.sum(-1)                    # per-example log|det J|
    def inverse(self, y):                      # base -> data
        ya = y * self.mask
        s, t = self._st(ya)                    # same nets, forward only
        return ya + (1 - self.mask) * ((y - t) * torch.exp(-s))

class RealNVP(nn.Module):
    def __init__(self, n=8):
        super().__init__()
        self.layers = nn.ModuleList(
            Coupling(torch.tensor([1., 0.] if i % 2 == 0 else [0., 1.]))
            for i in range(n))
    def forward(self, x):
        ld = torch.zeros(len(x), device=x.device)
        for l in self.layers:
            x, d = l(x); ld = ld + d
        return x, ld
    def inverse(self, z):
        for l in reversed(self.layers):
            z = l.inverse(z)
        return z
    def log_prob(self, x):                     # exact log-likelihood
        z, ld = self(x)
        return -0.5 * (z.pow(2) + math.log(2*math.pi)).sum(-1) + ld

flow = RealNVP().cuda()
opt = torch.optim.Adam(flow.parameters(), lr=1e-3)
for step in range(6000):
    x = sample_moons(1024)                     # (1024, 2), standardized
    loss = -flow.log_prob(x).mean()
    opt.zero_grad(); loss.backward(); opt.step()
# test NLL: 1.183 nats/point (single-Gaussian baseline: 2.838)

# --- check 1: invertibility ----------------------------------------
z, _ = flow(x_test)
print((flow.inverse(z) - x_test).abs().max())  # 8.6e-6 in float32

# --- check 2: log-det vs autograd jacobian -------------------------
xi = x_test[:1]                                # one point at a time
J = torch.autograd.functional.jacobian(
    lambda v: flow(v)[0], xi).reshape(2, 2)
ld_auto = torch.linalg.slogdet(J)[1]
ld_net = flow(xi)[1]
print((ld_auto - ld_net).abs())                # max 4.6e-5 over 200 pts

import math, jax, jax.numpy as jnp, optax

def init_coupling(key, d_h=128):
    k1, k2 = jax.random.split(key)
    g = lambda k, i, o: jax.random.normal(k, (i, o)) * jnp.sqrt(1.0 / i)
    return {"w1": g(k1, 2, d_h), "b1": jnp.zeros(d_h),
            "w2": g(k2, d_h, d_h), "b2": jnp.zeros(d_h),
            "w3": jnp.zeros((d_h, 4)), "b3": jnp.zeros(4)}  # identity init

def st(p, xa, mask):
    h = jnp.tanh(xa @ p["w1"] + p["b1"])
    h = jnp.tanh(h @ p["w2"] + p["b2"])
    s, t = jnp.split(h @ p["w3"] + p["b3"], 2, axis=-1)
    return 2.0 * jnp.tanh(s) * (1 - mask), t * (1 - mask)

def fwd_layer(p, x, mask):                     # data -> base
    xa = x * mask
    s, t = st(p, xa, mask)
    return xa + (1 - mask) * (x * jnp.exp(s) + t), s.sum(-1)

def inv_layer(p, y, mask):                     # base -> data
    ya = y * mask
    s, t = st(p, ya, mask)
    return ya + (1 - mask) * ((y - t) * jnp.exp(-s))

MASKS = jnp.array([[1., 0.] if i % 2 == 0 else [0., 1.]
                   for i in range(8)])         # (8, 2)

def fwd(params, x):
    ld = jnp.zeros(x.shape[0])
    for p, m in zip(params, MASKS):
        x, d = fwd_layer(p, x, m); ld = ld + d
    return x, ld

def log_prob(params, x):                       # exact log-likelihood
    z, ld = fwd(params, x)
    return -0.5 * jnp.sum(z**2 + math.log(2*math.pi), -1) + ld

loss_fn = lambda params, x: -jnp.mean(log_prob(params, x))

# --- check 1: invertibility ----------------------------------------
def inv(params, z):
    for p, m in zip(reversed(params), MASKS[::-1]):
        z = inv_layer(p, z, m)
    return z
# jnp.abs(inv(params, fwd(params, x)[0]) - x).max()  ~ 1e-5 (float32)

# --- check 2: log-det vs jax.jacfwd --------------------------------
def check_logdet(params, xi):                  # xi: (2,)
    J = jax.jacfwd(lambda v: fwd(params, v[None])[0][0])(xi)  # (2, 2)
    ld_auto = jnp.linalg.slogdet(J)[1]
    ld_net = fwd(params, xi[None])[1][0]
    return jnp.abs(ld_auto - ld_net)           # ~1e-5

MADE masks, and the Jacobian as proof

The mask construction is fifteen lines, and the property it buys, that output \( i \) depends on no input \( j \ge i \), is exactly verifiable. The \( 784 \times 784 \) Jacobian of logits with respect to inputs has maximum absolute entry 0.0 on and above the diagonal, a structural zero from the masked weights, not a small number. Trained 40 epochs on dynamically binarized MNIST (17 seconds), the model reaches 97.36 nats, 0.179 bits per pixel. The sampling loop at the bottom is the family's cost made concrete, 784 sequential forward passes, 172.3 ms for a batch of 64 where a single (training or likelihood) pass takes 0.43 ms, a 404× gap on the same weights and hardware.

import numpy as np, torch, torch.nn as nn, torch.nn.functional as F

def made_masks(d_in, d_hidden):
    """Degree-based masks (Germain et al., 2015), natural ordering."""
    rng = np.random.RandomState(42)
    degrees = [np.arange(1, d_in + 1)]         # inputs get 1..D
    for dh in d_hidden:                        # hidden degrees in [prev_min, D-1]
        degrees.append(rng.randint(degrees[-1].min(), d_in, size=dh))
    masks = [torch.tensor((degrees[l+1][:, None] >= degrees[l][None, :])
                          .astype(np.float32))
             for l in range(len(d_hidden))]
    # output k models p(x_k | x_<k): strict inequality kills self-loops
    masks.append(torch.tensor((degrees[0][:, None] > degrees[-1][None, :])
                              .astype(np.float32)))
    return masks

class MaskedLinear(nn.Linear):
    def __init__(self, d_in, d_out, mask):
        super().__init__(d_in, d_out)
        self.register_buffer("mask", mask)     # (d_out, d_in)
    def forward(self, x):
        return F.linear(x, self.weight * self.mask, self.bias)

class MADE(nn.Module):
    def __init__(self, d=784, hidden=(1024, 1024)):
        super().__init__()
        ms, dims = made_masks(d, list(hidden)), [d, *hidden]
        layers = []
        for l in range(len(hidden)):
            layers += [MaskedLinear(dims[l], dims[l+1], ms[l]), nn.ReLU()]
        layers += [MaskedLinear(dims[-1], d, ms[-1])]
        self.net = nn.Sequential(*layers)
    def forward(self, x):                      # (B, 784) -> logits (B, 784)
        return self.net(x)

made = MADE().cuda()   # train: BCE-with-logits summed over pixels
# test NLL after 40 epochs: 97.36 nats = 0.179 bits/pixel

# --- proof of autoregressiveness: the Jacobian is strictly lower-tri
x0 = torch.bernoulli(x_test[:1])
J = torch.autograd.functional.jacobian(made, x0).reshape(784, 784)
print(J.triu(0).abs().max())                   # exactly 0.0

# --- sampling: one forward pass per pixel --------------------------
x = torch.zeros(64, 784, device="cuda")
for i in range(784):                           # 172.3 ms total
    p_i = torch.sigmoid(made(x)[:, i])         # vs 0.43 ms for one pass
    x[:, i] = torch.bernoulli(p_i)             # -> 404x slower to sample

import numpy as np, jax, jax.numpy as jnp

def made_masks(d_in, d_hidden, seed=42):
    rng = np.random.RandomState(seed)
    degrees = [np.arange(1, d_in + 1)]
    for dh in d_hidden:
        degrees.append(rng.randint(degrees[-1].min(), d_in, size=dh))
    masks = [jnp.asarray((degrees[l+1][:, None] >= degrees[l][None, :]),
                         dtype=jnp.float32)
             for l in range(len(d_hidden))]
    masks.append(jnp.asarray((degrees[0][:, None] > degrees[-1][None, :]),
                             dtype=jnp.float32))
    return masks

def init_made(key, d=784, hidden=(1024, 1024)):
    dims = [d, *hidden, d]
    ks = jax.random.split(key, len(dims) - 1)
    ws = [jax.random.normal(k, (o, i)) * jnp.sqrt(2.0 / i)
          for k, i, o in zip(ks, dims[:-1], dims[1:])]
    bs = [jnp.zeros(o) for o in dims[1:]]
    return {"w": ws, "b": bs, "m": made_masks(d, list(hidden))}

def made_fwd(p, x):                            # (B, 784) -> logits
    h = x
    for l in range(len(p["w"]) - 1):
        h = jax.nn.relu(h @ (p["w"][l] * p["m"][l]).T + p["b"][l])
    return h @ (p["w"][-1] * p["m"][-1]).T + p["b"][-1]

def nll(p, x):                                 # nats per image
    logits = made_fwd(p, x)
    bce = jnp.maximum(logits, 0) - logits * x \
          + jnp.log1p(jnp.exp(-jnp.abs(logits)))
    return jnp.mean(bce.sum(-1))

# --- autoregressive property check ---------------------------------
J = jax.jacrev(lambda v: made_fwd(p, v[None])[0])(x0)   # (784, 784)
assert jnp.abs(jnp.triu(J, 0)).max() == 0.0    # structural zero

# --- sequential sampling (one jit-compiled step per pixel) ---------
def sample(p, key, n=64):
    x = jnp.zeros((n, 784))
    for i in range(784):                       # inherently sequential
        key, sub = jax.random.split(key)
        prob = jax.nn.sigmoid(made_fwd(p, x)[:, i])
        x = x.at[:, i].set(jax.random.bernoulli(sub, prob).astype(x.dtype))
    return x

VQ-VAE quantization with straight-through

The quantizer is the part worth reading closely, with nearest-code lookup by expanded squared distance, the three-term loss with its two stop-gradients placed exactly as the derivation demands, and the one-line straight-through identity. Trained 15 epochs on MNIST with a 7×7 grid over a 64-entry codebook, reconstruction MSE is 0.0035 per pixel through a 294-bit bottleneck (21× smaller than the 6272-bit input), with all 64 codes active at usage perplexity 30.2. Perplexity is the health metric to watch, since a collapsing codebook shows up as perplexity sliding toward 1 long before reconstructions visibly degrade.

import torch, torch.nn as nn, torch.nn.functional as F

class VectorQuantizer(nn.Module):
    def __init__(self, K=64, D=16, beta=0.25):
        super().__init__()
        self.emb = nn.Embedding(K, D)          # codebook e: (K, D)
        self.emb.weight.data.uniform_(-1/K, 1/K)
        self.beta = beta
    def forward(self, z_e):                    # z_e: (B, D, H, W)
        B, D, H, W = z_e.shape
        flat = z_e.permute(0, 2, 3, 1).reshape(-1, D)      # (BHW, D)
        # ||z - e||^2 = ||z||^2 - 2 z.e + ||e||^2  -> (BHW, K)
        d2 = (flat.pow(2).sum(1, keepdim=True)
              - 2 * flat @ self.emb.weight.t()
              + self.emb.weight.pow(2).sum(1))
        idx = d2.argmin(1)                     # (BHW,) code indices
        z_q = self.emb(idx).view(B, H, W, D).permute(0, 3, 1, 2)
        # codebook term moves codes; commitment term holds the encoder
        vq_loss = F.mse_loss(z_q, z_e.detach()) \
                + self.beta * F.mse_loss(z_e, z_q.detach())
        z_q = z_e + (z_q - z_e).detach()       # straight-through: d z_q = d z_e
        counts = F.one_hot(idx, self.emb.num_embeddings).float().mean(0)
        perplexity = torch.exp(-(counts * (counts + 1e-10).log()).sum())
        return z_q, vq_loss, perplexity, idx

# training step: total = F.mse_loss(decoder(z_q), x) + vq_loss
# measured (MNIST, 7x7 grid, K=64): recon MSE 0.0035/pixel,
# perplexity 30.2, 64/64 codes used, 294-bit latent vs 6272-bit input

import jax, jax.numpy as jnp

def quantize(codebook, z_e, beta=0.25):
    """codebook: (K, D); z_e: (B, H, W, D) channels-last."""
    K, D = codebook.shape
    flat = z_e.reshape(-1, D)                              # (BHW, D)
    d2 = (jnp.sum(flat**2, 1, keepdims=True)
          - 2.0 * flat @ codebook.T
          + jnp.sum(codebook**2, 1))                       # (BHW, K)
    idx = jnp.argmin(d2, axis=1)                           # (BHW,)
    z_q = codebook[idx].reshape(z_e.shape)                 # (B, H, W, D)
    codebook_loss = jnp.mean((z_q - jax.lax.stop_gradient(z_e))**2)
    commit_loss = jnp.mean((z_e - jax.lax.stop_gradient(z_q))**2)
    vq_loss = codebook_loss + beta * commit_loss
    # straight-through: forward value z_q, gradient flows to z_e
    z_q = z_e + jax.lax.stop_gradient(z_q - z_e)
    counts = jnp.mean(jax.nn.one_hot(idx, K), axis=0)      # (K,)
    perplexity = jnp.exp(-jnp.sum(counts * jnp.log(counts + 1e-10)))
    return z_q, vq_loss, perplexity, idx

# loss(params, x):
#   z_e = encoder(params, x)                  # (B, 7, 7, 16)
#   z_q, vq_loss, perp, _ = quantize(params["codebook"], z_e)
#   recon = decoder(params, z_q)
#   return jnp.mean((recon - x)**2) + vq_loss

Score matching on a 2-D density with a known answer

Denoising score matching on a two-Gaussian mixture whose score has a closed form, which makes the evaluation exact rather than aesthetic. After 8000 steps at \( \sigma = 0.05 \), the training loss settles at 835 against the irreducible floor \( \E\|\varepsilon\|^2/\sigma^2 = 800 \) (the target's variance, which no network can remove, so watching DSM loss go "only" to 835 and knowing that is correct requires having computed the floor). The learned field matches the analytic score to 1.6 percent relative squared error on-distribution, and 2000 steps of Langevin from \( \N(0, 9) \) noise land 49.79 percent of 20,000 chains in the right mode with marginal standard deviations within 2.5 percent of truth. This is the EBM sampling problem, solved end to end at toy scale by learning the score instead of the energy.

import math, torch, torch.nn as nn

mus = torch.tensor([[-1.5, 0.], [1.5, 0.]], device="cuda")
SIG2 = 0.35**2                                 # component variance

def sample_mix(n):                             # equal-weight 2-Gaussian mix
    comp = torch.randint(0, 2, (n,), device="cuda")
    return mus[comp] + math.sqrt(SIG2) * torch.randn(n, 2, device="cuda")

def true_score(x):                             # analytic, for evaluation
    d = torch.stack([-((x - m).pow(2).sum(-1)) / (2*SIG2) for m in mus], -1)
    w = torch.softmax(d, -1)                   # responsibilities (B, 2)
    return (w[:, :1] * (mus[0] - x) + w[:, 1:] * (mus[1] - x)) / SIG2

score = nn.Sequential(nn.Linear(2, 128), nn.SiLU(),
                      nn.Linear(128, 128), nn.SiLU(),
                      nn.Linear(128, 2)).cuda()
opt = torch.optim.Adam(score.parameters(), lr=1e-3)
SIGMA = 0.05                                   # DSM noise level
for step in range(8000):
    x = sample_mix(2048)
    eps = torch.randn_like(x)
    x_t = x + SIGMA * eps
    target = -eps / SIGMA                      # = grad log q(x_t | x)
    loss = (score(x_t) - target).pow(2).sum(-1).mean()
    opt.zero_grad(); loss.backward(); opt.step()
# loss -> 835 vs irreducible floor 2/sigma^2 = 800
# relative score error vs analytic: 1.6% on 50k mixture samples

# --- Langevin dynamics from pure noise -----------------------------
x = 3.0 * torch.randn(20_000, 2, device="cuda")
eta = 2e-3
for t in range(2000):
    x = x + eta * score(x) + math.sqrt(2 * eta) * torch.randn_like(x)
# right-mode fraction 0.4979 (true 0.5); marginal stds within 2.5%

import math, jax, jax.numpy as jnp, optax

MUS = jnp.array([[-1.5, 0.], [1.5, 0.]])
SIG2, SIGMA = 0.35**2, 0.05

def sample_mix(key, n):
    k1, k2 = jax.random.split(key)
    comp = jax.random.randint(k1, (n,), 0, 2)
    return MUS[comp] + math.sqrt(SIG2) * jax.random.normal(k2, (n, 2))

def true_score(x):                             # (B, 2) -> (B, 2)
    d = -jnp.stack([jnp.sum((x - m)**2, -1) for m in MUS], -1) / (2*SIG2)
    w = jax.nn.softmax(d, -1)
    return (w[:, :1] * (MUS[0] - x) + w[:, 1:] * (MUS[1] - x)) / SIG2

def score_net(p, x):
    h = jax.nn.silu(x @ p["w1"] + p["b1"])
    h = jax.nn.silu(h @ p["w2"] + p["b2"])
    return h @ p["w3"] + p["b3"]

def dsm_loss(p, key):
    k1, k2 = jax.random.split(key)
    x = sample_mix(k1, 2048)
    eps = jax.random.normal(k2, x.shape)
    x_t = x + SIGMA * eps
    return jnp.mean(jnp.sum((score_net(p, x_t) + eps / SIGMA)**2, -1))

@jax.jit
def langevin_step(p, x, key, eta=2e-3):
    noise = jax.random.normal(key, x.shape)
    return x + eta * score_net(p, x) + math.sqrt(2 * eta) * noise

def sample(p, key, n=20_000, steps=2000):
    key, k0 = jax.random.split(key)
    x = 3.0 * jax.random.normal(k0, (n, 2))    # start far off-manifold
    for _ in range(steps):
        key, k = jax.random.split(key)
        x = langevin_step(p, x, k)
    return x

How it is done in practice

The gap between these derivations and a production system is mostly about scale-induced failure modes and pipeline structure. Production VAE practice is dominated by the latent autoencoders inside image and video generation stacks, trained with a KL weight orders of magnitude below \( \beta = 1 \) (the point is a well-conditioned latent space, not a tight bound), a perceptual loss, and a light adversarial term, then frozen forever while diffusion or a transformer trains on top. The engineering issues are unglamorous and decisive. Variance of the latent space must be normalized before the second stage. Decoder artifacts cap final quality no matter how good the prior model is. And tokenizer or autoencoder changes invalidate every downstream checkpoint, so they version like database schemas.

Discrete-token pipelines center on the codebook, whose failure modes at scale, dead codes and index collapse, are fought with EMA updates, code resets, lower-dimensional code projection, and increasingly by replacing the learned codebook outright with finite scalar quantization or lookup-free binary quantization, which trade expressiveness per code for guaranteed usage. Autoregressive serving is its own discipline (KV caching, speculative decoding, batching), and the 404× sample-versus-score gap measured above is the toy version of exactly the asymmetry those systems engineer around. Diffusion serving optimizes the other axis, steps times cost-per-step, with latent spaces to shrink the step and distillation and consistency training to shrink the count, with the H100-class arithmetic on the diffusion page. Flow deployments in the sciences live or die on exactness details this page's checks foreshadow. An invertibility error of \( 10^{-6} \) per layer compounds over hundreds of layers, and samplers used for free-energy estimates must pair the flow with an importance-weight or MCMC correction so that model bias becomes statistical error rather than silent physics error.

Evaluation practice is the part most worth internalizing because it is where systems quietly go wrong. It means fixed evaluation noise and prompts to make FID-class metrics comparable across checkpoints, precision-recall alongside any scalar metric, likelihoods only within matched data treatments, and human or preference-model evaluation as the final arbiter for anything user-facing, since every proxy metric above has a documented way to be gamed.

The current research frontier

Four threads dominate the last few years. First is the consolidation of continuous generation around transport. Flow matching (Lipman et al. at Meta AI, 2023) and rectified flow (Liu et al., 2023) train continuous flows by regressing velocity fields along prescribed noise-to-data paths, with stochastic interpolants (Albergo & Vanden-Eijnden at NYU, 2023) as the general theory containing diffusion and deterministic flows as special cases, while consistency models (Song et al., OpenAI, 2023) and their training-without-distillation refinements compress the transport to one or two steps. Production image and video systems from Stability, Black Forest Labs, Google DeepMind, and OpenAI moved onto these formulations within roughly two years of the papers.

Second is discrete generation beyond left-to-right. Masked/absorbing discrete diffusion (the D3PM lineage at Google, MaskGIT-style parallel decoding, and MDLM and score-entropy discrete diffusion from academic groups including work at Cornell and elsewhere) now trains language models whose likelihoods approach autoregressive ones while decoding in parallel, and any-order autoregressive models blur the same boundary from the other side. Whether parallel discrete generation displaces next-token prediction anywhere that matters is one of the live questions of 2026.

Third is tokenizer research as a first-class problem. Finite scalar quantization (Google), lookup-free quantization in MAGVIT-v2, residual and multi-scale quantization (the VAR line from Peking and Bytedance showing next-scale autoregression beating diffusion baselines on ImageNet), and continuous-token autoregression with diffusion heads (MAR, MIT-affiliated authors) all attack the same question of what interface should sit between continuous media and sequence models. Unified multimodal models, Chameleon-style early-fusion token models at Meta, transfusion-style hybrids interleaving diffusion and next-token objectives, and their counterparts at DeepSeek, Alibaba's Qwen, and Tsinghua-linked Zhipu, are the deployment ground for whichever answer wins.

Fourth is the quieter theory-and-science thread, Boltzmann generators and transferable equilibrium samplers (Noé's group and DeepMind), lattice field theory flows (MIT and DeepMind collaborations), simulation-based inference stacks built on MAF/NSF-style flows (the Tübingen-led sbi ecosystem, with heavy use in astrophysics and neuroscience, and related tooling from ETH Zurich and EPFL groups), and diffusion models as priors for inverse problems in imaging and structural biology, where the score-based view earns its keep because posteriors, not samples, are the deliverable. Across all four threads the evaluation problem from this page recurs at larger stakes, and papers increasingly report precision-recall-style coverage metrics and downstream-task numbers instead of a single scalar.

Open source to read

Ordered roughly from pedagogical to production. For each, the first file worth opening.

• pytorch/examples, the canonical minimal VAE, the same model trained above. Open vae/main.py. The whole model and ELBO fit on one screen.
• openai/pixel-cnn, PixelCNN++ reference code. Open pixel_cnn_pp/nn.py and read discretized_mix_logistic_loss against the formula in the autoregressive section, edge bins and all.
• google-deepmind/distrax, JAX distributions and bijectors with clean flow composition. Open distrax/_src/bijectors/masked_coupling.py for a production coupling layer to compare with the one above.
• pyro-ppl/pyro, deep probabilistic programming, with the ELBO as reusable infrastructure rather than a hand-derived loss. Open examples/vae/vae.py to see the same VAE written as a model/guide pair.
• VincentStimper/normalizing-flows, a well-organized PyTorch flow zoo (RealNVP, Glow, MAF, residual flows). Open normflows/flows/affine/coupling.py.
• rtqichen/torchdiffeq, the neural-ODE solvers behind CNFs, with adjoint backpropagation. Open examples/ode_demo.py first, then torchdiffeq/_impl/odeint.py.
• CompVis/taming-transformers, VQ-GAN as shipped, the tokenizer recipe latent diffusion was built on. Open taming/modules/vqvae/quantize.py. The quantizer above is its minimal core.
• lucidrains/vector-quantize-pytorch, every modern quantizer variant (EMA, FSQ, LFQ, residual VQ) behind one interface, the fastest way to survey the design space. Open vector_quantize_pytorch/vector_quantize_pytorch.py.
• openai/consistency_models, the one-step end of the transport convergence story. Open cm/karras_diffusion.py, where the consistency training and distillation losses live.

Common misconceptions

"The VAE's KL term is the regularizer you tune like weight decay." At \( \beta = 1 \) the KL term is not a regularizer at all. It is a component of a single derived quantity, the ELBO, and the sum is a bound on \( \log p(x) \). Reweighting it is sometimes the right engineering move, but it changes the objective from a likelihood bound to a rate-distortion trade-off, and claims about likelihood made from a \( \beta \ne 1 \) model are unfounded.

"Blurry VAE samples prove likelihood is the wrong objective." The blur has specific, separable causes. One is a factorized Gaussian or Bernoulli decoder likelihood whose optimum given an uncertain \( z \) is a pixelwise mean. The other is posterior/prior mismatch at generation time. Hierarchical VAEs with the same maximum-likelihood objective produce sharp samples, and diffusion models, trained on a likelihood bound too, produce the sharpest of all. The objective was never the problem. The decoder's conditional independence was.

"Flows compute a determinant, so they must be slow." This is backwards. The entire flow design program is choosing layers whose log-determinant is a byproduct of the forward pass, \( O(D) \) for coupling and autoregressive layers, read off a triangular diagonal. Nobody ever materializes a general Jacobian except in a unit test, which is exactly what the autograd check above is.

"Higher held-out likelihood means better samples." Theis et al.'s mixture argument makes the decoupling quantitative. Contamination by 99 percent noise costs at most \( \ln 100 \approx 4.6 \) nats per image, about 0.002 bpd on CIFAR-10, invisible next to run-to-run variance, while destroying sample quality. The converse fails too. A training-set memorizer has superb samples and abysmal held-out likelihood. Likelihood measures coverage in high-dimensional detail space, while perception measures typicality of structure.

"Score matching approximates maximum likelihood." It is a different divergence, the Fisher divergence, with different weighting. Score error is averaged under the data density, so regions the data rarely visits are unconstrained regardless of how much probability the model puts there. That is not a defect of implementation but of the objective, and the multi-noise-scale construction exists precisely to re-constrain those regions by smoothing mass into them.

"The straight-through estimator is a hack with no analysis." It is biased, and its bias is controlled. The estimator treats the quantizer as identity, so the error shrinks with the quantization residual, which the commitment loss explicitly minimizes. The pairing of straight-through gradients with a loss term that keeps the operating point in the low-bias regime is a design, not an accident, and it is load-bearing in every deployed VQ tokenizer.

"GANs learn the data distribution." The minimax derivation says the JS divergence is minimized only at discriminator optimality, which finite training never provides. Nothing in practical GAN training certifies coverage, and measured recall confirms modes are dropped. The honest statement is that a GAN learns a sampler whose outputs a bounded critic cannot distinguish from data, a two-sample-test guarantee, not a density-estimation one.

"Diffusion made the other families obsolete." Diffusion won one regime, continuous media where iterative refinement is affordable. Language remains autoregressive. Multimodal token models are VQ plus autoregression. The latent space diffusion runs in is a VAE. Scientific samplers needing exact densities are flows. One-step distilled models are converging on the GAN's operating point with an adversarial term frequently reintroduced. The families merged. They were not eliminated.

Self-check

References

1. Goodfellow, I., Bengio, Y., Courville, A. Deep Learning. MIT Press, 2016. Chapters 16-20. deeplearningbook.org
2. Murphy, K. Probabilistic Machine Learning: Advanced Topics. MIT Press, 2023. Parts IV-V. probml.github.io/pml-book/book2
3. Tomczak, J. Deep Generative Modeling. Springer, 2nd ed. 2024. doi:10.1007/978-3-031-64087-2
4. Kingma, D. P., Welling, M. Auto-Encoding Variational Bayes. ICLR 2014. arXiv:1312.6114. An Introduction to Variational Autoencoders. Foundations and Trends in ML, 2019. arXiv:1906.02691
5. Rezende, D. J., Mohamed, S., Wierstra, D. Stochastic Backpropagation and Approximate Inference in Deep Generative Models. ICML 2014. arXiv:1401.4082
6. Higgins, I., et al. beta-VAE: Learning Basic Visual Concepts with a Constrained Variational Framework. ICLR 2017. openreview.net/forum?id=Sy2fzU9gl
7. van den Oord, A., Kalchbrenner, N., Kavukcuoglu, K. Pixel Recurrent Neural Networks. ICML 2016. arXiv:1601.06759. Conditional Image Generation with PixelCNN Decoders. NeurIPS 2016. arXiv:1606.05328. WaveNet: A Generative Model for Raw Audio, 2016. arXiv:1609.03499
8. Salimans, T., Karpathy, A., Chen, X., Kingma, D. P. PixelCNN++. ICLR 2017. arXiv:1701.05517
9. van den Oord, A., Vinyals, O., Kavukcuoglu, K. Neural Discrete Representation Learning (VQ-VAE). NeurIPS 2017. arXiv:1711.00937. Razavi, A., van den Oord, A., Vinyals, O. Generating Diverse High-Fidelity Images with VQ-VAE-2. NeurIPS 2019. arXiv:1906.00446
10. Esser, P., Rombach, R., Ommer, B. Taming Transformers for High-Resolution Image Synthesis (VQ-GAN). CVPR 2021. arXiv:2012.09841
11. Germain, M., Gregor, K., Murray, I., Larochelle, H. MADE: Masked Autoencoder for Distribution Estimation. ICML 2015. arXiv:1502.03509
12. Dinh, L., Krueger, D., Bengio, Y. NICE: Non-linear Independent Components Estimation. ICLR workshop 2015. arXiv:1410.8516. Dinh, L., Sohl-Dickstein, J., Bengio, S. Density Estimation Using Real NVP. ICLR 2017. arXiv:1605.08803
13. Kingma, D. P., Dhariwal, P. Glow: Generative Flow with Invertible 1x1 Convolutions. NeurIPS 2018. arXiv:1807.03039
14. Papamakarios, G., Pavlakou, T., Murray, I. Masked Autoregressive Flow for Density Estimation. NeurIPS 2017. arXiv:1705.07057. Papamakarios, G., et al. Normalizing Flows for Probabilistic Modeling and Inference. JMLR 2021. arXiv:1912.02762
15. Chen, R. T. Q., Rubanova, Y., Bettencourt, J., Duvenaud, D. Neural Ordinary Differential Equations. NeurIPS 2018. arXiv:1806.07366. Grathwohl, W., et al. FFJORD. ICLR 2019. arXiv:1810.01367
16. Hinton, G. E. Training Products of Experts by Minimizing Contrastive Divergence. Neural Computation 14(8), 2002. doi:10.1162/089976602760128018
17. Hyvärinen, A. Estimation of Non-Normalized Statistical Models by Score Matching. JMLR 6, 2005. jmlr.org/papers/v6/hyvarinen05a
18. Vincent, P. A Connection Between Score Matching and Denoising Autoencoders. Neural Computation 23(7), 2011. doi:10.1162/NECO_a_00142
19. Song, Y., Ermon, S. Generative Modeling by Estimating Gradients of the Data Distribution. NeurIPS 2019. arXiv:1907.05600
20. Gutmann, M., Hyvärinen, A. Noise-Contrastive Estimation. AISTATS 2010. proceedings.mlr.press/v9/gutmann10a
21. Jang, E., Gu, S., Poole, B. Categorical Reparameterization with Gumbel-Softmax. ICLR 2017. arXiv:1611.01144. Maddison, C. J., Mnih, A., Teh, Y. W. The Concrete Distribution. ICLR 2017. arXiv:1611.00712
22. Theis, L., van den Oord, A., Bethge, M. A Note on the Evaluation of Generative Models. ICLR 2016. arXiv:1511.01844
23. Locatello, F., et al. Challenging Common Assumptions in the Unsupervised Learning of Disentangled Representations. ICML 2019. arXiv:1811.12359
24. Goodfellow, I., et al. Generative Adversarial Networks. NeurIPS 2014. arXiv:1406.2661
25. Vahdat, A., Kautz, J. NVAE: A Deep Hierarchical Variational Autoencoder. NeurIPS 2020. arXiv:2007.03898. Child, R. Very Deep VAEs Generalize Autoregressive Models. ICLR 2021. arXiv:2011.10650