Expectation-Maximization Algorithm

The Expectation-Maximization algorithm (EM algorithm) is an iterative numerical optimization method in computation and algorithms and mathematical statistics for finding maximum-likelihood (ML) or maximum a posteriori (MAP) estimates of parameters θ in probabilistic models that depend on unobserved latent variables Z or incomplete data. Formalized and named in a landmark 1977 paper in the Journal of the Royal Statistical Society, Series B by Arthur P. Dempster, Nan M. Laird, and Donald B. Rubin, the algorithm unified dozens of previously isolated domain methods—including Cedric A. B. Smith's 1957 gene-counting estimator for human ABO blood-group allele frequencies, H. O. Hartley's 1958 missing-data analysis of variance, and Leonard E. Baum and Ted Petrie's 1966–1970 Baum–Welch forward-backward algorithm for hidden Markov models.

When a dataset contains observed variables X alongside hidden cluster assignments, missing measurements, or latent state paths Z, directly maximizing the marginal log-likelihood ℓ(θ; X) = ln ∑_Z p(X, Z | θ) requires differentiating the logarithm of a sum, which couples all model parameters together and prevents closed-form analytic solutions. The EM algorithm circumvents this coupling through minorization-maximization: at each iteration t, the E-step (Expectation step) uses the current parameter vector θ^(t) to compute the posterior distribution p(Z | X, θ^(t)) of the hidden variables and constructs a tight lower-bound surrogate function Q(θ | θ^(t)), and the M-step (Maximization step) maximizes that surrogate to obtain θ^(t+1), provably increasing (or leaving unchanged) the observed-data likelihood at every step without tuning learning-rate hyperparameters.

Why Marginal Log-Likelihood Is Intractable While Complete Data Is Simple

Let X denote the observed dataset, let Z denote the unobserved latent variables (or missing values), and let (X, Z) denote the complete data governed by joint probability density or mass function p(X, Z | θ) parameterized by θ ∈ Θ. The goal of maximum-likelihood estimation is to maximize the incomplete-data (marginal) log-likelihood:

ℓ(θ; X) = ln p(X | θ) = ln ∑_Z p(X, Z | θ)

(with summation replaced by integration ∫ p(X, Z | θ) dZ for continuous latent variables). Even when the complete-data distribution p(X, Z | θ) belongs to an exponential-family probability distribution of the form exp(η(θ)^TT(X, Z) − A(θ))—where maximizing the complete-data log-likelihood ln p(X, Z | θ) would reduce to a single linear closed-form equation in the sufficient statistics T(X, Z)—the presence of the sum ∑_Z inside the logarithm in ℓ(θ; X) destroys that linear structure.

Jensen's Inequality and the Evidence Lower Bound (ELBO)

Following the variational free-energy formulation of Radford M. Neal and Geoffrey E. Hinton (1998), introduce any arbitrary auxiliary probability distribution q(Z) over the latent variables (q(Z) > 0, ∑_Z q(Z) = 1). Multiplying and dividing by q(Z) inside the marginal log-likelihood and applying Jensen's inequality to the concave logarithm function yields:

ln p(X | θ) = ln ∑_Z q(Z) · [ p(X, Z | θ) / q(Z) ] ≥ ∑_Z q(Z) ln [ p(X, Z | θ) / q(Z) ] ≡ ℒ(q, θ)

The functional ℒ(q, θ) is the Evidence Lower Bound (ELBO) (or negative variational free energy). What is the difference between the true marginal log-likelihood ln p(X | θ) and the lower bound ℒ(q, θ)? Factoring the joint distribution via conditional probability as p(X, Z | θ) = p(Z | X, θ) · p(X | θ) gives the fundamental identity:

ln p(X | θ) = ℒ(q, θ) + D_KL( q(Z) ‖ p(Z | X, θ) )

where D_KL(q ‖ p) = ∑_Z q(Z) ln[q(Z) / p(Z | X, θ)] ≥ 0 is the Kullback–Leibler divergence between the auxiliary distribution q(Z) and the true Bayesian inference posterior distribution p(Z | X, θ) of the latent variables.

Alternating Coordinate Ascent on ℒ(q, θ)

Because D_KL(q ‖ p) ≥ 0 for all distributions and equals 0 if and only if q(Z) = p(Z | X, θ), the EM algorithm maximizes ℓ(θ; X) by alternating coordinate ascent on the Evidence Lower Bound ℒ(q, θ):

1. E-Step (Expectation Step — Maximize ℒ with respect to q at fixed θ^(t)):

Set q^(t+1)(Z) = p(Z | X, θ^(t)), the conditional posterior distribution of the latent variables given the observed data X and current parameters θ^(t). At this choice of q, the KL divergence vanishes (D_KL = 0), making the lower bound tangent to the true log-likelihood at θ^(t): ℒ(q^(t+1), θ^(t)) = ln p(X | θ^(t)) Expanding ℒ(q^(t+1), θ) = ∑_Z p(Z | X, θ^(t)) ln p(X, Z | θ) + H(q^(t+1)) shows that the only term depending on θ is Dempster, Laird, and Rubin's auxiliary Q-function (the expected complete-data log-likelihood): Q(θ | θ^(t)) = E_Z|X,θ^(t)[ ln p(X, Z | θ) ] = ∑_Z p(Z | X, θ^(t)) ln p(X, Z | θ)

1. M-Step (Maximization Step — Maximize ℒ with respect to θ at fixed q^(t+1)):

Choose the updated parameter vector θ^(t+1) to maximize Q(θ | θ^(t)): θ^(t+1) = arg max_θ∈Θ Q(θ | θ^(t))

Proof of Monotonic Likelihood Improvement

Why can the EM algorithm never decrease the observed-data likelihood? Using the KL decomposition at θ^(t+1) and θ^(t) with q^(t+1)(Z) = p(Z | X, θ^(t)):

• At θ^(t): ℓ(θ^(t); X) = Q(θ^(t) | θ^(t)) + H(q^(t+1)) + 0
• At θ^(t+1): ℓ(θ^(t+1); X) = Q(θ^(t+1) | θ^(t)) + H(q^(t+1)) + D_KL(p(Z | X, θ^(t)) ‖ p(Z | X, θ^(t+1)))

Subtracting the first equation from the second yields:

ℓ(θ^(t+1); X) − ℓ(θ^(t); X) = [ Q(θ^(t+1) | θ^(t)) − Q(θ^(t) | θ^(t)) ] + D_KL( p(Z | X, θ^(t)) ‖ p(Z | X, θ^(t+1)) ) ≥ 0

Because the M-step ensures Q(θ^(t+1) | θ^(t)) ≥ Q(θ^(t) | θ^(t)) and Gibbs' inequality ensures D_KL ≥ 0, the marginal log-likelihood increases monotonically at every iteration unless θ^(t) is already a stationary point of ℓ(θ; X). In 1983, C. F. Jeff Wu proved that if Q(θ | θ′) is continuous in both arguments and the level sets of ℓ(θ; X) are compact, every limit point of an EM sequence {θ^(t)} is a stationary point of ℓ(θ; X). When closed-form maximization of Q is unavailable, Generalized EM (GEM) performs a single gradient descent or Newton step that merely increases Q(θ | θ^(t)), preserving the monotonic ascent guarantee.

In a K-component Gaussian Mixture Model over N observed d-dimensional data points X = {x₁, ..., x_N}, each observation x_i ∈ R^d is generated by selecting a hidden cluster label z_i ∈ {1, ..., K} with prior mixing weight π_k = Pr(z_i = k) (∑_k=1^K π_k = 1) and drawing x_i ~ 𝒩(x_i | μ_k, Σ_k).

In the limit where all covariance matrices are fixed to isotropic spheres Σ_k = εI with ε → 0⁺, the soft posterior responsibilities γ_ik collapse to hard binary indicator assignments ∈ {0, 1}, recovering Lloyd's K-means clustering algorithm as a zero-noise limit of EM.

In bioinformatics and computational biology, the EM algorithm powers three foundational pipelines:

1. Baum–Welch Training of Hidden Markov Models: For profile hidden Markov models (used in Pfam and HMMER protein domain databases) and speech recognition, the E-step runs the forward-backward dynamic programming algorithm to compute the posterior expected number of transitions ξ_t(i, j) and symbol emissions γ_t(i) across all hidden state paths, and the M-step normalizes those expected counts into updated transition and emission probabilities.
2. De Novo Transcription Factor Motif Discovery (MEME): Introduced by Timothy L. Bailey and Charles Elkan in 1994, MEME (Multiple EM for Motif Elicitation) treats the unknown starting coordinate of a regulatory binding motif of width W inside a set of unaligned promoter DNA sequences as a latent indicator variable Z_ij, alternating between computing the posterior probability that a motif starts at offset j (E-step) and updating the 4 × W position weight matrix (PWM) of nucleotide frequencies (M-step).
3. RNA-Seq Isoform Quantification and Haplotype Phasing: Because short cDNA sequencing reads frequently map ambiguously to multiple alternatively spliced transcript isoforms of the same gene, quantification engines (RSEM, Kallisto, and Salmon) model the true transcript of origin of each read as a latent categorical variable and run accelerated EM (SQUAREM) to estimate transcript-Per-Million (TPM) abundances. Analogously, the Excoffier–Slatkin (1995) EM algorithm resolves diploid genotype data into phased gametic haplotype frequencies under Hardy–Weinberg equilibrium.

Does the Expectation-Maximization algorithm always find the global maximum-likelihood estimate?

No. Because the marginal log-likelihood surface ℓ(θ; X) of latent-variable models is generally non-concave and invariant under permutations of cluster labels, the EM algorithm converges monotonically to a local maximum (or saddle point) determined by the starting initialization θ⁽⁰⁾. In practice, EM is initialized using K-means++ seeding or run from multiple random starting points, selecting the solution with the highest converged log-likelihood or best Akaike information criterion (AIC) / Bayesian Information Criterion (BIC).

What happens when the E-step posterior p(Z | X, θ) cannot be computed in closed form?

When the latent variable space is too complex for exact marginalization (for example, in deep hierarchical graphical models or continuous phylogenetic trees), two major extensions replace the exact E-step: Monte Carlo EM (MCEM) approximates Q(θ | θ^(t)) by drawing Markov Chain Monte Carlo samples from p(Z | X, θ^(t)) via the Metropolis–Hastings algorithm, while Variational EM restricts q(Z) to a tractable factorized family (such as mean-field distributions q(Z) = ∏ q_i(Z_i)) that maximizes the Evidence Lower Bound ℒ(q, θ).

1. Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1), 1–22. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x
2. Wu, C. F. J. (1983). On the convergence properties of the EM algorithm. The Annals of Statistics, 11(1), 95–103. https://doi.org/10.1214/aos/1176346060
3. Neal, R. M., & Hinton, G. E. (1998). A view of the EM algorithm that justifies incremental, sparse, and other variants. In M. I. Jordan (Ed.), Learning in Graphical Models (pp. 355–368). Springer. https://doi.org/10.1007/978-94-011-5014-9_12
4. Bailey, T. L., & Elkan, C. (1994). Fitting a mixture model by expectation maximization to discover motifs in biopolymers. Proceedings of the Second International Conference on Intelligent Systems for Molecular Biology (ISMB), 2, 28–36. https://doi.org/10.1093/bioinformatics/10.6.665