DR-003: EP site approximations: keep EP.logit's rank-one Gaussian sites, identity start and undamped sweeps

Status: Accepted
Date recorded: 2026-10-06
Decision in one line: The expectation-propagation port keeps every site choice of the original EP.logit (rank-one Gaussian sites along each design row, identity-precision starting sites, sequential undamped updates, tilted moments by adaptive quadrature) because tests show the answer does not depend on the start or on damping, and EP lands within 0.4% of the best KL divergence any Gaussian can achieve.

Context

EP.logit, adapted from the MAST90125 lecture code, approximates the posterior by a product of Gaussian "sites", one per observation. Several choices sit inside that function, and the 2023 report did not discuss any of them.

Site family. Each site is Gaussian in natural parameters, and its precision is rank one along the design row $x_i$ . The likelihood depends on β only through $\eta_i = x_i^\top \beta$ , so nothing is lost by this restriction.
Starting sites. Every site starts with precision equal to the identity and zero shift, so the first sweep begins from a Gaussian with precision $N I$ .
Update schedule. Sites are updated one at a time in data order, with no damping.
Tilted moments. The mean and variance of each one-dimensional tilted distribution come from three integrate() calls over the cavity mean ± 10 cavity sd.
Prior. The prior is flat, so there is no prior site. The seven likelihood sites alone give a full-rank precision.
Stopping rule. Sweeps stop when the largest relative change in the global natural parameters is below $10^{-6}$ , or after 100 sweeps.

EP gives no general guarantee of convergence. Undamped EP can oscillate, and sites can acquire negative precision, so each of these choices could plausibly matter.

Decision

Keep all of them unchanged in the port (which matches R to about $10^{-12}$ ), and test the two that most often cause trouble in practice, the starting sites and the lack of damping. Optional initSitePrecision and damping arguments were added for the tests, with defaults that reproduce the original exactly.

Options considered

Keep the original choices (chosen). This preserves parity with the 2023 output, provided the evidence shows the choices are harmless here.
Damp every update (for example by 0.5). This is the standard safeguard against oscillation, but it slows convergence and changes the 2023 iteration counts.
Start from flat sites (zero precision). This is a more principled start, but the first cavity would be improper for some sites with a flat prior, so it needs a small positive start anyway.
Replace adaptive quadrature with Gauss-Hermite quadrature or a probit approximation of the logistic. This is faster, but the port deliberately mirrors R's integrate() routine for routine, because EP's iteration count depends on its exact tolerance behaviour.

Why

The evidence below shows that, for this model, the fixed point is the same whichever way the sites start, damping only slows the iteration down, and the Gaussian EP finds is as good as a Gaussian can be. Changing the choices would cost parity with the original and buy nothing measurable. The logistic likelihood is log-concave, so every tilted distribution is narrower than its cavity and every site precision stays positive (the test suite checks this at every update).

What happened

Starting sites do not matter. Starting every site at 0.01 I or at 100 I instead of I gives the same answer to within $6 \times 10^{-10}$ , in the same 5 sweeps.
Damping is unnecessary. Damping each update by 0.5 needs 23 sweeps instead of 5 and lands within $2.8 \times 10^{-7}$ of the undamped answer, inside the $10^{-6}$ stopping tolerance.
EP is as close to the exact posterior as a Gaussian can be. By quadrature on a 401 × 1,601 grid, the KL divergence from the exact posterior to EP's Gaussian is 0.0027. The best possible Gaussian, the one with the exact mean and covariance, also scores 0.0027, so EP is within 0.4% of it. The Laplace approximation scores 0.0085, about three times worse.
The weak point is shape, not location. EP's slope mean (0.3142) is within 0.0002 of the exact mean (0.3143), and its sd is 0.3% smaller. The exact posterior of the slope is right-skewed, however, and no Gaussian can represent that. EP's symmetric 95% interval [0.056, 0.573] sits 0.007 to 0.009 to the left of the exact interval [0.063, 0.582] at its two ends.

What I'd change

Report the tilted-moment integrals' accuracy and cost. Three adaptive integrations per site per sweep are far more than a smooth one-dimensional integrand needs, and 20-point Gauss-Hermite quadrature would do.
Make damping and the starting sites explicit arguments in the original function, with a sentence on why the defaults are safe, rather than leaving them implicit in the code.
Treat the prior as a site even when it is flat, so the same code handles weakly informative priors (see the prior-sensitivity analysis on /diagnostics).
When skewness matters for the conclusion, as it does for interval endpoints here, report a skew-aware approximation or a sampler alongside EP.