The Evidence Lower Bound as a Latent Variable Objective

A latent variable model describes data through a two-stage generative story: draw a code z from a prior p(z), then draw the observation x from a conditional p(x | z) parameterised by θ. The parameters are what training is supposed to recover. The obstacle is that the quantity training needs, p(x) = ∫ p(x | z) p(z) dz, is an integral over the entire latent space, and for any conditional more expressive than a linear-Gaussian pair it has no closed form. The integral is also over a space whose dimension is a design choice, commonly between 8 and 512 for image work, so quadrature is not an option either.

The evidence lower bound is the standard resolution. It converts an intractable integral into an expectation that can be estimated from a single sample and differentiated. It is worth being precise about what the bound does and does not guarantee, because almost every pathology of latent variable training follows from the gap between the bound and the quantity it bounds.

Why the marginal likelihood cannot be evaluated directly

Write the joint as p(x, z) = p(z) p(x | z) and suppose the prior is a standard normal in d dimensions. If the conditional is a Gaussian whose mean and variance come from a neural network, the marginal is a mixture of infinitely many Gaussians with continuously varying weights. Mixture models are tractable only when the mixing distribution is finite and discrete, so this integral has no elementary antiderivative.

Monte Carlo estimation of log p(x) is also unusable as a training signal. The quantity (1/K) Σ p(x | z_k) is an unbiased estimator of p(x) but a biased, high-variance estimator of its logarithm, and the bias grows worse exactly when the model is well fit and the samples concentrate. Jensen’s inequality puts the log of the average above the average of the logs, so the estimator overstates the log likelihood on average, and the overstatement is largest where the conditional is sharply peaked.

The way out is to stop asking for the marginal and ask instead for a distribution over z that explains a given x. That is inference, and in the exact posterior p(z | x) = p(x, z) / p(x) the intractable integral reappears in the denominator. Variational inference replaces the exact posterior with a tractable approximation q(z | x) parameterised by φ, and turns inference into optimisation.

Splitting the log evidence into two computable pieces

For any conditional distribution q(z | x) with support wherever p(z | x) has support, the identity

log p(x) = E_q[ log( p(x, z) / q(z | x) ) ] + KL( q(z | x) || p(z | x) )

holds exactly. The second term is non-negative, so the first term is a lower bound on the log evidence. Expanding the first term by splitting the joint gives the form that is actually computed:

L(θ, φ; x) = E_q[ log p(x | z) ] − KL( q(z | x) || p(z) )

The two terms have clear jobs. The first rewards codes that let the decoder reconstruct x; it is a negative reconstruction loss in expectation. The second penalises any divergence between the encoder’s distribution and the prior, which is what makes sampling from the prior at generation time meaningful. Without the KL term the encoder could use an arbitrary codebook and prior samples would generate nothing useful. With the KL term alone the model would ignore x entirely and drive the divergence to zero at the cost of a useless decoder.

The bound is tight when q(z | x) equals the true posterior. It is loose otherwise, and the looseness is the KL between approximation and truth, which is not observable during training because the true posterior is the thing that could not be computed. A model can therefore have a high bound and a poor marginal likelihood, and the gap is largest where the posterior is multimodal or strongly correlated and the approximation is a diagonal Gaussian.

The Gaussian case in closed form

The usual choice is a diagonal Gaussian q against a standard normal prior. The KL between two diagonal Gaussians, with the prior at zero mean and unit variance, decomposes into one independent summand per latent dimension:

KL = ½ Σ ( μ_i² + σ_i² − log σ_i² − 1 )

Each summand is non-negative, with equality at μ_i = 0 and σ_i = 1. This is why the term behaves as pressure toward the prior: it is a sum of d independent penalties, so widening the latent space increases the total penalty linearly unless the encoder commits harder to the prior. A 256-dimensional latent with every μ_i at 0.5 and every σ_i at 1 already pays roughly 32 nats per example, which is comparable to the reconstruction cost of a moderately detailed image under a unit-variance Gaussian decoder.

The reconstruction term needs a distributional assumption on x. For continuous data a fixed-variance Gaussian gives a log density equal to the negative squared error divided by twice the variance, plus a constant. The constant matters: the relative weight of the two terms is set by the variance, not chosen freely, and fixing it at 1 silently decides how many nats of reconstruction each nat of KL is worth. For binary or quantised data a Bernoulli or categorical decoder gives a cross-entropy term instead, whose scale is fixed by the data range rather than by a free parameter.

import numpy as np

def gaussian_kl(mu, log_var):
    # KL between N(mu, diag(exp(log_var))) and N(0, I), one value per example
    return 0.5 * np.sum(np.exp(log_var) + mu ** 2 - log_var - 1.0, axis=1)

def elbo_per_example(x, x_hat, mu, log_var, sigma2=1.0):
    # Gaussian decoder with fixed variance; both terms are in nats
    recon = -0.5 * np.sum((x - x_hat) ** 2, axis=1) / sigma2
    recon -= 0.5 * x.shape[1] * np.log(2.0 * np.pi * sigma2)
    return recon - gaussian_kl(mu, log_var)

The code computes the bound per example rather than per batch. Summing over a batch and then dividing by its size changes the effective weight of the KL by a factor of the batch size relative to a summed reconstruction term. That discrepancy shows up as a model whose latents are either ignored or over-regularised depending on how the loss was reduced, and it is invisible in the total loss value because both terms scale together.

What the gradient of the bound actually rewards

Differentiating with respect to θ is straightforward: the expectation is over q, which does not depend on θ, so the gradient passes inside and hits log p(x | z). Differentiating with respect to φ is the hard direction, because φ appears in the distribution the expectation is taken over. The score-function estimator, which multiplies the integrand by the gradient of the log density, is unbiased but has variance that grows quickly with the latent dimension, and it is rarely used for continuous latents. The reparameterisation trick, covered in the path from sampling to gradients, rewrites the sample as a deterministic function of φ and an independent noise draw, which removes the distribution from the differentiation path entirely.

With reparameterisation, the gradient of the KL term is available in closed form and the gradient of the reconstruction term is a single-sample estimate. One sample per example per step is standard. The estimator is unbiased for the gradient of the bound evaluated at the current parameters, though the bound itself is estimated, so the training signal has variance that depends on how noisy the decoder is.

There is a second-order effect that matters for scheduling. Early in training the encoder is far from the prior, so the KL term is large and its gradient dominates; late in training the reconstruction term dominates because the KL has been driven down. Any objective that reweights the two terms, such as a warm-up coefficient on the KL, is exploiting this ordering rather than changing the bound.

Where the bound is loose and what that costs

Three failure patterns follow directly from the structure of the bound rather than from optimisation bugs.

The first is that the KL term can win. If the decoder is weak, the cheapest way to reduce the divergence is to make q(z | x) equal the prior for every x, at which point z carries no information and the decoder falls back on the marginal mean of the data. This is posterior collapse, and it is a property of the objective, not a defect in the optimiser.

The second is that a diagonal Gaussian approximation cannot represent a posterior with several well-separated modes. The bound then forces the approximation to cover all of them, which pulls decoded outputs toward an average of the corresponding reconstructions. The model is not failing to fit; it is fitting the average because that is what maximises the bound under the constraint it was given.

The third is that the bound is a bound on log p(x), not on the geometry of samples. A model can have a tight bound and poor samples if the prior and the aggregate posterior disagree in ways the bound tolerates. Tighter estimators exist, notably the importance-weighted bound, which averages K samples of the ratio p(x, z)/q(z | x) inside a logarithm. It is a lower bound on the marginal likelihood and becomes tighter as K grows, but each additional sample costs a full decoder pass, and the gradient estimator’s variance at large K needs care because the logarithm of a sample average is dominated by its largest term.

The bound remains the default because it is the only construction that gives a differentiable, single-pass training signal for a model with an intractable marginal. Every alternative either needs an adversarial objective with its own stability problems, or needs a likelihood that is exactly computable, which restricts the architecture family. The diffusion models described in forward corruption as a fixed linear process avoid the intractable integral by changing what the model is asked to represent, rather than by bounding the integral.