View source for EM algorithm

''expectation-maximization algorithm''

An iterative optimization procedure [[#References|[a4]]] for computing

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201201.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a1)</td></tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201202.png" /> is a non-negative real valued function which is integrable as a function of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201203.png" /> for each <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201204.png" />. The EM algorithm was developed for statistical inference in problems with incomplete data or problems that can be formulated as such (e.g., with latent-variable modeling) and is a very popular method of computing maximum likelihood, restricted maximum likelihood, penalized maximum likelihood, and maximum posterior estimates. The name, EM algorithm, stems from this context, where  "E"  stands for the expectation (i.e., integration) step and  "M"  stands for the maximization step. In particular, the algorithm starts with an initial value <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201205.png" /> and iterates until convergence the two following steps for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201206.png" />:''''''<table border="0" cellpadding="0" cellspacing="0" style="background-color:black;"> <tr><td> <table border="0" cellspacing="1" cellpadding="4" style="background-color:black;"> <tbody> <tr> <td colname="1" style="background-color:white;" colspan="1">E-step:</td> </tr> <tr> <td colname="1" style="background-color:white;" colspan="1">Compute <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201207.png" />,</td> </tr> <tr> <td colname="1" style="background-color:white;" colspan="1">where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201208.png" />;</td> </tr> <tr> <td colname="1" style="background-color:white;" colspan="1">M-step:</td> </tr> <tr> <td colname="1" style="background-color:white;" colspan="1">Determine <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e1201209.png" /> by maximizing <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012010.png" />,</td> </tr> <tr> <td colname="1" style="background-color:white;" colspan="1">that is, find <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012011.png" /> so that</td> </tr> <tr> <td colname="1" style="background-color:white;" colspan="1"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012012.png" /> for all <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012013.png" />.</td> </tr> </tbody> </table>

</td></tr> </table>

The usefulness of the EM algorithm is apparent when both of these steps can be accomplished with minimal analytic and computation effort but either the direct maximization or the integration in (a1) is difficult. The attractive convergence properties of the algorithm along with its many generalizations and extensions will be discussed below. First, however, an illustration of the algorithm via several applications is given.

==Applications.==


===Maximum-likelihood estimation with missing observations.===
The standard notation and terminology for the EM algorithm stem from this important class of applications. In this case the function <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012014.png" /> is the complete-data likelihood function, which is written <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012015.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012016.png" /> is the complete data which consists of the observed data, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012017.png" />, and the missing data, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012018.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012019.png" /> is a model parameter, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012020.png" />. One is interested in computing the value of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012021.png" /> that maximizes <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012022.png" />, the maximum-likelihood estimate. Here it is assumed that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012023.png" /> is missing at random [[#References|[a26]]]; more complicated missing data mechanisms can also be accounted for in the formulation of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012024.png" />. More generally, one can augment <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012025.png" /> to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012026.png" />, the augmented data, via a many-to-one mapping, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012027.png" />, with the understanding that the augmented-data likelihood, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012028.png" />, is linked to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012029.png" /> via <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012030.png" />.

The EM algorithm builds on the intuitive notion that:

i) if there were no missing data, maximum-likelihood estimation would be easy; and

ii) if the model parameters were known, the missing data could easily be imputed (i.e., predicted) by its conditional expectation. These two observations correspond to the M-step and the E-step, respectively, with the proviso that not the missing data, but rather the complete-data log-likelihood should be imputed by its conditional expectation. In particular, the E-step reduces to computing

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012031.png" /></td> </tr></table>

the conditional expectation of the complete-data log-likelihood. If the complete-data model is from an exponential family, then <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012032.png" /> is linear in a set of complete-data sufficient statistics. In this common case which includes multivariate normal, Poisson, multinomial, and exponential models (among others), computing <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012033.png" /> involves routine calculations. The M-step then involves computing the maximum-likelihood estimates as if there were no missing data, by using the imputed complete-data sufficient statistics from the E-step as inputs.

===Linear inverse problems with positivity restrictions.===
Consider the system of equations

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012034.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012035.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012036.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012037.png" /> are known, and one wishes to solve for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012038.png" />. Such problems are omnipresent in scientific and engineering applications and can be solved with the EM algorithm [[#References|[a30]]]. First one notes that, without loss of generality, one may assume that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012039.png" />, that is, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012040.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012041.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012042.png" /> for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012043.png" /> are discrete probability measures (see [[#References|[a30]]] for details). Thus, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012044.png" /> can be viewed as the mixing weights in a finite mixture model, i.e., <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012045.png" />. A finite mixture model is a hierarchical model in that one can suppose that an integer <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012046.png" /> is chosen according to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012047.png" />, and then data are generated conditional on <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012048.png" /> according to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012049.png" />. The marginal distribution of the data generated is <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012050.png" />. If one considers <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012051.png" /> to be missing data, one can derive an EM algorithm of the form

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012052.png" /></td> </tr></table>

Y. Vardi and D. Lee [[#References|[a30]]] demonstrate that this iteration converges to the value that minimizes the [[Kullback–Leibler information|Kullback–Leibler information]] divergence between <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012053.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012054.png" /> over all non-negative <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012055.png" />, which is the desired non-negative solution if it exists.

===Multivariate <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012056.png" />-model, a latent variable model.===
The previous example illustrates an important principle in the application of the EM algorithm. Namely, one is often initially interested in optimizing a function <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012057.png" /> and purposely embeds <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012058.png" /> in a function, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012059.png" />, with a larger domain, such that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012060.png" />, in order to use the EM algorithm. In the previous example, this was accomplished by identifying <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012061.png" /> with <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012062.png" />. Similar strategies are often fruitful in various latent variable models. Consider, for example, the multivariate <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012063.png" />-model, which is useful for robust estimation (cf. [[Robust statistics|Robust statistics]]; see also, e.g., [[#References|[a9]]], [[#References|[a11]]]),

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012064.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a2)</td></tr></table>

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012065.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012066.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012067.png" /> is a positive-definite <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012068.png" />-matrix, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012069.png" />, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012070.png" /> is the known degree of freedom. Given a random sample <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012071.png" />, assumed to be from (a2), one wishes to find the maximizer of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012072.png" />, which is known to have no general closed-form solution. In order to apply the EM algorithm, one embeds <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012073.png" /> into a larger model <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012074.png" />. This is accomplished via the well-known representation

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012075.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a3)</td></tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012076.png" /> follows the [[T-distribution|<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012077.png" />-distribution]] in (a2), <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012078.png" />, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012079.png" /> is a mean chi-square variable with <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012080.png" /> degrees of freedom independent of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012081.png" />. In (a3), one sees that the distribution of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012082.png" /> conditional on <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012083.png" /> is multivariate normal. Thus, if one had observed <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012084.png" />, maximum-likelihood estimation would be easy. One therefore defines <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012085.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012086.png" />, which is easy to maximize as a function of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012087.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012088.png" /> due to the conditional normality of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012089.png" />. Thus, since <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012090.png" /> is linear in <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012091.png" />, the <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012092.png" />st iteration of the EM algorithm has a simple E-step that computes

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012093.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012094.png" />. The M-step then maximizes <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012095.png" /> to obtain

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012096.png" /></td> </tr></table>

and

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012097.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a4)</td></tr></table>

Given the weights <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012098.png" />, the M-step corresponds to weighted least squares, thus this EM algorithm, which often goes by the name iteratively reweighted least squares (IRLS), is trivial to program and use, and, as will be seen in the next section, exhibits unusually stable convergence.

==Properties of the EM algorithm.==
The EM algorithm has enjoyed wide popularity in many scientific fields from the 1970s onwards (e.g., [[#References|[a18]]], [[#References|[a23]]]). This is primarily due to easy implementation and stable convergence. As the previous paragraph illustrates, the EM algorithm is often easy to program and use. Although the algorithm may take many iterations to converge relative to other optimization routines (e.g., Newton–Raphson), each iteration is often easy to program and quick to compute. Moreover, the EM algorithm is less sensitive to poor starting values, and can be easier to use with many parameters since the iterations necessarily remain in the parameter space and no second derivatives are required. Finally, the EM algorithm has the very important property that the objective function is increased at each iteration. That is, by the definition of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e12012099.png" />,

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120100.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a5)</td></tr></table>

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120101.png" /></td> </tr></table>

where the second equality follows by averaging over <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120102.png" /> according to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120103.png" />. Since the first term of (a5) is maximized by <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120104.png" />, and the second is minimized by <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120105.png" /> (under the assumption that the support of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120106.png" /> does not depend on <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120107.png" />), one obtains

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120108.png" /></td> </tr></table>

for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120109.png" />. This property not only contributes to the stability of the algorithm, but also is very valuable for diagnosing implementation errors.

Although the EM algorithm is not guaranteed to converge to even a local mode (it can converge to a saddle point [[#References|[a25]]] or even a local minimum [[#References|[a1]]]), this can easily be avoided in practice by using several  "overdispersed"  starting values. (Details of convergence properties are developed in, a.o., [[#References|[a4]]], [[#References|[a32]]], [[#References|[a2]]], [[#References|[a21]]].) Running the EM algorithm with several starting values is also recommended because it can help one to find multiple local modes of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120110.png" />, an important advantage for statistical analysis.

==Extensions and enhancements.==
The advantages of the EM algorithm are diminished when either the E-step or the M-step are difficult to compute or the algorithm is very slow to converge. There are a number of extensions to the EM algorithm that can be successful in dealing with these difficulties.

For example, the ECM algorithm ([[#References|[a20]]], [[#References|[a17]]]) is useful when <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120111.png" /> is difficult to maximize as a function of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120112.png" />, but can be easily maximized when <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120113.png" /> is constrained to one of several subspaces of the original parameter space. One assumes that the aggregation of these subspaces is space-filling in the sense that the sequence of conditional maximizations searches the entire parameter space. Thus, ECM replaces the M-step with a sequence of CM-steps (i.e., conditional maximizations) while maintaining the convergence properties of the EM algorithm, including monotone convergence.

The situation is somewhat more difficult when the E-step is difficult to compute, since numerical integration can be very expensive computationally. The Monte-Carlo EM algorithm [[#References|[a31]]] suggests using Monte-Carlo integration in the E-step. In some cases, Markov chain Monte-Carlo integration methods have been used successfully (e.g., [[#References|[a15]]], [[#References|[a16]]], [[#References|[a22]]], [[#References|[a3]]]). The nested EM algorithm [[#References|[a27]]] offers a general strategy for efficient implementation of Markov chain Monte Carlo within the E-step of the EM algorithm.

Much emphasis in recent research is on speeding up the EM algorithm without sacrificing its stability or simplicity. A very fruitful method is through efficient data augmentation, that is, simple and  "small"  augmentation schemes, where  "small"  is measured by the so-called Fisher information — see [[#References|[a23]]] and the accompanying discussion for details. Briefly, the smaller the augmentation, the better <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120114.png" /> approximates <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120115.png" /> and thus the faster the algorithm; here, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120116.png" /> is the limit of the EM sequence — note that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120117.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120118.png" /> share the same stationary point(s) [[#References|[a4]]]. Thus, the SAGE [[#References|[a5]]], the ECME [[#References|[a11]]], and the more general the AECM [[#References|[a23]]] algorithms build on the ECM algorithm by allowing this approximation to vary from CM-step to CM-step, in order to obtain a better approximation for some CM-steps as long as this does not complicate the resulting CM-step. X.L. Meng and D.A. van Dyk [[#References|[a23]]] also introduce a working parameter into <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120119.png" />, that is, they define <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120120.png" /> such that

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120121.png" /></td> </tr></table>

for all <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120122.png" /> in some class <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120123.png" />. They then choose <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120124.png" /> in order to optimize the approximation in terms of the rate of convergence of the resulting algorithm while maintaining the simplicity and stability of the EM algorithm. This results in simple, stable, and fast algorithms in a number of statistical applications. For example, in the <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120125.png" />-model this procedure results in replacing the update for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120126.png" /> given in (a4) with

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120127.png" /></td> </tr></table>

The rest of the algorithm remains the same and this simple change clearly maintains the simplicity of the algorithm but can result in a dramatic reduction in computational time compared to the standard IRLS algorithm.

The parameter-expanded EM algorithm or PXEM algorithm [[#References|[a12]]] also works with a working parameter, but rather than conditioning on the optimal <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120128.png" /> it marginalizes out <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120129.png" /> by fitting it in each iteration. Detailed discussion and comparison of such conditional augmentation and marginal augmentation in the more general context of stochastic simulations (e.g., the Gibbs sampler) can be found in [[#References|[a24]]], [[#References|[a13]]], and [[#References|[a28]]].

Other acceleration techniques have been developed by combining the EM algorithm with various numerical methods, such as Aitken acceleration (e.g., [[#References|[a14]]]), Newton–Raphson [[#References|[a7]]], quasi-Newton methods (e.g., [[#References|[a8]]]), and conjugate-gradient acceleration [[#References|[a6]]]. These methods, however, typically sacrifice monotone convergence and therefore require extra programming and special monitoring.

A final pair of extensions are designed to compute

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120130.png" /></td> </tr></table>

whose inverse gives a large-sample variance of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120131.png" /> in statistical analysis. The supplemented EM [[#References|[a19]]] and supplemented ECM [[#References|[a29]]] algorithms combine the analytically computable

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120132.png" /></td> </tr></table>

with numerical differentiation of the EM mapping in order to compute <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/e/e120/e120120/e120120133.png" /> via a simple matrix identity.

====References====
<table><TR><TD valign="top">[a1]</TD> <TD valign="top">  O. Arslan,   P.D.L. Constable,   J.T. Kent,   "Domains of convergence for the EM algorithm: A cautionary tale in a location estimation problem"  ''Statist. Comput.'' , '''3'''  (1993)  pp. 103–108</TD></TR><TR><TD valign="top">[a2]</TD> <TD valign="top">  R.A. Boyles,   "On the convergence of the EM algorithm"  ''J.R. Statist. Soc.'' , '''B45'''  (1983)  pp. 47–50</TD></TR><TR><TD valign="top">[a3]</TD> <TD valign="top">  J.S.K. Chan,   A.Y.C. Kuk,   "Maximum likelihood estimation for probit-linear mixed models with correlated random effects"  ''Biometrics'' , '''53'''  (1997)  pp. 86–97</TD></TR><TR><TD valign="top">[a4]</TD> <TD valign="top">  A.P. Dempster,   N.M. Laird,   D.B. Rubin,   "Maximum likelihood estimation from incomplete-data via the EM algorithm (with discussion)"  ''J. R. Statist. Soc.'' , '''B39'''  (1977)  pp. 1–38</TD></TR><TR><TD valign="top">[a5]</TD> <TD valign="top">  J.A. Fessler,   A.O. Hero,   "Space-alternating generalized expectation-maximization algorithm"  ''IEEE Trans. Signal Processing'' , '''42'''  (1994)  pp. 2664–77</TD></TR><TR><TD valign="top">[a6]</TD> <TD valign="top">  M. Jamshidian,   R.I. Jennrich,   "Conjugate gradient acceleration of the EM algorithm"  ''J. Amer. Statist. Assoc.'' , '''88'''  (1993)  pp. 221–228</TD></TR><TR><TD valign="top">[a7]</TD> <TD valign="top">  K. Lange,   "A gradient algorithm locally equivalent to the EM algorithm"  ''J. R. Statist. Soc.'' , '''B57'''  (1995)  pp. 425–438</TD></TR><TR><TD valign="top">[a8]</TD> <TD valign="top">  K. Lange,   "A quasi-Newtonian acceleration of the EM algorithm"  ''Statistica Sinica'' , '''5'''  (1995)  pp. 1–18</TD></TR><TR><TD valign="top">[a9]</TD> <TD valign="top">  K. Lange,   R.J.A. Little,   J.M.G. Taylor,   "Robust statistical modeling using the t-distribution"  ''J. Amer. Statist. Assoc.'' , '''84'''  (1989)  pp. 881–896</TD></TR><TR><TD valign="top">[a10]</TD> <TD valign="top">  C. Liu,   D.B. Rubin,   "The ECME algorithm: a simple extension of EM and ECM with fast monotone convergence"  ''Biometrika'' , '''81'''  (1994)  pp. 633–48</TD></TR><TR><TD valign="top">[a11]</TD> <TD valign="top">  C. Liu,   D.B. Rubin,   "ML estimation of the t-distribution using EM and its extensions, ECM and ECME"  ''Statistica Sinica'' , '''5'''  (1995)  pp. 19–40</TD></TR><TR><TD valign="top">[a12]</TD> <TD valign="top">  C. Liu,   D.B. Rubin,   Y.N. Wu,   "Parameter expansion to accelerate EM: the PX-EM algorithm"  ''Biometrika''  (1998)  pp. 755–770</TD></TR><TR><TD valign="top">[a13]</TD> <TD valign="top">  J.S. Liu,   Y. Wu,   "Parameter expansion for data augmentation"  ''J. Amer. Statist. Assoc.'' , '''to appear'''  (1999)</TD></TR><TR><TD valign="top">[a14]</TD> <TD valign="top">  T.A. Louis,   "Finding the observed information matrix when using the EM algorithm"  ''J. R. Statist. Soc.'' , '''B44'''  (1982)  pp. 226–233</TD></TR><TR><TD valign="top">[a15]</TD> <TD valign="top">  C.E. McCulloch,   "Maximum likelihood variance components estimation for binary data"  ''J. Amer. Statist. Assoc.'' , '''89'''  (1994)  pp. 330–335</TD></TR><TR><TD valign="top">[a16]</TD> <TD valign="top">  C.E. McCulloch,   "Maximum likelihood algorithms for generalized linear mixed models"  ''J. Amer. Statist. Assoc.'' , '''92'''  (1997)  pp. 162–170</TD></TR><TR><TD valign="top">[a17]</TD> <TD valign="top">  X.L. Meng,   "On the rate of convergence of the ECM algorithm"  ''Ann. Math. Stat.'' , '''22'''  (1994)  pp. 326–339</TD></TR><TR><TD valign="top">[a18]</TD> <TD valign="top">  X.L. Meng,   S. Pedlow,   "EM: A bibliographic review with missing articles" , ''Proc. Statist. Comput. Sect. 24-27 Washington, D.C.'' , Amer. Statist. Assoc.  (1992)</TD></TR><TR><TD valign="top">[a19]</TD> <TD valign="top">  X.L. Meng,   D.B. Rubin,   "Using EM to obtain asymptotic variance-covariance matrices: the SEM algorithm"  ''J. Amer. Statist. Assoc.'' , '''86'''  (1991)  pp. 899–909</TD></TR><TR><TD valign="top">[a20]</TD> <TD valign="top">  X.L. Meng,   D.B. Rubin,   "Maximum likelihood estimation via the ECM algorithm: a general framework"  ''Biometrika'' , '''80'''  (1993)  pp. 267–78</TD></TR><TR><TD valign="top">[a21]</TD> <TD valign="top">  X.L. Meng,   D.B. Rubin,   "On the global and componentwise rates of convergence of the EM algorithm"  ''LALG'' , '''199'''  (1994)  pp. 413–425</TD></TR><TR><TD valign="top">[a22]</TD> <TD valign="top">  X.L. Meng,   S. Schilling,   "Fitting full-information item factor models and an empirical investigation of bridge sampling"  ''J. Amer. Statist. Assoc.'' , '''91'''  (1996)  pp. 1254–1267</TD></TR><TR><TD valign="top">[a23]</TD> <TD valign="top">  X.L. Meng,   D.A. van Dyk,   "The EM algorithm — an old folk song sung to a fast new tune (with discussion)"  ''J. R. Statist. Soc.'' , '''B59'''  (1997)  pp. 511–567</TD></TR><TR><TD valign="top">[a24]</TD> <TD valign="top">  X.L. Meng,   D.A. van Dyk,   "Seeking efficient data augmentation schemes via conditional and marginal augmentation"  ''Biometrika''  (1999)  pp. 301–320</TD></TR><TR><TD valign="top">[a25]</TD> <TD valign="top">  G.D. Murray,   "Comments on the paper by A.P. Dempster, N.M. Laird and D.B. Rubin"  ''J. R. Statist. Soc.'' , '''B39'''  (1977)  pp. 27–28</TD></TR><TR><TD valign="top">[a26]</TD> <TD valign="top">  D.B. Rubin,   "Inference and missing data"  ''Biometrika'' , '''63'''  (1976)  pp. 581–592</TD></TR><TR><TD valign="top">[a27]</TD> <TD valign="top">  D.A. van Dyk,   "Nesting EM Algorithms for computational efficiency"  ''Statistica Sinica'' , '''to appear'''  (2000)</TD></TR><TR><TD valign="top">[a28]</TD> <TD valign="top">  D.A. van Dyk,   X.L. Meng,   "The art of data augmentation"  ''Statist. Sci.'' , '''Submitted'''  (1999)</TD></TR><TR><TD valign="top">[a29]</TD> <TD valign="top">  D.A. van Dyk,   X.L. Meng,   D.B. Rubin,   "Maximum likelihood estimation via the ECM algorithm: computing the asymptotic variance"  ''Statistica Sinica'' , '''5'''  (1995)  pp. 55–75</TD></TR><TR><TD valign="top">[a30]</TD> <TD valign="top">  Y. Vardi,   D. Lee,   "From image deblurring to optimal investments: maximum likelihood solutions for positive linear inverse problems (with discussion)"  ''J. R. Statist. Soc.'' , '''B55'''  (1993)  pp. 569–598</TD></TR><TR><TD valign="top">[a31]</TD> <TD valign="top">  G.C.G. Wei,   M.A. Tanner,   "A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithm"  ''J. Amer. Statist. Assoc.'' , '''85'''  (1990)  pp. 699–704</TD></TR><TR><TD valign="top">[a32]</TD> <TD valign="top">  C.F.J. Wu,   "On the convergence properties of the EM algorithm"  ''Ann. Math. Stat.'' , '''11'''  (1983)  pp. 95–103</TD></TR></table>