9
Mixture Models and EM
Section 9.1
If we define a joint distribution over observed and latent variables, the correspond- ing distribution of the observed variables alone is obtained by marginalization. This allows relatively complex marginal distributions over observed variables to be ex- pressed in terms of more tractable joint distributions over the expanded space of observed and latent variables. The introduction of latent variables thereby allows complicated distributions to be formed from simpler components. In this chapter, we shall see that mixture distributions, such as the Gaussian mixture discussed in Section 2.3.9, can be interpreted in terms of discrete latent variables. Continuous latent variables will form the subject of Chapter 12.
As well as providing a framework for building more complex probability dis- tributions, mixture models can also be used to cluster data. We therefore begin our discussion of mixture distributions by considering the problem of finding clusters in a set of data points, which we approach first using a nonprobabilistic technique called the K-means algorithm (Lloyd, 1982). Then we introduce the latent variable
423
424
9. MIXTURE MODELS AND EM
Section 9.2
Section 9.3 Section 9.4
view of mixture distributions in which the discrete latent variables can be interpreted as defining assignments of data points to specific components of the mixture. A gen- eral technique for finding maximum likelihood estimators in latent variable models is the expectation-maximization (EM) algorithm. We first of all use the Gaussian mixture distribution to motivate the EM algorithm in a fairly informal way, and then we give a more careful treatment based on the latent variable viewpoint. We shall see that the K-means algorithm corresponds to a particular nonprobabilistic limit of EM applied to mixtures of Gaussians. Finally, we discuss EM in some generality.
Gaussian mixture models are widely used in data mining, pattern recognition, machine learning, and statistical analysis. In many applications, their parameters are determined by maximum likelihood, typically using the EM algorithm. However, as we shall see there are some significant limitations to the maximum likelihood ap- proach, and in Chapter 10 we shall show that an elegant Bayesian treatment can be given using the framework of variational inference. This requires little additional computation compared with EM, and it resolves the principal difficulties of maxi- mum likelihood while also allowing the number of components in the mixture to be inferred automatically from the data.
9.1. K-means Clustering
We begin by considering the problem of identifying groups, or clusters, of data points in a multidimensional space. Suppose we have a data set {x1 , . . . , xN } consisting of N observations of a random D-dimensional Euclidean variable x. Our goal is to partition the data set into some number K of clusters, where we shall suppose for the moment that the value of K is given. Intuitively, we might think of a cluster as comprising a group of data points whose inter-point distances are small compared with the distances to points outside of the cluster. We can formalize this notion by first introducing a set of D-dimensional vectors μk, where k = 1, . . . , K, in which μk is a prototype associated with the kth cluster. As we shall see shortly, we can think of the μk as representing the centres of the clusters. Our goal is then to find an assignment of data points to clusters, as well as a set of vectors {μk}, such that the sum of the squares of the distances of each data point to its closest vector μk, is a minimum.
It is convenient at this point to define some notation to describe the assignment of data points to clusters. For each data point xn, we introduce a corresponding set of binary indicator variables rnk ∈ {0, 1}, where k = 1, . . . , K describing which of the K clusters the data point xn is assigned to, so that if data point xn is assigned to cluster k then rnk = 1, and rnj = 0 for j à= k. This is known as the 1-of-K coding scheme. We can then define an objective function, sometimes called a distortion measure, given by
N K
rnk xn−μk 2 (9.1) which represents the sum of the squares of the distances of each data point to its
J=
n=1 k=1
Section 9.4
assigned vector μk. Our goal is to find values for the {rnk} and the {μk} so as to minimize J. We can do this through an iterative procedure in which each iteration involves two successive steps corresponding to successive optimizations with respect to the rnk and the μk . First we choose some initial values for the μk . Then in the first phase we minimize J with respect to the rnk, keeping the μk fixed. In the second phase we minimize J with respect to the μk, keeping rnk fixed. This two-stage optimization is then repeated until convergence. We shall see that these two stages of updating rnk and updating μk correspond respectively to the E (expectation) and M (maximization) steps of the EM algorithm, and to emphasize this we shall use the terms E step and M step in the context of the K-means algorithm.
Consider first the determination of the rnk. Because J in (9.1) is a linear func- tion of rnk, this optimization can be performed easily to give a closed form solution. The terms involving different n are independent and so we can optimize for each n separately by choosing rnk to be 1 for whichever value of k gives the minimum value of xn − μk 2. In other words, we simply assign the nth data point to the closest cluster centre. More formally, this can be expressed as
1 ifk=argminj xn−μj 2
rnk = 0 otherwise. (9.2)
Now consider the optimization of the μk with the rnk held fixed. The objective function J is a quadratic function of μk, and it can be minimized by setting its derivative with respect to μk to zero giving
(9.3)
(9.4)
The denominator in this expression is equal to the number of points assigned to cluster k, and so this result has a simple interpretation, namely set μk equal to the mean of all of the data points xn assigned to cluster k. For this reason, the procedure is known as the K-means algorithm.
The two phases of re-assigning data points to clusters and re-computing the clus- ter means are repeated in turn until there is no further change in the assignments (or until some maximum number of iterations is exceeded). Because each phase reduces the value of the objective function J, convergence of the algorithm is assured. How- ever, it may converge to a local rather than global minimum of J. The convergence properties of the K-means algorithm were studied by MacQueen (1967).
The K-means algorithm is illustrated using the Old Faithful data set in Fig- ure 9.1. For the purposes of this example, we have made a linear re-scaling of the data, known as standardizing, such that each of the variables has zero mean and unit standard deviation. For this example, we have chosen K = 2, and so in this
rnk(xn − μk) = 0 which we can easily solve for μk to give
2
N n=1
rnkxn n
μk=r . n nk
9.1. K-means Clustering 425
Exercise 9.1 Appendix A
426
9. MIXTURE MODELS AND EM
(a)
2
0
−2
−2 0 2
(b)
2
0
−2
−2 0 2
(c)
2
0
−2
−2 0 2
(d)
2
0
−2
−2 0 2
(e)
2
0
−2
−2 0 2
(f)
2
0
−2
−2 0 2
(g)
2
0
−2
−2 0 2
(h)
2
0
−2
−2 0 2
Illustration of the K-means algorithm using the re-scaled Old Faithful data set. (a) Green points denote the data set in a two-dimensional Euclidean space. The initial choices for centres μ1 and μ2 are shown by the red and blue crosses, respectively. (b) In the initial E step, each data point is assigned either to the red cluster or to the blue cluster, according to which cluster centre is nearer. This is equivalent to classifying the points according to which side of the perpendicular bisector of the two cluster centres, shown by the magenta line, they lie on. (c) In the subsequent M step, each cluster centre is re-computed to be the mean of the points assigned to the corresponding cluster. (d)–(i) show successive E and M steps through to final convergence of the algorithm.
Figure 9.1
(i)
2
0
−2
−2 0 2
Section 9.2.2
Section 2.3.5 Exercise 9.2
Figure 9.2
Plot of the cost function J given by (9.1) after each E step (blue points) and M step (red points) of the K- means algorithm for the example shown in Figure 9.1. The algo- rithm has converged after the third M step, and the final EM cycle pro- duces no changes in either the as- signments or the prototype vectors.
case, the assignment of each data point to the nearest cluster centre is equivalent to a classification of the data points according to which side they lie of the perpendicular bisector of the two cluster centres. A plot of the cost function J given by (9.1) for the Old Faithful example is shown in Figure 9.2.
Note that we have deliberately chosen poor initial values for the cluster centres so that the algorithm takes several steps before convergence. In practice, a better initialization procedure would be to choose the cluster centres μk to be equal to a random subset of K data points. It is also worth noting that the K-means algorithm itself is often used to initialize the parameters in a Gaussian mixture model before applying the EM algorithm.
A direct implementation of the K-means algorithm as discussed here can be relatively slow, because in each E step it is necessary to compute the Euclidean dis- tance between every prototype vector and every data point. Various schemes have been proposed for speeding up the K-means algorithm, some of which are based on precomputing a data structure such as a tree such that nearby points are in the same subtree (Ramasubramanian and Paliwal, 1990; Moore, 2000). Other approaches make use of the triangle inequality for distances, thereby avoiding unnecessary dis- tance calculations (Hodgson, 1998; Elkan, 2003).
So far, we have considered a batch version of K-means in which the whole data set is used together to update the prototype vectors. We can also derive an on-line stochastic algorithm (MacQueen, 1967) by applying the Robbins-Monro procedure to the problem of finding the roots of the regression function given by the derivatives of J in (9.1) with respect to μk. This leads to a sequential update in which, for each data point xn in turn, we update the nearest prototype μk using
μnew = μold + η (x − μold) (9.5) kknnk
where ηn is the learning rate parameter, which is typically made to decrease mono- tonically as more data points are considered.
The K-means algorithm is based on the use of squared Euclidean distance as the measure of dissimilarity between a data point and a prototype vector. Not only does this limit the type of data variables that can be considered (it would be inappropriate for cases where some or all of the variables represent categorical labels for instance),
9.1. K-means Clustering 427
1000
J
500
0
1234
428 9. MIXTURE MODELS AND EM
Section 2.3.7
but it can also make the determination of the cluster means nonrobust to outliers. We can generalize the K-means algorithm by introducing a more general dissimilarity measure V(x,x) between two vectors x and x and then minimizing the following distortion measure
N K J =
n=1 k=1
which gives the K-medoids algorithm. The E step again involves, for given cluster prototypes μk, assigning each data point to the cluster for which the dissimilarity to the corresponding prototype is smallest. The computational cost of this is O(KN), as is the case for the standard K-means algorithm. For a general choice of dissimi- larity measure, the M step is potentially more complex than for K-means, and so it is common to restrict each cluster prototype to be equal to one of the data vectors as- signed to that cluster, as this allows the algorithm to be implemented for any choice of dissimilarity measure V(·,·) so long as it can be readily evaluated. Thus the M step involves, for each cluster k, a discrete search over the Nk points assigned to that cluster, which requires O(Nk2) evaluations of V(·, ·).
One notable feature of the K-means algorithm is that at each iteration, every data point is assigned uniquely to one, and only one, of the clusters. Whereas some data points will be much closer to a particular centre μk than to any other centre, there may be other data points that lie roughly midway between cluster centres. In the latter case, it is not clear that the hard assignment to the nearest cluster is the most appropriate. We shall see in the next section that by adopting a probabilistic approach, we obtain ‘soft’ assignments of data points to clusters in a way that reflects the level of uncertainty over the most appropriate assignment. This probabilistic formulation brings with it numerous benefits.
9.1.1 Image segmentation and compression
As an illustration of the application of the K-means algorithm, we consider the related problems of image segmentation and image compression. The goal of segmentation is to partition an image into regions each of which has a reasonably homogeneous visual appearance or which corresponds to objects or parts of objects (Forsyth and Ponce, 2003). Each pixel in an image is a point in a 3-dimensional space comprising the intensities of the red, blue, and green channels, and our segmentation algorithm simply treats each pixel in the image as a separate data point. Note that strictly this space is not Euclidean because the channel intensities are bounded by the interval [0, 1]. Nevertheless, we can apply the K -means algorithm without diffi- culty. We illustrate the result of running K-means to convergence, for any particular value of K, by re-drawing the image replacing each pixel vector with the {R, G, B} intensity triplet given by the centre μk to which that pixel has been assigned. Results for various values of K are shown in Figure 9.3. We see that for a given value of K, the algorithm is representing the image using a palette of only K colours. It should be emphasized that this use of K-means is not a particularly sophisticated approach to image segmentation, not least because it takes no account of the spatial proximity of different pixels. The image segmentation problem is in general extremely difficult
rnkV(xn,μk) (9.6)
9.1. K-means Clustering 429
Original image
K=2
K=3
K = 10
Figure 9.3 Two examples of the application of the K-means clustering algorithm to image segmentation show- ing the initial images together with their K-means segmentations obtained using various values of K. This also illustrates of the use of vector quantization for data compression, in which smaller values of K give higher compression at the expense of poorer image quality.
and remains the subject of active research and is introduced here simply to illustrate the behaviour of the K-means algorithm.
We can also use the result of a clustering algorithm to perform data compres- sion. It is important to distinguish between lossless data compression, in which the goal is to be able to reconstruct the original data exactly from the compressed representation, and lossy data compression, in which we accept some errors in the reconstruction in return for higher levels of compression than can be achieved in the lossless case. We can apply the K-means algorithm to the problem of lossy data compression as follows. For each of the N data points, we store only the identity k of the cluster to which it is assigned. We also store the values of the K clus- ter centres μk, which typically requires significantly less data, provided we choose K à N. Each data point is then approximated by its nearest centre μk. New data points can similarly be compressed by first finding the nearest μk and then storing the label k instead of the original data vector. This framework is often called vector quantization, and the vectors μk are called code-book vectors.
430 9. MIXTURE MODELS AND EM
The image segmentation problem discussed above also provides an illustration of the use of clustering for data compression. Suppose the original image has N pixels comprising {R, G, B} values each of which is stored with 8 bits of precision. Then to transmit the whole image directly would cost 24N bits. Now suppose we first run K-means on the image data, and then instead of transmitting the original pixel intensity vectors we transmit the identity of the nearest vector μk. Because there are K such vectors, this requires log2 K bits per pixel. We must also transmit the K code book vectors μk, which requires 24K bits, and so the total number of bits required to transmit the image is 24K + N log2 K (rounding up to the nearest integer). The original image shown in Figure 9.3 has 240 × 180 = 43, 200 pixels and so requires 24 × 43, 200 = 1, 036, 800 bits to transmit directly. By comparison, the compressed images require 43, 248 bits (K = 2), 86, 472 bits (K = 3), and 173, 040 bits (K = 10), respectively, to transmit. These represent compression ratios compared to the original image of 4.2%, 8.3%, and 16.7%, respectively. We see that there is a trade-off between degree of compression and image quality. Note that our aim in this example is to illustrate the K-means algorithm. If we had been aiming to produce a good image compressor, then it would be more fruitful to consider small blocks of adjacent pixels, for instance 5 × 5, and thereby exploit the correlations that exist in natural images between nearby pixels.
9.2. Mixtures of Gaussians
In Section 2.3.9 we motivated the Gaussian mixture model as a simple linear super- position of Gaussian components, aimed at providing a richer class of density mod- els than the single Gaussian. We now turn to a formulation of Gaussian mixtures in terms of discrete latent variables. This will provide us with a deeper insight into this important distribution, and will also serve to motivate the expectation-maximization algorithm.
Recall from (2.188) that the Gaussian mixture distribution can be written as a linear superposition of Gaussians in the form
K
πkN(x|μk,Σk). (9.7) Let us introduce a K-dimensional binary random variable z having a 1-of-K repre-
sentation in which a particular element z
equal to 0. The values of zk therefore satisfy zk ∈ {0, 1} and k zk = 1, and we see that there are K possible states for the vector z according to which element is nonzero. We shall define the joint distribution p(x, z) in terms of a marginal dis- tribution p(z) and a conditional distribution p(x|z), corresponding to the graphical model in Figure 9.4. The marginal distribution over z is specified in terms of the mixing coefficients πk, such that
p(x) =
k=1
k
is equal to 1 and all other elements are
p(zk =1)=πk
Exercise 9.3
Figure 9.4
Graphical representation of a mixture model, in which the joint distribution is expressed in the form p(x,z) = p(z)p(x|z).
where the parameters {πk} must satisfy
0 à πk à 1
z
x
together with
p(x|z) =
K k=1
πk = 1
in order to be valid probabilities. Because z uses a 1-of-K representation, we can
also write this distribution in the form
p(z) =
9.2. Mixtures of Gaussians
431
K
k
k=1
Similarly, the conditional distribution of x given a particular value for z is a Gaussian
p(x|zk = 1) = N(x|μk,Σk) which can also be written in the form
then obtained by summing the joint distribution over all possible states of z to give K
p(x) = p(z)p(x|z) = πkN(x|μk,Σk) (9.12) z k=1
where we have made use of (9.10) and (9.11). Thus the marginal distribution of x is
a Gaussian mixture of the form (9.7). If we have several observations x1 , . . . , xN ,
K k=1
N (x|μk , Σk )zk . (9.11) The joint distribution is given by p(z)p(x|z), and the marginal distribution of x is
then, because we have represented the marginal distribution in the form p(x) =
p(x, z), it follows that for every observed data point xn there is a corresponding z
latent variable zn.
We have therefore found an equivalent formulation of the Gaussian mixture in-
volving an explicit latent variable. It might seem that we have not gained much by doing so. However, we are now able to work with the joint distribution p(x, z)
πzk . (9.10)
(9.8)
(9.9)
432 9. MIXTURE MODELS AND EM
instead of the marginal distribution p(x), and this will lead to significant simplifica- tions, most notably through the introduction of the expectation-maximization (EM) algorithm.
Another quantity that will play an important role is the conditional probability of z given x. We shall use γ(zk) to denote p(zk = 1|x), whose value can be found using Bayes’ theorem
γ(zk) ≡ p(zk = 1|x) =
p(zk = 1)p(x|zk = 1)
K
p(zj = 1)p(x|zj = 1)
j=1
= πkN(x|μk,Σk) .
K πjN(x|μj,Σj)
j=1
(9.13)
Section 8.1.2
We shall view πk as the prior probability of zk = 1, and the quantity γ(zk) as the corresponding posterior probability once we have observed x. As we shall see later, γ(zk) can also be viewed as the responsibility that component k takes for ‘explain- ing’ the observation x.
We can use the technique of ancestral sampling to generate random samples distributed according to the Gaussian mixture model. To do this, we first generate a value for z, which we denote z, from the marginal distribution p(z) and then generate a value for x from the conditional distribution p(x|z). Techniques for sampling from standard distributions are discussed in Chapter 11. We can depict samples from the joint distribution p(x,z) by plotting points at the corresponding values of x and then colouring them according to the value of z, in other words according to which Gaussian component was responsible for generating them, as shown in Figure 9.5(a). Similarly samples from the marginal distribution p(x) are obtained by taking the samples from the joint distribution and ignoring the values of z. These are illustrated in Figure 9.5(b) by plotting the x values without any coloured labels.
We can also use this synthetic data set to illustrate the ‘responsibilities’ by eval- uating, for every data point, the posterior probability for each component in the mixture distribution from which this data set was generated. In particular, we can represent the value of the responsibilities γ(znk) associated with data point xn by plotting the corresponding point using proportions of red, blue, and green ink given by γ(znk) for k = 1,2,3, respectively, as shown in Figure 9.5(c). So, for instance, a data point for which γ(zn1) = 1 will be coloured red, whereas one for which γ(zn2) = γ(zn3) = 0.5 will be coloured with equal proportions of blue and green ink and so will appear cyan. This should be compared with Figure 9.5(a) in which the data points were labelled using the true identity of the component from which they were generated.
9.2.1 Maximum likelihood
Suppose we have a data set of observations {x1 , . . . , xN }, and we wish to model this data using a mixture of Gaussians. We can represent this data set as an N × D
Figure 9.6
Graphical representation of a Gaussian mixture model
for a set of N i.i.d. data points {xn}, with corresponding π latent points {zn}, where n = 1,…,N.
zn
N lnp(X|π,μ,Σ) =
n=1
ln
K k=1
πkN(xn|μk,Σk) . (9.14)
Before discussing how to maximize this function, it is worth emphasizing that there is a significant problem associated with the maximum likelihood framework applied to Gaussian mixture models, due to the presence of singularities. For sim- plicity, consider a Gaussian mixture whose components have covariance matrices given by Σk = σk2I, where I is the unit matrix, although the conclusions will hold for general covariance matrices. Suppose that one of the components of the mixture model, let us say the jth component, has its mean μj exactly equal to one of the data
9.2. Mixtures of Gaussians 433
1
0.5
0
0 0.5 1
(a)
1
0.5
0
0 0.5 1
(b)
1
0.5
0
0 0.5 1
(c)
Figure 9.5 Example of 500 points drawn from the mixture of 3 Gaussians shown in Figure 2.23. (a) Samples from the joint distribution p(z)p(x|z) in which the three states of z, corresponding to the three components of the mixture, are depicted in red, green, and blue, and (b) the corresponding samples from the marginal distribution p(x), which is obtained by simply ignoring the values of z and just plotting the x values. The data set in (a) is said to be complete, whereas that in (b) is incomplete. (c) The same samples in which the colours represent the value of the responsibilities γ(znk) associated with data point xn, obtained by plotting the corresponding point using proportions of red, blue, and green ink given by γ(znk) for k = 1,2,3, respectively
matrix X in which the nth row is given by xTn . Similarly, the corresponding latent variables will be denoted by an N × K matrix Z with rows zTn . If we assume that the data points are drawn independently from the distribution, then we can express the Gaussian mixture model for this i.i.d. data set using the graphical representation shown in Figure 9.6. From (9.7) the log of the likelihood function is given by
xn μ
Σ
N
434
9. MIXTURE MODELS AND EM
Section 10.1
points so that μj = xn for some value of n. This data point will then contribute a term in the likelihood function of the form
N(xn|xn,σj2I) = 1 1 . (9.15) (2π)1/2 σj
If we consider the limit σj → 0, then we see that this term goes to infinity and so the log likelihood function will also go to infinity. Thus the maximization of the log likelihood function is not a well posed problem because such singularities will always be present and will occur whenever one of the Gaussian components ‘collapses’ onto a specific data point. Recall that this problem did not arise in the case of a single Gaussian distribution. To understand the difference, note that if a single Gaussian collapses onto a data point it will contribute multiplicative factors to the likelihood function arising from the other data points and these factors will go to zero exponentially fast, giving an overall likelihood that goes to zero rather than infinity. However, once we have (at least) two components in the mixture, one of the components can have a finite variance and therefore assign finite probability to all of the data points while the other component can shrink onto one specific data point and thereby contribute an ever increasing additive value to the log likelihood. This is illustrated in Figure 9.7. These singularities provide another example of the severe over-fitting that can occur in a maximum likelihood approach. We shall see that this difficulty does not occur if we adopt a Bayesian approach. For the moment, however, we simply note that in applying maximum likelihood to Gaussian mixture models we must take steps to avoid finding such pathological solutions and instead seek local maxima of the likelihood function that are well behaved. We can hope to avoid the singularities by using suitable heuristics, for instance by detecting when a Gaussian component is collapsing and resetting its mean to a randomly chosen value while also resetting its covariance to some large value, and then continuing with the optimization.
A further issue in finding maximum likelihood solutions arises from the fact that for any given maximum likelihood solution, a K-component mixture will have a total of K! equivalent solutions corresponding to the K! ways of assigning K sets of parameters to K components. In other words, for any given (nondegenerate) point in the space of parameter values there will be a further K!−1 additional points all of which give rise to exactly the same distribution. This problem is known as
Figure 9.7
Illustration of how singularities in the likelihood function arise with mixtures
of Gaussians. This should be com- p(x) pared with the case of a single Gaus-
sian shown in Figure 1.14 for which no singularities arise.
x
Section 10.1
identifiability (Casella and Berger, 2002) and is an important issue when we wish to interpret the parameter values discovered by a model. Identifiability will also arise when we discuss models having continuous latent variables in Chapter 12. However, for the purposes of finding a good density model, it is irrelevant because any of the equivalent solutions is as good as any other.
Maximizing the log likelihood function (9.14) for a Gaussian mixture model turns out to be a more complex problem than for the case of a single Gaussian. The difficulty arises from the presence of the summation over k that appears inside the logarithm in (9.14), so that the logarithm function no longer acts directly on the Gaussian. If we set the derivatives of the log likelihood to zero, we will no longer obtain a closed form solution, as we shall see shortly.
One approach is to apply gradient-based optimization techniques (Fletcher, 1987; Nocedal and Wright, 1999; Bishop and Nabney, 2008). Although gradient-based techniques are feasible, and indeed will play an important role when we discuss mixture density networks in Chapter 5, we now consider an alternative approach known as the EM algorithm which has broad applicability and which will lay the foundations for a discussion of variational inference techniques in Chapter 10.
9.2.2 EM for Gaussian mixtures
An elegant and powerful method for finding maximum likelihood solutions for models with latent variables is called the expectation-maximization algorithm, or EM algorithm (Dempster et al., 1977; McLachlan and Krishnan, 1997). Later we shall give a general treatment of EM, and we shall also show how EM can be generalized to obtain the variational inference framework. Initially, we shall motivate the EM algorithm by giving a relatively informal treatment in the context of the Gaussian mixture model. We emphasize, however, that EM has broad applicability, and indeed it will be encountered in the context of a variety of different models in this book.
Let us begin by writing down the conditions that must be satisfied at a maximum of the likelihood function. Setting the derivatives of ln p(X|π, μ, Σ) in (9.14) with respect to the means μk of the Gaussian components to zero, we obtain
N 0=−
n=1
π k N ( x n | μ k , Σ k )
πjN(xn|μj,Σj)Σk(xn −μk) (9.16)
j
μk = N Nk =
(9.17)
γ(znk)
where we have made use of the form (2.43) for the Gaussian distribution. Note that
the posterior probabilities, or responsibilities, given by (9.13) appear naturally on
the right-hand side. Multiplying by Σ−1 (which we assume to be nonsingular) and
rearranging we obtain
where we have defined
k
1 N
γ(znk)xn
k n=1 N
n=1
9.2. Mixtures of Gaussians 435
γ(znk). (9.18)
436 9. MIXTURE MODELS AND EM
Section 2.3.4
We can interpret Nk as the effective number of points assigned to cluster k. Note carefully the form of this solution. We see that the mean μk for the kth Gaussian component is obtained by taking a weighted mean of all of the points in the data set, in which the weighting factor for data point xn is given by the posterior probability γ(znk) that component k was responsible for generating xn.
If we set the derivative of ln p(X|π, μ, Σ) with respect to Σk to zero, and follow a similar line of reasoning, making use of the result for the maximum likelihood solution for the covariance matrix of a single Gaussian, we obtain
Appendix E
1 N
Σk = N
which has the same form as the corresponding result for a single Gaussian fitted to the data set, but again with each data point weighted by the corresponding poste- rior probability and with the denominator given by the effective number of points associated with the corresponding component.
Finally, we maximize ln p(X|π, μ, Σ) with respect to the mixing coefficients πk. Here we must take account of the constraint (9.9), which requires the mixing coefficients to sum to one. This can be achieved using a Lagrange multiplier and maximizing the following quantity
which gives
k n=1
γ(znk)(xn − μk)(xn − μk)T (9.19)
lnp(X|π,μ,Σ)+λ
K k=1
πk −1
(9.20)
(9.21)
N 0=
n=1
N ( x n | μ k , Σ k )
πjN(xn|μj,Σj)+λ
j
where again we see the appearance of the responsibilities. If we now multiply both sides by πk and sum over k making use of the constraint (9.9), we find λ = −N. Using this to eliminate λ and rearranging we obtain
πk = Nk (9.22) N
so that the mixing coefficient for the kth component is given by the average respon- sibility which that component takes for explaining the data points.
It is worth emphasizing that the results (9.17), (9.19), and (9.22) do not con- stitute a closed-form solution for the parameters of the mixture model because the responsibilities γ(znk) depend on those parameters in a complex way through (9.13). However, these results do suggest a simple iterative scheme for finding a solution to the maximum likelihood problem, which as we shall see turns out to be an instance of the EM algorithm for the particular case of the Gaussian mixture model. We first choose some initial values for the means, covariances, and mixing coefficients. Then we alternate between the following two updates that we shall call the E step
9.2. Mixtures of Gaussians 437
2
0
−2
−2 0 (a) 2
2
0
−2
−2 0 (b) 2
2
0
−2
−2 0 (c) 2
L=1
2
0
−2
−2 0 (d) 2
L=2
2
0
−2
−2 0 (e) 2
L=5
2
0
−2
−2 0 (f) 2
L = 20
Figure 9.8 Illustration of the EM algorithm using the Old Faithful set as used for the illustration of the K-means algorithm in Figure 9.1. See the text for details.
Section 9.4
and the M step, for reasons that will become apparent shortly. In the expectation step, or E step, we use the current values for the parameters to evaluate the posterior probabilities, or responsibilities, given by (9.13). We then use these probabilities in the maximization step, or M step, to re-estimate the means, covariances, and mix- ing coefficients using the results (9.17), (9.19), and (9.22). Note that in so doing we first evaluate the new means using (9.17) and then use these new values to find the covariances using (9.19), in keeping with the corresponding result for a single Gaussian distribution. We shall show that each update to the parameters resulting from an E step followed by an M step is guaranteed to increase the log likelihood function. In practice, the algorithm is deemed to have converged when the change in the log likelihood function, or alternatively in the parameters, falls below some threshold. We illustrate the EM algorithm for a mixture of two Gaussians applied to the rescaled Old Faithful data set in Figure 9.8. Here a mixture of two Gaussians is used, with centres initialized using the same values as for the K-means algorithm in Figure 9.1, and with precision matrices initialized to be proportional to the unit matrix. Plot (a) shows the data points in green, together with the initial configura- tion of the mixture model in which the one standard-deviation contours for the two
438 9. MIXTURE MODELS AND EM
Gaussian components are shown as blue and red circles. Plot (b) shows the result of the initial E step, in which each data point is depicted using a proportion of blue ink equal to the posterior probability of having been generated from the blue com- ponent, and a corresponding proportion of red ink given by the posterior probability of having been generated by the red component. Thus, points that have a significant probability for belonging to either cluster appear purple. The situation after the first M step is shown in plot (c), in which the mean of the blue Gaussian has moved to the mean of the data set, weighted by the probabilities of each data point belonging to the blue cluster, in other words it has moved to the centre of mass of the blue ink. Similarly, the covariance of the blue Gaussian is set equal to the covariance of the blue ink. Analogous results hold for the red component. Plots (d), (e), and (f) show the results after 2, 5, and 20 complete cycles of EM, respectively. In plot (f) the algorithm is close to convergence.
Note that the EM algorithm takes many more iterations to reach (approximate) convergence compared with the K-means algorithm, and that each cycle requires significantly more computation. It is therefore common to run the K-means algo- rithm in order to find a suitable initialization for a Gaussian mixture model that is subsequently adapted using EM. The covariance matrices can conveniently be ini- tialized to the sample covariances of the clusters found by the K-means algorithm, and the mixing coefficients can be set to the fractions of data points assigned to the respective clusters. As with gradient-based approaches for maximizing the log like- lihood, techniques must be employed to avoid singularities of the likelihood function in which a Gaussian component collapses onto a particular data point. It should be emphasized that there will generally be multiple local maxima of the log likelihood function, and that EM is not guaranteed to find the largest of these maxima. Because the EM algorithm for Gaussian mixtures plays such an important role, we summarize it below.
EM for Gaussian Mixtures
Given a Gaussian mixture model, the goal is to maximize the likelihood function with respect to the parameters (comprising the means and covariances of the components and the mixing coefficients).
1. Initialize the means μk, covariances Σk and mixing coefficients πk, and evaluate the initial value of the log likelihood.
2. E step. Evaluate the responsibilities using the current parameter values γ(znk) = πkN(xn|μk,Σk) . (9.23)
j=1
K πjN(xn|μj,Σj)
9.3. An Alternative View of EM 439 3. M step. Re-estimate the parameters using the current responsibilities
(9.24)
(9.25) (9.26)
(9.27)
where
πnew = Nk
Nk =
N n=1
new 1N
μk = N
γ(znk)xn
γ(z )(x −μnew)(x −μnew)T
n=1 kN
k n=1 1 N
Σnew =
k Nk nk n k n k
4. Evaluate the log likelihood
n=1
the convergence criterion is not satisfied return to step 2.
9.3. An Alternative View of EM
In this section, we present a complementary view of the EM algorithm that recog- nizes the key role played by latent variables. We discuss this approach first of all in an abstract setting, and then for illustration we consider once again the case of Gaussian mixtures.
The goal of the EM algorithm is to find maximum likelihood solutions for mod- els having latent variables. We denote the set of all observed data by X, in which the nth row represents xTn , and similarly we denote the set of all latent variables by Z, with a corresponding row zTn . The set of all model parameters is denoted by θ, and so the log likelihood function is given by
Z
Note that our discussion will apply equally well to continuous latent variables simply by replacing the sum over Z with an integral.
A key observation is that the summation over the latent variables appears inside the logarithm. Even if the joint distribution p(X, Z|θ) belongs to the exponential
N lnp(X|μ,Σ,π) =
K k=1
πkN(xn|μk,Σk)
and check for convergence of either the parameters or the log likelihood. If
ln
(9.28)
γ(znk).
ln p(X|θ) = ln
p(X, Z|θ) . (9.29)
440
9. MIXTURE MODELS AND EM
Section 9.4
family, the marginal distribution p(X|θ) typically does not as a result of this sum- mation. The presence of the sum prevents the logarithm from acting directly on the joint distribution, resulting in complicated expressions for the maximum likelihood solution.
Now suppose that, for each observation in X, we were told the corresponding value of the latent variable Z. We shall call {X, Z} the complete data set, and we shall refer to the actual observed data X as incomplete, as illustrated in Figure 9.5. The likelihood function for the complete data set simply takes the form ln p(X, Z|θ), and we shall suppose that maximization of this complete-data log likelihood function is straightforward.
In practice, however, we are not given the complete data set {X, Z}, but only the incomplete data X. Our state of knowledge of the values of the latent variables in Z is given only by the posterior distribution p(Z|X, θ). Because we cannot use the complete-data log likelihood, we consider instead its expected value under the posterior distribution of the latent variable, which corresponds (as we shall see) to the E step of the EM algorithm. In the subsequent M step, we maximize this expectation. If the current estimate for the parameters is denoted θold, then a pair of successive E and M steps gives rise to a revised estimate θnew. The algorithm is initialized by choosing some starting value for the parameters θ0. The use of the expectation may seem somewhat arbitrary. However, we shall see the motivation for this choice when we give a deeper treatment of EM in Section 9.4.
In the E step, we use the current parameter values θold to find the posterior distribution of the latent variables given by p(Z|X, θold). We then use this posterior distribution to find the expectation of the complete-data log likelihood evaluated for some general parameter value θ. This expectation, denoted Q(θ, θold), is given by
p(Z|X, θold) ln p(X, Z|θ). (9.30) In the M step, we determine the revised parameter estimate θnew by maximizing this
function
θnew =argmaxQ(θ,θold). (9.31) θ
Note that in the definition of Q(θ,θold), the logarithm acts directly on the joint distribution p(X, Z|θ), and so the corresponding M-step maximization will, by sup- position, be tractable.
The general EM algorithm is summarized below. It has the property, as we shall show later, that each cycle of EM will increase the incomplete-data log likelihood (unless it is already at a local maximum).
The General EM Algorithm
Given a joint distribution p(X, Z|θ) over observed variables X and latent vari- ables Z, governed by parameters θ, the goal is to maximize the likelihood func- tion p(X|θ) with respect to θ.
Q(θ, θold) =
Z
1. Choose an initial setting for the parameters θold.
Exercise 9.4
2. E step Evaluate p(Z|X, θold).
3. M step Evaluate θnew given by
where
9.3. An Alternative View of EM
441
θnew =argmaxQ(θ,θold) θ
(9.32)
Q(θ, θold) =
4. Check for convergence of either the log likelihood or the parameter values.
θold ← θnew (9.34)
and return to step 2.
The EM algorithm can also be used to find MAP (maximum posterior) solutions for models in which a prior p(θ) is defined over the parameters. In this case the E step remains the same as in the maximum likelihood case, whereas in the M step the quantity to be maximized is given by Q(θ, θold) + ln p(θ). Suitable choices for the prior will remove the singularities of the kind illustrated in Figure 9.7.
Here we have considered the use of the EM algorithm to maximize a likelihood function when there are discrete latent variables. However, it can also be applied when the unobserved variables correspond to missing values in the data set. The distribution of the observed values is obtained by taking the joint distribution of all the variables and then marginalizing over the missing ones. EM can then be used to maximize the corresponding likelihood function. We shall show an example of the application of this technique in the context of principal component analysis in Figure 12.11. This will be a valid procedure if the data values are missing at random, meaning that the mechanism causing values to be missing does not depend on the unobserved values. In many situations this will not be the case, for instance if a sensor fails to return a value whenever the quantity it is measuring exceeds some threshold.
9.3.1 Gaussian mixtures revisited
We now consider the application of this latent variable view of EM to the spe- cific case of a Gaussian mixture model. Recall that our goal is to maximize the log likelihood function (9.14), which is computed using the observed data set X, and we saw that this was more difficult than for the case of a single Gaussian distribution due to the presence of the summation over k that occurs inside the logarithm. Sup- pose then that in addition to the observed data set X, we were also given the values of the corresponding discrete variables Z. Recall that Figure 9.5(a) shows a ‘com- plete’ data set (i.e., one that includes labels showing which component generated each data point) while Figure 9.5(b) shows the corresponding ‘incomplete’ data set. The graphical model for the complete data is shown in Figure 9.9.
p(Z|X, θold) ln p(X, Z|θ). If the convergence criterion is not satisfied, then let
(9.33)
Z
442
9. MIXTURE MODELS AND EM
Figure 9.9
This shows the same graph as in Figure 9.6 except that
we now suppose that the discrete variables zn are ob- π served, as well as the data variables xn.
zn
xn μ
N
Now consider the problem of maximizing the likelihood for the complete data set {X, Z}. From (9.10) and (9.11), this likelihood function takes the form
(9.35)
(9.36)
p(X,Z|μ,Σ,π) =
πznkN(xn|μk,Σk)znk
N K
k
n=1 k=1
where znk denotes the kth component of zn. Taking the logarithm, we obtain
N K n=1 k=1
lnp(X,Z|μ,Σ,π)=
znk{lnπk +lnN(xn|μk,Σk)}.
Comparison with the log likelihood function (9.14) for the incomplete data shows that the summation over k and the logarithm have been interchanged. The loga- rithm now acts directly on the Gaussian distribution, which itself is a member of the exponential family. Not surprisingly, this leads to a much simpler solution to the maximum likelihood problem, as we now show. Consider first the maximization with respect to the means and covariances. Because zn is a K-dimensional vec- tor with all elements equal to 0 except for a single element having the value 1, the complete-data log likelihood function is simply a sum of K independent contribu- tions, one for each mixture component. Thus the maximization with respect to a mean or a covariance is exactly as for a single Gaussian, except that it involves only the subset of data points that are ‘assigned’ to that component. For the maximization with respect to the mixing coefficients, we note that these are coupled for different values of k by virtue of the summation constraint (9.9). Again, this can be enforced using a Lagrange multiplier as before, and leads to the result
1 N
πk = N
so that the mixing coefficients are equal to the fractions of data points assigned to the corresponding components.
Thus we see that the complete-data log likelihood function can be maximized trivially in closed form. In practice, however, we do not have values for the latent variables so, as discussed earlier, we consider the expectation, with respect to the posterior distribution of the latent variables, of the complete-data log likelihood.
n=1
znk (9.37)
Σ
Exercise 9.5 Section 8.2
Exercise 9.8
9.3. An Alternative View of EM 443 Using (9.10) and (9.11) together with Bayes’ theorem, we see that this posterior
distribution takes the form
p(Z|X,μ,Σ,π) ∝
and hence factorizes over n so that under the posterior distribution the {zn} are independent. This is easily verified by inspection of the directed graph in Figure 9.6 and making use of the d-separation criterion. The expected value of the indicator variable znk under this posterior distribution is then given by
znk znk [πkN(xn|μk,Σk)]
EZ[lnp(X,Z|μ,Σ,π)]=
γ(znk){lnπk +lnN(xn|μk,Σk)}. (9.40)
znk E[znk]=
znj
=
πkN(xn|μk,Σk) = γ(znk) (9.39)
N K n=1 k=1
[πkN(xn|μk,Σk)]znk . (9.38)
πjN(xn|μj,Σj) znj
K πjN(xn|μj,Σj)
j=1
which is just the responsibility of component k for data point xn. The expected value
of the complete-data log likelihood function is therefore given by
N K n=1 k=1
We can now proceed as follows. First we choose some initial values for the param- eters μold, Σold and πold, and use these to evaluate the responsibilities (the E step). We then keep the responsibilities fixed and maximize (9.40) with respect to μk, Σk and πk (the M step). This leads to closed form solutions for μnew, Σnew and πnew given by (9.17), (9.19), and (9.22) as before. This is precisely the EM algorithm for Gaussian mixtures as derived earlier. We shall gain more insight into the role of the expected complete-data log likelihood function when we give a proof of convergence of the EM algorithm in Section 9.4.
9.3.2 Relation to K-means
Comparison of the K-means algorithm with the EM algorithm for Gaussian mixtures shows that there is a close similarity. Whereas the K-means algorithm performs a hard assignment of data points to clusters, in which each data point is associated uniquely with one cluster, the EM algorithm makes a soft assignment based on the posterior probabilities. In fact, we can derive the K-means algorithm as a particular limit of EM for Gaussian mixtures as follows.
Consider a Gaussian mixture model in which the covariance matrices of the mixture components are given by I, where is a variance parameter that is shared
444 9. MIXTURE MODELS AND EM
by all of the components, and I is the identity matrix, so that
Exercise 9.11
p(x|μk,Σk)= 1 exp −1 x−μk 2 . (9.41) (2π)1/2 2
We now consider the EM algorithm for a mixture of K Gaussians of this form in which we treat as a fixed constant, instead of a parameter to be re-estimated. From (9.13) the posterior probabilities, or responsibilities, for a particular data point xn, are given by
πkexp{−xn−μk 2/2}
γ(znk)=πjexp −xn−μj 2/2 . (9.42)
j
If we consider the limit → 0, we see that in the denominator the term for which xn − μj 2 is smallest will go to zero most slowly, and hence the responsibilities γ(znk) for the data point xn all go to zero except for term j, for which the responsi- bility γ(znj) will go to unity. Note that this holds independently of the values of the πk so long as none of the πk is zero. Thus, in this limit, we obtain a hard assignment of data points to clusters, just as in the K-means algorithm, so that γ(znk) → rnk where rnk is defined by (9.2). Each data point is thereby assigned to the cluster
having the closest mean.
The EM re-estimation equation for the μk, given by (9.17), then reduces to the
K-means result (9.4). Note that the re-estimation formula for the mixing coefficients (9.22) simply re-sets the value of πk to be equal to the fraction of data points assigned to cluster k, although these parameters no longer play an active role in the algorithm.
Finally, in the limit → 0 the expected complete-data log likelihood, given by (9.40), becomes
Section 13.2
1 N K
EZ[lnp(X,Z|μ,Σ,π)]→−2
Thus we see that in this limit, maximizing the expected complete-data log likelihood is equivalent to minimizing the distortion measure J for the K-means algorithm given by (9.1).
Note that the K-means algorithm does not estimate the covariances of the clus- ters but only the cluster means. A hard-assignment version of the Gaussian mixture model with general covariance matrices, known as the elliptical K-means algorithm, has been considered by Sung and Poggio (1994).
9.3.3 Mixtures of Bernoulli distributions
So far in this chapter, we have focussed on distributions over continuous vari- ables described by mixtures of Gaussians. As a further example of mixture mod- elling, and to illustrate the EM algorithm in a different context, we now discuss mix- tures of discrete binary variables described by Bernoulli distributions. This model is also known as latent class analysis (Lazarsfeld and Henry, 1968; McLachlan and Peel, 2000). As well as being of practical importance in its own right, our discus- sion of Bernoulli mixtures will also lay the foundation for a consideration of hidden Markov models over discrete variables.
n=1 k=1
rnk xn −μk 2 +const. (9.43)
Exercise 9.12
9.3. An Alternative View of EM 445 Consider a set of D binary variables xi, where i = 1,…,D, each of which is
governed by a Bernoulli distribution with parameter μi, so that
i=1
where x = (x1,…,xD)T and μ = (μ1,…,μD)T. We see that the individual variables xi are independent, given μ. The mean and covariance of this distribution are easily seen to be
p(x|μ) =
μxi (1 − μi)(1−xi) (9.44)
E[x] = μ
cov[x] = diag{μi(1 − μi)}.
Now let us consider a finite mixture of these distributions given by
K
πkp(x|μk) where μ = {μ1,…,μK}, π = {π1,…,πK}, and
i=1
K k=1
K k=1
(9.45) (9.46)
(9.47)
(9.48)
(9.49) (9.50)
p(x|μ, π) =
)(1−xi).
The mean and covariance of this mixture distribution are given by
E[x] = cov[x] =
−E[x]E[x]T
p(x|μ ) =
D
i
D
k ki ki
k=1
μxi (1 − μ
πkμk
πk Σk +μkμTk
= diag {μki(1 − μki)}. Because the covariance matrix cov[x] is no longer diagonal, the mixture distribution can capture correlations between the vari- ables, unlike a single Bernoulli distribution.
where Σk
If we are given a data set X = {x1,…,xN} then the log likelihood function
for this model is given by
N ln p(X|μ, π) =
n=1
ln
K k=1
πkp(xn|μk) . (9.51)
Again we see the appearance of the summation inside the logarithm, so that the maximum likelihood solution no longer has closed form.
We now derive the EM algorithm for maximizing the likelihood function for the mixture of Bernoulli distributions. To do this, we first introduce an explicit latent
446 9. MIXTURE MODELS AND EM
Exercise 9.14
If we form the product of p(x|z, μ) and p(z|π) and then marginalize over z, then we recover (9.47).
In order to derive the EM algorithm, we first write down the complete-data log likelihood function, which is given by
variable z associated with each instance of x. As in the case of the Gaussian mixture, z = (z1,…,zK)T is a binary K-dimensional variable having a single component equal to 1, with all other components equal to 0. We can then write the conditional distribution of x, given the latent variable, as
ln p(X, Z|μ, π) =
N K
+
(9.54)
D
K k=1
p(x|μk )zk (9.52)
p(x|z, μ) =
while the prior distribution for the latent variables is the same as for the mixture of
Gaussians model, so that
p(z|π) =
πzk . (9.53)
K
k
k=1
znk ln πk
[xni ln μki + (1 − xni) ln(1 − μki)]
n=1 k=1
i=1
where X = {xn} and Z = {zn}. Next we take the expectation of the complete-data
log likelihood with respect to the posterior distribution of the latent variables to give
EZ[ln p(X, Z|μ, π)] = D
+
N K
i=1
γ(znk) ln πk [xni ln μki + (1 − xni) ln(1 − μki)]
n=1 k=1
where γ(znk) = E[znk] is the posterior probability, or responsibility, of component k given data point xn. In the E step, these responsibilities are evaluated using Bayes’
theorem, which takes the form
γ(znk) = E[znk]
znk znk [πkp(xn|μk)]
znk
= znj
(9.55)
πj p(xn |μj ) znj
= πk p(xn |μk ) .
K
πj p(xn |μj )
j=1
(9.56)
Exercise 9.15
Exercise 9.16
where Nk is the effective number of data points associated with component k. In the M step, we maximize the expected complete-data log likelihood with respect to the parameters μk and π. If we set the derivative of (9.55) with respect to μk equal to zero and rearrange the terms, we obtain
μk = xk. (9.59)
We see that this sets the mean of component k equal to a weighted mean of the
data, with weighting coefficients given by the responsibilities that component k takes
Exercise 9.17
Section 9.4
steps to those used for the mixture of Gaussians, we then obtain
πk = Nk (9.60)
N
which represents the intuitively reasonable result that the mixing coefficient for com- ponent k is given by the effective fraction of points in the data set explained by that component.
Note that in contrast to the mixture of Gaussians, there are no singularities in
which the likelihood function goes to infinity. This can be seen by noting that the
likelihood function is bounded above because 0 à p(xn|μk) à 1. There exist
singularities at which the likelihood function goes to zero, but these will not be
found by EM provided it is not initialized to a pathological starting point, because
the EM algorithm always increases the value of the likelihood function, until a local
maximum is found. We illustrate the Bernoulli mixture model in Figure 9.10 by
using it to model handwritten digits. Here the digit images have been turned into
binary vectors by setting all elements whose values exceed 0.5 to 1 and setting the
remaining elements to 0. We now fit a data set of N = 600 such digits, comprising
the digits ‘2’, ‘3’, and ‘4’, with a mixture of K = 3 Bernoulli distributions by
running 10 iterations of the EM algorithm. The mixing coefficients were initialized
the data set corresponding to the different digits.
The conjugate prior for the parameters of a Bernoulli distribution is given by
the beta distribution, and we have seen that a beta prior is equivalent to introducing
9.3. An Alternative View of EM 447 If we consider the sum over n in (9.55), we see that the responsibilities enter
only through two terms, which can be written as
N
Nk = γ(znk)
(9.57)
(9.58)
n=1
1 N
xk = N
γ(znk)xn
k n=1
for data points. For the maximization with respect to π , we need to introduce a
k
Lagrange multiplier to enforce the constraint k πk = 1. Following analogous
to π = 1/K, and the parameters μ were set to random values chosen uniformly in
k kj
j μkj = 1. We see that a mixture of 3 Bernoulli distributions is able to find the three clusters in
the range (0.25, 0.75) and then normalized to satisfy the constraint that
448 9. MIXTURE MODELS AND EM
Figure 9.10 Illustration of the Bernoulli mixture model in which the top row shows examples from the digits data set after converting the pixel values from grey scale to binary using a threshold of 0.5. On the bottom row the first three images show the parameters μki for each of the three components in the mixture model. As a comparison, we also fit the same data set using a single multivariate Bernoulli distribution, again using maximum likelihood. This amounts to simply averaging the counts in each pixel and is shown by the right-most image on the bottom row.
Section 2.1.1 Exercise 9.18 Exercise 9.19
additional effective observations of x. We can similarly introduce priors into the Bernoulli mixture model, and use EM to maximize the posterior probability distri- butions.
It is straightforward to extend the analysis of Bernoulli mixtures to the case of multinomial binary variables having M > 2 states by making use of the discrete dis- tribution (2.26). Again, we can introduce Dirichlet priors over the model parameters if desired.
9.3.4 EM for Bayesian linear regression
As a third example of the application of EM, we return to the evidence ap- proximation for Bayesian linear regression. In Section 3.5.2, we obtained the re- estimation equations for the hyperparameters α and β by evaluation of the evidence and then setting the derivatives of the resulting expression to zero. We now turn to an alternative approach for finding α and β based on the EM algorithm. Recall that our goal is to maximize the evidence function p(t|α, β) given by (3.77) with respect to α and β. Because the parameter vector w is marginalized out, we can regard it as a latent variable, and hence we can optimize this marginal likelihood function using EM. In the E step, we compute the posterior distribution of w given the current set- ting of the parameters α and β and then use this to find the expected complete-data log likelihood. In the M step, we maximize this quantity with respect to α and β. We have already derived the posterior distribution of w because this is given by (3.49). The complete-data log likelihood function is then given by
ln p(t, w|α, β) = ln p(t|w, β) + ln p(w|α) (9.61)
Exercise 9.20
Exercise 9.21
equation
9.3. An Alternative View of EM 449
where the likelihood p(t|w, β) and the prior p(w|α) are given by (3.10) and (3.52), respectively, and y(x, w) is given by (3.3). Taking the expectation with respect to the posterior distribution of w then gives
E[lnp(t,w|α,β)] = Mln α −αEwTw+Nln β 2 2π 2 2 2π
β N
E (tn−wTφn)2 .
Setting the derivatives with respect to α to zero, we obtain the M step re-estimation
−2
α= M = M . (9.63)
An analogous result holds for β.
Note that this re-estimation equation takes a slightly different form from the
corresponding result (3.92) derived by direct evaluation of the evidence function. However, they each involve computation and inversion (or eigen decomposition) of an M × M matrix and hence will have comparable computational cost per iteration.
These two approaches to determining α should of course converge to the same result (assuming they find the same local maximum of the evidence function). This can be verified by first noting that the quantity γ is defined by
M 1
λi + α = M − αTr(SN ). (9.64) At a stationary point of the evidence function, the re-estimation equation (3.92) will
be self-consistently satisfied, and hence we can substitute for γ to give
αmTNmN =γ=M−αTr(SN) (9.65)
and solving for α we obtain (9.63), which is precisely the EM re-estimation equation. As a final example, we consider a closely related model, namely the relevance vector machine for regression discussed in Section 7.2.1. There we used direct max- imization of the marginal likelihood to derive re-estimation equations for the hyper- parameters α and β. Here we consider an alternative approach in which we view the weight vector w as a latent variable and apply the EM algorithm. The E step involves finding the posterior distribution over the weights, and this is given by (7.81). In the M step we maximize the expected complete-data log likelihood, which is defined by
Ew [ln p(t|X, w, β)p(w|α)] (9.66)
where the expectation is taken with respect to the posterior distribution computed using the ‘old’ parameter values. To compute the new parameter values we maximize with respect to α and β to give
n=1
(9.62)
Exercise 9.22
γ = M − α
E[wTw] mTNmN +Tr(SN)
i=1
450 9. MIXTURE MODELS AND EM
αnew =
1
(9.67)
i
m2i +Σii
t − ΦmN
(βnew)−1 =
These re-estimation equations are formally equivalent to those obtained by direct
maxmization.
2 + β−1
i γi (9.68)
Exercise 9.23
N
Section 10.1
9.4. The EM Algorithm in General
The expectation maximization algorithm, or EM algorithm, is a general technique for finding maximum likelihood solutions for probabilistic models having latent vari- ables (Dempster et al., 1977; McLachlan and Krishnan, 1997). Here we give a very general treatment of the EM algorithm and in the process provide a proof that the EM algorithm derived heuristically in Sections 9.2 and 9.3 for Gaussian mixtures does indeed maximize the likelihood function (Csisza`r and Tusna`dy, 1984; Hath- away, 1986; Neal and Hinton, 1999). Our discussion will also form the basis for the derivation of the variational inference framework.
Consider a probabilistic model in which we collectively denote all of the ob- served variables by X and all of the hidden variables by Z. The joint distribution p(X, Z|θ) is governed by a set of parameters denoted θ. Our goal is to maximize the likelihood function that is given by
Z
Here we are assuming Z is discrete, although the discussion is identical if Z com- prises continuous variables or a combination of discrete and continuous variables, with summation replaced by integration as appropriate.
We shall suppose that direct optimization of p(X|θ) is difficult, but that opti- mization of the complete-data likelihood function p(X, Z|θ) is significantly easier. Next we introduce a distribution q(Z) defined over the latent variables, and we ob- serve that, for any choice of q(Z), the following decomposition holds
where we have defined
L(q, θ) = KL(q p) = −
q(Z) ln q(Z)
p(Z|X, θ)
p(X|θ) =
p(X, Z|θ). (9.69)
lnp(X|θ)=L(q,θ)+KL(q p)
(9.70)
(9.71)
(9.72)
p(X, Z|θ)
Z
Z
q(Z)ln q(Z) .
Note that L(q, θ) is a functional (see Appendix D for a discussion of functionals) of the distribution q(Z), and a function of the parameters θ. It is worth studying
Figure 9.11
Illustration of the decomposition given by (9.70), which holds for any choice of distribution q(Z). Because the Kullback-Leibler divergence satisfies KL(q p) 0, we see that the quan- tity L(q, θ) is a lower bound on the log likelihood function ln p(X|θ).
9.4. The EM Algorithm in General
451
KL(q||p)
L(q, θ)
ln p(X|θ)
Exercise 9.24
Section 1.6.1
carefully the forms of the expressions (9.71) and (9.72), and in particular noting that they differ in sign and also that L(q, θ) contains the joint distribution of X and Z while KL(q p) contains the conditional distribution of Z given X. To verify the decomposition (9.70), we first make use of the product rule of probability to give
ln p(X, Z|θ) = ln p(Z|X, θ) + ln p(X|θ) (9.73)
which we then substitute into the expression for L(q, θ). This gives rise to two terms, one of which cancels KL(q p) while the other gives the required log likelihood ln p(X|θ) after noting that q(Z) is a normalized distribution that sums to 1.
From (9.72), we see that KL(q p) is the Kullback-Leibler divergence between q(Z) and the posterior distribution p(Z|X, θ). Recall that the Kullback-Leibler di- vergence satisfies KL(q p) 0, with equality if, and only if, q(Z) = p(Z|X, θ). It therefore follows from (9.70) that L(q, θ) à ln p(X|θ), in other words that L(q, θ) is a lower bound on lnp(X|θ). The decomposition (9.70) is illustrated in Fig- ure 9.11.
The EM algorithm is a two-stage iterative optimization technique for finding maximum likelihood solutions. We can use the decomposition (9.70) to define the EM algorithm and to demonstrate that it does indeed maximize the log likelihood. Suppose that the current value of the parameter vector is θold. In the E step, the lower bound L(q, θold) is maximized with respect to q(Z) while holding θold fixed. The solution to this maximization problem is easily seen by noting that the value of ln p(X|θold) does not depend on q(Z) and so the largest value of L(q, θold) will occur when the Kullback-Leibler divergence vanishes, in other words when q(Z) is equal to the posterior distribution p(Z|X,θold). In this case, the lower bound will equal the log likelihood, as illustrated in Figure 9.12.
In the subsequent M step, the distribution q(Z) is held fixed and the lower bound L(q,θ) is maximized with respect to θ to give some new value θnew. This will cause the lower bound L to increase (unless it is already at a maximum), which will necessarily cause the corresponding log likelihood function to increase. Because the distribution q is determined using the old parameter values rather than the new values and is held fixed during the M step, it will not equal the new posterior distribution p(Z|X,θnew), and hence there will be a nonzero KL divergence. The increase in the log likelihood function is therefore greater than the increase in the lower bound, as
452
9. MIXTURE MODELS AND EM
Figure 9.12
Illustration of the E step of the EM algorithm. The q distribution is set equal to the posterior distribution for the current parameter val- ues θold, causing the lower bound to move up to the same value as the log like- lihood function, with the KL divergence vanishing.
KL(q||p) = 0
L(q, θold)
ln p(X|θold)
shown in Figure 9.13. If we substitute q(Z) = p(Z|X, θold) into (9.71), we see that, after the E step, the lower bound takes the form
L(q, θ) = p(Z|X, θold) ln p(X, Z|θ) − p(Z|X, θold) ln p(Z|X, θold) ZZ
= Q(θ, θold) + const (9.74)
where the constant is simply the negative entropy of the q distribution and is there- fore independent of θ. Thus in the M step, the quantity that is being maximized is the expectation of the complete-data log likelihood, as we saw earlier in the case of mix- tures of Gaussians. Note that the variable θ over which we are optimizing appears only inside the logarithm. If the joint distribution p(Z, X|θ) comprises a member of the exponential family, or a product of such members, then we see that the logarithm will cancel the exponential and lead to an M step that will be typically much simpler than the maximization of the corresponding incomplete-data log likelihood function p(X|θ).
The operation of the EM algorithm can also be viewed in the space of parame- ters, as illustrated schematically in Figure 9.14. Here the red curve depicts the (in-
Figure 9.13
Illustration of the M step of the EM
algorithm. The distribution q(Z) KL(q||p) is held fixed and the lower bound
L(q,θ) is maximized with respect
to the parameter vector θ to give
a revised value θnew. Because the
KL divergence is nonnegative, this
causes the log likelihood lnp(X|θ)
to increase by at least as much as
the lower bound does.
L(q, θnew)
ln p(X|θnew)
n
that is evaluated in the E step takes the form
9.4. The EM Algorithm in General 453
lnp(X|θ)
L(q,θ)
θold θnew
Exercise 9.25
complete data) log likelihood function whose value we wish to maximize. We start with some initial parameter value θold, and in the first E step we evaluate the poste- rior distribution over latent variables, which gives rise to a lower bound L(θ, θ(old)) whose value equals the log likelihood at θ(old), as shown by the blue curve. Note that the bound makes a tangential contact with the log likelihood at θ(old), so that both curves have the same gradient. This bound is a convex function having a unique maximum (for mixture components from the exponential family). In the M step, the bound is maximized giving the value θ(new), which gives a larger value of log likeli- hood than θ(old). The subsequent E step then constructs a bound that is tangential at θ(new) as shown by the green curve.
For the particular case of an independent, identically distributed data set, X
will comprise N data points {xn} while Z will comprise N corresponding latent
variables {z }, where n = 1, . . . , N . From the independence assumption, we have n
Figure 9.14
The EM algorithm involves alter- nately computing a lower bound on the log likelihood for the cur- rent parameter values and then maximizing this bound to obtain the new parameter values. See the text for a full discussion.
p(X, Z) = p(x , z ) and, by marginalizing over the {z } we have p(X) =
nnn
n p(xn). Using the sum and product rules, we see that the posterior probability
p(Z|X, θ) =
p(X, Z|θ)
Z
= N Z n=1
=
p(zn|xn, θ)
(9.75)
p(X, Z|θ)
N n=1
p(xn, zn|θ)
p(xn, zn|θ)
N n=1
and so the posterior distribution also factorizes with respect to n. In the case of the Gaussian mixture model this simply says that the responsibility that each of the mixture components takes for a particular data point xn depends only on the value of xn and on the parameters θ of the mixture components, not on the values of the other data points.
We have seen that both the E and the M steps of the EM algorithm are increas- ing the value of a well-defined bound on the log likelihood function and that the
454 9. MIXTURE MODELS AND EM
complete EM cycle will change the model parameters in such a way as to cause the log likelihood to increase (unless it is already at a maximum, in which case the parameters remain unchanged).
We can also use the EM algorithm to maximize the posterior distribution p(θ|X) for models in which we have introduced a prior p(θ) over the parameters. To see this, we note that as a function of θ, we have p(θ|X) = p(θ, X)/p(X) and so
ln p(θ|X) = ln p(θ, X) − ln p(X). (9.76) Making use of the decomposition (9.70), we have
lnp(θ|X) = L(q,θ)+KL(q p)+lnp(θ)−lnp(X)
L(q, θ) + ln p(θ) − ln p(X). (9.77)
where ln p(X) is a constant. We can again optimize the right-hand side alternately with respect to q and θ. The optimization with respect to q gives rise to the same E- step equations as for the standard EM algorithm, because q only appears in L(q, θ). The M-step equations are modified through the introduction of the prior term ln p(θ), which typically requires only a small modification to the standard maximum likeli- hood M-step equations.
The EM algorithm breaks down the potentially difficult problem of maximizing the likelihood function into two stages, the E step and the M step, each of which will often prove simpler to implement. Nevertheless, for complex models it may be the case that either the E step or the M step, or indeed both, remain intractable. This leads to two possible extensions of the EM algorithm, as follows.
The generalized EM, or GEM, algorithm addresses the problem of an intractable M step. Instead of aiming to maximize L(q, θ) with respect to θ, it seeks instead to change the parameters in such a way as to increase its value. Again, because L(q, θ) is a lower bound on the log likelihood function, each complete EM cycle of the GEM algorithm is guaranteed to increase the value of the log likelihood (unless the parameters already correspond to a local maximum). One way to exploit the GEM approach would be to use one of the nonlinear optimization strategies, such as the conjugate gradients algorithm, during the M step. Another form of GEM algorithm, known as the expectation conditional maximization, or ECM, algorithm, involves making several constrained optimizations within each M step (Meng and Rubin, 1993). For instance, the parameters might be partitioned into groups, and the M step is broken down into multiple steps each of which involves optimizing one of the subset with the remainder held fixed.
We can similarly generalize the E step of the EM algorithm by performing a partial, rather than complete, optimization of L(q, θ) with respect to q(Z) (Neal and Hinton, 1999). As we have seen, for any given value of θ there is a unique maximum of L(q, θ) with respect to q(Z) that corresponds to the posterior distribution qθ (Z) = p(Z|X,θ) and that for this choice of q(Z) the bound L(q,θ) is equal to the log likelihood function lnp(X|θ). It follows that any algorithm that converges to the global maximum of L(q,θ) will find a value of θ that is also a global maximum of the log likelihood ln p(X|θ). Provided p(X, Z|θ) is a continuous function of θ
Exercise 9.26
then, by continuity, any local maximum of L(q, θ) will also be a local maximum of ln p(X|θ).
Consider the case of N independent data points x1 , . . . , xN with corresponding latent variables z1, . . . , zN . The joint distribution p(X, Z|θ) factorizes over the data points, and this structure can be exploited in an incremental form of EM in which at each EM cycle only one data point is processed at a time. In the E step, instead of recomputing the responsibilities for all of the data points, we just re-evaluate the responsibilities for one data point. It might appear that the subsequent M step would require computation involving the responsibilities for all of the data points. How- ever, if the mixture components are members of the exponential family, then the responsibilities enter only through simple sufficient statistics, and these can be up- dated efficiently. Consider, for instance, the case of a Gaussian mixture, and suppose we perform an update for data point m in which the corresponding old and new values of the responsibilities are denoted γold(zmk) and γnew(zmk). In the M step, the required sufficient statistics can be updated incrementally. For instance, for the means the sufficient statistics are defined by (9.17) and (9.18) from which we obtain
μnew = μold + γnew(zmk) − γold(zmk) x − μold kk Nnew mk
(9.78)
(9.79)
Exercises 455
Exercises
9.1
9.2
The corresponding results for the covariances and the mixing coefficients are analo- gous.
Thus both the E step and the M step take fixed time that is independent of the total number of data points. Because the parameters are revised after each data point, rather than waiting until after the whole data set is processed, this incremental ver- sion can converge faster than the batch version. Each E or M step in this incremental algorithm is increasing the value of L(q,θ) and, as we have shown above, if the algorithm converges to a local (or global) maximum of L(q, θ), this will correspond to a local (or global) maximum of the log likelihood function ln p(X|θ).
( ) www Consider the K-means algorithm discussed in Section 9.1. Show that as a consequence of there being a finite number of possible assignments for the set of discrete indicator variables rnk, and that for each such assignment there is a unique optimum for the {μk}, the K-means algorithm must converge after a finite number of iterations.
( ) Apply the Robbins-Monro sequential estimation procedure described in Sec- tion 2.3.5 to the problem of finding the roots of the regression function given by the derivatives of J in (9.1) with respect to μk. Show that this leads to a stochastic K-means algorithm in which, for each data point xn, the nearest prototype μk is updated using (9.5).
together with
k
Nnew = Nold + γnew(z ) − γold(z ). k k mk mk
456 9. MIXTURE MODELS AND EM
9.3 ( ) www Consider a Gaussian mixture model in which the marginal distribution p(z) for the latent variable is given by (9.10), and the conditional distribution p(x|z) for the observed variable is given by (9.11). Show that the marginal distribution p(x), obtained by summing p(z)p(x|z) over all possible values of z, is a Gaussian mixture of the form (9.7).
9.4 ( ) Suppose we wish to use the EM algorithm to maximize the posterior distri- bution over parameters p(θ|X) for a model containing latent variables, where X is the observed data set. Show that the E step remains the same as in the maximum likelihood case, whereas in the M step the quantity to be maximized is given by Q(θ, θold) + ln p(θ) where Q(θ, θold) is defined by (9.30).
9.5 ( ) Consider the directed graph for a Gaussian mixture model shown in Figure 9.6. By making use of the d-separation criterion discussed in Section 8.2, show that the posterior distribution of the latent variables factorizes with respect to the different data points so that
N
n=1
9.6 ( ) Consider a special case of a Gaussian mixture model in which the covari- ance matrices Σk of the components are all constrained to have a common value Σ. Derive the EM equations for maximizing the likelihood function under such a model.
9.7 ( ) www Verify that maximization of the complete-data log likelihood (9.36) for a Gaussian mixture model leads to the result that the means and covariances of each component are fitted independently to the corresponding group of data points, and the mixing coefficients are given by the fractions of points in each group.
9.8 ( ) www Show that if we maximize (9.40) with respect to μk while keeping the responsibilities γ(znk) fixed, we obtain the closed form solution given by (9.17).
9.9 ( ) Show that if we maximize (9.40) with respect to Σk and πk while keeping the responsibilities γ(znk) fixed, we obtain the closed form solutions given by (9.19) and (9.22).
9.10 ( ) Consider a density model given by a mixture distribution
K k=1
and suppose that we partition the vector x into two parts so that x = (xa,xb). Show that the conditional density p(xb|xa) is itself a mixture distribution and find expressions for the mixing coefficients and for the component densities.
p(Z|X, μ, Σ, π) =
p(zn|xn, μ, Σ, π). (9.80)
p(x) =
πkp(x|k) (9.81)
9.11 ( ) In Section 9.3.2, we obtained a relationship between K means and EM for Gaussian mixtures by considering a mixture model in which all components have covariance I. Show that in the limit → 0, maximizing the expected complete- data log likelihood for this model, given by (9.40), is equivalent to minimizing the distortion measure J for the K-means algorithm given by (9.1).
9.12 ( ) www Consider a mixture distribution of the form
K k=1
where the elements of x could be discrete or continuous or a combination of these. Denote the mean and covariance of p(x|k) by μk and Σk, respectively. Show that the mean and covariance of the mixture distribution are given by (9.49) and (9.50).
9.13 ( ) Using the re-estimation equations for the EM algorithm, show that a mix- ture of Bernoulli distributions, with its parameters set to values corresponding to a maximum of the likelihood function, has the property that
1 N
E[x] = N
Hence show that if the parameters of this model are initialized such that all compo- nents have the same mean μk = μ for k = 1,…,K, then the EM algorithm will converge after one iteration, for any choice of the initial mixing coefficients, and that this solution has the property μk = x. Note that this represents a degenerate case of the mixture model in which all of the components are identical, and in practice we try to avoid such solutions by using an appropriate initialization.
9.14 ( ) Consider the joint distribution of latent and observed variables for the Bernoulli distribution obtained by forming the product of p(x|z, μ) given by (9.52) and p(z|π) given by (9.53). Show that if we marginalize this joint distribution with respect to z, then we obtain (9.47).
9.15 ( ) www Show that if we maximize the expected complete-data log likelihood function (9.55) for a mixture of Bernoulli distributions with respect to μk , we obtain the M step equation (9.59).
9.16 ( ) Show that if we maximize the expected complete-data log likelihood function (9.55) for a mixture of Bernoulli distributions with respect to the mixing coefficients πk, using a Lagrange multiplier to enforce the summation constraint, we obtain the M step equation (9.60).
9.17 ( ) www Show that as a consequence of the constraint 0 à p(xn|μk) à 1 for the discrete variable xn, the incomplete-data log likelihood function for a mixture of Bernoulli distributions is bounded above, and hence that there are no singularities for which the likelihood goes to infinity.
n=1
Exercises 457
p(x) =
πkp(x|k) (9.82)
xn ≡ x. (9.83)
458 9. MIXTURE MODELS AND EM
9.18 ( ) Consider a Bernoulli mixture model as discussed in Section 9.3.3, together with a prior distribution p(μk|ak,bk) over each of the parameter vectors μk given by the beta distribution (2.13), and a Dirichlet prior p(π|α) given by (2.38). Derive the EM algorithm for maximizing the posterior probability p(μ, π|X).
9.19 ( ) Consider a D-dimensional variable x each of whose components i is itself a multinomial variable of degree M so that x is a binary vector with components x
ij where i = 1,…,D and j = 1,…,M, subject to the constraint that j xij = 1 for all i. Suppose that the distribution of these variables is described by a mixture of the
discrete multinomial distributions considered in Section 2.2 so that
K k=1
D M i=1 j=1
where
p(x) =
πkp(x|μk)
(9.84)
p(x|μ ) = k
μxij . kij
(9.85) = 1|μ ) and must satisfy
The parameters μ
kij ijk
represent the probabilities p(x
0 à μkij à 1 together with the constraint j μkij = 1 for all values of k and i.
Given an observed data set {xn}, where n = 1,…,N, derive the E and M step equations of the EM algorithm for optimizing the mixing coefficients πk and the component parameters μkij of this distribution by maximum likelihood.
9.20 ( ) www Show that maximization of the expected complete-data log likelihood function (9.62) for the Bayesian linear regression model leads to the M step re- estimation result (9.63) for α.
9.21 ( ) Using the evidence framework of Section 3.5, derive the M-step re-estimation equations for the parameter β in the Bayesian linear regression model, analogous to the result (9.63) for α.
9.22 ( ) By maximization of the expected complete-data log likelihood defined by (9.66), derive the M step equations (9.67) and (9.68) for re-estimating the hyperpa- rameters of the relevance vector machine for regression.
9.23 ( ) www In Section 7.2.1 we used direct maximization of the marginal like- lihood to derive the re-estimation equations (7.87) and (7.88) for finding values of the hyperparameters α and β for the regression RVM. Similarly, in Section 9.3.4 we used the EM algorithm to maximize the same marginal likelihood, giving the re-estimation equations (9.67) and (9.68). Show that these two sets of re-estimation equations are formally equivalent.
9.24 ( ) Verify the relation (9.70) in which L(q, θ) and KL(q p) are defined by (9.71) and (9.72), respectively.
9.25 ( ) www Show that the lower bound L(q,θ) given by (9.71), with q(Z) = p(Z|X, θ(old)), has the same gradient with respect to θ as the log likelihood function ln p(X|θ) at the point θ = θ(old).
9.26 ( ) www Consider the incremental form of the EM algorithm for a mixture of Gaussians, in which the responsibilities are recomputed only for a specific data point xm. Starting from the M-step formulae (9.17) and (9.18), derive the results (9.78) and (9.79) for updating the component means.
9.27 ( ) Derive M-step formulae for updating the covariance matrices and mixing coefficients in a Gaussian mixture model when the responsibilities are updated in- crementally, analogous to the result (9.78) for updating the means.
Exercises 459