Laplace Approximation From Local Gaussian Curvature

From a Marginal Integral to Local Curvature

In a nonlinear mixed-effects model, a subject’s marginal likelihood contribution averages the conditional data likelihood over that subject’s random effects. This averaging is an integral, and it is often difficult to evaluate exactly when the individual predictions depend nonlinearly on the random effects. The Laplace approximation replaces the integrand near its dominant region with a Gaussian function whose width is determined by local curvature.

The key idea is local: find the conditional mode, measure how sharply the log integrand falls away from that mode, and integrate the resulting quadratic approximation. This tutorial derives that step. It does not make the approximation exact for a general nonlinear model; its reliability depends on how well the local quadratic shape represents the part of the integrand that contributes appreciably to the integral.

Definition: For subject \(i\), let \(\eta_i\) be a \(q\)-dimensional random-effect vector and let \(L_i(\eta_i)=p(y_i\mid \eta_i,\theta)p(\eta_i\mid\theta)\) be the conditional likelihood multiplied by the random-effect density. The marginal likelihood contribution is \(L_i^{\mathrm{marg}}(\theta)=\int_{\mathbb{R}^q}L_i(\eta_i)\,d\eta_i\). Write \(g_i(\eta_i)=\log L_i(\eta_i)\). An interior conditional mode \(\hat{\eta}_i\) is a point maximizing \(g_i\), and its negative Hessian is \(H_i=-\nabla^2g_i(\hat{\eta}_i)\).

The random-effect density belongs in the integrand. In particular, when \(\eta_i\sim N(0,\Omega)\), its log density contributes both a quadratic penalty and a normalization term to \(g_i\). Omitting the density would approximate a different integral. The earlier material on marginal likelihoods and conditional random-effect modes establishes this integral and the role of the mode; the step here is to approximate its integral using the mode’s curvature.

Deriving the Local Gaussian Function

Assume \(g_i\) is twice continuously differentiable in a neighborhood of an interior mode \(\hat{\eta}_i\). Because the mode is interior and differentiable, its gradient is zero. A second-order Taylor expansion, with \(h=\eta_i-\hat{\eta}_i\), gives

g_i(\hat{\eta}_i+h) = g_i(\hat{\eta}_i) +\frac{1}{2}h^\mathsf{T}\nabla^2g_i(\hat{\eta}_i)h +o(\|h\|^2) = g_i(\hat{\eta}_i)-\frac{1}{2}h^\mathsf{T}H_i h+o(\|h\|^2).

The remainder is \(o(\|h\|^2)\): divided by \(\|h\|^2\), it tends to zero as \(h\) tends to zero. This statement alone does not guarantee that the quadratic term is negative in every direction. For that conclusion, assume \(H_i\) is positive definite. Then \(h^\mathsf{T}H_i h>0\) for every nonzero \(h\), so the quadratic approximation decreases away from the mode.

Exponentiating the quadratic part gives a constant times a Gaussian-shaped function:

L_i(\hat{\eta}_i+h) \ \approx\ \exp\{g_i(\hat{\eta}_i)\} \exp\left(-\frac{1}{2}h^\mathsf{T}H_i h\right).

If \(H_i\) is positive definite, the integral of this quadratic approximation over all \(q\) dimensions is finite. Diagonalizing \(H_i\), or changing variables to \(z=H_i^{1/2}h\), transforms the quadratic form into \(z^\mathsf{T}z\). The change of variables has volume factor \(|H_i|^{-1/2}\), where \(|H_i|\) is the determinant. The standard \(q\)-dimensional Gaussian integral is \((2\pi)^{q/2}\), yielding the Laplace approximation:

L_i^{\mathrm{marg}}(\theta) \ \approx\ \exp\{g_i(\hat{\eta}_i)\} (2\pi)^{q/2}|H_i|^{-1/2}.

The exponential term is the integrand’s height at the mode. The remaining factor is its local Gaussian volume: flatter curvature means a wider region and a larger volume, while sharper curvature means a smaller volume. For a log integrand that is exactly quadratic over the integration domain, this Gaussian integral is exact. Otherwise, it is a local approximation. Its usual asymptotic justification is strongest when the integrand becomes increasingly concentrated around a regular, isolated mode, while the mode dimension remains fixed.

Worked Calculations

Worked Example: A Scalar Integrand With a Quartic Tail

Consider a hypothetical scalar integrand with log \(g(\eta)=-1-\frac{9}{2}(\eta-0.2)^2-0.2(\eta-0.2)^4\). It is integrable, and its mode is at \(\hat{\eta}=0.2\): the terms after \(-1\) are nonpositive and vanish only at that point. The quartic term makes the log integrand nonquadratic away from the mode.

Set \(h=\eta-0.2\). The second derivative at the mode is \(-9\), so the negative curvature is \(H=9>0\). The local Gaussian approximation is therefore

\int_{-\infty}^{\infty}\exp\{g(\eta)\}\,d\eta \ \approx\ \exp(-1)\sqrt{\frac{2\pi}{9}} = 0.3679(0.8355) \approx 0.3074.

This is an approximation because the quartic term is not included in the local quadratic function. The conditions for the local calculation hold: the mode is interior, the log integrand is smooth, the integral is finite, and the negative curvature is strictly positive. The calculation does not, by itself, show how close 0.3074 is to the exact integral.

Worked Example: One Nonlinear Random Effect

Suppose a hypothetical subject has one observation \(y=0\), conditional model \(Y\mid\eta\sim N(\eta+0.2\eta^2,0.5^2)\), and random effect \(\eta\sim N(0,1)\). Both densities are included in \(L(\eta)\). With \(f(\eta)=\eta+0.2\eta^2\), the log integrand is

g(\eta) = -\log(\pi)-\frac{\eta^2}{2}-2f(\eta)^2.

The normalizing constant is \(-\log(\pi)\), since the likelihood and prior density constants multiply to \(1/\pi\). The exponent terms are nonpositive and can both equal zero only at \(\eta=0\), so \(\hat{\eta}=0\) is the global interior mode. Differentiating \(f\) gives \(f'(0)=1\). The second derivative of the log integrand at zero is \(-1-4\{f'(0)\}^2=-5\), and thus \(H=5>0\).

Here \(q=1\), so the approximation is

L^{\mathrm{marg}}(\theta) \ \approx\ \exp\{-\log(\pi)\}\sqrt{\frac{2\pi}{5}} = \frac{1}{\pi}(1.1210) \approx 0.3568.

This is an approximate marginal density for the observation, in inverse units of \(y\). The mode, smoothness, finite integral, and positive curvature conditions are all satisfied. Because the conditional mean is nonlinear, the log integrand is not globally quadratic, so the result need not equal the exact marginal density.

Worked Example: Two Random Effects and the Determinant

For a hypothetical two-dimensional random effect, suppose a smooth, integrable log integrand has mode \(\hat{\eta}\), value \(g(\hat{\eta})=-2\), and local expansion \(g(\hat{\eta}+h)=-2-\frac12h^\mathsf{T}Hh-0.1(h_1^4+h_2^4)\), where

H= \begin{pmatrix} 4&1\\ 1&3 \end{pmatrix}.

The determinant is \(|H|=(4)(3)-(1)(1)=11>0\), and the first diagonal entry is positive; these conditions make this symmetric matrix positive definite. Thus the quadratic curvature is positive in every nonzero direction. The quartic terms are zero with zero first and second derivatives at the mode, so they do not change its value or curvature there.

Substitute \(q=2\), \(g(\hat{\eta})=-2\), and \(|H|=11\):

L^{\mathrm{marg}}(\theta) \ \approx\ \exp(-2)\frac{2\pi}{\sqrt{11}} = 0.1353\frac{6.2832}{3.3166} \approx 0.2564.

The determinant accounts for the combined spread and orientation of the local Gaussian in both dimensions. Replacing it with a product of the diagonal entries would give \(12\), not \(11\), and would ignore the off-diagonal curvature.

When the Approximation Is Trustworthy

Positive-definite curvature is necessary for this local Gaussian integral, but it does not guarantee that the local region captures most of the full integral. For example, a broad shoulder, a distant region of substantial integrand, or strong asymmetry can make the mass away from the mode important. The approximation uses only the height and second-order local shape; it does not directly account for those features.

A useful mathematical justification for the local shape follows from positive definiteness. Let \(\lambda_{\min}>0\) be the smallest eigenvalue of \(H_i\). Then \(h^\mathsf{T}H_i h\geq\lambda_{\min}\|h\|^2\). Since the Taylor remainder is \(o(\|h\|^2)\), sufficiently close to the mode its absolute value is at most \((\lambda_{\min}/4)\|h\|^2\). Consequently,

-\frac12h^\mathsf{T}H_i h+o(\|h\|^2) \ \leq\ -\frac{\lambda_{\min}}{2}\|h\|^2 +\frac{\lambda_{\min}}{4}\|h\|^2 = -\frac{\lambda_{\min}}{4}\|h\|^2.

Thus, near the mode, the log integrand falls at least quadratically under these conditions. This is a local statement; it does not establish that the tail or other regions of the full integrand are negligible.

In practice, warning signs include a mode that changes substantially with numerical starting values, a nearly singular curvature matrix, or a calculated mode close to a parameter-space boundary. These can indicate optimization or approximation problems. A nearly singular \(H_i\) implies a very large Gaussian-volume factor, and a non-positive-definite \(H_i\) does not define the required finite Gaussian integral. Neither case should be handled by blindly substituting a determinant into the formula.

Laplace approximation differs from first-order and FOCE approximations in what it uses to approximate the subject contribution. First-order methods linearize individual predictions around a reference random effect; FOCE uses a conditional mode and conditional linearization. Laplace uses the log integrand’s local second-order curvature at the mode and integrates the resulting Gaussian. Exact integration is preferable when available, as in the linear-Gaussian benchmark, but nonlinear mixed-effects models often do not permit it. Laplace can be effective when the integrand is concentrated and locally well represented by a Gaussian; it is not a guarantee against bias or poor fit when the integrand is skewed, broad, or otherwise poorly represented near its mode.

Common Mistakes and Practical Checks

Using only the conditional likelihood. The marginal integral includes the random-effect density as well as the conditional data likelihood. Leaving out the prior changes both the mode and the curvature.
Treating a strict mode as proof of positive curvature. A strict maximum can have a singular Hessian. For example, \(g(x)=-x^4\) has a strict maximum at zero but zero second derivative there. Check positive definiteness rather than inferring it from the word “maximum.”
Forgetting the Gaussian volume. The approximation is not just \(\exp\{g_i(\hat{\eta}_i)\}\); it also includes \((2\pi)^{q/2}|H_i|^{-1/2}\).
Assuming local accuracy implies global accuracy. Curvature describes the integrand near the mode. It does not rule out important contributions elsewhere.

Key takeaway: Under smoothness, an interior mode, and positive-definite negative curvature, a second-order expansion produces a local Gaussian whose integral is \(\exp\{g_i(\hat{\eta}_i)\}(2\pi)^{q/2}|H_i|^{-1/2}\). The formula is an approximation unless the log integrand is quadratic, and local conditions alone do not ensure global accuracy.

Check Your Understanding

Try each question, then compare with the solutions below.

Why must the random-effect density be part of the subject-level integrand?
Suppose \(q=1\), \(g(\hat{\eta})=-1\), and \(H=9\). Calculate the Laplace approximation.
Does a strict interior maximum alone guarantee that the negative Hessian is positive definite? Explain.
For a two-dimensional positive-definite curvature matrix with determinant 16 and \(g(\hat{\eta})=0\), write the Laplace approximation.

Solutions. (1) The marginal likelihood integrates the conditional data likelihood against the random-effect distribution; omitting that density changes the integral being approximated. (2) The approximation is \(\exp(-1)\sqrt{2\pi/9}=0.3679(0.8355)\approx 0.3074\). (3) No. A strict maximum may have a singular Hessian, as \(g(x)=-x^4\) at zero demonstrates; positive definiteness must be checked separately. (4) With \(q=2\), the result is \(\exp(0)(2\pi)/\sqrt{16}=2\pi/4=\pi/2\approx1.571\).