Mixture Models and EM

A mixture model explains each observation by a hidden label saying which simple component produced it. This chapter turns that idea into two algorithms, K-means and expectation-maximisation (EM), proves why EM can only climb the likelihood, and applies it to Gaussian and Bernoulli mixtures and to Bayesian regression. Every derivation the book leaves to the reader is written out below.

TrungTin Nguyen · Following Bishop (2006), Pattern Recognition and Machine Learning, Chapter 9 · 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

9.1 K-means clustering

pp. 424–428 We have \(N\) points \(\mathbf{x}_1,\dots,\mathbf{x}_N\in\R^D\) and want to split them into \(K\) groups, with \(K\) fixed in advance. Each cluster gets a prototype \(\bm\mu_k\in\R^D\) (it will turn out to be the cluster's centre), and each point gets a 1-of-K indicator vector \((r_{n1},\dots,r_{nK})\): \(r_{nk}\in\{0,1\}\) and \(\sum_k r_{nk}=1\), with \(r_{nk}=1\) meaning "point \(n\) is in cluster \(k\)".

Minimising \(J\) jointly is a hard combinatorial problem: finding the global minimum is NP-hard even for \(K=2\) (Aloise et al., 2009). K-means (Lloyd, 1982) settles for block coordinate descent: starting from some prototypes, it alternately minimises \(J\) over the assignments with the prototypes held fixed, and over the prototypes with the assignments held fixed. Both minimisations have closed forms (9.2)–(9.4). Bishop calls them the E step and the M step, because Section 9.3.2 shows that they are a limiting case of the E and M steps of EM.

\[ r_{nk}=\begin{cases}1 & \text{if } k=\argmin_j \lVert\mathbf{x}_n-\bm\mu_j\rVert^2,\\ 0 & \text{otherwise,}\end{cases} \tag{9.2} \] \[ 2\sum_{n=1}^{N} r_{nk}\,(\mathbf{x}_n-\bm\mu_k)=\mathbf{0}, \tag{9.3} \] \[ \bm\mu_k=\frac{\sum_n r_{nk}\,\mathbf{x}_n}{\sum_n r_{nk}} . \tag{9.4} \]

Proof

(a) Each step minimises \(J\) exactly over its own block of variables (gap box above); the current values of that block are feasible, so the minimum cannot exceed the current value of \(J\).

(b) For an assignment \(r\) let \(F(r)=\min_{\bm\mu}J(r,\bm\mu)\), the within-cluster scatter of the partition (empty clusters contribute 0). Number the iterations so that iteration \(t\) produces \(r^{(t)}\) (E step) and then \(\bm\mu^{(t)}\) (M step). After the M step, \(J(r^{(t)},\bm\mu^{(t)})=F(r^{(t)})\) whatever convention is used for empty clusters. Suppose \(r^{(t+1)}\ne r^{(t)}\). By the tie rule a point changes cluster only if its new prototype is strictly nearer, so the E step strictly lowers \(J\):

\[ F(r^{(t+1)})=J(r^{(t+1)},\bm\mu^{(t+1)})\le J(r^{(t+1)},\bm\mu^{(t)})\lt J(r^{(t)},\bm\mu^{(t)})=F(r^{(t)}). \]

So \(F\) strictly decreases along the sequence of assignments for as long as it changes, and an assignment cannot reappear (it would have to have a smaller \(F\) than itself). There are at most \(K^N\) assignments, so after at most \(K^N\) changes the E step reproduces the previous assignment. The M step then reproduces the same means (and empty clusters keep their prototypes), so nothing changes in later iterations either. (The bound \(K^N\) is very crude; in practice a few dozen iterations suffice.)

(c) The first statement says that each block is optimal given the other, which is what the stopping rule checks. For the local minimum, let \(\bm\mu^\star\) be the final prototypes and \(r^\star\) the final assignment. If every point has a unique nearest prototype, the strict inequalities \(\lVert\mathbf{x}_n-\bm\mu^\star_{k(n)}\rVert\)\({}\lt\lVert\mathbf{x}_n-\bm\mu^\star_j\rVert\) (\(j\ne k(n)\)) persist under small perturbations, because distances are continuous: there is \(\delta\gt0\) such that every \(\bm\mu\) with \(\max_k\lVert\bm\mu_k-\bm\mu^\star_k\rVert\)\({}\lt\delta\) gives the same nearest prototypes. On that neighbourhood \(J^\star(\bm\mu)\)\({}=J(r^\star,\bm\mu)\)\({}=F(r^\star)+\sum_kN_k\lVert\bm\mu_k-\bm\mu^\star_k\rVert^2\) by the identity of the gap box, which is larger than \(F(r^\star)=J^\star(\bm\mu^\star)\) unless \(\bm\mu=\bm\mu^\star\).

Not global. Take the four points \((\pm1,\pm a)\) with \(0\lt a\lt1\) and \(K=2\). Splitting them into top and bottom, with prototypes \((0,\pm a)\), is a fixed point: each point is at distance 1 from its own prototype and \(\sqrt{1+4a^2}\gt1\) from the other, and each prototype is the mean of its two points. Its distortion is \(J=4\). The left/right split, with prototypes \((\pm1,0)\), has \(J=4a^2\lt4\). Started near the first configuration, K-means stays there. ∎

Practical remarks p. 427. Because of Result 9.1(c), K-means is usually run from several starts and the run with the smallest \(J\) is kept. Starting from \(K\) data points chosen at random is the simplest choice; spreading the initial points out (as in the "k-means++" seeding of Arthur and Vassilvitskii, 2007) is usually better. Since \(J\) uses Euclidean distance, the result depends on the units of the variables, which is why the Old Faithful data are standardised first. K-means is also the standard way to initialise EM for Gaussian mixtures (Section 9.2.2). Each E step compares every point with every prototype; tree data structures and the triangle inequality can skip most of these distance computations (the book cites Moore, 2000, and Elkan, 2003).

Sequential (on-line) K-means

When data arrive one at a time, MacQueen's (1967) version updates only the prototype nearest to the new point \(\mathbf{x}_n\):

\[ \bm\mu_k^{\rm new}=\bm\mu_k^{\rm old}+\eta_n\,(\mathbf{x}_n-\bm\mu_k^{\rm old}), \tag{9.5} \]

where \(\eta_n\) is a learning rate that decreases as more data are seen.

K-medoids

Squared Euclidean distance is unsuitable for categorical variables and makes the means sensitive to outliers (one far-away point can drag a mean anywhere). Replacing it by a general dissimilarity \(\mathcal{V}(\mathbf{x},\mathbf{x}')\) gives the distortion

\[ \widetilde J=\sum_{n=1}^{N}\sum_{k=1}^{K} r_{nk}\,\mathcal{V}(\mathbf{x}_n,\bm\mu_k), \tag{9.6} \]

minimised by the K-medoids algorithm. Its E step is again "assign each point to the prototype with the smallest dissimilarity". In the M step the prototype of each cluster is usually restricted to be one of the cluster's own points (its medoid), so that only evaluations of \(\mathcal{V}\) are needed (Kaufman and Rousseeuw, 1990).

Hard versus soft assignments p. 428. K-means assigns every point to exactly one cluster, even a point that lies almost halfway between two centres. The probabilistic treatment of Section 9.2 replaces the 0/1 indicators by posterior probabilities, which express how uncertain each assignment is.

9.1.1 Image segmentation and compression

pp. 428–430 Treat every pixel of a colour image as a point \((R,G,B)\in[0,1]^3\) and run K-means on these \(N\) points, ignoring where the pixels are in the image. Redrawing each pixel in the colour of its cluster centre represents the image with a palette of \(K\) colours (Figure 9.3). As a segmentation method this is crude, because neighbouring pixels are not encouraged to share a cluster; it simply shows what K-means does. (The colour cube is bounded rather than all of \(\R^3\), but nothing in the algorithm requires an unbounded space: means of colours are colours.)

Vector quantisation. The same output is a lossy compression scheme. Store, for every point, only the label \(k\) of its cluster, plus the \(K\) code-book vectors \(\bm\mu_k\); decode each point as its code-book vector. A new point is compressed by finding its nearest code-book vector. The mean squared reconstruction error is \(J/N\), so K-means is precisely the algorithm that designs a code book of size \(K\) to minimise this error for the given data; in quantisation theory it is Lloyd's algorithm (Lloyd, 1982; Gersho and Gray, 1992).

9.2 Mixtures of Gaussians

pp. 430–432 Section 2.3.9 introduced the Gaussian mixture (2.188) as a weighted sum of Gaussian densities,

\[ p(\mathbf{x})=\sum_{k=1}^{K}\pi_k\,\N(\mathbf{x}\given\bm\mu_k,\bm\Sigma_k). \tag{9.7} \]

The same distribution arises from a two-stage story: first pick a component at random, then draw \(\mathbf{x}\) from that component. Making the hidden choice an explicit random variable is what leads to EM.

Proof

Everything rests on one observation: if \(\mathbf{z}=\mathbf{e}_j\) then \(z_k=1\) for \(k=j\) and \(z_k=0\) otherwise, so for any numbers \(a_1,\dots,a_K\), \(\prod_ka_k^{z_k}=a_j\) (every factor with \(z_k=0\) equals \(a_k^0=1\)).

(a) Multiply (9.10) and (9.11) and combine the two products over \(k\): \(\prod_k\pi_k^{z_k}\prod_k\N_k^{z_k}\)\({}=\prod_k(\pi_k\N_k)^{z_k}\), writing \(\N_k\)\({}=\N(\mathbf{x}\given\bm\mu_k,\bm\Sigma_k)\).

(b) By the sum rule, \(p(\mathbf{x})\) is the sum of (a) over the \(K\) possible values \(\mathbf{z}=\mathbf{e}_j\), and by the observation above the term for \(\mathbf{e}_j\) is \(\pi_j\N_j\). Summing over \(j\) gives (9.12).

(c) Bayes' theorem: \(p(z_k=1\given\mathbf{x})\)\({}=p(\mathbf{x},\mathbf{z}=\mathbf{e}_k)/p(\mathbf{x})\)\({}=\pi_k\N_k/\sum_j\pi_j\N_j\). ∎

Reading the responsibilities. \(\pi_k\) is the prior probability of component \(k\) and \(\gamma(z_k)\) its posterior probability after seeing \(\mathbf{x}\). If all \(\pi_k\gt0\) then \(0\lt\gamma(z_k)\lt1\) (Gaussian densities are positive) and \(\sum_k\gamma(z_k)=1\). Writing \(a_k\)\({}=\ln\pi_k+\ln\N(\mathbf{x}\given\bm\mu_k,\bm\Sigma_k)\), (9.13) is the softmax \(\gamma(z_k)=e^{a_k}/\sum_je^{a_j}\), which is how it should be computed (log-sum-exp). For \(K=2\) it is the logistic sigmoid of the log posterior odds, \(\gamma(z_1)=\sigma(a_1-a_2)\), as for the generative classifiers of Section 4.2; and when all components share one covariance matrix, the quadratic terms cancel in \(a_k-a_j\), so the sets where two responsibilities are equal are hyperplanes.

Sampling. p. 432 The two-stage story gives a direct way to simulate from the mixture, an instance of ancestral sampling (Section 8.1.2): sample the parent \(\mathbf{z}\), then the child \(\mathbf{x}\).

In the language of Section 9.3, the pairs \(\{(\mathbf{x}_n,\mathbf{z}_n)\}\) form a complete data set and the \(\mathbf{x}_n\) alone an incomplete one. Figure 9.5 shows both, and also the responsibilities, which are the best we can say about the labels once they have been thrown away.

9.2.1 Maximum likelihood

pp. 432–435 Collect the observations in the \(N\times D\) matrix \(\mathbf{X}\) with rows \(\mathbf{x}_n^\T\), and the latent vectors in the \(N\times K\) matrix \(\mathbf{Z}\) with rows \(\mathbf{z}_n^\T\). If the pairs \((\mathbf{x}_n,\mathbf{z}_n)\) are independent and identically distributed (the plate diagram of Figure 9.6), the log likelihood of the observed data is, from (9.7),

\[ \ln p(\mathbf{X}\given\bm\pi,\bm\mu,\bm\Sigma)=\sum_{n=1}^{N}\ln\Big\{\sum_{k=1}^{K}\pi_k\,\N(\mathbf{x}_n\given\bm\mu_k,\bm\Sigma_k)\Big\}. \tag{9.14} \]

Before maximising (9.14) we must face two features that a single Gaussian does not have: the likelihood can be made arbitrarily large, and every solution comes in \(K!\) copies.

Proof

The singular term. A component centred on a data point contributes, at that point, the density

\[ \N(\mathbf{x}_n\given\mathbf{x}_n,\sigma_j^2\mathbf{I})=\frac{1}{(2\pi)^{D/2}}\,\frac{1}{\sigma_j^{D}}. \tag{9.15} \]

Erratum. The book prints \(\frac{1}{(2\pi)^{1/2}}\frac{1}{\sigma_j}\), which is the case \(D=1\); the official errata replace \(\sigma_j\) by \(\sigma_j^D\), and for \(D\gt1\) the constant is also \((2\pi)^{D/2}\). The conclusion is the same: the term grows without bound as \(\sigma_j\to0\).

(a) Every term of (9.14) is the log of a sum of nonnegative numbers, so it is at least the log of any one of them. For \(n\ne n_0\) keep the term of component \(i\): \(\ln\sum_k\pi_k\N_{nk}\)\({}\ge\ln\{\pi_i\N(\mathbf{x}_n\given\bm\mu_i,\bm\Sigma_i)\}\)\({}=:c_n\), a finite number that does not depend on \(\sigma\). For \(n=n_0\) keep the term of component \(j\) and use (9.15). Adding up,

\[ \ln p(\mathbf{X}\given\bm\theta)\ \ge\ \sum_{n\ne n_0}c_n+\ln\pi_j-\frac D2\ln(2\pi)-D\ln\sigma\ \longrightarrow\ +\infty\quad(\sigma\to0). \]

(b) For fixed \(\sigma^2\), \(\ln p\)\({}=-\frac{ND}{2}\ln(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_n\lVert\mathbf{x}_n-\bm\mu\rVert^2\), and \(\sum_n\lVert\mathbf{x}_n-\bm\mu\rVert^2\)\({}=NDs^2+N\lVert\bar{\mathbf{x}}-\bm\mu\rVert^2\)\({}\ge NDs^2\) (the decomposition of the first gap box). So \(\ln p\)\({}\le g(\sigma^2)\)\({}=-\frac{ND}{2}\ln(2\pi\sigma^2)-\frac{NDs^2}{2\sigma^2}\). If \(s^2\gt0\), \(g'(\sigma^2)\)\({}=-\frac{ND}{2\sigma^2}+\frac{NDs^2}{2\sigma^4}\) vanishes only at \(\sigma^2=s^2\), where \(g\) is largest (it increases before and decreases after), with value \(-\frac{ND}{2}\{\ln(2\pi s^2)+1\}\).

Why the collapse fails here: with \(\bm\mu=\mathbf{x}_{n_0}\), \(\ln p\)\({}=-ND\ln\sigma-\frac{ND}{2}\ln(2\pi)-\frac{c}{\sigma^2}\) with \(c\)\({}=\frac12\sum_{n\ne n_0}\lVert\mathbf{x}_n-\mathbf{x}_{n_0}\rVert^2\)\({}\gt0\), and \(c/\sigma^2\) grows faster than \(ND\ln(1/\sigma)\). Each of the other points contributes a factor proportional to \(\sigma^{-D}e^{-\lVert\mathbf{x}_n-\mathbf{x}_{n_0}\rVert^2/(2\sigma^2)}\), which goes to zero faster than any power of \(\sigma\): this is the book's "multiplicative factors ... go to zero exponentially fast". In a mixture, component \(i\) keeps the density at those points bounded below, which is why (a) works.

Full covariance: \(\ln p\)\({}=-\frac N2\{D\ln2\pi+\ln|\bm\Sigma|+\Tr(\bm\Sigma^{-1}\mathbf{S}_{\bm\mu})\}\) with \(\mathbf{S}_{\bm\mu}\)\({}=\mathbf{S}+(\bar{\mathbf{x}}-\bm\mu)(\bar{\mathbf{x}}-\bm\mu)^\T\) (trace identities), maximised at \(\bm\mu=\bar{\mathbf{x}}\), \(\bm\Sigma=\mathbf{S}\) by the argument of the covariance gap box with all \(\gamma_{nk}=1\), giving the stated value. If \(\mathbf{S}\) is singular (for instance \(N\le D\)), choose a unit vector \(\mathbf{v}\) with \(\mathbf{S}\mathbf{v}=\mathbf{0}\) and \(\bm\Sigma\)\({}=\mathbf{I}+(\epsilon-1)\mathbf{v}\mathbf{v}^\T\), \(\bm\mu=\bar{\mathbf{x}}\): then \(\ln|\bm\Sigma|\)\({}=\ln\epsilon\)\({}\to-\infty\) while \(\Tr(\bm\Sigma^{-1}\mathbf{S})\)\({}=\Tr(\mathbf{S})\) stays fixed, so even one Gaussian has an unbounded likelihood. ∎

So maximum likelihood for Gaussian mixtures is, strictly speaking, ill-posed: the "best" fit puts a spike on a data point. Collapse onto lower-dimensional sets is also possible: with full covariance matrices a component can flatten onto any \(D\) or fewer data points, which always lie in a hyperplane (in two dimensions, onto the line through two points). In practice we look for a well-behaved local maximum instead, and guard against collapse: by detecting a shrinking component and resetting it (the book's suggestion), by a floor on the variances, by a prior (a MAP estimate, Section 9.3), or by the Bayesian treatment of Chapter 10, where the problem disappears. The phenomenon goes back at least to Kiefer and Wolfowitz (1956); bounding the ratio of the component variances makes the constrained likelihood bounded (Hathaway, 1985).

Proof

\(p(\mathbf{x}\given\bm\theta_\tau)\)\({}=\sum_k\pi_{\tau(k)}\N(\mathbf{x}\given\bm\mu_{\tau(k)},\bm\Sigma_{\tau(k)})\) is the same finite sum as (9.7) with its terms reordered, hence equal to it; the log likelihood (9.14) is a function of these densities only. If \(\bm\theta_\tau=\bm\theta_{\tau'}\), then the triple in position \(k\) is the same for both, \((\pi,\bm\mu,\bm\Sigma)_{\tau(k)}\)\({}=(\pi,\bm\mu,\bm\Sigma)_{\tau'(k)}\), and since the triples are pairwise different, \(\tau(k)=\tau'(k)\) for every \(k\). So \(\tau\mapsto\bm\theta_\tau\) is one-to-one and gives \(K!\) distinct parameter points. A local maximiser stays a local maximiser after relabelling, because relabelling is a bijection of the parameter space that preserves the likelihood and is continuous with a continuous inverse. ∎

Consequences. The likelihood surface of a mixture with \(K\ge2\) is multimodal by symmetry alone. The parameters are at best identifiable up to relabelling: for finite Gaussian mixtures with distinct components and positive weights, equality of the densities does imply equality of the parameter sets up to a permutation (Teicher, 1963; Yakowitz and Spragins, 1968). At degenerate points (two equal components, or some \(\pi_k=0\)) whole continua of parameters give the same density. For density estimation none of this matters; for interpreting fitted components one reports label-free summaries or fixes an ordering. In a Bayesian treatment the posterior has \(K!\) symmetric modes, which is why the variational lower bound used for comparing values of \(K\) in Chapter 10 gets a \(\ln K!\) correction (Section 10.2.4).

Why there is no closed form p. 435. For a single Gaussian the logarithm acts directly on an exponential and the likelihood equations are linear in the mean. In (9.14) the logarithm acts on a sum, and the stationarity conditions derived in the next section involve the responsibilities, which depend on all the parameters through the ratio (9.13). The result is a coupled system of transcendental equations. No general closed-form solution is known, and none could be a single formula, because the system typically has many solutions (at least \(K!\) maxima by Result 9.4, plus saddle points such as the one where all components coincide). One option is a general-purpose gradient method (as for the mixture density networks of Chapter 5); EM is the other.

9.2.2 EM for Gaussian mixtures

pp. 435–439 At an interior maximum of (9.14) (all \(\pi_k\gt0\), all \(\bm\Sigma_k\) positive definite) the derivatives with respect to every parameter vanish. Working them out gives three conditions in which the responsibilities appear by themselves. For the means:

\[ 0=\sum_{n=1}^{N}\underbrace{\frac{\pi_k\,\N(\mathbf{x}_n\given\bm\mu_k,\bm\Sigma_k)}{\sum_j\pi_j\,\N(\mathbf{x}_n\given\bm\mu_j,\bm\Sigma_j)}}_{\gamma(z_{nk})}\,\bm\Sigma_k^{-1}(\mathbf{x}_n-\bm\mu_k), \tag{9.16} \] \[ \bm\mu_k=\frac{1}{N_k}\sum_{n=1}^{N}\gamma(z_{nk})\,\mathbf{x}_n, \tag{9.17} \] \[ N_k=\sum_{n=1}^{N}\gamma(z_{nk}). \tag{9.18} \]

\(N_k\) is the effective number of points assigned to component \(k\), and (9.17) is a weighted mean in which point \(n\) has weight \(\gamma(z_{nk})\).

For the covariances the same reasoning, together with the maximum likelihood covariance of a single Gaussian (Section 2.3.4), gives

\[ \bm\Sigma_k=\frac{1}{N_k}\sum_{n=1}^{N}\gamma(z_{nk})\,(\mathbf{x}_n-\bm\mu_k)(\mathbf{x}_n-\bm\mu_k)^\T, \tag{9.19} \]

a weighted sample covariance with the effective count \(N_k\) in the denominator.

For the mixing coefficients the constraint (9.9) must be respected. Following Appendix E, maximise the Lagrangian

\[ \ln p(\mathbf{X}\given\bm\pi,\bm\mu,\bm\Sigma)+\lambda\Big(\sum_{k=1}^{K}\pi_k-1\Big), \tag{9.20} \]

whose derivative with respect to \(\pi_k\) gives

\[ 0=\sum_{n=1}^{N}\frac{\N(\mathbf{x}_n\given\bm\mu_k,\bm\Sigma_k)}{\sum_j\pi_j\,\N(\mathbf{x}_n\given\bm\mu_j,\bm\Sigma_j)}+\lambda, \tag{9.21} \]

and then, with \(\lambda=-N\),

\[ \pi_k=\frac{N_k}{N}: \tag{9.22} \]

the mixing coefficient of a component is its average responsibility.

Proof

Each step of the three gap boxes is an equivalence. The gradient with respect to \(\bm\mu_k\) vanishes if and only if (9.17) holds (multiplication by the invertible \(\bm\Sigma_k\) is reversible). Given that, the derivative with respect to \(\bm\Lambda_k\) vanishes if and only if \(\mathbf{G}_k=\mathbf{0}\), i.e. (9.19). The constrained stationarity condition for \(\bm\pi\) (gradient of the Lagrangian zero for some \(\lambda\)) holds if and only if (9.22). In all three, \(\gamma_{nk}\) is evaluated at the same \(\bm\theta\), which is exactly the statement that \(\bm\theta\) is mapped to itself. ∎

Two warnings. First, (9.17), (9.19) and (9.22) are not a closed-form solution: the responsibilities on the right depend on the unknowns. Second, fixed points need not be maxima: if all components are identical (\(\bm\mu_k=\bar{\mathbf{x}}\), \(\bm\Sigma_k=\mathbf{S}\), any \(\bm\pi\)), then every \(\gamma_{nk}=\pi_k\) and one EM cycle returns the same parameters; this symmetric point is usually a saddle (compare Exercise 9.13). What the equations do suggest is to iterate them, and Section 9.4 proves that the iteration never decreases the likelihood.

9.3 An alternative view of EM

pp. 439–441 Now let \(\mathbf{X}\) be all the observed data, \(\mathbf{Z}\) all the latent variables (row \(\mathbf{z}_n^\T\) for observation \(n\)) and \(\bm\theta\) all the parameters. The log likelihood is

\[ \ln p(\mathbf{X}\given\bm\theta)=\ln\Big\{\sum_{\mathbf{Z}}p(\mathbf{X},\mathbf{Z}\given\bm\theta)\Big\}, \tag{9.29} \]

with the sum replaced by an integral for continuous latent variables. The sum inside the logarithm is the whole difficulty: even when the joint distribution \(p(\mathbf{X},\mathbf{Z}\given\bm\theta)\) is a simple exponential-family distribution, its marginal usually is not, and the logarithm can no longer undo the exponential.

The M step is tractable because, in (9.30), the logarithm acts on the joint distribution. When the complete-data model is an exponential family this has a simple and general form.

Proof

Take the logarithm: \(\ln p(\mathbf{X},\mathbf{Z}\given\bm\eta)\)\({}=\ln h+\ln g(\bm\eta)+\bm\eta^\T\mathbf{u}\). Its expectation under \(p(\mathbf{Z}\given\mathbf{X},\bm\eta^{\rm old})\) is \(\E[\ln h]+\ln g(\bm\eta)+\bm\eta^\T\bar{\mathbf{u}}\), by linearity, and \(\E[\ln h]\) does not depend on \(\bm\eta\). By the exponential-family identities, \(-\nabla\ln g(\bm\eta)\)\({}=\E_{\bm\eta}[\mathbf{u}]\) and \(-\nabla\nabla\ln g(\bm\eta)\)\({}=\cov_{\bm\eta}[\mathbf{u}]\), which is positive semidefinite. So \(\nabla_{\bm\eta}Q\)\({}=\bar{\mathbf{u}}-\E_{\bm\eta}[\mathbf{u}]\), and \(Q\) is concave, which makes any solution of \(\E_{\bm\eta}[\mathbf{u}]\)\({}=\bar{\mathbf{u}}\) a global maximiser. ∎

For the Gaussian mixture, the complete-data log likelihood (9.36) below is, for each \(k\), linear in \(\sum_nz_{nk}\), \(\sum_nz_{nk}\mathbf{x}_n\) and \(\sum_nz_{nk}\mathbf{x}_n\mathbf{x}_n^\T\) (expand the quadratic form). Their expectations are \(N_k\), \(\sum_n\gamma_{nk}\mathbf{x}_n\) and \(\sum_n\gamma_{nk}\mathbf{x}_n\mathbf{x}_n^\T\), which are exactly the quantities from which (9.24)–(9.27) are computed. (The mixing weights satisfy the constraint (9.9), so this is a constrained version of the result, handled by the Lagrange argument.)

MAP estimation. p. 441 With a prior \(p(\bm\theta)\) the same algorithm finds a mode of the posterior \(p(\bm\theta\given\mathbf{X})\): the E step is unchanged and the M step maximises \(Q(\bm\theta,\bm\theta^{\rm old})+\ln p(\bm\theta)\). A suitable prior removes the singularities of Section 9.2.1.

Missing data. p. 441 The unobserved variables need not be labels: EM also handles data sets in which some values are missing, by treating the missing values as \(\mathbf{Z}\). This maximises the likelihood of the observed values, which is the right thing to do when the data are missing at random: whether a value is missing may depend on the observed values, but not on the missing value itself. A sensor that fails to report whenever the true value exceeds a threshold violates this, and ignoring the mechanism then biases the estimates (here downwards).

9.3.1 Gaussian mixtures revisited

pp. 441–443 Suppose that, besides \(\mathbf{X}\), we were told which component generated each point (the complete data of Figure 9.5(a), graph in Figure 9.9). From (9.10) and (9.11) the complete-data likelihood is

\[ p(\mathbf{X},\mathbf{Z}\given\bm\mu,\bm\Sigma,\bm\pi)=\prod_{n=1}^{N}\prod_{k=1}^{K}\pi_k^{z_{nk}}\,\N(\mathbf{x}_n\given\bm\mu_k,\bm\Sigma_k)^{z_{nk}}, \tag{9.35} \] \[ \ln p(\mathbf{X},\mathbf{Z}\given\bm\mu,\bm\Sigma,\bm\pi)=\sum_{n=1}^{N}\sum_{k=1}^{K}z_{nk}\big\{\ln\pi_k+\ln\N(\mathbf{x}_n\given\bm\mu_k,\bm\Sigma_k)\big\}. \tag{9.36} \]

Compared with (9.14), the sum over \(k\) and the logarithm have changed places, and the logarithm now acts directly on each Gaussian. The maximisation splits into \(K\) single-Gaussian fits and

\[ \pi_k=\frac1N\sum_{n=1}^{N}z_{nk}. \tag{9.37} \]

We do not observe \(\mathbf{Z}\), so we average (9.36) over its posterior. From (9.35) and Bayes' theorem,

\[ p(\mathbf{Z}\given\mathbf{X},\bm\mu,\bm\Sigma,\bm\pi)\propto\prod_{n=1}^{N}\prod_{k=1}^{K}\big[\pi_k\,\N(\mathbf{x}_n\given\bm\mu_k,\bm\Sigma_k)\big]^{z_{nk}}. \tag{9.38} \]

Proof

Factorisation. \(p(\mathbf{Z}\given\mathbf{X},\bm\theta)\)\({}=p(\mathbf{X},\mathbf{Z}\given\bm\theta)/p(\mathbf{X}\given\bm\theta)\). The numerator is the product over \(n\) of \(f_n(\mathbf{z}_n)\)\({}=\prod_k[\pi_k\N_{nk}]^{z_{nk}}\), and the denominator is \(\prod_np(\mathbf{x}_n\given\bm\theta)\)\({}=\prod_n\sum_{\mathbf{z}_n}f_n(\mathbf{z}_n)\) (independent observations). So \(p(\mathbf{Z}\given\mathbf{X},\bm\theta)\)\({}=\prod_n\{f_n(\mathbf{z}_n)/\sum_{\mathbf{z}'_n}f_n(\mathbf{z}'_n)\}\), a product of normalised distributions, one for each \(\mathbf{z}_n\). This is the general fact (9.75) of Section 9.4.

The same by d-separation. In Figure 9.6, treat \(\bm\pi\), \(\bm\mu\), \(\bm\Sigma\) as fixed (conditioned on). Every path from \(\mathbf{z}_n\) to \(\mathbf{z}_m\), \(m\ne n\), passes through one of these parameter nodes with both edges pointing away from it (tail to tail: \(\mathbf{z}_n\leftarrow\bm\pi\)\({}\to\mathbf{z}_m\), or \(\mathbf{x}_n\leftarrow\bm\mu\)\({}\to\mathbf{x}_m\) inside the longer path \(\mathbf{z}_n\)\({}\to\mathbf{x}_n\leftarrow\bm\mu\)\({}\to\mathbf{x}_m\leftarrow\mathbf{z}_m\)), and a tail-to-tail node that is conditioned on blocks the path. All paths are blocked, so \(\mathbf{z}_n\) and \(\mathbf{z}_m\) are independent given \(\mathbf{X}\) and the parameters (Section 8.2).

The expectation. By the factorisation, \(\E[z_{nk}]\) involves only the distribution of \(\mathbf{z}_n\): the factors for \(m\ne n\) sum to the same constant in the numerator and the denominator and cancel, which leaves the first expression in (9.39). The sums run over the \(K\) values \(\mathbf{z}_n=\mathbf{e}_j\); in the numerator only \(\mathbf{z}_n=\mathbf{e}_k\) has \(z_{nk}\ne0\), and its product equals \(\pi_k\N_{nk}\) (the 1-of-K observation in the proof of Result 9.2); the denominator is \(\sum_j\pi_j\N_{nj}\). (Early printings write the first line of (9.39) differently; the form above is the corrected one.) ∎

Since (9.36) is linear in the \(z_{nk}\), its expectation is obtained by replacing each \(z_{nk}\) by \(\E[z_{nk}]=\gamma(z_{nk})\):

\[ \E_{\mathbf{Z}}\big[\ln p(\mathbf{X},\mathbf{Z}\given\bm\mu,\bm\Sigma,\bm\pi)\big]=\sum_{n=1}^{N}\sum_{k=1}^{K}\gamma(z_{nk})\big\{\ln\pi_k+\ln\N(\mathbf{x}_n\given\bm\mu_k,\bm\Sigma_k)\big\}. \tag{9.40} \]

Here the responsibilities are computed from the old parameters (E step), while the parameters inside the braces are free. Maximising over them (M step) gives back the updates (9.17), (9.19), (9.22), but now in a stronger sense.

Proof

Write (9.40) as \(\sum_kN_k\ln\pi_k+\sum_kQ_k(\bm\mu_k,\bm\Lambda_k)\), where, with \(\bm\Lambda_k=\bm\Sigma_k^{-1}\) and up to constants,

\[ Q_k(\bm\mu,\bm\Lambda)=\frac{N_k}{2}\ln|\bm\Lambda|-\frac12\sum_n\gamma_{nk}(\mathbf{x}_n-\bm\mu)^\T\bm\Lambda(\mathbf{x}_n-\bm\mu). \]

The three groups of parameters can be maximised separately.

(i) Mixing coefficients. By the Gibbs-inequality argument of the previous gap box, with \(p_k=N_k/N\), the unique maximiser is \(\pi_k=N_k/N\).

(ii) Means (Exercise 9.8). Fix \(\bm\Lambda\) positive definite. The gradient of \(Q_k\) in \(\bm\mu\) is \(\bm\Lambda\sum_n\gamma_{nk}(\mathbf{x}_n-\bm\mu)\) and the Hessian is \(-N_k\bm\Lambda\), negative definite. So \(Q_k(\cdot,\bm\Lambda)\) is strictly concave with unique maximiser \(\bm\mu_k^{\rm new}\)\({}=N_k^{-1}\sum_n\gamma_{nk}\mathbf{x}_n\), and this maximiser does not depend on \(\bm\Lambda\). Hence \(Q_k(\bm\mu,\bm\Lambda)\)\({}\le Q_k(\bm\mu_k^{\rm new},\bm\Lambda)\) for every \(\bm\mu\) and \(\bm\Lambda\).

(iii) Precisions (Exercise 9.9). By the trace trick, \(f(\bm\Lambda)\)\({}=Q_k(\bm\mu_k^{\rm new},\bm\Lambda)\)\({}=\frac{N_k}{2}\ln|\bm\Lambda|-\frac12\Tr(\bm\Lambda\widetilde{\mathbf{S}}_k)\). It is strictly concave on positive definite matrices: along a line \(\bm\Lambda+t\mathbf{H}\) (\(\mathbf{H}\) symmetric, nonzero) the trace term is linear, and \(\frac{\d^2}{\d t^2}\ln|\bm\Lambda+t\mathbf{H}|\)\({}=-\Tr(\bm\Lambda^{-1}\mathbf{H}\bm\Lambda^{-1}\mathbf{H})\)\({}=-\lVert\bm\Lambda^{-1/2}\mathbf{H}\bm\Lambda^{-1/2}\rVert_F^2\)\({}\lt0\) at \(t=0\) (differentiate \(\d\ln|\bm\Lambda|\)\({}=\Tr(\bm\Lambda^{-1}\d\bm\Lambda)\) once more, using \(\d\bm\Lambda^{-1}\)\({}=-\bm\Lambda^{-1}\d\bm\Lambda\,\bm\Lambda^{-1}\)). Its differential vanishes, as in the covariance gap box, exactly at \(\bm\Lambda^{-1}\)\({}=\widetilde{\mathbf{S}}_k/N_k\), which is positive definite when \(\widetilde{\mathbf{S}}_k\) is nonsingular; a stationary point of a strictly concave function is its unique global maximiser. This is (9.25).

Unboundedness. If \(\widetilde{\mathbf{S}}_k\mathbf{v}\)\({}=\mathbf{0}\) for a unit vector \(\mathbf{v}\), take \(\bm\Lambda\)\({}=\mathbf{I}+t\mathbf{v}\mathbf{v}^\T\): then \(\ln|\bm\Lambda|=\ln(1+t)\to\infty\) while \(\Tr(\bm\Lambda\widetilde{\mathbf{S}}_k)\)\({}=\Tr(\widetilde{\mathbf{S}}_k)\) stays fixed.

Combining (i)–(iii): for every admissible parameter value, (9.40) is at most its value at the EM update, with equality only there. Contrast with Section 9.2.2, where the same formulas were only necessary conditions, because there the responsibilities moved with the parameters. ∎

9.3.2 Relation to K-means

pp. 443–444 K-means makes hard assignments, EM soft ones. To connect them, give every component the same covariance \(\epsilon\mathbf{I}\), with \(\epsilon\) a fixed constant (not re-estimated):

\[ p(\mathbf{x}\given\bm\mu_k,\bm\Sigma_k)=\frac{1}{(2\pi\epsilon)^{D/2}}\exp\Big\{-\frac{1}{2\epsilon}\lVert\mathbf{x}-\bm\mu_k\rVert^2\Big\}. \tag{9.41} \]

Erratum. The book prints \((2\pi\epsilon)^{M/2}\); the dimension of \(\mathbf{x}\) is \(D\). The responsibilities (9.13) become

\[ \gamma(z_{nk})=\frac{\pi_k\exp\{-\lVert\mathbf{x}_n-\bm\mu_k\rVert^2/2\epsilon\}}{\sum_j\pi_j\exp\{-\lVert\mathbf{x}_n-\bm\mu_j\rVert^2/2\epsilon\}}, \tag{9.42} \]

a softmax of the negative squared distances at "temperature" \(2\epsilon\). As \(\epsilon\to0\) they harden into the K-means assignments (9.2), and the expected complete-data log likelihood becomes, up to a scale factor and constants,

\[ \E_{\mathbf{Z}}\big[\ln p(\mathbf{X},\mathbf{Z}\given\bm\mu,\bm\Sigma,\bm\pi)\big]\ \to\ -\frac12\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}\,\lVert\mathbf{x}_n-\bm\mu_k\rVert^2+\const. \tag{9.43} \]

The next result says precisely in what sense (9.43) holds.

Proof

(a) Write \(d_k\)\({}=\lVert\mathbf{x}_n-\bm\mu_k\rVert^2\). Divide the numerator and the denominator of (9.42) by \(\pi_je^{-d_j/2\epsilon}\):

\[ \gamma(z_{nk})=\frac{(\pi_k/\pi_j)\,e^{-(d_k-d_j)/2\epsilon}}{1+\sum_{i\ne j}(\pi_i/\pi_j)\,e^{-(d_i-d_j)/2\epsilon}}. \]

For \(k\ne j\), \(d_k-d_j\ge\Delta_n\gt0\), so the numerator is at most \((\pi_k/\pi_j)e^{-\Delta_n/2\epsilon}\)\({}\to0\) while the denominator is at least 1; hence \(\gamma(z_{nk})\to0\), and \(1-\gamma(z_{nj})\)\({}=\sum_{k\ne j}\gamma(z_{nk})\)\({}\le\sum_{k\ne j}(\pi_k/\pi_j)e^{-\Delta_n/2\epsilon}\)\({}=\frac{1-\pi_j}{\pi_j}e^{-\Delta_n/2\epsilon}\). The positivity of \(\pi_j\) is essential: a nearest component with \(\pi_j=0\) gets no responsibility. Points close to a bisector (small \(\Delta_n\)) need a very small \(\epsilon\) before their assignment hardens. With ties, the tied terms have \(d_i-d_j=0\) and keep the weights \(\pi_i/\pi_j\), while all others vanish, giving the stated split.

(b) \(N_k\)\({}=\sum_n\gamma(z_{nk})\)\({}\to\sum_nr_{nk}\), and (9.17) is a continuous function of the responsibilities wherever \(N_k\gt0\), so it converges to (9.4); likewise \(N_k/N\) converges to the cluster fraction.

(c) From (9.40) and (9.41),

\[ \E_{\mathbf{Z}}[\ln p]=\sum_{n,k}\gamma(z_{nk})\ln\pi_k-\frac{ND}{2}\ln(2\pi\epsilon)-\frac{1}{2\epsilon}\sum_{n,k}\gamma(z_{nk})\lVert\mathbf{x}_n-\bm\mu'_k\rVert^2 , \]

using \(\sum_k\gamma(z_{nk})=1\). Multiply by \(\epsilon\). The first term is bounded (by \(N\max_k|\ln\pi_k|\)), so times \(\epsilon\) it tends to 0; the second tends to 0 because \(\epsilon\ln\epsilon\to0\); the third is \(-\frac12\sum\gamma(z_{nk})\lVert\mathbf{x}_n-\bm\mu'_k\rVert^2\)\({}\to-\frac12\sum r_{nk}\lVert\mathbf{x}_n-\bm\mu'_k\rVert^2\) by (a). The distances are bounded for \(\bm\mu'\) in a bounded set, which makes the convergence uniform there. Unscaled, \(\E_{\mathbf{Z}}[\ln p]\)\({}=-\frac{1}{2\epsilon}\{J(r,\bm\mu')+o(1)\}+(\text{terms free of }\bm\mu')\), so for small \(\epsilon\) maximising it over the means is minimising \(J\). ∎

In the limit the mixing coefficients are still re-estimated as cluster fractions, but they no longer influence the assignments (except to split ties). K-means estimates no covariances at all. A hard-assignment algorithm with general covariances, which maximises the likelihood of the complete data jointly over the labels and the parameters, is known as the classification EM algorithm (Celeux and Govaert, 1992); with Mahalanobis distances it is the "elliptical K-means" that the book attributes to Sung and Poggio (1994).

9.3.3 Mixtures of Bernoulli distributions

pp. 444–448 Mixtures also model discrete data. Let \(\mathbf{x}=(x_1,\dots,x_D)^\T\) be a vector of \(D\) binary variables, independent Bernoulli given \(\bm\mu=(\mu_1,\dots,\mu_D)^\T\):

\[ p(\mathbf{x}\given\bm\mu)=\prod_{i=1}^{D}\mu_i^{x_i}(1-\mu_i)^{(1-x_i)}, \tag{9.44} \] \[ \E[\mathbf{x}]=\bm\mu, \tag{9.45} \] \[ \cov[\mathbf{x}]=\diag\{\mu_i(1-\mu_i)\}. \tag{9.46} \]

(Each \(x_i\) is Bernoulli: \(\E[x_i]=\mu_i\), and \(x_i^2=x_i\) gives \(\var[x_i]=\mu_i-\mu_i^2\); independence makes the off-diagonal covariances zero.) A single such distribution cannot represent any correlation between the variables. A mixture of them, also called latent class analysis (Lazarsfeld and Henry, 1968), can:

\[ p(\mathbf{x}\given\bm\mu,\bm\pi)=\sum_{k=1}^{K}\pi_k\,p(\mathbf{x}\given\bm\mu_k), \tag{9.47} \] \[ p(\mathbf{x}\given\bm\mu_k)=\prod_{i=1}^{D}\mu_{ki}^{x_i}(1-\mu_{ki})^{(1-x_i)}, \tag{9.48} \]

with \(\bm\mu=\{\bm\mu_1,\dots,\bm\mu_K\}\) and \(\bm\pi=\{\pi_1,\dots,\pi_K\}\). Its moments are

\[ \E[\mathbf{x}]=\sum_{k=1}^{K}\pi_k\bm\mu_k, \tag{9.49} \] \[ \cov[\mathbf{x}]=\sum_{k=1}^{K}\pi_k\big\{\bm\Sigma_k+\bm\mu_k\bm\mu_k^\T\big\}-\E[\mathbf{x}]\E[\mathbf{x}]^\T, \tag{9.50} \]

with \(\bm\Sigma_k\)\({}=\diag\{\mu_{ki}(1-\mu_{ki})\}\).

Proof

Expectations under the mixture are \(\pi\)-weighted averages of expectations under the components: \(\E[f(\mathbf{x})]\)\({}=\sum_k\pi_k\E_k[f(\mathbf{x})]\) (a sum or an integral of \(f\) against \(\sum_k\pi_kp(\mathbf{x}\given k)\)). With \(f(\mathbf{x})=\mathbf{x}\) this is (9.49). With \(f(\mathbf{x})\)\({}=\mathbf{x}\mathbf{x}^\T\) and \(\E_k[\mathbf{x}\mathbf{x}^\T]\)\({}=\bm\Sigma_k+\bm\mu_k\bm\mu_k^\T\) (the definition of a covariance), \(\E[\mathbf{x}\mathbf{x}^\T]\)\({}=\sum_k\pi_k(\bm\Sigma_k+\bm\mu_k\bm\mu_k^\T)\), and \(\cov[\mathbf{x}]\)\({}=\E[\mathbf{x}\mathbf{x}^\T]-\E[\mathbf{x}]\E[\mathbf{x}]^\T\) (1.42) gives (9.50). For the second form, expand \(\sum_k\pi_k(\bm\mu_k-\bar{\bm\mu})(\bm\mu_k-\bar{\bm\mu})^\T\)\({}=\sum_k\pi_k\bm\mu_k\bm\mu_k^\T-\bar{\bm\mu}\bar{\bm\mu}^\T-\bar{\bm\mu}\bar{\bm\mu}^\T+\bar{\bm\mu}\bar{\bm\mu}^\T\), using \(\sum_k\pi_k=1\) and \(\sum_k\pi_k\bm\mu_k=\bar{\bm\mu}\). (This is the law of total covariance, \(\cov[\mathbf{x}]\)\({}=\E[\cov[\mathbf{x}\given k]]+\cov[\E[\mathbf{x}\given k]]\).) For Bernoulli components \(\bm\Sigma_k\) is diagonal, so the off-diagonal entries come from the between-components term only. Example: \(K=2\), \(\bm\pi=(\frac12,\frac12)\), \(\bm\mu_1=(0.9,0.9)\), \(\bm\mu_2=(0.1,0.1)\): \(\bar{\bm\mu}=(0.5,0.5)\) and \(\cov[x_1,x_2]\)\({}=\frac12(0.4)(0.4)+\frac12(-0.4)(-0.4)\)\({}=0.16\), a correlation of \(0.16/0.25=0.64\) between two variables that are independent within each component. ∎

Given data \(\mathbf{X}\)\({}=\{\mathbf{x}_1,\dots,\mathbf{x}_N\}\), the log likelihood again has a sum inside the logarithm,

\[ \ln p(\mathbf{X}\given\bm\mu,\bm\pi)=\sum_{n=1}^{N}\ln\Big\{\sum_{k=1}^{K}\pi_k\,p(\mathbf{x}_n\given\bm\mu_k)\Big\}, \tag{9.51} \]

so we introduce a 1-of-K latent \(\mathbf{z}\) exactly as before:

\[ p(\mathbf{x}\given\mathbf{z},\bm\mu)=\prod_{k=1}^{K}p(\mathbf{x}\given\bm\mu_k)^{z_k}, \tag{9.52} \] \[ p(\mathbf{z}\given\bm\pi)=\prod_{k=1}^{K}\pi_k^{z_k}. \tag{9.53} \]

Summing their product over \(\mathbf{z}\) gives back (9.47) (Exercise 9.14: the argument of Result 9.2(b) word for word). The complete-data log likelihood, its expectation and the E step are

\[ \ln p(\mathbf{X},\mathbf{Z}\given\bm\mu,\bm\pi)=\sum_{n=1}^{N}\sum_{k=1}^{K}z_{nk}\Big\{\ln\pi_k+\sum_{i=1}^{D}\big[x_{ni}\ln\mu_{ki}+(1-x_{ni})\ln(1-\mu_{ki})\big]\Big\}, \tag{9.54} \] \[ \E_{\mathbf{Z}}\big[\ln p(\mathbf{X},\mathbf{Z}\given\bm\mu,\bm\pi)\big]=\sum_{n=1}^{N}\sum_{k=1}^{K}\gamma(z_{nk})\Big\{\ln\pi_k+\sum_{i=1}^{D}\big[x_{ni}\ln\mu_{ki}+(1-x_{ni})\ln(1-\mu_{ki})\big]\Big\}, \tag{9.55} \] \[ \gamma(z_{nk})=\E[z_{nk}]=\frac{\sum_{\mathbf{z}_n}z_{nk}\prod_{k'}[\pi_{k'}p(\mathbf{x}_n\given\bm\mu_{k'})]^{z_{nk'}}}{\sum_{\mathbf{z}_n}\prod_j[\pi_jp(\mathbf{x}_n\given\bm\mu_j)]^{z_{nj}}}=\frac{\pi_k\,p(\mathbf{x}_n\given\bm\mu_k)}{\sum_{j=1}^{K}\pi_j\,p(\mathbf{x}_n\given\bm\mu_j)}. \tag{9.56} \]

The responsibilities enter (9.55) only through

\[ N_k=\sum_{n=1}^{N}\gamma(z_{nk}) \tag{9.57} \] \[ \bar{\mathbf{x}}_k=\frac{1}{N_k}\sum_{n=1}^{N}\gamma(z_{nk})\,\mathbf{x}_n, \tag{9.58} \]

and the M step is

\[ \bm\mu_k=\bar{\mathbf{x}}_k, \tag{9.59} \] \[ \pi_k=\frac{N_k}{N}. \tag{9.60} \]

Proof

Each factor \(\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{1-x_{ni}}\) equals \(\mu_{ki}\) or \(1-\mu_{ki}\), a number in \([0,1]\), so \(0\)\({}\le p(\mathbf{x}_n\given\bm\mu_k)\)\({}\le1\). Then \(0\)\({}\le\sum_k\pi_kp(\mathbf{x}_n\given\bm\mu_k)\)\({}\le\sum_k\pi_k\)\({}=1\), each term of (9.51) is at most \(\ln1=0\), and so is their sum. For attainment: the likelihood \(\prod_n\sum_k\pi_kp(\mathbf{x}_n\given\bm\mu_k)\) is a polynomial in the parameters, hence continuous, and the parameter set is a product of closed intervals and a simplex, hence compact; a continuous function on a compact set attains its maximum. At \(\mu_{ki}=\frac12\), \(\pi_k=\frac1K\) the likelihood is \(2^{-ND}\gt0\), so the maximum is positive and its logarithm finite. ∎

The Gaussian case differs on both counts: there the densities are not bounded by 1, and the parameter set is not compact (a covariance can approach a singular matrix). Bernoulli mixtures do have points where the likelihood is zero (for instance \(\mu_{ki}=0\) for all \(k\) while some \(x_{ni}=1\)), but EM, which only moves uphill (Section 9.4), never reaches them from a start with positive likelihood. Maximum likelihood estimates on the boundary (\(\mu_{ki}=0\) or 1, a pixel that is always off in a cluster) are common, and they give probability zero to any new image with that pixel on; the Beta prior below fixes this.

Handwritten digits p. 447. Figure 9.10 fits \(K=3\) components to 600 binarised images of the digits 2, 3 and 4 by 10 EM iterations, without using the labels; each component's mean \(\bm\mu_k\), drawn as an image, looks like one of the three digits. The book initialises \(\pi_k=1/K\) and draws the \(\mu_{ki}\) uniformly in \((0.25,0.75)\), and adds that they were "normalized to satisfy the constraint that \(\sum_j\mu_{kj}=1\)". Bernoulli parameters satisfy no such constraint: each \(\mu_{ki}\) is a separate probability, and a constraint \(\sum_i\mu_{ki}=1\) belongs to the multinomial model of Exercise 9.19. The normalisation is harmless, though: it only changes the starting point (every \(\mu_{ki}\) becomes about \(1/784\), with random variation), and the first M step replaces each \(\bm\mu_k\) by the weighted mean \(\bar{\mathbf{x}}_k\). The reproduction below follows the book's recipe literally.

Priors and more than two states p. 448. Conjugate Beta priors on the \(\mu_{ki}\) act like extra pseudo-observations (Section 2.1.1), and EM then maximises the posterior (Exercise 9.18). Variables with \(M\gt2\) states are handled with the discrete distribution (2.26) in place of the Bernoulli, optionally with Dirichlet priors (Exercise 9.19).

9.3.4 EM for Bayesian linear regression

pp. 448–450 The latent variable can also be continuous. In the evidence framework of Section 3.5, the hyperparameters \(\alpha\) (prior precision of the weights) and \(\beta\) (noise precision) are chosen to maximise the marginal likelihood \(p(\mathbf{t}\given\alpha,\beta)\)\({}=\int p(\mathbf{t}\given\mathbf{w},\beta)\,p(\mathbf{w}\given\alpha)\,\d\mathbf{w}\) (3.77). Since \(\mathbf{w}\) is integrated out, it can play the role of \(\mathbf{Z}\). Recall the model: targets \(\mathbf{t}=(t_1,\dots,t_N)^\T\), design matrix \(\bm\Phi\) (\(N\times M\)) with rows \(\bm\phi_n^\T\)\({}=\bm\phi(\mathbf{x}_n)^\T\), likelihood \(p(\mathbf{t}\given\mathbf{w},\beta)\)\({}=\prod_n\N(t_n\given\mathbf{w}^\T\bm\phi_n,\beta^{-1})\) (3.10) and prior \(p(\mathbf{w}\given\alpha)\)\({}=\N(\mathbf{w}\given\mathbf{0},\alpha^{-1}\mathbf{I})\) (3.52). The E step needs the posterior \(\N(\mathbf{w}\given\mathbf{m}_N,\mathbf{S}_N)\) (3.49), (3.53)–(3.54); the complete-data log likelihood is

\[ \ln p(\mathbf{t},\mathbf{w}\given\alpha,\beta)=\ln p(\mathbf{t}\given\mathbf{w},\beta)+\ln p(\mathbf{w}\given\alpha), \tag{9.61} \]

and its posterior expectation is

\[ \E\big[\ln p(\mathbf{t},\mathbf{w}\given\alpha,\beta)\big]=\frac M2\ln\Big(\frac{\alpha}{2\pi}\Big)-\frac\alpha2\E\big[\mathbf{w}^\T\mathbf{w}\big]+\frac N2\ln\Big(\frac{\beta}{2\pi}\Big)-\frac\beta2\sum_{n=1}^{N}\E\big[(t_n-\mathbf{w}^\T\bm\phi_n)^2\big]. \tag{9.62} \]

Maximising (9.62) over \(\alpha\) gives (Exercise 9.20)

\[ \alpha=\frac{M}{\E[\mathbf{w}^\T\mathbf{w}]}=\frac{M}{\mathbf{m}_N^\T\mathbf{m}_N+\Tr(\mathbf{S}_N)}, \tag{9.63} \]

and an analogous update holds for \(\beta\) (Exercise 9.21). This differs in form from the update \(\alpha\)\({}=\gamma/\mathbf{m}_N^\T\mathbf{m}_N\) (3.92) obtained by differentiating the evidence directly, but both need an \(M\times M\) inverse (or eigendecomposition) per iteration. Their fixed points agree: with the eigenvalues \(\lambda_i\) of \(\beta\bm\Phi^\T\bm\Phi\) (3.87),

\[ \gamma=M-\alpha\sum_{i=1}^{M}\frac{1}{\lambda_i+\alpha}=M-\alpha\Tr(\mathbf{S}_N), \tag{9.64} \]

and at a stationary point of the evidence, where (3.92) holds,

\[ \alpha\,\mathbf{m}_N^\T\mathbf{m}_N=\gamma=M-\alpha\Tr(\mathbf{S}_N), \tag{9.65} \]

which, solved for \(\alpha\), is (9.63).

Proof

M step. \(\partial(9.62)/\partial\alpha\)\({}=\frac{M}{2\alpha}-\frac12\E[\mathbf{w}^\T\mathbf{w}]\), zero at (9.63); the second derivative \(-M/(2\alpha^2)\) is negative, so this is the maximum. Likewise \(\partial(9.62)/\partial\beta\)\({}=\frac{N}{2\beta}-\frac12\sum_n\E[(t_n-\mathbf{w}^\T\bm\phi_n)^2]\), zero at the stated \(\beta^{\rm new}\) by the previous gap box, with second derivative \(-N/(2\beta^2)\lt0\).

Two trace identities. Let \(\beta\bm\Phi^\T\bm\Phi\mathbf{u}_i\)\({}=\lambda_i\mathbf{u}_i\) with orthonormal \(\mathbf{u}_i\) (a symmetric positive semidefinite matrix). Then \(\mathbf{S}_N^{-1}\)\({}=\alpha\mathbf{I}+\beta\bm\Phi^\T\bm\Phi\) has the same eigenvectors and eigenvalues \(\alpha+\lambda_i\), so \(\Tr(\mathbf{S}_N)\)\({}=\sum_i\frac{1}{\alpha+\lambda_i}\) and

\[ \gamma=\sum_i\frac{\lambda_i}{\alpha+\lambda_i}=\sum_i\Big(1-\frac{\alpha}{\alpha+\lambda_i}\Big)=M-\alpha\Tr(\mathbf{S}_N),\qquad \Tr(\bm\Phi^\T\bm\Phi\mathbf{S}_N)=\frac1\beta\sum_i\frac{\lambda_i}{\alpha+\lambda_i}=\frac\gamma\beta . \]

The first is (9.64) (with (3.91)); the second holds because \(\beta\bm\Phi^\T\bm\Phi\mathbf{S}_N\) has eigenvalues \(\lambda_i/(\alpha+\lambda_i)\).

Fixed points. \(\alpha\)\({}=M/(\mathbf{m}_N^\T\mathbf{m}_N+\Tr\mathbf{S}_N)\)\({}\iff\alpha\mathbf{m}_N^\T\mathbf{m}_N\)\({}=M-\alpha\Tr\mathbf{S}_N\)\({}=\gamma\)\({}\iff\) (3.92); this is (9.65). And \(N/\beta\)\({}=\lVert\mathbf{t}-\bm\Phi\mathbf{m}_N\rVert^2+\gamma/\beta\)\({}\iff(N-\gamma)/\beta\)\({}=\lVert\mathbf{t}-\bm\Phi\mathbf{m}_N\rVert^2\)\({}\iff\) (3.95).

Stationarity of the evidence. Equations (3.92) and (3.95) are rearrangements of \(\partial\ln p(\mathbf{t}\given\alpha,\beta)/\partial\alpha\)\({}=0\) and \(\partial/\partial\beta=0\), (3.89) and (3.94). The book derives (3.89) "neglecting" the dependence of \(\mathbf{m}_N\) on \(\alpha\), but nothing is lost: \(E(\mathbf{m}_N)\)\({}=\min_{\mathbf{w}}E(\mathbf{w})\) with \(E(\mathbf{w})\)\({}=\frac\beta2\lVert\mathbf{t}-\bm\Phi\mathbf{w}\rVert^2+\frac\alpha2\mathbf{w}^\T\mathbf{w}\) (3.79), and since \(\nabla_{\mathbf{w}}E(\mathbf{m}_N)\)\({}=\mathbf{0}\), the chain rule gives \(\frac{\d}{\d\alpha}E(\mathbf{m}_N)\)\({}=\frac{\partial E}{\partial\alpha}(\mathbf{m}_N)\)\({}=\frac12\mathbf{m}_N^\T\mathbf{m}_N\) exactly (the envelope theorem), and similarly for \(\beta\). So the fixed points of EM are exactly the stationary points of the evidence, as Result 9.16 below guarantees in general. (A numerical check on random design matrices: the EM iterations converge to a point where both (3.92) and (3.95) hold and the finite-difference gradient of the log evidence vanishes.) ∎

The relevance vector machine pp. 449–450. In the RVM of Section 7.2.1 each weight has its own prior precision, \(p(\mathbf{w}\given\bm\alpha)\)\({}=\prod_i\N(w_i\given0,\alpha_i^{-1})\) (7.80), and the posterior is \(\N(\mathbf{w}\given\mathbf{m},\bm\Sigma)\) with \(\bm\Sigma\)\({}=(\mathbf{A}+\beta\bm\Phi^\T\bm\Phi)^{-1}\), \(\mathbf{A}=\diag(\alpha_i)\), \(\mathbf{m}\)\({}=\beta\bm\Sigma\bm\Phi^\T\mathbf{t}\) (7.81)–(7.83). The M step maximises

\[ \E_{\mathbf{w}}\big[\ln\{p(\mathbf{t}\given\mathbf{X},\mathbf{w},\beta)\,p(\mathbf{w}\given\bm\alpha)\}\big] \tag{9.66} \]

(the expectation under the posterior at the old hyperparameters; early printings lack the braces), giving

\[ \alpha_i^{\rm new}=\frac{1}{m_i^2+\Sigma_{ii}}, \tag{9.67} \] \[ (\beta^{\rm new})^{-1}=\frac{\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2+\beta^{-1}\sum_i\gamma_i}{N}, \tag{9.68} \]

with \(\gamma_i=1-\alpha_i\Sigma_{ii}\) (7.89). Erratum. Early printings have \(\mathbf{m}_N\) in (9.68); it is the RVM posterior mean \(\mathbf{m}\) of (7.82).

Proof

(a) As in the previous gap box, (9.66) equals, up to constants,

\[ \sum_i\Big[\frac12\ln\alpha_i-\frac{\alpha_i}{2}(m_i^2+\Sigma_{ii})\Big]+\frac N2\ln\beta-\frac\beta2\Big\{\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2+\Tr(\bm\Phi^\T\bm\Phi\bm\Sigma)\Big\}, \]

using \(\E[w_i^2]=m_i^2+\Sigma_{ii}\). Each \(\alpha_i\) appears in one concave term, \(\frac{1}{2\alpha_i}-\frac12(m_i^2+\Sigma_{ii})\)\({}=0\) gives (9.67), and the \(\beta\) derivative gives \((\beta^{\rm new})^{-1}\)\({}=\{\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2+\Tr(\bm\Phi^\T\bm\Phi\bm\Sigma)\}/N\). Finally, \(\beta\bm\Phi^\T\bm\Phi\)\({}=\bm\Sigma^{-1}-\mathbf{A}\) (old values), so \(\beta\bm\Phi^\T\bm\Phi\bm\Sigma\)\({}=\mathbf{I}-\mathbf{A}\bm\Sigma\) and \(\Tr(\bm\Phi^\T\bm\Phi\bm\Sigma)\)\({}=\beta^{-1}\sum_i(1-\alpha_i\Sigma_{ii})\)\({}=\beta^{-1}\sum_i\gamma_i\), which turns the \(\beta\) update into (9.68).

(b) At a fixed point old and new values coincide. Then \(\alpha_i(m_i^2+\Sigma_{ii})\)\({}=1\)\({}\iff\alpha_im_i^2\)\({}=1-\alpha_i\Sigma_{ii}\)\({}=\gamma_i\)\({}\iff\alpha_i\)\({}=\gamma_i/m_i^2\) (when \(m_i\ne0\)), which is (7.87); and \(N\beta^{-1}\)\({}=\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2+\beta^{-1}\sum_i\gamma_i\)\({}\iff\beta^{-1}(N-\sum_i\gamma_i)\)\({}=\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2\), which is (7.88). Away from fixed points the two iterations differ; the book reports (p. 348) that the direct updates converge somewhat faster (Tipping, 2001). ∎

9.4 The EM algorithm in general

pp. 450–455 We now prove that EM works, for any model with observed variables \(\mathbf{X}\), hidden variables \(\mathbf{Z}\) and a joint distribution \(p(\mathbf{X},\mathbf{Z}\given\bm\theta)\) (Dempster et al., 1977; the view below is that of Csiszár and Tusnády, 1984, Hathaway, 1986, and Neal and Hinton, 1998). The goal is to maximise

\[ p(\mathbf{X}\given\bm\theta)=\sum_{\mathbf{Z}}p(\mathbf{X},\mathbf{Z}\given\bm\theta), \tag{9.69} \]

assuming that \(p(\mathbf{X},\mathbf{Z}\given\bm\theta)\) is much easier to work with. For discrete \(\mathbf{Z}\) we write sums; for continuous or mixed \(\mathbf{Z}\) replace them by integrals. For any distribution \(q(\mathbf{Z})\) over the latent variables,

\[ \ln p(\mathbf{X}\given\bm\theta)=\LB(q,\bm\theta)+\KL(q\|p), \tag{9.70} \] \[ \LB(q,\bm\theta)=\sum_{\mathbf{Z}}q(\mathbf{Z})\ln\Big\{\frac{p(\mathbf{X},\mathbf{Z}\given\bm\theta)}{q(\mathbf{Z})}\Big\}, \tag{9.71} \] \[ \KL(q\|p)=-\sum_{\mathbf{Z}}q(\mathbf{Z})\ln\Big\{\frac{p(\mathbf{Z}\given\mathbf{X},\bm\theta)}{q(\mathbf{Z})}\Big\}. \tag{9.72} \]

\(\LB(q,\bm\theta)\) is a functional of the distribution \(q\) (Appendix D) and a function of \(\bm\theta\). Note the differences between (9.71) and (9.72): opposite signs, and the joint distribution in \(\LB\) against the posterior in \(\KL\).

Proof

Line 1 (product rule). For every \(\mathbf{Z}\) with \(p(\mathbf{X},\mathbf{Z}\given\bm\theta)\)\({}\gt0\), \(p(\mathbf{X},\mathbf{Z}\given\bm\theta)\)\({}=p(\mathbf{Z}\given\mathbf{X},\bm\theta)\,p(\mathbf{X}\given\bm\theta)\), so

\[ \ln p(\mathbf{X},\mathbf{Z}\given\bm\theta)=\ln p(\mathbf{Z}\given\mathbf{X},\bm\theta)+\ln p(\mathbf{X}\given\bm\theta). \tag{9.73} \]

Line 2 (divide by \(q\)). Let \(\mathbf{Z}\) have \(q(\mathbf{Z})\gt0\); by assumption its posterior, hence its joint probability, is positive. Subtracting \(\ln q(\mathbf{Z})\) from both sides of (9.73):

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

Line 3 (average over \(q\)). Multiply by \(q(\mathbf{Z})\) and sum over the \(\mathbf{Z}\) with \(q(\mathbf{Z})\gt0\) (the others contribute nothing to (9.71) or (9.72), with \(0\ln0=0\)). The left side becomes \(\LB(q,\bm\theta)\), the first term on the right becomes \(-\KL(q\|p)\), and the last becomes \(\ln p(\mathbf{X}\given\bm\theta)\sum_{\mathbf{Z}}q(\mathbf{Z})\)\({}=\ln p(\mathbf{X}\given\bm\theta)\), because \(q\) is normalised and \(\ln p(\mathbf{X}\given\bm\theta)\) does not depend on \(\mathbf{Z}\). Rearranging gives (9.70). All sums are finite sums of finite numbers.

Line 4 (sign of KL). By Gibbs' inequality (KL divergence, Section 1.6.1), using \(\ln t\le t-1\) with equality only at \(t=1\): \(-\KL(q\|p)\)\({}=\sum_{q\gt0}q\ln\frac pq\)\({}\le\sum_{q\gt0}q\big(\frac pq-1\big)\)\({}=\sum_{q\gt0}p-1\)\({}\le0\), with equality if and only if \(p=q\) wherever \(q\gt0\) and \(p\) puts no mass elsewhere, i.e. \(q=p\). Then \(\LB\)\({}=\ln p(\mathbf{X}\given\bm\theta)-\KL\)\({}\le\ln p(\mathbf{X}\given\bm\theta)\).

(If \(q\) puts mass where the posterior is zero, then \(\KL=+\infty\) and \(\LB=-\infty\), and (9.70) still holds in the extended sense. For mixtures with all \(\pi_k\gt0\) the posterior is positive everywhere and the condition is empty.) ∎

The two steps as coordinate ascent on \(\LB\). pp. 451–452 In the E step, \(\LB(q,\bm\theta^{\rm old})\) is maximised over \(q\) with \(\bm\theta^{\rm old}\) fixed. Since \(\ln p(\mathbf{X}\given\bm\theta^{\rm old})\) does not depend on \(q\), (9.70) says that maximising \(\LB\) is minimising the KL term, so the unique maximiser is the posterior \(q(\mathbf{Z})\)\({}=p(\mathbf{Z}\given\mathbf{X},\bm\theta^{\rm old})\), at which the bound equals the log likelihood (Result 9.14). In the M step, \(q\) is held fixed and \(\LB(q,\bm\theta)\) is maximised over \(\bm\theta\). Substituting the posterior into (9.71),

\[ \LB(q,\bm\theta)=\sum_{\mathbf{Z}}p(\mathbf{Z}\given\mathbf{X},\bm\theta^{\rm old})\ln p(\mathbf{X},\mathbf{Z}\given\bm\theta)-\sum_{\mathbf{Z}}p(\mathbf{Z}\given\mathbf{X},\bm\theta^{\rm old})\ln p(\mathbf{Z}\given\mathbf{X},\bm\theta^{\rm old})=Q(\bm\theta,\bm\theta^{\rm old})+\const, \tag{9.74} \]

so the M step of Section 9.3 is exactly the maximisation of the bound.

Proof

Apply (9.70) at \(\bm\theta^{\rm new}\) with this \(q\), then (9.74):

\[ \ln p(\mathbf{X}\given\bm\theta^{\rm new})=\LB(q,\bm\theta^{\rm new})+\KL\big(q\|p(\cdot\given\mathbf{X},\bm\theta^{\rm new})\big)=Q(\bm\theta^{\rm new},\bm\theta^{\rm old})+\mathrm{H}[q]+\KL\big(q\|p(\cdot\given\mathbf{X},\bm\theta^{\rm new})\big). \]

At \(\bm\theta^{\rm old}\) the KL term is zero (E step), so \(\ln p(\mathbf{X}\given\bm\theta^{\rm old})\)\({}=Q(\bm\theta^{\rm old},\bm\theta^{\rm old})+\mathrm{H}[q]\). Subtract. Both brackets are nonnegative (the first by the choice of \(\bm\theta^{\rm new}\), the second by Result 9.14), and the total is zero only if both are zero, i.e. \(Q(\bm\theta^{\rm new},\bm\theta^{\rm old})\)\({}=Q(\bm\theta^{\rm old},\bm\theta^{\rm old})\) and \(q\)\({}=p(\cdot\given\mathbf{X},\bm\theta^{\rm new})\). ∎

What monotonicity does and does not give. The sequence \(\ln p(\mathbf{X}\given\bm\theta^{(t)})\) is nondecreasing, so it converges whenever the likelihood is bounded above (Bernoulli mixtures; Gaussian mixtures with a variance floor or a Wishart prior), and otherwise it may run off to \(+\infty\) along a singularity (Figure 9.7): monotonicity drives EM into a collapse rather than protecting against it. It does not by itself imply that the parameters converge, or that their limit is a maximum. Under mild regularity conditions (for example, \(Q\) continuous in both arguments and the set \(\{\bm\theta:\ln p(\mathbf{X}\given\bm\theta)\)\({}\ge\ln p(\mathbf{X}\given\bm\theta^{(0)})\}\) compact), every limit point of an EM sequence is a stationary point of the likelihood (Wu, 1983; Boyles, 1983). That stationary point is usually a local maximum but can be a saddle (EM started exactly at a symmetric point such as the one in Exercise 9.13 stays there), and which local maximum is reached depends on the start.

Shape of the bound. p. 453 The book says the bound is "a convex function having a unique maximum"; erratum: it should read concave. For exponential-family components, \(Q(\bm\eta,\bm\eta^{\rm old})\) is concave in the natural parameters by Result 9.6 (for instance \(N_1\ln\pi+N_2\ln(1-\pi)\) for a mixing weight, or a downward parabola for a mean), which is what the figure shows. The bound also touches the log likelihood tangentially:

Proof

First proof (the KL term is minimal). By (9.70), \(\LB(q^{\rm old},\bm\theta)\)\({}=\ln p(\mathbf{X}\given\bm\theta)-\KL(q^{\rm old}\|p(\cdot\given\mathbf{X},\bm\theta))\). As a function of \(\bm\theta\), the KL term is nonnegative and equals 0 at \(\bm\theta^{\rm old}\), so \(\bm\theta^{\rm old}\) is a minimum of it, and its gradient there is zero. The gradients of \(\LB\) and \(\ln p\) therefore agree at \(\bm\theta^{\rm old}\).

Second proof (Fisher's identity). \(\nabla_{\bm\theta}\LB(q^{\rm old},\bm\theta)\)\({}=\sum_{\mathbf{Z}}q^{\rm old}(\mathbf{Z})\nabla_{\bm\theta}\ln p(\mathbf{X},\mathbf{Z}\given\bm\theta)\). At \(\bm\theta=\bm\theta^{\rm old}\), write \(q^{\rm old}\)\({}=p(\mathbf{X},\mathbf{Z}\given\bm\theta^{\rm old})/p(\mathbf{X}\given\bm\theta^{\rm old})\) and \(\nabla\ln p(\mathbf{X},\mathbf{Z}\given\bm\theta)\)\({}=\nabla p(\mathbf{X},\mathbf{Z}\given\bm\theta)/p(\mathbf{X},\mathbf{Z}\given\bm\theta)\):

\[ \nabla_{\bm\theta}\LB\big|_{\bm\theta^{\rm old}}=\frac{\sum_{\mathbf{Z}}\nabla_{\bm\theta}p(\mathbf{X},\mathbf{Z}\given\bm\theta^{\rm old})}{p(\mathbf{X}\given\bm\theta^{\rm old})}=\frac{\nabla_{\bm\theta}p(\mathbf{X}\given\bm\theta^{\rm old})}{p(\mathbf{X}\given\bm\theta^{\rm old})}=\nabla_{\bm\theta}\ln p(\mathbf{X}\given\bm\theta^{\rm old}), \]

i.e. the observed-data score is the posterior mean of the complete-data score.

Fixed points. If \(\bm\theta^\star\) maximises \(Q(\cdot,\bm\theta^\star)\) at an interior point, then \(\nabla_{\bm\theta}Q(\bm\theta,\bm\theta^\star)\)\({}=\mathbf{0}\) at \(\bm\theta^\star\); since \(\LB(q^\star,\bm\theta)\)\({}=Q(\bm\theta,\bm\theta^\star)+\const\), the identity gives \(\nabla\ln p(\mathbf{X}\given\bm\theta^\star)\)\({}=\mathbf{0}\). The converse reads the same identity backwards. (A numerical check on a two-component mixture with five free parameters: finite-difference gradients of \(\LB\) and of \(\ln p\) at \(\bm\theta^{\rm old}\) agree to about \(10^{-5}\).) ∎

Going further: how fast does EM converge?

The book notes that EM needs many more iterations than K-means. Near a maximum, EM converges linearly, at a rate set by the proportion of information that the latent variables hide (Dempster et al., 1977).

Proof

Jacobian of the EM map. Let \(\mathbf{g}(\bm\theta,\bm\theta')\)\({}=\nabla_{\bm\theta}Q(\bm\theta,\bm\theta')\); the M step solves \(\mathbf{g}(M(\bm\theta'),\bm\theta')\)\({}=\mathbf{0}\), and \(M(\bm\theta^\star)=\bm\theta^\star\) since EM converged there. Differentiating with respect to \(\bm\theta'\) at \(\bm\theta^\star\) (implicit function theorem): \(\partial_{\bm\theta}\mathbf{g}\,\mathbf{J}+\partial_{\bm\theta'}\mathbf{g}\)\({}=\mathbf{0}\), with \(\partial_{\bm\theta}\mathbf{g}\)\({}=-\mathbf{I}_{\rm c}\). So \(\mathbf{J}\)\({}=\mathbf{I}_{\rm c}^{-1}\,\partial_{\bm\theta'}\mathbf{g}\).

The cross derivative. \(g_i(\bm\theta,\bm\theta')\)\({}=\sum_{\mathbf{Z}}p(\mathbf{Z}\given\mathbf{X},\bm\theta')\,\partial_i\ln p(\mathbf{X},\mathbf{Z}\given\bm\theta)\). Differentiating in \(\theta'_j\) and setting \(\bm\theta\)\({}=\bm\theta'\)\({}=\bm\theta^\star\), with \(s_j\)\({}=\partial_j\ln p(\mathbf{Z}\given\mathbf{X},\bm\theta^\star)\),

\[ \frac{\partial g_i}{\partial\theta'_j}=\E\big[s_j\,\partial_i\ln p(\mathbf{X},\mathbf{Z}\given\bm\theta^\star)\big]=\E[s_js_i]+\partial_i\ln p(\mathbf{X}\given\bm\theta^\star)\,\E[s_j]=\E[s_is_j], \]

using \(\ln p(\mathbf{X},\mathbf{Z})\)\({}=\ln p(\mathbf{Z}\given\mathbf{X})+\ln p(\mathbf{X})\) and \(\E[s_j]\)\({}=\sum_{\mathbf{Z}}\partial_jp(\mathbf{Z}\given\mathbf{X},\bm\theta)\)\({}=\partial_j1\)\({}=0\). Differentiating \(\sum_{\mathbf{Z}}p(\mathbf{Z}\given\mathbf{X},\bm\theta)\,\partial_i\ln p(\mathbf{Z}\given\mathbf{X},\bm\theta)\)\({}=0\) once more gives the information identity \(\E[s_is_j]\)\({}=-\E[\partial_i\partial_j\ln p(\mathbf{Z}\given\mathbf{X},\bm\theta^\star)]\), so \(\partial_{\bm\theta'}\mathbf{g}\)\({}=\mathbf{I}_{\rm mis}\), a positive semidefinite matrix.

Missing information principle. For every \(\mathbf{Z}\), \(\ln p(\mathbf{X}\given\bm\theta)\)\({}=\ln p(\mathbf{X},\mathbf{Z}\given\bm\theta)-\ln p(\mathbf{Z}\given\mathbf{X},\bm\theta)\). Differentiate twice, average over \(p(\mathbf{Z}\given\mathbf{X},\bm\theta^\star)\) and evaluate at \(\bm\theta^\star\): \(-\mathbf{I}_{\rm obs}\)\({}=-\mathbf{I}_{\rm c}+\mathbf{I}_{\rm mis}\) (the expected Hessian of the complete-data log likelihood is \(\nabla^2_{\bm\theta}Q\) at \(\bm\theta^\star\)). Hence \(\mathbf{J}\)\({}=\mathbf{I}_{\rm c}^{-1}(\mathbf{I}_{\rm c}-\mathbf{I}_{\rm obs})\)\({}=\mathbf{I}-\mathbf{I}_{\rm c}^{-1}\mathbf{I}_{\rm obs}\).

Eigenvalues. \(\mathbf{J}\) is similar to \(\mathbf{I}_{\rm c}^{-1/2}\mathbf{I}_{\rm mis}\mathbf{I}_{\rm c}^{-1/2}\), which is symmetric positive semidefinite, so its eigenvalues are real and \(\ge0\); and \(\mathbf{I}-\mathbf{J}\) is similar to \(\mathbf{I}_{\rm c}^{-1/2}\mathbf{I}_{\rm obs}\mathbf{I}_{\rm c}^{-1/2}\), positive definite, so they are \(\lt1\). Finally \(\bm\theta^{(t+1)}-\bm\theta^\star\)\({}=M(\bm\theta^{(t)})-M(\bm\theta^\star)\)\({}=\mathbf{J}(\bm\theta^{(t)}-\bm\theta^\star)+o(\cdot)\) by differentiability of \(M\). ∎

A numerical illustration (seeded, our own run): 300 points, half from \(\N(0,1)\) and half from \(\N(s,1)\), fitted by a two-component mixture with unit variances and free \((\pi,\mu_1,\mu_2)\). The largest eigenvalue of \(\mathbf{J}\) is 0.25 for \(s=4\) (19 iterations to reach a tolerance of \(10^{-13}\)), 0.89 for \(s=2\) (224 iterations) and 0.97 for \(s=1.2\) (863 iterations); a finite-difference Jacobian of the actual EM map agrees with \(\mathbf{I}-\mathbf{I}_{\rm c}^{-1}\mathbf{I}_{\rm obs}\) to about \(3\times10^{-6}\). Newton-type methods converge faster near a maximum but give up EM's monotonicity and simplicity.

Independent observations. p. 453 For an i.i.d. data set, \(\mathbf{X}=\{\mathbf{x}_n\}\), \(\mathbf{Z}=\{\mathbf{z}_n\}\) and \(p(\mathbf{X},\mathbf{Z})\)\({}=\prod_np(\mathbf{x}_n,\mathbf{z}_n)\). The posterior then factorises,

\[ p(\mathbf{Z}\given\mathbf{X},\bm\theta)=\frac{p(\mathbf{X},\mathbf{Z}\given\bm\theta)}{\sum_{\mathbf{Z}}p(\mathbf{X},\mathbf{Z}\given\bm\theta)}=\frac{\prod_{n=1}^{N}p(\mathbf{x}_n,\mathbf{z}_n\given\bm\theta)}{\sum_{\mathbf{Z}}\prod_{n=1}^{N}p(\mathbf{x}_n,\mathbf{z}_n\given\bm\theta)}=\prod_{n=1}^{N}p(\mathbf{z}_n\given\mathbf{x}_n,\bm\theta), \tag{9.75} \]

so the responsibility of a component for \(\mathbf{x}_n\) depends only on \(\mathbf{x}_n\) and \(\bm\theta\), not on the other data points.

MAP estimation. p. 454 For a prior \(p(\bm\theta)\), \(p(\bm\theta\given\mathbf{X})\)\({}=p(\bm\theta,\mathbf{X})/p(\mathbf{X})\), so

\[ \ln p(\bm\theta\given\mathbf{X})=\ln p(\bm\theta,\mathbf{X})-\ln p(\mathbf{X}), \tag{9.76} \]

and with the decomposition (9.70),

\[ \begin{aligned}\ln p(\bm\theta\given\mathbf{X})&=\LB(q,\bm\theta)+\KL(q\|p)+\ln p(\bm\theta)-\ln p(\mathbf{X})\\&\ge\LB(q,\bm\theta)+\ln p(\bm\theta)-\ln p(\mathbf{X}),\end{aligned} \tag{9.77} \]

where \(\ln p(\mathbf{X})\) is a constant.

Generalised EM and ECM

p. 454 If the M step cannot be carried out exactly, Result 9.15 still applies to any \(\bm\theta^{\rm new}\) with \(Q(\bm\theta^{\rm new},\bm\theta^{\rm old})\)\({}\ge Q(\bm\theta^{\rm old},\bm\theta^{\rm old})\): increasing the bound is enough. This is the generalised EM (GEM) algorithm; a few iterations of a gradient or conjugate-gradient method on \(Q\) are one way to implement it. The expectation conditional maximisation (ECM) algorithm (Meng and Rubin, 1993) partitions \(\bm\theta\) into blocks and replaces the M step by a sequence of conditional maximisations, each over one block with the others held fixed; each increases \(Q\), so a full sequence is a GEM step. The Gaussian-mixture M step, means first and then covariances at the new means, (9.24)–(9.25), is such a pair of conditional maximisations; here it happens to reach the joint maximum, because the maximiser over the means does not depend on the covariances (Result 9.8).

Partial E steps

pp. 454–455 The E step can be generalised in the same way: any change of \(q\) that increases \(\LB(q,\bm\theta)\) is allowed (Neal and Hinton, 1998). EM is then plain coordinate ascent on the single function \(\LB(q,\bm\theta)\), and maxima of \(\LB\) correspond to maxima of the likelihood.

Proof

(a) Result 9.14: \(\LB(q,\bm\theta)\)\({}=\ln p(\mathbf{X}\given\bm\theta)-\KL(q\|q_{\bm\theta})\)\({}\le\ln p(\mathbf{X}\given\bm\theta)\), with equality exactly at \(q=q_{\bm\theta}\). So \(\sup_{q,\bm\theta}\LB\)\({}=\sup_{\bm\theta}\ln p(\mathbf{X}\given\bm\theta)\), and the global statements follow.

(b) First part: for \((q,\bm\theta)\) near \((q_{\bm\theta^\star},\bm\theta^\star)\), \(\LB(q,\bm\theta)\)\({}\le\ln p(\mathbf{X}\given\bm\theta)\)\({}\le\ln p(\mathbf{X}\given\bm\theta^\star)\)\({}=\LB(q_{\bm\theta^\star},\bm\theta^\star)\), the middle inequality because \(\bm\theta\) is near \(\bm\theta^\star\). Converse: \(q\mapsto\LB(q,\bm\theta^\star)\) is concave (\(\KL(q\|p)=\sum q\ln q-\sum q\ln p\) is convex in \(q\)), so a local maximiser over \(q\) is global, and (a) gives \(q^\star=q_{\bm\theta^\star}\). If \(\bm\theta^\star\) were not a local maximiser of the likelihood, there would be \(\bm\theta_j\to\bm\theta^\star\) with \(\ln p(\mathbf{X}\given\bm\theta_j)\)\({}\gt\ln p(\mathbf{X}\given\bm\theta^\star)\). Then \((q_{\bm\theta_j},\bm\theta_j)\)\({}\to(q^\star,\bm\theta^\star)\) by continuity, and \(\LB(q_{\bm\theta_j},\bm\theta_j)\)\({}=\ln p(\mathbf{X}\given\bm\theta_j)\)\({}\gt\ln p(\mathbf{X}\given\bm\theta^\star)\)\({}=\LB(q^\star,\bm\theta^\star)\), contradicting the local maximality of \((q^\star,\bm\theta^\star)\). This is the book's "by continuity". ∎

So any algorithm that keeps increasing \(\LB\), by full or partial steps in \(q\) or \(\bm\theta\), is hunting for maxima of the likelihood. Chapter 10 keeps this objective but restricts \(q\) to a tractable family, which is variational inference.

Incremental EM

p. 455 For i.i.d. data, take \(q=\prod_nq_n(\mathbf{z}_n)\), so that \(\LB\) is a sum over points (see the gap box on (9.75)). Updating a single \(q_m\) to the posterior of point \(m\) is a partial E step, and it increases \(\LB\). For exponential-family components the following M step only needs sufficient statistics, which change by the contribution of one point. For a Gaussian mixture, writing \(\gamma^{\rm old}(z_{mk})\) and \(\gamma^{\rm new}(z_{mk})\) for the responsibilities of point \(m\) before and after its partial E step,

\[ \bm\mu_k^{\rm new}=\bm\mu_k^{\rm old}+\Big(\frac{\gamma^{\rm new}(z_{mk})-\gamma^{\rm old}(z_{mk})}{N_k^{\rm new}}\Big)\big(\mathbf{x}_m-\bm\mu_k^{\rm old}\big), \tag{9.78} \] \[ N_k^{\rm new}=N_k^{\rm old}+\gamma^{\rm new}(z_{mk})-\gamma^{\rm old}(z_{mk}). \tag{9.79} \]

Proof

Write \(N=N_k^{\rm old}\), \(N'=N_k^{\rm new}\), \(\bm\mu=\bm\mu_k^{\rm old}\), \(\bm\mu'=\bm\mu_k^{\rm new}\), \(\mathbf{x}=\mathbf{x}_m\), \(\Delta=\Delta_k\) (for this proof only, \(N\) is the effective count of component \(k\)).

Counts (9.79). \(N_k=\sum_n\gamma(z_{nk})\) and only the term \(n=m\) changes, by \(\Delta\).

Means (9.78), Exercise 9.26. \(N'\bm\mu'\)\({}=\sum_n\gamma^{\rm new}(z_{nk})\mathbf{x}_n\)\({}=N\bm\mu+\Delta\mathbf{x}\). With \(N=N'-\Delta\): \(\bm\mu'\)\({}=\{(N'-\Delta)\bm\mu+\Delta\mathbf{x}\}/N'\)\({}=\bm\mu+\frac{\Delta}{N'}(\mathbf{x}-\bm\mu)\). Two consequences used below: \(\bm\mu'-\bm\mu\)\({}=\frac{\Delta}{N'}(\mathbf{x}-\bm\mu)\) and \(\mathbf{x}-\bm\mu'\)\({}=(1-\frac{\Delta}{N'})(\mathbf{x}-\bm\mu)\)\({}=\frac{N}{N'}(\mathbf{x}-\bm\mu)\).

Covariances, Exercise 9.27. Expanding (9.25) with \(\mathbf{W}_k\)\({}=\sum_n\gamma(z_{nk})\mathbf{x}_n\mathbf{x}_n^\T\) gives \(N_k\bm\Sigma_k\)\({}=\mathbf{W}_k-N_k\bm\mu_k\bm\mu_k^\T\), and \(\mathbf{W}'\)\({}=\mathbf{W}+\Delta\mathbf{x}\mathbf{x}^\T\). Using \(N'\bm\mu'=N\bm\mu+\Delta\mathbf{x}\),

\[ \begin{aligned} N'\bm\Sigma'-N\bm\Sigma&=\Delta\mathbf{x}\mathbf{x}^\T-(N\bm\mu+\Delta\mathbf{x})\bm\mu'^\T+N\bm\mu\bm\mu^\T=\Delta\,\mathbf{x}(\mathbf{x}-\bm\mu')^\T-N\bm\mu(\bm\mu'-\bm\mu)^\T\\ &=\frac{\Delta N}{N'}\,\mathbf{x}(\mathbf{x}-\bm\mu)^\T-\frac{\Delta N}{N'}\,\bm\mu(\mathbf{x}-\bm\mu)^\T=\frac{\Delta N}{N'}(\mathbf{x}-\bm\mu)(\mathbf{x}-\bm\mu)^\T. \end{aligned} \]

Divide by \(N'\) and write \(N/N'=1-\Delta/N'\): \(\bm\Sigma'\)\({}=\bm\Sigma+\frac{\Delta}{N'}\big[\frac{N}{N'}(\mathbf{x}-\bm\mu)(\mathbf{x}-\bm\mu)^\T-\bm\Sigma\big]\), and \(\frac{N}{N'}(\mathbf{x}-\bm\mu)\)\({}=\mathbf{x}-\bm\mu'\) gives the stated form. (This is the weighted version of Welford's running-covariance update; it agrees with direct recomputation to machine precision.)

Mixing coefficients. \(\pi_k=N_k/N_{\rm tot}\) with the total \(N_{\rm tot}=\sum_kN_k=N\) unchanged, because both the old and the new responsibilities of point \(m\) sum to one over \(k\), i.e. \(\sum_k\Delta_k=0\). So \(\pi_k\) changes by \(\Delta_k/N\) and still sums to one.

Cost. Each update is a few vector and outer-product operations per component, \(O(KD^2)\); \(\bm\Sigma_k\) changes by a scaling plus a rank-one term, so its Cholesky factor can be updated in \(O(D^2)\) rather than recomputed. ∎

Exercises

All 27 exercises of Chapter 9 (Bishop, 2006, pp. 455–459), restated in my own words, with complete solutions (hidden). Stars give the book's difficulty and "www" marks the exercises with official solutions. Several solutions are compact versions of the gap boxes and results above, which are linked for more detail.

1. Exercise 9.1 (★, www). Using that there are finitely many possible assignments \(\{r_{nk}\}\), and that each assignment has a unique optimal set of prototypes, show that K-means converges after a finite number of iterations.

   Solution

   Let \(F(r)=\min_{\bm\mu}J(r,\bm\mu)\); for nonempty clusters the minimiser is unique, the cluster means (gap box 9.1), and empty clusters contribute nothing to \(J\). After every M step, \(J=F(r^{(t)})\). Break ties in the E step in favour of the current cluster. If the E step changes the assignment, some point moves to a strictly nearer prototype, so \(J(r^{(t+1)},\bm\mu^{(t)})\)\({}\lt J(r^{(t)},\bm\mu^{(t)})\)\({}=F(r^{(t)})\), and the M step cannot increase \(J\): \(F(r^{(t+1)})\lt F(r^{(t)})\). A strictly decreasing sequence of values of \(F\) cannot revisit an assignment, and there are at most \(K^N\) assignments, so after finitely many iterations the E step leaves the assignment unchanged; the M step then reproduces the same means, and the algorithm has converged. The limit is a configuration in which neither step can improve \(J\), not necessarily the global minimum (Result 9.1).

   ▶ See it in a demo: a K-means run that stops at a poor local minimum

2. Exercise 9.2 (★). Apply the Robbins–Monro procedure of Section 2.3.5 to find the roots of the derivatives of \(J\) with respect to the \(\bm\mu_k\), and show that it gives the sequential update (9.5) of the nearest prototype.

   Solution

   For data arriving from \(p(\mathbf{x})\), the per-point expected distortion \(\bar J(\bm\mu)\)\({}=\E[\min_k\lVert\mathbf{x}-\bm\mu_k\rVert^2]\) has \(\nabla_{\bm\mu_k}\bar J\)\({}=\E[-2r_k(\mathbf{x};\bm\mu)(\mathbf{x}-\bm\mu_k)]\), where \(r_k(\mathbf{x};\bm\mu)=1\) if \(\bm\mu_k\) is the prototype nearest to \(\mathbf{x}\) (the integrand is differentiable except on cell boundaries, a null set when \(p\) has a density and the prototypes are distinct). This is a regression function whose root we want, and \(\mathbf{z}\)\({}=-2r_k(\mathbf{x}_n;\bm\mu)(\mathbf{x}_n-\bm\mu_k)\) is an unbiased noisy observation of it. The Robbins–Monro step (2.129), \(\bm\mu_k\leftarrow\bm\mu_k-a_{n-1}\mathbf{z}\), is \(\bm\mu_k\leftarrow\bm\mu_k+2a_{n-1}r_k(\mathbf{x}_n;\bm\mu)(\mathbf{x}_n-\bm\mu_k)\): nothing moves except the nearest prototype, which moves by (9.5) with \(\eta_n=2a_{n-1}\). The step sizes must satisfy (2.130)–(2.132); the choice \(\eta=1/m_k\), with \(m_k\) the number of points assigned to prototype \(k\) so far, does, and makes \(\bm\mu_k\) the running mean of those points (MacQueen, 1967). Details in gap box 9.2.

3. Exercise 9.3 (★, www). With \(p(\mathbf{z})\) given by (9.10) and \(p(\mathbf{x}\given\mathbf{z})\) by (9.11), show that summing \(p(\mathbf{z})p(\mathbf{x}\given\mathbf{z})\) over the values of \(\mathbf{z}\) gives the Gaussian mixture (9.7).

   Solution

   \(p(\mathbf{z})p(\mathbf{x}\given\mathbf{z})\)\({}=\prod_k\{\pi_k\N(\mathbf{x}\given\bm\mu_k,\bm\Sigma_k)\}^{z_k}\). The variable \(\mathbf{z}\) takes the \(K\) values \(\mathbf{e}_j\), and for \(\mathbf{z}=\mathbf{e}_j\) every factor with \(k\ne j\) is raised to the power 0, so the product equals \(\pi_j\N(\mathbf{x}\given\bm\mu_j,\bm\Sigma_j)\). Summing over \(j\) gives \(\sum_j\pi_j\N(\mathbf{x}\given\bm\mu_j,\bm\Sigma_j)\), which is (9.7).

4. Exercise 9.4 (★). To maximise the posterior \(p(\bm\theta\given\mathbf{X})\) of a latent-variable model by EM, show that the E step is as for maximum likelihood and that the M step maximises \(Q(\bm\theta,\bm\theta^{\rm old})+\ln p(\bm\theta)\).

   Solution

   \(\ln p(\bm\theta\given\mathbf{X})\)\({}=\ln p(\mathbf{X}\given\bm\theta)+\ln p(\bm\theta)-\ln p(\mathbf{X})\), and the last term is constant. With \(q(\mathbf{Z})\)\({}=p(\mathbf{Z}\given\mathbf{X},\bm\theta^{\rm old})\), Jensen's inequality gives \(\ln p(\mathbf{X}\given\bm\theta)\)\({}\ge Q(\bm\theta,\bm\theta^{\rm old})+\mathrm{H}[q]\) with equality at \(\bm\theta^{\rm old}\) (gap box 9.8). Adding \(\ln p(\bm\theta)\): the function \(Q(\bm\theta,\bm\theta^{\rm old})+\ln p(\bm\theta)+\mathrm{H}[q]\) is a lower bound on \(\ln p(\mathbf{X}\given\bm\theta)+\ln p(\bm\theta)\) that touches it at \(\bm\theta^{\rm old}\), so any \(\bm\theta\) that increases \(Q+\ln p(\bm\theta)\) increases the log posterior. The E step, the computation of \(q\), involves only the posterior of \(\mathbf{Z}\) at \(\bm\theta^{\rm old}\), in which the prior plays no role, so it is unchanged; the M step maximises \(Q(\bm\theta,\bm\theta^{\rm old})+\ln p(\bm\theta)\). (Equivalently, in the language of Section 9.4: maximise the right-hand side of (9.77) alternately over \(q\) and \(\bm\theta\); see gap box 9.17.)

5. Exercise 9.5 (★). Using d-separation in the graph of Figure 9.6, show that the posterior of the latent variables factorises over data points, \(p(\mathbf{Z}\given\mathbf{X},\bm\mu,\bm\Sigma,\bm\pi)\)\({}=\prod_{n=1}^{N}p(\mathbf{z}_n\given\mathbf{x}_n,\bm\mu,\bm\Sigma,\bm\pi)\) (9.80).

   Solution

   Condition on the observed \(\mathbf{x}_1,\dots,\mathbf{x}_N\) and on the parameters \(\bm\pi,\bm\mu,\bm\Sigma\). Any path from \(\mathbf{z}_n\) to \(\mathbf{z}_m\) (\(m\ne n\)) must leave the plate copy \(n\) through a parameter node: either \(\mathbf{z}_n\leftarrow\bm\pi\)\({}\to\mathbf{z}_m\), or \(\mathbf{z}_n\)\({}\to\mathbf{x}_n\leftarrow\bm\mu\)\({}\to\mathbf{x}_m\leftarrow\mathbf{z}_m\) (or the same with \(\bm\Sigma\)). In each case the parameter node is tail to tail and conditioned on, which blocks the path (the head-to-head nodes \(\mathbf{x}_n\), \(\mathbf{x}_m\) are observed and would let the path through, but one blocked node suffices). So all \(\mathbf{z}_n\) are mutually independent given \(\mathbf{X}\) and the parameters, and since \(\mathbf{z}_n\) is independent of \(\mathbf{x}_m\) (\(m\ne n\)) given \(\mathbf{x}_n\) and the parameters by the same argument, each factor is \(p(\mathbf{z}_n\given\mathbf{x}_n,\bm\theta)\). Algebraically: \(p(\mathbf{Z}\given\mathbf{X})\)\({}=\prod_np(\mathbf{x}_n,\mathbf{z}_n)/\prod_np(\mathbf{x}_n)\) (Result 9.7).

6. Exercise 9.6 (★★). Derive the EM algorithm for a Gaussian mixture in which all components share one covariance matrix \(\bm\Sigma\).

   Solution

   E step: (9.23) with every \(\bm\Sigma_k=\bm\Sigma\). M step: with \(\bm\Lambda=\bm\Sigma^{-1}\) the expected complete-data log likelihood is

   \[ \sum_{n,k}\gamma_{nk}\Big\{\ln\pi_k+\frac12\ln|\bm\Lambda|-\frac12(\mathbf{x}_n-\bm\mu_k)^\T\bm\Lambda(\mathbf{x}_n-\bm\mu_k)\Big\}+\const . \]

   For fixed \(\bm\Lambda\) it is strictly concave in each \(\bm\mu_k\) with maximiser \(\bm\mu_k\)\({}=N_k^{-1}\sum_n\gamma_{nk}\mathbf{x}_n\), independent of \(\bm\Lambda\) (as in Result 9.8); \(\pi_k=N_k/N\) as before. With these means, and \(\sum_kN_k=N\), the \(\bm\Lambda\)-terms are \(\frac N2\ln|\bm\Lambda|-\frac12\Tr\big(\bm\Lambda\sum_k\widetilde{\mathbf{S}}_k\big)\), with \(\widetilde{\mathbf{S}}_k\)\({}=\sum_n\gamma_{nk}(\mathbf{x}_n-\bm\mu_k)(\mathbf{x}_n-\bm\mu_k)^\T\). The differential calculation of gap box 9.6 gives

   \[ \bm\Sigma=\frac1N\sum_{k=1}^{K}\sum_{n=1}^{N}\gamma_{nk}(\mathbf{x}_n-\bm\mu_k)(\mathbf{x}_n-\bm\mu_k)^\T=\sum_{k=1}^{K}\pi_k\,\bm\Sigma_k^{\rm sep}, \]

   where \(\bm\Sigma_k^{\rm sep}\)\({}=\widetilde{\mathbf{S}}_k/N_k\) is what (9.25) would give for component \(k\) alone: the shared covariance is the \(\pi\)-weighted average of the separate ones (the pooled within-component covariance). Because \(\bm\Sigma\) is estimated from all \(N\) points, it is nonsingular as soon as the data span \(\R^D\) around the means, so a single component can no longer collapse onto one point.

7. Exercise 9.7 (★, www). Show that maximising the complete-data log likelihood (9.36) fits each component's mean and covariance to its own group of points, and sets the mixing coefficients to the fractions of points in the groups.

   Solution

   With \(\mathcal{C}_k=\{n:z_{nk}=1\}\) and \(N_k=|\mathcal{C}_k|\), (9.36) equals \(\sum_kN_k\ln\pi_k+\sum_k\sum_{n\in\mathcal{C}_k}\ln\N(\mathbf{x}_n\given\bm\mu_k,\bm\Sigma_k)\). Each \((\bm\mu_k,\bm\Sigma_k)\) appears only in the \(k\)th inner sum, the single-Gaussian log likelihood of group \(k\), maximised by the group mean and the group covariance with divisor \(N_k\) (2.121)–(2.122). For \(\bm\pi\), the Lagrangian \(\sum_kN_k\ln\pi_k+\lambda(\sum_k\pi_k-1)\) gives \(N_k/\pi_k+\lambda=0\); multiplying by \(\pi_k\) and summing over \(k\) gives \(\lambda=-N\) and \(\pi_k=N_k/N\), which is (9.37) and a global maximum by Gibbs' inequality (gap box 9.11).

8. Exercise 9.8 (★, www). Show that maximising (9.40) over \(\bm\mu_k\) with the responsibilities fixed gives (9.17).

   Solution

   By (2.43), the \(\bm\mu_k\)-dependent part of (9.40) is \(-\frac12\sum_n\gamma_{nk}(\mathbf{x}_n-\bm\mu_k)^\T\bm\Sigma_k^{-1}(\mathbf{x}_n-\bm\mu_k)\). Its gradient is \(\bm\Sigma_k^{-1}\sum_n\gamma_{nk}(\mathbf{x}_n-\bm\mu_k)\), which vanishes if and only if \(\bm\mu_k\)\({}=N_k^{-1}\sum_n\gamma_{nk}\mathbf{x}_n\), i.e. (9.17). The Hessian \(-N_k\bm\Sigma_k^{-1}\) is negative definite (for \(N_k\gt0\)), so this is the unique maximum.

9. Exercise 9.9 (★). Show that maximising (9.40) over \(\bm\Sigma_k\) and \(\pi_k\) with the responsibilities fixed gives (9.19) and (9.22).

   Solution

   With \(\bm\Lambda_k=\bm\Sigma_k^{-1}\) and the trace trick, the \(\bm\Lambda_k\)-terms of (9.40) are \(\frac{N_k}{2}\ln|\bm\Lambda_k|-\frac12\Tr(\bm\Lambda_k\widetilde{\mathbf{S}}_k)\) with \(\widetilde{\mathbf{S}}_k\)\({}=\sum_n\gamma_{nk}(\mathbf{x}_n-\bm\mu_k)(\mathbf{x}_n-\bm\mu_k)^\T\). For a symmetric perturbation, its differential is \(\frac12\Tr\{(N_k\bm\Lambda_k^{-1}-\widetilde{\mathbf{S}}_k)\d\bm\Lambda_k\}\); the matrix in braces is symmetric, so the differential vanishes for all symmetric \(\d\bm\Lambda_k\) if and only if \(\bm\Sigma_k\)\({}=\widetilde{\mathbf{S}}_k/N_k\), which is (9.19). The function is strictly concave in \(\bm\Lambda_k\), so this is the maximum (Result 9.8). For \(\pi_k\), the \(\bm\pi\)-terms are \(\sum_kN_k\ln\pi_k\); the Lagrange argument with \(\sum_k\pi_k=1\) gives \(\lambda=-N\) and \(\pi_k=N_k/N\), which is (9.22) (gap box 9.7).

10. Exercise 9.10 (★★). For a mixture \(p(\mathbf{x})\)\({}=\sum_k\pi_kp(\mathbf{x}\given k)\) (9.81) and a partition \(\mathbf{x}\)\({}=(\mathbf{x}_a,\mathbf{x}_b)\), show that \(p(\mathbf{x}_b\given\mathbf{x}_a)\) is again a mixture, and find its mixing coefficients and components.

    Solution

    The marginal is \(p(\mathbf{x}_a)\)\({}=\sum_j\pi_jp(\mathbf{x}_a\given j)\), with \(p(\mathbf{x}_a\given j)\) the marginal of component \(j\). Then

    \[ p(\mathbf{x}_b\given\mathbf{x}_a)=\frac{\sum_k\pi_kp(\mathbf{x}_a,\mathbf{x}_b\given k)}{p(\mathbf{x}_a)}=\sum_{k=1}^{K}\underbrace{\frac{\pi_k\,p(\mathbf{x}_a\given k)}{\sum_j\pi_j\,p(\mathbf{x}_a\given j)}}_{w_k(\mathbf{x}_a)}\ p(\mathbf{x}_b\given\mathbf{x}_a,k), \]

    writing \(p(\mathbf{x}_a,\mathbf{x}_b\given k)\)\({}=p(\mathbf{x}_a\given k)\,p(\mathbf{x}_b\given\mathbf{x}_a,k)\). The new mixing coefficients \(w_k(\mathbf{x}_a)\) are the responsibilities of the components for the observed part \(\mathbf{x}_a\); they are nonnegative and sum to one, and they depend on \(\mathbf{x}_a\). The components are the conditionals of the original components. For Gaussian components (partitioned Gaussians), \(p(\mathbf{x}_a\given k)\)\({}=\N(\mathbf{x}_a\given\bm\mu_{k,a},\bm\Sigma_{k,aa})\) and

    \[ p(\mathbf{x}_b\given\mathbf{x}_a,k)=\N\big(\mathbf{x}_b\given\bm\mu_{k,b}+\bm\Sigma_{k,ba}\bm\Sigma_{k,aa}^{-1}(\mathbf{x}_a-\bm\mu_{k,a}),\ \bm\Sigma_{k,bb}-\bm\Sigma_{k,ba}\bm\Sigma_{k,aa}^{-1}\bm\Sigma_{k,ab}\big): \]

    a mixture of linear-Gaussian regressions with input-dependent weights, as in the mixtures of experts of Section 14.5.

11. Exercise 9.11 (★). For the mixture with covariances \(\epsilon\mathbf{I}\) of Section 9.3.2, show that as \(\epsilon\to0\) maximising the expected complete-data log likelihood (9.40) is equivalent to minimising the distortion \(J\) of (9.1).

    Solution

    By (9.40)–(9.41),

    \[ \E_{\mathbf{Z}}[\ln p]=\sum_{n,k}\gamma_{nk}\ln\pi_k-\frac{ND}{2}\ln(2\pi\epsilon)-\frac{1}{2\epsilon}\sum_{n,k}\gamma_{nk}\lVert\mathbf{x}_n-\bm\mu_k\rVert^2 . \]

    The middle term does not involve the parameters. As \(\epsilon\to0\), each \(\gamma_{nk}\to r_{nk}\) of (9.2) (for points with a unique nearest centre, at rate \(e^{-\Delta_n/2\epsilon}\), and assuming all \(\pi_k\gt0\)). Multiplying by \(\epsilon\gt0\) does not change the maximiser, and \(\epsilon\,\E_{\mathbf{Z}}[\ln p]\)\({}\to-\frac12\sum_{n,k}r_{nk}\lVert\mathbf{x}_n-\bm\mu_k\rVert^2\)\({}=-\frac12J\), uniformly over bounded sets of means (the \(\ln\pi_k\) term is bounded and \(\epsilon\ln\epsilon\to0\)). Hence maximising over the means becomes minimising \(J\), and the maximiser itself, the weighted mean (9.17), converges to the K-means update (9.4). Full details in Result 9.9.

    ▶ See it in a demo: EM with shared covariance \(\epsilon\mathbf{I}\) turning into K-means

12. Exercise 9.12 (★, www). For a mixture (9.82) whose components (discrete, continuous or mixed) have means \(\bm\mu_k\) and covariances \(\bm\Sigma_k\), show that the mixture's mean and covariance are (9.49) and (9.50).

    Solution

    For any function \(f\), \(\E[f(\mathbf{x})]\)\({}=\sum_k\pi_k\E_k[f(\mathbf{x})]\) by linearity of the sum or integral. With \(f(\mathbf{x})=\mathbf{x}\): \(\E[\mathbf{x}]=\sum_k\pi_k\bm\mu_k\), (9.49). With \(f(\mathbf{x})\)\({}=\mathbf{x}\mathbf{x}^\T\) and \(\E_k[\mathbf{x}\mathbf{x}^\T]\)\({}=\bm\Sigma_k+\bm\mu_k\bm\mu_k^\T\): \(\E[\mathbf{x}\mathbf{x}^\T]\)\({}=\sum_k\pi_k(\bm\Sigma_k+\bm\mu_k\bm\mu_k^\T)\), and \(\cov[\mathbf{x}]\)\({}=\E[\mathbf{x}\mathbf{x}^\T]-\E[\mathbf{x}]\E[\mathbf{x}]^\T\) gives (9.50). Equivalently \(\cov[\mathbf{x}]\)\({}=\sum_k\pi_k\bm\Sigma_k+\sum_k\pi_k(\bm\mu_k-\bar{\bm\mu})(\bm\mu_k-\bar{\bm\mu})^\T\) (Result 9.10).

13. Exercise 9.13 (★★). Show that a Bernoulli mixture whose parameters are at a maximum of the likelihood satisfies \(\E[\mathbf{x}]\)\({}=\frac1N\sum_n\mathbf{x}_n\equiv\bar{\mathbf{x}}\) (9.83). Show also that if all components start with the same mean \(\bm\mu_k=\widehat{\bm\mu}\), EM converges after one iteration, for any initial mixing coefficients, to \(\bm\mu_k=\bar{\mathbf{x}}\).

    Solution

    First part. At an interior maximum the likelihood is stationary, so by the argument of Result 9.5 (with the Bernoulli derivatives of gap box 9.12) the parameters satisfy \(\bm\mu_k=\bar{\mathbf{x}}_k\) and \(\pi_k=N_k/N\) with the responsibilities evaluated at the maximum. Then

    \[ \E[\mathbf{x}]=\sum_k\pi_k\bm\mu_k=\sum_k\frac{N_k}{N}\cdot\frac{1}{N_k}\sum_n\gamma_{nk}\mathbf{x}_n=\frac1N\sum_n\Big(\sum_k\gamma_{nk}\Big)\mathbf{x}_n=\bar{\mathbf{x}}, \]

    since \(\sum_k\gamma_{nk}=1\). (The same identity holds after every M step, maximum or not.)

    Second part. If \(\bm\mu_k=\widehat{\bm\mu}\) for all \(k\), then \(p(\mathbf{x}_n\given\bm\mu_k)\) is the same for every \(k\), and (9.56) gives \(\gamma_{nk}=\pi_k\) for all \(n\): the data carry no information about the labels. The M step gives \(N_k=N\pi_k\), \(\bm\mu_k^{\rm new}\)\({}=\frac{1}{N\pi_k}\sum_n\pi_k\mathbf{x}_n\)\({}=\bar{\mathbf{x}}\) and \(\pi_k^{\rm new}=\pi_k\). The components are again all equal (now to \(\bar{\mathbf{x}}\)), so the next E step again gives \(\gamma_{nk}=\pi_k\), and the next M step returns the same parameters: EM has converged after one iteration, for any \(\pi_k\gt0\), to a degenerate fixed point (a stationary point, generally not a maximum) in which all components equal the single-Bernoulli fit. Random, distinct initial means avoid it.

14. Exercise 9.14 (★). Show that marginalising \(p(\mathbf{x}\given\mathbf{z},\bm\mu)p(\mathbf{z}\given\bm\pi)\), (9.52)–(9.53), over \(\mathbf{z}\) gives (9.47).

    Solution

    \(p(\mathbf{x}\given\mathbf{z},\bm\mu)p(\mathbf{z}\given\bm\pi)\)\({}=\prod_k\{\pi_kp(\mathbf{x}\given\bm\mu_k)\}^{z_k}\). For \(\mathbf{z}=\mathbf{e}_j\) this is \(\pi_jp(\mathbf{x}\given\bm\mu_j)\), and summing over the \(K\) values of \(\mathbf{z}\) gives \(\sum_j\pi_jp(\mathbf{x}\given\bm\mu_j)\), which is (9.47). This is Exercise 9.3 with Bernoulli components; nothing about Gaussians was used there.

15. Exercise 9.15 (★, www). Show that maximising (9.55) over \(\bm\mu_k\) gives (9.59).

    Solution

    \(\partial(9.55)/\partial\mu_{ki}\)\({}=\sum_n\gamma_{nk}\{x_{ni}/\mu_{ki}-(1-x_{ni})/(1-\mu_{ki})\}\)\({}=\{\sum_n\gamma_{nk}x_{ni}-\mu_{ki}N_k\}/\{\mu_{ki}(1-\mu_{ki})\}\), which is zero at \(\mu_{ki}\)\({}=N_k^{-1}\sum_n\gamma_{nk}x_{ni}\); in vector form \(\bm\mu_k=\bar{\mathbf{x}}_k\), (9.59). The objective in \(\mu_{ki}\) is \(A\ln\mu_{ki}+B\ln(1-\mu_{ki})\) with \(A,B\ge0\), concave, so this stationary point is the maximum (gap box 9.12).

16. Exercise 9.16 (★). Show that maximising (9.55) over the mixing coefficients, with a Lagrange multiplier for \(\sum_k\pi_k=1\), gives (9.60).

    Solution

    The \(\bm\pi\)-terms of (9.55) are \(\sum_kN_k\ln\pi_k\). Setting the derivative of \(\sum_kN_k\ln\pi_k+\lambda(\sum_k\pi_k-1)\) to zero gives \(N_k/\pi_k+\lambda=0\), i.e. \(N_k=-\lambda\pi_k\). Summing over \(k\), with \(\sum_kN_k=N\) (responsibilities sum to one over \(k\)) and \(\sum_k\pi_k=1\), gives \(\lambda=-N\), hence \(\pi_k=N_k/N\), which is (9.60). The value \(\sum_kN_k\ln\pi_k\)\({}=N\sum_kp_k\ln p_k-N\KL(\mathbf{p}\|\bm\pi)\) with \(p_k=N_k/N\) shows it is the global maximum.

17. Exercise 9.17 (★, www). Using \(0\)\({}\le p(\mathbf{x}_n\given\bm\mu_k)\)\({}\le1\), show that the log likelihood of a Bernoulli mixture is bounded above, so it has no singularities at which it tends to infinity.

    Solution

    Each \(p(\mathbf{x}_n\given\bm\mu_k)\) is a product of numbers \(\mu_{ki}\) or \(1-\mu_{ki}\) in \([0,1]\), so it lies in \([0,1]\). Hence \(\sum_k\pi_kp(\mathbf{x}_n\given\bm\mu_k)\)\({}\le\sum_k\pi_k\)\({}=1\), each term of (9.51) is at most 0, and \(\ln p(\mathbf{X}\given\bm\mu,\bm\pi)\)\({}\le0\). The parameter set is compact and the likelihood continuous, so the maximum is even attained (Result 9.11).

18. Exercise 9.18 (★★). For a Bernoulli mixture with a Beta prior \(\mathrm{Beta}(\mu_{ki}\given a_k,b_k)\) on each parameter (2.13) and a Dirichlet prior \(\Dir(\bm\pi\given\bm\alpha)\) (2.38), derive the EM algorithm that maximises the posterior \(p(\bm\mu,\bm\pi\given\mathbf{X})\).

    Solution

    By Exercise 9.4, the E step is (9.56) and the M step maximises (9.55) plus the log prior

    \[ \sum_k(\alpha_k-1)\ln\pi_k+\sum_{k,i}\big\{(a_k-1)\ln\mu_{ki}+(b_k-1)\ln(1-\mu_{ki})\big\}. \]

    The \(\mu_{ki}\)-terms are \(A\ln\mu_{ki}+B\ln(1-\mu_{ki})\) with \(A=\sum_n\gamma_{nk}x_{ni}+a_k-1\), \(B\)\({}=\sum_n\gamma_{nk}(1-x_{ni})+b_k-1\); for \(a_k,b_k\ge1\) the maximiser is

    \[ \mu_{ki}=\frac{\sum_n\gamma_{nk}x_{ni}+a_k-1}{N_k+a_k+b_k-2}. \]

    The \(\bm\pi\)-terms are \(\sum_k(N_k+\alpha_k-1)\ln\pi_k\); with a Lagrange multiplier, \(\lambda=-(N+\hat\alpha-K)\) and

    \[ \pi_k=\frac{N_k+\alpha_k-1}{N+\hat\alpha-K},\qquad\hat\alpha=\sum_k\alpha_k, \]

    valid (a global maximum) when all \(\alpha_k\ge1\). With \(a_k,b_k\gt1\) every \(\mu_{ki}\) stays strictly inside \((0,1)\). Flat priors (\(a_k=b_k=\alpha_k=1\)) give back (9.59)–(9.60). Details and the case \(\alpha_k\lt1\) in gap box 9.13.

19. Exercise 9.19 (★★). Let \(\mathbf{x}\) consist of \(D\) multinomial variables, each a 1-of-M binary vector \((x_{i1},\dots,x_{iM})\), with mixture distribution (9.84)–(9.85): \(p(\mathbf{x}\given\bm\mu_k)\)\({}=\prod_i\prod_j\mu_{kij}^{x_{ij}}\), \(\sum_j\mu_{kij}=1\). Derive the E and M steps of maximum likelihood EM.

    Solution

    E step: \(\gamma_{nk}\)\({}=\pi_k\prod_{i,j}\mu_{kij}^{x_{nij}}\big/\sum_l\pi_l\prod_{i,j}\mu_{lij}^{x_{nij}}\) (in log space). M step: the expected complete-data log likelihood is \(\sum_{n,k}\gamma_{nk}\{\ln\pi_k+\sum_{i,j}x_{nij}\ln\mu_{kij}\}\). For each \((k,i)\), maximise \(\sum_j c_{kij}\ln\mu_{kij}\) with \(c_{kij}=\sum_n\gamma_{nk}x_{nij}\) subject to \(\sum_j\mu_{kij}=1\): the Lagrange condition \(c_{kij}/\mu_{kij}+\lambda_{ki}=0\), summed over \(j\) after multiplying by \(\mu_{kij}\), gives \(\lambda_{ki}=-\sum_jc_{kij}=-N_k\) (because \(\sum_jx_{nij}=1\)). Hence

    \[ \mu_{kij}=\frac{1}{N_k}\sum_{n=1}^{N}\gamma_{nk}\,x_{nij},\qquad\pi_k=\frac{N_k}{N}, \]

    global maxima by Gibbs' inequality. For \(M=2\) this is the Bernoulli mixture. With Dirichlet priors \(\Dir(\bm\mu_{ki}\given\bm\beta_{ki})\), \(\mu_{kij}\)\({}=(c_{kij}+\beta_{kij}-1)/(N_k+\sum_j\beta_{kij}-M)\) (gap box 9.13).

20. Exercise 9.20 (★, www). Show that maximising (9.62) over \(\alpha\) gives (9.63).

    Solution

    The \(\alpha\)-terms of (9.62) are \(\frac M2\ln\alpha-\frac\alpha2\E[\mathbf{w}^\T\mathbf{w}]\). The derivative \(\frac{M}{2\alpha}-\frac12\E[\mathbf{w}^\T\mathbf{w}]\) vanishes at \(\alpha\)\({}=M/\E[\mathbf{w}^\T\mathbf{w}]\), and the second derivative \(-M/(2\alpha^2)\) is negative. Under the posterior \(\N(\mathbf{m}_N,\mathbf{S}_N)\), \(\E[\mathbf{w}^\T\mathbf{w}]\)\({}=\Tr(\mathbf{S}_N+\mathbf{m}_N\mathbf{m}_N^\T)\)\({}=\mathbf{m}_N^\T\mathbf{m}_N+\Tr(\mathbf{S}_N)\), which gives (9.63).

21. Exercise 9.21 (★★). In the same setting (Section 3.5), derive the M-step update for the noise precision \(\beta\), analogous to (9.63).

    Solution

    The \(\beta\)-terms of (9.62) are \(\frac N2\ln\beta-\frac\beta2\sum_n\E[(t_n-\mathbf{w}^\T\bm\phi_n)^2]\). Setting the derivative to zero (the second derivative \(-N/(2\beta^2)\) is negative) and using \(\E[(t_n-\mathbf{w}^\T\bm\phi_n)^2]\)\({}=(t_n-\mathbf{m}_N^\T\bm\phi_n)^2+\bm\phi_n^\T\mathbf{S}_N\bm\phi_n\):

    \[ \frac{1}{\beta^{\rm new}}=\frac1N\Big\{\lVert\mathbf{t}-\bm\Phi\mathbf{m}_N\rVert^2+\Tr(\bm\Phi^\T\bm\Phi\mathbf{S}_N)\Big\}=\frac1N\Big\{\lVert\mathbf{t}-\bm\Phi\mathbf{m}_N\rVert^2+\frac{\gamma}{\beta}\Big\}, \]

    the last form because \(\Tr(\bm\Phi^\T\bm\Phi\mathbf{S}_N)\)\({}=\gamma/\beta\) with \(\gamma\) of (3.91) evaluated at the old values. At a fixed point this is \((N-\gamma)/\beta\)\({}=\lVert\mathbf{t}-\bm\Phi\mathbf{m}_N\rVert^2\), the evidence update (3.95) (Result 9.12).

22. Exercise 9.22 (★★). By maximising the expected complete-data log likelihood (9.66), derive the RVM updates (9.67) and (9.68).

    Solution

    Under the posterior \(\N(\mathbf{m},\bm\Sigma)\) at the old hyperparameters, (9.66) equals, up to constants,

    \[ \sum_i\Big[\frac12\ln\alpha_i-\frac{\alpha_i}{2}(m_i^2+\Sigma_{ii})\Big]+\frac N2\ln\beta-\frac\beta2\Big\{\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2+\Tr(\bm\Phi^\T\bm\Phi\bm\Sigma)\Big\}, \]

    using \(\E[w_i^2]=m_i^2+\Sigma_{ii}\) and the bias-plus-variance decomposition of the expected squared residuals. The derivative in \(\alpha_i\), \(\frac{1}{2\alpha_i}-\frac12(m_i^2+\Sigma_{ii})\), gives (9.67); the one in \(\beta\) gives \(N/\beta\)\({}=\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2+\Tr(\bm\Phi^\T\bm\Phi\bm\Sigma)\). Since \(\beta\bm\Phi^\T\bm\Phi\)\({}=\bm\Sigma^{-1}-\mathbf{A}\) at the old values, \(\Tr(\bm\Phi^\T\bm\Phi\bm\Sigma)\)\({}=\beta^{-1}\Tr(\mathbf{I}-\mathbf{A}\bm\Sigma)\)\({}=\beta^{-1}\sum_i\gamma_i\), which gives (9.68) (Result 9.13).

23. Exercise 9.23 (★★, www). Show that, at any stationary point, the RVM re-estimation equations obtained by direct maximisation of the evidence, (7.87)–(7.88), and those obtained by EM, (9.67)–(9.68), are equivalent. (Wording as corrected in the errata: the equivalence holds at stationary points only.)

    Solution

    At a stationary point the old and new hyperparameters coincide. Then (9.67) says \(\alpha_i(m_i^2+\Sigma_{ii})=1\), i.e. \(\alpha_im_i^2\)\({}=1-\alpha_i\Sigma_{ii}\)\({}=\gamma_i\) by (7.89), i.e. \(\alpha_i=\gamma_i/m_i^2\), which is (7.87). And (9.68) says \(N\beta^{-1}\)\({}=\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2+\beta^{-1}\sum_i\gamma_i\), i.e. \(\beta^{-1}(N-\sum_i\gamma_i)\)\({}=\lVert\mathbf{t}-\bm\Phi\mathbf{m}\rVert^2\), which is (7.88). Each step is reversible, so the two sets of equations have the same solutions. Away from stationary points the updates differ, which is why the original wording ("formally equivalent", without qualification) needed the correction.

24. Exercise 9.24 (★). Verify (9.70) with \(\LB(q,\bm\theta)\) and \(\KL(q\|p)\) defined by (9.71) and (9.72).

    Solution

    For \(\mathbf{Z}\) with \(q(\mathbf{Z})\gt0\), the product rule gives \(\ln\frac{p(\mathbf{X},\mathbf{Z}\given\bm\theta)}{q(\mathbf{Z})}\)\({}=\ln\frac{p(\mathbf{Z}\given\mathbf{X},\bm\theta)}{q(\mathbf{Z})}+\ln p(\mathbf{X}\given\bm\theta)\). Multiply by \(q(\mathbf{Z})\) and sum: \(\LB(q,\bm\theta)\)\({}=-\KL(q\|p)+\ln p(\mathbf{X}\given\bm\theta)\sum_{\mathbf{Z}}q(\mathbf{Z})\)\({}=-\KL(q\|p)+\ln p(\mathbf{X}\given\bm\theta)\), which rearranges to (9.70) (Result 9.14, which also treats \(q\) with zeros).

25. Exercise 9.25 (★, www). Show that the lower bound \(\LB(q,\bm\theta)\) with \(q(\mathbf{Z})\)\({}=p(\mathbf{Z}\given\mathbf{X},\bm\theta^{\rm old})\) has the same gradient with respect to \(\bm\theta\) as \(\ln p(\mathbf{X}\given\bm\theta)\) at \(\bm\theta=\bm\theta^{\rm old}\).

    Solution

    By (9.70), \(\LB(q,\bm\theta)\)\({}=\ln p(\mathbf{X}\given\bm\theta)-\KL(q\|p(\cdot\given\mathbf{X},\bm\theta))\). With \(q\) fixed at the posterior for \(\bm\theta^{\rm old}\), the function \(\bm\theta\mapsto\KL(q\|p(\cdot\given\mathbf{X},\bm\theta))\) is nonnegative and equals zero at \(\bm\theta^{\rm old}\), so \(\bm\theta^{\rm old}\) minimises it and (if it is differentiable there) its gradient vanishes. Hence \(\nabla\LB=\nabla\ln p\) at \(\bm\theta^{\rm old}\). Directly: \(\nabla\LB|_{\bm\theta^{\rm old}}\)\({}=\sum_{\mathbf{Z}}\frac{p(\mathbf{X},\mathbf{Z}\given\bm\theta^{\rm old})}{p(\mathbf{X}\given\bm\theta^{\rm old})}\frac{\nabla p(\mathbf{X},\mathbf{Z}\given\bm\theta^{\rm old})}{p(\mathbf{X},\mathbf{Z}\given\bm\theta^{\rm old})}\)\({}=\frac{\nabla p(\mathbf{X}\given\bm\theta^{\rm old})}{p(\mathbf{X}\given\bm\theta^{\rm old})}\) (Result 9.16).

    ▶ See it in a demo: the bound touching the log likelihood with a common tangent

26. Exercise 9.26 (★, www). In incremental EM for a Gaussian mixture, only the responsibilities of one point \(\mathbf{x}_m\) are recomputed. Starting from (9.17) and (9.18), derive (9.78) and (9.79).

    Solution

    Let \(\Delta_k\)\({}=\gamma^{\rm new}(z_{mk})-\gamma^{\rm old}(z_{mk})\). In \(N_k=\sum_n\gamma(z_{nk})\) only the \(m\)th term changes, so \(N_k^{\rm new}\)\({}=N_k^{\rm old}+\Delta_k\), (9.79). From (9.17), \(N_k^{\rm new}\bm\mu_k^{\rm new}\)\({}=\sum_{n\ne m}\gamma^{\rm old}(z_{nk})\mathbf{x}_n+\gamma^{\rm new}(z_{mk})\mathbf{x}_m\)\({}=N_k^{\rm old}\bm\mu_k^{\rm old}+\Delta_k\mathbf{x}_m\). Substituting \(N_k^{\rm old}\)\({}=N_k^{\rm new}-\Delta_k\) and dividing by \(N_k^{\rm new}\): \(\bm\mu_k^{\rm new}\)\({}=\bm\mu_k^{\rm old}+\frac{\Delta_k}{N_k^{\rm new}}(\mathbf{x}_m-\bm\mu_k^{\rm old})\), (9.78).

27. Exercise 9.27 (★★). Derive the incremental updates of the covariance matrices and mixing coefficients, analogous to (9.78).

    Solution

    With \(\mathbf{W}_k\)\({}=\sum_n\gamma(z_{nk})\mathbf{x}_n\mathbf{x}_n^\T\), (9.25) gives \(N_k\bm\Sigma_k\)\({}=\mathbf{W}_k-N_k\bm\mu_k\bm\mu_k^\T\), and \(\mathbf{W}_k^{\rm new}\)\({}=\mathbf{W}_k^{\rm old}+\Delta_k\mathbf{x}_m\mathbf{x}_m^\T\). Using the mean update, the algebra of Result 9.19 gives

    \[ N_k^{\rm new}\bm\Sigma_k^{\rm new}=N_k^{\rm old}\bm\Sigma_k^{\rm old}+\frac{\Delta_kN_k^{\rm old}}{N_k^{\rm new}}(\mathbf{x}_m-\bm\mu_k^{\rm old})(\mathbf{x}_m-\bm\mu_k^{\rm old})^\T, \]

    i.e.

    \[ \bm\Sigma_k^{\rm new}=\bm\Sigma_k^{\rm old}+\frac{\Delta_k}{N_k^{\rm new}}\Big[(\mathbf{x}_m-\bm\mu_k^{\rm new})(\mathbf{x}_m-\bm\mu_k^{\rm old})^\T-\bm\Sigma_k^{\rm old}\Big]. \]

    For the mixing coefficients, \(\pi_k=N_k/N\) with \(N\) fixed, so \(\pi_k^{\rm new}\)\({}=\pi_k^{\rm old}+\Delta_k/N\); they still sum to one because \(\sum_k\Delta_k=0\).

Demonstrations for this chapter

Each demonstration runs in the browser, with seeded random numbers 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 under discussion.

The Python code that produces every figure of this page, and checks the chapter's identities numerically, is code/ch09_mixture_em.py (see the code index).

References

1. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer, New York. Chapter 9, pp. 423–459 (and the official errata).
2. Aloise, D., Deshpande, A., Hansen, P. and Popat, P. (2009). NP-hardness of Euclidean sum-of-squares clustering. Machine Learning, 75(2), 245–248. doi:10.1007/s10994-009-5103-0.
3. Arthur, D. and Vassilvitskii, S. (2007). k-means++: the advantages of careful seeding. In Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithms, 1027–1035.
4. Boyles, R. A. (1983). On the convergence of the EM algorithm. Journal of the Royal Statistical Society: Series B, 45(1), 47–50. doi:10.1111/j.2517-6161.1983.tb01229.x.
5. Celeux, G. and Govaert, G. (1992). A classification EM algorithm for clustering and two stochastic versions. Computational Statistics & Data Analysis, 14(3), 315–332. doi:10.1016/0167-9473(92)90042-E.
6. Csiszár, I. and Tusnády, G. (1984). Information geometry and alternating minimization procedures. Statistics and Decisions, Supplement Issue 1, 205–237.
7. Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B, 39(1), 1–38. doi:10.1111/j.2517-6161.1977.tb01600.x.
8. Gersho, A. and Gray, R. M. (1992). Vector Quantization and Signal Compression. Kluwer, Boston. doi:10.1007/978-1-4615-3626-0.
9. Hathaway, R. J. (1985). A constrained formulation of maximum-likelihood estimation for normal mixture distributions. The Annals of Statistics, 13(2), 795–800. doi:10.1214/aos/1176349557.
10. Hathaway, R. J. (1986). Another interpretation of the EM algorithm for mixture distributions. Statistics & Probability Letters, 4(2), 53–56. doi:10.1016/0167-7152(86)90016-7.
11. Kaufman, L. and Rousseeuw, P. J. (1990). Finding Groups in Data: An Introduction to Cluster Analysis. Wiley, New York. doi:10.1002/9780470316801.
12. Kiefer, J. and Wolfowitz, J. (1956). Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters. The Annals of Mathematical Statistics, 27(4), 887–906. doi:10.1214/aoms/1177728066.
13. Lazarsfeld, P. F. and Henry, N. W. (1968). Latent Structure Analysis. Houghton Mifflin, Boston.
14. Little, R. J. A. and Rubin, D. B. (2019). Statistical Analysis with Missing Data, 3rd edition. Wiley, Hoboken, NJ. doi:10.1002/9781119482260.
15. Lloyd, S. P. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129–137. doi:10.1109/TIT.1982.1056489.
16. MacQueen, J. (1967). Some methods for classification and analysis of multivariate observations. In L. M. Le Cam and J. Neyman (eds), Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1, 281–297. University of California Press, Berkeley.
17. McLachlan, G. J. and Krishnan, T. (2008). The EM Algorithm and Extensions, 2nd edition. Wiley, Hoboken, NJ. doi:10.1002/9780470191613.
18. McLachlan, G. J. and Peel, D. (2000). Finite Mixture Models. Wiley, New York. doi:10.1002/0471721182.
19. Meng, X.-L. and Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: a general framework. Biometrika, 80(2), 267–278. doi:10.1093/biomet/80.2.267.
20. Neal, R. M. and Hinton, G. E. (1998). A view of the EM algorithm that justifies incremental, sparse, and other variants. In M. I. Jordan (ed.), Learning in Graphical Models, 355–368. Kluwer, Dordrecht. doi:10.1007/978-94-011-5014-9_12. (Bishop cites the 1999 MIT Press printing of the same volume.)
21. Robbins, H. and Monro, S. (1951). A stochastic approximation method. The Annals of Mathematical Statistics, 22(3), 400–407. doi:10.1214/aoms/1177729586.
22. Rubin, D. B. (1976). Inference and missing data. Biometrika, 63(3), 581–592. doi:10.1093/biomet/63.3.581.
23. Sung, K.-K. and Poggio, T. (1994). Example-based learning for view-based human face detection. A.I. Memo 1521, Massachusetts Institute of Technology.
24. Teicher, H. (1963). Identifiability of finite mixtures. The Annals of Mathematical Statistics, 34(4), 1265–1269. doi:10.1214/aoms/1177703862.
25. Tipping, M. E. (2001). Sparse Bayesian learning and the relevance vector machine. Journal of Machine Learning Research, 1, 211–244.
26. Wu, C. F. J. (1983). On the convergence properties of the EM algorithm. The Annals of Statistics, 11(1), 95–103. doi:10.1214/aos/1176346060.
27. Yakowitz, S. J. and Spragins, J. D. (1968). On the identifiability of finite mixtures. The Annals of Mathematical Statistics, 39(1), 209–214. doi:10.1214/aoms/1177698520.

DOIs were checked against Crossref. The numerical statements made in the text (the rates of convergence in Result 9.17, the checks of Results 9.12, 9.16 and 9.19) come from short Python computations written for these notes; the figures come from the course code.