VAE: Auto-Encoding Variational Bayes

An autoencoder you can sample from, trained by backpropagating through a random draw.

A latent-variable model makes data by decoding a random code, but the probability it assigns to a real datapoint is an integral over every possible code, and nobody can compute it. Kingma and Welling train the model anyway: they maximize a lower bound on that probability, give the job of guessing codes to a second network, and rewrite the random draw so that gradients pass through it.

Trained on 28×28 MNIST digits with a 20-number code and about 816,000 weights, the model turns any 20 draws from a standard Gaussian into a new digit.

The paper, Auto-Encoding Variational Bayes by Diederik Kingma and Max Welling (2013), gives a general recipe for fitting models with continuous hidden variables. Its neural-network instance, which the paper calls the variational auto-encoder, is what is now usually called the VAE. The page builds it from four ideas, in order:

1. the model and the integral that makes its likelihood uncomputable;
2. a lower bound on the likelihood that can be computed;
3. an encoder network that stands in for the unknown posterior, which turns the bound into an autoencoder;
4. a rewrite of the random sampling step, the reparameterization trick, so backprop reaches the encoder.

After that come the full training step with real shapes, the learned latent space, and the paper's experiments.

The model and its integral

A latent-variable model explains each datapoint $\mathbf{x}$ by an unobserved code $\mathbf{z}$. The paper's version has two parts. The code comes from a fixed prior, a standard Gaussian $p(\mathbf{z}) = \mathcal{N}(\mathbf{0}, \mathbf{I})$ with no learned parameters. A decoder network with weights $\boldsymbol{\theta}$ maps the code to a distribution over data, $p_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{z})$. To generate, draw $\mathbf{z}$ from the prior and then draw $\mathbf{x}$ from the decoder.

For binarized MNIST that means: $\mathbf{z}$ is a vector of 20 numbers, the decoder is a one-hidden-layer network (20 → 500 → 784), and its output is 784 probabilities, one per pixel, each the chance that pixel is on. Figure 1 shows the generative direction: codes from a round Gaussian on the left pass through the decoder and land on the thin, curved set where real data lives.

Training means choosing $\boldsymbol{\theta}$ so that real data is probable under the model. The probability of one datapoint has to account for every code that could have produced it, each weighted by its prior probability. This is the marginal likelihood, also called the evidence:

The integrand is cheap. For one pair $(\mathbf{x}, \mathbf{z})$, the joint $p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}) = p(\mathbf{z})\,p_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{z})$ costs one decoder pass and one Gaussian density. The integral is the expensive part. The decoder is a nonlinear network, so there is no closed form, and a grid over 20 dimensions with only 10 points per dimension already needs $10^{20}$ decoder passes for a single digit. Random codes from the prior do not help much either: almost all of them decode to something unlike this particular $\mathbf{x}$ and contribute nearly nothing, so a plain Monte Carlo average needs an enormous number of samples.

The same integral blocks inference, the question of which codes explain a given $\mathbf{x}$. By Bayes' rule the posterior is

and its denominator is the same marginal likelihood. The numerator can be evaluated at any $\mathbf{z}$; the normalizer cannot. That rules out the standard tool for latent-variable models, the EM (expectation-maximization) algorithm, whose E-step needs the exact posterior. It also rules out mean-field variational Bayes, the usual fallback, because its update equations need integrals that are just as intractable for a neural-network decoder. The paper states both failures in Section 2.1.

The paper gives up on the exact posterior. A second network, the recognition model or encoder $q_{\boldsymbol{\phi}}(\mathbf{z}\mid\mathbf{x})$ with weights $\boldsymbol{\phi}$, takes a datapoint and returns an approximate posterior, a Gaussian over the codes that might have produced it. The next section shows how that approximation yields a quantity that can be computed and optimized.

A lower bound on log p(x)

With any encoder in hand, $\log p_{\boldsymbol{\theta}}(\mathbf{x})$ splits exactly into two terms. This is the paper's equation (1):

The first term is a KL divergence, $D_{\mathrm{KL}}(q\,\|\,p) = \mathbb{E}_q[\log q - \log p]$, which measures how far the encoder's Gaussian is from the true posterior. It is never negative, and it is zero only when the two distributions are identical. The left side is a fixed number once $\boldsymbol{\theta}$ is fixed, so the second term $\mathcal{L}$ can never exceed it. The paper calls $\mathcal{L}$ the variational lower bound; later papers call it the evidence lower bound, or ELBO. It equals $\log p_{\boldsymbol{\theta}}(\mathbf{x})$ only when the encoder matches the true posterior exactly, and the KL term, called the gap below, is how much the bound falls short.

The bound itself can be written without the posterior, as the paper's equation (2):

Everything inside the expectation can be evaluated: the encoder's own density and the cheap joint. The expectation is over codes drawn from the encoder, which is easy to sample. The same bound falls out of Jensen's inequality. Write $p_{\boldsymbol{\theta}}(\mathbf{x})$ as an average under $q$, then use the fact that the log of an average is at least the average of the log:

and the difference between the two sides is exactly the KL gap in (1).

Raising $\mathcal{L}$ with respect to the encoder shrinks the gap, even though $\log p_{\boldsymbol{\theta}}(\mathbf{x})$ is never computed. With $\boldsymbol{\theta}$ held fixed, the left side of (1) does not move, so every nat the bound gains is a nat removed from the gap, and the encoder moves toward the true posterior. Raising it with respect to $\boldsymbol{\theta}$ pushes up a floor under $\log p_{\boldsymbol{\theta}}(\mathbf{x})$.

Figure 2 checks (1) with real numbers on a model small enough to integrate. The code is one number, $z \sim \mathcal{N}(0,1)$, and the decoder is $x \mid z \sim \mathcal{N}(z^2,\, 0.5^2)$. An observed $x = 2.5$ could have come from $z$ near $+1.6$ or near $-1.6$, so the true posterior has two peaks. Here $\log p(x) = -2.55$ nats, computed by numerical integration, and the ELBO of the encoder's Gaussian comes in closed form. At the starting Gaussian ($\mu = 0.4$, $\sigma = 1$) the ELBO is $-9.18$, so the gap is 6.63 nats. Press “fit q” to run gradient ascent on the ELBO and watch where the Gaussian settles.

From the starting point, gradient ascent moves the Gaussian onto the right-hand peak ($\mu = 1.51$, $\sigma = 0.165$) and stops with an ELBO of $-3.25$ and a gap of 0.70 nats. The gap does not close because no single Gaussian can match a two-peaked posterior. It settles on one peak instead of spreading over both because of which way round the KL is written. $D_{\mathrm{KL}}(q\,\|\,p)$ averages $\log q - \log p$ over samples from $q$, so it charges heavily for any place where $q$ puts mass and $p$ has almost none (the valley at $z = 0$), and charges nothing for a peak of $p$ that $q$ ignores. This is the reverse KL, and it is mode-seeking. Covering one of two equal, well-separated peaks costs about $\log 2 = 0.69$ nats, because the posterior density at that peak is half of what it would be if the peak held all the mass. Drag $x$ to 0 or below and the posterior has a single peak; the fitted gap falls to 0.03 nats at $x = 0$ and 0.002 at $x = -1$.

For training, the bound is rewritten once more. Splitting the joint in (2) as $p_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{z})\,p(\mathbf{z})$ and collecting the two terms that involve the prior gives the paper's equation (3), written here with the reconstruction first:

This KL has a different second argument from the one in (1). The gap in (1) compares the encoder with the true posterior $p_{\boldsymbol{\theta}}(\mathbf{z}\mid\mathbf{x})$, which is unknown. The regularizer in (3) compares it with the prior $\mathcal{N}(\mathbf{0},\mathbf{I})$, which is fixed and known, and a later section computes it in closed form. Both are reverse KLs with $q$ first.

Equation (3) also fixes the direction of optimization. The paper maximizes $\mathcal{L}$ by stochastic gradient ascent, so the regularizer enters with a minus sign. Most code minimizes a loss equal to the negative bound, where the same KL is a positive penalty:

The two describe the same training. On this page, every signed quantity is labelled as either the bound (maximized) or the loss (minimized).

The bound as an autoencoder

Read (3) as a computation. The encoder takes $\mathbf{x}$ and outputs a Gaussian over codes. A code $\mathbf{z}$ is drawn from it. The decoder takes $\mathbf{z}$ and assigns a probability to the original $\mathbf{x}$; the reconstruction term is the log of that probability, averaged over the draw. That is an autoencoder, an encoder and a decoder joined through a narrow code, with two differences from the ordinary kind: the code is a random draw from a distribution, and the KL term pulls each such distribution toward the prior. The paper calls the reconstruction term “a negative reconstruction error in auto-encoder parlance”. For binary pixels it is the negative of the binary cross-entropy loss.

The two terms pull in opposite directions. The reconstruction term prefers a narrow Gaussian far from its neighbours for each datapoint, so the decoder can tell datapoints apart. The KL term prefers every Gaussian to be the prior itself, $\mu = 0$ and $\sigma = 1$, which carries no information about $\mathbf{x}$ at all. Figure 4 makes the balance adjustable. It multiplies the KL term by a weight $\beta$; the paper's bound is $\beta = 1$. The decoder is held fixed at $x \mid z \sim \mathcal{N}(2z,\, 0.5^2)$ so that the best encoder for each of eight datapoints has a closed form: $\mu = 2x/(4 + 0.25\beta)$ and $\sigma^2 = 0.25\beta/(4 + 0.25\beta)$.

At $\beta = 1$ each Gaussian has $\sigma = 0.243$, which is the exact posterior for this model. The reconstruction averages $-0.72$ nats, the KL $1.35$, and the ELBO $-2.07$, the highest value any setting of the slider reaches. The Gaussians overlap, their average roughly traces the prior, and 4% of the prior's mass is left uncovered, all of it beyond $|z| = 1.9$.

Slide left toward $\beta = 0$. The Gaussians narrow into spikes and the reconstruction term approaches its best value, $-0.23$ nats, but at $\beta = 0.01$ the KL is already 3.64 nats and 71% of the prior's mass lies between codes. At $\beta = 0$ exactly, which is a plain autoencoder, every code is a point, the KL to the prior is infinite, and every prior draw lands in a hole. In a trained model the decoder would only ever see codes at the eight spikes, so a sample from the prior gives the decoder a code it was never trained on, and a plain autoencoder cannot be used as a generator.

To the right, at $\beta = 10$ the Gaussians have widened to $\sigma = 0.62$ and pulled in toward the origin, and the decoded means have shrunk from the data (the reconstruction falls to $-4.37$). At $\beta = 1000$ every Gaussian is the prior, the KL is 0.00, and every datapoint decodes to nearly 0. The code now carries no information about $\mathbf{x}$. Later work calls this posterior collapse, and it happens to individual latent dimensions in real VAEs when the decoder can model the data without them.

One encoder for every datapoint

The encoder also changes how inference is paid for. Classical variational inference fits a separate Gaussian for each datapoint: hold $\mathbf{x}$ fixed and run an inner optimization over that point's own $\boldsymbol{\mu}$ and $\boldsymbol{\sigma}$ until its bound stops improving. A thousand datapoints means a thousand optimizations, and a new datapoint at test time means one more, started from scratch. The paper instead trains one network that maps any $\mathbf{x}$ straight to its $\boldsymbol{\mu}$ and $\boldsymbol{\sigma}$. The optimization cost moves into training the weights $\boldsymbol{\phi}$, once, and inference on any datapoint, seen or unseen, is a single forward pass. Later papers call this amortized inference, since the cost is spread across every datapoint the encoder will ever see.

The shared network has a price. With a fixed number of weights it cannot, in general, hit each datapoint's best Gaussian exactly, and the shortfall is called the amortization gap (a later term, not the paper's). It adds to the gap from Figure 2, which comes from the Gaussian shape itself. Figure 5 sets the costs side by side.

Past the crossover the amortized line stays nearly flat while the classical line keeps rising in proportion to D; each datapoint encoded after training adds one forward pass and no optimization.

The KL term in closed form

The paper's encoder is a Gaussian with a diagonal covariance, its equation (9). For each datapoint the encoder network outputs a mean vector $\boldsymbol{\mu}$ and a vector of standard deviations $\boldsymbol{\sigma}$, one entry per latent dimension:

The prior is also Gaussian, so the KL between them has an exact formula and never needs sampling. Appendix B of the paper gives it in the form that is added to the maximized bound, summed over the $J$ latent dimensions:

Each dimension contributes independently. A dimension with $\mu_j = 0$ and $\sigma_j = 1$ contributes $\tfrac{1}{2}(1 + 0 - 0 - 1) = 0$, the largest possible value; every other setting contributes a negative number. The $-\mu_j^2$ term charges for a mean away from the origin. The pair $\log \sigma_j^2 - \sigma_j^2$ peaks at $\sigma_j = 1$ and falls on both sides, slowly for a large $\sigma_j$ and without limit as $\sigma_j \to 0$, where the log goes to $-\infty$. A loss that is minimized uses the positive KL, with all four signs inside the sum flipped:

Two worked values: one dimension with $\mu = -1.1$ and $\sigma = 1.8$ gives $\tfrac{1}{2}(1 + 1.176 - 1.21 - 3.24) = -1.14$ nats in the bound (the starting point of Figure 6), and a 20-dimensional code with every $\mu_j = 0.5$ and $\sigma_j = 0.5$ costs $0.443$ nats per dimension, 8.86 nats in total, in the loss.

With the KL handled exactly, one term of (3) still needs sampling: the reconstruction $\mathbb{E}_{q_{\boldsymbol{\phi}}}\!\big[\log p_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{z})\big]$. Estimating its value is easy: draw a few codes and average. Estimating its gradient with respect to the encoder is the problem the rest of the method solves.

Gradients through a random draw

Write the reconstruction's integrand as $f(\mathbf{z}) = \log p_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{z})$. The decoder's gradient is straightforward: draw a code, then $\boldsymbol{\theta}$ appears inside $f$ at that fixed code and backprop handles it. The encoder's weights $\boldsymbol{\phi}$ are different. They do not appear in $f$; they set the distribution the codes are drawn from. Changing $\boldsymbol{\phi}$ changes which codes come out, and a sampling operation has no derivative with respect to the parameters of its distribution, so backprop stops there.

One way around this is the score-function estimator (also called REINFORCE or the likelihood-ratio estimator), which moves the gradient onto the log-density:

Its sample average is an unbiased estimate of the gradient, and it works even for discrete codes. The paper dismisses it in one sentence as having “very high variance”, citing Blei, Jordan and Paisley (2012). The form shows why. Each sample multiplies the value of $f$ by a score term, so the estimate's spread scales with the size of $f$. Adding a constant to $f$ leaves the true gradient unchanged but makes this estimator noisier. It uses $f$ only through correlation: codes that happened to score well become slightly more likely.

The paper's alternative is the reparameterization trick. Its equation (4) writes the random code as a deterministic, differentiable function of the encoder's parameters and a separate noise variable whose distribution has no parameters:

For the diagonal Gaussian encoder, the instance of (4) the VAE uses is

where $\odot$ multiplies elementwise. The randomness now sits entirely in $\boldsymbol{\varepsilon}$, which does not depend on $\boldsymbol{\phi}$, and $\mathbf{z}$ is an ordinary formula of $\boldsymbol{\mu}$ and $\boldsymbol{\sigma}$. Expectations over $q$ become expectations over the fixed noise, which the paper's equation (5) estimates with $L$ samples:

The right side is differentiable in $\boldsymbol{\phi}$, so its gradient is the gradient of each term, computed by backprop. The first equality is a change of variables: the probability mass of a small region $q_{\boldsymbol{\phi}}(\mathbf{z}\mid\mathbf{x})\,d\mathbf{z}$ equals the mass $p(\boldsymbol{\epsilon})\,d\boldsymbol{\epsilon}$ of the noise region that maps onto it.

A one-coordinate example of what receives gradient: the encoder outputs $\mu = 0.5$ and $\log\sigma^2 = -1.386$ (so $\sigma = 0.5$), the noise draw is $\varepsilon = 1.2$, and the code is $z = 0.5 + 0.5 \times 1.2 = 1.1$. The decoder produces $f(1.1)$ and its slope $f'(1.1)$. By the chain rule the gradient reaching $\mu$ is $f'(1.1) \times 1$ and the gradient reaching $\sigma$ is $f'(1.1) \times 1.2$, and from there both flow back into $\boldsymbol{\phi}$. The noise $\varepsilon$ is an input with no parameters behind it, so no gradient flows into it. Figure 7 draws the two wirings.

Section 2.4 of the paper lists three families of distributions that can be written this way. Distributions with a tractable inverse CDF (exponential, Cauchy, logistic, Gumbel and others) use $\boldsymbol{\epsilon} \sim \mathcal{U}(0, 1)$ and apply the inverse CDF. Location-scale families (Gaussian, Laplace, Student's t, uniform and others) use location + scale × a standard draw, as above. Compositions build on those: a log-normal is the exponential of a Gaussian, a gamma is a sum of exponentials, a Dirichlet is a normalized set of gammas. All of them need a continuous $\mathbf{z}$.

The pathwise gradient uses the slope of $f$ at the drawn code, which keeps its variance low when $f$ is smooth. Figure 8 measures it on a one-dimensional case: $z \sim \mathcal{N}(\mu, \sigma^2)$ with $\mu = 0.6$, $\sigma = 1$, and $f(z) = (z - a)^2$ with $a = -0.9$. The true gradient is $\partial_\mu \mathbb{E}[f] = 2(\mu - a) = 3.0$. From one draw $\varepsilon$, with $z = \mu + \sigma\varepsilon$, the two estimates are

Both average to 3.0. Their standard deviations are 7.2 and 2.0 (exactly $\sqrt{d^4 + 14 d^2 + 15}$ with $d = \mu - a = 1.5$, against $2\sigma$), so the score-function estimate has about 13 times the variance and needs about 13 times as many samples to reach the same precision.

Both estimators target the same gradient, so the reparameterization changes the variance, not the answer. Lower variance is also a strong tendency rather than a theorem: for an $f$ that varies very sharply with $\mathbf{z}$, the pathwise estimate can be the noisier one. For the smooth log-likelihoods of a VAE decoder its variance is far lower (12.9 times lower in the toy of Figure 8), and that lets the paper train with a single code per datapoint.

The full algorithm, on MNIST

Applying (5) to the bound gives the paper's SGVB (Stochastic Gradient Variational Bayes) estimators. The generic one, $\tilde{\mathcal{L}}^{A}$ in equation (6), samples the whole integrand of (2):

with $\mathbf{z}^{(i,l)} = g_{\boldsymbol{\phi}}(\boldsymbol{\epsilon}^{(i,l)}, \mathbf{x}^{(i)})$. It works for any prior and encoder that can be reparameterized. When the KL has a closed form, the second estimator $\tilde{\mathcal{L}}^{B}$, equation (7), computes it exactly and samples only the reconstruction:

The paper says $\tilde{\mathcal{L}}^{B}$ “typically has less variance than the generic estimator”, since an exact term contributes no sampling noise. For the Gaussian encoder, substituting the Appendix B formula gives the paper's equation (10), the quantity a VAE maximizes per datapoint:

where $\mathbf{z}^{(i,l)} = \boldsymbol{\mu}^{(i)} + \boldsymbol{\sigma}^{(i)}\odot\boldsymbol{\varepsilon}^{(l)}$ and $\boldsymbol{\varepsilon}^{(l)} \sim \mathcal{N}(\mathbf{0},\mathbf{I})$. For a dataset of $N$ points, equation (8) estimates the full bound from a random minibatch of $M$:

The paper's Algorithm 1, called AEVB (Auto-Encoding Variational Bayes), repeats three steps: draw a minibatch, draw fresh noise, and step both $\boldsymbol{\theta}$ and $\boldsymbol{\phi}$ along the gradient of (8). The experiments use $M = 100$ and $L = 1$; the paper found one sample per datapoint enough as long as the minibatch is large enough, giving $M = 100$ as the example. The 100 independent draws in a minibatch average out much of the per-sample noise. Here is one step for binarized MNIST in the minimized-loss convention, with every shape written out:

# one AEVB step on a minibatch X of shape (M, 784), binary pixels, M = 100
# prior p(z) = N(0, I), J = 20 latent dims, 500 hidden units; minimizes -L
h      = tanh(X @ We + be)                  # (M, 500)  encoder hidden layer
mu     = h @ Wm + bm                        # (M, 20)
logvar = h @ Wv + bv                        # (M, 20)   log sigma^2 (App. C.2)
eps    = randn(M, 20)                       # outside noise, no parameters
z      = mu + exp(0.5 * logvar) * eps       # reparameterization, eq (4)
y      = sigmoid(tanh(z @ Wd + bd) @ Wo + bo)   # (M, 784) pixel probs, eq (11)
recon  = sum(X * log(y) + (1 - X) * log(1 - y), axis=1)       # (M,)
kl     = 0.5 * sum(mu**2 + exp(logvar) - logvar - 1, axis=1)  # (M,)  +KL
loss   = mean(kl - recon)                   # -L per datapoint
loss.backward()                             # gradients for encoder AND decoder
adagrad_step(params, lr=0.02)               # paper: Adagrad, lr in {0.01, 0.02, 0.1}

The shapes follow the paper's Appendix C: one hidden layer of 500 tanh units in each network, a sigmoid output for the Bernoulli decoder, and an encoder that outputs $\log\boldsymbol{\sigma}^2$ (the paper's equation (12) with $\mathbf{x}$ and $\mathbf{z}$ swapped), from which $\boldsymbol{\sigma} = \exp(\tfrac{1}{2}\log\boldsymbol{\sigma}^2)$. Outputting the log keeps the variance positive whatever the network produces. The encoder has 412,540 parameters ($784 \times 500$ weights and 500 biases, then two heads of $500 \times 20 + 20$), and the decoder 403,284 ($20 \times 500 + 500$, then $500 \times 784 + 784$), 815,824 in total. The listing uses a mean over the minibatch where (8) uses a sum scaled by $N/M$; the two differ by a constant factor that the learning rate absorbs. The paper also adds a small weight-decay term, which the listing leaves out.

The reconstruction line is the Bernoulli log-likelihood of the paper's equation (11), summed over the 784 pixels:

where $\mathbf{y}$ is the decoder's vector of pixel probabilities. A pixel that is on ($x_i = 1$) with $y_i = 0.9$ contributes $\log 0.9 = -0.105$ nats; a pixel that is off with the same $y_i$ contributes $\log 0.1 = -2.30$. This is binary cross-entropy with the sign flipped. For continuous data such as the Frey Face images, the decoder is instead a Gaussian with its own mean and log-variance heads, equation (12), and the reconstruction is a Gaussian log-density. The paper constrains that decoder's mean to $(0,1)$ with a sigmoid.

The single loss.backward() reaches the decoder through the reconstruction term and the encoder through both terms: the KL directly, and the reconstruction through the reparameterized z. After training, generation uses only the decoder and the prior:

# generation: no encoder and no data
z = randn(n, 20)                            # z ~ p(z) = N(0, I)
y = sigmoid(tanh(z @ Wd + bd) @ Wo + bo)    # (n, 784) pixel probabilities
x = bernoulli(y)                            # n new digits (or display y itself)

Walking the latent space

Because every datapoint's Gaussian sits near the origin and overlaps its neighbours, the decoder is trained on codes spread across the whole bulk of the prior, and nearby codes decode to similar outputs. The paper shows this with a two-dimensional latent space (its Figure 4). It takes a grid of codes, decodes each one, and lays the images out in grid order.

Figure 9 is a synthetic stand-in: its “decoder” blends four digit shapes across the plane, so it shows the layout and the grid recipe, not a trained model. The paper does not space its grid evenly in $\mathbf{z}$. It spaces the points evenly on the unit square and passes each coordinate through the inverse Gaussian CDF, so every row and column of the picture covers the same share of the prior's mass. With 9 columns spread evenly over $[-2.3, 2.3]$, the centre column holds 20.2% of the prior and each edge column 2.6%; with the inverse-CDF grid every column holds 11.1%.

In the inverse-CDF grid the ticks for $z = \pm 1$ sit at 16% and 84% of the way across, and $z = \pm 2$ sits near the edges, so most of the picture is spent on the region where prior draws actually land. Walking a straight line between two codes then gives a gradual change in the decoded output, which made VAE latent spaces a common tool for interpolation and visualization.

Experiments and limits

The experiments are small by current standards. The datasets are MNIST and the Frey Face video frames, and the comparison measures the variational lower bound itself (averaged per datapoint, in nats, since the paper uses natural logarithms) plus an estimate of the marginal likelihood for a small latent space. There is no FID, no Inception Score, and no CIFAR or ImageNet; FID dates from 2017.

The settings: minibatches of $M = 100$ with $L = 1$ sample each; all weights initialized from $\mathcal{N}(0, 0.01)$; 500 hidden units for MNIST and 200 for Frey Face, fewer for Frey Face because it is a much smaller dataset and more units overfit. Updates are stochastic gradient ascent on the bound with step sizes adapted by Adagrad, the global step size chosen from {0.01, 0.02, 0.1} by training performance in the first few iterations. Adam, which Kingma published with Ba the following year, did not exist yet. The objective also includes a small weight decay, equivalent to a prior $p(\boldsymbol{\theta}) = \mathcal{N}(\mathbf{0}, \mathbf{I})$ on the weights, which makes the training an approximate MAP estimate. The Figure 2 caption reports 20 to 40 minutes of computation per million training samples on an Intel Xeon CPU running at an effective 40 GFLOPS.

The paper compares against two baselines. The wake-sleep algorithm (Hinton et al., 1995) also trains an encoder and a decoder, in alternating phases; the paper calls it the only other on-line method for the same class of continuous latent-variable models, and notes that its two objectives do not together optimize a bound on the marginal likelihood. Monte Carlo EM replaces the E-step's exact posterior with samples from Hybrid Monte Carlo, a gradient-based sampler; the paper runs 10 leapfrog steps followed by 5 weight updates per sample.

On the lower bound, for every latent size tried, AEVB “converged considerably faster and reached a better solution” than wake-sleep on both datasets. On the estimated marginal likelihood, using 100 hidden units and 3 latent dimensions (beyond that the estimates became unreliable), AEVB is compared with both wake-sleep and Monte Carlo EM on a small and a large MNIST training set. Monte Carlo EM is not an on-line algorithm and cannot be applied efficiently to the full dataset, while AEVB and wake-sleep can. Adding latent dimensions did not cause overfitting, which the authors attribute to the regularizing effect of the bound; a dimension the decoder does not use costs KL, so the optimum sets it to the prior. Random samples from models with 2, 5, 10 and 20 latent dimensions appear in the paper's Figure 5.

The limits became their own research topics. The Gaussian shape and the shared encoder leave a gap between the bound and $\log p(\mathbf{x})$, as in Figures 2 and 5; importance-weighted bounds (IWAE, 2015) and richer encoders attack it. The balance between the two terms of (3) is fixed at $\beta = 1$; the $\beta$-VAE (Higgins et al., 2017) raises the weight on the KL term, usually to $\beta > 1$, to get latent dimensions that each track one factor of variation, at a cost in reconstruction. Latent dimensions can collapse onto the prior, as at the right end of Figure 4. The reparameterization needs a continuous code, so discrete latents waited for relaxations such as Gumbel-Softmax (2016). And the usual way to display a sample is the decoder's mean, which averages over the pixel noise the model assigns, part of why VAE samples look soft next to those of a GAN.

The core recipe is still in wide use: a bound in place of an intractable likelihood, an encoder in place of per-datapoint inference, and a reparameterized draw so that one backward pass trains both networks. The agent in World Models compresses frames with a VAE, DDPM derives its training objective from a variational bound of the same kind, and the autoencoder under latent diffusion (Stable Diffusion) carries a lightly weighted version of this page's KL term.