Approximate Inference: Deterministic Methods

At the heart of Bayesian inference lies a fundamental computational challenge: the posterior distribution p(θ|D) = p(D|θ)p(θ)/p(D) is often intractable. The marginal likelihood p(D) = ∫p(D|θ)p(θ)dθ requires integrating over the entire parameter space—an operation that is analytically impossible for most realistic models.

We've seen how sampling methods (MCMC, importance sampling) handle this by drawing samples from the posterior. But these methods can be slow, require careful tuning, and provide stochastic rather than deterministic answers. For many applications—especially those requiring fast predictions or online learning—we need deterministic approximations that provide closed-form, tractable representations of complex posteriors.

The Laplace approximation is the simplest and most elegant of these methods. It approximates any posterior with a Gaussian distribution centered at the posterior mode. Despite its simplicity, it remains surprisingly effective for many machine learning applications and forms the conceptual foundation for more sophisticated methods.

The Laplace approximation is named after Pierre-Simon Laplace, who used it in the 18th century for asymptotic expansions of integrals. The core insight is deceptively simple:

Any smooth, peaked probability density can be locally approximated by a Gaussian near its mode.

This isn't just a heuristic—it emerges from the mathematics of Taylor expansion. When a probability density p(θ) has a single sharp peak at θ̂, the behavior near that peak dominates all integrals involving p(θ). If we can approximate this local behavior with a Gaussian, we obtain a tractable representation that captures the essential structure of the distribution.

Why Gaussians are special:

Gaussians are the "natural" choice for approximating peaked distributions because:

1. Quadratic log-density: log p(θ) is a quadratic function of θ, which is exactly what a second-order Taylor expansion produces

2. Analytical tractability: All moments, marginals, and conditionals of Gaussians have closed forms

3. Maximum entropy: Among distributions with specified mean and covariance, the Gaussian has maximum entropy—it assumes nothing beyond these moments

4. Closure under marginalization: Marginals of Gaussians are Gaussian, enabling recursive approximation schemes

Let's derive the Laplace approximation rigorously. Consider a posterior distribution:

$$p(\theta|D) = \frac{p(D|\theta)p(\theta)}{p(D)} \propto p(D|\theta)p(\theta)$$

Define the negative log posterior (ignoring the normalization constant):

$$\Psi(\theta) = -\log p(D|\theta) - \log p(\theta)$$

The posterior can be written as:

$$p(\theta|D) \propto \exp(-\Psi(\theta))$$

Our goal is to approximate Ψ(θ) with a simpler function that yields a Gaussian when exponentiated.

Step 1: Find the mode

First, find the Maximum A Posteriori (MAP) estimate:

$$\hat{\theta} = \arg\min_\theta \Psi(\theta) = \arg\max_\theta p(\theta|D)$$

At the mode, the gradient vanishes:

$$\nabla\Psi(\hat{\theta}) = 0$$

Step 2: Second-order Taylor expansion

Expand Ψ(θ) around θ̂ using Taylor's theorem:

$$\Psi(\theta) \approx \Psi(\hat{\theta}) + \underbrace{\nabla\Psi(\hat{\theta})^T(\theta - \hat{\theta})}_{= 0 \text{ (at mode)}} + \frac{1}{2}(\theta - \hat{\theta})^T H (\theta - \hat{\theta})$$

where H is the Hessian matrix of Ψ evaluated at θ̂:

$$H = \nabla^2\Psi(\hat{\theta}) = -\nabla^2\log p(D|\hat{\theta}) - \nabla^2\log p(\hat{\theta})$$

Step 3: Exponentiate to get the Gaussian

Substituting back:

$$p(\theta|D) \propto \exp(-\Psi(\theta)) \approx \exp\left(-\Psi(\hat{\theta}) - \frac{1}{2}(\theta - \hat{\theta})^T H (\theta - \hat{\theta})\right)$$

This is proportional to a Gaussian:

$$q(\theta) = \mathcal{N}(\theta | \hat{\theta}, H^{-1})$$

where:

• Mean = θ̂ (the MAP estimate)
• Covariance = H⁻¹ (inverse Hessian of negative log posterior)

The complete Laplace approximation:

$$p(\theta|D) \approx \mathcal{N}(\theta | \hat{\theta}, (\nabla^2\Psi(\hat{\theta}))^{-1})$$

The practical challenge in Laplace approximation is computing the Hessian. For a model with d parameters, the Hessian is a d × d matrix containing all second-order partial derivatives:

$$H_{ij} = \frac{\partial^2 \Psi}{\partial \theta_i \partial \theta_j}$$

The Hessian decomposes into contributions from likelihood and prior:

$$H = -\nabla^2\log p(D|\theta) - \nabla^2\log p(\theta)$$

The Fisher Information Matrix:

For many models, the Hessian of the negative log-likelihood equals the Fisher Information Matrix:

$$\mathcal{I}(\theta) = -\mathbb{E}[\nabla^2\log p(D|\theta)]$$

The Fisher Information measures the curvature of the log-likelihood, which determines how much information the data provides about θ. Higher curvature means more concentrated (less uncertain) posteriors.

Observed vs. Expected Fisher Information:

• Observed Fisher: Use actual Hessian at θ̂ (what we compute)
• Expected Fisher: Use expectation over data (sometimes analytically simpler)

For large datasets, these converge. The observed Fisher is preferred when available because it uses actual data rather than expectations.

The Laplace approximation has an elegant geometric interpretation. Imagine the negative log posterior Ψ(θ) as a "energy landscape" where valleys correspond to high-probability regions.

The approximation process:

1. Find the lowest point (minimum of Ψ, maximum of posterior)
2. Fit a paraboloid tangent to the surface at this point
3. This paraboloid defines a Gaussian when exponentiated

The Hessian captures the curvature of this landscape:

• Large eigenvalues → steep curvature → tightly constrained parameters
• Small eigenvalues → gentle curvature → loosely constrained parameters
• Eigenvalue ratios → shape of the Gaussian (elongated vs. spherical)

The eigenvalue interpretation:

The covariance matrix Σ = H⁻¹ has spectral decomposition:

$$\Sigma = V \Lambda V^T$$

where V contains eigenvectors (principal directions) and Λ contains eigenvalues (variances along each direction).

• Eigenvectors point along the principal axes of uncertainty
• Eigenvalues measure uncertainty magnitude along each axis
• Condition number κ = λ_max/λ_min indicates how "elongated" the Gaussian is

A large condition number (κ >> 1) signals potential numerical issues and suggests parameters are on very different scales.

Beyond approximating the posterior, Laplace provides an approximation to the marginal likelihood (model evidence):

$$p(D) = \int p(D|\theta)p(\theta)d\theta$$

This integral is typically intractable, but Laplace approximation yields:

$$\log p(D) \approx \log p(D|\hat{\theta}) + \log p(\hat{\theta}) + \frac{d}{2}\log(2\pi) - \frac{1}{2}\log|H|$$

where d is the number of parameters.

Decomposition of the marginal likelihood:

$$\log p(D) \approx \underbrace{\log p(D|\hat{\theta})}{\text{fit to data}} + \underbrace{\log p(\hat{\theta})}{\text{prior consistency}} + \underbrace{\frac{d}{2}\log(2\pi) - \frac{1}{2}\log|H|}_{\text{Occam factor}}$$

The Occam factor automatically penalizes model complexity:

• More parameters → larger d → smaller Occam factor
• Sharper posterior → larger |H| → smaller Occam factor

This provides a principled trade-off between model fit and complexity, implementing Bayesian Occam's razor.

Computing log|H|:

The log-determinant can be computed efficiently via Cholesky decomposition:

$$\log|H| = 2\sum_{i=1}^d \log L_{ii}$$

where H = LLᵀ is the Cholesky factorization. This is O(d³) but numerically stable.

Implementing Laplace approximation in practice requires attention to several numerical and statistical details.

Scaling to high dimensions:

For models with many parameters (d >> 1000), storing and inverting the full Hessian becomes prohibitive. Several approximations help:

1. Diagonal approximation: Use only diagonal of Hessian, assuming parameter independence $$\Sigma \approx \text{diag}(1/H_{11}, ..., 1/H_{dd})$$

2. Block-diagonal approximation: Capture correlations within groups (e.g., per-layer in neural networks) but not across groups

3. Low-rank approximation: Represent covariance as Σ ≈ D + VVᵀ where D is diagonal and V is low-rank

4. Kronecker-factored approximation (K-FAC): For neural networks, approximate Fisher as Kronecker product of smaller matrices

Laplace approximation is widely used across machine learning, particularly in settings where speed or closed-form solutions are prioritized.

Laplace for Neural Networks (LLA):

Recent work has revived Laplace approximation for deep learning through Last-Layer Laplace Approximation:

1. Train a neural network to MAP estimate using standard methods
2. Apply Laplace only to the final layer weights
3. Propagate uncertainty through the fixed feature extractor

This provides calibrated uncertainty estimates with minimal computational overhead—the Hessian is only computed for the last layer, making it tractable even for large networks.

Laplace for continual learning:

Sequential Laplace approximation naturally supports continual learning:

• After seeing data D₁, posterior is N(θ̂₁, Σ₁)
• This becomes the prior for D₂
• The prior's precision (Σ₁⁻¹) acts as regularization, preventing catastrophic forgetting

Understanding when Laplace approximation fails is as important as knowing how to apply it. Several fundamental limitations restrict its applicability.

Extensions to address limitations:

1. Mixture of Laplace approximations: Fit Gaussian at each mode, combine as mixture $$q(\theta) = \sum_k \pi_k \mathcal{N}(\theta | \hat{\theta}_k, H_k^{-1})$$

2. Reparameterization: Transform parameters so posterior is more Gaussian (e.g., log-transform positive parameters)

3. Higher-order expansions: Include third/fourth-order terms for skewness and kurtosis correction (Edgeworth expansions)

4. Leave-one-out corrections: WAIC and LOO-CV can detect when Laplace marginal likelihood is unreliable

5. Variational refinement: Use Laplace as initialization for more flexible variational approximations

The Laplace approximation provides a remarkably simple yet powerful approach to approximate inference. By fitting a Gaussian to the posterior via second-order Taylor expansion, we obtain tractable representations that enable fast prediction, uncertainty quantification, and model comparison.

What's next:

Laplace approximation provides a single Gaussian centered at the mode, which may be too restrictive for complex posteriors. In the next page, we'll explore variational inference—a more flexible framework that optimizes over a family of distributions to find the best approximation, trading off fidelity and tractability through a principled objective function.