1. Big picture: what EP is doing

Goal: approximate the posterior $p(\theta \mid y)$

with something simpler, say a multivariate normal $g(\theta)$.

EP is:

  • A deterministic iterative algorithm (no MCMC randomness).
  • It works with a factorization over data points:
    • $p(\theta \mid y) \propto f(\theta) = \prod_{i=0}^n f_i(\theta)$
      • Think:
        • $f_0(\theta)$: prior.
        • $f_i(\theta)$ for $i \ge 1$: likelihood contribution from data point $i$.

EP approximates each $f_i(\theta)$ by a “site” factor $g_i(\theta)$ in a convenient parametric family (usually exponential family, often Normal). Then the global approximation is $g(\theta) = \prod_{i=0}^n g_i(\theta).$

At convergence, each $g_i(\theta)$ is supposed to be a good approximation to $f_i(\theta)$ in the context of all the other data (through the other factors).

Key difference vs variational Bayes (mean-field):

  • Variational Bayes often factorizes by parameters, e.g. $g(\theta) = \prod_{j=1}^J g_j(\theta_j).$
  • EP factorizes by data / likelihood terms, e.g. $g(\theta) = \prod_{i=0}^n g_i(\theta).$

So VB tries to make parameters independent; EP keeps parameters coupled but approximates each likelihood factor.

2. Exponential family, sufficient statistics, and natural parameters

EP usually chooses $g(\theta)$ from an exponential family. That means:

$g(\theta) = \frac{1}{Z(\eta)} \exp\{ \eta^\top T(\theta) \}$

  • $T(\theta)$: vector of sufficient statistics,
  • $\eta$: natural parameters,
  • $Z(\eta)$: normalizing constant.

Example: univariate normal $N(\theta \mid \mu, \sigma^2)$:

$N(\theta \mid \mu, \sigma^2) \propto \exp\left( -\frac{1}{2\sigma^2}(\theta – \mu)^2 \right) = \exp\left( -\frac{1}{2\sigma^2}\theta^2 + \frac{\mu}{\sigma^2}\theta \right)$

  • Sufficient statistics: $T(\theta) = (\theta, \theta^2)$
  • Natural parameters:
    • coeff of $\theta$: $\mu / \sigma^2$,
    • coeff of $\theta^2$: $-1/(2\sigma^2)$.

For a multivariate normal $N(\theta \mid \mu, \Sigma)$:

  • Sufficient statistics: $\theta$ and $\theta\theta^\top$.
  • Natural parameters:
    • vector: $\Sigma^{-1}\mu$,
    • matrix: $\Sigma^{-1}$.

EP works by matching expectations of sufficient statistics under a “tilted” distribution, so knowing these is important.

3. EP: cavity distribution and tilted distribution

We approximate the posterior by $g(\theta) = \prod_{i=0}^n g_i(\theta)$,

with each $g_i(\theta)$ in the chosen family.

At each EP update step, pick one index $i$ (one site) and:

  1. Remove the current site $g_i(\theta)$ from the global approximation to get the cavity distribution:
    • $g_{-i}(\theta) \propto \frac{g(\theta)}{g_i(\theta)}$
    • This is “what we think about θ using all data except data point i.”
  2. Insert the true factor $f_i(\theta)$ back in to get the tilted distribution:
    • $\text{tilted}_i(\theta) \propto g_{-i}(\theta) f_i(\theta)$
    • This is the “local true target” mixing:
      • our approximation from other data, and
      • the exact contribution of data point $i$.
  3. Approximate the tilted distribution by a new parametric distribution in the chosen family, by moment matching:
    • compute expectations of sufficient statistics $T(\theta)$ under the tilted distribution,
    • choose new parameters of $g(\theta)$ so that its sufficient-statistic expectations match those of the tilted distribution.
  4. From this updated global approximate $g(\theta)$, derive the updated site factor:
    • $g_i(\theta) = \frac{g(\theta)}{g_{-i}(\theta)}$

Repeat this for all i (cycling) until convergence.

Two variants:

  • Sequential EP: update $g(\theta)$ after each site update.
  • Parallel EP: compute all new site updates first, then update global $g(\theta)$ once; usually faster.

4. Moment matching: the core step

Suppose the global approximation is multivariate normal:

$g(\theta) = N(\theta \mid \mu, \Sigma)$

Sufficient statistics: $\theta$ and $\theta\theta^\top$.

In one EP update step for site i:

  1. We have cavity distribution:
    • $g_{-i}(\theta) = N(\theta \mid \mu_{-i}, \Sigma_{-i})$ (in practice we track this via natural parameters, not directly).
  2. Tilted distribution: $\tilde{p}_i(\theta) \propto g_{-i}(\theta) f_i(\theta)$
  3. Moment matching means:
    • Set the new global approximation $g^{\text{new}}(\theta)$ to be the normal distribution with:
      • $\mu^{\text{new}} = E_{\tilde{p}_i}[\theta], \quad \Sigma^{\text{new}} = \text{Var}_{\tilde{p}_i}(\theta)$
  4. Then compute the new site factor as:
    • $g_i^{\text{new}}(\theta) = \frac{g^{\text{new}}(\theta)}{g_{-i}(\theta)}$

So each EP step is:

  • compute the “true” (tilted) one-step posterior for that site,
  • find the best normal approximation by matching first and second moments,
  • update the site so that the product of all sites has those moments.

The hard part is computing $E_{\tilde{p}_i}[\theta]$ and $\text{Var}_{\tilde{p}_i}(\theta)$. But usually we exploit structure to reduce this to low-dimensional integrals.

5. Logistic regression example (general, with θ)

Model:

  • Data: $y_i \sim \text{Binomial}(m_i, \pi_i), \quad \pi_i = \text{logit}^{-1}(X_i \theta)$
    • where:
      • $X_i$​ is the i-th row of the $n \times k$ design matrix,
      • $\theta \in \mathbb{R}^k$ are regression coefficients.
  • Prior: $\theta \sim N(\mu_0, \Sigma_0)$

Posterior is not conjugate, but we can approximate it by a multivariate normal with EP.

5.1 Factorization

Write the unnormalized posterior as:

$p(\theta \mid y) \propto f(\theta) = \underbrace{f_0(\theta)}_{\text{prior}} \prod_{i=1}^n \underbrace{f_i(\theta)}_{\text{likelihood for } y_i}$

where:

  • $f_0(\theta) = N(\theta \mid \mu_0, \Sigma_0)$
  • $f_i(\theta) = \text{Bin}(y_i \mid m_i, \text{logit}^{-1}(X_i \theta))$

EP approximation:

$g(\theta) = \prod_{i=0}^n g_i(\theta)$

each $g_i(\theta)$ a multivariate normal factor. In practice:

  • $g_0(\theta)$ is just the prior (so $g_0 = f_0$​).
  • For each $i \ge 1$, $g_i(\theta)$ is a normal factor with its own natural parameters.

5.2 Cavity distribution

The global approximation is normal:

$g(\theta) = N(\theta \mid \mu, \Sigma)$

with natural parameters:

  • $\Sigma^{-1} \mu$ (vector),
  • $\Sigma^{-1}$ (matrix).

The i-th site factor is also normal:

$g_i(\theta) = N(\theta \mid \mu_i, \Sigma_i)$

with natural parameters:

  • $\Sigma_i^{-1} \mu_i$
  • $\Sigma_i^{-1}$

To get the cavity distribution $g_{-i}(\theta)$:

  • In natural parameters, $\Sigma_{-i}^{-1} = \Sigma^{-1} – \Sigma_i^{-1}$​, $\Sigma_{-i}^{-1}\mu_{-i} = \Sigma^{-1} \mu – \Sigma_i^{-1}\mu_i.$

So $g_{-i}(\theta) = N(\theta \mid \mu_{-i}, \Sigma_{-i})$.

5.3 Reducing to 1D: projection onto the linear predictor

For data point $i$, the likelihood factor depends on $\theta$ only through the linear predictor:

$\eta_i = X_i \theta \quad (\text{scalar})$

Under the cavity distribution

$N(\theta \mid \mu_{-i}, \Sigma_{-i})$, $\eta_i = X_i \theta \sim N(M_{-i}, V_{-i})$,

with:

  • $M_{-i} = X_i \mu_{-i}$
  • $V_{-i} = X_i \Sigma_{-i} X_i^\top$

Thus the tilted distribution in terms of $\eta_i$ is:

$\text{tilted}_i(\eta_i) \propto N(\eta_i \mid M_{-i}, V_{-i}) \times \text{Bin}(y_i \mid m_i, \text{logit}^{-1}(\eta_i))$

This is a one-dimensional distribution over $\eta_i$​. Therefore:

  • We can compute:
    • $E[\eta_i]$
    • $E[\eta_i^2]$
  • using 1D numerical integration (e.g. Gauss–Kronrod quadrature).

Let:

  • $E_0 = \int \tilde{p}_i(\eta_i)\, d\eta_i$ (normalizing constant),
  • $E_1 = \int \eta_i \tilde{p}_i(\eta_i)\, d\eta_i$
  • $E_2 = \int \eta_i^2 \tilde{p}_i(\eta_i)\, d\eta_i$

where

$\tilde{p}_i(\eta_i) = N(\eta_i \mid M_{-i}, V_{-i}) \times \text{Bin}(y_i \mid m_i, \text{logit}^{-1}(\eta_i))$

Then:

  • Mean of tilted $\eta_i$​: $M = \frac{E_1}{E_0}$
  • Variance of tilted $\eta_i$​: $V = \frac{E_2}{E_0} – \left(\frac{E_1}{E_0}\right)^2$.

The integration domain can be taken as e.g. $[M_{-i} – \delta \sqrt{V_{-i}},\ M_{-i} + \delta \sqrt{V_{-i}}]$ with a large $\delta$ (e.g. 10) to cover essentially all the mass.

5.4 Updating the site factor in 1D (for $\eta_i$​)

We want a new normal approximation $N(\eta_i \mid M_i, V_i)$ that, combined with the cavity distribution, gives the tilted mean and variance.

Intuitively:

  • Tilted distribution ∝ cavity × true factor.
  • New site factor should be the “difference” between tilted and cavity in exponential-family terms.

Algebraically, one can express $M_i, V_i$​ so that:

$\underbrace{N(\eta_i \mid M_{-i}, V_{-i})}_{\text{cavity}} \times \underbrace{N(\eta_i \mid M_i, V_i)}_{\text{site}} \propto \underbrace{N(\eta_i \mid M, V)}_{\text{tilted approximation}}$

In natural parameters (for a normal), means and variances combine additively in precision space, so one can solve for $M_i, V_i$ from $M, V$ and $M_{-i}, V_{-i}$.

5.5 Lifting back to θ and updating the global approximation

Once we have the updated scalar site for $\eta_i$​, we convert it back to a multivariate normal factor $g_i(\theta)$:

  • The dependence is only through $X_i \theta$, so the natural parameters of the site in θ-space are:
    • precision contribution: $\Sigma_i^{-1} = \frac{1}{V_i} X_i^\top X_i$​,
    • mean contribution: $\Sigma_i^{-1} \mu_i = \frac{M_i}{V_i} X_i^\top$
    • (dimensions: $X_i$​ is $1 \times k$, so $X_i^\top X_i$​ is $k \times k$, $X_i^\top$ is $k \times 1$).

Then update the global normal approximation by adding natural parameters:

$\Sigma^{-1}_{\text{new}} = \Sigma^{-1}_{-i} + \Sigma_i^{-1}$

$\Sigma^{-1}_{\text{new}} \mu_{\text{new}} = \Sigma^{-1}_{-i} \mu_{-i} + \Sigma_i^{-1} \mu_i$

Cycle over i until convergence.

This is feasible because:

  • For each i, the expensive computation is only a 1D integral in $\eta_i$,
  • The rest is linear algebra with dimension $k$ (number of parameters).

6. Bioassay logistic regression example (2D)

Special case:

  • Simple logistic regression:
    • $y_i \sim \text{Bin}(m_i, \text{logit}^{-1}(\theta_1 + \theta_2 x_i)), \quad i = 1,\dots,4,$,
    • with:
      • $\theta = (\theta_1, \theta_2)$: intercept and slope,
      • a flat prior for $\theta$.

There are only:

  • 2 parameters,
  • 4 likelihood factors.

The exact posterior in this example is:

  • Skewed,
  • Has a long tail for large $\theta_1, \theta_2$​.

A simple Laplace (mode-based normal) approximation:

  • Centers at the posterior mode,
  • Uses the local curvature (Hessian) at that point,
  • Gives a symmetric ellipsoid that does not fully reflect the skew.

EP:

  • Starts with a rough normal,
  • Repeatedly updates site factors for the 4 observations,
  • Adjusts the global $N(\mu, \Sigma)$ to better match the tilted distributions.

Result:

  • The EP approximation ellipse shifts toward the actual mass of the posterior.
  • It better captures location and orientation than the mode-based Laplace approximation.
  • Still normal, so it can’t capture the full skew, but it’s closer to the truth.

Plots in the text:

  • Show the evolution of mean, std devs, and correlation for $\theta_1, \theta_2$​ over iterations.
  • Show the final EP ellipse vs:
    • Laplace ellipse,
    • true posterior contour.

EP’s ellipse lies closer to the true posterior contour.

7. Properties and limitations of EP

  • EP approximates the joint distribution $p(\theta \mid y)$ with a single multivariate normal, not just the marginals.
  • Each step can be seen as moment matching for a tilted distribution.
  • Moment matching corresponds to minimizing a certain local Kullback–Leibler divergence from the tilted distribution to the approximating family.
  • However:
    • The global objective (KL between full posterior and global approximation) is not guaranteed to be minimized overall.
    • EP is not guaranteed to converge in general.

But in many models:

  • Where the factors $f_i(\theta)$ are log-concave,
  • And starting from a reasonable prior-based initialization,

EP performs well in practice.

8. Extensions and improvements

Several extensions are mentioned:

  1. Double-loop EP
    • A theoretically convergent but slower variant.
    • Can be used as a “safety net” if standard EP doesn’t converge.
  2. Damping
    • Do not fully replace old $g_i$ with new $g_i$.
    • Instead, mix old and new (like under-relaxation):
      • $g_i^{\text{new}} \leftarrow (g_i^{\text{old}})^{1-\lambda} (g_i^{\text{EP}})^\lambda,\quad 0<\lambda<1.$
    • Helps avoid oscillations and instability.
  3. Fractional EP (power EP)
    • Raise each factor to a power $\alpha$ to control influence:
      • $f_i(\theta) \to f_i(\theta)^\alpha$
    • This corresponds to minimizing an α-divergence instead of standard KL.
    • Includes various divergences:
      • α → 1: recovers standard EP (KL in one direction),
      • other α: closer to KL in the other direction, or Hellinger distance.
    • Helps when:
      • the approximating family is too rigid,
      • priors are vague,
      • information spreads poorly across the graph.
  4. Improved marginals using ratio tricks
    • Even if $g(\theta)$ is a joint approximation, one can sometimes improve individual marginal posteriors for components $\theta_i$​ using tricks like:
      • $p(\phi \mid y) \approx \frac{p(\gamma,\phi \mid y)}{p_{\text{approx}}(\gamma \mid \phi, y)}$​
    • This can give more accurate 1D marginals than the raw multivariate normal approximation.