IEOR E4570: Machine Learning for OR&FE Spring 2015 …

1 ieor e4570 : Machine Learning for OR&FESpring 2015c 2015 by Martin HaughThe EM AlgorithmThe EM algorithm is used for obtaining maximum likelihood estimates of parameters when some of the data ismissing. More generally, however, the EM algorithm can also be applied when there islatent, unobserved,data which was never intended to be observed in the first place. In that case, we simply assume that the latentdata is missing and proceed to apply the EM algorithm. The EM algorithm has many applications throughoutstatistics. It is often used for example, in Machine Learning and data mining applications, and in Bayesianstatistics where it is often used to obtain the mode of the posterior marginal distributions of The Classical EM AlgorithmWe begin by assuming that the complete data-set consists ofZ= (X,Y)but that onlyXis observed.

2 Thecomplete-data log likelihood is then denoted byl( ;X,Y)where is the unknown parameter vector for whichwe wish to find the :The E-step of the EM algorithm computes the expected value ofl( ;X,Y)given the observed data,X, and the current parameter estimate, oldsay. In particular, we defineQ( ; old) :=E[l( ;X,Y)|X, old]= l( ;X,y)p(y|X, old)dy(1)wherep( |X, old)is the conditional density ofYgiven the observed data,X, and assuming = :The M-step consists of maximizing over the expectation computed in (1). That is, we set new:= max Q( ; old).We then set old= two steps are repeated as necessary until the sequence of new s converges. Indeed under very generalcircumstances convergence to a local maximum can be guaranteed and we explain why this is the case below. Ifit is suspected that the log-likelihood function has multiple local maximums then the EM algorithm should berun many times, using a different starting value of oldon each occasion.

3 The ML estimate of is then taken tobe the best of the set of local maximums obtained from the various runs of the EM Does the EM Algorithm Work?We usep( | )to denote a generic conditional PDF. Now observe thatl( ;X) = lnp(X | ) = ln p(X,y| )dy= ln p(X,y| )p(y|X, old)p(y|X, old)dyThe EM Algorithm2= lnE[p(X,Y | )p(Y |X, old)|X, old] E[ln(p(X,Y | )p(Y |X, old))|X, old](2)=E[lnp(X,Y | )|X, old] E[lnp(Y |X, old)|X, old]=Q( ; old) E[lnp(Y |X, old)|X, old](3)where (2) follows from Jensen s Inequality and since thelnfunction is concave. It is also clear (because theterm inside the expectation becomes a constant) that the inequality in (2) becomes an equality if we take = old. Lettingg( | old)denote the right-hand-side of (3), we therefore havel( ;X) g( | old)for all with equality when = old. Therefore any value of that increasesg( | old)beyondg( old| old)must also increasel( ;X)beyondl( old;X).

4 The M-step finds such a by maximizingQ( ; old)over whichis equivalent (why?) to maximizingg( | old)over . It is also worth mentioning that in many applications thefunctionQ( ; old)will be a convex function of and therefore easy to ExamplesExample 1(Missing Data in a Multinomial Model)Supposex:= (x1,x2,x3,x4)is a sample from a Mult(n, )distribution where =(12+14 ,14(1 ),14(1 ),14 ).The likelihood,L( ;x), is then given byL( ;x) =n!x1!x2!x3!x4!(12+14 )x1(14(1 ))x2(14(1 ))x3(14 )x4so that the log-likelihoodl( ;x)isl( ;x) =C+x1ln(12+14 )+ (x2+x3) ln (1 ) +x4ln ( )whereCis a constant that does not depend on . We could try to maximizel( ;x)over directly usingstandard non-linear optimization algorithms. However, in this example we will perform the optimization insteadusing the EM algorithm. To do this we assume that the complete data is given byy:= (y1,y2,y3,y4,y5)andthatyhas a Mult(n, )distribution where =(12,14 ,14(1 ),14(1 ),14 ).

5 However, instead of observingywe only observe(y1+y2,y3,y4,y5), , we only observex. We therefore takeX= (y1+y2,y3,y4,y5)and takeY=y2. The log-likelihood of the complete data is then given byl( ;X,Y) =C+y2ln ( ) + (y3+y4) ln (1 ) +y5ln ( )where againCis a constant containing all terms that do no depend on . It is also clear that the conditionaldensity ofYsatisfiesf(Y |X, ) =Bin(y1+y2, /41/2 + /4).The EM Algorithm3We can now implement the E-step and :Recalling thatQ( ; old) :=E[l( ;X,Y)|X, old], we haveQ( ; old) :=C+E[y2ln ( )|X, old] + (y3+y4) ln (1 ) +y5ln ( )=C+ (y1+y2)poldln ( ) + (y3+y4) ln (1 ) +y5ln ( )wherepold:= old/41/2 + old/4.(4)M-Step:We now maximizeQ( ; old)to find new. Taking the derivative we obtaindQd =(y1+y2) pold (y3+y4)1 +y5 which is zero when we take = newwhere new:=y5+pold(y1+y2)y3+y4+y5+pold(y1+y2)( 5)Equations (4) and (5) now define the EM iteration which begins with some (judiciously) chosen value of 2(A Simple Normal-Mixture Model)An extremely common application of the EM algorithm is to estimate the MLE of normal mixture models.

6 Thisis often used, for example, in clustering algorithms. Suppose for example thatX= (X1,..,Xn)are IIDrandom variables each with PDFfx(x) =m j=1pje (x j)2/2 2j 2 2jwherepj 0for alljand where pj= 1. The parameters in this model are thepj s, the j s and the j of trying to finding the maximum likelihood estimates of these parameters directly via numericaloptimization, we can use the EM algorithm. We do this by assuming the presence of an additional randomvariable,Ysay, whereP(Y=j) =pjforj= 1,..,m. The realized value ofYthen determines which of themnormal distributions generates the corresponding value ofX. There arensuch random variables,(Y1,..,Yn) :=Y. Note thatfx|y(xi|yi=j, ) =1 2 2je (xi j)2/2 2j(6)where := (p1,..,pm, 1,.., m, 1,.., m)is the unknown parameter vector and that the likelihood isgiven byL( ;X,Y) =n i=1pyi1 2 2yie (xi yi)2/2 EM algorithm starts with an initial guess, old, and then iterates the E-step and M-step as described belowuntil :We need to computeQ( ; old) :=E[l( ;X,Y)|X, old].

7 Indeed it is straightforward to show thatQ( ; old) =n i=1m j=1P(Yi=j|xi, old) ln(fx|y(xi|yi=j, )P(Yi=j| )).(7)The EM Algorithm4 Note thatfx|y(xi|yi=j, )is given by (6) and thatP(Yi=j| old) =pj,old. Finally, sinceP(Yi=j|xi, old) =P(Yi=j, Xi=xi| old)P(Xi=xi| old)=fx|y(xi|yi=j, old)P(Yi=j| old) mk=1fx|y(xi|yi=k, old)P(Yi=k| old)it is clear that we can compute (7) :We can now maximizeQ( ; old)by setting the vector of partial derivatives, Q/ , equal to0andsolving for new. After some algebra, we obtain j,new= ni=1xiP(Yi=j|xi, old) ni=1P(Yi=j|xi, old)(8) 2j,new= ni=1(xi j)2P(Yi=j|xi, old) ni=1P(Yi=j|xi, old)(9)pj,new=1nn i=1P(Yi=j|xi, old).(10)Given an initial estimate, old, the EM algorithm cycles through (8) to (10) repeatedly, setting old= newaftereach cycle, until the estimates A More General Version of the EM AlgorithmThe EM algorithm is often stated more generally using the language of information theory.

8 In this section1wewill describe this more general formulation and relate it back to EM algorithm as described in Section 1. Asbefore the goal is to maximize the likelihood function,L( ;X), which is given by2L( ;X) =p(X | ) = yp(X,y| )dy.(11)The implicit assumption underlying the EM algorithm is that it is difficult to optimizep(X | )with respect to but that it is much easier to optimizep(X,Y| ). We first introduce an arbitrary distribution,q(Y), over thelatent variables,Y, and note that we can decompose the log-likelihood,l( ;X), into two terms (the first ofwhich is sometimes called the energy term) according tol( ;X) := lnp(X | ) =L(q, ) energy +KL(q||pY|X)(12)whereL(q, )and KL(q||pY|X)are the likelihood and Kullback-Leibler (KL) divergence3which are given byL(q, ) = Yq(Y) ln(p(X,Y| )q(Y))(13)KL(q||pY|X) = Yq(Y) ln(p(Y|X, )q(Y)).

9 It is well-known (see the Appendix) that the KL divergence satisfies KL(q||pY|X) 0and equals0if and only ifq(Y) =pY|X. It therefore follows thatL(q, ) l( ;X)for all distributions,q( ). We can now use thedecomposition of (12) to define the EM algorithm. We begin with an initial parameter estimate, :The E-step maximizes the lower bound,L(q, old), with respect toq( )while keeping oldfixed. Inprinciple this is a variational problem since we are optimizing a functional, but the solution is easily note thatl( old;X)does not depend onq( ). It then follows from (12) (with = old) that maximizingL(q, old)is equivalent to minimizing KL(q||pY|X). Since this latter term is always non-negative we see that1 The material in this section is drawn fromPattern Recognition and Machine Learning (2006) by Chris discrete then we replace the integral in (11) with a KL divergence is also often called therelative EM Algorithm5L(q, old)is optimized when KL(q||pY|X) = 0which, by our earlier observation, is the case when we takeq(Y) =p(Y|X, old).

10 At this point we see that the lower bound,L(q, old), will now equal the current value ofthe log-likelihood,l( old;X).M-Step:In the M-step we keepq(Y)fixed and maximizeL(q, )over to obtain new. This will thereforecause the lower bound to increase (if it is not already at a maximum) which in turn means that the log-likelihoodmust also increase. Moreover, at this new value newit will no longer be the case that KL(q||pY|X) = 0and soby (12) the increase in the log-likelihood will be greater than the increase in the lower the General EM Algorithm with the Classical EM AlgorithmIt is instructive to compare the E-step and M-step of the general EM algorithm with the corresponding steps ofSection 1. To do this, first substituteq(Y) =p(Y|X, old)into (13) to obtainL(q, ) =Q( ; old) +constant(14)whereQ( ; old)is the expected complete-date log-likelihood as defined in (1) where the expectation is takenassuming = old.