Laplace and variational approximations

What do you do when the posterior has no formula?

Part 16 built a posterior by multiplying a prior by a likelihood, and it worked because the two were conjugate: the product landed back inside a family you already knew, and the normalising constant fell out of a table. Conjugacy is a luxury. A logistic regression or a hierarchical model produces a product whose shape you can evaluate but whose integral you cannot do. The posterior exists and you cannot write it down.

One response is to sample, as Part 19 did: build a Markov chain whose stationary distribution is the posterior and let its draws stand in. Sampling is asymptotically exact, but stochastic, slow, and awkward at scale. The other response, and the subject of this part, is to approximate the posterior by a distribution you can write down: commit to a simple family — almost always Gaussians — and choose the member closest to the truth. The Laplace approximation does this locally, by expanding the log-posterior around its mode; variational inference does it globally, by defining a family of candidates and optimising a divergence. Both turn an integral into an optimisation, and both fit a Gaussian to something that is not one.

The shape is free; the constant is not

Everything begins with Bayes' rule in its unnormalised form. For data $D$ and parameters $\theta$,

The numerator is a function you can evaluate anywhere: give me a value of $\theta$ and I can compute the prior and the likelihood in a few operations. The denominator is the evidence, and it is where the trouble lives. Without it the numerator is just a shape; dividing by it is what makes the shape a density, and that division is the step we cannot perform. In one or two dimensions this is a nuisance — grid, sum, normalise — but the curse of dimensionality ends that: a grid of ten points per dimension is $10^{30}$ points at thirty parameters, and the posterior sits in a thin shell a coarse grid will miss.

The probability chapter of the calculus guide built expectations out of exactly this integral, and every summary you want — a mean, a variance, a credible interval — is another intractable integral against the posterior. There is a practical motive too: a Gaussian posterior can be fused with the next measurement by the precision-weighted rule of Part 17 and stored as a matrix, whereas a cloud of samples is much harder to hand to the next stage of a pipeline.

Variational inference, and the asymmetry that shapes it

Variational inference fixes a family of candidates $q$ — here, every Gaussian — and searches for the member closest to the posterior $p$. "Closest" needs a definition, and the standard one is the Kullback–Leibler divergence. It is asymmetric, and that asymmetry is the single most consequential design choice in the method:

The first direction, $\mathrm{KL}(q \| p)$, averages the log-ratio under the approximation. Wherever $q$ places mass and $p$ places almost none, the integrand explodes, so the optimiser is punished for putting $q$ where $p$ is not. The fit is mode-seeking: it locks onto one mode and shrinks until it avoids every region of low posterior density — an effect called zero-forcing — and the price is that it understates uncertainty and can ignore a secondary mode. The second direction, $\mathrm{KL}(p \| q)$, averages under the posterior, so the punishment falls where $p$ has mass and $q$ does not, and $q$ stretches to cover everything the posterior cares about. It is mass-covering, and for a Gaussian $q$ minimising it matches the posterior's mean and covariance exactly.

Variational inference almost always minimises the first direction, for computability. The evidence decomposes as

and since the evidence is fixed, maximising the evidence lower bound is the same as minimising $\mathrm{KL}(q\|p)$. The ELBO needs no evidence integral, only expectations under $q$ of quantities you can evaluate, which for a Gaussian $q$ are often closed-form. $\mathrm{KL}(p\|q)$ has no such free lunch: its expectation is under the very posterior you cannot sample from. The asymmetry is therefore not a matter of elegance — it is the direction the computation permits, and its zero-forcing bias is the price paid for tractability.

Set the skew high and watch the two objectives disagree. Minimising $\mathrm{KL}(q\|p)$ parks a narrow Gaussian on the peak, shading the right tail it has decided to ignore, and underestimates the variance badly. Minimising $\mathrm{KL}(p\|q)$ produces a wider Gaussian that straddles the posterior's bulk, matching its mean and variance at the cost of assigning density where the posterior dislikes it. One asks for the most probable region, the other for the best average representation — which is why a variational autoencoder's latent code can look sharp and slightly overconfident, and why calibration matters when an approximate posterior reaches production (the evaluation chapter).

Mode plus curvature, in one Taylor expansion

Laplace takes the opposite tack: it does not search over shapes, it looks very closely at the peak. Let $\theta^{*}$ be the posterior mode, found by maximising the log-posterior as Part 7 maximised a log-likelihood, and expand the log-posterior in a Taylor series there. The first derivative vanishes, so the first non-trivial term is the quadratic:

Exponentiating gives a Gaussian whose normalising constant is free, so the approximation is explicit:

The mean is the mode and the covariance is the inverse of the curvature. $H$ is the observed information, the negative Hessian at the peak: exactly the object Part 8 called curvature, and exactly the matrix a Newton or Gauss–Newton optimiser assembles on its way to the maximum. This is why the method is so cheap — if you found the mode by Newton's method you already factorised the Hessian, and the optimisers chapter develops the mechanics. In code, Stats.bayes.laplace(logpost, mode) returns the Gaussian by differencing the log-posterior twice.

The demo fits that Gaussian to a skewed posterior. At low skew the posterior is nearly symmetric and the fit is excellent, with mean, variance, and the first two posterior moments almost coinciding. As skew grows the mode drifts away from the mean and the curvature at the mode describes only the sharp inner shoulder, so the fitted Gaussian is centred in the wrong place and far too narrow: the long right tail contributes mass but almost no curvature.

Errors have a comprehensible structure. The mean is the mode, so the fit inherits whatever gap separates mode from mean; the variance is the reciprocal of local curvature, so it ignores tails that carry mass without bending the log-density. The error is governed by the third derivative — the term the quadratic threw away — and shrinks as data accumulate, so for small samples, strong skew, or several modes the approximation is a rough sketch rather than a portrait.

Two caveats follow. The mode must be found, and for a nonconvex log-posterior that is itself a global optimisation problem; a local maximum yields a Gaussian around the wrong point. And the curvature must be positive definite — a negative eigenvalue means the point was not a maximum and $H^{-1}$ is not a covariance.

Local and cheap, or global and flexible

The two approximations trade different currencies. Laplace is local: one optimisation and a Hessian, no iteration over distributions and no samples, giving a Gaussian whose fidelity is a question about the log-posterior near one point. Variational inference is global: it commits to a family and optimises a bound over many gradient steps, and can represent skew, correlation, and even multimodality if the family is rich enough — at the cost of far more work.

The decision rule follows. When the posterior is close to Gaussian — a well-identified model, plenty of data, no boundary pathologies — Laplace is hard to beat. When it is skewed, heavy-tailed, or multimodal, a single quadratic will mislead and a flexible family is the better investment, provided you remember that zero-forcing will still make it too confident. When the parameter space is enormous, a full Hessian is out of reach, but a diagonal or low-rank Laplace fit, or a mean-field posterior, remains tractable.

One diagnostic is worth remembering: when the two disagree about the variance, that disagreement is a signal. A Laplace covariance much smaller than a variational one suggests the log-posterior curves hard at the mode and then flattens — the signature of skew or heavy tails. The numerics chapter supplies the conditioning vocabulary and the scaling chapter explains why low-rank curvature matters. The same ideas reappear under other names: a recursive filter runs a Laplace-style update in a loop, the pose-graph chapter assembles the same information matrix, and a variational autoencoder is variational inference with an amortised family.

The Gaussian belief and the optimisation it rides on

The pattern generalises. Whenever a system reports uncertainty about a fitted quantity — a pose, a calibration, a regression coefficient, a benchmark score — it reports the inverse curvature of some negative log-posterior at an optimum, and is therefore reporting a Laplace approximation whether or not it uses the name. Whenever a generative model learns a latent space by maximising a lower bound rather than an exact likelihood, it is doing variational inference, and zero-forcing shapes how much uncertainty that latent space admits. The serving metrics chapter inherits the vocabulary: such an estimate should be read as exactly that, with the failure modes above attached.

What unites the two is the question this part began with. The normalising constant is unavailable, so we replace the posterior with something we can write: one answer looks locally at the peak, the other searches globally for the closest member of a family. Knowing which you have, and which way its errors bend, is the difference between using an approximation and being used by it.

Further reading

These references treat approximate inference as a first-class subject rather than a numerical afterthought. If you take away one thing, take away the picture of a posterior replaced by a Gaussian, and of the KL divergence deciding where that Gaussian may put its mass.

• David MacKay, Information Theory, Inference, and Learning Algorithms, chapters 27–28 — the Laplace method and variational free energy, derived from the evidence decomposition.
• Christopher Bishop, Pattern Recognition and Machine Learning, chapters 4 and 10 — the Laplace approximation and variational inference, with the ELBO developed carefully.
• David Blei, Alp Kucukelbir and Jon McAuliffe, "Variational Inference: A Review for Statisticians", JASA 2017 — the modern survey, including the consequences of the KL direction.
• Luke Tierney and Joseph Kadane, "Accurate Approximations for Posterior Moments", JASA 1986 — when Laplace is accurate, and how to correct it.

Cheat sheet