📊 statistics

Corrected Integrated Laplace Approximation for Bayesian Inference in Latent Gaussian Models

This paper proposes an importance sampling scheme to correct the errors introduced by the integrated Laplace approximation in Bayesian inference for latent Gaussian models, enabling convergence to the correct posterior through techniques like pseudo-marginalization and randomized quasi-Monte Carlo within an automatic differentiation framework.

Original authors: Jinlin Lai, Charles C. Margossian, Daniel R. Sheldon

Published 2026-05-21

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Jinlin Lai, Charles C. Margossian, Daniel R. Sheldon

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to solve a massive, multi-layered puzzle. In the world of statistics, this puzzle is called a Latent Gaussian Model (LGM). It's a way of understanding data that has hidden patterns (the "latent" parts) and observable results.

The problem is, the hidden parts of the puzzle are so numerous and complex that trying to solve the whole thing at once is like trying to drink from a firehose. It's computationally impossible for standard tools.

The Old Shortcut: The "Best Guess" Map

To make this manageable, statisticians use a clever shortcut called the Integrated Laplace Approximation (ILA).

Think of the hidden parts of your puzzle as a foggy mountain range. To navigate it, you don't map every single tree and rock. Instead, you find the highest peak (the "mode") and assume the terrain around it looks like a smooth, perfect hill. You draw a map based only on that peak.

This is fast and efficient. However, the paper points out a flaw: The map is wrong. The real mountain isn't a perfect smooth hill; it has bumps, valleys, and weird shapes. Because the map is an oversimplification, the path you take based on it leads you to the wrong destination. In statistics, this means your final conclusions (the "posterior") are biased and inaccurate.

The New Solution: The "Correction Crew"

The authors of this paper propose a way to keep the speed of the shortcut but fix the errors. They use a technique called Importance Sampling, which they describe as sending out a "correction crew" to check the map.

Here is how their three new methods work, using analogies:

1. The "Pseudo-Marginal" Method (PM-ADLA)

The Analogy: Imagine you are drawing a map, but instead of just looking at the peak, you send out a team of explorers to randomly check different spots on the mountain. You then average their reports to get a better picture.
How it works: The math gets a bit heavy here. They add "noise" (random variables) to the system. By averaging many random samples, the error cancels out, and the map becomes mathematically perfect in the long run.
The Trade-off: It's very accurate, but because you have to send out a whole team of explorers for every step, it slows down the process significantly. It's like trading speed for absolute precision.

2. The "Quasi-Monte Carlo" Method (QMC-ADLA)

The Analogy: Instead of sending explorers randomly (who might all clump together in one area), you send them out in a perfectly organized grid, like a chessboard, to ensure they cover every inch of the mountain evenly.
How it works: They use a special, non-random sequence of points (called a low-discrepancy sequence) to sample the mountain. This is much more efficient than random sampling.
The Trade-off: It's faster than the first method. However, the paper found that sometimes, even with a perfect grid, the map can still have a "blind spot" if the mountain has a very strange shape. The error gets smaller as you add more grid points, but it might not disappear completely for difficult problems.

3. The "Randomized Quasi-Monte Carlo" Method (RQMC-ADLA)

The Analogy: This is the "Goldilocks" solution. You take your organized chessboard grid, but you give the whole board a tiny, random shake before you start.
How it works: They combine the organization of the grid with a tiny bit of randomness. This keeps the efficiency of the grid but removes the "blind spots" that happened when the grid was too rigid.
The Trade-off: This method is the star of the show. It is fast, it doesn't make the puzzle bigger (unlike the first method), and it fixes the errors better than the others. However, because of the "shake" (the modulo operation), the terrain becomes slightly "jagged," so the explorers have to take smaller, more careful steps to navigate it.

The Results: Why It Matters

The authors tested these methods on three types of puzzles:

Synthetic Gaussian Processes: A made-up dataset designed to be tricky.
Sparse Kernel Interaction Models: A model used for real-world data with complex interactions.
Mixed-Effects Models: Used for clinical trial data (like counting seizures in patients).

What they found:

The Old Way (Standard Laplace): Fast, but the answers were consistently wrong (biased).
The "Base" Way (No Shortcut): The most accurate, but so slow and unstable that it often crashed or got stuck (divergent transitions).
The New Ways: They found that their new methods (especially the Randomized one) gave answers that were much closer to the truth than the old shortcut, without crashing like the "Base" way. They corrected the "wrong map" problem while keeping the speed benefits.

In a Nutshell

The paper says: "We found a way to fix the broken maps statisticians have been using for years. We didn't throw away the shortcut; we just added a smart correction system that makes the shortcut accurate again. This lets us solve complex statistical puzzles faster and more correctly than before."

They implemented this in modern software tools (using JAX and BlackJAX) so that other researchers can use these "corrected maps" immediately.

Technical Summary: Corrected Integrated Laplace Approximation for Bayesian Inference in Latent Gaussian Models

Problem Statement

Latent Gaussian Models (LGMs) are a ubiquitous class of Bayesian hierarchical models, including Gaussian processes, spatial models, and mixed-effect models. They are defined by a joint density $\pi(\theta, z, y) = \pi(\theta)\pi(z|\theta)\pi(y|\theta, z)$ , where $\theta$ represents hyperparameters, $z$ represents latent Gaussian variables, and $y$ represents observations.

A primary challenge in Bayesian inference for LGMs is the high dimensionality of $z$ and the complex geometry of the joint posterior $\pi(\theta, z|y)$ , which often hinders efficient sampling. A common strategy is to marginalize out $z$ to infer $\pi(\theta|y)$ , a lower-dimensional target. When the likelihood $\pi(y|\theta, z)$ is non-Gaussian, exact marginalization is intractable. The standard approach, Integrated Nested Laplace Approximation (INLA), uses a Laplace approximation to approximate the conditional posterior $\pi(z|\theta, y)$ as a Gaussian centered at the mode $\hat{z}$ .

While efficient, this approximation introduces a systematic "posterior error." When the likelihood is non-Gaussian, the approximate marginal posterior $\hat{\pi}(\theta|y)$ differs from the true posterior $\pi(\theta|y)$ , even asymptotically. Existing methods like the Adjoint-Differentiated Laplace Approximation (ADLA) allow for gradient-based sampling (e.g., Hamiltonian Monte Carlo) on this approximate posterior but do not correct the underlying bias introduced by the Laplace approximation. Running MCMC on the full joint posterior $\pi(\theta, z|y)$ is asymptotically exact but often computationally infeasible due to the high dimensionality of $z$ .

Methodology

The authors propose an importance sampling framework to correct the error introduced by the Laplace approximation while retaining the computational benefits of marginalization. The core insight is that the integrated Laplace approximation can be reformulated as an importance sampling estimator with a single sample. By increasing the number of samples, the estimator converges to the true marginal.

The paper details three specific implementations of this idea, all integrated into an automatic differentiation framework (JAX) to support gradient-based inference:

Pseudo-Marginal ADLA (PM-ADLA):
- Mechanism: This approach targets an extended distribution over the parameters $\theta$ and auxiliary random variables $\epsilon$ used to generate samples $z_i$ from the Laplace approximation. Specifically, $z_i = \hat{z} + \sqrt{\Sigma}\epsilon_i$ , where $\epsilon_i \sim \mathcal{N}(0, I)$ .
- Properties: It provides an unbiased estimator of the marginal likelihood. As the number of samples $n \to \infty$ , the posterior converges to the correct one.
- Trade-off: It increases the dimensionality of the inference problem by $(n-1)d_z$ (where $d_z$ is the dimension of $z$ ) and requires $n$ density evaluations per HMC step.
Quasi-Monte Carlo ADLA (QMC-ADLA):
- Mechanism: Instead of random auxiliary variables, this method fixes the randomness using a low-discrepancy sequence (LDS), such as the Sobol sequence. The auxiliary variables $u_i$ are fixed points in the unit cube, and $z_i$ are generated via inverse-CDF transformation.
- Properties: This reduces the estimation variance of the importance sampling to $O(n^{-2})$ under proper conditions, offering faster convergence than standard Monte Carlo. It maintains a lower-dimensional parameter space (same as ADLA).
- Trade-off: While asymptotically correct, the paper notes that for finite $n$ , the bias can remain significant, and the fixed nature of the points can lead to large errors in practice even with moderate $n$ .
Randomized Quasi-Monte Carlo ADLA (RQMC-ADLA):
- Mechanism: This unifies the previous two approaches. It uses a single random shift variable $U$ applied to a fixed LDS to generate the auxiliary variables $u_i = (v_i + U) \pmod 1$ .
- Properties: It provides an unbiased estimator of the marginal likelihood and does not increase the model dimension (the auxiliary variable $U$ has the same dimension as $z$ , but is shared across samples).
- Trade-off: The modulo operation introduces discontinuities in the target density with respect to $U$ , making standard HMC inapplicable. The authors propose a Metropolis-within-Gibbs (MwG) sampler: HMC updates $\theta$ conditional on $U$ , and a Metropolis-Hastings step updates $U$ conditional on $\theta$ .

Recovering Latent Variables:
The paper also describes a method to recover samples of the latent variables $z$ that are closer to the true posterior $\pi(z|\theta, y)$ than the standard Laplace approximation. By computing importance weights for the $n$ samples $z_1, \dots, z_n$ , one can construct a discrete approximation of the true conditional posterior, effectively correcting the latent variable inference as well.

Key Contributions

Theoretical Reformulation: The authors demonstrate that the integrated Laplace approximation is a specific case of an importance sampling estimator with one sample. They prove that increasing the sample size in this framework yields an asymptotically correct posterior.
Algorithmic Variants: They propose and implement three distinct algorithms (PM-ADLA, QMC-ADLA, RQMC-ADLA) that balance computational efficiency, dimensionality, and asymptotic correctness.
Gradient Computation: The work extends the ADLA framework to compute gradients for these importance sampling estimators, enabling the use of Hamiltonian Monte Carlo (HMC) and other gradient-based samplers.
Latent Variable Recovery: A procedure is provided to obtain improved posterior samples for the latent variables $z$ using the importance weights derived during the marginal inference.

Experimental Results

The methods were evaluated on synthetic Gaussian processes with Poisson likelihoods, sparse kernel interaction models (SKIMs) on real data, and mixed-effects models (Epil2 dataset).

Error Reduction: In all experiments, the proposed methods (PM, QMC, RQMC) significantly reduced the estimation error of parameter means compared to standard ADLA.
Convergence:
- PM-ADLA: Corrected the posterior error effectively but was computationally slower due to the increased dimensionality and multiple density evaluations.
- QMC-ADLA: Showed reduced variance but exhibited persistent bias in some cases (e.g., estimating means of specific parameters in SKIMs) even with large $n$ .
- RQMC-ADLA: Achieved a balance, correcting the bias of ADLA with performance comparable to standard HMC in terms of error, though the Gibbs sampling step incurred higher computational costs per iteration.
Divergences: In models where standard HMC on the full posterior suffered from divergent transitions (indicating poor geometry), all marginalization methods (including ADLA and the proposed corrections) eliminated divergences, demonstrating that the approximate marginal geometry is well-behaved.
Efficiency: While ADLA was the fastest, it suffered from bias. The proposed methods offered a tunable trade-off: increasing $n$ reduced bias at the cost of computation time.

Significance and Claims

The paper claims to provide a "corrected" version of the integrated Laplace approximation that retains the scalability of marginalization while offering asymptotic guarantees of correctness. The significance lies in bridging the gap between the efficiency of approximate methods (like INLA/ADLA) and the accuracy of exact MCMC.

The authors emphasize that their methods allow users to leverage the flexibility of probabilistic programming languages (via automatic differentiation) without being constrained by the accuracy limitations of the Laplace approximation. They position RQMC-ADLA as a particularly promising approach that avoids increasing model dimensionality while correcting the posterior error, albeit with a more complex sampling scheme (MwG) to handle discontinuities. The work suggests that for many applications, the computational cost of these corrections is justified by the reduction in posterior bias, enabling more reliable downstream applications.