AppendixA Expectation-Maximization (EM) Method

Appendix AExpectation-Maximization (EM) MethodExpectation-Maximization (EM) is a method for deriving algorithms to maximizelikelihood functions—using it with different likelihood functions results in differentML (or MAP, as the case may be) estimators. It is, on first encounter, a mysteriousmethod with strong appeal to those who love to see elaborate mathematical constructions telescope into simple forms.EM is reviewed in this appendix. The distinctions between the named randomvariables—incomplete, missing, and complete—are discussed, together with theauxiliary function called Q, and the E- and M-steps. It is shown that increasing theauxiliary function increases the likelihood function. This means that all EM basedalgorithms monotonically increase the data likelihood function p Z (z ; θ ). Iterativemaximization is a strictly numerical analysis method that interprets the auxiliaryfunction Q as a two parameter family of curves whose envelop is the data likelihoodfunction.EM works amazingly well in some applications and is seemingly useless inothers. It is often a hand-in-glove fit for likelihood functions that involve sums(integrals), for in these cases the natural missing data are the summation indices(variables of integration) and they simplify the likelihood function to a product ofthe summands (integrands). Such products are often much simpler to manipulate.However, other missing data are possible, and different choices of missing data canlead to different EM convergence rates.A.1 FormulationThe observed data z are called the incomplete data. As is seen shortly, the namemakes sense in the context of the EM method. The measurements z are a realizationof the random variable Z , and Z takes values in the space Z. The pdf of the dataz is specified parametrically by p Z (z ; θ ), where θ Θ is an unknown parametervector. Here, Θ is the set of all valid parameter values.The likelihood function of θ is the pdf p Z (z ; θ ) thought of as a function of θ fora given z. It is assumed that the likelihood of the data for any θ Θ is finite. It isalso assumed that the likelihood function of θ is uniformly bounded above, that is,p Z (z ; θ ) B for all θ Θ. The latter assumption is important.R.L. Streit, Poisson Point Processes, DOI 10.1007/978-1-4419-6923-1, C Springer Science Business Media, LLC 2010223

224Appendix A: Expectation-Maximization (EM) MethodThe maximum likelihood estimate of θ isθ̂ M L arg max p Z (z ; θ ) .θ ΘFor pdfs differentiable with respect to θ , the natural way to compute θ̂ M L is solvethe so-called necessary conditions θ p Z (z ; θ ) 0using, say, a safeguarded Newton-Raphson algorithm. This approach can encounterdifficulties in some nonlinear problems, especially when the dimension of θ is large.The EM method is an alternative way to compute θ̂ M L . Let K be a random variable whose realizations k occur in the set K, that is, k K. For ease of exposition,the discussion in this appendix assumes that K is a continuous variable. If K isdiscrete, the discussion below is modified merely by replacing integrals over K bysums over K.Now, K is called missing data in the sense of EM if p Z (z ; θ ) Kp Z K (z, k ; θ ) dk ,(A.1)where p Z K (z, k ; θ ) is the joint pdf of Z and K . The pair of variables (Z , K )are called the complete data. In words, (A.1) says that the pdf of the data Z is themarginal of the complete data pdf over the missing data K . The user of the EMmethod is at liberty to choose K and the joint pdf p Z K ( · ) as fancy or need dictates,as long as the (A.1) holds. This is crucial. In many applications, however, the choiceof missing data is very natural.The pdf of K conditioned on Z is, by the definition of conditioning,p K X (k x ; θ ) p Z K (z, k ; θ ) p Z (z ; θ )p Z K (z, k ; θ ).K p Z K (z, k ; θ ) dk(A.2)Let n 0 denote the EM iteration index. Let θ (0) Θ be an initial (valid) value ofthe parameter.A.1.1 E-stepFor n 1, the EM auxiliary function is defined by the conditional expectation: Q θ ; θ (n 1) E K X ;θ (n 1) log p Z K (z, k ; θ ) ) (log p Z K (z, k ; θ )) p K Z k ) z ; θ (n 1) dk .K(A.3)

A.1Formulation225 The sole purpose of the E-step is to evaluate Q θ ; θ (n 1) . In practice, this is veryoften a purely symbolic step, that is, the expectation (A.3) is manipulated so that ittakes a simple analytic form. When manipulation does not yield an analytic form, itis necessary to replace the expectation with a suitably chosen discrete sum (defined,e.g., via Monte Carlo integration).A.1.2 M-stepThe EM updated parameter is θ (n) arg max Q θ ; θ (n 1) .θ Θ(A.4)The maximum in (A.4) can be computed by any available method. It is sometimesthought that it is necessary to solve the maximization step explicitly, but this isemphatically not so—a strictly numerical method is quite acceptable.It is unnecessary to solve for the maximum in (A.4). Instead it is only necessaryto find a value θ (n) such that Q θ (n) ; θ (n 1) Q θ (n 1) ; θ (n 1) .(A.5)If the update θ (n) is chosen in this way, the method is called the generalized EM(GEM) method. The GEM method is very useful in many problems. It is the startingpoint of other closely related EM-based methods. A prominent example is the SAGE(Space Alternating Generalized EM) algorithm [35].A.1.3 ConvergenceUnder mild conditions, convergence is guaranteed to a critical point of the likelihood function, that is, θ (n) θ such that θ p Z (z θ ) θ θ 0 .Experience shows that in practice, θ is almost certainly a local maximum and nota saddle point, so that θ θ̂ M L .One of the mild conditions that cannot be overlooked is the assumption thatthe likelihood function of θ is uniformly bounded above. If it is not, then the EMiteration will converge only if the initial value of θ is by chance in the domain ofattraction of a point of local maximum likelihood. Otherwise, it will diverge, meaning that the likelihood function of the iterates will grow unbounded. Unboundednesswould not be an issue were it not for the fact that the only valid values of θ havefinite likelihoods, so the iterates—if they converge to a point—converge to an invalidparameter value, that is, to a point not in Θ. Estimation of heteroscedastic Gaussian

226Appendix A: Expectation-Maximization (EM) Methodsums is notorious for this behavior (due to covariance matrix collapse). See Section3.4 for further discussion.The details of the EM convergence proof are given in too many places to repeatall the details here. The book [80] gives a full and careful treatment of convergence,as well as examples where convergence fails in various ways.For present purposes, it is enough to establish two facts. One is that if the updateθ (n) satisfies (A.5), then it also increases the likelihood function. The other is thata critical point of the auxiliary function is also a critical point of the likelihoodfunction, and conversely.To see that each EM step monotonically increases the likelihood function, recallthat log x x 1 for all x 0 with equality if and only if x 1. Then, usingonly the above definitions, 0 Q θ ; θ (n 1) Q θ (n 1) ; θ (n 1) ) (log p Z K (z, k ; θ )) p K Z k ) z ; θ (n 1) dkK ) log p Z K z, k ; θ (n 1)p K Z k ) z ; θ (n 1) dkK,2' & )p Z K (z, k ; θ )) z ; θ (n 1) dk klogpK Zp Z K z, k ; θ (n 1)K,2 p Z K z, k ; θ (n 1)p Z K (z, k ; θ ) 1 dkp Z (z ; θ (n 1) )p Z K z, k ; θ (n 1)K * 1(n 1) pz,k;θdk k;θ p(z,)ZKZKp Z z ; θ (n 1) K* 1 p Z (z ; θ ) p Z z ; θ (n 1) . p Z z ; θ (n 1)Clearly, any increase in Q will increase the likelihood function.To see that critical points of the likelihood function are critical points of theauxiliary function, and conversely, simply note that0 θ0 p Z (z ; θ0 ) θ0 p Z K (z, k ; θ0 ) dk K p Z K (z, k ; θ0 ) θ0 log p Z K (z, k ; θ0 ) dkK p Z K (z, k ; θ0 ) θ0 log p Z K (z, k ; θ0 ) dk p Z (z ; θ0 )p Z (z ; θ0 )K

A.2Iterative Majorization227 p Z (z ; θ0 )K p K Z (k z ; θ0 ) θ log p Z K (z, k θ ) θ θ dk0 p Z (z ; θ0 ) [ θ Q (θ ; θ0 )]θ θ0 .For further discussion and references to the literature, see [80].A.2 Iterative MajorizationThe widely referenced paper [23] is the synthesis of several earlier discoveries ofEM in a statistical setting. However, EM is not essentially statistical in character, butis rather only one member of a larger class of strictly numerical methods called iterative majorization [19–21]. This connection is not often mentioned in the literature,but the insight was noticed almost immediately after the publication of [23].The observed data likelihood function p Z (z ; θ ) is the envelop of a two parameter family of functions. This family is defined in the EM method by the auxiliaryfunction, Q(θ ; φ). As seen from the results in [23], the data likelihood functionmajorizes every function in this family. For each specified parameter φ, the function Q(θ ; φ) is tangent to the p Z (z ; θ ) at the point φ. This situation is depictedin Fig. A.1 for the sequence φ θ0 , θ1 , θ2 , . . . . It is now intuitively clear thatEM based algorithms monotonically increase the data likelihood function, and thatthey converge with high probability to a local maximum of the likelihood functionp Z (z ; θ ) that depends on the starting point.It is very often observed in practice that EM based algorithms make large stridestoward the solution in the early iterations, but that progress toward the solutionFig. A.1 Iterative majorization interpretation of the observed (incomplete) data likelihood functionas the envelop of the EM auxiliary function, Q(θ, φ)

228Appendix A: Expectation-Maximization (EM) Methodslows significantly as the iteration progresses. The iterative majorization interpretation of EM provides an intuitive insight into both phenomena. See [80] for furtherdiscussion and, in particular, for a proof that the rate of convergence is ultimatelyonly linear. This convergence behavior explains the many efforts proposed in theliterature to speed up EM convergence.A.3 Observed InformationThe observed information matrix (OIM) corresponding to the ML estimate θ̂ M L isdefined by O I M(θ̂ M L ) θ ( θ log p Z (z ; θ ))Tθ θ̂ M L.(A.6)The matrix O I M(θ ) is an OIM only if θ is an ML estimate. Evaluating the OIMis often a relatively straightforward procedure. Moreover, EM based ML estimation algorithms are easily adapted to compute the OIM as a by-product of the EMiteration at very little additional computational expense. See [69] for details.The OIM and FIM differ in important ways: The FIM computes the expectation of the negative Hessian of the loglikelihoodfunction with respect to the data, while the OIM does not. The OIM is evaluated at the ML estimate θ̂ M L , while the FIM is evaluated at thetrue value of the parameter.The expected value of the OIM is sometimes said to be the FIM, but this is notprecisely correct.The OIM and FIM are alike in that they both evaluate the negative Hessian ofthe loglikelihood function. Because of this and other similarities between them, theOIM is often used as a surrogate for the FIM when the FIM is unknown or otherwiseunavailable. See [26] for further discussion of the OIM, its origins, and potentialutility.