Variational Inference and the Evidence Lower Bound

Every quantity that makes Bayesian computation hard comes from the same integral. The evidence is one, and posteriors of latent variables are others. Variational inference converts the integral into an optimization by introducing a family of candidate distributions and asking which member of the family is closest to the posterior in a divergence that requires no normalizing constant. The objective it maximizes is a lower bound on the log evidence, so the answer comes with an explicit measure of how much it may be understating the target.

An identity that produces an objective

Write the model as a joint over observed data $x$ and latent variables $z$, and take any distribution $q$ over $z$ with full support relative to the posterior. Two applications of the definition of expectation give

$$\log p(x) = \mathrm{KL}\bigl(q(z) ,|, p(z \mid x)\bigr) + \mathbb{E}_q\bigl[\log p(x, z) - \log q(z)\bigr]$$

The decomposition is exact, and the expectation term is the evidence lower bound, or ELBO. Because a Kullback-Leibler divergence is non-negative, the ELBO never exceeds $\log p(x)$, and the gap between them is exactly the divergence from $q$ to the posterior. Tightness has a single cause: the family must contain the posterior. Maximizing the ELBO over a family that cannot contain it returns the closest member rather than the right answer, and the residual gap is not reported by the optimizer, which sees only the bound.

The practical attraction is that the ELBO contains no posterior normalizer. Its first term is an expectation of the log joint under $q$, and its second term is the entropy of $q$; for a model whose complete conditionals are in a tractable family, the first term reduces to expectations of the model’s sufficient statistics, which are known functions of the variational parameters. No integration over the posterior appears anywhere in the objective, and the gradient of the bound with respect to those parameters can be estimated by drawing samples from $q$ rather than from the posterior.

import math

def kl_normal(mu_q, sigma_q, mu_p, sigma_p):
    # KL(N(mu_q, sigma_q**2) || N(mu_p, sigma_p**2)) in closed form
    return (math.log(sigma_p / sigma_q)
            + (sigma_q ** 2 + (mu_q - mu_p) ** 2) / (2 * sigma_p ** 2)
            - 0.5)

The quantity being bounded is the normalizing constant that makes the partition function of an undirected model hard, as treated in Markov Random Fields and the Intractability of the Partition Function. Bounding it is one of two responses. The other is to estimate it from the states of a chain, which returns an answer that is correct in the limit of an unbounded run — Metropolis-Hastings: Detailed Balance and Proposal Design sets out that construction — while the bound costs a fixed number of iterations per evaluation and carries a bias that the objective itself does not measure.

The divergence direction decides which error is made

The gap is $\mathrm{KL}(q | p)$, not $\mathrm{KL}(p | q)$, and the two directions fail differently. Terms of the first kind are large when $q$ places mass where $p$ has little, and they are small when $q$ is near zero where $p$ is large. Minimizing that divergence therefore has no incentive to cover all of $p$‘s mass: the optimizer narrows $q$ until it sits inside the region where $p$ is substantial, which makes the approximation understate variance and, when the posterior has separated modes, drop all but the ones it happens to be near. Choosing the opposite direction would spread $q$ across the whole posterior, including the low-density valleys between modes, and overstate variance instead. The direction is fixed by the identity above; the ELBO inherits it, so a variational posterior that looks sharply confident is the expected behaviour of the objective rather than a sign that the data are informative.

For the Gaussian case the consequence is quantitative. If the family is Gaussian and the target is approximately Gaussian around a mode, the optimal member scales its covariance to the local curvature of the log density at that mode, which matches the target’s variance only when the curvature is constant. Curvature that varies across the posterior’s mass produces a variance between the two extremes, and the ratio between the fitted and true standard deviation grows with the amount of curvature variation.

Mean-field updates are a coordinate-wise optimum

The most common family factorizes the latents into independent groups, $q(z) = \prod_i q_i(z_i)$, which sets every correlation between groups to zero by assumption. To maximize the bound over one factor with the others held fixed, isolate the terms in the ELBO that depend on $q_j$: they are $\mathbb{E}{q_j}\bigl[\mathbb{E}{q_{-j}}[\log p(x, z)]\bigr]$ minus the entropy of $q_j$. The optimum over all distributions of $z_j$ is reached by

$$q_j^*(z_j) ;\propto; \exp\Bigl(\mathbb{E}{q{-j}}\bigl[\log p(x, z)\bigr]\Bigr)$$

which is a statement about the complete conditional of $z_j$: the optimum has the same form as $p(z_j \mid z_{-j}, x)$ with the conditioning variables replaced by their expectations under the current approximation. When the complete conditionals belong to an exponential family, so does every $q_j$, and the update becomes a rule for the natural parameters in terms of expected sufficient statistics. Two properties follow. Each coordinate update cannot decrease the ELBO, because it solves an exact maximization over one block; and the iteration converges to a local optimum, since the bound is concave in each block separately for conditionally conjugate models but not jointly. Every reported posterior is therefore one of possibly several fixed points, selected by initialization.

The per-iteration cost is set by the expectations, not by the coordinate count. Updating $q_j$ requires $\mathbb{E}{q{-j}}[\log p(x, z)]$, which for a latent-variable model with $n$ observations and $K$ components requires touching each observation against each component, so the cost grows linearly in $n$ and $K$ with a factor of $d^2$ or $d^3$ per component depending on how the covariance is parameterized. That linear dependence on $n$ is the term that a full-batch iteration cannot avoid, and it is entirely absent from the bound’s structure.

Minibatches and amortized variational parameters

The ELBO is a sum of per-observation contributions plus the entropy terms, so a uniform minibatch gives an unbiased estimate of the summed part when observations are exchangeable, and an escalated gradient step on that estimate ascends the same objective in expectation. The practical effect is that each iteration costs work proportional to the minibatch rather than to the data set, while the number of iterations required grows, so the trade is the familiar one between per-iteration cost and variance of the step direction. A second option removes the per-datapoint optimization altogether: make the variational parameters outputs of a shared function of the data point, so fitting a new observation costs one function evaluation instead of an inner optimization loop. Estimating the gradient of the bound through such a function requires either a reparameterization that moves the randomness outside the parameters or a score-function estimator whose variance is typically larger.

What the bound never reports

Four limitations follow from the objective rather than from any implementation. First, the ELBO is a lower bound on one model’s evidence, and comparing bounds across different models compares two gaps; a model with a tighter bound can be ranked above a model with higher evidence. Only comparisons within one family, where the gap is similar, have a supporting argument. Second, no second-order information comes for free: a fitted covariance is a property of the objective’s optimum, not an estimate of posterior curvature, and mean-field covariances are diagonal by construction.

Third, exchangeable structure produces symmetric solutions. In a mixture with $K$ components, the posterior is invariant under permutation of component labels, so the bound is too, and a local optimum selects one arbitrary labelling; the fitted parameters depend on initialization rather than on the data. Fourth, the objective at convergence quantifies the fit only relative to the family. Running several initializations and recording the bound for each is the cheapest available check on whether the optimizer landed in a poor local optimum, and a large spread across runs is a signal about the family rather than about the data.

The natural next question is where the bound’s maximizing distribution comes from when the model has parameters to fit as well as latents to infer. Maximizing the same bound over a factorized family and over model parameters, with the family member replaced by the exact posterior at each parameter value, is precisely the algorithm treated in Expectation-Maximization as Inference with Latent Variables, which makes the bound’s behavior in the presence of intractable posteriors concrete.