Approximate Inference

When the posterior cannot be computed, we replace it by the closest member of a tractable family. This chapter builds the two deterministic ways of doing so, variational inference (minimise \(\KL(q\|p)\) and obtain a lower bound on the evidence) and expectation propagation (match moments factor by factor, a local \(\KL(p\|q)\)), and works every example of the book through to the last line.

TrungTin Nguyen · Following Bishop (2006), Pattern Recognition and Machine Learning, Chapter 10 · Course home · Mathematical toolkit · Python code for the figures

Ch. 1 Introduction Ch. 2 Probability distributions Ch. 3 Linear regression Ch. 4 Linear classification Ch. 5 Neural networks Ch. 6 Kernel methods Ch. 7 Sparse kernel machines Ch. 8 Graphical models Ch. 9 Mixture models and EM Ch. 10 Approximate inference Ch. 11 Sampling methods Ch. 12 Continuous latent variables Ch. 13 Sequential data Ch. 14 Combining models

Throughout, \(\mathbf{X}\) denotes all observed variables and \(\mathbf{Z}\) all unobserved ones: latent variables and, in a fully Bayesian model, the parameters too. The model is the joint density \(p(\mathbf{X},\mathbf{Z})\), which we can evaluate. What we want, the posterior \(p(\mathbf{Z}\given\mathbf{X})\) and the evidence \(p(\mathbf{X})\)\({}=\int p(\mathbf{X},\mathbf{Z})\,\d\mathbf{Z}\), is what we cannot compute: the integral may have no closed form, numerical integration fails in high dimension, and for discrete \(\mathbf{Z}\) the sum may have exponentially many terms (Bishop, pp. 461–462).

10.1 Variational inference

A functional maps a function to a number. The example of the book is the entropy

\[ \mathrm{H}[p]=-\int p(x)\ln p(x)\,\d x . \tag{10.1} \]

Erratum. The first printing omits the minus sign in (10.1).

The functional derivative \(\delta F/\delta q(x)\) says how \(F[q]\) changes when \(q\) is perturbed near \(x\): \(F[q+\epsilon\eta]\)\({}=F[q]+\epsilon\int \frac{\delta F}{\delta q(x)}\eta(x)\,\d x\)\({}+O(\epsilon^2)\) (Appendix D; toolkit: optimising over distributions). Nothing in a variational method is approximate until we restrict the set of functions over which we optimise; the approximation comes entirely from that restriction.

Proof

By the product rule \(p(\mathbf{X},\mathbf{Z})\)\({}=p(\mathbf{Z}\given\mathbf{X})\,p(\mathbf{X})\), so for every \(\mathbf{Z}\) with \(q(\mathbf{Z})\gt0\),

\[ \ln\frac{p(\mathbf{X},\mathbf{Z})}{q(\mathbf{Z})}=\ln\frac{p(\mathbf{Z}\given\mathbf{X})}{q(\mathbf{Z})}+\ln p(\mathbf{X}). \]

Multiply by \(q(\mathbf{Z})\) and integrate over \(\mathbf{Z}\). The last term does not depend on \(\mathbf{Z}\) and \(\int q\,\d\mathbf{Z}=1\), so it contributes \(\ln p(\mathbf{X})\); the other two terms are \(\LB(q)\) and \(-\KL(q\|p)\). This is (10.2).

By Gibbs' inequality (toolkit: KL divergence), \(\KL(q\|p)\ge0\) with equality if and only if \(q=p(\cdot\given\mathbf{X})\) almost everywhere. In one line, by Jensen's inequality for the concave logarithm, \(-\KL(q\|p)\)\({}=\E_q[\ln(p(\mathbf{Z}\given\mathbf{X})/q(\mathbf{Z}))]\)\({}\le\ln\E_q[p(\mathbf{Z}\given\mathbf{X})/q(\mathbf{Z})]\)\({}\)\({}=\ln\int_{q\gt0}p(\mathbf{Z}\given\mathbf{X})\,\d\mathbf{Z}\)\({}\le0\), with equality only when the ratio is constant. Since \(\ln p(\mathbf{X})\) does not depend on \(q\), maximising \(\LB\) and minimising \(\KL(q\|p)\) are the same problem, solved by the posterior, where \(\LB=\ln p(\mathbf{X})\). ∎

The point of (10.2) is that \(\LB(q)\) involves only the joint \(p(\mathbf{X},\mathbf{Z})\), which we can evaluate, while the KL term involves the unknown posterior. So we can make the KL small by making \(\LB\) large without ever computing \(p(\mathbf{Z}\given\mathbf{X})\). The identity is the same as (9.70) in the EM chapter, except that the parameters are now random and live inside \(\mathbf{Z}\). Two further remarks. First, \(\LB\) is a concave functional of \(q\) (the entropy \(-\int q\ln q\) is strictly concave and \(\int q\ln p\) is linear), so over any convex family of densities it has a unique maximiser. Second, the family is restricted purely for tractability: a richer family can only get closer to the posterior, so there is no over-fitting in choosing a flexible \(q\).

The simplest restriction is a parametric family \(q(\mathbf{Z}\given\bm\omega)\), for instance a Gaussian, with \(\LB\) maximised over \(\bm\omega\) by ordinary optimisation. Figure 10.1 does this for the skewed density of Figure 4.14, \(p(z)\propto e^{-z^2/2}\sigma(20z+4)\). The next result explains what the optimal Gaussian is.

Proof

The entropy of \(\N(m,s^2)\) is \(\frac12\ln(2\pi e s^2)\) (toolkit: Gaussian facts), and \(\int q\ln(\tilde p/q)\)\({}=\int q\ln(p/q)+\ln Z\). Write \(z=m+s\varepsilon\) with \(\varepsilon\sim\N(0,1)\) and \(f=\ln\tilde p\). Differentiating under the expectation (allowed by the growth condition), \(\partial_m\E[f(m+s\varepsilon)]\)\({}=\E[f'(z)]\) and \(\partial_s\E[f(m+s\varepsilon)]\)\({}=\E[\varepsilon f'(m+s\varepsilon)]\). Gaussian integration by parts (Stein's identity) gives \(\E[\varepsilon g(\varepsilon)]\)\({}=\E[g'(\varepsilon)]\) for the standard normal: indeed \(\int\varepsilon\phi(\varepsilon)g(\varepsilon)\d\varepsilon\)\({}=-\int\phi'(\varepsilon)g(\varepsilon)\d\varepsilon\)\({}=\int\phi(\varepsilon)g'(\varepsilon)\d\varepsilon\). With \(g(\varepsilon)=f'(m+s\varepsilon)\), \(g'=s f''\). Hence

\[ \frac{\partial\LB}{\partial m}=\E_q[f'(z)],\qquad \frac{\partial\LB}{\partial s}=s\,\E_q[f''(z)]+\frac1s , \]

and setting both to zero gives the two conditions. ∎

For the density of Figure 10.1 I computed the three Gaussians by numerical integration (the numbers below were computed for these notes and are not from the book):

Each criterion wins on its own divergence. The variational Gaussian sits in the bulk of the mass but is narrower than \(p\); the Laplace Gaussian is centred at the mode, near the steep left edge, and so leaks mass into the region where \(p\) is almost zero, which \(\KL(q\|p)\) punishes heavily.

10.1.1 Factorised distributions

The second way to restrict \(q\) leaves the functional form free but imposes independence between groups. Partition \(\mathbf{Z}\) into disjoint groups \(\mathbf{Z}_1,\dots,\mathbf{Z}_M\) and take

\[ q(\mathbf{Z})=\prod_{i=1}^{M}q_i(\mathbf{Z}_i). \tag{10.5} \]

This is the mean-field family of statistical physics. Nothing is assumed about the individual factors \(q_i\); their forms will come out of the optimisation. Write \(q_j\) for \(q_j(\mathbf{Z}_j)\), and \(\E_{i\ne j}[\cdot]\) for the expectation over all factors except the \(j\)th:

\[ \E_{i\ne j}[\ln p(\mathbf{X},\mathbf{Z})]=\int\ln p(\mathbf{X},\mathbf{Z})\prod_{i\ne j}q_i\,\d\mathbf{Z}_i . \tag{10.8} \]

In practice we never compute \(\mathcal{Z}_j\) first. We write \(\ln p(\mathbf{X},\mathbf{Z})\), keep only the terms that involve \(\mathbf{Z}_j\) (everything else is absorbed into "const"), take expectations over the other factors, and recognise the result as the logarithm of a standard density; the normaliser is then read off. Every example below follows this recipe.

Proof

Each update maximises \(\LB\) over \(q_j\) with the others held fixed (Result 10.3), so it cannot decrease \(\LB\). By Result 10.1, \(\LB\le\ln p(\mathbf{X})\). A non-decreasing sequence bounded above converges. ∎

Note on p. 466. The book says convergence holds "because bound is convex with respect to each of the factors". The property used is that \(\LB\) is concave in each factor (Result 10.3), so that each update is a global maximum in that factor. Convergence is of the values of \(\LB\); the limit is in general only a local maximum of \(\LB\) over the product family, because \(\LB\) is not jointly concave in \((q_1,\dots,q_M)\), and different starts can end at different solutions (Figure 10.3).

10.1.2 Properties of factorised approximations

To see what the factorisation costs, approximate a correlated Gaussian \(p(\mathbf{z})\)\({}=\N(\mathbf{z}\given\bm\mu,\bm\Lambda^{-1})\), \(\mathbf{z}=(z_1,z_2)\), by \(q(\mathbf{z})=q_1(z_1)q_2(z_2)\), where

\[ \bm\mu=\begin{pmatrix}\mu_1\\ \mu_2\end{pmatrix},\qquad \bm\Lambda=\begin{pmatrix}\Lambda_{11} & \Lambda_{12}\\ \Lambda_{21} & \Lambda_{22}\end{pmatrix},\qquad \Lambda_{21}=\Lambda_{12}. \tag{10.10} \]

Proof

(a) Put \(e_i=\E[z_i]-\mu_i\). Using \(\E[z_i]=m_i\), (10.13) and (10.15) read \(e_1\)\({}=-(\Lambda_{12}/\Lambda_{11})e_2\) and \(e_2\)\({}=-(\Lambda_{12}/\Lambda_{22})e_1\). Substituting the second into the first, \(e_1=\rho^2e_1\) with \(\rho^2\)\({}=\Lambda_{12}^2/(\Lambda_{11}\Lambda_{22})\). Since \(\det\bm\Lambda\)\({}=\Lambda_{11}\Lambda_{22}\)\({}-\Lambda_{12}^2\)\({}\gt0\), \(\rho^2\lt1\), so \(e_1=0\) and then \(e_2=0\).

(b) For a \(2\times2\) matrix, \(\bm\Sigma\)\({}=\bm\Lambda^{-1}\)\({}=\frac{1}{\det\bm\Lambda}\begin{pmatrix}\Lambda_{22} & -\Lambda_{12}\\ -\Lambda_{12} & \Lambda_{11}\end{pmatrix}\). So \(\Sigma_{11}\)\({}=\Lambda_{22}/\det\bm\Lambda\) and the correlation satisfies \(\rho^2\)\({}=\Sigma_{12}^2/(\Sigma_{11}\Sigma_{22})\)\({}=\Lambda_{12}^2/(\Lambda_{11}\Lambda_{22})\). Then \(\Sigma_{11}(1-\rho^2)\)\({}=\frac{\Lambda_{22}}{\det\bm\Lambda}\cdot\frac{\det\bm\Lambda}{\Lambda_{11}\Lambda_{22}}\)\({}=\Lambda_{11}^{-1}\).

(c) One sweep updates \(e_1'\)\({}=-(\Lambda_{12}/\Lambda_{11})e_2\) and then \(e_2'\)\({}=-(\Lambda_{12}/\Lambda_{22})e_1'\)\({}=\rho^2e_2\). ∎

The same holds in \(D\) dimensions: each mean-field factor of a Gaussian has variance \(1/\Lambda_{ii}\)\({}=\var_p[z_i\given\mathbf{z}_{-i}]\), which is smaller than the marginal variance whenever \(z_i\) is correlated with the rest. The approximation is controlled by the direction of smallest variance of \(p\) and is too compact along the others (Figure 10.2a). With \(\rho=0.9\), for instance, the factorised \(q\) has only \(19\%\) of the true marginal variance, and the cyclic updates shrink the error by only \(0.81\) per sweep.

The reverse divergence. Expectation propagation (Section 10.7) minimises \(\KL(p\|q)\) instead. For a factorised \(q\),

\[ \KL(p\|q)=-\int p(\mathbf{Z})\Big[\sum_{i=1}^{M}\ln q_i(\mathbf{Z}_i)\Big]\d\mathbf{Z}+\const \tag{10.16} \]

where the constant is \(\int p\ln p=-\mathrm{H}[p]\) (the book calls it the entropy; it is minus the entropy).

Proof

In (10.16), \(\int p(\mathbf{Z})\ln q_i(\mathbf{Z}_i)\,\d\mathbf{Z}\)\({}=\int p_i\ln q_i\,\d\mathbf{Z}_i\) (integrate out the other groups first). Add and subtract \(\sum_i\int p_i\ln p_i\):

\[ \KL\Big(p\Big\|\prod_iq_i\Big)=\Big[\int p\ln p-\sum_i\int p_i\ln p_i\Big]+\sum_i\int p_i\ln\frac{p_i}{q_i}. \]

The bracket equals \(\int p\ln\big(p/\prod_ip_i\big)\)\({}=\KL(p\|\prod_ip_i)\), because \(\int p\ln p_i\,\d\mathbf{Z}\)\({}=\int p_i\ln p_i\,\d\mathbf{Z}_i\). Each \(\KL(p_i\|q_i)\ge0\), with equality iff \(q_i=p_i\). The Lagrange-multiplier route asked for in Exercise 10.3 is in the solution of Exercise 10.3. ∎

For the Gaussian, the marginals are \(\N(z_i\given\mu_i,\Sigma_{ii})\) (2.98): again the right means, but now the full marginal variances, so \(q\) puts substantial mass where \(p\) is small, outside the narrow correlated ellipse (Figure 10.2b). The first term \(\KL(p\|\prod_ip_i)\) is the multi-information of \(p\): the irreducible error of any factorised approximation under this divergence.

Zero-forcing and zero-avoiding. In

\[ \KL(q\|p)=-\int q(\mathbf{Z})\ln\left\{\frac{p(\mathbf{Z})}{q(\mathbf{Z})}\right\}\d\mathbf{Z} \tag{10.18} \]

the integrand \(q\ln(q/p)\) is huge wherever \(p\approx0\) but \(q\) is not, so the minimiser keeps \(q\approx0\) wherever \(p\approx0\): it is zero-forcing, it under-estimates the support of \(p\), and on a multimodal \(p\) it settles on one mode. In \(\KL(p\|q)=\int p\ln(p/q)\) the penalty is large wherever \(p\gt0\) but \(q\approx0\), so \(q\) must cover all of the mass of \(p\): it is zero-avoiding. For a single Gaussian \(q\) the minimiser of \(\KL(p\|q)\) matches the overall mean and covariance of \(p\) (Result 10.27 below), which for a bimodal \(p\) puts the peak of \(q\) between the modes (Figure 10.3a), while \(\KL(q\|p)\) has one local minimum per mode (Figures 10.3b, c). In a mixture model the modes are the \(K!\) relabellings (Section 10.2.4) and averaging over them would be useless, which is why variational Bayes, not expectation propagation, is used for mixtures.

Erratum. In the first printing the integrand of (10.20) is not squared.

Proof

Put \(a=(1+\alpha)/2\), so \((1-\alpha)/2=1-a\), \(1-\alpha^2=4a(1-a)\) and \(\mathrm{D}_\alpha\)\({}=\{1-I(a)\}/\{a(1-a)\}\) with \(I(a)=\int q\,(p/q)^a\).

(a) If \(0\lt a\lt1\) (\(|\alpha|\lt1\)), \(t\mapsto t^a\) is concave, so by Jensen \(I(a)\)\({}=\E_q[(p/q)^a]\)\({}\le(\E_q[p/q])^a\)\({}\le1\), and \(a(1-a)\gt0\). If \(a\gt1\) or \(a\lt0\) (\(|\alpha|\gt1\)), \(t^a\) is convex and \(a(1-a)\lt0\). Either the integral is infinite (when \(p\gt0=q\) on a set of positive mass for \(a\gt1\), or \(q\gt0=p\) for \(a\lt0\)), or the supports agree and Jensen gives \(I(a)\ge(\E_q[p/q])^a=1\); in both cases \(\mathrm{D}_\alpha\ge0\). Jensen's inequality is an equality only if \(p/q\) is constant, i.e. \(p=q\).

(b) (The route suggested in Exercise 10.6.) As \(\alpha\to1\), \(\epsilon=1-a\to0\) and \(I\)\({}=\int p\,(q/p)^{\epsilon}\)\({}=\int p\,e^{\epsilon\ln(q/p)}\)\({}=1+\epsilon\int p\ln(q/p)\)\({}+O(\epsilon^2)\), so \(\mathrm{D}_\alpha\)\({}=\{-\epsilon\int p\ln(q/p)\)\({}+O(\epsilon^2)\}/\{(1-\epsilon)\epsilon\}\)\({}\to\int p\ln(p/q)\)\({}=\KL(p\|q)\). As \(\alpha\to-1\), \(a\to0\) and \(I=1+a\int q\ln(p/q)+O(a^2)\), so \(\mathrm{D}_\alpha\)\({}\to-\int q\ln(p/q)\)\({}=\KL(q\|p)\). (Term-by-term integration is justified, for example, when \(p/q\) is bounded above and below on the common support.)

(c) At \(\alpha=0\), \(\mathrm{D}_0=4(1-\int\sqrt{pq})\), while \(\mathrm{D}_{\mathrm H}\)\({}=\int p+\int q-2\int\sqrt{pq}\)\({}=2(1-\int\sqrt{pq})\). And \(\sqrt{\mathrm{D}_{\mathrm H}}\)\({}=\|\sqrt p-\sqrt q\|_{L^2}\) is a norm distance, hence a metric.

(d) For \(a\lt0\) the integrand \(p^aq^{1-a}\) is infinite where \(p=0\lt q\), so any \(q\) with \(\mathrm{D}_\alpha\lt\infty\) vanishes where \(p\) does (and \(a=0\) is the limit \(\KL(q\|p)\), which behaves the same way). For \(a\gt1\) the factor \(q^{1-a}\) is infinite where \(q=0\lt p\), so \(q\) must be positive wherever \(p\) is (and \(a=1\) is \(\KL(p\|q)\)). ∎

Minka (2005) shows that variational message passing (\(\alpha=-1\)), loopy belief propagation and expectation propagation (\(\alpha=1\)), and fractional and power EP (other \(\alpha\)) are all local minimisations of members of this family (Bishop, p. 517).

10.1.3 Example: the univariate Gaussian

Data \(\mathcal{D}=\{x_1,\dots,x_N\}\) are iid \(\N(\mu,\tau^{-1})\) with unknown mean \(\mu\) and precision \(\tau\) (MacKay 2003):

\[ p(\mathcal{D}\given\mu,\tau)=\Big(\frac{\tau}{2\pi}\Big)^{N/2}\exp\Big\{-\frac\tau2\sum_{n=1}^{N}(x_n-\mu)^2\Big\} \tag{10.21} \] \[ p(\mu\given\tau)=\N\big(\mu\given\mu_0,(\lambda_0\tau)^{-1}\big) \tag{10.22} \] \[ p(\tau)=\Gam(\tau\given a_0,b_0). \tag{10.23} \]

This Gaussian–Gamma prior is conjugate, so the exact posterior is again Gaussian–Gamma (Exercise 2.44) and we can compare. The mean-field approximation is

\[ q(\mu,\tau)=q_\mu(\mu)\,q_\tau(\tau). \tag{10.24} \]

The true posterior does not factorise, because the precision of \(\mu\) given \(\tau\) is proportional to \(\tau\).

For this model the fixed point can be found in closed form, and the comparison with the exact posterior is instructive.

Proof

Fixed point. Substitute \(\E[\tau]=a_N/b_N\) into \(b_N=b_0+\frac12(C+1/\E[\tau])\): \(b_N=b_0+\frac C2+\frac{b_N}{2a_N}\), so \(b_N(1-\frac1{2a_N})=b_0+\frac C2\) and \(b_N=\frac{(2b_0+C)\,a_N}{2a_N-1}\). With \(2a_N-1=2a_0+N\) this gives \(b_N\) and \(\E_q[\tau]\)\({}=a_N/b_N\)\({}=(2a_0+N)/(2b_0+C)\).

Exact posterior. Multiplying (10.21)–(10.23) and completing the square in \(\mu\) gives the Gaussian–Gamma posterior with \(a^\ast=a_0+N/2\) and \(b^\ast\)\({}=b_0+\frac12\big[\sum_n(x_n-\bar x)^2+\frac{\lambda_0N}{\lambda_0+N}(\bar x-\mu_0)^2\big]\) (Exercise 2.44). To see that \(2b^\ast=2b_0+C\), write \(\sum_n(x_n-\mu_N)^2\)\({}=\sum_n(x_n-\bar x)^2\)\({}+N(\bar x-\mu_N)^2\) and use \(\bar x-\mu_N\)\({}=\frac{\lambda_0(\bar x-\mu_0)}{\lambda_0+N}\), \(\mu_N-\mu_0\)\({}=\frac{N(\bar x-\mu_0)}{\lambda_0+N}\): then \(N(\bar x-\mu_N)^2+\lambda_0(\mu_N-\mu_0)^2\)\({}=\frac{N\lambda_0^2+\lambda_0N^2}{(\lambda_0+N)^2}(\bar x-\mu_0)^2\)\({}=\frac{\lambda_0N}{\lambda_0+N}(\bar x-\mu_0)^2\).

(i) Hence \(\E_q[\tau]\)\({}=\frac{2a_0+N}{2b^\ast}\)\({}=\frac{a^\ast}{b^\ast}\)\({}=\E[\tau\given\mathcal{D}]\). The mean of \(\mu\) is \(\mu_N\) under both.

(ii) \(\var_q[\tau]\)\({}=a_N/b_N^2\)\({}=\E_q[\tau]^2/a_N\) and \(\var[\tau\given\mathcal{D}]\)\({}=a^\ast/b^{\ast2}\)\({}=\E[\tau\given\mathcal{D}]^2/a^\ast\); the means agree and \(a_N=a^\ast+\frac12\).

(iii) \(\var_q[\mu]\)\({}=\lambda_N^{-1}\)\({}=\frac{b^\ast}{(\lambda_0+N)a^\ast}\). Integrating \(\tau\) out of the exact posterior gives a Student t density with \(2a^\ast\) degrees of freedom, location \(\mu_N\) and precision \((\lambda_0+N)a^\ast/b^\ast\) (toolkit: Student t), whose variance is \(\frac{2a^\ast}{2a^\ast-2}\cdot\frac{b^\ast}{(\lambda_0+N)a^\ast}\)\({}=\frac{b^\ast}{(\lambda_0+N)(a^\ast-1)}\). ∎

I checked these formulas numerically against the iteration and against direct numerical integration of the exact posterior. With the non-informative choice below (\(a^\ast=N/2\)) and \(N=5\), \(q_\mu\) has only \(60\%\) of the posterior variance of \(\mu\); with \(N=100\), \(98\%\).

Non-informative priors. Take \(\mu_0=a_0=b_0=\lambda_0=0\). The priors are improper but the posterior is proper (for \(N\ge2\)). Then \(\mu_N=\bar x\), \(\lambda_N=N\E[\tau]\), \(a_N=\frac{N+1}{2}\), \(b_N\)\({}=\frac12\E_\mu[\sum_n(x_n-\mu)^2]\), and with \(\overline{x^2}=\frac1N\sum_nx_n^2\),

\[ \frac{1}{\E[\tau]}=\frac{b_N}{a_N}=\E\Big[\frac{1}{N+1}\sum_{n=1}^{N}(x_n-\mu)^2\Big]=\frac{N}{N+1}\Big(\overline{x^2}-2\bar x\,\E[\mu]+\E[\mu^2]\Big) \tag{10.31} \] \[ \E[\mu]=\bar x,\qquad \E[\mu^2]=\bar x^2+\frac{1}{N\E[\tau]}. \tag{10.32} \]

10.1.4 Model comparison

Let \(m\) index candidate models with prior probabilities \(p(m)\). Different models can have different latent spaces, so \(q(\mathbf{Z})q(m)\) makes no sense; instead take \(q(\mathbf{Z},m)\)\({}=q(\mathbf{Z}\given m)q(m)\). Then (with sums over a discrete \(\mathbf{Z}\); integrals otherwise)

\[ \ln p(\mathbf{X})=\LB-\sum_m\sum_{\mathbf{Z}}q(\mathbf{Z}\given m)q(m)\ln\left\{\frac{p(\mathbf{Z},m\given\mathbf{X})}{q(\mathbf{Z}\given m)q(m)}\right\} \tag{10.34} \] \[ \LB=\sum_m\sum_{\mathbf{Z}}q(\mathbf{Z}\given m)q(m)\ln\left\{\frac{p(\mathbf{Z},\mathbf{X},m)}{q(\mathbf{Z}\given m)q(m)}\right\}. \tag{10.35} \]

Proof

(10.34) is Result 10.1 applied to the pair \((\mathbf{Z},m)\) as the unobserved variables: \(p(\mathbf{Z},\mathbf{X},m)\)\({}=p(\mathbf{Z},m\given\mathbf{X})p(\mathbf{X})\). For (10.36), use \(p(\mathbf{Z},\mathbf{X},m)\)\({}=p(\mathbf{Z},\mathbf{X}\given m)p(m)\) in (10.35) and sum over \(\mathbf{Z}\) first:

\[ \begin{aligned}\LB&=\sum_mq(m)\Big\{\ln p(m)-\ln q(m)+\LB_m\Big\}=\sum_mq(m)\ln\frac{p(m)e^{\LB_m}}{q(m)}\\ &=-\KL\big(q(m)\,\|\,q^\star(m)\big)+\ln\sum_{m'}p(m')e^{\LB_{m'}},\end{aligned} \]

with \(q^\star(m)\)\({}=p(m)e^{\LB_m}/\sum_{m'}p(m')e^{\LB_{m'}}\). The KL term vanishes iff \(q=q^\star\). (The Lagrange-multiplier version is the solution of Exercise 10.11.) The terms of \(\LB\) involving \(q(\mathbf{Z}\given m)\) are \(q(m)\LB_m\), so for each \(m\) they are maximised by maximising \(\LB_m\), independently of the other models. Finally \(\LB_m\le\ln p(\mathbf{X}\given m)\) gives \(\ln\sum_mp(m)e^{\LB_m}\)\({}\le\ln\sum_mp(m)p(\mathbf{X}\given m)\)\({}=\ln p(\mathbf{X})\). ∎

So the posterior over models is approximated by replacing each log evidence \(\ln p(\mathbf{X}\given m)\) with its lower bound \(\LB_m\). This is how the bounds of Figures 10.7 and 10.9 are used to choose the number of components and the polynomial degree.

Erratum. The first printing writes \(\LB_m\) for \(\LB\) in (10.34), (10.35) and the line between them, and does not define \(\LB_m\) after (10.36). The statements above follow the corrected text.

10.2 Illustration: variational mixture of Gaussians

The Gaussian mixture of Chapter 9 is the model on which variational Bayes is best understood (Attias 1999), and the book asks the reader to work through it in detail. Each \(\mathbf{x}_n\in\R^D\) has a one-of-\(K\) latent vector \(\mathbf{z}_n\); \(\mathbf{X}=\{\mathbf{x}_n\}\), \(\mathbf{Z}=\{\mathbf{z}_n\}\). Precision matrices \(\bm\Lambda_k\) are used instead of covariances because they simplify the algebra.

Both priors are conjugate, which is what will make every factor of the variational posterior a standard distribution. For \(\alpha_0\lt1\) the Dirichlet density is largest near the faces of the simplex (Figure 2.5), so the prior favours solutions with some \(\pi_k\) close to zero; this is the source of the pruning seen in Figure 10.6.

10.2.1 Variational distribution

The joint distribution of all the variables is

\[ p(\mathbf{X},\mathbf{Z},\bm\pi,\bm\mu,\bm\Lambda)=p(\mathbf{X}\given\mathbf{Z},\bm\mu,\bm\Lambda)\,p(\mathbf{Z}\given\bm\pi)\,p(\bm\pi)\,p(\bm\mu\given\bm\Lambda)\,p(\bm\Lambda), \tag{10.41} \]

one factor per node of Figure 10.5, each conditioned on its parents. The only assumption we make is a split between latent variables and parameters:

\[ q(\mathbf{Z},\bm\pi,\bm\mu,\bm\Lambda)=q(\mathbf{Z})\,q(\bm\pi,\bm\mu,\bm\Lambda). \tag{10.42} \]

All further structure (a factor per data point, a separate Dirichlet, a Gaussian–Wishart per component) will be derived. We drop the subscripts on \(q\) and let the arguments say which factor is meant.

Compare the maximum likelihood responsibilities (9.13), written with precisions,

\[ r_{nk}\propto\pi_k|\bm\Lambda_k|^{1/2}\exp\Big\{-\frac12(\mathbf{x}_n-\bm\mu_k)^\T\bm\Lambda_k(\mathbf{x}_n-\bm\mu_k)\Big\}. \tag{10.68} \]

Point values are replaced by geometric means (\(\tilde\pi_k=e^{\E[\ln\pi_k]}\), \(\tilde\Lambda_k\)\({}=e^{\E\ln|\bm\Lambda_k|}\)) and the quadratic form gains the term \(D/\beta_k\), which accounts for the uncertainty in \(\bm\mu_k\). The alternation is a variational EM algorithm: a variational E step computes the expectations and responsibilities, a variational M step re-computes \(q(\bm\pi)\) and \(q(\bm\mu_k,\bm\Lambda_k)\) from the same weighted sums (10.51)–(10.53) as maximum likelihood EM. Every factor of the variational posterior has the same form as the corresponding factor of the joint (10.41): a general consequence of conjugacy (Section 10.4.1).

The posterior mean of the mixing proportions is (Exercise 10.15; the Dirichlet mean \(\alpha_k/\hat\alpha\) with \(\hat\alpha\)\({}=\sum_k(\alpha_0+N_k)\)\({}=K\alpha_0+N\))

\[ \E[\pi_k]=\frac{\alpha_0+N_k}{K\alpha_0+N}. \tag{10.69} \]

Erratum. The first printing has \(\alpha_k\) instead of \(\alpha_0\) in the numerator of (10.69).

Proof

From \(\Gamma(x+1)=x\Gamma(x)\), \(\psi(x+1)=\psi(x)+1/x\), so \(\psi(\alpha_0)\)\({}=\psi(1+\alpha_0)-1/\alpha_0\)\({}=\psi(1)+O(\alpha_0)-1/\alpha_0\) with \(\psi(1)=-\gamma\). For large \(x\), \(\psi(x)=\ln x-\frac1{2x}+O(x^{-2})\), so \(e^{\psi(x)}\)\({}=x\,e^{-1/(2x)+O(x^{-2})}\)\({}=x-\frac12+O(1/x)\). The common factor \(e^{-\psi(\hat\alpha)}\) cancels in the normalisation (10.49). The statements about \(\E[\pi_k]\) are (10.69). ∎

This is the mechanical side of the "automatic trade-off between fitting the data and the complexity of the model" (Section 3.4). The statistical side is visible in the bound of Section 10.2.2, which can be rearranged as an expected log-likelihood minus \(\KL(q(\bm\pi)\|p(\bm\pi))\) and minus \(\sum_k\KL(q(\bm\mu_k,\bm\Lambda_k)\|p(\bm\mu_k,\bm\Lambda_k))\): a component pays a KL price for moving away from its prior, and a component that merely duplicates another improves the fit too little to pay it. Maximum likelihood has no such price, so EM with \(K=6\) keeps all six components.

Three practical advantages over maximum likelihood come at almost no extra cost (the dominant work, the responsibilities and the weighted covariances, is the same as in EM). The singularities of Section 9.2.1, where a component collapses onto one data point, are absent: \(\mathbf{W}_k^{-1}\succeq\mathbf{W}_0^{-1}\) keeps every precision bounded (they are absent from MAP estimation too, Exercise 10.24). Large \(K\) does not over-fit (Figure 10.6). And the bound can be used to choose \(K\) (Section 10.2.4).

10.2.2 Variational lower bound

Monitoring \(\LB\) is the best test of both the algebra and the code: it must not decrease at any step, and finite differences can check that each update maximises it (Svensén and Bishop 2005). For the mixture,

\[ \begin{aligned}\LB&=\sum_{\mathbf{Z}}\iiint q(\mathbf{Z},\bm\pi,\bm\mu,\bm\Lambda)\ln\left\{\frac{p(\mathbf{X},\mathbf{Z},\bm\pi,\bm\mu,\bm\Lambda)}{q(\mathbf{Z},\bm\pi,\bm\mu,\bm\Lambda)}\right\}\d\bm\pi\,\d\bm\mu\,\d\bm\Lambda\\ &=\E[\ln p(\mathbf{X}\given\mathbf{Z},\bm\mu,\bm\Lambda)]+\E[\ln p(\mathbf{Z}\given\bm\pi)]+\E[\ln p(\bm\pi)]+\E[\ln p(\bm\mu,\bm\Lambda)]\\ &\quad-\E[\ln q(\mathbf{Z})]-\E[\ln q(\bm\pi)]-\E[\ln q(\bm\mu,\bm\Lambda)]\end{aligned} \tag{10.70} \]

where each expectation is over all the variables in its argument, under the current \(q\).

Proof

(a) By (10.6), \(\LB\)\({}=\E_{q_j}\{\E_{i\ne j}[\ln p]\}+\mathrm{H}[q_j]+\const\). With \(\bm\mu(\bm\eta)\)\({}=\E_{\bm\eta}[\mathbf{u}]\)\({}=\nabla A(\bm\eta)\), the first term is \(\tilde{\bm\eta}^\T\bm\mu(\bm\eta)+\E[\ln h]\) and the entropy is \(-\E[\ln h]-\bm\eta^\T\bm\mu(\bm\eta)+A(\bm\eta)\), so \(\LB\)\({}=(\tilde{\bm\eta}-\bm\eta)^\T\bm\mu(\bm\eta)\)\({}+A(\bm\eta)+\const\). Differentiating, \(\nabla\LB\)\({}=\nabla^2A\,(\tilde{\bm\eta}-\bm\eta)\)\({}-\bm\mu+\nabla A\)\({}=\nabla^2A\,(\tilde{\bm\eta}-\bm\eta)\), and \(\nabla^2A=\cov[\mathbf{u}]\) is positive definite for a minimal family (toolkit: exponential families).

(b) The explicit derivatives, from Result 10.12:

Responsibilities. (10.71) and (10.72) together equal \(\sum_{n,k}r_{nk}\ln\rho_{nk}\) with \(\ln\rho_{nk}\) of (10.46) (expand \(N_k,\bar{\mathbf{x}}_k,\mathbf{S}_k\) back into sums over \(n\)), and (10.75) contributes \(-\sum_{n,k}r_{nk}\ln r_{nk}\). With a multiplier \(\lambda_n\) for \(\sum_kr_{nk}=1\): \(\ln\rho_{nk}-\ln r_{nk}-1+\lambda_n\)\({}=0\), so \(r_{nk}\propto\rho_{nk}\), which is (10.49); the entropy term makes this a maximum.

Dirichlet. The terms with \(\bm\alpha\) are \(\sum_k(\alpha_0+N_k-\alpha_k)\ln\tilde\pi_k\)\({}-\ln C(\bm\alpha)+\const\), with \(\ln\tilde\pi_k\)\({}=\psi(\alpha_k)-\psi(\hat\alpha)\). Since \(\ln C(\bm\alpha)\)\({}=\ln\Gamma(\hat\alpha)\)\({}-\sum_k\ln\Gamma(\alpha_k)\), \(\partial\ln C/\partial\alpha_j\)\({}=\psi(\hat\alpha)-\psi(\alpha_j)\)\({}=-\ln\tilde\pi_j\), which cancels the term \(-\ln\tilde\pi_j\) from differentiating the coefficient. What is left is \(\sum_k(\alpha_0+N_k-\alpha_k)\{\psi'(\alpha_k)\delta_{kj}\)\({}-\psi'(\hat\alpha)\}\)\({}=0\) for all \(j\). The matrix \(\{\psi'(\alpha_k)\delta_{kj}\)\({}-\psi'(\hat\alpha)\}\) is the covariance of \((\ln\pi_1,\dots,\ln\pi_K)\) under the Dirichlet, positive definite (no non-zero combination \(\sum_kv_k\ln\pi_k\) is constant on the simplex), so \(\alpha_k=\alpha_0+N_k\).

Means. The terms with \(\mathbf{m}_k\) are \(-\frac{\nu_k}2\{N_k(\bar{\mathbf{x}}_k-\mathbf{m}_k)^\T\mathbf{W}_k(\bar{\mathbf{x}}_k-\mathbf{m}_k)\)\({}+\beta_0(\mathbf{m}_k-\mathbf{m}_0)^\T\mathbf{W}_k(\mathbf{m}_k-\mathbf{m}_0)\}\); the gradient \(\nu_k\mathbf{W}_k\{N_k(\bar{\mathbf{x}}_k-\mathbf{m}_k)\)\({}-\beta_0(\mathbf{m}_k-\mathbf{m}_0)\}\) vanishes iff (10.61) holds.

\(\beta_k\). The terms are \(-\frac{D(N_k+\beta_0)}{2\beta_k}\)\({}-\frac D2\ln\beta_k\); the derivative \(\frac{D(N_k+\beta_0)}{2\beta_k^2}\)\({}-\frac{D}{2\beta_k}\) vanishes iff \(\beta_k=\beta_0+N_k\), and the second derivative there is \(-D/(2\beta_k^2)\lt0\).

Wishart parameters. Collecting the terms with \(\mathbf{W}=\mathbf{W}_k\), \(\nu=\nu_k\) (the coefficient of \(\ln\tilde\Lambda_k\) was found in Gap 10.11):

\[ f(\mathbf{W},\nu)=\frac{\nu_0+N_k-\nu}{2}\ln\tilde\Lambda(\mathbf{W},\nu)-\frac\nu2\Tr(\mathbf{W}\mathbf{M}_k)-\ln B(\mathbf{W},\nu)+\frac{\nu D}{2}, \]

with \(\mathbf{M}_k\)\({}=\mathbf{W}_0^{-1}\)\({}+N_k\mathbf{S}_k+N_k(\bar{\mathbf{x}}_k-\mathbf{m}_k)(\cdot)^\T\)\({}+\beta_0(\mathbf{m}_k-\mathbf{m}_0)(\cdot)^\T\), which equals the right side of (10.62) at (10.61) (Gap 10.8). Using \(\partial\ln\tilde\Lambda/\partial\mathbf{W}\)\({}=\mathbf{W}^{-1}\), \(\partial\ln B/\partial\mathbf{W}\)\({}=-\frac\nu2\mathbf{W}^{-1}\), \(\partial\ln B/\partial\nu\)\({}=-\frac12\ln\tilde\Lambda\) and \(\partial\ln\tilde\Lambda/\partial\nu\)\({}=\frac12\sum_i\psi'(\frac{\nu+1-i}{2})\):

\[ \frac{\partial f}{\partial\mathbf{W}}=\frac{\nu_0+N_k}{2}\mathbf{W}^{-1}-\frac\nu2\mathbf{M}_k,\qquad \frac{\partial f}{\partial\nu}=\frac{\nu_0+N_k-\nu}{2}\frac{\partial\ln\tilde\Lambda}{\partial\nu}-\frac12\Tr(\mathbf{W}\mathbf{M}_k)+\frac D2. \]

The first vanishes iff \(\mathbf{W}\)\({}=\frac{\nu_0+N_k}{\nu}\mathbf{M}_k^{-1}\); substituting into the second gives \(\frac{\nu_0+N_k-\nu}{2}\big\{\frac{\partial\ln\tilde\Lambda}{\partial\nu}\)\({}-\frac D\nu\big\}\). Since \(\psi'(x)\)\({}=\sum_{j\ge0}(x+j)^{-2}\)\({}\gt\int_0^\infty(x+t)^{-2}\d t\)\({}=1/x\), the braces are positive, so the only solution is \(\nu=\nu_0+N_k\), and then \(\mathbf{W}^{-1}=\mathbf{M}_k\). ∎

So the bound offers a second route to the update equations: fix the functional forms that conjugacy guarantees (discrete, Dirichlet, Gaussian–Wishart), write \(\LB\) as a function of their parameters and maximise. In code, the finite-difference version of this check catches most errors.

10.2.3 Predictive density

For a new point \(\hat{\mathbf{x}}\) with latent \(\hat{\mathbf{z}}\),

\[ p(\hat{\mathbf{x}}\given\mathbf{X})=\sum_{\hat{\mathbf{z}}}\iiint p(\hat{\mathbf{x}}\given\hat{\mathbf{z}},\bm\mu,\bm\Lambda)\,p(\hat{\mathbf{z}}\given\bm\pi)\,p(\bm\pi,\bm\mu,\bm\Lambda\given\mathbf{X})\,\d\bm\pi\,\d\bm\mu\,\d\bm\Lambda. \tag{10.78} \]

Summing over the \(K\) values of \(\hat{\mathbf{z}}\) (the summand is \(\prod_k[\pi_k\N(\hat{\mathbf{x}}\given\bm\mu_k,\bm\Lambda_k^{-1})]^{\hat z_k}\)),

\[ p(\hat{\mathbf{x}}\given\mathbf{X})=\sum_{k=1}^{K}\iiint\pi_k\,\N\big(\hat{\mathbf{x}}\given\bm\mu_k,\bm\Lambda_k^{-1}\big)\,p(\bm\pi,\bm\mu,\bm\Lambda\given\mathbf{X})\,\d\bm\pi\,\d\bm\mu\,\d\bm\Lambda \tag{10.79} \]

and replacing the unknown posterior by \(q(\bm\pi)\prod_jq(\bm\mu_j,\bm\Lambda_j)\), in each term the factors \(j\ne k\) integrate to one:

\[ p(\hat{\mathbf{x}}\given\mathbf{X})\simeq\sum_{k=1}^{K}\iiint\pi_k\,\N\big(\hat{\mathbf{x}}\given\bm\mu_k,\bm\Lambda_k^{-1}\big)\,q(\bm\pi)\,q(\bm\mu_k,\bm\Lambda_k)\,\d\bm\pi\,\d\bm\mu_k\,\d\bm\Lambda_k. \tag{10.80} \]

Erratum. The first printing has "=" instead of "\(\simeq\)" in (10.80) and (10.81).

10.2.4 Determining the number of components

For generic parameter values, relabelling the components (a permutation \(\sigma\) of \(\{1,\dots,K\}\) applied to \((\pi_k,\bm\mu_k,\bm\Lambda_k)\) and to the coordinates of every \(\mathbf{z}_n\)) leaves \(p(\mathbf{X})\) and the prior unchanged; there are \(K!\) such permutations, all giving different parameter vectors (Exercise 10.21). So every mode of the posterior comes with \(K!-1\) equivalent copies. Maximum likelihood does not care, since EM finds one of them. A variational \(q\) minimising \(\KL(q\|p)\) also sits on one copy (Figure 10.3, not Figure 10.2 as in the first printing), which is harmless for prediction with a fixed \(K\), but the bound then misses most of the posterior mass, and the amount it misses depends on \(K\).

Proof

Write \(\bm\theta\) for all unobserved variables (\(\mathbf{Z}\) and the parameters). \(\LB(q_{\rm mix})\)\({}=\E_{q_{\rm mix}}[\ln p(\mathbf{X},\bm\theta)]\)\({}-\E_{q_{\rm mix}}[\ln q_{\rm mix}]\). Energy: \(p(\mathbf{X},\bm\theta)\) is invariant under relabelling (the likelihood is, and the priors (10.39), (10.40) are symmetric in \(k\)), so \(\E_{q\circ\sigma}[\ln p]\)\({}=\E_q[\ln p]\) for every \(\sigma\), and the mixture has the same energy as \(q\). Entropy: on the region where copy \(\sigma\) lives, the other copies are negligible, so \(q_{\rm mix}\)\({}\approx\frac1{K!}q\circ\sigma\) and \(\ln q_{\rm mix}\)\({}\approx\ln(q\circ\sigma)-\ln K!\). Hence \(\E_{q_{\rm mix}}[\ln q_{\rm mix}]\)\({}\approx\frac1{K!}\sum_\sigma\{\E_{q\circ\sigma}[\ln q\circ\sigma]-\ln K!\}\)\({}=\E_q[\ln q]-\ln K!\). Subtracting gives the claim. The entropy of the mixture cannot exceed that of \(q\) by more than \(\ln K!\) (the entropy of the mixture label), with equality in the non-overlapping limit, so the correction is an upper limit on the gain. ∎

Figure 10.7 shows the corrected bound against \(K\) for the Old Faithful data: a clear peak at \(K=2\), whereas the maximised likelihood keeps increasing with \(K\) and cannot choose the model (Section 3.4). Two cautions. First (Exercise 10.25), the factorised \(q\) under-estimates the posterior spread (Section 10.1.2), so \(\LB\) under-estimates \(\ln p(\mathbf{X}\given K)\), and the shortfall tends to grow with the number of parameters: this biases the comparison towards fewer components. Second, the \(\ln K!\) correction assumes that the \(K!\) copies do not overlap; when some components coincide (for example several switched-off components all sitting at the prior) there are fewer distinct copies, and adding \(\ln K!\) over-states the bound for large \(K\).

Point estimates of \(\bm\pi\). Instead of comparing several \(K\), one can drop the prior on \(\bm\pi\), treat \(\bm\pi\) as a parameter and maximise \(\LB\) over it (Corduneanu and Bishop 2001). Only (10.72), now \(\sum_{n,k}r_{nk}\ln\pi_k\), depends on \(\bm\pi\); maximising with a multiplier for \(\sum_k\pi_k=1\), exactly as for (9.22), gives (Exercise 10.23)

\[ \pi_k=\frac1N\sum_{n=1}^{N}r_{nk}, \tag{10.83} \]

interleaved with the variational updates of the other factors. Components that explain too little have \(\pi_k\) driven to zero and drop out: a single run from a generous \(K\) prunes itself, by the same automatic relevance determination mechanism as in the relevance vector machine (Section 7.2.2).

10.2.5 Induced factorisations

We assumed only (10.42), yet \(q^\star\) factorised further: over \(n\) in \(q(\mathbf{Z})\) and into \(q(\bm\pi)\prod_kq(\bm\mu_k,\bm\Lambda_k)\). (Not over \(k\) within \(\mathbf{z}_n\): its components are tied by \(\sum_kz_{nk}=1\).) These induced factorisations come from the conditional independences of the true model, and an implementation should exploit them, for example by never storing a full precision matrix for a Gaussian that is known to factorise. Partition the latent variables into \(\mathbf{A},\mathbf{B},\mathbf{C}\) and assume

\[ q(\mathbf{A},\mathbf{B},\mathbf{C})=q(\mathbf{A},\mathbf{B})\,q(\mathbf{C}). \tag{10.84} \]

Proof

(10.85) is (10.9) followed by the product rule \(\ln p(\mathbf{X},\mathbf{A},\mathbf{B},\mathbf{C})\)\({}=\ln p(\mathbf{A},\mathbf{B}\given\mathbf{X},\mathbf{C})\)\({}+\ln p(\mathbf{X},\mathbf{C})\), whose second term does not involve \(\mathbf{A},\mathbf{B}\). If (10.86) holds, \(\ln p(\mathbf{A},\mathbf{B}\given\mathbf{X},\mathbf{C})\)\({}=\ln p(\mathbf{A}\given\mathbf{X},\mathbf{C})\)\({}+\ln p(\mathbf{B}\given\mathbf{X},\mathbf{C})\) for every \(\mathbf{C}\); averaging over \(\mathbf{C}\) keeps the sum of a function of \(\mathbf{A}\) and a function of \(\mathbf{B}\), whose exponential is a product. Conversely, take \(q(\mathbf{C})\) concentrated at any value \(\mathbf{C}_0\): factorisation of \(q^\star\) then means \(p(\mathbf{A},\mathbf{B}\given\mathbf{X},\mathbf{C}_0)\) factorises, which for all \(\mathbf{C}_0\) is (10.86). ∎

The mixture. With \(\mathbf{C}=\mathbf{Z}\), every path from \(\bm\pi\) to \(\bm\mu\) or \(\bm\Lambda\) in Figure 10.5 passes through some \(\mathbf{z}_n\) as \(\bm\pi\)\({}\to\mathbf{z}_n\)\({}\to\mathbf{x}_n\), head-to-tail at the conditioned node \(\mathbf{z}_n\), hence blocked: \(\bm\pi\perp\!\!\!\perp(\bm\mu,\bm\Lambda)\given\mathbf{X},\mathbf{Z}\), and \(q(\bm\pi)\) separates. With \(\mathbf{C}\) the parameters, \(\mathbf{z}_n\) and \(\mathbf{z}_m\) are connected only through conditioned tail-to-tail nodes (\(\bm\pi\), \(\bm\mu\), \(\bm\Lambda\)), so \(q(\mathbf{Z})\)\({}=\prod_nq(\mathbf{z}_n)\) (also seen directly in Gap 10.6). The factorisation over \(k\) of \(q(\bm\mu,\bm\Lambda)\) is a subtler case: given \(\mathbf{X}\) and \(\mathbf{Z}\) the components are independent, because

\[ p(\bm\mu,\bm\Lambda\given\mathbf{X},\mathbf{Z})\propto\prod_k\Big[p(\bm\mu_k,\bm\Lambda_k)\prod_n\N(\mathbf{x}_n\given\bm\mu_k,\bm\Lambda_k^{-1})^{z_{nk}}\Big], \]

but the graph cannot show it, since \(\mathbf{x}_n\) has every \(\bm\mu_k\) as a parent and conditioning on it couples them. The independence holds for each value of \(\mathbf{Z}\) separately ("context-specific" independence). d-separation is a sufficient test, not a necessary one.

10.3 Variational linear regression

In the evidence framework of Section 3.5 the hyperparameters of Bayesian linear regression were fixed at the values maximising the marginal likelihood. A fully Bayesian treatment integrates over them; exact integration is intractable, but a variational approximation is easy. Keep the noise precision \(\beta\) fixed at its true value (Exercise 10.26 adds a Gamma prior on \(\beta\)) and put a Gamma prior on the weight precision \(\alpha\). With \(\bm\phi_n\)\({}=\bm\phi(\mathbf{x}_n)\in\R^M\) and design matrix \(\bm\Phi\) (rows \(\bm\phi_n^\T\)),

\[ p(\mathbf{t}\given\mathbf{w})=\prod_{n=1}^{N}\N(t_n\given\mathbf{w}^\T\bm\phi_n,\beta^{-1}) \tag{10.87} \] \[ p(\mathbf{w}\given\alpha)=\N(\mathbf{w}\given\mathbf{0},\alpha^{-1}\mathbf{I}) \tag{10.88} \] \[ p(\alpha)=\Gam(\alpha\given a_0,b_0) \tag{10.89} \] \[ p(\mathbf{t},\mathbf{w},\alpha)=p(\mathbf{t}\given\mathbf{w})\,p(\mathbf{w}\given\alpha)\,p(\alpha). \tag{10.90} \]

The graph (Figure 10.8) is the chain \(\alpha\to\mathbf{w}\to t_n\), with the inputs \(\bm\phi_n\) and \(\beta\) as fixed quantities. This example is simple, but it sets up the variational logistic regression of Section 10.6.

10.3.1 Variational distribution

Take \(q(\mathbf{w},\alpha)\)\({}=q(\mathbf{w})\,q(\alpha)\) (10.91).

Relation to the evidence framework. With \(a_0=b_0=0\) (an infinitely broad prior),

\[ \E[\alpha]=\frac{a_N}{b_N}=\frac{M/2}{\E[\mathbf{w}^\T\mathbf{w}]/2}=\frac{M}{\mathbf{m}_N^\T\mathbf{m}_N+\Tr(\mathbf{S}_N)}. \tag{10.104} \]

This is the EM re-estimate (9.63) of \(\alpha\) for the evidence (Section 9.3.4), with the point estimate replaced by \(\E[\alpha]\). Because \(q(\mathbf{w})\) depends on \(q(\alpha)\) only through \(\E[\alpha]\), the variational iteration and the EM iteration are the same map, so they converge to the same \(\mathbf{m}_N,\mathbf{S}_N\), with \(\E[\alpha]\) equal to the evidence maximiser. The variational approach adds a whole distribution for \(\alpha\), \(\Gam(M/2,b_N)\), whose coefficient of variation \(\sqrt{2/M}\) shows how poorly \(\alpha\) is determined when there are few weights.

10.3.2 Predictive distribution

\[ \begin{aligned}p(t\given\mathbf{x},\mathbf{t})&=\int p(t\given\mathbf{x},\mathbf{w})p(\mathbf{w}\given\mathbf{t})\,\d\mathbf{w}\\ &\simeq\int\N(t\given\mathbf{w}^\T\bm\phi(\mathbf{x}),\beta^{-1})\N(\mathbf{w}\given\mathbf{m}_N,\mathbf{S}_N)\,\d\mathbf{w}\\ &=\N\big(t\given\mathbf{m}_N^\T\bm\phi(\mathbf{x}),\sigma^2(\mathbf{x})\big)\end{aligned} \tag{10.105} \] \[ \sigma^2(\mathbf{x})=\frac1\beta+\bm\phi(\mathbf{x})^\T\mathbf{S}_N\bm\phi(\mathbf{x}). \tag{10.106} \]

10.3.3 Lower bound

\[ \begin{aligned}\LB(q)&=\E[\ln p(\mathbf{w},\alpha,\mathbf{t})]-\E[\ln q(\mathbf{w},\alpha)]\\ &=\E_{\mathbf{w}}[\ln p(\mathbf{t}\given\mathbf{w})]+\E_{\mathbf{w},\alpha}[\ln p(\mathbf{w}\given\alpha)]+\E_\alpha[\ln p(\alpha)]-\E_{\mathbf{w}}[\ln q(\mathbf{w})]-\E_\alpha[\ln q(\alpha)].\end{aligned} \tag{10.107} \]

Erratum. The first printing has \(\ln\Gamma(a_N)\) as the last term of (10.110); it should be \(\ln\Gamma(a_0)\) (the normaliser of the prior). The last term of (10.107) is written "\(-\E_\alpha[\ln q(\mathbf{w})]_{\mathbf{w}}\)" in the book; it is \(-\E_{\mathbf{w}}[\ln q(\mathbf{w})]\).

Figure 10.9 plots \(\LB\) against the polynomial order for ten points from a cubic. By Section 10.1.4, with equal prior probabilities for the models, \(q(M)\propto e^{\LB_M}\): the bound acts as an approximate log evidence \(\ln p(\mathbf{t}\given M)\), and it peaks at the true order, whereas the training error keeps decreasing with the order. The figure uses \(a_0=b_0=0\), the prior \(p(\alpha)\propto1/\alpha\) (uniform in \(\ln\alpha\)). Then the term \(-\ln\Gamma(a_0)\) of (10.110) is infinite, but it is the same for every model, so differences of \(\LB\) between models, which are all that model comparison uses, are finite; in practice it is dropped (or \(a_0=b_0\) is set to a tiny positive number).

Erratum. On p. 489, line −7, "\(p(\mathbf{t}\given M)\)" should be "\(\ln p(\mathbf{t}\given M)\)".

10.4 Exponential family distributions

In many models the complete-data likelihood \(p(\mathbf{X},\mathbf{Z}\given\bm\eta)\) is in the exponential family even though the marginal likelihood is not (a Gaussian mixture is not an exponential family, but the joint of \(\mathbf{x}_n\) and its label \(\mathbf{z}_n\) is). Distinguish latent variables \(\mathbf{Z}\), which are extensive (one per data point), from parameters \(\bm\eta\), which are intensive (fixed in number). For iid data,

\[ p(\mathbf{X},\mathbf{Z}\given\bm\eta)=\prod_{n=1}^{N}h(\mathbf{x}_n,\mathbf{z}_n)\,g(\bm\eta)\exp\big\{\bm\eta^\T\mathbf{u}(\mathbf{x}_n,\mathbf{z}_n)\big\} \tag{10.113} \]

with the conjugate prior, interpretable as \(\nu_0\) pseudo-observations whose sufficient statistic has average \(\bm\chi_0\),

\[ p(\bm\eta\given\nu_0,\bm\chi_0)=f(\nu_0,\bm\chi_0)\,g(\bm\eta)^{\nu_0}\exp\big\{\nu_0\bm\eta^\T\bm\chi_0\big\}. \tag{10.114} \]

Erratum. The first printing has \(p(\bm\eta\given\nu_0,\mathbf{v}_0)\) on the left of (10.114), and omits the factors \(\nu_0\) and \(\nu_N\) in front of \(\bm\eta^\T\bm\chi_0\) in (10.118) and \(\bm\eta^\T\bm\chi_N\) in (10.119), and on the left of (10.121). The corrected forms are used below.

The two updates are coupled and are iterated: one stage computes the expected sufficient statistics \(\E[\mathbf{u}(\mathbf{x}_n,\mathbf{z}_n)]\) under \(q(\mathbf{z}_n)\) and with them \(q(\bm\eta)\); the other computes the expected natural parameters \(\E[\bm\eta]\) under \(q(\bm\eta)\) and with them \(q(\mathbf{z}_n)\). The book calls the first the variational E step and the second the variational M step. Exercise 10.28 recovers the Gaussian-mixture results of Section 10.2 from this general form.

10.4.1 Variational message passing

For a directed graph with nodes \(\mathbf{x}_i\) (latent or observed) and parent sets \(\mathrm{pa}_i\),

\[ p(\mathbf{x})=\prod_ip(\mathbf{x}_i\given\mathrm{pa}_i) \tag{10.122} \]

and a variational distribution that factorises over the unobserved nodes (observed nodes get no factor),

\[ q(\mathbf{x})=\prod_iq_i(\mathbf{x}_i), \tag{10.123} \]

the general result (10.9) reads

\[ \ln q_j^\star(\mathbf{x}_j)=\E_{i\ne j}\Big[\sum_i\ln p(\mathbf{x}_i\given\mathrm{pa}_i)\Big]+\const. \tag{10.124} \]

Erratum. On p. 492, line 7 after (10.124), \(q^\star(\mathbf{x}_j)\) should be \(q_j^\star(\mathbf{x}_j)\).

Proof

(a) In (10.124), every term \(\ln p(\mathbf{x}_i\given\mathrm{pa}_i)\) that does not contain \(\mathbf{x}_j\) is a constant after averaging over the other factors. The terms containing \(\mathbf{x}_j\) are \(i=j\) and the children \(i=c\) (for which \(\mathbf{x}_j\in\mathrm{pa}_c\)). Their expectations involve the parents of \(j\), the children and the co-parents, which together form the Markov blanket (Figure 8.26).

(b) Substitute the two forms: \(\ln q_j^\star(\mathbf{x}_j)\)\({}=\E[\bm\eta_j(\mathrm{pa}_j)]^\T\mathbf{u}_j(\mathbf{x}_j)\)\({}+f_j(\mathbf{x}_j)\)\({}+\sum_c\E[\bm\phi_{cj}]^\T\mathbf{u}_j(\mathbf{x}_j)+\const\), which is the stated exponential family. In a conjugate-exponential model \(\bm\eta_j(\mathrm{pa}_j)\) and \(\bm\phi_{cj}\) are multilinear in the sufficient statistics of the variables they depend on, and \(q\) is factorised, so each expectation is obtained by substituting the \(\langle\mathbf{u}_k\rangle\). ∎

Example. In Section 10.1.3, \(\mathbf{u}_\mu(\mu)=(\mu,\mu^2)\). Each data node sends \(\E[\tau](x_n,-\frac12)\) (from \(\ln\N(x_n\given\mu,\tau^{-1})\)\({}=\tau x_n\mu-\frac\tau2\mu^2+\cdots\)), and the prior node sends \(\E[\tau](\lambda_0\mu_0,-\frac{\lambda_0}2)\). Their sum \(\E[\tau](\lambda_0\mu_0+N\bar x,\,-\frac{\lambda_0+N}{2})\) is the natural parameter of \(\N(\mu_N,\lambda_N^{-1})\), which is (10.26)–(10.27). A node can be updated once it has heard from all its parents and children, the children having first heard from their co-parents; the bound is also a sum of local terms, so general-purpose software can run variational inference on a model described only by its graph (Bishop et al. 2003; Winn and Bishop 2005).

10.5 Local variational methods

Sections 10.1–10.4 approximate the whole posterior at once (a "global" method). A local method instead bounds one awkward factor of the model, for example a sigmoid likelihood term, by a simpler function with its own adjustable parameter, and repeats this until the model becomes tractable. Convexity is the tool: a convex function lies above each of its tangent lines (and a concave one below), so tangents give bounds. (Everything below holds for concave functions with "max" and "min", and lower and upper bounds, interchanged.)

Erratum. In the first printing the variational parameter of this section is called \(\lambda\) from Figure 10.10 up to the line before (10.141); it should be \(\eta\), which is the notation used here. The symbol \(\lambda\) is then reserved for \(\lambda(\xi)\) of (10.141), consistent with Jaakkola and Jordan (2000).

A first example. \(f(x)=e^{-x}\) is convex. Its tangent at \(x=\xi\) (first-order Taylor expansion)

\[ y(x)=f(\xi)+f'(\xi)(x-\xi) \tag{10.125} \]

satisfies \(y(x)\le f(x)\), with equality at \(x=\xi\). For this \(f\),

\[ y(x)=e^{-\xi}-e^{-\xi}(x-\xi). \tag{10.126} \]

Re-parameterise by the slope \(\eta=-e^{-\xi}\lt0\):

\[ y(x,\eta)=\eta x-\eta+\eta\ln(-\eta) \tag{10.127} \] \[ f(x)=\max_\eta\{\eta x-\eta+\eta\ln(-\eta)\}. \tag{10.128} \]

If \(f\) is neither convex nor concave, look for an invertible transformation of the function or of its argument that makes it so, bound there, and transform back. The key example is the logistic sigmoid

\[ \sigma(x)=\frac1{1+e^{-x}}, \tag{10.134} \]

which is neither: \(\sigma''\)\({}=\sigma(1-\sigma)(1-2\sigma)\) changes sign at \(x=0\).

Proof

\(\frac{\d}{\d x}\ln\sigma(x)\)\({}=1-\sigma(x)\) (since \(\sigma'=\sigma(1-\sigma)\)) and \(\frac{\d^2}{\d x^2}\ln\sigma(x)\)\({}=-\sigma(x)(1-\sigma(x))\)\({}\lt0\): concave. In (10.133), the derivative of \(\eta x-\ln\sigma(x)\) is \(\eta-1+\sigma(x)\), which vanishes where \(\sigma(x)=1-\eta\), i.e. \(x^\ast=\ln\frac{1-\eta}{\eta}\) (a minimum, since the second derivative \(\sigma(1-\sigma)\) is positive). There \(\ln\sigma(x^\ast)=\ln(1-\eta)\), so

\[ g(\eta)=\eta\ln\frac{1-\eta}{\eta}-\ln(1-\eta)=-\eta\ln\eta-(1-\eta)\ln(1-\eta). \]

(10.136) is (10.132) and (10.137) follows because \(\exp\) is increasing. Taylor route (Exercise 10.30): a concave function lies below its tangent, so \(\ln\sigma(x)\)\({}\le\ln\sigma(\xi)\)\({}+(1-\sigma(\xi))(x-\xi)\); with \(\eta=1-\sigma(\xi)\) the intercept is \(\ln\sigma(\xi)-\eta\xi\)\({}=\ln(1-\eta)-\eta\ln\frac{1-\eta}{\eta}\)\({}=-g(\eta)\), the same bound. ∎

A lower bound of Gaussian form is more useful, because it combines with Gaussian priors. Following Jaakkola and Jordan (2000), transform both the function and its argument:

\[ \ln\sigma(x)=-\ln(1+e^{-x})=-\ln\big\{e^{-x/2}(e^{x/2}+e^{-x/2})\big\}=\frac x2-\ln\big(e^{x/2}+e^{-x/2}\big). \tag{10.138} \]

Proof

Concavity in \(x\). \(f(x)=-\ln(2\cosh\frac x2)\), \(f'(x)=-\frac12\tanh\frac x2\), \(f''(x)\)\({}=-\frac14\operatorname{sech}^2\frac x2\)\({}\lt0\).

Convexity in \(t=x^2\). \(f\) is even, so it is a function of \(t\ge0\): with \(u=\sqrt t\), \(f=-\ln(2\cosh\frac u2)\). By the chain rule (\(\d u/\d t=1/(2u)\)),

\[ \frac{\d f}{\d t}=-\frac{\tanh(u/2)}{4u},\qquad \frac{\d^2f}{\d t^2}=\frac{\tanh(u/2)-\frac u2\operatorname{sech}^2(u/2)}{8u^3}. \]

Let \(h(v)\)\({}=\tanh v-v\operatorname{sech}^2v\). Then \(h(0)=0\) and \(h'(v)\)\({}=\operatorname{sech}^2v\)\({}-\operatorname{sech}^2v\)\({}+2v\operatorname{sech}^2v\tanh v\)\({}\ge0\), so \(h\ge0\) for \(v\ge0\), and \(f\) is convex in \(t\) (strictly for \(t\gt0\)).

The tangent in \(t\). A convex function of \(t\) lies above its tangent at \(t_0=\xi^2\): \(f(x)\)\({}\ge f(\xi)+\frac{\d f}{\d t}\big|_{t=\xi^2}(x^2-\xi^2)\)\({}=f(\xi)-\lambda(\xi)(x^2-\xi^2)\), since \(\frac{\d f}{\d t}\big|_{\xi^2}\)\({}=-\frac{\tanh(\xi/2)}{4\xi}\)\({}=-\lambda(\xi)\). The identity \(\tanh(\xi/2)\)\({}=\frac{1-e^{-\xi}}{1+e^{-\xi}}\)\({}=2\sigma(\xi)-1\) gives the second form of \(\lambda\).

Back to the sigmoid. By (10.138) at \(x=\xi\), \(f(\xi)=\ln\sigma(\xi)-\xi/2\). So \(\ln\sigma(x)\)\({}=\frac x2+f(x)\)\({}\ge\ln\sigma(\xi)\)\({}+\frac{x-\xi}2-\lambda(\xi)(x^2-\xi^2)\); exponentiate. Equality holds when \(x^2=\xi^2\), and both sides are invariant under \(\xi\to-\xi\) (\(\sigma(-\xi)e^{\xi/2}\)\({}=\sigma(\xi)e^{-\xi/2}\)), so we may take \(\xi\ge0\).

Properties of \(\lambda\). \(\tanh(\xi/2)/\xi\) is even; \(\tanh(\xi/2)\approx\xi/2\) near 0 gives \(\lambda(0)=\frac18\); and \(\lambda'(\xi)\)\({}=\frac{2\xi\operatorname{sech}^2(\xi/2)-4\tanh(\xi/2)}{16\xi^2}\)\({}=-\frac{h(\xi/2)}{4\xi^2}\)\({}\lt0\) for \(\xi\gt0\). ∎

The book's conjugate-function route. With \(\eta\) the slope in \(x^2\), (10.129) reads

\[ g(\eta)=\max_{x^2}\big\{\eta x^2-f(\sqrt{x^2})\big\}, \tag{10.139} \]

with stationarity condition

\[ 0=\eta-\frac{\d x}{\d x^2}\frac{\d}{\d x}f(x)=\eta+\frac1{4x}\tanh\Big(\frac x2\Big). \tag{10.140} \]

Writing \(\xi\) for the contact point,

\[ \eta=-\frac1{4\xi}\tanh\Big(\frac\xi2\Big)=-\frac1{2\xi}\Big[\sigma(\xi)-\frac12\Big]=-\lambda(\xi), \tag{10.141} \] \[ g(\lambda(\xi))=-\lambda(\xi)\xi^2-f(\xi)=-\lambda(\xi)\xi^2+\ln\big(e^{\xi/2}+e^{-\xi/2}\big) \tag{10.142} \] \[ f(x)\ge-\lambda(\xi)x^2-g(\lambda(\xi))=-\lambda(\xi)x^2+\lambda(\xi)\xi^2-\ln\big(e^{\xi/2}+e^{-\xi/2}\big), \tag{10.143} \]

which is the tangent inequality above.

Erratum. In the first printing (10.141) defines the variational parameter as \(\lambda\) (with sign conventions that clash with Jaakkola and Jordan), and (10.142)–(10.143) are correspondingly different. The corrected forms, used above, write \(\eta=-\lambda(\xi)\) and \(g(\lambda(\xi))\).

I checked (10.144) numerically on a fine grid of \(x\in[-30,30]\) for several \(\xi\): the gap \(\sigma(x)-\text{bound}\) is never negative and vanishes at \(\pm\xi\). The bound is only for two classes; for the softmax, Gibbs (1997) proposed a Gaussian-form construction that is conjectured, but not proved, to be a bound.

Bounds inside integrals. The typical use is an integral such as

\[ I=\int\sigma(a)\,p(a)\,\d a \tag{10.145} \]

with \(p(a)\) Gaussian (a predictive probability, say). Writing (10.144) as \(\sigma(a)\ge f(a,\xi)\), the integrand of

\[ I\ge\int f(a,\xi)\,p(a)\,\d a=F(\xi) \tag{10.146} \]

is a product of two exponential-quadratic functions, so \(F\) is available in closed form, and we choose \(\xi^\star\) to maximise it. The bound is exact for one value of \(a\) only, so \(\xi^\star\) is a compromise weighted by \(p(a)\), and \(F(\xi^\star)\lt I\) in general.

Proof

\(F(\xi)\)\({}=\sigma(\xi)e^{-\xi/2+\lambda\xi^2}\int e^{a/2-\lambda a^2}\N(a\given\mu,s^2)\d a\), and the exponent of the integrand is \(-\frac P2a^2+ha-\frac{\mu^2}{2s^2}\). The Gaussian integral \(\int e^{-\frac P2a^2+ha}\d a\)\({}=\sqrt{2\pi/P}\,e^{h^2/(2P)}\), times \((2\pi s^2)^{-1/2}\), gives \((s^2P)^{-1/2}\)\({}=(1+2\lambda s^2)^{-1/2}\).

For the maximiser, treat \(a\) as a latent variable and use the EM argument of Section 9.4: for \(\tilde q(a)\propto f(a,\xi^{\rm old})p(a)\), \(\ln F(\xi)\)\({}\ge\E_{\tilde q}[\ln f(a,\xi)]+\const\) with equality at \(\xi=\xi^{\rm old}\), and \(\E_{\tilde q}[\ln f(a,\xi)]\)\({}=\ln\sigma(\xi)-\frac\xi2\)\({}-\lambda(\xi)(\E_{\tilde q}[a^2]-\xi^2)+\const\), whose \(\xi\)-derivative is \(-\lambda'(\xi)(\E_{\tilde q}[a^2]-\xi^2)\) (the same computation as (10.162) below). Since the bound touches \(\ln F\) at \(\xi^{\rm old}\), the derivatives agree there, so \(\frac{\d}{\d\xi}\ln F\)\({}=-\lambda'(\xi)(\E_{\tilde q}[a^2]-\xi^2)\) with \(\tilde q\) computed at \(\xi\): zero iff \(\xi^2=\E_{\tilde q}[a^2]\). ∎

For example, with \(\mu=0.7\) and \(s^2=2.3\) I computed \(I=0.6210\) by quadrature, while the optimised bound is \(F(\xi^\star)=0.5974\) at \(\xi^\star=1.784\), found both by the fixed-point iteration and by direct maximisation. The bound is useful mostly inside larger computations, as in the next section.

10.6 Variational logistic regression

Return to Bayesian logistic regression (Section 4.5), where the Laplace approximation gave a Gaussian posterior. The variational treatment of Jaakkola and Jordan (2000) also gives a Gaussian, but it is fitted by maximising a rigorous lower bound on the evidence, with extra flexibility (one parameter per data point), and is more accurate than Laplace.

10.6.1 Variational posterior distribution

Take the Gaussian prior \(p(\mathbf{w})\)\({}=\N(\mathbf{w}\given\mathbf{m}_0,\mathbf{S}_0)\) (4.140) with fixed \(\mathbf{m}_0,\mathbf{S}_0\) (Section 10.6.3 learns the prior precision). The evidence is

\[ p(\mathbf{t})=\int p(\mathbf{t}\given\mathbf{w})p(\mathbf{w})\,\d\mathbf{w}=\int\Big[\prod_{n=1}^{N}p(t_n\given\mathbf{w})\Big]p(\mathbf{w})\,\d\mathbf{w}, \tag{10.147} \]

and with \(a=\mathbf{w}^\T\bm\phi\) each likelihood term can be written as

\[ \begin{aligned}p(t\given\mathbf{w})&=\sigma(a)^t\{1-\sigma(a)\}^{1-t}=\Big(\frac1{1+e^{-a}}\Big)^t\Big(1-\frac1{1+e^{-a}}\Big)^{1-t}\\ &=e^{at}\frac{e^{-a}}{1+e^{-a}}=e^{at}\sigma(-a).\end{aligned} \tag{10.148} \]

(The third equality: \(1-\sigma(a)=e^{-a}/(1+e^{-a})\), so the product is \(e^{-a(1-t)}/(1+e^{-a})\)\({}=e^{at}\cdot e^{-a}/(1+e^{-a})\), and \(e^{-a}/(1+e^{-a})\)\({}=1/(e^a+1)\)\({}=\sigma(-a)\).) Restating (10.144) for convenience,

\[ \sigma(z)\ge\sigma(\xi)\exp\big\{(z-\xi)/2-\lambda(\xi)(z^2-\xi^2)\big\} \tag{10.149} \] \[ \lambda(\xi)=\frac1{2\xi}\Big[\sigma(\xi)-\frac12\Big]. \tag{10.150} \]

Apply (10.149) with \(z=-a\) (and \((-a)^2=a^2\)):

\[ p(t\given\mathbf{w})=e^{at}\sigma(-a)\ge e^{at}\sigma(\xi)\exp\big\{-(a+\xi)/2-\lambda(\xi)(a^2-\xi^2)\big\}. \tag{10.151} \]

The best \(\xi\) depends on \(a=\mathbf{w}^\T\bm\phi_n\), so each data point gets its own \(\xi_n\). Multiplying by the prior,

\[ p(\mathbf{t},\mathbf{w})=p(\mathbf{t}\given\mathbf{w})p(\mathbf{w})\ge h(\mathbf{w},\bm\xi)p(\mathbf{w}) \tag{10.152} \] \[ h(\mathbf{w},\bm\xi)=\prod_{n=1}^{N}\sigma(\xi_n)\exp\big\{\mathbf{w}^\T\bm\phi_nt_n-(\mathbf{w}^\T\bm\phi_n+\xi_n)/2-\lambda(\xi_n)([\mathbf{w}^\T\bm\phi_n]^2-\xi_n^2)\big\}. \tag{10.153} \]

Normalising the left side of (10.152) would give the exact posterior, which is intractable; we normalise the right side instead. The resulting \(q(\mathbf{w})\) is a density, but no longer a bound. Taking logarithms (an increasing function preserves the inequality),

\[ \begin{aligned}\ln\{p(\mathbf{t}\given\mathbf{w})p(\mathbf{w})\}\ge\ln p(\mathbf{w})+\sum_{n=1}^{N}\big\{&\ln\sigma(\xi_n)+\mathbf{w}^\T\bm\phi_nt_n\\ &-(\mathbf{w}^\T\bm\phi_n+\xi_n)/2-\lambda(\xi_n)([\mathbf{w}^\T\bm\phi_n]^2-\xi_n^2)\big\},\end{aligned} \tag{10.154} \]

and as a function of \(\mathbf{w}\) the right side is

\[ -\frac12(\mathbf{w}-\mathbf{m}_0)^\T\mathbf{S}_0^{-1}(\mathbf{w}-\mathbf{m}_0)+\sum_{n=1}^{N}\big\{\mathbf{w}^\T\bm\phi_n(t_n-1/2)-\lambda(\xi_n)\mathbf{w}^\T(\bm\phi_n\bm\phi_n^\T)\mathbf{w}\big\}+\const. \tag{10.155} \]

Proof

Expand (10.155): the quadratic terms are \(-\frac12\mathbf{w}^\T\{\mathbf{S}_0^{-1}\)\({}+2\sum_n\lambda(\xi_n)\bm\phi_n\bm\phi_n^\T\}\mathbf{w}\) (the factor 2 because \(-\lambda\,\mathbf{w}^\T\mathbf{A}\mathbf{w}\)\({}=-\frac12\mathbf{w}^\T(2\lambda\mathbf{A})\mathbf{w}\)), and the linear terms are \(\mathbf{w}^\T\{\mathbf{S}_0^{-1}\mathbf{m}_0\)\({}+\sum_n(t_n-\frac12)\bm\phi_n\}\); the terms \(\mathbf{w}^\T\bm\phi_nt_n\)\({}-\mathbf{w}^\T\bm\phi_n/2\) combine to \(\mathbf{w}^\T\bm\phi_n(t_n-\frac12)\). The matrix is positive definite (\(\lambda\gt0\)), so completing the square gives (10.156)–(10.158). ∎

Each data point acts as a Gaussian pseudo-observation of \(a_n=\mathbf{w}^\T\bm\phi_n\) with precision \(2\lambda(\xi_n)\) and pseudo-target \((t_n-\frac12)/(2\lambda(\xi_n))\). The Laplace approximation also adds a rank-one precision per point, \(\sigma_n(1-\sigma_n)\bm\phi_n\bm\phi_n^\T\) at the mode; here the curvature \(2\lambda(\xi_n)\) is chosen by optimising the evidence bound instead. Because the bound (10.154) is a sum of per-point terms, the posterior can also be updated one point at a time, \(\mathbf{S}_n^{-1}\)\({}=\mathbf{S}_{n-1}^{-1}\)\({}+2\lambda(\xi_n)\bm\phi_n\bm\phi_n^\T\) and \(\mathbf{S}_n^{-1}\mathbf{m}_n\)\({}=\mathbf{S}_{n-1}^{-1}\mathbf{m}_{n-1}\)\({}+(t_n-\frac12)\bm\phi_n\), with \(\xi_n\) optimised when its point arrives (Exercise 10.32). Predictions use the same probit approximation as the Laplace method (Section 4.5.2): \(p(t=1\given\bm\phi)\simeq\sigma\big(\kappa(\sigma_a^2)\mu_a\big)\) with \(\mu_a=\mathbf{m}_N^\T\bm\phi\), \(\sigma_a^2\)\({}=\bm\phi^\T\mathbf{S}_N\bm\phi\) and \(\kappa(s^2)=(1+\pi s^2/8)^{-1/2}\).

10.6.2 Optimising the variational parameters

The parameters \(\xi_n\) are chosen to maximise the bound on the log evidence obtained by substituting (10.152) into (10.147):

\[ \ln p(\mathbf{t})=\ln\int p(\mathbf{t}\given\mathbf{w})p(\mathbf{w})\,\d\mathbf{w}\ge\ln\int h(\mathbf{w},\bm\xi)p(\mathbf{w})\,\d\mathbf{w}=\LB(\bm\xi). \tag{10.159} \]

As with the hyperparameters of linear regression (Section 3.5) there are two routes: treat \(\mathbf{w}\) as a latent variable and use EM, or do the Gaussian integral in closed form and maximise directly. They give the same fixed point.

EM. The E step computes \(q(\mathbf{w})\) of (10.156) at \(\bm\xi^{\rm old}\). The M step maximises the expected complete-data log likelihood

\[ Q(\bm\xi,\bm\xi^{\rm old})=\E\big[\ln\{h(\mathbf{w},\bm\xi)p(\mathbf{w})\}\big] \tag{10.160} \] \[ Q(\bm\xi,\bm\xi^{\rm old})=\sum_{n=1}^{N}\big\{\ln\sigma(\xi_n)-\xi_n/2-\lambda(\xi_n)(\bm\phi_n^\T\E[\mathbf{w}\mathbf{w}^\T]\bm\phi_n-\xi_n^2)\big\}+\const \tag{10.161} \]

(the expectation is under \(q(\mathbf{w})\) at \(\bm\xi^{\rm old}\); \(p(\mathbf{w})\) and the terms linear in \(\mathbf{w}\) do not involve \(\bm\xi\)). Setting \(\partial Q/\partial\xi_n=0\):

\[ 0=\lambda'(\xi_n)\big(\bm\phi_n^\T\E[\mathbf{w}\mathbf{w}^\T]\bm\phi_n-\xi_n^2\big) \tag{10.162} \] \[ (\xi_n^{\rm new})^2=\bm\phi_n^\T\E[\mathbf{w}\mathbf{w}^\T]\bm\phi_n=\bm\phi_n^\T\big(\mathbf{S}_N+\mathbf{m}_N\mathbf{m}_N^\T\big)\bm\phi_n. \tag{10.163} \]

Erratum. The first printing writes \(\ln h(\mathbf{w},\bm\xi)p(\mathbf{w})\) in (10.160), without the braces.

Direct maximisation. The integral in (10.159) can be done in closed form:

\[ \begin{aligned}\LB(\bm\xi)&=\frac12\ln\frac{|\mathbf{S}_N|}{|\mathbf{S}_0|}+\frac12\mathbf{m}_N^\T\mathbf{S}_N^{-1}\mathbf{m}_N-\frac12\mathbf{m}_0^\T\mathbf{S}_0^{-1}\mathbf{m}_0\\ &\quad+\sum_{n=1}^{N}\Big\{\ln\sigma(\xi_n)-\frac12\xi_n+\lambda(\xi_n)\xi_n^2\Big\}.\end{aligned} \tag{10.164} \]

Erratum. The first printing has several sign errors in (10.164); the form above is the corrected one.

10.6.3 Inference of hyperparameters

Now let the prior precision be unknown, combining the global and local methods in one bound (as Bishop and Svensén 2003 did for hierarchical mixtures of experts):

\[ p(\mathbf{w}\given\alpha)=\N(\mathbf{w}\given\mathbf{0},\alpha^{-1}\mathbf{I}) \tag{10.165} \] \[ p(\alpha)=\Gam(\alpha\given a_0,b_0) \tag{10.166} \] \[ p(\mathbf{t})=\iint p(\mathbf{w},\alpha,\mathbf{t})\,\d\mathbf{w}\,\d\alpha \tag{10.167} \] \[ p(\mathbf{w},\alpha,\mathbf{t})=p(\mathbf{t}\given\mathbf{w})\,p(\mathbf{w}\given\alpha)\,p(\alpha). \tag{10.168} \]

For any \(q(\mathbf{w},\alpha)\), Result 10.1 gives

\[ \ln p(\mathbf{t})=\LB(q)+\KL(q\|p) \tag{10.169} \] \[ \LB(q)=\iint q(\mathbf{w},\alpha)\ln\left\{\frac{p(\mathbf{w},\alpha,\mathbf{t})}{q(\mathbf{w},\alpha)}\right\}\d\mathbf{w}\,\d\alpha \tag{10.170} \] \[ \KL(q\|p)=-\iint q(\mathbf{w},\alpha)\ln\left\{\frac{p(\mathbf{w},\alpha\given\mathbf{t})}{q(\mathbf{w},\alpha)}\right\}\d\mathbf{w}\,\d\alpha. \tag{10.171} \]

\(\LB(q)\) is still intractable because of \(p(\mathbf{t}\given\mathbf{w})\), so replace it by its lower bound \(h(\mathbf{w},\bm\xi)\) of (10.152)–(10.153); inside the logarithm this can only lower \(\LB\):

\[ \ln p(\mathbf{t})\ge\LB(q)\ge\tilde\LB(q,\bm\xi)=\iint q(\mathbf{w},\alpha)\ln\left\{\frac{h(\mathbf{w},\bm\xi)p(\mathbf{w}\given\alpha)p(\alpha)}{q(\mathbf{w},\alpha)}\right\}\d\mathbf{w}\,\d\alpha \tag{10.172} \]

and assume \(q(\mathbf{w},\alpha)\)\({}=q(\mathbf{w})q(\alpha)\) (10.173).

Erratum. The first printing has \(a_0,b_0\) instead of \(a_N,b_N\) on the right of (10.177), the transposes misplaced in (10.183), and refers to (10.159) instead of (10.160) when introducing (10.181).

10.7 Expectation propagation

Expectation propagation (EP; Minka 2001a, 2001b) is also based on a KL divergence, but the reverse one, \(\KL(p\|q)\), and only locally. Its basic step is the following fact about exponential families. Let

\[ q(\mathbf{z})=h(\mathbf{z})\,g(\bm\eta)\exp\big\{\bm\eta^\T\mathbf{u}(\mathbf{z})\big\} \tag{10.184} \]

(the form (2.194)). As a function of \(\bm\eta\), for a fixed \(p\),

\[ \KL(p\|q)=-\ln g(\bm\eta)-\bm\eta^\T\E_{p(\mathbf{z})}[\mathbf{u}(\mathbf{z})]+\const. \tag{10.185} \]

Proof

\(\KL(p\|q)\)\({}=\int p\ln p-\int p\ln h\)\({}-\ln g(\bm\eta)-\bm\eta^\T\E_p[\mathbf{u}]\), and the first two terms do not depend on \(\bm\eta\): this is (10.185). Its gradient is \(-\nabla\ln g(\bm\eta)-\E_p[\mathbf{u}]\), giving (10.186). By (2.226), \(-\nabla\ln g(\bm\eta)\)\({}=\E_q[\mathbf{u}]\), giving (10.187); and \(-\nabla\nabla\ln g(\bm\eta)\)\({}=\cov_q[\mathbf{u}]\succeq\mathbf{0}\) (toolkit: exponential families), so (10.185) is convex in \(\bm\eta\) and every stationary point is a global minimum (unique when the family is minimal). For the Gaussian, matching \(\E[\mathbf{z}]\) and \(\E[\mathbf{z}\mathbf{z}^\T]\) is matching the mean and the covariance. A direct proof by differentiation for the Gaussian is the solution of Exercise 10.4. ∎

This is the "moment matching" seen in Figure 10.3(a) and in the first row of the table after Result 10.2.

The EP approximation. Many joint densities are products of factors,

\[ p(\mathcal{D},\bm\theta)=\prod_if_i(\bm\theta), \tag{10.188} \]

for example \(f_0(\bm\theta)=p(\bm\theta)\) and \(f_n(\bm\theta)\)\({}=p(\mathbf{x}_n\given\bm\theta)\) for iid data, or the conditionals of a directed graph, or the potentials of an undirected one. Then

\[ p(\bm\theta\given\mathcal{D})=\frac1{p(\mathcal{D})}\prod_if_i(\bm\theta) \tag{10.189} \] \[ p(\mathcal{D})=\int\prod_if_i(\bm\theta)\,\d\bm\theta. \tag{10.190} \]

EP approximates the posterior by a product of the same shape, with each factor replaced by a member of an exponential family (so that \(q\) itself is in that family, e.g. Gaussian):

\[ q(\bm\theta)=\frac1Z\prod_i\tilde f_i(\bm\theta). \tag{10.191} \]

Ideally we would minimise

\[ \KL(p\|q)=\KL\Big(\frac1{p(\mathcal{D})}\prod_if_i(\bm\theta)\,\Big\|\,\frac1Z\prod_i\tilde f_i(\bm\theta)\Big), \tag{10.192} \]

but that needs expectations under the true posterior. Fitting each \(\tilde f_i\) to \(f_i\) in isolation is easy and non-iterative but can be poor, because a factor need only be accurate where the posterior has mass. EP fits each factor in the context of the others: to refine \(\tilde f_j\), make

\[ q^{\rm new}(\bm\theta)\propto\tilde f_j(\bm\theta)\prod_{i\ne j}\tilde f_i(\bm\theta) \tag{10.193} \]

as close as possible to

\[ f_j(\bm\theta)\prod_{i\ne j}\tilde f_i(\bm\theta) \tag{10.194} \]

with the other factors fixed. Concretely, remove \(\tilde f_j\) to get the cavity distribution

\[ q^{\setminus j}(\bm\theta)=\frac{q(\bm\theta)}{\tilde f_j(\bm\theta)}, \tag{10.195} \]

multiply by the true factor to get the tilted distribution

\[ \frac1{Z_j}f_j(\bm\theta)\,q^{\setminus j}(\bm\theta) \tag{10.196} \] \[ Z_j=\int f_j(\bm\theta)\,q^{\setminus j}(\bm\theta)\,\d\bm\theta, \tag{10.197} \]

and choose \(q^{\rm new}\) in the family to minimise

\[ \KL\Big(\frac{f_j(\bm\theta)q^{\setminus j}(\bm\theta)}{Z_j}\,\Big\|\,q^{\rm new}(\bm\theta)\Big), \tag{10.198} \]

which by Result 10.27 means matching the moments of the tilted distribution (for any exponential family whose moments can be related to the derivatives of its normaliser, (2.226)). The new factor is then obtained by division,

\[ \tilde f_j(\bm\theta)=K\frac{q^{\rm new}(\bm\theta)}{q^{\setminus j}(\bm\theta)}. \tag{10.199} \]

At a fixed point, \(q\) has exactly the moments of every tilted distribution \(q^{\setminus j}f_j\): this "moment consistency" is what EP achieves instead of the intractable global minimisation of (10.192).

Assumed density filtering (ADF, also called moment matching; see the references on p. 510 of the book) is the special case in which all factors except the first are initialised to 1 and each is updated once, in a single pass: suitable for on-line learning, but in batch use its answer depends on the arbitrary order of the data. Its evidence estimate is \(\prod_jZ_j\) (Exercise 10.36). EP revisits the factors until they agree, which removes the order dependence and improves accuracy.

Properties. EP is not guaranteed to converge. When \(q\) is in an exponential family and it does converge, the fixed point is a stationary point of a certain energy function (Minka 2001a), though the iterations need not decrease it; algorithms that minimise that energy directly converge, but are slower and more complex. Variational Bayes, by contrast, increases its bound at every step. Because EP minimises \(\KL(p\|q)\) locally, it averages over modes, which makes it a poor choice for mixture posteriors; for logistic-type models it is often more accurate than both local variational bounds and Laplace (Kuss and Rasmussen 2005).

Erratum. In the caption of Figure 10.14, "obtained by" should be inserted between "that" and "variational".

10.7.1 Example: the clutter problem

Following Minka (2001b): infer the mean \(\bm\theta\) of a \(D\)-dimensional Gaussian from observations embedded in background clutter,

\[ p(\mathbf{x}\given\bm\theta)=(1-w)\N(\mathbf{x}\given\bm\theta,\mathbf{I})+w\N(\mathbf{x}\given\mathbf{0},a\mathbf{I}) \tag{10.209} \] \[ p(\bm\theta)=\N(\bm\theta\given\mathbf{0},b\mathbf{I}) \tag{10.210} \] \[ p(\mathcal{D},\bm\theta)=p(\bm\theta)\prod_{n=1}^{N}p(\mathbf{x}_n\given\bm\theta) \tag{10.211} \]

with the clutter fraction \(w\) known; Minka uses \(a=10\), \(b=100\), \(w=0.5\). Expanding the product, the posterior is a mixture of \(2^N\) Gaussians (each data point is either signal or clutter), so exact inference costs \(O(2^N)\). For EP take \(f_0=p(\bm\theta)\), \(f_n=p(\mathbf{x}_n\given\bm\theta)\), and a spherical Gaussian

\[ q(\bm\theta)=\N(\bm\theta\given\mathbf{m},v\mathbf{I}), \tag{10.212} \]

so the factor approximations are exponential-quadratic,

\[ \tilde f_n(\bm\theta)=s_n\N(\bm\theta\given\mathbf{m}_n,v_n\mathbf{I}), \tag{10.213} \]

for \(n=1,\dots,N\), with \(\tilde f_0=p(\bm\theta)\). Here \(\N\) is only notation: \(v_n\) may be negative or infinite. Initialising \(\tilde f_n=1\) (\(s_n=(2\pi v_n)^{D/2}\), \(v_n\to\infty\), \(\mathbf{m}_n=\mathbf{0}\)) makes the initial \(q\) equal to the prior; the prior factor is already in the family and is never changed by an EP update (Exercise 10.37).

Erratum. In (10.217) and (10.218) the left sides read \(\mathbf{m}\) and \(v\) in the book (also in the printing used here); they are \(\mathbf{m}^{\rm new}\) and \(v^{\rm new}\). In (10.245) of Exercise 10.39 the term \(v^{\setminus n}D\) is missing in early printings.

EP has no guarantee of reaching a good fixed point: in a single sweep from the prior, two outlying clutter points processed first can pull the cavity towards themselves and steer the iteration to a poor solution, so it is worth comparing runs with different orders (or damping the updates).

Erratum. On p. 513, the first line after (10.224) should read "Examples of factor approximations".

10.7.2 Expectation propagation on graphs

Now let each factor depend on a subset of the variables, as in a factor graph (Section 8.4.3). We show that EP with a fully factorised \(q\) is loopy belief propagation (Minka 2001a). Start with the factor graph of Figure 8.51,

\[ p(\mathbf{x})=f_a(x_1,x_2)\,f_b(x_2,x_3)\,f_c(x_2,x_4). \tag{10.225} \]

Seek \(q\) with the same factorisation,

\[ q(\mathbf{x})\propto\tilde f_a(x_1,x_2)\,\tilde f_b(x_2,x_3)\,\tilde f_c(x_2,x_4), \tag{10.226} \]

with normalisation restored locally at the end, and restrict each approximate factor to factorise over its variables (Figure 10.18, right):

\[ q(\mathbf{x})\propto\tilde f_{a1}(x_1)\tilde f_{a2}(x_2)\tilde f_{b2}(x_2)\tilde f_{b3}(x_3)\tilde f_{c2}(x_2)\tilde f_{c4}(x_4). \tag{10.227} \]

Then \(q\) is fully factorised. To refine \(\tilde f_b(x_2,x_3)\)\({}=\tilde f_{b2}(x_2)\tilde f_{b3}(x_3)\), remove it and multiply by the exact factor:

\[ q^{\setminus b}(\mathbf{x})\propto\tilde f_{a1}(x_1)\tilde f_{a2}(x_2)\tilde f_{c2}(x_2)\tilde f_{c4}(x_4) \tag{10.228} \] \[ \hat p(\mathbf{x})=q^{\setminus b}(\mathbf{x})f_b(x_2,x_3)=\tilde f_{a1}(x_1)\tilde f_{a2}(x_2)\tilde f_{c2}(x_2)\tilde f_{c4}(x_4)f_b(x_2,x_3). \tag{10.229} \]

Erratum. The first printing has "=" instead of "\(\propto\)" in (10.228), and \(\tilde f\) is mistyped in the labels of Figure 10.18 (right).

By Result 10.6, minimising \(\KL(\hat p\|q^{\rm new})\) over fully factorised \(q^{\rm new}\) gives the product of the marginals of \(\hat p\):

\[ \hat p(x_1)\propto\tilde f_{a1}(x_1) \tag{10.230} \] \[ \hat p(x_2)\propto\tilde f_{a2}(x_2)\tilde f_{c2}(x_2)\sum_{x_3}f_b(x_2,x_3) \tag{10.231} \] \[ \hat p(x_3)\propto\sum_{x_2}\big\{f_b(x_2,x_3)\tilde f_{a2}(x_2)\tilde f_{c2}(x_2)\big\} \tag{10.232} \] \[ \hat p(x_4)\propto\tilde f_{c4}(x_4). \tag{10.233} \]

Dividing \(q^{\rm new}\) by \(q^{\setminus b}\) gives the refined factor:

\[ \tilde f_{b2}(x_2)\propto\sum_{x_3}f_b(x_2,x_3) \tag{10.234} \] \[ \tilde f_{b3}(x_3)\propto\sum_{x_2}\big\{f_b(x_2,x_3)\tilde f_{a2}(x_2)\tilde f_{c2}(x_2)\big\}. \tag{10.235} \]

These are exactly the belief-propagation messages with the variable-to-factor messages folded into the factor-to-variable messages (Section 8.4.4): \(\tilde f_{b2}(x_2)\) is \(\mu_{f_b\to x_2}(x_2)\) of (8.81), and substituting (8.78) into (8.79) gives (10.235), with \(\tilde f_{a2}=\mu_{f_a\to x_2}\) and \(\tilde f_{c2}=\mu_{f_c\to x_2}\), so that \(\tilde f_{b3}=\mu_{f_b\to x_3}\). EP sends the messages in both directions at once; refining one sub-factor at a time (for example only \(\tilde f_{b3}\)) gives the usual sum-product schedule, and on a tree the two-pass schedule gives exact marginals whatever the initialisation.

General factor graphs. For

\[ p(\bm\theta)=\prod_if_i(\bm\theta_i), \tag{10.236} \]

with \(\bm\theta_i\) the variables of factor \(i\), use the fully factorised approximation

\[ q(\bm\theta)\propto\prod_i\prod_k\tilde f_{ik}(\theta_k), \tag{10.237} \]

where \(\theta_k\) runs over the variable nodes (with \(\tilde f_{ik}\equiv1\) if \(\theta_k\notin\bm\theta_i\)). To refine \(\tilde f_{jl}(\theta_l)\), remove \(\tilde f_j(\bm\theta_j)\),

\[ q^{\setminus j}(\bm\theta)\propto\prod_{i\ne j}\prod_k\tilde f_{ik}(\theta_k), \tag{10.238} \]

and take the marginal on \(\theta_l\) of

\[ q^{\setminus j}(\bm\theta)f_j(\bm\theta_j). \tag{10.239} \]

Proof

Marginalise (10.239) onto \(\theta_l\). Split the factors \(\tilde f_{ik}(\theta_k)\), \(i\ne j\), of \(q^{\setminus j}\) into three groups: those with \(k=l\), which are functions of \(\theta_l\) and come out of the sum as \(\prod_{i\ne j}\tilde f_{il}(\theta_l)\); those with \(\theta_k\in\bm\theta_j\), \(k\ne l\), which stay inside the sum over \(\bm\theta_j\setminus\theta_l\) together with \(f_j\); and those with \(\theta_k\notin\bm\theta_j\), which sum to constants. Moment matching for a fully factorised \(q\) sets the new \(\theta_l\)-factor of \(q\) to this marginal (Result 10.6). Dividing by the \(\theta_l\)-factor of \(q^{\setminus j}\), which is \(\prod_{i\ne j}\tilde f_{il}(\theta_l)\), leaves (10.240): the terms of the other factors cancel between numerator and denominator. ∎

So the sum-product algorithm is the special case of EP with the most factorised \(q\). Less factorised approximations (partially disconnected graphs) and refining groups of factors together can both be more accurate (Minka 2001b); choosing them well is an open problem. More generally, Minka (2005) shows that variational message passing, loopy belief propagation and EP, as well as tree-reweighted message passing (Wainwright et al. 2005), fractional belief propagation (Wiegerinck and Heskes 2003) and power EP (Minka 2004), all minimise local alpha-divergences (10.19) for different \(\alpha\) (Result 10.7).

Exercises

All 39 exercises of Chapter 10 (pp. 517–522), restated briefly, with full solutions (hidden). The www exercises are those with official solutions (Svensén and Bishop 2009); the others are solved here independently. Where the main text already contains the complete derivation, the solution gives the argument in short form and links to it.

1. Exercise 10.1 (www). Show that \(\ln p(\mathbf{X})=\LB(q)+\KL(q\|p)\) with \(\LB\) and \(\KL\) as in (10.3) and (10.4).

   Solution

   By the product rule, \(\ln\{p(\mathbf{X},\mathbf{Z})/q(\mathbf{Z})\}\)\({}=\ln\{p(\mathbf{Z}\given\mathbf{X})/q(\mathbf{Z})\}+\ln p(\mathbf{X})\). Multiply by \(q(\mathbf{Z})\) and integrate; since \(\int q=1\), \(\LB(q)\)\({}=-\KL(q\|p)+\ln p(\mathbf{X})\). Details and consequences: Result 10.1.

2. Exercise 10.2. Using \(\E[z_1]=m_1\), \(\E[z_2]=m_2\), solve (10.13) and (10.15) and show that for non-singular \(p\) the unique solution is \(\E[z_i]=\mu_i\).

   Solution

   With \(e_i=m_i-\mu_i\) the equations are \(e_1\)\({}=-(\Lambda_{12}/\Lambda_{11})e_2\) and \(e_2\)\({}=-(\Lambda_{21}/\Lambda_{22})e_1\), so \(e_1(1-\rho^2)=0\) with \(\rho^2\)\({}=\Lambda_{12}^2/(\Lambda_{11}\Lambda_{22})\). Non-singularity of the positive definite \(\bm\Lambda\) means \(\Lambda_{11}\Lambda_{22}-\Lambda_{12}^2\)\({}\gt0\), i.e. \(\rho^2\lt1\); hence \(e_1=0\) and then \(e_2=0\). If \(\bm\Lambda\) were singular (\(\rho^2=1\)), every point of the line \(e_1\)\({}=-(\Lambda_{12}/\Lambda_{11})e_2\) would solve both equations, so uniqueness fails. See Result 10.5.

3. Exercise 10.3 (www). For factorised \(q\), use a Lagrange multiplier to show that minimising \(\KL(p\|q)\) over one factor \(q_j\) gives (10.17).

   Solution

   By (10.16), the terms of \(\KL(p\|q)\) that involve \(q_j\) are \(-\int p(\mathbf{Z})\ln q_j(\mathbf{Z}_j)\,\d\mathbf{Z}\)\({}=-\int F_j(\mathbf{Z}_j)\ln q_j(\mathbf{Z}_j)\,\d\mathbf{Z}_j\) with \(F_j(\mathbf{Z}_j)\)\({}=\int p(\mathbf{Z})\prod_{i\ne j}\d\mathbf{Z}_i\). Minimise \(-\int F_j\ln q_j+\lambda(\int q_j-1)\). The functional derivative is \(-F_j/q_j+\lambda=0\), so \(q_j=F_j/\lambda\), and integrating, \(\lambda=\int F_j\,\d\mathbf{Z}_j=1\). Hence \(q_j^\star=F_j=p_j\), the marginal. The second variation \(\int F_j(\delta q_j)^2/q_j^2\gt0\) shows it is a minimum. (A shorter route is the decomposition in Result 10.6.)

4. Exercise 10.4. For fixed \(p(\mathbf{x})\) and Gaussian \(q(\mathbf{x})\)\({}=\N(\mathbf{x}\given\bm\mu,\bm\Sigma)\), show by differentiation that \(\KL(p\|q)\) is minimised when \(\bm\mu\) and \(\bm\Sigma\) are the mean and covariance of \(p\).

   Solution

   Let \(\mathbf{m}=\E_p[\mathbf{x}]\), \(\mathbf{C}=\cov_p[\mathbf{x}]\) and \(\bm\Lambda=\bm\Sigma^{-1}\). Up to terms not involving \(q\),

   \[ \begin{aligned}\KL(p\|q)&=-\E_p[\ln\N(\mathbf{x}\given\bm\mu,\bm\Sigma)]=-\tfrac12\ln|\bm\Lambda|+\tfrac12\E_p[(\mathbf{x}-\bm\mu)^\T\bm\Lambda(\mathbf{x}-\bm\mu)]\\ &=-\tfrac12\ln|\bm\Lambda|+\tfrac12\Tr\big(\bm\Lambda\{\mathbf{C}+(\mathbf{m}-\bm\mu)(\mathbf{m}-\bm\mu)^\T\}\big).\end{aligned} \]

   The gradient in \(\bm\mu\) is \(-\bm\Lambda(\mathbf{m}-\bm\mu)\), zero iff \(\bm\mu=\mathbf{m}\). Then the gradient in \(\bm\Lambda\) is \(-\frac12\bm\Lambda^{-1}\)\({}+\frac12\mathbf{C}\) (toolkit: matrix calculus), zero iff \(\bm\Sigma=\mathbf{C}\). The function is convex in \((\bm\mu,\bm\Lambda)\) separately (a positive definite quadratic in \(\bm\mu\); \(-\ln|\bm\Lambda|\) is convex and the trace is linear), so this stationary point is the minimum. This is the Gaussian case of Result 10.27.

5. Exercise 10.5 (www). With \(q(\mathbf{z},\bm\theta)\)\({}=q_{\mathbf{z}}(\mathbf{z})q_{\bm\theta}(\bm\theta)\) and \(q_{\bm\theta}\)\({}=\delta(\bm\theta-\bm\theta_0)\), show that variational optimisation is EM: the E step optimises \(q_{\mathbf{z}}\), the M step maximises the expected complete-data log posterior over \(\bm\theta_0\).

   Solution

   E step. By (10.9), \(\ln q_{\mathbf{z}}(\mathbf{z})\)\({}=\E_{\bm\theta}[\ln p(\mathbf{X},\mathbf{z},\bm\theta)]+\const\)\({}=\ln p(\mathbf{X},\mathbf{z},\bm\theta_0)+\const\), so \(q_{\mathbf{z}}(\mathbf{z})\)\({}=p(\mathbf{z}\given\mathbf{X},\bm\theta_0)\), the E step. M step. \(\LB\)\({}=\iint q_{\mathbf{z}}q_{\bm\theta}\ln p(\mathbf{X},\mathbf{z},\bm\theta)\)\({}+\mathrm{H}[q_{\mathbf{z}}]\)\({}+\mathrm{H}[q_{\bm\theta}]\). The entropy of a point mass is \(-\infty\), but it does not depend on the location \(\bm\theta_0\); to be careful, take \(q_{\bm\theta}\)\({}=\N(\bm\theta_0,\epsilon\mathbf{I})\), whose entropy \(\frac d2\ln(2\pi e\epsilon)\) is free of \(\bm\theta_0\), and let \(\epsilon\to0\). So maximising \(\LB\) over \(\bm\theta_0\) maximises \(\E_{q_{\mathbf{z}}}[\ln p(\mathbf{X},\mathbf{z},\bm\theta_0)]\)\({}=\E_{q_{\mathbf{z}}}[\ln p(\mathbf{X},\mathbf{z}\given\bm\theta_0)]\)\({}+\ln p(\bm\theta_0)\): the M step of EM for the MAP estimate (Section 9.3).

6. Exercise 10.6. Show that \(\mathrm{D}_\alpha(p\|q)\)\({}\to\KL(p\|q)\) as \(\alpha\to1\) and \(\to\KL(q\|p)\) as \(\alpha\to-1\), by writing \(p^\epsilon\)\({}=\exp(\epsilon\ln p)\)\({}=1+\epsilon\ln p+O(\epsilon^2)\).

   Solution

   With \(\epsilon=(1-\alpha)/2\to0\), \(p^{(1+\alpha)/2}q^{(1-\alpha)/2}\)\({}=p\,(q/p)^\epsilon\)\({}=p\{1+\epsilon\ln(q/p)+O(\epsilon^2)\}\) and \(1-\alpha^2=4\epsilon(1-\epsilon)\). So \(\mathrm{D}_\alpha\)\({}=\frac{1}{\epsilon(1-\epsilon)}\{-\epsilon\int p\ln(q/p)\)\({}+O(\epsilon^2)\}\)\({}\to\int p\ln(p/q)\). The case \(\alpha\to-1\) is the same with the roles of \(p\) and \(q\) exchanged. Full statement and the other properties: Result 10.7.

7. Exercise 10.7. For the univariate Gaussian of Section 10.1.3, show that \(q_\mu\)\({}=\N(\mu\given\mu_N,\lambda_N^{-1})\) with (10.26)–(10.27) and \(q_\tau=\Gam(\tau\given a_N,b_N)\) with (10.29)–(10.30).

   Solution

   For \(q_\mu\): keep the terms of \(\E_\tau[\ln p]\) involving \(\mu\), \(-\frac{\E[\tau]}2\{\lambda_0(\mu-\mu_0)^2\)\({}+\sum_n(x_n-\mu)^2\}\), a quadratic with precision \((\lambda_0+N)\E[\tau]\) and mean \((\lambda_0\mu_0+N\bar x)/(\lambda_0+N)\). For \(q_\tau\): the terms in \(\tau\) are

   \[ \big(a_0-1+\tfrac N2+\tfrac12\big)\ln\tau-\tau\Big\{b_0+\tfrac12\E_\mu\Big[\sum_n(x_n-\mu)^2+\lambda_0(\mu-\mu_0)^2\Big]\Big\}, \]

   a Gamma kernel; the \(\frac12\ln\tau\) is the normaliser of \(p(\mu\given\tau)\). Full derivations: Gap 10.3 and Gap 10.4.

8. Exercise 10.8. Show that as \(N\to\infty\), \(q_\tau\) has mean equal to the inverse of the maximum likelihood variance and variance tending to zero.

   Solution

   From Gap 10.4, \(a_N=a_0+\frac{N+1}2\) and \(b_N\)\({}=b_0+\frac12\{\sum_n(x_n-\mu_N)^2\)\({}+\lambda_0(\mu_N-\mu_0)^2\)\({}+1/\E[\tau]\}\). As \(N\to\infty\), \(\mu_N\to\bar x\) and every term except \(\sum_n(x_n-\mu_N)^2=Ns^2+O(1)\) is \(O(1)\), where \(s^2=\frac1N\sum_n(x_n-\bar x)^2\) is the maximum likelihood variance. By (B.27), (B.28), \(\E[\tau]\)\({}=a_N/b_N\)\({}\to(N/2)/(Ns^2/2)\)\({}=1/s^2\) and \(\var[\tau]\)\({}=a_N/b_N^2\)\({}\approx2/(Ns^4)\)\({}\to0\).

9. Exercise 10.9. Derive (10.33) from \(\E[\tau]=a_N/b_N\), (10.26), (10.27), (10.29), (10.30).

   Solution

   With \(\mu_0=a_0=b_0=\lambda_0=0\), \(1/\E[\tau]\)\({}=b_N/a_N\)\({}=\frac{N}{N+1}\{\overline{x^2}\)\({}-\bar x^2+\frac{1}{N\E[\tau]}\}\); collecting the \(1/\E[\tau]\) terms gives \(1/\E[\tau]\)\({}=\overline{x^2}-\bar x^2\)\({}=\frac1N\sum_n(x_n-\bar x)^2\). See Gap 10.5 and, for the general prior and the 1/N versus 1/(N−1) question, Result 10.8 and the box after Gap 10.5.

10. Exercise 10.10 (www). Derive the decomposition (10.34).

    Solution

    Substitute \(p(\mathbf{Z},\mathbf{X},m)\)\({}=p(\mathbf{Z},m\given\mathbf{X})p(\mathbf{X})\) into (10.35): \(\LB\)\({}=\sum_m\sum_{\mathbf{Z}}q(\mathbf{Z}\given m)q(m)\ln\frac{p(\mathbf{Z},m\given\mathbf{X})}{q(\mathbf{Z}\given m)q(m)}\)\({}+\ln p(\mathbf{X})\), because \(\sum_m\sum_{\mathbf{Z}}q(\mathbf{Z}\given m)q(m)\)\({}=1\). Rearranging gives (10.34). (In the first printing \(\LB_m\) appears where \(\LB\) is meant; see the erratum in Section 10.1.4.)

11. Exercise 10.11 (www). With a Lagrange multiplier for \(\sum_mq(m)=1\), show that the maximum of (10.35) over \(q(m)\) is (10.36).

    Solution

    Summing over \(\mathbf{Z}\) first, \(\LB\)\({}=\sum_mq(m)\{\ln p(m)+\LB_m-\ln q(m)\}\). Add \(\lambda(\sum_mq(m)-1)\) and differentiate with respect to \(q(m)\): \(\ln p(m)+\LB_m-\ln q(m)-1+\lambda\)\({}=0\), so \(q(m)\)\({}=e^{\lambda-1}p(m)e^{\LB_m}\propto p(m)\exp\{\LB_m\}\), with \(\lambda\) fixed by normalisation. The objective is concave in \(q(m)\) (entropy plus linear terms), so this is the maximum. See Result 10.9.

12. Exercise 10.12. Starting from (10.41) and (10.9), verify that the optimal \(q^\star(\mathbf{Z})\) for the Bayesian mixture is (10.48).

    Solution

    Only \(\ln p(\mathbf{Z}\given\bm\pi)\) and \(\ln p(\mathbf{X}\given\mathbf{Z},\bm\mu,\bm\Lambda)\) involve \(\mathbf{Z}\); both are linear in \(z_{nk}\), so \(\ln q^\star(\mathbf{Z})\)\({}=\sum_{n,k}z_{nk}\ln\rho_{nk}+\const\) with (10.46). The normaliser factorises over \(n\), and each \(\mathbf{z}_n\) is one-of-\(K\), giving \(r_{nk}=\rho_{nk}/\sum_j\rho_{nj}\). All steps: Gap 10.6.

13. Exercise 10.13 (www). From (10.54), derive (10.59) and the parameters (10.60)–(10.63).

    Solution

    Collect the terms of (10.54) in \((\bm\mu_k,\bm\Lambda_k)\). Completing the square in \(\bm\mu_k\) gives \(q^\star(\bm\mu_k\given\bm\Lambda_k)\)\({}=\N(\mathbf{m}_k,(\beta_k\bm\Lambda_k)^{-1})\) with (10.60), (10.61). Subtracting its logarithm leaves \(\frac{\nu_0+N_k-D-1}{2}\ln|\bm\Lambda_k|\)\({}-\frac12\Tr(\mathbf{M}\bm\Lambda_k)\) with \(\mathbf{M}\) free of \(\bm\mu_k\) and equal to (10.62), a Wishart with \(\nu_k=\nu_0+N_k\). All steps, including why \(\mathbf{M}\) does not depend on \(\bm\mu_k\): Gap 10.8. Note that this gives \(\nu_k=\nu_0+N_k\), not the \(\nu_0+N_k+1\) printed in some copies of (10.63).

14. Exercise 10.14. Using (10.59), verify (10.64).

    Solution

    Condition on \(\bm\Lambda_k\): \(\E[(\mathbf{x}-\bm\mu_k)^\T\bm\Lambda_k(\mathbf{x}-\bm\mu_k)\given\bm\Lambda_k]\)\({}=(\mathbf{x}-\mathbf{m}_k)^\T\bm\Lambda_k(\mathbf{x}-\mathbf{m}_k)\)\({}+\Tr(\bm\Lambda_k(\beta_k\bm\Lambda_k)^{-1})\), and the trace is \(D/\beta_k\). Averaging over \(\bm\Lambda_k\) with \(\E[\bm\Lambda_k]=\nu_k\mathbf{W}_k\) gives \(D\beta_k^{-1}+\nu_k(\mathbf{x}-\mathbf{m}_k)^\T\mathbf{W}_k(\mathbf{x}-\mathbf{m}_k)\). See Gap 10.9.

15. Exercise 10.15. Using (B.17), show (10.69).

    Solution

    (B.17) gives \(\E[\pi_k]=\alpha_k/\hat\alpha\) for \(\Dir(\bm\pi\given\bm\alpha)\). Here \(\alpha_k=\alpha_0+N_k\) and \(\hat\alpha\)\({}=\sum_k(\alpha_0+N_k)\)\({}=K\alpha_0+N\), because \(\sum_kN_k\)\({}=\sum_n\sum_kr_{nk}\)\({}=\sum_n1\)\({}=N\). So \(\E[\pi_k]\)\({}=(\alpha_0+N_k)/(K\alpha_0+N)\).

16. Exercise 10.16 (www). Verify (10.71) and (10.72).

    Solution

    We have

    \[ \E[\ln p(\mathbf{X}\given\mathbf{Z},\bm\mu,\bm\Lambda)]=\tfrac12\sum_{n,k}r_{nk}\big\{\ln\tilde\Lambda_k-D\ln2\pi-D\beta_k^{-1}-\nu_k(\mathbf{x}_n-\mathbf{m}_k)^\T\mathbf{W}_k(\mathbf{x}_n-\mathbf{m}_k)\big\} \]

    by independence of \(q(\mathbf{Z})\) and \(q(\bm\mu,\bm\Lambda)\) and (10.64)–(10.65); then \(\sum_nr_{nk}(\mathbf{x}_n-\mathbf{m}_k)^\T\mathbf{W}_k(\mathbf{x}_n-\mathbf{m}_k)\)\({}=N_k\Tr(\mathbf{S}_k\mathbf{W}_k)\)\({}+N_k(\bar{\mathbf{x}}_k-\mathbf{m}_k)^\T\mathbf{W}_k(\bar{\mathbf{x}}_k-\mathbf{m}_k)\). And \(\E[\ln p(\mathbf{Z}\given\bm\pi)]\)\({}=\sum_{n,k}r_{nk}\E[\ln\pi_k]\). See Gap 10.11.

17. Exercise 10.17. Verify (10.73)–(10.77).

    Solution

    Each is the expectation of the logarithm of a Dirichlet, Gaussian or Wishart density under \(q\), evaluated with \(\E[\ln\pi_k]=\ln\tilde\pi_k\), \(\E[\ln|\bm\Lambda_k|]\)\({}=\ln\tilde\Lambda_k\), \(\E[\bm\Lambda_k]=\nu_k\mathbf{W}_k\) and the expected quadratic forms \(\E[(\bm\mu_k-\mathbf{m}_0)^\T\bm\Lambda_k(\bm\mu_k-\mathbf{m}_0)]\)\({}=D/\beta_k+\nu_k(\mathbf{m}_k-\mathbf{m}_0)^\T\mathbf{W}_k(\mathbf{m}_k-\mathbf{m}_0)\) and \(\E[(\bm\mu_k-\mathbf{m}_k)^\T\beta_k\bm\Lambda_k(\bm\mu_k-\mathbf{m}_k)]\)\({}=D\). Term by term: Gap 10.11.

18. Exercise 10.18. Write the bound (10.70) as a function of the parameters of the factors (10.48), (10.57), (10.59) and show that maximising it gives the same re-estimation equations as Section 10.2.1.

    Solution

    The bound is concave in the responsibilities of each \(n\) (a linear term plus an entropy), maximised at (10.49) by a Lagrange multiplier. For the parameter factors, the general identity \(\nabla_{\bm\eta}\LB\)\({}=\cov[\mathbf{u}](\tilde{\bm\eta}-\bm\eta)\) shows that the only stationary point is the natural parameter given by (10.9). Explicitly: the \(\bm\alpha\)-gradient is \(\sum_k(\alpha_0+N_k-\alpha_k)\{\psi'(\alpha_k)\delta_{kj}\)\({}-\psi'(\hat\alpha)\}\); the \(\mathbf{m}_k\)-gradient is \(\nu_k\mathbf{W}_k\{N_k(\bar{\mathbf{x}}_k-\mathbf{m}_k)\)\({}-\beta_0(\mathbf{m}_k-\mathbf{m}_0)\}\); the \(\beta_k\)-derivative is \(\frac D{2\beta_k^2}(N_k+\beta_0-\beta_k)\); and the \((\mathbf{W}_k,\nu_k)\) derivatives vanish only at (10.62) with \(\nu_k=\nu_0+N_k\). All details: Result 10.13.

19. Exercise 10.19. Derive the predictive density (10.81).

    Solution

    Integrate \(\bm\pi\) (giving \(\alpha_k/\hat\alpha\)), then \(\bm\mu_k\) (a linear-Gaussian marginal, \(\N(\hat{\mathbf{x}}\given\mathbf{m}_k,(\frac{\beta_k}{1+\beta_k}\bm\Lambda_k)^{-1})\)), then \(\bm\Lambda_k\) against the Wishart: the integrand is an unnormalised Wishart with \(\nu_k+1\) degrees of freedom, and the ratio of normalisers is a Student density with \(\nu_k+1-D\) degrees of freedom and precision (10.82). See Gap 10.12.

20. Exercise 10.20 (www). Show that for large \(N\) the variational posterior concentrates at the maximum likelihood solution, the responsibilities become the EM ones (use \(\psi(x)=\ln x+O(1/x)\), (10.241)), and the predictive density becomes a mixture of Gaussians.

    Solution

    The parameter posteriors have means tending to \(N_k/N\), \(\bar{\mathbf{x}}_k\), \(\mathbf{S}_k^{-1}\) and variances of order \(1/N_k\); the digamma asymptotics turn \(\ln\tilde\pi_k\) and \(\ln\tilde\Lambda_k\) into \(\ln(N_k/N)\) and \(\ln|\mathbf{S}_k^{-1}|\), and the Student components become Gaussians as their degrees of freedom grow. Full argument: Gap 10.10.

21. Exercise 10.21. Show that a \(K\)-component mixture has \(K!\) equivalent parameter settings.

    Solution

    A relabelling is a bijection \(\sigma\) of \(\{1,\dots,K\}\); there are \(K\) choices for \(\sigma(1)\), then \(K-1\) for \(\sigma(2)\), and so on: \(K!\) in total. Replacing \((\pi_k,\bm\mu_k,\bm\Lambda_k)\) by \((\pi_{\sigma(k)},\bm\mu_{\sigma(k)},\bm\Lambda_{\sigma(k)})\) permutes the terms of \(\sum_k\pi_k\N(\mathbf{x}\given\bm\mu_k,\bm\Lambda_k^{-1})\) and so leaves the density, hence the likelihood, unchanged. When the component parameters are all different (the generic case), different permutations give different parameter vectors, so each setting has exactly \(K!\) equivalent versions; when some components coincide, there are fewer.

22. Exercise 10.22. If \(q\) is localised at one mode and we replace it by the equal mixture of its \(K!\) relabelled copies, show that with negligible overlap the bound increases by \(\ln K!\).

    Solution

    The energy term \(\E[\ln p(\mathbf{X},\bm\theta)]\) is the same for every copy, because the joint is invariant under relabelling; the entropy of the mixture is the entropy of \(q\) plus the entropy \(\ln K!\) of the uniformly chosen copy when the copies do not overlap. So \(\LB(q_{\rm mix})=\LB(q)+\ln K!\). Details, and why \(\ln K!\) is an upper limit when copies overlap: Result 10.15.

23. Exercise 10.23 (www). With \(\bm\pi\) a parameter (no prior, no \(q(\bm\pi)\)), maximise the bound over \(\bm\pi\) subject to \(\sum_k\pi_k=1\) and obtain (10.83).

    Solution

    Only \(\E_{q(\mathbf{Z})}[\ln p(\mathbf{Z}\given\bm\pi)]\)\({}=\sum_{n,k}r_{nk}\ln\pi_k\)\({}=\sum_kN_k\ln\pi_k\) depends on \(\bm\pi\). With \(\lambda(\sum_k\pi_k-1)\) added, \(N_k/\pi_k+\lambda=0\), so \(\pi_k=-N_k/\lambda\); summing over \(k\), \(-\lambda=\sum_kN_k=N\). Hence \(\pi_k=N_k/N=\frac1N\sum_nr_{nk}\), (10.83). The objective is concave in \(\bm\pi\), so this is the maximum.

24. Exercise 10.24 (www). Would the singularities of maximum likelihood (Section 9.2.1) occur if the Bayesian mixture were fitted by MAP?

    Solution

    No. A singularity needs a component with \(\bm\mu_k=\mathbf{x}_n\) and \(|\bm\Lambda_k|\to\infty\): the log likelihood then grows like \(\frac12\ln|\bm\Lambda_k|\). The Wishart prior contributes \(-\frac12\Tr(\mathbf{W}_0^{-1}\bm\Lambda_k)\)\({}\le-\frac12\lambda_{\min}(\mathbf{W}_0^{-1})\Tr\bm\Lambda_k\), which goes to \(-\infty\) linearly, so the log posterior stays bounded above. Explicitly, for fixed responsibilities and \(\bm\mu_k\), the terms of the log posterior in \(\bm\Lambda_k\) are \(\frac{\nu_0+N_k-D}{2}\ln|\bm\Lambda_k|\)\({}-\frac12\Tr(\mathbf{M}_k\bm\Lambda_k)\) with \(\mathbf{M}_k\)\({}=\mathbf{W}_0^{-1}\)\({}+\beta_0(\bm\mu_k-\mathbf{m}_0)(\cdot)^\T\)\({}+\sum_nr_{nk}(\mathbf{x}_n-\bm\mu_k)(\cdot)^\T\succeq\mathbf{W}_0^{-1}\) (the power \(\frac{\nu_0+N_k-D}{2}\) collects \(\frac12\) from \(p(\bm\mu_k\given\bm\Lambda_k)\), \(\frac{\nu_0-D-1}{2}\) from the Wishart and \(\frac{N_k}{2}\) from the data). The maximiser is \(\bm\Lambda_k^{-1}\)\({}=\mathbf{M}_k/(\nu_0+N_k-D)\), which can never be singular because \(\mathbf{W}_0^{-1}\) is positive definite. (MAP is not invariant under reparameterisation, but it does remove the collapse.)

25. Exercise 10.25. The factorised approximation under-estimates posterior variance in some directions. Discuss the effect on the variational approximation to the evidence, how it varies with \(K\), and whether the variational mixture will tend to under- or over-estimate the number of components.

    Solution

    The bound falls short of \(\ln p(\mathbf{X}\given K)\) by \(\KL(q\|p)\), which is positive whenever the posterior has correlations that \(q(\mathbf{Z})q(\bm\pi)\prod_kq(\bm\mu_k,\bm\Lambda_k)\) cannot represent; a too-compact \(q\) (Section 10.1.2) misses posterior volume, so the bound under-estimates the evidence. With more components there are more parameters, more neglected dependencies (between the \(z_{nk}\) and the component parameters, and between components competing for the same points), and more equivalent modes, so the shortfall tends to grow with \(K\). Model comparison with \(\LB\) therefore penalises large \(K\) too much: the variational mixture tends to under-estimate the number of components. (The \(\ln K!\) correction of Result 10.15 removes only the part of the shortfall that is due to label switching.)

26. Exercise 10.26. Extend variational linear regression with a prior \(\Gam(\beta\given c_0,d_0)\) and \(q(\mathbf{w})q(\alpha)q(\beta)\): derive the three updates, the bound and the predictive distribution.

    Solution

    Updates. By (10.9):

    • \(\ln q(\mathbf{w})\)\({}=-\frac{\E[\beta]}2\|\mathbf{t}\)\({}-\bm\Phi\mathbf{w}\|^2\)\({}-\frac{\E[\alpha]}2\mathbf{w}^\T\mathbf{w}+\const\), so \(q(\mathbf{w})\)\({}=\N(\mathbf{m}_N,\mathbf{S}_N)\) with \(\mathbf{S}_N^{-1}\)\({}=\E[\alpha]\mathbf{I}\)\({}+\E[\beta]\bm\Phi^\T\bm\Phi\) and \(\mathbf{m}_N\)\({}=\E[\beta]\mathbf{S}_N\bm\Phi^\T\mathbf{t}\);
    • \(q(\alpha)=\Gam(a_N,b_N)\) exactly as in (10.94)–(10.95);
    • \(\ln q(\beta)\)\({}=(c_0-1+\frac N2)\ln\beta\)\({}-\beta\{d_0+\frac12\E\|\mathbf{t}\)\({}-\bm\Phi\mathbf{w}\|^2\}+\const\), so \(q(\beta)=\Gam(c_N,d_N)\) with \(c_N=c_0+\frac N2\) and \(d_N\)\({}=d_0+\frac12\{\|\mathbf{t}\)\({}-\bm\Phi\mathbf{m}_N\|^2\)\({}+\Tr(\bm\Phi^\T\bm\Phi\mathbf{S}_N)\}\).

    Bound. \(\LB\)\({}=\E[\ln p(\mathbf{t}\given\mathbf{w},\beta)]\)\({}+\E[\ln p(\mathbf{w}\given\alpha)]\)\({}+\E[\ln p(\alpha)]\)\({}+\E[\ln p(\beta)]+\mathrm{H}[q(\mathbf{w})]\)\({}+\mathrm{H}[q(\alpha)]\)\({}+\mathrm{H}[q(\beta)]\), where \(\E[\ln p(\mathbf{t}\given\mathbf{w},\beta)]\)\({}=\frac N2\{\psi(c_N)\)\({}-\ln d_N-\ln2\pi\}-\frac{c_N}{2d_N}\{\|\mathbf{t}\)\({}-\bm\Phi\mathbf{m}_N\|^2\)\({}+\Tr(\bm\Phi^\T\bm\Phi\mathbf{S}_N)\}\); the \(\beta\) prior and entropy terms are (10.110) and (10.112) with \((a_0,b_0,a_N,b_N)\) replaced by \((c_0,d_0,c_N,d_N)\); the others are (10.109)–(10.111).

    Predictive.

    \[ \begin{aligned}p(t\given\mathbf{x},\mathbf{t})&\simeq\iint\N(t\given\mathbf{w}^\T\bm\phi,\beta^{-1})q(\mathbf{w})q(\beta)\,\d\mathbf{w}\,\d\beta\\ &=\int\N\big(t\given\mathbf{m}_N^\T\bm\phi,\beta^{-1}+\bm\phi^\T\mathbf{S}_N\bm\phi\big)\Gam(\beta\given c_N,d_N)\,\d\beta.\end{aligned} \]

    Because \(q(\mathbf{w})\) does not scale with \(\beta\), this one-dimensional integral has no closed form; it can be done numerically, or approximated by \(\N(t\given\mathbf{m}_N^\T\bm\phi,\,d_N/c_N+\bm\phi^\T\mathbf{S}_N\bm\phi)\), which is accurate when \(q(\beta)\) is sharp (large \(N\)). (With the conjugate prior \(\mathbf{w}\given\alpha,\beta\sim\N(\mathbf{0},(\alpha\beta)^{-1}\mathbf{I})\) and an unfactorised \(q(\mathbf{w},\beta)\), the predictive would be a Student t density.)

27. Exercise 10.27. Show that the bound for the linear basis function model has the form (10.107) with the terms (10.108)–(10.112).

    Solution

    Each term is an expectation of a Gaussian or Gamma log density under \(q(\mathbf{w})q(\alpha)\), using \(\E[\mathbf{w}]=\mathbf{m}_N\), \(\E[\mathbf{w}\mathbf{w}^\T]\)\({}=\mathbf{m}_N\mathbf{m}_N^\T\)\({}+\mathbf{S}_N\), \(\E[\alpha]=a_N/b_N\), \(\E[\ln\alpha]=\psi(a_N)-\ln b_N\), and the Gaussian and Gamma entropies. Term by term (with the erratum \(\ln\Gamma(a_0)\) in (10.110)): Gap 10.15.

28. Exercise 10.28. Write the Bayesian mixture of Gaussians as a conjugate exponential-family model and recover (10.48), (10.57) and (10.59) from (10.115) and (10.119).

    Solution

    Complete-data likelihood. For one point,

    \[ \begin{aligned}\ln p(\mathbf{x}_n,\mathbf{z}_n\given\bm\pi,\bm\mu,\bm\Lambda)&=\sum_kz_{nk}\Big\{\ln\pi_k+\tfrac12\ln|\bm\Lambda_k|-\tfrac12\bm\mu_k^\T\bm\Lambda_k\bm\mu_k\Big\}\\ &\quad+\sum_kz_{nk}\mathbf{x}_n^\T(\bm\Lambda_k\bm\mu_k)+\sum_k\Tr\Big(z_{nk}\mathbf{x}_n\mathbf{x}_n^\T\big(-\tfrac12\bm\Lambda_k\big)\Big)-\tfrac D2\ln2\pi,\end{aligned} \]

    which is \(\bm\eta^\T\mathbf{u}(\mathbf{x}_n,\mathbf{z}_n)+\ln h\) with sufficient statistics \(\mathbf{u}\)\({}=\{z_{nk},\,z_{nk}\mathbf{x}_n,\,z_{nk}\mathbf{x}_n\mathbf{x}_n^\T\}_{k}\), natural parameters \(\bm\eta\)\({}=\{\ln\pi_k+\frac12\ln|\bm\Lambda_k|\)\({}-\frac12\bm\mu_k^\T\bm\Lambda_k\bm\mu_k,\ \bm\Lambda_k\bm\mu_k,\ -\frac12\bm\Lambda_k\}_k\), \(h=(2\pi)^{-D/2}\), and \(g(\bm\eta)=1\) (the density is normalised automatically because \(\sum_k\pi_k=1\); it is a curved exponential family).

    (10.115) gives (10.48). \(\ln q^\star(\mathbf{Z})\)\({}=\sum_{n,k}z_{nk}\{\E[\ln\pi_k]\)\({}+\frac12\E\ln|\bm\Lambda_k|\)\({}-\frac12\E[\bm\mu_k^\T\bm\Lambda_k\bm\mu_k]\)\({}+\mathbf{x}_n^\T\E[\bm\Lambda_k\bm\mu_k]\)\({}-\frac12\mathbf{x}_n^\T\E[\bm\Lambda_k]\mathbf{x}_n\}-\frac D2\ln2\pi+\const\). Since \(\E[(\mathbf{x}-\bm\mu)^\T\bm\Lambda(\mathbf{x}-\bm\mu)]\)\({}=\mathbf{x}^\T\E[\bm\Lambda]\mathbf{x}\)\({}-2\mathbf{x}^\T\E[\bm\Lambda\bm\mu]\)\({}+\E[\bm\mu^\T\bm\Lambda\bm\mu]\), the braces equal \(\ln\rho_{nk}\) of (10.46), and normalising over each \(\mathbf{z}_n\) gives (10.48).

    (10.118)–(10.121) give (10.57) and (10.59). The prior \(\Dir(\bm\pi\given\bm\alpha_0)\prod_k\N(\bm\mu_k\given\mathbf{m}_0,(\beta_0\bm\Lambda_k)^{-1})\)\(\Wish(\bm\Lambda_k\given\mathbf{W}_0,\nu_0)\) is conjugate: its logarithm, \(\sum_k(\alpha_0-1)\ln\pi_k\)\({}+\sum_k\{\frac{\nu_0-D}{2}\ln|\bm\Lambda_k|\)\({}-\frac{\beta_0}{2}\bm\mu_k^\T\bm\Lambda_k\bm\mu_k\)\({}+\beta_0\mathbf{m}_0^\T\bm\Lambda_k\bm\mu_k\)\({}-\frac12\Tr[(\mathbf{W}_0^{-1}+\beta_0\mathbf{m}_0\mathbf{m}_0^\T)\bm\Lambda_k]\}+\const\), is a linear combination of the same functions of the parameters that appear in \(\bm\eta\). (It has a separate pseudo-count for each block: \(\alpha_0\), \(\beta_0\), \(\nu_0\), rather than the single \(\nu_0\) of (10.114), but the argument of (10.117)–(10.118) is unchanged.) Adding \(\bm\eta^\T\sum_n\E[\mathbf{u}_n]\), with \(\sum_n\E[\mathbf{u}_n]\)\({}=\{N_k,\,N_k\bar{\mathbf{x}}_k,\,\sum_nr_{nk}\mathbf{x}_n\mathbf{x}_n^\T\}_k\), adds \(N_k\ln\pi_k\) to the \(\bm\pi\) part, giving \(\Dir(\alpha_0+N_k)\), (10.57); and adds \(\frac{N_k}{2}\ln|\bm\Lambda_k|\)\({}-\frac{N_k}{2}\bm\mu_k^\T\bm\Lambda_k\bm\mu_k\)\({}+N_k\bar{\mathbf{x}}_k^\T\bm\Lambda_k\bm\mu_k\)\({}-\frac12\Tr(\sum_nr_{nk}\mathbf{x}_n\mathbf{x}_n^\T\bm\Lambda_k)\) to each component, which are exactly the terms collected in Gap 10.8, giving the Gaussian–Wishart (10.59)–(10.63).

29. Exercise 10.29 (www). Show that \(\ln x\) is concave on \((0,\infty)\), find its dual function (10.133), and verify that (10.132) recovers \(\ln x\).

    Solution

    \((\ln x)''=-1/x^2\lt0\). By (10.133), \(g(\eta)=\min_x\{\eta x-\ln x\}\); for \(\eta\gt0\) the derivative \(\eta-1/x\) vanishes at \(x=1/\eta\) (a minimum, second derivative \(1/x^2\gt0\)), so \(g(\eta)=1-\ln(1/\eta)=1+\ln\eta\). (For \(\eta\le0\), \(g=-\infty\).) Then (10.132): \(\min_\eta\{\eta x-1-\ln\eta\}\), with derivative \(x-1/\eta=0\) at \(\eta=1/x\), gives \(1-1-\ln(1/x)=\ln x\).

30. Exercise 10.30. Show that \(\ln\sigma(x)\) is concave and derive the upper bound (10.137) from a first-order Taylor expansion.

    Solution

    \((\ln\sigma)'=1-\sigma\) and \((\ln\sigma)''\)\({}=-\sigma(1-\sigma)\)\({}\lt0\). A concave function lies below its tangent: \(\ln\sigma(x)\)\({}\le\ln\sigma(\xi)\)\({}+(1-\sigma(\xi))(x-\xi)\). With \(\eta=1-\sigma(\xi)\) the intercept is \(\ln\sigma(\xi)-\eta\xi=-g(\eta)\), the binary entropy, giving (10.136) and, exponentiating, (10.137). See Result 10.21. (The first printing says "second order"; it is first order.)

31. Exercise 10.31. Show that \(f(x)=-\ln(e^{x/2}+e^{-x/2})\) is concave in \(x\) and convex in \(x^2\); sketch it against \(x\) and \(x^2\); derive (10.144) from a first-order Taylor expansion in \(x^2\) at \(\xi^2\).

    Solution

    \(f''(x)\)\({}=-\frac14\operatorname{sech}^2(x/2)\)\({}\lt0\); as a function of \(t=x^2\), \(\d^2f/\d t^2\)\({}=\{\tanh(u/2)-\frac u2\operatorname{sech}^2(u/2)\}/(8u^3)\)\({}\ge0\) with \(u=\sqrt t\). Sketch: against \(x\), \(f\) is even and concave, with maximum \(f(0)=-\ln2\) and \(f(x)\approx-|x|/2\) for large \(|x|\) (a rounded inverted V); against \(t=x^2\) it decreases like \(-\sqrt t/2\), a convex curve starting at \(-\ln2\) with slope \(-\frac18\). The tangent in \(t\) at \(\xi^2\) gives \(f(x)\)\({}\ge f(\xi)-\lambda(\xi)(x^2-\xi^2)\), and \(\ln\sigma(x)=x/2+f(x)\) gives (10.144). Full proof: Result 10.22. The demo plots the resulting bound.

32. Exercise 10.32 (www). For sequential learning in variational logistic regression, show that a Gaussian posterior can be maintained with the bound (10.151), starting from the prior and optimising each \(\xi_n\) as its point is absorbed.

    Solution

    The log of the bound (10.154) is the log prior plus a sum of per-point terms, each quadratic in \(\mathbf{w}\). After \(n-1\) points the accumulated prior-plus-data terms are \(\ln\N(\mathbf{w}\given\mathbf{m}_{n-1},\mathbf{S}_{n-1})\) up to a constant, and acting as the prior for point \(n\):

    \[ \mathbf{S}_n^{-1}=\mathbf{S}_{n-1}^{-1}+2\lambda(\xi_n)\bm\phi_n\bm\phi_n^\T,\qquad \mathbf{m}_n=\mathbf{S}_n\big\{\mathbf{S}_{n-1}^{-1}\mathbf{m}_{n-1}+(t_n-\tfrac12)\bm\phi_n\big\}, \]

    which unrolls to (10.157)–(10.158). For the new point, alternate this update with \(\xi_n^2\)\({}=\bm\phi_n^\T(\mathbf{S}_n+\mathbf{m}_n\mathbf{m}_n^\T)\bm\phi_n\) (the one-point version of (10.163)) until \(\xi_n\) settles, then discard the point. Unlike the batch scheme, earlier \(\xi\)'s are never revisited.

33. Exercise 10.33. Differentiate \(Q(\bm\xi,\bm\xi^{\rm old})\) of (10.161) to obtain (10.163).

    Solution

    Using \((\ln\sigma)'=1-\sigma\) and \(2\xi\lambda(\xi)\)\({}=\sigma(\xi)-\frac12\), \(\partial Q/\partial\xi_n\)\({}=-\lambda'(\xi_n)(\bm\phi_n^\T\E[\mathbf{w}\mathbf{w}^\T]\bm\phi_n-\xi_n^2)\). Since \(\lambda'\lt0\) on \(\xi\gt0\), the unique maximiser is \(\xi_n^2\)\({}=\bm\phi_n^\T(\mathbf{S}_N+\mathbf{m}_N\mathbf{m}_N^\T)\bm\phi_n\). See Gap 10.18.

34. Exercise 10.34. Re-derive the \(\xi\) updates by setting \(\partial\LB(\bm\xi)/\partial\xi_n\)\({}=0\) in (10.164), using (3.117), (10.157) and (10.158).

    Solution

    Only \(\mathbf{S}_N\) depends on \(\bm\xi\), with \(\partial\mathbf{S}_N^{-1}/\partial\xi_n\)\({}=2\lambda'(\xi_n)\bm\phi_n\bm\phi_n^\T\). With (3.117), \(\partial\frac12\ln|\mathbf{S}_N|/\partial\xi_n\)\({}=-\lambda'\bm\phi_n^\T\mathbf{S}_N\bm\phi_n\); writing \(\mathbf{m}_N^\T\mathbf{S}_N^{-1}\mathbf{m}_N\)\({}=\mathbf{b}^\T\mathbf{S}_N\mathbf{b}\) with \(\mathbf{b}\) independent of \(\bm\xi\), its half-derivative is \(-\lambda'(\bm\phi_n^\T\mathbf{m}_N)^2\); and the last sum contributes \(\lambda'\xi_n^2\). Hence \(\partial\LB/\partial\xi_n\)\({}=\lambda'(\xi_n)\{\xi_n^2\)\({}-\bm\phi_n^\T(\mathbf{S}_N+\mathbf{m}_N\mathbf{m}_N^\T)\bm\phi_n\}\), zero at (10.163). See Gap 10.18.

35. Exercise 10.35. Derive (10.164) by substituting the prior and \(h(\mathbf{w},\bm\xi)\) into (10.159), completing the square and using the Gaussian normaliser.

    Solution

    \(h(\mathbf{w},\bm\xi)p(\mathbf{w})\)\({}=\exp\{-\frac12\mathbf{w}^\T\mathbf{S}_N^{-1}\mathbf{w}\)\({}+\mathbf{w}^\T\mathbf{S}_N^{-1}\mathbf{m}_N+c\}\) with \(c\)\({}=\sum_n\{\ln\sigma(\xi_n)\)\({}-\frac{\xi_n}2+\lambda(\xi_n)\xi_n^2\}-\frac12\mathbf{m}_0^\T\mathbf{S}_0^{-1}\mathbf{m}_0\)\({}-\frac12\ln|2\pi\mathbf{S}_0|\); the Gaussian integral is \(|2\pi\mathbf{S}_N|^{1/2}\exp\{\frac12\mathbf{m}_N^\T\mathbf{S}_N^{-1}\mathbf{m}_N\}\). Taking logs gives (10.164) in its corrected form. See Gap 10.18.

36. Exercise 10.36. For assumed density filtering, show that including factor \(f_j\) updates the evidence as \(p_j(\mathcal{D})\simeq p_{j-1}(\mathcal{D})Z_j\), and hence \(p(\mathcal{D})\simeq\prod_jZ_j\).

    Solution

    In ADF the approximate factors not yet processed equal 1, so after processing factors \(1,\dots,j-1\) the product of the approximate factors is \(p_{j-1}(\mathcal{D})\,q_{j-1}(\bm\theta)\) with \(q_{j-1}\) normalised (induction hypothesis; at the start, \(p_0=1\) and \(q_0\) is the first factor, e.g. the prior, if it is normalised). When \(f_j\) is processed, the cavity is \(q^{\setminus j}=q_{j-1}\) (because \(\tilde f_j=1\)), \(Z_j=\int f_jq_{j-1}\), and (10.207) gives \(\tilde f_j=Z_jq_j/q_{j-1}\). Then \(\prod_{i\le j}\tilde f_i\)\({}=p_{j-1}q_{j-1}\cdot Z_jq_j/q_{j-1}\)\({}=p_{j-1}Z_j\,q_j\), and integrating, \(p_j=p_{j-1}Z_j\) (10.242). By induction \(p(\mathcal{D})\simeq\prod_jZ_j\) (10.243). This is the chain rule \(p(\mathcal{D})\)\({}=\prod_jp(\mathbf{x}_j\given\mathbf{x}_{\lt j})\), with each one-step predictive \(\int f_j\,p(\bm\theta\given\mathbf{x}_{\lt j})\d\bm\theta\) approximated by \(Z_j\).

37. Exercise 10.37 (www). If \(f_0\) has the same exponential-family form as \(q\) and \(\tilde f_0\) is initialised to \(f_0\), show that an EP update leaves \(\tilde f_0\) unchanged.

    Solution

    Write \(q\)\({}=\frac1Zf_0\prod_{i\ne0}\tilde f_i\). Then \(q^{\setminus0}\)\({}=q/\tilde f_0\)\({}=\frac1Z\prod_{i\ne0}\tilde f_i\), and the tilted function is \(f_0q^{\setminus0}=q\), which is already in the family and normalised. Moment matching therefore returns \(q^{\rm new}=q\), and \(Z_0\)\({}=\int f_0q^{\setminus0}\)\({}=\int q\)\({}=1\). By (10.207), \(\tilde f_0\)\({}=Z_0q^{\rm new}/q^{\setminus0}\)\({}=q/q^{\setminus0}\)\({}=f_0\): unchanged. So a prior in the approximating family is included exactly once and never refined.

38. Exercise 10.38. For the clutter problem, derive (10.214)–(10.215) from the division (10.205) by completing the square, and show that \(Z_n\) is (10.216) using (2.115).

    Solution

    Dividing the Gaussians subtracts precisions and precision-times-means: \((v^{\setminus n})^{-1}\)\({}=v^{-1}-v_n^{-1}\) and \(\mathbf{m}^{\setminus n}\)\({}=v^{\setminus n}(\mathbf{m}/v-\mathbf{m}_n/v_n)\)\({}=\mathbf{m}+v^{\setminus n}v_n^{-1}(\mathbf{m}-\mathbf{m}_n)\). Then \(Z_n\) is a linear-Gaussian marginal plus a constant clutter term. See Gap 10.21.

39. Exercise 10.39. Prove (10.244)–(10.245), deduce (10.217)–(10.218), then derive (10.220)–(10.222) from (10.207) and (10.223) from (10.208).

    Solution

    Differentiating \(\ln Z\) with respect to the cavity mean and variance gives the tilted moments (10.244)–(10.245) (with the term \(v^{\setminus n}D\) that early printings omit); applying them to (10.216) gives (10.217)–(10.219). The site parameters follow from dividing \(q^{\rm new}\) by the cavity, and the evidence from integrating the product of all sites, including the prior. Complete derivations: Gap 10.22 and Gap 10.23; for the conventions behind (10.222) and (10.223) see the note on the evidence.

Demonstrations for this chapter

Each demonstration runs in your browser, uses a fixed random seed so that a class sees the same pictures, and links back to these notes. The ▶ Demo boxes above open them already set to the case being discussed.

References

1. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer, New York. Chapter 10, pp. 461–522, and the official errata list.
2. Ali, S. M. and Silvey, S. D. (1966). A general class of coefficients of divergence of one distribution from another. Journal of the Royal Statistical Society: Series B, 28(1), 131–142. doi:10.1111/j.2517-6161.1966.tb00626.x.
3. Amari, S. (1985). Differential-Geometrical Methods in Statistics. Lecture Notes in Statistics 28, Springer, New York. doi:10.1007/978-1-4612-5056-2.
4. Attias, H. (1999). Inferring parameters and structure of latent variable models by variational Bayes. Proceedings of the 15th Conference on Uncertainty in Artificial Intelligence, 21–30. arXiv:1301.6676.
5. Bishop, C. M., Spiegelhalter, D. and Winn, J. (2003). VIBES: a variational inference engine for Bayesian networks. Advances in Neural Information Processing Systems 15.
6. Bishop, C. M. and Svensén, M. (2003). Bayesian hierarchical mixtures of experts. Proceedings of the 19th Conference on Uncertainty in Artificial Intelligence, 57–64. arXiv:1212.2447.
7. Blei, D. M., Kucukelbir, A. and McAuliffe, J. D. (2017). Variational inference: a review for statisticians. Journal of the American Statistical Association, 112(518), 859–877. doi:10.1080/01621459.2017.1285773. (A modern introduction to the material of Sections 10.1–10.4.)
8. Corduneanu, A. and Bishop, C. M. (2001). Variational Bayesian model selection for mixture distributions. Proceedings of the 8th International Conference on Artificial Intelligence and Statistics.
9. Gibbs, M. N. (1997). Bayesian Gaussian Processes for Regression and Classification. PhD thesis, University of Cambridge.
10. Jaakkola, T. S. and Jordan, M. I. (2000). Bayesian parameter estimation via variational methods. Statistics and Computing, 10(1), 25–37. doi:10.1023/A:1008932416310.
11. Jordan, M. I., Ghahramani, Z., Jaakkola, T. S. and Saul, L. K. (1999). An introduction to variational methods for graphical models. Machine Learning, 37(2), 183–233. doi:10.1023/A:1007665907178.
12. Kuss, M. and Rasmussen, C. E. (2005). Assessing approximate inference for binary Gaussian process classification. Journal of Machine Learning Research, 6, 1679–1704. (Cited by Bishop as Kuss and Rasmussen 2006.)
13. MacKay, D. J. C. (2003). Information Theory, Inference, and Learning Algorithms. Cambridge University Press.
14. Minka, T. P. (2001a). A Family of Algorithms for Approximate Bayesian Inference. PhD thesis, Massachusetts Institute of Technology.
15. Minka, T. P. (2001b). Expectation propagation for approximate Bayesian inference. Proceedings of the 17th Conference on Uncertainty in Artificial Intelligence, 362–369. arXiv:1301.2294.
16. Minka, T. P. (2004). Power EP. Technical report MSR-TR-2004-149, Microsoft Research.
17. Minka, T. P. (2005). Divergence measures and message passing. Technical report MSR-TR-2005-173, Microsoft Research.
18. Rockafellar, R. T. (1970). Convex Analysis. Princeton University Press.
19. Svensén, M. and Bishop, C. M. (2005). Robust Bayesian mixture modelling. Neurocomputing, 64, 235–252. doi:10.1016/j.neucom.2004.11.018. (Cited by Bishop as Svensén and Bishop 2004.)
20. Svensén, M. and Bishop, C. M. (2009). Pattern Recognition and Machine Learning: Solutions to the Exercises (web edition). Microsoft Research.
21. Wainwright, M. J., Jaakkola, T. S. and Willsky, A. S. (2005). A new class of upper bounds on the log partition function. IEEE Transactions on Information Theory, 51(7), 2313–2335. doi:10.1109/TIT.2005.850091.
22. Wainwright, M. J. and Jordan, M. I. (2008). Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, 1(1–2), 1–305. doi:10.1561/2200000001. (Further reading on Sections 10.4 and 10.7.)
23. Wiegerinck, W. and Heskes, T. (2003). Fractional belief propagation. Advances in Neural Information Processing Systems 15.
24. Winn, J. and Bishop, C. M. (2005). Variational message passing. Journal of Machine Learning Research, 6, 661–694.

Numbers quoted in these notes (the table after Result 10.2, the integral example after Result 10.23, and the numerical checks mentioned in the derivations) were computed for these notes in Python with numpy and scipy. Figures are reproduced by the course code, not copied from the book.