Laplace Approximation

Laplace Approximation

The Laplace approximation approximates a posterior $p(\theta \mid y)$ by the Gaussian obtained from a second-order Taylor expansion of the log-posterior around its mode. Let $\theta \in {\mathbb{R}}^d$, let $\ell (\theta )=\log [p(\theta )p(y\mid \theta )]$ be the log of the unnormalized posterior, and let $\hat{\theta }=\text{MAP}=\arg {\max }_{\theta }\ell (\theta )$. A Taylor expansion at the mode, using $\nabla \ell (\hat{\theta })=0$, gives
$\ell (\theta )\approx \ell (\hat{\theta })+\frac{1}{2}(\theta -\hat{\theta }{)}^{\top }H(\hat{\theta })(\theta -\hat{\theta }),$
where $H(\hat{\theta })={\nabla }^2\ell (\hat{\theta })$ is the Hessian at the mode (negative definite, since $\hat{\theta }$ is a maximum). Exponentiating and recognizing the resulting density as a Gaussian kernel gives
$p(\theta \mid y)\approx N\left(\hat{\theta },{\Sigma }_L\right),{\Sigma }_L=[-H(\hat{\theta }){]}^{-1}.$
Finding the mode is typically done with Newton’s method, which uses the same Hessian and converges quadratically near $\hat{\theta }$.

Integrating the Gaussian kernel gives an approximation to the marginal likelihood
$p(y)=\int p(\theta )p(y\mid \theta )d\theta \approx (2\pi {)}^{d/2}|-H(\hat{\theta }){|}^{-1/2}p(\hat{\theta })p(y\mid \hat{\theta }),$
which is the standard shortcut for evidence computation and Bayes factors when no conjugate form is available.

Three main uses:

1. Direct posterior approximation. Report $N(\hat{\theta },{\Sigma }_L)$ as the posterior. Credible intervals come from Gaussian quantiles: ${\hat{\theta }}_j\pm 1.96\sqrt{{\Sigma }_L[j,j]}$ for a $95\%$ marginal interval. Fastest available Bayesian analysis; accurate for smooth, unimodal, high-$n$ posteriors.
2. MCMC proposal covariance. Use ${\Sigma }_L$ as the covariance of a random-walk Metropolis proposal, typically scaled by $c^2\cdot {\Sigma }_L$ with $c\approx 2.38/\sqrt{d}$ (Roberts-Gelman-Gilks) to target the optimal acceptance rate.
3. Marginal likelihood approximation. Plug into the Laplace evidence formula above for model comparison via Bayes factors.

Validity. The approximation is exact when the posterior is Gaussian (e.g., Normal-Normal conjugacy with known variance). Its error is $O(n^{-1})$ for smooth parametric models by the Bernstein-von Mises theorem; in the large-data limit the true posterior converges to exactly this Gaussian. It fails for multi-modal posteriors (the approximation captures only one mode), bounded parameters near the boundary (the Gaussian extends past the support), and heavy-tailed posteriors (the Gaussian under-weights the tails).

Intuition

The log-posterior of a regular model looks roughly quadratic near its mode: zero gradient, negative-definite Hessian, higher-order terms small for most of the posterior’s mass. A Gaussian has exactly this log-form everywhere, so approximating the log-posterior by its second-order Taylor expansion amounts to replacing the posterior by the Gaussian that best matches its curvature at the peak. For large $n$, the likelihood sharpens faster than the prior shifts the mode, so the approximation gets better — the Bernstein-von Mises theorem makes this precise.

The practical appeal is pragmatic: Newton’s method to find $\hat{\theta }$ and a Hessian evaluation at $\hat{\theta }$ together give you a Bayesian posterior approximation in the same compute budget as computing the MAP. When the posterior is smooth and unimodal, the Laplace approximation is essentially free and good enough for many purposes. When it isn’t enough, its MAP and covariance become the natural starting point for MCMC — the Hessian-scaled proposal is almost always the right thing to use for a Metropolis random walk.

Properties

(Newton’s Method Finds the MAP) The same Hessian that defines the approximation drives the optimization.

Newton’s step for maximizing $\ell$ is ${\theta }^{(k+1)}={\theta }^{(k)}-H({\theta }^{(k)}{)}^{-1}\nabla \ell ({\theta }^{(k)})$, using the Hessian of the log-posterior. Near the mode the method converges quadratically: the error $\Vert {\theta }^{(k)}-\hat{\theta }\Vert$ squares every iteration, so a dozen iterations from a reasonable start give machine precision. Once $\hat{\theta }$ is found, the same $H(\hat{\theta })$ defines the covariance ${\Sigma }_L=-H^{-1}$. In practice the cost of the Laplace approximation is dominated by the Hessian factorization, $O(d^3)$, which is fine for $d$ up to a few thousand.

(Fisher Information Connection) The Hessian of the log-likelihood equals the observed Fisher information.

Writing $\ell (\theta )=\log p(\theta )+\log p(y\mid \theta )$, the Hessian decomposes as $H=H_{\text{prior}}+H_{\text{lik}}$. The likelihood part $-H_{\text{lik}}(\hat{\theta })$ is the observed Fisher information $\hat{\mathcal{I}}(\hat{\theta })$. For a flat prior and large $n$, $-H(\hat{\theta })\approx \hat{\mathcal{I}}(\hat{\theta })\approx n\mathcal{I}({\theta }_0)$ where $\mathcal{I}$ is the Fisher information per observation and ${\theta }_0$ is the true parameter. The Laplace covariance ${\Sigma }_L=[-H{]}^{-1}$ then matches the frequentist asymptotic covariance $[n\mathcal{I}({\theta }_0){]}^{-1}$ of the MLE, which is the concrete form of Bernstein-von Mises.

(Asymptotic Exactness — Bernstein-von Mises) Under regularity, the Laplace approximation converges to the true posterior as $n\to \infty$.

For smooth parametric models with identifiability, interior true parameter, and positive prior density at ${\theta }_0$, the total variation distance between the true posterior and the Laplace Gaussian is $O(n^{-1/2})$, and the KL divergence is $O(n^{-1})$. Asymptotically, every regular Bayesian posterior is Gaussian centered at the MLE with covariance $[n\mathcal{I}({\theta }_0){]}^{-1}$, independent of the prior. This justifies Laplace as a default fast approximation when $n$ is large, and it implies that reporting MAP $\pm$ Laplace standard errors is asymptotically equivalent to reporting MLE $\pm$ standard Wald errors.

(Marginal Likelihood Estimate) The Laplace evidence formula is $O(n^{-1})$-accurate for smooth unimodal posteriors.

Writing the evidence as $p(y)=\int e^{\ell (\theta )}d\theta$ and substituting the Taylor expansion gives
$p(y)\approx e^{\ell (\hat{\theta })}(2\pi {)}^{d/2}|-H(\hat{\theta }){|}^{-1/2}.$
Taking logs, $\log p(y)\approx \ell (\hat{\theta })+\frac{d}{2}\log (2\pi )-\frac{1}{2}\log |-H(\hat{\theta })|$. The approximation is $O(n^{-1})$ for smooth likelihoods; comparing models by plugging this into log-Bayes-factors $\log p(y\mid M_1)-\log p(y\mid M_2)$ gives a practical model-comparison tool. Two cautions: the Laplace evidence is sensitive to prior tails (since the normalizing constant depends on the prior’s shape away from $\hat{\theta }$), and the approximation degrades when the dimension $d$ grows with $n$.

(Failure Modes — Where the Approximation Breaks) Laplace misses multi-modality, boundary effects, and heavy tails.

Three canonical failures: (a) multi-modal posteriors, where the Gaussian at one mode ignores the others entirely; a mixture-of-Laplace approximations or MCMC is needed. (b) boundary effects, where a parameter lives on $[0,\infty )$ or $[0,1]$ and the mode lies near the boundary; the Gaussian extends past the support and assigns mass to impossible regions. Reparameterizing (e.g., $\log \theta$ for rate parameters, $\log [\theta /(1-\theta )]$ for probabilities) usually fixes this by moving the boundary to $\pm \infty$. (c) heavy tails, where the true posterior has tails fatter than any Gaussian; the Laplace approximation under-weights the tails and produces credible intervals that are too short. For these cases, a $t$-distribution approximation or direct MCMC is preferred.

(Proposal Covariance for MH) Scaled Laplace covariance is the default random-walk Metropolis proposal.

The optimal scaling for a $d$-dimensional random-walk Metropolis on a target with posterior covariance $\Sigma$ is $c^2\Sigma$ with $c=2.38/\sqrt{d}$ (Roberts-Gelman-Gilks, 1997), targeting an acceptance rate near $0.234$. Using the Laplace covariance ${\Sigma }_L$ as the plug-in estimate of $\Sigma$ is practical and nearly optimal whenever the posterior is approximately Gaussian. When the posterior is strongly non-Gaussian, ${\Sigma }_L$ is still a far better proposal shape than the isotropic ${\sigma }^{2}I$ default — at minimum, it rescales the coordinates by their posterior widths and aligns the proposal with the posterior’s principal axes.

Examples

Bayesian logistic regression via Laplace

With $y_i\in \{0,1\}$ and $P(y_i=1\mid x_i,\beta )=\sigma (x_i^{\top }\beta )$ (where $\sigma$ is the logistic function) and a flat prior on $\beta \in {\mathbb{R}}^p$, the posterior has no closed form. The log-posterior is
$\ell (\beta )={\sum }_{i=1}^n[y_i{x}_i^{\top }\beta -\log (1+e^{x_i^{\top }\beta })],$
with gradient $\nabla \ell (\beta )=X^{\top }(y-p(\beta ))$ (where $p_i(\beta )=\sigma (x_i^{\top }\beta )$) and Hessian $H(\beta )=-X^{\top }W(\beta )X$ (where $W(\beta )=\text{diag}(p_i(1-p_i))$). Both are negative-definite, so the log-posterior is concave and has a unique MAP.

Newton’s iteration: ${\beta }^{(k+1)}={\beta }^{(k)}+(X^{\top }{W}^{(k)}X{)}^{-1}{X}^{\top }(y-p^{(k)})$ — this is exactly the IRLS update for the MLE. Starting from ${\beta }^{(0)}=0$, typically $5$–$15$ iterations reach machine precision.

Laplace approximation: $\beta \mid y\sim N(\hat{\beta },(X^{\top }\hat{W}X{)}^{-1})$ where $\hat{W}=W(\hat{\beta })$. A $95\%$ credible interval for coefficient ${\beta }_j$ is ${\hat{\beta }}_j\pm 1.96\sqrt{[(X^{\top }\hat{W}X{)}^{-1}{]}_{jj}}$ — identical in form to the frequentist Wald interval, which is no coincidence (Bernstein-von Mises again).

Posterior predictive for a new $x_{\ast }$: exact computation requires marginalizing a logistic function against a Gaussian, which has no closed form. The standard approximation uses either Monte Carlo draws from $N(\hat{\beta },{\Sigma }_L)$ or the MacKay probit-approximation $\sigma (x_{\ast }^{\top }\hat{\beta }/\sqrt{1+\pi x_{\ast }^{\top }{\Sigma }_L{x}_{\ast }/8})$.

Laplace as MH proposal for Poisson regression

With $y_i\sim \text{Poisson}(e^{x_i^{\top }\beta })$ and a flat prior, the log-posterior is $\ell (\beta )={\sum }_i[y_i{x}_i^{\top }\beta -e^{x_i^{\top }\beta }]$ with Hessian $H(\beta )=-{\sum }_i{e}^{x_i^{\top }\beta }{x}_i{x}_i^{\top }=-X^{\top }D(\beta )X$, where $D(\beta )=\text{diag}(e^{x_i^{\top }\beta })$. Newton’s method gives the MAP $\hat{\beta }$; the Laplace covariance is ${\Sigma }_L=(X^{\top }\hat{D}X{)}^{-1}$.

A random-walk Metropolis chain with proposal ${\beta }^{\prime }\sim N(\beta ,c^2{\Sigma }_L)$ mixes dramatically faster than the same chain with ${\sigma }^{2}I$. For $p=10$ covariates, isotropic proposals might take $10^5$ iterations to produce $100$ effective samples; the Laplace-shaped proposal reaches the same effective sample size in roughly $10^3$ iterations. This is the concrete payoff of using a Hessian-informed proposal: mixing speed scales with how well the proposal shape matches the posterior shape, and ${\Sigma }_L$ is the best second-order match available without further work.

Evidence for model comparison

Compare two Poisson regression models with $p_1=3$ and $p_2=5$ covariates using Bayes factors. Under each model, compute the MAP ${\hat{\beta }}^{(k)}$ and Hessian $H^{(k)}({\hat{\beta }}^{(k)})$. The Laplace log-evidence is
$\log p(y\mid M_k)\approx {\ell }^{(k)}({\hat{\beta }}^{(k)})+\frac{p_k}{2}\log (2\pi )-\frac{1}{2}\log |-H^{(k)}({\hat{\beta }}^{(k)})|.$
The log-Bayes factor is the difference. The first term rewards fit (higher maximized log-likelihood); the combination $\frac{p_k}{2}\log (2\pi )-\frac{1}{2}\log |-H|$ penalizes model complexity (determinant of the Hessian grows with $p$), giving a built-in Occam penalty. This is the Laplace-BIC connection: expanding $\log |-H({\hat{\beta }}^{(k)})|\approx p_k\log n+O(1)$ for large $n$ recovers the Bayesian Information Criterion $\text{BIC}=-2\ell (\hat{\beta })+p\log n$. BIC is therefore a crude Laplace approximation in which the Hessian determinant is approximated by its leading $n$-dependence.

Derivation

Goal. Derive the Gaussian approximation $p(\theta \mid y)\approx N(\hat{\theta },-H(\hat{\theta }{)}^{-1})$ from a Taylor expansion of the log-posterior.

Let $\ell (\theta )=\log p(\theta )+\log p(y\mid \theta )$ be the log of the unnormalized posterior. Assume $\ell$ is twice differentiable and has a unique global maximum at $\hat{\theta }$ with negative-definite Hessian. A second-order Taylor expansion of $\ell$ at $\hat{\theta }$ gives
$\ell (\theta )=\ell (\hat{\theta })+\nabla \ell (\hat{\theta }{)}^{\top }(\theta -\hat{\theta })+\frac{1}{2}(\theta -\hat{\theta }{)}^{\top }H(\hat{\theta })(\theta -\hat{\theta })+O(\Vert \theta -\hat{\theta }{\Vert }^3).$
The gradient vanishes at the mode, so the linear term drops out. Exponentiating and dropping higher-order terms,
$p(\theta )p(y\mid \theta )\approx p(\hat{\theta })p(y\mid \hat{\theta })\exp (-\frac{1}{2}(\theta -\hat{\theta }{)}^{\top }[-H(\hat{\theta })](\theta -\hat{\theta })).$
The right-hand side, as a function of $\theta$, is the kernel of a Gaussian density with mean $\hat{\theta }$ and precision matrix $-H(\hat{\theta })$. Normalizing gives
$p(\theta \mid y)\approx N\left(\theta ;\hat{\theta },[-H(\hat{\theta }){]}^{-1}\right).$
The normalizing constant absorbs the factor $p(\hat{\theta })p(y\mid \hat{\theta })(2\pi {)}^{d/2}|-H(\hat{\theta }){|}^{-1/2}$, which — equated with the definition of the marginal likelihood — yields the Laplace evidence formula
$p(y)\approx p(\hat{\theta })p(y\mid \hat{\theta })(2\pi {)}^{d/2}|-H(\hat{\theta }){|}^{-1/2}.$

The error of the approximation is controlled by the $O(\Vert \theta -\hat{\theta }{\Vert }^3)$ term. For smooth parametric likelihoods with $n$ observations, the posterior concentrates on a $O(n^{-1/2})$ neighborhood of $\hat{\theta }$, so the cubic remainder is $O(n^{-3/2})$ relative to the quadratic term of order $O(n^{-1})$. The relative error is therefore $O(n^{-1/2})$ in total-variation distance and $O(n^{-1})$ in KL divergence, giving the Bernstein-von Mises rate.

See Also

• 3.1 Bayesian Framework — prior, likelihood, posterior, MAP; Laplace is a Gaussian approximation of this machinery.
• 1.2.2 Maximum a Posteriori Estimation — the MAP that Newton’s method targets.
• 3.3 Bayes Factors — Laplace evidence formula is the standard shortcut for Bayes factors.
• 3.7 Marginal Likelihood — exact and approximate evidence computations, including Laplace.
• 3.11.1 Metropolis-Hastings — uses Laplace covariance as the standard proposal shape.
• 3.8.3 Bayesian Logistic Regression — canonical application of Laplace to a non-conjugate likelihood.
• 1.2.1 Maximum Likelihood Estimation — MLE = MAP under a flat prior; Fisher information = $-$Hessian of log-likelihood.