Integrated Laplace Approximation (ILA)

• ILA is a family of approximation methods that employs local Gaussian expansions around dominant modes to marginalize latent variables efficiently in complex models.
• It transforms high-dimensional integrals into tractable approximations by leveraging conditional modes and curvature, facilitating hyperparameter inference in latent Gaussian settings.
• Recent refinements incorporate adjoint differentiation, importance sampling corrections, and nested strategies (INLA) to enhance accuracy and computational performance.

Integrated Laplace Approximation (ILA) denotes a family of approximation strategies for replacing an intractable integral by a local Gaussian expansion around a dominant mode and then integrating out latent variables, nuisance parameters, or intermediate states. In latent Gaussian models, the central use of ILA is to marginalize a high-dimensional Gaussian latent field so that inference proceeds on hyperparameters and posterior marginals rather than on the full joint posterior. In spatial statistics and econometrics, this role is often realized through Integrated Nested Laplace Approximation (INLA), which targets marginal posterior distributions of latent components and hyperparameters by combining Laplace approximations with numerical integration over a low-dimensional hyperparameter space (Gomez-Rubio et al., 2017, Gómez-Rubio et al., 2019, Margossian, 2023, Lai et al., 19 May 2026).

1. Terminological scope and major variants

The literature surveyed here suggests that the label “Integrated Laplace Approximation” is used heterogeneously. In some papers it refers directly to approximate marginalization of a latent Gaussian variable in a latent Gaussian model; in others, closely related methodology appears under the more specific label INLA, emphasizing the nested structure of the approximation; and in yet other settings, Laplace-based integration is adapted to path integrals, marginal likelihoods, or bespoke model classes (Gómez-Rubio et al., 2019, Ruli et al., 2015, Thygesen, 27 Mar 2025).

Variant	Integrated quantity	Characteristic setting
ILA in latent Gaussian models	Latent Gaussian variable $z$ or $\theta$	Hyperparameter inference after approximate marginalization
INLA	Latent marginals and hyperparameters via nested approximations	Latent GMRF models and R-INLA workflows
Improved Laplace approximation	Normalizing constant $I=\int e^{-h(x)}\,dx$	General marginal likelihoods and related integrals
Continuous-time Laplace approximation	Intermediate states or Brownian paths	Transition densities for stochastic differential equations

For latent Gaussian models, the basic hierarchy is written as

$$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$

or equivalently

$$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$

The exact marginal density after integrating out the latent field is

$$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$

When the likelihood is Gaussian, this integral is exact; for non-Gaussian likelihoods, ILA replaces it by a Laplace approximation built from the conditional mode and curvature of the latent posterior (Margossian, 2023, Lai et al., 19 May 2026).

A distinct but related line of work studies general multidimensional integrals of the form

$I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$

including Bayesian marginal likelihoods, posterior normalizing constants, and GLMM random-effects integrals. There the “integrated” idea appears through sequential approximation of normalized marginal and conditional densities, followed by re-normalization (Ruli et al., 2015).

2. Core mathematical construction in latent Gaussian models

In the latent Gaussian setting, ILA begins from the conditional posterior of the latent field. The conditional mode is

$$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$

and the negative Hessian of the log likelihood is

$$W \triangleq -\nabla_\theta^2 \log \pi(y\mid \theta,\eta).$$

The approximate marginal likelihood is then obtained by Laplace approximation around $\hat\theta$. One formulation writes

$\theta$0

with

$\theta$1

Using Woodbury identities, the determinant term can be rewritten as

$\theta$2

so the same approximation can be expressed through a convenient $\theta$3-matrix factorization (Margossian, 2023).

This marginalization is motivated by posterior geometry. The hierarchical prior in latent Gaussian models induces a posterior geometry that is “prone to frustrate inference algorithms”; marginalizing out the latent Gaussian variable by ILA “removes the offending geometry,” making gradient-based inference on hyperparameters substantially more tractable (Margossian, 2023). A similar motivation appears in corrected ILA work: direct inference on $\theta$4 is difficult because $\theta$5 is often large and strongly coupled with $\theta$6, so the usual strategy is to integrate out $\theta$7 and perform inference on $\theta$8 (Lai et al., 19 May 2026).

A convenient expression for standard ILA in this setting is the ratio

$\theta$9

evaluated at the posterior mode $I=\int e^{-h(x)}\,dx$0, where $I=\int e^{-h(x)}\,dx$1 is a Gaussian approximation constructed from the mode and curvature of the conditional posterior (Lai et al., 19 May 2026). This shows directly where the approximation enters: the denominator is only an approximation to the true conditional posterior.

3. Nested Laplace approximation and conditional averaging

INLA is the most developed nested version of this idea for latent Gaussian Markov random fields. Its generic model class is

$I=\int e^{-h(x)}\,dx$2

with latent effects

$I=\int e^{-h(x)}\,dx$3

and hyperparameters $I=\int e^{-h(x)}\,dx$4. The crucial computational feature is that the latent field is a GMRF, so the associated precision matrix is sparse (Gómez-Rubio et al., 2019).

INLA does not primarily target the full joint posterior. Instead, it targets the posterior marginals

$I=\int e^{-h(x)}\,dx$5

The approximation proceeds in two nested stages. First, INLA computes an approximation $I=\int e^{-h(x)}\,dx$6 to the hyperparameter posterior. It then recovers the latent marginals by numerical integration: $I=\int e^{-h(x)}\,dx$7 approximated by a weighted sum over grid points $I=\int e^{-h(x)}\,dx$8: $I=\int e^{-h(x)}\,dx$9 The conditional marginals $$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$0 are themselves obtained via Laplace-type approximations (Gomez-Rubio et al., 2017).

A further extension conditions on a subset of hyperparameters. If

$$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$1

one can fit the model with INLA conditional on $$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$2, obtaining $$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$3 and $$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$4. Bayesian model averaging then restores uncertainty in $$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$5 through

$$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$6

approximated on a grid by

$$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$7

with weights

$$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$8

This construction extends INLA beyond “pure” latent GMRF models, but it works best when the dimension of $$\phi \sim \pi(\phi), \qquad \theta \sim \mathrm{Normal}(0,K(\phi)), \qquad y \sim \pi(y \mid \theta,\eta),$$9 is low (Gómez-Rubio et al., 2019).

4. Algorithmic refinements, error correction, and higher-order variants

A major recent direction replaces hand-derived derivatives and restrictive Hessian assumptions by more general differentiation and linear-algebra machinery. The general adjoint-differentiated Laplace approximation supports three Newton solvers unified by $$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$0-matrix formulations, works with block-diagonal Hessians rather than only diagonal ones, uses automatic differentiation instead of analytical derivatives of the likelihood, and reuses the $$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$1-matrix factorization already computed in the Newton solve. The method organizes gradients of the approximate marginal likelihood into three components: an explicit derivative term, a log-determinant derivative term, and an implicit term from the dependence of the latent mode $$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$2 on the hyperparameters (Margossian, 2023).

This refinement addresses explicit limitations of earlier adjoint-differentiated Laplace methods: the diagonal Hessian restriction, the need for analytical first-, second-, and third-order likelihood derivatives, and limited support for non-log-concave or unconventional likelihoods. Numerical experiments reported in the same work indicate that, on a standard latent Gaussian model, the new method is slightly faster than the earlier adjoint-differentiated Laplace approximation, while remaining applicable to a broader class of likelihoods (Margossian, 2023).

A separate refinement corrects the approximation error of standard ILA by reinterpreting it as importance sampling. If

$$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$3

then

$$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$4

is an unbiased estimator of $$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$5, since

$$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$6

With $$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$7 samples from the Laplace approximation, the corrected estimator remains unbiased and has variance $$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$8. This observation leads to pseudo-marginal, quasi-Monte Carlo, and randomized quasi-Monte Carlo corrected schemes. PM-ADLA and RQMC-ADLA are asymptotically exact for $$\theta \sim \pi(\theta), \qquad z \sim \mathrm{normal}(0,K(\theta)), \qquad y \sim \pi(y\mid \theta,z).$$9; QMC-ADLA converges to the true posterior in total variation under a domination condition, but finite-$$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$0 bias can still be substantial (Lai et al., 19 May 2026).

Another branch of the literature pursues higher-order accuracy for general integrals. The improved Laplace approximation starts from the identity

$$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$1

chooses $$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$2, factorizes the normalized density into marginal and conditional densities, approximates each factor by Laplace-type formulas, and then re-normalizes numerically. The resulting approximation is

$$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$3

Under the usual regularity assumptions,

$$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$4

whereas the standard Laplace approximation has relative error $$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$5. The paper therefore describes the method as having third-order accuracy (Ruli et al., 2015).

5. Representative model classes and disciplinary uses

The practical impact of ILA-type methods is most visible in settings where repeated high-dimensional integration or repeated model fitting would make MCMC or dense Gaussian-process computation expensive. The following examples are representative rather than exhaustive.

Domain	Contribution	Paper
Spatial econometrics	slm latent class in R-INLA for spatial lag / SAR models	(Gomez-Rubio et al., 2017)
Conditional latent GMRF spatial econometrics	BMA with INLA for SAC models	(Gómez-Rubio et al., 2019)
Dirichlet regression	Gaussian pseudo-observations and dirinla	(Martínez-Minaya et al., 2019)
Spatial hierarchical modelling	LGCP-SPDE and BYM inference with INLA	(D'Angelo et al., 2020)
Health economic EVPPI	SPDE-INLA with low-dimensional projection	(Heath et al., 2015)
Stochastic differential equations	Continuous-time Laplace approximation for transition densities	(Thygesen, 27 Mar 2025)
Structural equation models	INLA-inspired profiling with skew-normal marginals	(Jamil et al., 26 Mar 2026)

In spatial econometrics, the standard spatial lag model

$$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$6

or equivalently

$$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$7

was implemented in R-INLA through a new latent model class called slm. A central step is to cast the model as a GMRF with sparse precision matrix, which makes it compatible with INLA’s sparse-matrix machinery. This implementation makes spatial econometrics models native to R-INLA and enables model fitting, model selection, posterior marginals for coefficients and spatial dependence, marginal likelihood, DIC and CPO, variable selection, prediction of missing responses, and Bayesian model averaging across alternative spatial structures (Gomez-Rubio et al., 2017).

The conditional-averaging extension appears in a spatial autoregressive combined model,

$$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$8

where $$\pi(\theta,y)=\pi(\theta)\int \pi(z\mid \theta)\pi(y\mid \theta,z)\,dz.$$9 are treated as the conditioning hyperparameters. Once $I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$0 are fixed, the model becomes a Gaussian mixed model with known precision structure, so INLA can fit it conditionally and BMA can average over the grid of conditioning values (Gómez-Rubio et al., 2019).

For compositional data, base R-INLA cannot directly fit Dirichlet likelihoods with more than two categories. The dirinla approach approximates the Dirichlet log-likelihood by a quadratic expansion in the linear predictor, defines Gaussian pseudo-data

$I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$1

and obtains the pseudo-observation model

$I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$2

Stacking all observations converts the multivariate Dirichlet likelihood into a Gaussian observation model that R-INLA can handle. In the Arabidopsis thaliana case study, the reported computation times were about $I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$3 s for dirinla, about $I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$4 s for standard JAGS, and about $I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$5 s for long JAGS (Martínez-Minaya et al., 2019).

In spatial hierarchical modelling more broadly, INLA is used for latent Gaussian fields in both continuous and discrete spatial domains. The surveyed applications include an inhomogeneous Log-Gaussian Cox Process with SPDE representation for a seismic sequence in Greece and a Besag–York–Mollié disease-mapping model for Covid-19 infection in North Italy. In both cases, the posterior is not available in closed form because of non-Gaussian observation models and latent spatial effects, so inference proceeds through latent Gaussian modelling, sparse precision matrices, and nested Laplace approximations (D'Angelo et al., 2020).

In health economics, SPDE-INLA has been used to accelerate Gaussian-process regression for Expected Value of Partial Perfect Information. There the computational gain comes from two approximations: replacing dense GP covariance matrices by a sparse SPDE/GMRF representation and projecting a high-dimensional parameter subset into a low-dimensional, ideally two-dimensional, space by Principal Fitted Components. Reported timings show standard GP regression between $I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$6 and $I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$7 seconds in the Vaccine example and $I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$8 to $I=\int_{\mathbb R^d} e^{-h(x)}\,dx,$9 seconds in the SAVI example, versus about $$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$0 to $$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$1 seconds and $$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$2 to $$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$3 seconds, respectively, for SPDE-INLA; the method was often up to about $$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$4 times faster (Heath et al., 2015).

Outside latent Gaussian hierarchical models, closely related integrated Laplace constructions appear in path-space problems. For stochastic differential equations, a discrete-time approximation introduces intermediate states and Brownian increments, then applies Laplace approximation to the resulting high-dimensional integral. In the continuous-time limit, the dominant path becomes the solution of a minimum-energy optimal control problem, and the final approximation is expressed through canonical equations, a Riccati equation, and a Lyapunov equation (Thygesen, 27 Mar 2025). In Bayesian structural equation modelling, an INLA-inspired method integrates out latent factors exactly in the linear-Gaussian case, then uses a joint Laplace approximation, deterministic marginal profiling, a volume correction, skew-normal marginal fitting, and a variational Bayes location shift (Jamil et al., 26 Mar 2026).

6. Accuracy, failure modes, and methodological boundaries

The strengths of ILA-type methods are tied to identifiable boundaries. INLA is fundamentally a method for latent Gaussian models, and its strongest guarantees concern accurate approximation of marginal posterior quantities rather than the full joint posterior. Conditional INLA plus BMA extends the scope, but numerical integration over the conditioning parameters is limited to low dimensions and may not scale well (Gomez-Rubio et al., 2017, Gómez-Rubio et al., 2019).

Approximation quality remains model dependent. In latent Gaussian models, Laplace approximation may fail badly for multimodal or strongly non-Gaussian conditional posteriors of the latent field, and there is no cheap universal diagnostic. Suggested diagnostics include importance sampling, PSIS, leave-one-out cross-validation, and simulation-based calibration (Margossian, 2023). Corrected ILA addresses this by replacing the one-sample Laplace ratio by multi-sample importance estimators, but the trade-off is computational: PM is exact but can be expensive, QMC can still be biased at finite $$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$5, and RQMC offers a compromise between exactness and variance reduction (Lai et al., 19 May 2026).

Application-specific caveats are also recurrent. In spatial econometrics, the new slm class makes SAR-type models native to R-INLA, but more complex models with multiple spatial parameters or nonstandard likelihoods may still require conditional modelling strategies, approximation, or model averaging. The reported discussion also emphasizes practical numerical care: covariates may need rescaling, lagged covariates can be highly collinear, and binary or other highly parameterized likelihoods require attention to identifiability (Gomez-Rubio et al., 2017). In the Dirichlet case, the approximation is local and depends on the posterior mode, while hyperparameter approximation becomes harder in larger latent structures (Martínez-Minaya et al., 2019).

Path-space approximations display a related boundary in a different form. For the double-well stochastic differential equation,

$$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$6

the most probable path can be trivial even though the true transition density receives substantial contributions from trajectories that go to one well and return. This shows that non-near-path contributions can impair the approximation when the path integral is effectively multimodal. Numerical results for the Cox–Ingersoll–Ross process also indicate that intrinsic Laplace error scales approximately like

$$\hat\theta = \arg\max_\theta \pi(\theta\mid y,\phi,\eta),$$7

which is consistent with improvement in the weak-noise or short-time regime and deterioration as trajectories wander farther from the minimum-effort path (Thygesen, 27 Mar 2025).

Taken together, these results delineate a coherent methodological picture. ILA is most effective when the posterior or integrand is dominated by a single mode, when local curvature captures most of the relevant geometry, and when sparse or low-dimensional structure makes deterministic approximation preferable to full sampling. Where those conditions weaken, the contemporary literature increasingly supplements the basic Laplace construction with conditional averaging, automatic differentiation, importance-sampling correction, re-normalization, or path-space control formulations rather than treating the standard approximation as universally sufficient (Ruli et al., 2015, Margossian, 2023, Lai et al., 19 May 2026).