Auto-Encoding Variational Bayes, annotated

How to read this page

Nothing on this page assumes you already know the jargon. Three things help:

• Any dotted word explains itself when you hover it, tab to it, or tap it.
• Every symbol inside an equation does the same, and each equation is followed by a table of its symbols.
• The diagrams are live: hover or tap any part to see what it does. The sliders recompute everything as you drag.

Each idea is explained in the same order: an everyday picture, a tiny example you could check by hand, a diagram, then the math, and finally why it still matters. The whole page leans on one tiny world, small enough to solve exactly, so you can see every quantity the paper can only approximate:

The code that builds a VAE from scratch, in NumPy, lives in the autoencoders lesson.

Abstract

“How can we perform efficient inference and learning in directed probabilistic models, in the presence of continuous latent variables with intractable posterior distributions, and large datasets?”Kingma and Welling (2013), Abstract. Read the original

Everyday picture

A detective sees a footprint and wants to know who made it. Going forwards is easy: pick a person, and you can say what footprint they would leave. Going backwards (footprint to person) is the hard direction. This paper trains two networks together: a forward one that turns hidden causes into data (a decoder) and a backward one that makes a quick educated guess of the causes from the data (an encoder). The guesser is trained so that its guesses make the forward model as good as possible.

What the paper claims

• A reparameterization of the variational lower bound gives an estimate of it that ordinary stochastic gradient methods can optimise directly.
• For a dataset of independent examples, each with its own hidden code, training an approximate inference model (a network that guesses the code from the example) with that estimate makes inference fast. With neural networks on both sides, this is the variational autoencoder.
• The experiments (on handwritten digits and faces) show it trains faster and reaches a better bound than the older wake-sleep algorithm.

Why it matters today

The reparameterization trick, z = μ + σ·ε, is now a standard tool for sending gradients through randomness. VAEs themselves became the compressors inside image generators: the autoencoder at both ends of latent diffusion is a VAE with a small KL weight.

1 Introduction · original

“When a neural network is used for the recognition model, we arrive at the variational auto-encoder.”Kingma and Welling (2013), §1

Everyday picture

Before this paper, the standard way to guess hidden causes (variational inference) worked like solving a fresh algebra problem for every single photo: slow, and only possible when the algebra happened to come out neatly. The paper's move is to replace all those individual solutions with one trained guesser that answers for any photo in a single pass. The paper calls it a recognition model; today everyone calls it the encoder.

Tiny example

In the one-number world, the exact backward answer for x = 2 is a bell curve centred at 1 with variance 0.5. For x = 4 it would be centred at 2. A recognition model learns the rule “centre = x / 2, variance = 0.5” once, then applies it to every x for free. In a real model no such neat rule exists, so a neural network learns an approximate one.

Two names the paper introduces

• SGVB (Stochastic Gradient Variational Bayes): the estimator of the bound, built with the reparameterization trick (§2.3).
• AEVB (Auto-Encoding Variational Bayes): the training algorithm that uses SGVB to fit an encoder and a decoder together (Algorithm 1). With neural networks, AEVB trains a VAE.

Why it matters

Replacing a per-example calculation with a trained network is called amortising the cost, and the idea spread far beyond VAEs. It also avoids slow sampling loops such as MCMC for every example, which is what made these models trainable on large datasets.

2 Method · original

Everyday picture

Think of a story about how every example was made, told in two steps. First, nature rolls some hidden dice: that is the code z. Then it turns the code into the thing you see: that is x. The paper draws this story as a small graphical model, and adds a second, dashed story going backwards: a guess of the dice from what you see.

Reading it: each circle is a quantity. z (open) is hidden and x (shaded) is observed. Solid arrows are the generative story: θ sets how z is drawn and how z becomes x. The dashed arrows are the guess going the other way: φ sets how the encoder reads x and proposes values of z. The rounded box with N in its corner is a plate: everything inside happens once per example, while θ and φ, outside it, are shared by all N examples. The paper's key point is that θ and φ are learned together.

2.1 Problem scenario · original

Everyday picture

Suppose a factory makes cakes from a secret recipe card. You only ever see cakes. To judge how good your own guess at the factory's process is, you would like to know: “how likely is this exact cake, averaged over every recipe card the factory might have drawn?” That average over all hidden causes is the marginal likelihood. For a real model it means summing over every possible code, which is hopeless.

Tiny example

In the one-number world, the chance of seeing x = 2 is the average, over every possible z, of “how likely is z” times “how likely is x = 2 given that z”. Chop the z line into steps of 0.01 from −8 to 8, add up the 1,601 products times the step width, and you get 0.1038. That is fine for one number. With a 20-number code and just 100 grid points per number you would need 100²⁰ = 10⁴⁰ evaluations of a neural network, for each example.

In words: “the probability the model gives to an example is the chance of each possible code times the chance that code produces the example, added up over every code.”

With the numbers: p(x = 2) = ∫ N(z; 0, 1) · N(2; z, 1) dz = 0.1038, which is exactly the bell curve of spread √2 evaluated at 2 (adding two independent spreads of 1 gives variance 2).

In Python:

import math
def N(x, mean, var):
    # the bell-curve density at x
    return math.exp(-(x - mean) ** 2 / (2 * var)) / math.sqrt(2 * math.pi * var)
x, dz = 2.0, 0.01
zs = [-8 + i * dz for i in range(1601)]
# ∫ p(z) p(x | z) dz, as a sum of thin slices
p_x = sum(N(z, 0, 1) * N(x, z, 1) * dz for z in zs)
round(p_x, 4)  # → 0.1038
# the exact answer: x is the sum of two spreads of 1, so its variance is 2
round(N(x, 0, 2), 4)  # → 0.1038
# log p(x), used on the rest of the page
round(math.log(p_x), 4)  # → -2.2655

The backward question is answered by Bayes' rule, and the marginal likelihood sits in its denominator, which is why the backward question is just as hopeless:

In words: “how likely a code is, now that we have seen the example, equals how well that code explains the example, times how likely the code was beforehand, divided by how likely the example is overall.”

With the numbers: at z = 1, p(z = 1 | x = 2) = N(2; 1, 1) · N(1; 0, 1) / 0.1038 = 0.2420 × 0.2420 / 0.1038 = 0.5642. That is the height at its centre of the bell curve with mean 1 and variance 0.5: the exact posterior of this world.

In Python:

import math
def N(x, mean, var):
    return math.exp(-(x - mean) ** 2 / (2 * var)) / math.sqrt(2 * math.pi * var)
x, z = 2.0, 1.0
p_x = N(x, 0, 2)
# p(x | z) and p(z)
round(N(x, z, 1), 4), round(N(z, 0, 1), 4)  # → (0.242, 0.242)
# Bayes' rule
round(N(x, z, 1) * N(z, 0, 1) / p_x, 4)  # → 0.5642
# the posterior is the bell curve N(1, 0.5): same height at z = 1
round(N(z, 1, 0.5), 4)  # → 0.5642

What the paper wants from one method

1. Good values of θ (maximum likelihood or MAP), so the model can generate data that looks real.
2. Fast approximate inference of z from x, useful as a compact code.
3. Approximate marginal inference of x, useful for denoising, inpainting and super-resolution.

And it wants these when the integral is intractable and the dataset is too large for anything but small random batches. Its answer is a recognition model q_φ(z | x), which the paper also calls a probabilistic encoder; p_θ(x | z) is the probabilistic decoder.

Why it matters today

This two-step story, “draw a simple code, then decode it with a network”, is the template for most generative models since. A diffusion model is the same story told over a thousand small steps; see the DDPM companion.

2.2 The variational bound · original

“Since this KL-divergence is non-negative, the second RHS term is called the (variational) lower bound on the marginal likelihood of datapoint i”Kingma and Welling (2013), §2.2

Everyday picture

You cannot weigh a ship directly, but you can weigh everything you loaded onto it and know the ship weighs at least that much. Push that lower figure up and you are pushing towards the true weight. The evidence lower bound (ELBO) is that lower figure for log p(x). The shortfall is exactly how far the encoder's guess q is from the true posterior, measured by the KL divergence.

Tiny example

In the one-number world, log p(x = 2) = −2.2655. Let the encoder guess q = the standard bell curve N(0, 1), a lazy guess that ignores x. Its bound is −3.4189, short by 1.1534. Let it guess the true posterior N(1, 0.5) instead: the bound is −2.2655, short by exactly 0.

Reading it: the chart shows two bell curves over the hidden number z. The fixed one is the true posterior for x = 2; the one you move is the encoder's guess q. Below it, the bars show sizes, so longer means more negative. The grey bar is log p(x) = −2.2655, which never moves. The solid bar is the ELBO, always at least as long, and the striped bar is the gap: exactly the difference in length between the two, so ELBO + gap = log p(x) whatever you choose. Drag the centre towards 1 and the variance towards 0.5: the striped gap shrinks to nothing and the ELBO rises to meet log p(x). Everything here is computed exactly from the formulas below; nothing is sampled.

The math

In words: “the log-probability of an example splits into two pieces: how far the encoder's guess is from the true posterior, plus the bound. The first piece is never negative, so the bound is never above the log-probability.”

With the numbers: for the lazy guess N(0, 1): −2.2655 = 1.1534 + (−3.4189). For the perfect guess N(1, 0.5): −2.2655 = 0 + (−2.2655).

In Python:

import math
x = 2.0
log_p_x = -0.5 * math.log(2 * math.pi * 2) - x ** 2 / 4
def elbo(m, v):
    # E_q[log p(x | z)] minus KL(q ‖ p(z)), both in closed form
    reconstruction = -0.5 * math.log(2 * math.pi) - ((x - m) ** 2 + v) / 2
    kl_to_prior = 0.5 * (m ** 2 + v - math.log(v) - 1)
    return reconstruction - kl_to_prior
def kl_to_posterior(m, v):
    # KL(N(m, v) ‖ N(1, 0.5)), the gap
    return 0.5 * (v / 0.5 + (m - 1) ** 2 / 0.5 - 1 - math.log(v / 0.5))
round(log_p_x, 4)  # → -2.2655
# the lazy guess q = N(0, 1): gap plus bound
round(kl_to_posterior(0, 1), 4), round(elbo(0, 1), 4)  # → (1.1534, -3.4189)
round(kl_to_posterior(0, 1) + elbo(0, 1), 4)  # → -2.2655
# the perfect guess q = N(1, 0.5): no gap
round(kl_to_posterior(1, 0.5), 4), round(elbo(1, 0.5), 4)  # → (0.0, -2.2655)

The paper then rewrites the bound in the form every VAE trains on:

In words: “the bound is how well codes drawn from the guess rebuild the example, on average, minus how far the guess strays from the prior.”

With the numbers: for the perfect guess N(1, 0.5): the rebuild term is −0.9189 − (1² + 0.5)/2 = −1.6689, the KL to the prior is ½(1 + 0.5 − log 0.5 − 1) = 0.5966, and −0.5966 + (−1.6689) = −2.2655.

In Python:

import math
x, m, v = 2.0, 1.0, 0.5
# E_q[log p(x | z)]: the average of −½ log 2π − (x − z)² / 2 over z ~ N(m, v)
reconstruction = -0.5 * math.log(2 * math.pi) - ((x - m) ** 2 + v) / 2
round(reconstruction, 4)  # → -1.6689
# D_KL(q ‖ p(z)) for q = N(m, v) and p(z) = N(0, 1)
kl = 0.5 * (m ** 2 + v - math.log(v) - 1)
round(kl, 4)  # → 0.5966
round(-kl + reconstruction, 4)  # → -2.2655

Why not just sample the gradient?

To train the encoder you need the gradient of the bound with respect to φ, and φ shapes the very distribution you are averaging over. The textbook trick (the score-function estimator) works, but the paper notes that it “exhibits very high variance” and is impractical here. Section 2.4 shows the fix, with a live comparison.

Why it matters today

The same bound, written the same way, trains every VAE. The autoencoders lesson uses its negative as the loss, rebuild error plus β times KL, and shows what happens as β varies. The diffusion paper writes an ELBO of exactly this kind over a thousand steps.

2.3 The SGVB estimator and the AEVB algorithm · original

Everyday picture

To estimate the average height in a city you do not measure everyone: you measure a random handful and scale up. SGVB estimates the bound the same way: the KL part has an exact formula, and the rebuild part is estimated from a few random codes, often just one. Then AEVB repeats: grab a small random batch of examples, estimate, nudge both networks, repeat.

Tiny example

Take the perfect guess N(1, 0.5) and draw one noise number, ε = 0.5. The code is z = 1 + √0.5 × 0.5 = 1.3536. Its rebuild score is log N(2; 1.3536, 1) = −1.1279. The estimate of the bound is −0.5966 + (−1.1279) = −1.7245. One draw is noisy: the true bound is −2.2655. Averaged over 100,000 draws, the estimate comes to −2.26.

In words: “keep the exact KL, and replace the average rebuild score with the average over L codes, each made by pushing a fresh noise draw through the encoder's recipe.”

With the numbers: L = 1 and ε = 0.5: −0.5966 + log N(2; 1.3536, 1) = −0.5966 − 1.1279 = −1.7245.

In Python:

import math, random
x, m, v = 2.0, 1.0, 0.5
kl = 0.5 * (m ** 2 + v - math.log(v) - 1)
def log_p_x_given_z(z):
    return -0.5 * math.log(2 * math.pi) - (x - z) ** 2 / 2
# z = g_φ(ε, x) = m + √v · ε, with one noise draw (L = 1)
z = m + math.sqrt(v) * 0.5
round(z, 4)  # → 1.3536
round(-kl + log_p_x_given_z(z), 4)  # → -1.7245
# (1/L) Σ over many draws: the estimate is right on average
random.seed(0)
L = 100_000
round(-kl + sum(log_p_x_given_z(m + math.sqrt(v) * random.gauss(0, 1)) for _ in range(L)) / L, 2)  # → -2.26

The paper gives two versions. Estimator A (its equation 6) samples the whole bound; estimator B (above, equation 7) computes the KL exactly and so “typically has less variance”. For a whole dataset of N examples, a random minibatch of M is scaled up:

In words: “estimate the bound for the whole dataset by adding up the estimates for a random batch and multiplying by how many batches would fit in the dataset.”

With the numbers (illustrative: a dataset of N = 1,000 examples, all equal to our x = 2, and the paper's batch size M = 100): if the batch's per-example estimates average −2.26, the batch sum is −226 and the dataset estimate is 10 × (−226) = −2,260, close to the true 1,000 × (−2.2655) = −2,265.5.

In Python:

N, M = 1_000, 100
batch = [-2.26] * M
# (N / M) Σᵢ L̃(x⁽ⁱ⁾)
round(N / M * sum(batch), 1)  # → -2260.0
# the true total, 1,000 × log p(x)
round(N * -2.2655, 1)  # → -2265.5

Reading it: read from the top. After one initialisation, the loop on the left repeats four steps: pick a random minibatch, draw fresh noise for every example, differentiate the minibatch estimate of the bound with respect to both networks' weights, and take a step. There is no inner loop per example, no sampling chain and no separate phase for the encoder. That single, ordinary training loop is the practical contribution: the paper uses M = 100 and L = 1.

The connection to autoencoders

Read estimator B again. Its second term asks: encode x to a code, decode the code, and score how well x comes back. That is an autoencoder's rebuild error, with a minus sign. The first term, the KL, acts as a regulariser that the bound itself dictates, rather than one added by hand.

In code: VAE.loss_and_gradients computes both terms for a batch with one noise draw per example, and train is Algorithm 1 with Adam as the update.

Why it matters today

“One sample per example, a big enough batch, and backpropagation” is still how VAEs are trained. The whole method fits in a standard training loop, which is why it spread so fast.

2.4 The reparameterization trick · original

“The essential parameterization trick is quite simple.”Kingma and Welling (2013), §2.4

Everyday picture

Two ways to pick a random seat in a row. You can close your eyes and point: if someone slides the row one seat left, there is no saying how your pick changes. Or you can roll a die for “how many seats from the middle” and count from wherever the middle now is: slide the row, and your seat slides with it by exactly the same amount. Same randomness, but now the result responds smoothly to the row's position. That is the trick: roll the dice first, as an ordinary input, then build the random code with arithmetic the gradient can pass through.

Tiny example

The encoder says μ = 1 and σ = 1. The die gives ε = 0.5, so z = 1 + 1 × 0.5 = 1.5. Nudge μ to 1.001 and z becomes 1.501: z moved by exactly the nudge, so ∂z/∂μ = 1. Nudge σ instead and z moves by ε times the nudge, so ∂z/∂σ = ε = 0.5.

In words: “averaging f over codes drawn from a bell curve with centre μ and spread σ is the same as averaging f(μ + σε) over standard noise ε, and a handful of noise draws estimates it.”

With the numbers: take f(z) = z² with μ = 1 and σ = 1. The exact average is μ² + σ² = 2. The single draw ε = 0.5 gives f(1.5) = 2.25; 100,000 draws give 2.01. Its gradient with respect to μ is 2μ = 2, and each draw's own gradient, 2(μ + σε), averages to 2 as well.

In Python:

import random
mu, sigma = 1.0, 1.0
def f(z):
    return z ** 2
# one draw, ε = 0.5
f(mu + sigma * 0.5)  # → 2.25
# (1/L) Σ f(μ + σ ε⁽ˡ⁾): many draws approach μ² + σ² = 2
random.seed(0)
eps = [random.gauss(0, 1) for _ in range(100_000)]
round(sum(f(mu + sigma * e) for e in eps) / len(eps), 2)  # → 2.01
# each draw's gradient with respect to μ is f'(μ + σε) = 2(μ + σε)
round(sum(2 * (mu + sigma * e) for e in eps) / len(eps), 2)  # → 2.01

Try it: two ways to estimate one gradient

Both estimators below aim at the same number, the gradient of the average of z² with respect to μ, which is exactly 2μ. The reparameterized one differentiates through z = μ + σε. The score-function one (the “naïve” estimator of §2.2) multiplies each f(z) by (z − μ)/σ², the slope of the log bell curve. Drag the sliders and watch how fast each running average settles.

In words: “one draw's guess at the gradient is either the slope of f at the drawn code, or f's value at the drawn code times how much the drawn code favours a larger μ.”

With the numbers: with μ = σ = 1 and ε = 0.5 (so z = 1.5): the reparameterized guess is 2 × 1.5 = 3 and the score-function guess is 2.25 × 0.5 = 1.125. Both average to 2 over many draws, but their spreads per draw are 2 and 5.48: the score function is almost three times as noisy, and far worse as σ shrinks (5.48 becomes 10.7 at σ = 0.1, while the reparameterized spread falls to 0.2).

In Python:

import math
mu, sigma, eps = 1.0, 1.0, 0.5
z = mu + sigma * eps
# reparameterized: f'(z) = 2z
2 * z  # → 3.0
# score function: f(z) (z − μ) / σ²
z ** 2 * (z - mu) / sigma ** 2  # → 1.125
# spread per draw, from the bell curve's moments E[ε²] = 1, E[ε⁴] = 3, E[ε⁶] = 15
def spreads(a, s):
    score_var = (a ** 4 + 18 * a ** 2 * s ** 2 + 15 * s ** 4) / s ** 2 - 4 * a ** 2
    return round(2 * s, 2), round(math.sqrt(score_var), 2)
spreads(1.0, 1.0)  # → (2.0, 5.48)
spreads(1.0, 0.1)  # → (0.2, 10.68)

Reading it: the x-axis counts noise draws; each line is the running average of one estimator's guesses, using the same draws. The flat dotted line is the true gradient 2μ. The page opens at μ = 1 and σ = 0.3, a small, confident region like the ones a VAE encoder learns. There the reparameterized line (solid) sits on the target almost from the first draw, while the score-function line (dashed) swings and is still visibly off after hundreds of draws: one draw's guess is about 9 times as noisy. Widen σ towards 1 and the two get closer (spreads 2 and 5.48); shrink it to 0.1 and the gap grows to about 50 times. Lines are clipped to ±3 around the target so a wild early swing cannot squash the picture. The draws come from a fixed seed; “Draw new noise” picks another.

Which distributions can be reparameterized?

The paper lists three recipes for writing a draw as g_φ(ε, x) with ε from a fixed distribution:

• Inverse CDF: draw ε uniformly between 0 and 1 and push it through the inverse CDF. The exponential, Cauchy, logistic and several others work this way.
• Location and scale: take the standard version of the distribution and set g = location + scale · ε. The Gaussian, Laplace, Student's t and others.
• Composition: build the draw from other reparameterizable draws, as a log-normal is e raised to a Gaussian draw.

In words: “to draw a waiting time with rate λ, draw a uniform number and push it through the exponential distribution's inverse CDF.”

With the numbers: λ = 2 and u = 0.5 give z = ln 2 / 2 = 0.3466. The gradient with respect to λ exists: ∂z/∂λ = ln(1 − u)/λ² = −0.1733, so a faster rate means a shorter wait.

In Python:

import math
lam, u = 2.0, 0.5
# z = −ln(1 − u) / λ
z = -math.log(1 - u) / lam
round(z, 4)  # → 0.3466
# ∂z/∂λ = ln(1 − u) / λ²
round(math.log(1 - u) / lam ** 2, 4)  # → -0.1733

Why it matters today

This trick is why VAEs could be trained with plain backpropagation. The same move (keep the noise as an input, make everything else differentiable) shows up in diffusion models, whose noisy training inputs are written exactly as √ᾱ·x + √(1 − ᾱ)·ε, and in reinforcement learning with continuous actions. Discrete choices cannot be reparameterized this way, which is why training with them stays harder. See the reparameterization section of the autoencoders lesson.

In code: reparameterize is z = μ + σ ⊙ ε, and a test checks the two gradients through it against finite differences.

3 Example: the variational auto-encoder · original

Everyday picture

Now put neural networks on both sides. The encoder reads a picture and, instead of pinning it to one exact spot on a map, draws a small fuzzy circle: “somewhere around here”. The decoder is handed a random point from inside that circle and must still redraw the picture. Every circle pays rent (the KL) for sitting far from the middle of the map or for shrinking to a pin, so circles crowd together and overlap, and the whole middle of the map ends up decoding to plausible pictures.

Tiny example

A 2-number code. The encoder says μ = (0.5, 0) and log σ² = (log 0.25, 0), so σ = (0.5, 1). Noise ε = (1.2, −0.4) gives z = (0.5 + 0.5 × 1.2, 0 + 1 × (−0.4)) = (1.1, −0.4). The decoder turns z into pixel probabilities; say, for three pixels, y = (0.9, 0.2, 0.6) (illustrative values), and the real pixels are x = (1, 0, 1).

Reading it: data flows upwards. The encoder reads the picture and has two heads: a centre μ and a size, given as log σ² so that it can be any number. The sampling box turns the fuzzy region into one concrete code, with noise ε coming in from the side as an ordinary input (the trick of §2.4). The decoder turns the code into a chance for every pixel. Two wires reach the top: the long one on the left carries the real picture to be scored against the decoder's chances, and the one on the right carries μ and log σ² straight to the KL, which needs no sampling at all. The bound adds the two.

The math

The paper chooses the prior p_θ(z) = N(0, I), which has no parameters, and a Gaussian guess with a separate spread for each code number:

In words: “the encoder's guess for example i is a bell curve for each code number, with a centre and a spread that the encoder network computes from the example.”

With the numbers: μ = (0.5, 0) and σ² = (0.25, 1). The log-density of the guess at its own centre is Σ_j −½ log(2πσ_j²) = −½ log(2π × 0.25) − ½ log(2π) = −0.2258 − 0.9189 = −1.1447.

In Python:

import math
mu, var = [0.5, 0.0], [0.25, 1.0]
def log_N(z, m, v):
    return -0.5 * math.log(2 * math.pi * v) - (z - m) ** 2 / (2 * v)
# log q(z | x) at z = μ: one term per code number, added (σ² I means independent numbers)
round(sum(log_N(m, m, v) for m, v in zip(mu, var)), 4)  # → -1.1447

With both prior and guess Gaussian, the KL has a closed form (derived in Appendix B), and the whole estimator for one example is:

In words: “for each code number, add one, plus the log of its variance, minus its squared centre, minus its variance, and halve the total (that is minus the KL); then add the rebuild score of a code made by stretching fresh noise by σ and shifting it by μ.”

With the numbers: the first code number gives ½(1 − 1.3863 − 0.25 − 0.25) = −0.4431 and the second ½(1 + 0 − 0 − 1) = 0, so minus the KL is −0.4431. The code is z = (1.1, −0.4), and with the illustrative pixel chances y = (0.9, 0.2, 0.6) the rebuild score is −0.8393 (Appendix C shows how). The estimate is −0.4431 − 0.8393 = −1.2825; a VAE's loss is its negative, 1.2825.

In Python:

import math
mu, var, eps = [0.5, 0.0], [0.25, 1.0], [1.2, -0.4]
# ½ Σ_j (1 + log σ_j² − μ_j² − σ_j²): minus the KL
neg_kl = 0.5 * sum(1 + math.log(v) - m ** 2 - v for m, v in zip(mu, var))
round(neg_kl, 4)  # → -0.4431
# z = μ + σ ⊙ ε
z = [m + math.sqrt(v) * e for m, v, e in zip(mu, var, eps)]
[round(z_j, 2) for z_j in z]  # → [1.1, -0.4]
# the rebuild score log p(x | z) for x = (1, 0, 1) and illustrative y = (0.9, 0.2, 0.6)
x, y = [1, 0, 1], [0.9, 0.2, 0.6]
log_p = sum(x_i * math.log(y_i) + (1 - x_i) * math.log(1 - y_i) for x_i, y_i in zip(x, y))
round(neg_kl + log_p, 4)  # → -1.2825

In code: VAE is this model with a two-headed encoder, VAE.encode_distribution returns μ and log σ², and kl_to_standard_normal is the KL above, summed over the code numbers. The lesson uses squared error for the rebuild score, which is the Gaussian decoder of Appendix C with the noise level folded into a weight β.

Why it matters today

This section is the recipe every VAE still follows: a standard normal prior, a diagonal Gaussian encoder with μ and log σ² heads, one reparameterized sample, and a closed-form KL. Modern image autoencoders swap the small MLPs for deep convolutional networks and weight the KL lightly, but the equation is the same.

4 Related work · original

Everyday picture

Before AEVB, the closest method was the wake-sleep algorithm: two students correcting each other's homework in turns, each with its own marking scheme. AEVB gives both students one shared exam, the bound, so improving either one can only improve the same score.

What the paper connects

• Wake-sleep also trains a recognition model and costs about the same per example, but it optimises two objectives that together are not a bound on the likelihood. It does handle discrete codes, which AEVB does not.
• Linear models: PCA is the maximum-likelihood solution of a linear-Gaussian model with prior N(0, I) and decoder N(x; Wz, εI), in the limit of vanishing noise ε (a result of Roweis, 1998). The VAE is the nonlinear, variational generalisation of that picture.
• Regularised autoencoders (denoising, contractive, sparse) need a hand-tuned penalty to learn useful codes. In the VAE the penalty, the KL, falls out of the bound.
• Rezende, Mohamed and Wierstra (2014) found the same connection, with the same reparameterization, independently.

Tiny example

Why does a plain rebuild score not suffice? An autoencoder whose code is as wide as its input can learn to copy, scoring a perfect rebuild while learning nothing. The paper cites the result that the rebuild criterion alone “is in itself not sufficient for learning useful representations”. The KL term stops the copying: every bit of information stored in the code costs rent.

Why it matters today

The linear case is worth remembering: the autoencoders lesson shows a linear autoencoder learning exactly what PCA learns, then what the bends of a nonlinear network add.

5 Experiments · original

Everyday picture

A race between three ways of training the same kind of model, judged on two scoreboards: the bound itself, and an independent estimate of the true log-likelihood made by a slow but careful sampling method.

The setup

What they found

• The bound (Figure 2): compared with wake-sleep at several code sizes, AEVB “converged considerably faster and reached a better solution in all experiments”.
• No overfitting from extra code numbers: adding superfluous latent variables did not make things worse, which the paper credits to the regularising KL term.
• The marginal likelihood (Figure 3): with a 3-number code and 100 hidden units (where the careful estimate is still reliable), AEVB was compared with wake-sleep and Monte Carlo EM. Monte Carlo EM, which runs a sampling chain per example, cannot be applied efficiently to the full MNIST dataset at all.

Tiny example: why the careful estimate only works in low dimensions

The estimate of log p(x) that judges Figure 3 (Appendix D) fits a density to posterior samples, and the paper notes it is only good below about five code dimensions. The same curse as the grid of §2.1: a grid of 100 points per dimension is 10⁶ points in 3 dimensions but 10²⁰ in 10.

Why it matters today

The benchmarks were small by today's standards, but the claim that mattered held up: one ordinary gradient loop beat bespoke inference schemes. The figures' curves are not redrawn here because the paper publishes them only as plots; read them in the original.

6–7 Conclusion and future work · original

What the paper concludes

SGVB is an estimator of the variational bound that ordinary stochastic gradient methods can optimise; AEVB uses it to learn an inference model for datasets with a continuous code per example.

What it proposed next, and what happened

• Deep, convolutional encoders and decoders: now the norm. The autoencoder inside latent diffusion is a deep convolutional VAE.
• Time-series models: sequential and recurrent VAEs followed.
• SGVB for the global parameters (Appendix F): variational treatment of the weights themselves.
• Supervised models with latent variables: conditional VAEs, which generate given a label.

Why it matters today

The two things to take away are the bound and the trick. Almost every later latent-variable model, from β-VAE to VQ-VAE to diffusion, is either a new bound, a new way to push gradients through randomness, or both.

Appendices A–F · original

The appendices hold the derivations and details the main text points to. A, B and C get the full treatment; D, E and F are summarised.

from statistics import NormalDist
u = [0.1, 0.3, 0.5, 0.7, 0.9]
# Φ⁻¹(u_k): the inverse of the standard normal CDF
[round(NormalDist().inv_cdf(u_k), 4) for u_k in u]  # → [-1.2816, -0.5244, 0.0, 0.5244, 1.2816]
# check: 90% of the mass lies below 1.2816
round(NormalDist().cdf(1.2816), 4)  # → 0.9

import math
def N(z, m, v):
    return math.exp(-(z - m) ** 2 / (2 * v)) / math.sqrt(2 * math.pi * v)
mu, var, dz = 0.5, 0.25, 0.001
zs = [-10 + i * dz for i in range(20_001)]
# ∫ q log p dz and ∫ q log q dz, as sums of thin slices
q_log_p = sum(N(z, mu, var) * math.log(N(z, 0, 1)) * dz for z in zs)
q_log_q = sum(N(z, mu, var) * math.log(N(z, mu, var)) * dz for z in zs)
round(q_log_p, 4), round(q_log_q, 4)  # → (-1.1689, -0.7258)
round(q_log_p - q_log_q, 4)  # → -0.4431
# the closed form ½ Σ (1 + log σ² − μ² − σ²)
round(0.5 * (1 + math.log(var) - mu ** 2 - var), 4)  # → -0.4431

import math
x, y = [1, 0, 1], [0.9, 0.2, 0.6]
# Σ x_i log y_i + (1 − x_i) log(1 − y_i)
round(sum(x_i * math.log(y_i) + (1 - x_i) * math.log(1 - y_i) for x_i, y_i in zip(x, y)), 4)  # → -0.8393
# the same as the log of the picture's probability
round(math.log(0.9 * 0.8 * 0.6), 4)  # → -0.8393
# one pixel's y from a one-number code, with illustrative weights
def sigmoid(a):
    return 1 / (1 + math.exp(-a))
round(sigmoid(3 * math.tanh(1.0 * 1.1 + 0) - 0.2), 2)  # → 0.9

What changed since 2013

The bound and the trick are unchanged. Almost everything around them has moved:

Glossary

Every term with hover guidance on this page, in one place.