The Laplace approximation, and where it stops working

“Bayesian logistic regression” sounds expensive. In BEACON it costs one extra matrix inverse over the fit you were already doing, because of an approximation that is almost two and a half centuries old and still the right first move.

It is also an approximation with a specific shape of failure, and the failure is easier to reason about than most people assume.

The problem

You have a posterior over parameters:

p(θ | D) ∝ p(D | θ) p(θ)

For logistic regression with a Gaussian prior, the normalising constant has no closed form. You cannot integrate it, so you cannot directly compute the things you want — posterior means, credible intervals, predictive distributions.

The options are to sample it (MCMC, expensive), to optimise a surrogate for it (variational inference, more machinery), or to approximate it locally. The third is Laplace’s.

The move

Work in log space, where the posterior is a smooth function with a single peak for this model. Take a second-order Taylor expansion of the log posterior about its mode θ*, the MAP estimate:

log p(θ | D)  ≈  log p(θ* | D)  −  ½ (θ − θ*)ᵀ H (θ − θ*)

The first-order term vanishes, which is the whole trick: at a mode the gradient is zero, so there is no linear term to carry. What remains is a constant plus a quadratic form, with H the Hessian of the negative log posterior evaluated at the mode.

Exponentiate that and you have the kernel of a Gaussian:

p(θ | D)  ≈  N(θ*, H⁻¹)

A quadratic in log space is a Gaussian in probability space. So the claim being made is precise and modest: near its peak, this posterior curves like a Gaussian. The inverse Hessian is the covariance because curvature and uncertainty are the same quantity — a sharply peaked log posterior is a confident one.

What it looks like for logistic regression

With a Gaussian prior N(0, σ²I), the negative log posterior is

L(θ) = Σᵢ log(1 + exp(−yᵢ xᵢᵀθ))  +  ‖θ‖² / (2σ²)

The first term is the logistic loss and is convex. The second is a quadratic penalty and is strictly convex. Their sum is strictly convex, which means one mode, no local optima, and the Taylor expansion has an unambiguous point to expand about.

Differentiate twice and the Hessian is

H = XᵀSX + I/σ²        where  S = diag( pᵢ(1 − pᵢ) )

with pᵢ the fitted probability for observation i. That diagonal is worth staring at. The weight an observation contributes to the curvature is p(1−p), maximised at p = 0.5 and vanishing as predictions approach 0 or 1.

Confidently classified points contribute almost nothing to the posterior precision. All the information about parameter uncertainty comes from points near the decision boundary. That is not an artefact of the approximation; it is a real property of the likelihood that the approximation makes visible.

It also has a direct consequence for rare-event problems. If your positives are few and the model is confident about most of your negatives, the effective sample size governing your parameter uncertainty is much smaller than your row count suggests.

Where it stops working

The approximation is local and symmetric. Everything it gets wrong follows from those two words.

Skew. A Gaussian is symmetric about its mode. A skewed posterior is not, and Laplace will match the peak while misplacing the mass. Under separation — when a feature perfectly splits the classes — the logistic posterior develops a long tail toward large coefficients that a symmetric approximation cannot represent.

Multimodality. Expanding about one mode ignores the others entirely. Logistic regression with a proper prior is safe here by convexity, but the moment you approximate something non-convex — a mixture, a neural network, anything with symmetries — you are describing one basin and calling it the distribution.

Mode is not mean. The approximation is centred at the MAP estimate, which for a skewed posterior is not the posterior mean. If you want a mean, Laplace hands you the wrong summary confidently.

Boundaries. Parameters constrained to an interval — variances, probabilities, anything non-negative — have posteriors that pile up against the constraint. A Gaussian puts mass outside the feasible region. Reparameterise onto an unbounded scale before approximating, or expect nonsense.

Small data. The justification is asymptotic: posteriors tend toward Gaussian as data accumulates. With a dozen positives, that argument has not started working yet.

Why BEACON uses it anyway

The conditions here are the favourable ones. The posterior is strictly convex and unimodal, the parameters are unbounded, and the model is a baseline rather than the headline result. Laplace gives a genuine posterior at roughly the cost of the fit, with failure modes that are understood rather than mysterious.

Its job here is not to win; the gradient boosting ensemble ranks better. It exists so the study contains one genuine posterior to hold the ensemble’s Bayesian-inspired spread up against — a reference point rather than a contender.

Which is the useful framing. An approximation with known, checkable failure modes is worth more as a reference point than a better score from a method whose uncertainty nobody can interrogate.