Gaussian integrals and Laplace's method

Before this lesson: lesson 03 (the Hessian and the second-order Taylor model), lesson 10 (Fubini's theorem, polar coordinates $dx\,dy = r\,dr\,d\theta$, change of variables with $|\det DT|$), linear algebra lesson 06 (the spectral theorem $A = Q\Lambda Q^\top$) and linear algebra lesson 07 (positive definite matrices, Cholesky, $\log\det$). It helps to know the Gaussian distribution and Bayes' rule (prior, likelihood, posterior) from the probability and statistics track; we repeat what we need. Time: about 85 minutes, plus the notebook.

---

Why this matters #

Almost every closed form in probabilistic ML is one Gaussian integral: the normalizing constant of $N(\mu, \Sigma)$, the KL divergence between two Gaussians (in the loss of a variational autoencoder), the marginal likelihood of Bayesian linear regression and of a Gaussian process. When the integral is not Gaussian, two tools remain before you fall back to Monte Carlo.

• Gauss–Hermite quadrature computes a low-dimensional expectation under a Gaussian with a handful of function evaluations. Gaussian-process libraries use it for non-Gaussian likelihoods.
• Laplace's method replaces a peaked integrand by a Gaussian with the same peak and the same curvature. It gives Stirling's formula, the Bayesian information criterion (BIC), and the Laplace approximation, which turns a trained network into a Bayesian one for the price of a Hessian (Laplace Redux, 2021).

In science, the same idea gives uncertainty for inverse problems. In 4D-Var data assimilation (weather forecasting) and in seismic inversion, the inverse Hessian of the cost at its minimum is the covariance of the Laplace (Gaussian) posterior.

1. Predict first #

2. The picture #

One bump: height times width #

The function $e^{-a x^2/2}$ with $a > 0$ has height 1 at $x = 0$. Its curvature at the peak is $a$: the second derivative of the exponent $-ax^2/2$ is $-a$. A large $a$ gives a narrow bump. The width is of order $1/\sqrt{a}$, and the area is exactly $\sqrt{2\pi/a}$ (proved in section 3). Area is height times width: four times the curvature gives half the area.

Two copies make a round bump #

Let $I = \int_{\R} e^{-x^2}\,dx$. No formula built from the usual functions has the derivative $e^{-x^2}$, so the fundamental theorem of calculus does not help directly. The trick is to compute $I^2$ instead:

This is the volume under a round bump in the plane. A round bump is easy to cut into rings. A thin ring at radius $r$ with width $\Delta r$ is almost a long thin rectangle: length $2\pi r$, width $\Delta r$, so its area is about $2\pi r\,\Delta r$. Over the ring, the bump has height $e^{-r^2}$. So the ring holds the volume

The factor $r$ saves us: $2\pi r\,e^{-r^2}$ is the derivative of $-\pi e^{-r^2}$. So the volume of all rings is $\big[-\pi e^{-r^2}\big]_0^\infty = \pi$, and $I = \sqrt\pi$. This factor $r$ is the factor in $dx\,dy = r\,dr\,d\theta$ from lesson 10.

Interactive figure: loading…

Many dimensions: rotate, then factor #

No coupling. If the integrand is a product, the integral is a product (Fubini's theorem, lesson 10):

Here $ac$ is the determinant of $\operatorname{diag}(a, c)$.

Coupling is a rotation. A general quadratic form $x^\top A x$ has tilted elliptical level sets. Rotate the coordinates to the eigenvectors of $A$. A rotation does not change volume, and in the new coordinates there is no coupling. So the volume is $2\pi/\sqrt{\lambda_1\lambda_2} = 2\pi/\sqrt{\det A}$, where $\lambda_1, \lambda_2$ are the eigenvalues of $A$. Along the eigenvector $v_i$ the bump has width of order $1/\sqrt{\lambda_i}$.

A linear term only moves the bump. Adding $b^\top x$ to the exponent moves the peak to $m = A^{-1}b$ and changes its height to $e^{\frac12 b^\top A^{-1}b}$. The width, and so the determinant factor, stays the same.

Interactive figure: loading…

Laplace's idea: a sharp peak looks Gaussian #

Now take an integrand that is not Gaussian: $e^{-n\varphi(x)}$ for a function $\varphi$ with a single minimum at $x^*$ and a large number $n$. At $x^*$ the integrand is largest. Elsewhere it is smaller by the factor $e^{-n(\varphi(x) - \varphi(x^*))}$, which is exponentially small in $n$. So only a small region around $x^*$ matters, of width of order $1/\sqrt{n}$. In that region, $\varphi$ is close to its second-order Taylor model (lesson 03):

There is no linear term, because $\varphi'(x^*) = 0$ at a minimum.

So the integrand is close to a Gaussian bump with the same peak and the same curvature, and we know its area: $e^{-n\varphi(x^*)}\sqrt{2\pi/(n\varphi''(x^*))}$. This is Laplace's method.

Interactive figure: loading…

3. Definitions and main results #

3.1 The Gaussian integral in one dimension #

Theorem (Gaussian integral in 1-D).

Proof.

1. Let $I = \int e^{-x^2}\,dx$. It is finite: from $(|x| - \tfrac12)^2 \ge 0$ we get $x^2 \ge |x| - \tfrac14$, so $e^{-x^2} \le e^{1/4}\,e^{-|x|}$, and $\int e^{-|x|}dx = 2$.
2. $I^2 = \int e^{-x^2}\,dx \int e^{-y^2}\,dy = \iint e^{-(x^2 + y^2)}\,dx\,dy$. This is Tonelli's theorem (the integrand is positive) in the product case (lesson 10, section 3.3).
3. Change to polar coordinates, $x = r\cos\theta$, $y = r\sin\theta$, with $dx\,dy = r\,dr\,d\theta$. The map from $(0, \infty) \times (0, 2\pi)$ covers the plane except one ray, a set of area zero.
4. $I^2 = \int_0^{2\pi}\!\int_0^\infty e^{-r^2}\,r\,dr\,d\theta = 2\pi \cdot \big[-\tfrac12 e^{-r^2}\big]_0^\infty = 2\pi \cdot \tfrac12 = \pi$.
5. $I > 0$, so $I = \sqrt{\pi}$.
6. For (b), substitute $x = \sqrt{2/a}\,t$, so $dx = \sqrt{2/a}\,dt$ and $a x^2/2 = t^2$: $\int e^{-a x^2/2}\,dx = \sqrt{2/a}\,\sqrt{\pi} = \sqrt{2\pi/a}$. $\blacksquare$

In step 4 the factor $r$ from the Jacobian is what makes $r e^{-r^2}$ have an elementary antiderivative. With $a = 1$, (b) says $\int e^{-x^2/2}\,dx = \sqrt{2\pi}$: this is the $\sqrt{2\pi}$ in the standard normal density $e^{-x^2/2}/\sqrt{2\pi}$.

3.2 The Gaussian integral in $d$ dimensions #

A symmetric matrix $A$ is positive definite, written $A \succ 0$, if $x^\top A x > 0$ for every $x \ne 0$. Equivalently, all its eigenvalues are positive (linear algebra lesson 07).

Theorem (Gaussian integral). Let $A$ be a symmetric $d \times d$ matrix with $A \succ 0$, and let $b \in \R^d$. Then

Proof.

1. Complete the square. Let $m = A^{-1} b$. Then
   $$-\tfrac12 x^\top A x + b^\top x = -\tfrac12 (x - m)^\top A (x - m) + \tfrac12 b^\top A^{-1} b .$$
   To check, expand the right side: $-\tfrac12 x^\top A x + \tfrac12 x^\top A m + \tfrac12 m^\top A x - \tfrac12 m^\top A m + \tfrac12 b^\top A^{-1} b$. The two cross terms add up to $x^\top A m = x^\top b$ (this uses $A = A^\top$). The constant $-\tfrac12 m^\top A m = -\tfrac12 b^\top A^{-1} b$ cancels the last term.
2. Shift. Put $y = x - m$. A shift has Jacobian determinant 1, so $dx = dy$.
3. Rotate. By the spectral theorem, $A = Q\Lambda Q^\top$ with $Q$ orthogonal and $\Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_d)$, all $\lambda_i > 0$ because $A \succ 0$. Put $u = Q^\top y$. Then $|\det Q| = 1$, so $dy = du$, and $y^\top A y = u^\top \Lambda u = \sum_i \lambda_i u_i^2$.
4. Factor (Fubini). The integrand is now a product of $d$ one-dimensional bumps: $\int e^{-\frac12 \sum_i \lambda_i u_i^2}\,du = \prod_i \sqrt{2\pi/\lambda_i} = (2\pi)^{d/2} \big(\prod_i \lambda_i\big)^{-1/2} = (2\pi)^{d/2} \det(A)^{-1/2}$, because the determinant is the product of the eigenvalues.
5. Multiply by the constant $e^{\frac12 b^\top A^{-1} b}$ from step 1. $\blacksquare$

Why $\det^{-1/2}$. Each eigen-direction contributes its width $\sqrt{2\pi/\lambda_i}$. The product of the widths is a volume, and the product of the $\lambda_i$ is $\det A$. So a large curvature (a large $\det A$) means a small volume.

The condition $A \succ 0$ is necessary. If one $\lambda_i \le 0$, the $u_i$-integral in step 4 is $\int e^{-\lambda_i u_i^2/2}\,du_i = \infty$. If $A$ is not symmetric, note that $x^\top A x$ only sees the symmetric part $(A + A^\top)/2$: use that matrix in the formula.

In code: Cholesky, not eigenvectors. With $A = LL^\top$ ($L$ lower triangular, linear algebra lesson 07) and $u = L^\top y$, you get $y^\top A y = \|u\|^2$ and $dy = du/\det L$, with $\det L = \det(A)^{1/2}$. So $\log\det A = 2\sum_i \log L_{ii}$, and $m = A^{-1}b$ comes from two triangular solves.

The Gaussian normalizer. Take $A = \Sigma^{-1}$ and $b = 0$, and shift by the mean $\mu$: $\int e^{-\frac12 (x - \mu)^\top \Sigma^{-1}(x - \mu)}\,dx = (2\pi)^{d/2}\det(\Sigma)^{1/2}$. Dividing by this number gives the density of $N(\mu, \Sigma)$:

3.3 Moments and marginal likelihoods from $Z$ #

Fix $A$ and view $Z(b)$ as a function of $b$. The normalized density is $p_b(x) = e^{-\frac12 x^\top A x + b^\top x}/Z(b)$.

1. Mean. Differentiate under the integral sign. This is allowed here, because the integrand and its $b$-derivatives are bounded by integrable Gaussians. $\nabla_b Z = \int x\, e^{-\frac12 x^\top A x + b^\top x}\,dx$, so $\nabla_b \log Z = \nabla_b Z / Z = \mathbb E_{p_b}[x]$.
2. Covariance. Differentiate once more: $\nabla_b^2 \log Z = \mathbb E_{p_b}[x x^\top] - \mathbb E_{p_b}[x]\,\mathbb E_{p_b}[x]^\top = \operatorname{Cov}_{p_b}[x]$.
3. Read them off. From the closed form, $\log Z = \text{const} + \tfrac12 b^\top A^{-1} b$, so $\nabla_b \log Z = A^{-1}b$ and $\nabla_b^2 \log Z = A^{-1}$. So $p_b = N(A^{-1}b, A^{-1})$.

So $A$ is the inverse covariance: the precision matrix. The function $\log Z(b)$ is the log-partition function; "its gradient is the mean and its Hessian is the covariance" holds for every exponential family. For example, the log-sum-exp in the softmax cross-entropy of lesson 06 is the log-partition function of the categorical distribution: its gradient is the vector of softmax probabilities. The same trick in a scale parameter gives a second moment: $\frac{d}{da}\int e^{-a x^2/2}dx = -\frac12\int x^2 e^{-a x^2/2}dx$ and $\frac{d}{da}\sqrt{2\pi}\,a^{-1/2} = -\frac12\sqrt{2\pi}\,a^{-3/2}$, so $\mathbb E[X^2] = 1/a$ for $X \sim N(0, 1/a)$.

Marginal likelihood of Bayesian linear regression. Take $n$ data points with features in a design matrix $\Phi \in \R^{n \times d}$ and targets $y \in \R^n$. The model has the prior $w \sim N(0, \alpha^{-1} I_d)$ and the likelihood $y \mid w \sim N(\Phi w, \sigma^2 I_n)$. The marginal likelihood or evidence is the probability of the data with the weights integrated out: $p(y) = \int p(y \mid w)\,p(w)\,dw$.

1. The exponent of $p(y \mid w)\,p(w)$ is quadratic in $w$: $-\tfrac12 w^\top A w + w^\top b - \|y\|^2/(2\sigma^2)$, with $A = \alpha I + \Phi^\top\Phi/\sigma^2$ and $b = \Phi^\top y/\sigma^2$.
2. The constants in front are $(2\pi\sigma^2)^{-n/2}$ from the likelihood and $(\alpha/2\pi)^{d/2}$ from the prior.
3. The Gaussian integral over $w$ gives $(2\pi)^{d/2}\det(A)^{-1/2}e^{\frac12 b^\top A^{-1} b}$. The two factors $(2\pi)^{\pm d/2}$ cancel, so

The matrix $A$ is the posterior precision, and $m = A^{-1}b$ is the posterior mean (section 3.3, step 3). This value equals $\log N(y; 0, \sigma^2 I + \alpha^{-1}\Phi\Phi^\top)$, the Gaussian-process form; section 5 checks this numerically. The term $-\tfrac12 \log\det A$ penalizes a posterior that is much narrower than the prior: we come back to this "Occam factor" in section 3.7.

3.4 Gauss–Hermite quadrature #

A quadrature rule replaces an integral by a weighted sum of function values. A $K$-node Gauss–Hermite rule is a set of nodes $t_1, \dots, t_K$ and weights $w_1, \dots, w_K > 0$ with

that is exact for every polynomial $g$ of degree $\le 2K - 1$. The nodes are the roots of the Hermite polynomial $H_K$. The count makes sense: $K$ nodes and $K$ weights are $2K$ free numbers, enough to match the $2K$ monomials $1, t, \dots, t^{2K-1}$. NumPy gives them as np.polynomial.hermite.hermgauss(K). This is the "physicists'" weight $e^{-t^2}$, not $e^{-t^2/2}$.

Gaussian expectations. Let $X \sim N(\mu, \sigma^2)$. Substitute $x = \mu + \sqrt{2}\,\sigma t$, so that $(x - \mu)^2/(2\sigma^2) = t^2$ and $dx = \sqrt{2}\,\sigma\,dt$:

with $\tilde w_k = w_k/\sqrt\pi$. Check: with $f = 1$ the sum must be 1, and indeed $\sum_k w_k = \int e^{-t^2}dt = \sqrt{\pi}$.

• For smooth $f$, the error falls very fast with $K$. For $f$ with a kink (ReLU, $\max$, a threshold) it falls slowly, because no polynomial fits a kink well.
• In $D$ dimensions, the tensor-product rule needs $K^D$ nodes. That is fine for $D \le 3$ or so; beyond that, use Monte Carlo (lesson 10) or sparse grids.

Interactive figure: loading…

Why K nodes are exact up to degree 2K − 1

Write $\langle p, q\rangle = \int e^{-t^2} p(t)\,q(t)\,dt$. The Hermite polynomials $H_0, H_1, \dots$ are orthogonal for this inner product, and $H_K$ has degree $K$ and $K$ distinct real roots $t_1, \dots, t_K$.

1. Interpolatory weights. Let $\ell_k$ be the Lagrange polynomial of degree $K - 1$ with $\ell_k(t_j) = 1$ if $j = k$ and $0$ otherwise. Set $w_k = \int e^{-t^2}\ell_k(t)\,dt$. Every polynomial $r$ of degree $\le K - 1$ equals $\sum_k r(t_k)\ell_k$, so the rule is exact for it.
2. Divide by $H_K$. Let $p$ have degree $\le 2K - 1$. Polynomial division gives $p = qH_K + r$ with $\deg q \le K - 1$ and $\deg r \le K - 1$.
3. The $qH_K$ part. The exact integral is $\langle q, H_K\rangle = 0$, because $q$ is a combination of $H_0, \dots, H_{K-1}$. The rule gives $\sum_k w_k\,q(t_k)H_K(t_k) = 0$, because the nodes are roots of $H_K$.
4. The $r$ part is exact by step 1. So the rule is exact for $p$.
5. Degree $2K$ fails. For $p(t) = \prod_k (t - t_k)^2$, the integral is positive, but the rule gives 0.
6. Positive weights. Apply the rule to $\ell_k^2$ (degree $2K - 2$): $w_k = \int e^{-t^2}\ell_k^2\,dt > 0$.

In practice, the nodes and weights come from the Golub–Welsch method: the nodes are the eigenvalues of a small symmetric tridiagonal matrix built from the three-term recurrence of the Hermite polynomials, and the weights come from the first entries of its eigenvectors. The figure above computes them this way.

3.5 Laplace's method #

Theorem (Laplace's method). Let $\varphi: \R^d \to \R$ have a unique global minimizer $x^*$. Assume $\varphi$ has continuous derivatives up to order 4 near $x^*$, with $H = \nabla^2\varphi(x^*) \succ 0$; $\varphi$ stays at least some $\delta > 0$ above $\varphi(x^*)$ outside every ball around $x^*$ (with $\delta$ depending on the ball); and $\int e^{-\varphi} < \infty$. Then, as $n \to \infty$,

Why it is true, in three steps. (1) Factor out the peak value $e^{-n\varphi(x^*)}$. (2) Near $x^*$, replace $\varphi$ by its Taylor model $\varphi(x^*) + \tfrac12 (x - x^*)^\top H (x - x^*)$; far from $x^*$ the integrand is exponentially small. (3) The Gaussian integral of section 3.2 with $A = nH$ and $b = 0$ gives $(2\pi)^{d/2}\det(nH)^{-1/2} = (2\pi/n)^{d/2}\det(H)^{-1/2}$.

Proof

1. Factor out the peak. $\int e^{-n\varphi} = e^{-n\varphi(x^*)} \int e^{-n[\varphi(x) - \varphi(x^*)]}\,dx$.
2. Taylor at the minimum. $\nabla\varphi(x^*) = 0$, so $\varphi(x) - \varphi(x^*) = \tfrac12 (x - x^*)^\top H (x - x^*) + R(x)$ with $|R(x)| \le C\|x - x^*\|^3$ near $x^*$.
3. Zoom in. Put $u = \sqrt{n}\,(x - x^*)$, so $dx = n^{-d/2}\,du$. The integral becomes $n^{-d/2}\int \exp\!\big(-\tfrac12 u^\top H u - n R(x^* + u/\sqrt{n})\big)\,du$.
4. The remainder vanishes. For fixed $u$, $n R(x^* + u/\sqrt{n}) = O(\|u\|^3/\sqrt{n}) \to 0$. Choose a ball around $x^*$ so small that $|R(x)| \le \tfrac14\lambda_{\min}(H)\|x - x^*\|^2$ in it. Inside the ball, the integrand is at most $e^{-\frac14\lambda_{\min}(H)\|u\|^2}$, which is integrable, so the dominated convergence theorem lets the limit pass inside the integral. Outside the ball, the integrand of step 1 is at most $e^{-(n-1)\delta}\,e^{-[\varphi(x) - \varphi(x^*)]}$: exponentially small.
5. Gaussian integral. The limit is $\int e^{-\frac12 u^\top H u}\,du = (2\pi)^{d/2}\det(H)^{-1/2}$. With the factor $n^{-d/2}$ this is $(2\pi/n)^{d/2}\det(H)^{-1/2}$.
6. Why the error is $O(1/n)$. Expand $e^{-nR} \approx 1 - nR + \tfrac12 n^2 R^2$. The cubic part of $nR$ is odd in $u$ and of size $n^{-1/2}$; its integral against the symmetric Gaussian is zero. The next terms, the quartic part of $nR$ and the square of the cubic part, are of size $1/n$. (Steps 1–5 need only three derivatives; this step uses the fourth.) $\blacksquare$

The first correction in 1-D. With $\varphi_k = \varphi^{(k)}(x^*)$, the exact integral is the Laplace value times $1 + \frac{c}{n} + O(n^{-2})$, where

The Laplace figure of section 2 draws $|c|/n$ next to the true error.

Without an explicit $n$. In practice you approximate $\int e^{-\Psi(\theta)}d\theta$ for a function $\Psi$ that grows with the amount of data. The same steps give

where $\hat\theta$ is the minimizer of $\Psi$. The factor $(2\pi/n)^{d/2}$ of the theorem is this formula with $\Psi = n\varphi$, because $\det(nH) = n^d\det H$.

3.6 Stirling's formula #

1. Write $n! = \Gamma(n+1) = \int_0^\infty x^n e^{-x}\,dx = \int_0^\infty e^{n\log x - x}\,dx$.
2. The peak is at $x = n$, so rescale: $x = nt$, $dx = n\,dt$. Then $n\log x - x = n\log n - n(t - \log t)$, and $n! = n^{n+1}\int_0^\infty e^{-n\varphi(t)}\,dt$ with $\varphi(t) = t - \log t$.
3. $\varphi'(t) = 1 - 1/t = 0$ at $t^* = 1$; $\varphi(1) = 1$; $\varphi''(t) = 1/t^2$, so $\varphi''(1) = 1$.
4. Laplace: $\int_0^\infty e^{-n\varphi(t)}dt \approx e^{-n}\sqrt{2\pi/n}$. So
   $$n! \approx n^{n+1} e^{-n}\sqrt{2\pi/n} = \sqrt{2\pi n}\,\Big(\frac{n}{e}\Big)^n .$$
5. With $\varphi_3 = \varphi'''(1) = -2$ and $\varphi_4 = \varphi''''(1) = 6$, the correction formula gives $c = \frac{5 \cdot 4}{24} - \frac{6}{8} = \frac{1}{12}$. So $n! = \sqrt{2\pi n}\,(n/e)^n\,\big(1 + \frac{1}{12n} + O(n^{-2})\big)$.

$n$	$n!$	Stirling	$n!/\text{Stirling} - 1$	$1/(12n)$
1	1	0.9221	0.0844	0.0833
10	3 628 800	3 598 696	0.0084	0.0083
100	$9.333 \cdot 10^{157}$	$9.325 \cdot 10^{157}$	0.00083	0.00083

3.7 The Laplace approximation of a posterior and of the evidence #

Let $\theta \in \R^d$ be the parameters of a model, $p(\theta)$ the prior and $p(D \mid \theta)$ the likelihood of the data $D$. Bayes' rule gives the posterior $p(\theta \mid D) = p(D \mid \theta)\,p(\theta)/p(D)$, where the evidence $p(D) = \int p(D \mid \theta)\,p(\theta)\,d\theta$ is the normalizing constant. Define the negative log joint

In deep learning, $\Psi$ is the training loss with weight decay, up to a constant: if the loss is the summed negative log-likelihood, weight decay $\frac{\gamma}{2}\|\theta\|^2$ is a Gaussian prior with precision $\gamma$.

1. MAP. The maximum a posteriori (MAP) estimate is $\hat\theta = \arg\min_\theta \Psi(\theta)$. You find it with Newton's method, L-BFGS or ordinary training.
2. Curvature. $H = \nabla^2\Psi(\hat\theta) = -\nabla^2\log p(D \mid \hat\theta) - \nabla^2\log p(\hat\theta)$. A prior $N(0, v I)$ adds $I/v$. For logistic regression with inputs $x_i$ (the rows of $X$) and predicted probabilities $p_i$, $H = X^\top \operatorname{diag}\big(p_i(1 - p_i)\big) X + I/v$ (lesson 06).
3. Posterior. $\log p(\theta \mid D) = -\Psi(\theta) + \text{const}$. Its second-order Taylor model at $\hat\theta$ has no linear term, so it is the log of a Gaussian density:
   $$p(\theta \mid D) \approx N(\hat\theta, H^{-1}).$$
   This is the Laplace approximation. The inverse Hessian $H^{-1}$ is its covariance: directions of high curvature are well determined by the data, flat directions are uncertain.
4. Evidence. $p(D) = \int e^{-\Psi(\theta)}\,d\theta$, so Laplace's method gives
   $$\log p(D) \approx \log p(D \mid \hat\theta) + \log p(\hat\theta) + \tfrac{d}{2}\log 2\pi - \tfrac12\log\det H .$$

MacKay reads the evidence as best-fit likelihood × Occam factor. The Occam factor $p(\hat\theta)\,(2\pi)^{d/2}\det(H)^{-1/2}$ is roughly the posterior volume divided by the prior volume. It is small for a flexible model whose data pin down many parameters. So the evidence prefers the simplest model that explains the data, and you can use it to choose hyperparameters, such as the prior variance $v$ (the weight decay), without a validation set.

The BIC. With $n$ independent data points, $H \approx n\,\bar{I}$, where $\bar{I}$ is the average curvature per data point. So $\log\det H = d\log n + \log\det\bar{I}$. Keep only the terms that grow with $n$: $\log p(D) \approx \log p(D \mid \hat\theta) - \tfrac{d}{2}\log n$. The usual form is the Bayesian information criterion $\mathrm{BIC} = -2\log p(D \mid \hat\theta) + d\log n$; smaller is better.

If $\Psi$ is exactly quadratic, steps 3 and 4 are exact. This is the case for a Gaussian prior and a Gaussian likelihood whose mean is linear in $\theta$, such as Bayesian linear regression (section 3.3).

3.8 When the Laplace approximation fails #

1. Several modes (mixtures; the permutations of hidden units in networks). One Gaussian sees one mode and underestimates the uncertainty. Remedy: an ensemble of Laplace approximations, or sampling.
2. Skewed or bounded posteriors (a rate, a variance, a probability near 0 or 1, little data). The Gaussian is centred at the mode, not the mean, puts mass outside the support, and the $O(1/n)$ error is large. Remedy: an unconstrained parameter ($\log\lambda$, $\operatorname{logit} p$). The density then picks up a Jacobian factor (lesson 10), so the result depends on the parameterization.
3. Not a smooth interior minimum. A mode on the boundary, an active constraint, an $L_1$ prior or a ReLU kink at $\hat\theta$ breaks the Taylor model. Early stopping leaves $\nabla\Psi(\hat\theta) \ne 0$.
4. Singular or indefinite Hessian (flat and negative directions of large networks): $H^{-1}$ is not a covariance. Remedy: the generalized Gauss–Newton matrix or the Fisher matrix (both positive semidefinite, lesson 08) plus the prior precision.
5. Size. A $d \times d$ Hessian for $d = 10^7$ cannot be stored. Remedy: diagonal, Kronecker-factored (K-FAC) or low-rank curvature, or Laplace on the last layer only.
6. High dimension. The $O(1/n)$ rate is for a fixed $d$; the correction terms grow with $d$.

Interactive figure: loading…

Answers to section 1. (1) $\pi$. The $y$-direction has curvature 4, so its width halves: $\sqrt{2\pi}\cdot\sqrt{2\pi/4} = \pi$, or $2\pi/\sqrt{\det\operatorname{diag}(1, 4)} = 2\pi/2$. (2) Only about 8%: Stirling gives $\sqrt{2\pi}/e = 0.9221$ instead of $1! = 1$. The ratio $n!/\text{Stirling}$ is about $1 + 1/(12n)$, so the error is 0.83% at $n = 10$ and 0.083% at $n = 100$. Laplace's method is often good long before $n$ is large. (3) Below. The Gaussian is centred at the mode, and for a distribution skewed to the right, the mode is below the mean. In worked example (c) the mode is 1.5 and the mean is 1.75.

4. Worked example #

(a) A 2-D Gaussian integral. Compute $Z$ for $A = \begin{pmatrix} 3 & 1 \\ 1 & 2 \end{pmatrix}$ and $b = (1, 1)$.

1. $A \succ 0$: the trace is $5 > 0$ and $\det A = 6 - 1 = 5 > 0$, so both eigenvalues are positive.
2. $A^{-1} = \tfrac15\begin{pmatrix} 2 & -1 \\ -1 & 3 \end{pmatrix}$, $m = A^{-1}b = \tfrac15(2 - 1, -1 + 3) = (0.2, 0.4)$, and $b^\top A^{-1} b = b^\top m = 0.2 + 0.4 = 0.6$.
3. $Z = (2\pi)^{2/2}\cdot 5^{-1/2}\cdot e^{0.3} = 2.8099 \times 1.3499 = 3.7930$, and $\log Z = 1.8379 - 0.8047 + 0.3 = 1.3332$.
4. The normalized density is $N(m, A^{-1})$: mean $(0.2, 0.4)$ and covariance $\begin{pmatrix} 0.4 & -0.2 \\ -0.2 & 0.6 \end{pmatrix}$. Both come from $\log Z(b) = \text{const} + \tfrac12 b^\top A^{-1} b$.

(The "worked example" button of the second figure shows this bump.)

(b) Gauss–Hermite by hand. Compute $\mathbb E[e^X]$ for $X \sim N(0, 1)$. The exact value is $e^{1/2} = 1.6487$: it is $Z(1, 1)/Z(1, 0)$ in 1-D, that is $e^{\frac12 \cdot 1^2}$.

• $K = 2$: hermgauss(2) gives $t = \pm 1/\sqrt2$ and $w = \sqrt{\pi}/2$ each. The mapped nodes are $x = \sqrt2\,t = \pm 1$, with weights $\tilde w = w/\sqrt\pi = \tfrac12$. So $\mathbb E[e^X] \approx \tfrac12(e + e^{-1}) = 1.5431$: 6.4% too low.
• $K = 3$: $t = 0, \pm\sqrt{3/2}$ with $w = 2\sqrt\pi/3, \sqrt\pi/6, \sqrt\pi/6$. Mapped: $x = 0, \pm\sqrt3$ with weights $\tfrac23, \tfrac16, \tfrac16$. So $\mathbb E[e^X] \approx \tfrac23 + \tfrac16(5.6522 + 0.1769) = 1.6382$: 0.64% too low.
• $K = 5$ gives $1.648679$ (relative error $2.5 \cdot 10^{-5}$), and $K = 10$ agrees with $e^{1/2}$ to about 12 digits.

The 2-node rule is the "sigma point" rule $\mu \pm \sigma$. It is exact for polynomials up to degree 3.

(c) Laplace for a Poisson rate. Three equal time windows show $y = (2, 1, 3)$ events. Model: $y_i \mid \lambda \sim \mathrm{Poisson}(\lambda)$, so $p(y_i \mid \lambda) = \lambda^{y_i}e^{-\lambda}/y_i!$, and the prior is $p(\lambda) = e^{-\lambda}$ on $\lambda > 0$.

1. Joint. $p(y, \lambda) = \prod_i \frac{\lambda^{y_i}e^{-\lambda}}{y_i!}\cdot e^{-\lambda} = \frac{\lambda^6 e^{-4\lambda}}{12}$, because $\sum_i y_i = 6$ and $2!\,1!\,3! = 12$.
2. Exact answers. $p(y) = \frac{1}{12}\int_0^\infty \lambda^6 e^{-4\lambda}d\lambda = \frac{6!}{12 \cdot 4^7} = 0.0036621$ (substitute $x = 4\lambda$ and use $\int_0^\infty x^6 e^{-x}dx = 6!$), so $\log p(y) = -5.6097$. The posterior is $\propto \lambda^6 e^{-4\lambda}$: a Gamma distribution with shape 7 and rate 4, with mean $7/4 = 1.75$ and standard deviation $\sqrt7/4 = 0.661$.
3. MAP. $\Psi(\lambda) = 4\lambda - 6\log\lambda + \log 12$. $\Psi'(\lambda) = 4 - 6/\lambda = 0$ gives $\hat\lambda = 1.5$.
4. Curvature. $H = \Psi''(\hat\lambda) = 6/\hat\lambda^2 = 2.667$. Laplace posterior: $N(1.5,\ 0.375)$, standard deviation $\sqrt{0.375} = 0.612$.
5. Evidence. $\log p(y \mid \hat\lambda) = 6\log 1.5 - 3 \cdot 1.5 - \log 12 = -4.5521$, $\log p(\hat\lambda) = -1.5$, and $\tfrac12\log(2\pi/H) = 0.4285$. Sum: $\log p(y) \approx -5.6236$, so $p(y) \approx 0.0036116$: 1.38% too low.
6. Check of the error. The 1-D correction of section 3.5 (with $\Psi$ in place of $n\varphi$, so $n = 1$) uses $\Psi_2 = 6/\hat\lambda^2 = 8/3$, $\Psi_3 = -12/\hat\lambda^3 = -32/9$ and $\Psi_4 = 36/\hat\lambda^4 = 64/9$. It gives $c = \frac{5 (32/9)^2}{24 (8/3)^3} - \frac{64/9}{8 (8/3)^2} = \frac{5}{36} - \frac18 = \frac{1}{72} = 1.39\%$, as observed.

The defects of section 3.8 are visible. The Laplace mean is the mode 1.5, not the mean 1.75, because the posterior is skewed. The Laplace standard deviation (0.612) is too small (0.661). And the Gaussian puts 0.7% of its mass on $\lambda < 0$. This is the "$y = (2, 1, 3)$" setting of the last figure.

5. In code #

NumPy and SciPy check the numbers above with numerical integration and library densities.

import math
import numpy as np
from scipy import integrate, stats

# (a) sqrt(pi), and the 2-D example against the hand value 2*pi/sqrt(5)*e^0.3
print(round(integrate.quad(lambda x: np.exp(-x * x), -np.inf, np.inf)[0], 6), round(math.sqrt(math.pi), 6))
f = lambda y, x: np.exp(-0.5 * (3 * x * x + 2 * x * y + 2 * y * y) + x + y)
Z = integrate.dblquad(f, -15, 15, -15, 15)[0]
mean_x = integrate.dblquad(lambda y, x: x * f(y, x), -15, 15, -15, 15)[0] / Z
print(round(Z, 4), round(2 * math.pi / math.sqrt(5) * math.exp(0.3), 4), round(mean_x, 4))
# 1.772454 1.772454
# 3.793 3.793 0.2

# Bayesian linear regression: the evidence (an integral over w) is the density of y
rng = np.random.default_rng(0)
Phi, y, alpha, s2 = rng.normal(size=(5, 2)), rng.normal(size=5), 2.0, 0.5
def joint(w2, w1):
    r = y - Phi @ np.array([w1, w2])
    lik = np.exp(-r @ r / (2 * s2)) / (2 * math.pi * s2) ** 2.5
    return lik * alpha / (2 * math.pi) * np.exp(-alpha * (w1 * w1 + w2 * w2) / 2)
evidence = integrate.dblquad(joint, -10, 10, -10, 10, epsabs=1e-14)[0]
A, b = alpha * np.eye(2) + Phi.T @ Phi / s2, Phi.T @ y / s2
formula = (-5 / 2 * math.log(2 * math.pi * s2) + math.log(alpha) - 0.5 * np.linalg.slogdet(A)[1]
           - y @ y / (2 * s2) + 0.5 * b @ np.linalg.solve(A, b))
gp_form = stats.multivariate_normal(np.zeros(5), s2 * np.eye(5) + Phi @ Phi.T / alpha).logpdf(y)
print(round(math.log(evidence), 6), round(formula, 6), round(gp_form, 6))
# -10.580234 -10.580234 -10.580234

# (b) Gauss-Hermite nodes and weights after x = sqrt(2) t, for N(0, 1)
for K in (2, 3):
    t, w = np.polynomial.hermite.hermgauss(K)
    print(np.round(np.sqrt(2) * t, 4), np.round(w / np.sqrt(np.pi), 4))
# [-1.  1.] [0.5 0.5]
# [-1.7321  0.      1.7321] [0.1667 0.6667 0.1667]
exact = integrate.quad(lambda x: np.exp(x) * stats.norm.pdf(x), -30, 30)[0]     # = e^0.5
print(round(math.cosh(1), 4), round(2 / 3 + math.cosh(math.sqrt(3)) / 3, 4), round(exact, 4))
# 1.5431 1.6382 1.6487

import math
from scipy import special, stats

# Stirling from Laplace: n!/Stirling - 1 is close to 1/(12n)
for n in (1, 10, 100):
    log_stirling = 0.5 * math.log(2 * math.pi * n) + n * math.log(n) - n
    print(n, round(math.expm1(math.lgamma(n + 1) - log_stirling), 6), round(1 / (12 * n), 6))
# 1 0.084438 0.083333
# 10 0.008365 0.008333
# 100 0.000834 0.000833

# (c) the Poisson rate: exact evidence and posterior against the Laplace numbers of section 4
log_ev_exact = special.gammaln(7) - 7 * math.log(4) - math.log(12)
log_ev_laplace = (6 * math.log(1.5) - 4 * 1.5 - math.log(12)) + 0.5 * math.log(2 * math.pi / (6 / 1.5**2))
post = stats.gamma(a=7, scale=1 / 4)
print(round(log_ev_exact, 4), round(log_ev_laplace, 4))              # -5.6097 -5.6236
print(post.mean(), round(post.std(), 4), round(math.sqrt(1.5**2 / 6), 4))   # 1.75 0.6614 0.6124
print(round(stats.norm(1.5, math.sqrt(0.375)).cdf(0), 4))            # 0.0072: mass on lambda < 0

For a large $d$, never compute det and then log: the determinant overflows or underflows. Use np.linalg.slogdet, or $2\sum_i \log L_{ii}$ from a Cholesky factor, and solve with $A$ instead of forming $A^{-1}$.

6. Where it shows up #

• Gaussian log-densities. $\log N(x; \mu, \Sigma) = -\tfrac12(x - \mu)^\top\Sigma^{-1}(x - \mu) - \tfrac12\log\det\Sigma - \tfrac{d}{2}\log 2\pi$ is in every Gaussian likelihood and in the KL term of a variational autoencoder. It is evaluated with a Cholesky factor of $\Sigma$.
• Marginal likelihood of Gaussian processes. Gaussian-process regression tunes the kernel hyperparameters by maximizing $\log p(y) = -\tfrac12 y^\top K_y^{-1} y - \tfrac12\log\det K_y - \tfrac{n}{2}\log 2\pi$, where $K_y$ is the kernel matrix plus the noise variance times $I$. The $-\tfrac12\log\det K_y$ term penalizes kernels that are too flexible.
• Gauss–Hermite in practice. GPflow and GPyTorch compute the 1-D expectations $\mathbb E_{q(f)}[\log p(y \mid f)]$ of variational Gaussian-process classification by Gauss–Hermite quadrature. For a Bayesian logistic model whose score is $X \sim N(1, 2^2)$, $\mathbb E[\sigma(X)] = 0.6477$; the probit approximation $\sigma\big(\mu/\sqrt{1 + \pi\sigma^2/8}\big)$ gives 0.6511, and the plug-in $\sigma(\mu)$ gives 0.7311 (overconfident). In 1-D, the unscented transform of the unscented Kalman filter with $\kappa = 2$ is exactly the 3-node rule (nodes $\mu, \mu \pm \sqrt3\,\sigma$, weights $\tfrac23, \tfrac16, \tfrac16$).
• Bayesian deep learning. Laplace on the last layer or with K-FAC curvature, with the generalized Gauss–Newton matrix in place of the Hessian, the evidence to tune the weight decay, and a probit or linearized predictive (Laplace Redux and the laplace-torch library). How to compute curvature products without forming the matrix: lesson 08.
• ML for science. In 4D-Var data assimilation and in seismic full-waveform inversion, the minimizer of the cost is the MAP estimate, and the inverse Gauss–Newton Hessian there is the covariance of the Laplace posterior (exact for a linear forward model). INLA (integrated nested Laplace approximations) uses nested Laplace approximations for latent Gaussian fields in spatial statistics.
• Model selection and Stirling. The BIC chooses the number of components of a Gaussian mixture (see the ML fundamentals track, coming later). PyTorch's poisson_nll_loss(..., full=True) adds the Stirling term $y\log y - y + \tfrac12\log(2\pi y)$ for $\log y!$.

7. Common mistakes #

• Forgetting the $r$ in $r\,dr\,d\theta$. Without it, the ring areas are wrong and the integral $\int_0^\infty e^{-r^2}dr$ has no elementary antiderivative.
• Precision and covariance mixed up. $\int e^{-\frac12 x^\top A x}dx$ has $\det(A)^{-1/2}$, and $A$ is the precision. The Gaussian normalizer has $\det(\Sigma)^{+1/2}$ in the integral and $\det(\Sigma)^{-1/2}$ in the density.
• A non-symmetric $A$. $x^\top A x$ only sees $(A + A^\top)/2$. Symmetrize before you use $\det A$ and $A^{-1}b$.
• Forgetting the square root. Each direction contributes $\sqrt{2\pi/\lambda_i}$, so the integral has $\det(A)^{-1/2}$, not $\det(A)^{-1}$.
• Gauss–Hermite constants. hermgauss uses the weight $e^{-t^2}$. You need both the $\sqrt2$ in $x = \mu + \sqrt2\,\sigma t$ and the factor $1/\sqrt\pi$. Check that the weights sum to 1 for $f = 1$.
• Gauss–Hermite for a kink. ReLU, $\max$ and thresholds converge slowly; split the integral at the kink or use more nodes.
• Dropping the prior in $H$. The posterior precision is the likelihood Hessian plus the prior precision. Without the prior term, $H$ can be singular and the evidence is wrong.
• Treating Laplace as parameterization-free. The mode, the Gaussian and the evidence estimate all change when you change variables, for example from $\lambda$ to $\log\lambda$.
• Evidence and likelihood mixed up. The maximum likelihood $p(D \mid \hat\theta)$ always grows with model size; the evidence $p(D) = \int p(D \mid \theta)p(\theta)\,d\theta$ includes the Occam factor and does not.

8. Check yourself #

1. Compute $\int_{\R^3}\exp\big(-\tfrac12(x^2 + 2y^2 + 8z^2)\big)\,dx\,dy\,dz$.
   Answer

   $A = \operatorname{diag}(1, 2, 8)$, $\det A = 16$, so the value is $(2\pi)^{3/2}/\sqrt{16} = 15.7496/4 = 3.9374$.
2. Let $p(x) \propto \exp(-\tfrac a2 x^2 + bx)$ with $a > 0$. Find its mean and variance from $\log Z(b)$.
   Answer

   $\log Z(b) = \tfrac12\log(2\pi/a) + \frac{b^2}{2a}$. Its derivative $b/a$ is the mean, and its second derivative $1/a$ is the variance: $p = N(b/a, 1/a)$.
3. The 3-node Gauss–Hermite rule uses $x = 0, \pm\sqrt3$ with weights $\tfrac23, \tfrac16, \tfrac16$. How many nodes make $\mathbb E[X^6]$ exact for $X \sim N(0, 1)$, and why?
   Answer

   $K$ nodes are exact up to degree $2K - 1$. You need $2K - 1 \ge 6$, so $K = 4$. The 3-node rule gives $2 \cdot \tfrac16 \cdot 27 = 9$ instead of 15.
4. Redo the evidence of worked example (c) with Laplace in $\theta = \log\lambda$. Remember the Jacobian.
   Answer

   $\lambda = e^\theta$, $d\lambda = e^\theta d\theta$, so the integrand in $\theta$ is $e^{7\theta - 4e^\theta}/12$. The mode is $e^{\hat\theta} = 7/4$; the curvature is $4e^{\hat\theta} = 7$. Laplace: $1.75^7 e^{-7}/12 \cdot \sqrt{2\pi/7} = 0.0038197 \times 0.9474 = 0.0036188$, 1.18% too low. This is better than 1.38% in $\lambda$, and it puts no mass on invalid values. Two parameterizations give two answers.
5. With $n = 1000$ data points, model 1 has $d = 3$ parameters and maximum log-likelihood $-520.0$; model 2 has $d = 5$ and $-515.0$. Which one does the BIC prefer?
   Answer

   $\log p(D) \approx \hat\ell - \tfrac{d}{2}\log n$ with $\log 1000 = 6.9078$. Model 1: $-520.0 - 10.36 = -530.36$. Model 2: $-515.0 - 17.27 = -532.27$. So model 1: the gain of 5 nats does not pay for two more parameters, which cost $\log 1000 = 6.9$ nats.
6. Why does a Laplace approximation of a deep network usually use the generalized Gauss–Newton matrix instead of the true Hessian?
   Answer

   At a trained network, $\nabla^2\Psi$ is often indefinite (SGD does not stop at an exact minimum, and there are flat and negative directions), so its inverse is not a covariance. The GGN $J^\top\nabla^2\ell\,J$ is positive semidefinite for convex losses; adding the prior precision makes it positive definite; and it has cheap matrix-free and Kronecker-factored forms.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is $\int_{\R} e^{-x^2}\,dx$, and what is the key step of the proof?

   Show answer

2. Give $\int_{\R^d}\exp(-\tfrac12 x^\top A x + b^\top x)\,dx$ for symmetric $A \succ 0$.

   Show answer

3. What are the mean and covariance of the density $\propto \exp(-\tfrac12 x^\top A x + b^\top x)$?

   Show answer

4. Why is $\det(A)^{-1/2}$ in the Gaussian integral?

   Show answer

5. How do you compute $\mathbb E[f(X)]$, $X \sim N(\mu, \sigma^2)$, with Gauss–Hermite nodes $t_k$ and weights $w_k$ (weight $e^{-t^2}$)?

   Show answer

6. For which functions is a $K$-node Gauss–Hermite rule exact?

   Show answer

7. State Laplace's method for $\int_{\R^d} e^{-n\varphi(x)}\,dx$.

   Show answer

8. State Stirling's formula and the size of its relative error.

   Show answer

9. What is the Laplace approximation of a posterior $p(\theta \mid D)$?

   Show answer

10. Why does the Laplace approximation fail for a bimodal posterior?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Compute $\int e^{-x^2}dx$." — Square it, use polar coordinates and the factor $r$: $\sqrt\pi$.
• "Where does $\det(\Sigma)^{1/2}$ in the Gaussian normalizer come from?" — Rotate to the eigenvectors; the integral factors into $d$ 1-D integrals with widths $\sqrt{2\pi\lambda_i(\Sigma)}$, and their product is $(2\pi)^{d/2}\det(\Sigma)^{1/2}$.
• "Write the log marginal likelihood of a GP and say how you compute it." — Cholesky of $K_y$: $\log\det K_y = 2\sum_i\log L_{ii}$, and two triangular solves for $y^\top K_y^{-1} y$.
• "What is the Laplace approximation? When is it exact, and when does it fail?" — A Gaussian at the MAP with covariance $H^{-1}$; exact for linear-Gaussian models; it fails for multimodal, skewed or bounded posteriors, kinks, and indefinite Hessians.
• "Where does the BIC come from?" — From the Laplace evidence, keeping only the terms that grow with $n$.
• "How do you get cheap uncertainty for a trained network?" — Last-layer or K-FAC Laplace with the GGN, the prior precision tuned on the evidence, and a probit or linearized predictive.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. Write log_gaussian_integral(A, b) with a Cholesky factor, and check it against dblquad and in $d = 200$, where det overflows.
2. Write gauss_hermite_expectation(f, mu, sigma, K), and compare its error with Monte Carlo for a smooth function and for ReLU.
3. Write laplace_log_integral with Newton's method, and check it on Stirling's formula and on an exact Gaussian.
4. Write laplace_logistic_evidence for Bayesian logistic regression on the breast-cancer data, check it against importance sampling, and choose the prior variance by the evidence.
5. Compute the Laplace evidence of the Poisson-rate example in $\lambda$ and in $\log\lambda$, and show that the error depends only on the total count.

11. Watch and read #

Why π is in the normal distribution (beyond integral tricks)3Blue1Brown · 25 min

The volume under the round 2-D bump, cut into rings, and why this gives √π.

Laplace's Method (Quadratic Approximation)Meerkat Statistics · 16 min

Laplace's method step by step: second-order Taylor expansion at the peak, then the Gaussian integral.

4.4 The Laplace Approximation - Pattern Recognition and Machine LearningSina Tootoonian · Bishop's PRML, section 4.4 · 13 min

The Laplace approximation of a distribution and of its normalizing constant, as in Bishop's book.

BDA 2019 Lecture 11.1 Normal approximation, Laplace approximation.Aki Vehtari · Bayesian Data Analysis course, Aalto University · 44 min

When a normal approximation of a posterior works, how the parameterization matters, and how to diagnose it.

NeurIPS 2021: Laplace Redux -- Effortless Bayesian Deep LearningTübingen Machine Learning · 11 min

Laplace approximations for neural networks: which weights, which curvature (GGN, K-FAC) and the evidence for hyperparameters.

Preview: The Magic of Gaussian Quadrature - A Billion Times Better than the Next Best ThingMathTheBeautiful · 10 min

Why choosing the nodes as well as the weights makes Gaussian quadrature so accurate.

• MacKay, Information Theory, Inference, and Learning Algorithms, free PDF: ch. 27 "Laplace's Method" (exercise 27.1 is a Poisson rate in $\lambda$ and in $\log\lambda$) and §28.1 (evidence ≈ best-fit likelihood × Occam factor).
• Deisenroth, Faisal and Ong, Mathematics for Machine Learning, free PDF: §6.5 (Gaussian distribution), §8.6 (Bayesian model selection) and §9.3.5 (marginal likelihood of Bayesian linear regression).
• Rasmussen and Williams, Gaussian Processes for Machine Learning, free PDF: §2.2 (the marginal likelihood as a Gaussian integral) and §3.4 (the Laplace approximation for Gaussian-process classification).
• Gelman et al., Bayesian Data Analysis, 3rd ed., free PDF: §4.1 (normal approximations to the posterior), §13.1 (finding posterior modes) and §13.3 (normal and related mixture approximations).
• Daxberger et al., Laplace Redux — Effortless Bayesian Deep Learning (NeurIPS 2021): §2 (the Laplace posterior and evidence), then subsets of weights, curvature approximations and predictives.
• NIST Digital Library of Mathematical Functions, §3.5(v): Gauss quadrature, including the Gauss–Hermite rule.
• SciPy, Integration tutorial: "General multiple integration (dblquad, tplquad, nquad)" and "Gaussian quadrature".