Enhanced Laplace Approximations

• Enhanced Laplace approximations are advanced techniques that extend the classical Gaussian approximation by using iterative mixtures, skew adjustments, and higher-order corrections to accurately capture multimodality and asymmetry.
• They employ algorithmic refinements such as iterated Laplace and blockwise normalization to improve moment recovery and computational efficiency in settings where standard Laplace fails.
• These methods are widely applied in Bayesian inference, neural networks, and hierarchical models to robustly estimate integrals for complex posteriors with heavy tails and non-ellipticity.

Improved and enhanced Laplace approximations refer to a spectrum of algorithmic and theoretical developments that extend the classical Laplace method, addressing its well-known limitations—such as its inability to represent skewness, multimodality, or heavy tails in integrands/posteriors, and its lack of third-order accuracy in high dimensions or for complex, non-elliptical densities. These advances encompass mixture refinements, higher-order and skew corrections, modifications for Bayesian model selection, and adaptations for structured and hierarchical models, among others. Enhanced Laplace methods have become crucial for modern Bayesian and likelihood-based computation, finding widespread use in statistics, machine learning, and applied sciences.

1. Iterated and Mixture-Based Laplace Approximations

The classical Laplace approximation constructs a single local Gaussian around the mode of the target density $\pi(\theta)$, yielding accurate approximations when $\pi$ is unimodal and nearly elliptically contoured. To overcome its failure in multimodal, asymmetric, or heavy-tailed settings, Bornkamp (2011) introduced the iterated Laplace approximation (iterLap), which iteratively improves a Gaussian-mixture approximation by successively fitting new Gaussian (“Laplace”) components to the residual, i.e., the difference between the current mixture and the true density (Bornkamp, 2011). At each step, a local mode of the residual is found, a local quadratic approximation is formed, and the resulting Gaussian is added to the mixture. After $M$ steps, the approximation $$\hat{\pi}(\theta) = \sum_{m=1}^M w_m \phi(\theta; \mu_m, \Sigma_m)$$ captures non-Gaussian structure far beyond any single Gaussian.

Bornkamp’s iterLap further employs adaptive least-squares fitting of weights over a sampled grid to ensure global $L_2$ optimality. The algorithm demonstrates substantial gains in normalized effective sample size and moment recovery for highly skewed, multimodal, or nonlinear posteriors, with modest computational cost even in moderate dimension ($p \lesssim 20$).

Subsequent work further refines iterLap. Key modifications include omitting premature stopping based on normalization-constant stabilization, introducing absolute residual-based optimization, symmetric residuals for local over- and under-fitting, and Hessian scaling to control component sharpness. See (Mai et al., 2015) for a systematic investigation of these practical enhancements and their impact on accuracy-runtime trade-offs.

2. Higher-Order, Skew, and Asymmetry Corrections

Classical Laplace is second-order: it matches curvature but not skewness or kurtosis. Recent progress establishes analytic and algorithmic frameworks for systematically incorporating third (skewness) and fourth (kurtosis) order information:

• Skew-adjusted Laplace (Katsevich, 2023): In the high-dimensional regime, the leading error of Laplace is $O(d/\sqrt{n})$. A cubic correction using the third derivative tensor provides a signed-density approximation

$$\text{dLA}_S(\theta) = \mathrm{LA}(\theta) \left[1 - \frac{1}{6} \nabla^3 V(\hat{\theta}) : (\theta - \hat{\theta})^{\otimes 3}\right]$$

which provably reduces error to $O((d/\sqrt{n})^2)$, an order-of-magnitude improvement (Katsevich, 2023). Explicit mean and covariance corrections, as well as closed-form formulas for shifted means, are available and computationally trivial to apply in practice.

• Skew-Laplace for nonparametric mixtures (Pozza, Durante & Szabó, 2026): For multimodal or asymmetric posteriors, “skew-Laplace” constructs a signed mixture

$$f_{\mathrm{SL}}(\theta) = 2 f_L(\theta) w(\theta), \quad w(\theta) = \frac{p(\theta | y)}{p(\theta | y) + p(2\hat{\theta} - \theta | y)}$$

which captures higher-order asymmetry and delivers 30%+ reduction in total variation error relative to Gaussian Laplace for Dirichlet process mixtures and related models (Franzolini et al., 28 Apr 2026).

• Refined univariate approximations: The extended simplified Laplace approach (ESLA, (Chiuchiolo et al., 2022)) fits an extended skew-normal (ESN) density to the marginal, matching not only curvature and asymmetry but also excess kurtosis via fourth-derivative expansions. This achieves up to 75% reduction in mode error in cases of pronounced skew, without additional computational cost.
• Blockwise asymmetry corrections in hierarchical models (LPS; (Lambert et al., 2022)): In Bayesian P-spline models, splitting the parameter vector into penalized (well-informed, approximately Gaussian by Bernstein–von Mises) and unpenalized (potentially poorly identified, skewed) subgroups, and applying skew-normal corrections only to the latter, yields near-MCMC accuracy for intervals and means with significantly lower computation.

3. Improved and Enhanced Marginal Likelihood Estimation

Laplace approximations are foundational in integrated likelihoods, but may be inaccurate in high-dimension or complex settings:

• Blockwise and normalized improved Laplace (Ruli, Sartori & Ventura, (Ruli et al., 2015)): By factorizing the target density and enforcing exact normalization on each block, this third-order “improved Laplace” reduces the relative error of the Laplace approximation for integrals from $\pi$0 to $\pi$1. Each marginal or conditional is re-approximated and renormalized, yielding much higher accuracy, even in dimensions up to $\pi$2 for $\pi$3-skewed and spatial random effects posteriors. This approach is implemented in the R package iLaplace.
• Enhanced Laplace Approximation (ELA) (Han et al., 2022): In latent variable models (GLMMs), the ELA substitutes a Monte Carlo moment-corrected average (over the local Gaussian at the mode) for the single-point plug-in, delivering asymptotically unbiased MLE and consistent variance estimators even when classical Laplace yields bias and underestimates uncertainty (e.g., binary/correlated models). ELA seamlessly generalizes to REMLE and is provably consistent as the Monte Carlo sample size increases.
• Approximate Laplace Approximation (ALA) (Rossell et al., 2020): For massive model selection (e.g., high-dimensional variable selection), ALA removes the need for repeated optimization per model, instead performing a single Taylor expansion at a cheap fixed point, and then solves a least squares problem per model. Theoretical results establish that ALA achieves model selection consistency at the same rates as LA, with two to three orders-of-magnitude reduction in computation. Implemented in the mombf R package.

4. Structured, Online, and Subnetwork Laplace for Neural Networks

Scaling Laplace to high-dimensional neural posteriors demands additional structure:

• Online structured Laplace for continual learning (Ritter et al., 2018): Maintains a recursive Gaussian posterior update after each task under a Bayesian online learning framework, with curvature matrices stored as block-diagonal Kronecker products per layer. This allows scalable approximate Bayesian continual learning, overcoming catastrophic forgetting and outperforming diagonal and simple EWC-style schemes.
• Subnetwork Laplace and optimal parameter selection (Raha et al., 9 May 2026): Subnetwork Laplace restricts curvature approximations to a subset of parameters, but always underestimates predictive uncertainty. Gradient-Laplace—selecting parameters whose weights have largest mean squared gradient of the output—and Greedy-Laplace—which sequentially accounts for cross-parameter precision—are theoretically optimal or near-optimal and consistently outperform heuristic diagonal/last-layer schemes in predictive uncertainty and credible interval quality.
• Mixture-of-Laplace Approximations (MoLA) (Eschenhagen et al., 2021): Post hoc combination of Laplace approximations fitted at the MAP of each independently trained model in a deep ensemble; this Gaussian mixture posterior leverages both global multimodality and local curvature, yielding sharper uncertainty—especially off-manifold and under distribution shift—than either classical Laplace or plain ensembles.
• Riemannian Laplace (Bergamin et al., 2023): Samples are drawn along geodesics in a loss-gradient-inducing Riemannian metric, concentrating posteriors in regions of low loss and yielding superior calibration with reduced sensitivity to prior hyperparameters.
• Quadratic Laplace Approximation (QLA) (Jiménez et al., 3 Feb 2026): For deep nets, local Hessian curvature is approximated per data-point via rank-one updates using power iterations on the second derivative, instead of the linearization used in LLA. QLA yields modest, consistent improvements in uncertainty estimation metrics over LLA with minimal extra overhead.
• Learnable Uncertainty via Uncertainty Units (Kristiadi et al., 2020): “Uncertainty units”—hidden units with trainable weights that are inactive at the MAP—can be added to pretrained nets to explicitly increase the marginal likelihood (Laplace evidence) with respect to the Hessian of the posterior. This modular approach improves calibration without impacting point predictions.

5. Improved Laplace for Non-Smooth and Non-Gaussian Problems

Several classes of problems are addressed where standard Laplace breaks down:

• Quantile regression with non-smooth likelihoods (Nava et al., 20 May 2026): The Hessian of the pinball loss is everywhere zero, so Laplace fails. Replacing the observed Hessian with the Fisher information (correct case) or population curvature (misspecification case), possibly estimated nonparametrically via the Triangular Kernel Curvature (TKC) estimator, restores the validity of the Laplace method for mixed effects and GP quantile regression.
• Latent variable models with multimodal or highly skewed marginals (Brown et al., 2019): Integration/sampling over hyperparameters or latent variables is accelerated and made more robust in INLA-type approximations by use of low-discrepancy (quasi-Monte Carlo) sequences combined with polynomial corrections (LDS-CX), achieving full multimodal recovery at a fraction of the computational cost of dense grids.

6. Empirical Benchmarks and Application Domains

Extensive empirical studies across the literature establish that improved Laplace procedures can:

• Drastically increase normalized effective sample size, reduce marginal moment errors, and accurately estimate normalizing constants in synthetic and real Bayesian inference problems, compared to classical Laplace (Bornkamp, 2011, Ruli et al., 2015).
• Substantially improve estimation in high-noise or complex systems identification, such as dynamical systems with discontinuities, high-order derivatives, or rapid growth terms (Laplace-Enhanced SINDy (Zheng et al., 2024)).
• Maintain or improve estimation relative to MCMC and variational methods in high-dimensional, non-Gaussian, or misspecified settings, while yielding dramatic speedups (orders of magnitude) (Nava et al., 20 May 2026, Bergamin et al., 2023).
• Achieve closer calibration and reduced bias in hierarchical regression, spatial GLMMs, and Bayesian penalized spline models (Lambert et al., 2022, Han et al., 2022).

7. Practical Recommendations and Limitations

• Use mixture or iterated Laplace when multimodality or non-ellipticity is evident, especially for multimodal posteriors or heavy-tailed models.
• Apply skew or higher-order corrections when pronounced marginal skewness is present or $\pi$4 is not negligible (Katsevich, 2023, Chiuchiolo et al., 2022).
• In high-dimensional and computational settings, prefer structured, blockwise, or subset Laplace approximations tailored to computational constraints (Ritter et al., 2018, Raha et al., 9 May 2026).
• For quantile/generic non-smooth likelihoods, replace the observed Hessian by Fisher or population curvature, using plug-in or TKC estimation as appropriate (Nava et al., 20 May 2026).
• While improved Laplace methods increase accuracy and robustness, they may incur extra computation or storage, especially as the number of mixture components or dimensions increase. In very high dimensions, careful parallelization and exploitation of structure (sparsity, Kronecker, or gradient selection) are essential.

The proliferation of improved and enhanced Laplace methods has decisively expanded the range of models, applications, and dimensionalities for which deterministic, efficient, and theoretically well-understood approximation is possible, often matching the accuracy of full MCMC at a fraction of the cost (Bornkamp, 2011, Ruli et al., 2015, Katsevich, 2023, Lambert et al., 2022, Chiuchiolo et al., 2022, Raha et al., 9 May 2026).