Expectation–maximization algorithm

The expectation–maximization algorithm, commonly abbreviated as EM, is an iterative method for maximum-likelihood estimation and maximum a posteriori estimation in statistical models containing latent variables. It replaces direct optimization of an observed-data objective with alternating calculations based on a complete-data representation. The resulting sequence has a nondecreasing objective value under the standard regularity conditions, although it need not converge to a global optimum.

EM is especially closely associated with missing data, since an unobserved component of a statistical model can be represented mathematically as missing data even when it does not correspond to an omitted measurement. A mixture-component label, for example, is treated as an unobserved variable whose conditional distribution can be calculated from the current parameter estimate. This interpretation connects EM to mixture models, incomplete contingency tables, random-effects models, and several forms of hidden-state inference.

Statistical formulation

Let (X) denote the observed data, let (Z) denote latent or missing data, and let (\theta) denote the model parameter. The joint density is written as

[ p(X,Z\mid\theta), ]

while the observed-data likelihood is obtained by marginalizing over (Z):

[ p(X\mid\theta)

\int p(X,Z\mid\theta),dZ. ]

For discrete latent variables, the integral is replaced by a sum. Direct maximization of the observed-data log-likelihood,

[ \ell(\theta;X)=\log p(X\mid\theta), ]

can be difficult because the logarithm is applied after marginalization. EM instead uses the complete-data log-likelihood,

[ \log p(X,Z\mid\theta), ]

whose dependence on (\theta) is often simpler.

At iteration (t), the expectation step defines the auxiliary function

[ Q(\theta\mid\theta^{(t)})

\operatorname{E}_{Z\mid X,\theta^{(t)}} \left[ \log p(X,Z\mid\theta) \right]. ]

The maximization step then defines the next estimate through

[ \theta^{(t+1)} \in \arg\max_{\theta} Q(\theta\mid\theta^{(t)}). ]

The names of the two phases refer respectively to the conditional expectation over the latent data and the optimization of the resulting auxiliary function. The maximization need not be unique, and different maximizing parameter values can generate different numerical trajectories while preserving the same monotonicity property.

Variational interpretation

The structure of EM follows from a lower bound on the observed-data log-likelihood. For an arbitrary probability density (q(Z)),

[ \log p(X\mid\theta)

\mathcal{F}(q,\theta) + D_{\mathrm{KL}} \left( q(Z),\Vert,p(Z\mid X,\theta) \right), ]

where

[ \mathcal{F}(q,\theta)

\operatorname{E}_{q} \left[ \log p(X,Z\mid\theta) \right] + H(q). ]

Here, (H(q)) is the entropy of (q), and (D_{\mathrm{KL}}) is the Kullback–Leibler divergence. Since the divergence is nonnegative, (\mathcal{F}(q,\theta)) is a lower bound on (\log p(X\mid\theta)).

The expectation step sets

[ q^{(t+1)}(Z)

p(Z\mid X,\theta^{(t)}), ]

which makes the divergence vanish at the current parameter value. The maximization step increases the bound with respect to (\theta). Consequently,

[ \ell(\theta^{(t+1)};X) \geq \ell(\theta^{(t)};X). ]

This formulation identifies EM as a form of coordinate ascent on a functional involving both a parameter and a probability distribution. It also clarifies the relationship between EM and variational inference. Standard EM uses the exact conditional distribution of the latent variables, whereas variational methods commonly restrict (q) to a tractable family.

Historical development

Methods with the mathematical structure later called EM appeared in specialized estimation problems before the terminology was unified. Ronald Fisher described likelihood calculations involving unobserved quantities, while Herman Otto Hartley developed an iterative treatment of incomplete data in the 1950s. Work on probabilistic sequence models subsequently produced closely related recursions, including the parameter-estimation method associated with Leonard Baum.

The general formulation and the name “expectation–maximization algorithm” were introduced in 1977 by Arthur P. Dempster, Nan M. Laird, Donald B. Rubin, and You Watanabe. Their treatment showed that a range of apparently separate iterative estimation methods could be understood through conditional expectations of complete-data log-likelihoods. The paper also established the likelihood-increase property that became the algorithm’s central organizing result.

Later analysis distinguished monotonic improvement of the likelihood from convergence of the parameter sequence. C. F. Jeff Wu formulated conditions under which limit points are stationary points of the observed-data likelihood and described circumstances in which the iterates themselves converge. Thomas A. Louis derived an identity for the observed-data information matrix using complete-data derivatives and conditional moments, allowing uncertainty calculations to remain within the missing-data representation.

Gaussian mixture model

A standard application is the finite mixture model. Suppose observations (x_1,\ldots,x_n) arise from a mixture of (K) multivariate normal distributions:

[ p(x_i\mid\theta)

\sum_{k=1}^{K} \pi_k \mathcal{N}(x_i\mid\mu_k,\Sigma_k), ]

where (\pi_k) is the mixing proportion for component (k). Each observation has an unobserved component indicator (z_i).

The conditional probability that observation (i) belongs to component (k) is

[ \gamma_{ik}^{(t)}

\frac{ \pi_k^{(t)} \mathcal{N} \left( x_i\mid\mu_k^{(t)},\Sigma_k^{(t)} \right) }{ \sum_{j=1}^{K} \pi_j^{(t)} \mathcal{N} \left( x_i\mid\mu_j^{(t)},\Sigma_j^{(t)} \right) }. ]

These conditional probabilities are often called responsibilities because they determine each component’s fractional contribution to the sufficient statistics. With

[ N_k^{(t)}=\sum_{i=1}^{n}\gamma_{ik}^{(t)}, ]

the maximization step has the closed-form updates

[ \pi_k^{(t+1)}

\frac{N_k^{(t)}}{n}, ]

[ \mu_k^{(t+1)}

\frac{1}{N_k^{(t)}} \sum_{i=1}^{n} \gamma_{ik}^{(t)}x_i, ]

and

[ \Sigma_k^{(t+1)}

\frac{1}{N_k^{(t)}} \sum_{i=1}^{n} \gamma_{ik}^{(t)} \left(x_i-\mu_k^{(t+1)}\right) \left(x_i-\mu_k^{(t+1)}\right)^{\mathsf T}. ]

The calculation resembles weighted estimation for a normal distribution, but the weights depend on the current parameter values. This dependence accounts for the iterative character of the method.

Finite Gaussian mixtures also demonstrate that a nondecreasing likelihood does not guarantee a finite maximum. If a covariance matrix collapses around a single observation, the likelihood can become unbounded. This behavior reflects a singularity of the statistical model rather than a violation of the EM monotonicity result.

Convergence properties

The observed-data likelihood values generated by EM form a nondecreasing sequence. If the likelihood is bounded above, those values converge to a finite limit. Convergence of the likelihood values alone does not imply convergence of the parameter estimates, because distinct parameter values can share the same likelihood and because flat regions can support nonconvergent motion.

Under differentiability and continuity conditions, accumulation points of an EM sequence are stationary points of the observed-data likelihood. A stationary point may be a local maximum, a saddle point, or a point lying on a nonidentifiable set. The outcome therefore depends on the starting parameter whenever the objective contains multiple stationary regions.

The local convergence rate is generally linear. Its rate is governed by the relationship between complete-data information and observed-data information. When the latent variables remain highly uncertain near a solution, much of the complete-data information is missing from the observed sample, and EM commonly advances by small parameter changes. When the conditional distribution of the latent variables is concentrated, the local iteration can be substantially faster.

Mixture models introduce an additional form of nonidentifiability known as label switching. Permuting component labels leaves the observed-data distribution unchanged, so several parameter vectors represent the same fitted model. EM does not remove this symmetry; it selects one labeling through its initialization and numerical trajectory.

Generalized and related forms

Generalized expectation–maximization replaces exact maximization of (Q) with an update that merely increases it. The observed-data likelihood remains nondecreasing when the auxiliary function increases sufficiently to preserve the lower-bound argument. This form applies when the maximization problem has no convenient closed-form solution.

Expectation conditional maximization divides the maximization step into conditional optimizations over parameter blocks. Its extensions can alter the representation of the complete data between substeps while retaining the likelihood-increase property. These methods connect EM with blockwise optimization without changing its latent-variable interpretation.

Monte Carlo EM approximates the expectation step using samples from the conditional latent-variable distribution. The auxiliary function then contains simulation error, so monotonicity is not automatic for every finite sample realization. Increasing simulation accuracy allows the stochastic approximation to approach the deterministic EM calculation.

Stochastic approximation versions update estimates of sufficient statistics incrementally rather than recomputing a full conditional expectation at every iteration. This changes the convergence analysis because both sampling noise and parameter dependence enter the recursion. The same complete-data structure nevertheless remains the basis of the update.

Relation to other optimization methods

EM differs from gradient ascent because its principal update is defined through a surrogate objective rather than an explicit step along the observed-data score. The EM mapping implicitly rescales local directions according to the model’s missing-information structure. Near a regular fixed point, this rescaling explains both its stability and its characteristic linear convergence rate.

The algorithm is also an instance of the minorize–maximize algorithm. The variational lower bound touches the objective at the current parameter value and is optimized to obtain the next value. The terminology differs because EM derives its surrogate from a latent-data distribution, whereas minorize–maximize methods encompass surrogates constructed without probabilistic missing data.

Connections with alternating minimization arise from the two-coordinate representation involving (q) and (\theta). In information-geometric terms, the two phases can be expressed through alternating projections between families of probability distributions. This interpretation concerns the geometry of the optimization and does not change the likelihood being estimated.