The EM Algorithm

A latent variable model assumes each observed data point $x$ is generated together with an unobserved variable $z$ (a cluster label, a hidden state, a missing feature). The joint density is $p(x, z \mid \theta)$, and the observed-data likelihood is obtained by marginalizing out $z$:

For a $K$-component Gaussian mixture, $z \in \{1,\dots,K\}$ indicates which component generated $x$, and $p(x \mid \theta) = \sum_{k=1}^K \pi_k \, \mathcal N(x \mid \mu_k, \Sigma_k)$. If we knew $z$ for every point, maximum likelihood estimation would be trivial: compute per-cluster means, covariances, and mixing weights in closed form. The difficulty is entirely due to the sum inside the logarithm of the log-likelihood,

which couples all components together and has no closed-form maximizer. EM sidesteps this by introducing, for each data point, a distribution $q(z_i)$ over the hidden variable and using Jensen's inequality to build a lower bound on $\ell(\theta)$ called the Evidence Lower Bound (ELBO):

The gap between the true log-likelihood and the ELBO is exactly the Kullback–Leibler divergence between $q(z_i)$ and the true posterior $p(z_i \mid x_i, \theta)$:

This identity is the crux of the whole method. It says the bound is tight exactly when $q(z_i) = p(z_i \mid x_i, \theta)$, and it is this fact that EM exploits by alternating two moves: set $q$ to the exact posterior (closing the gap — the E-step), then maximize the ELBO over $\theta$ with $q$ fixed (climbing the bound — the M-step). Because the ELBO touches $\ell(\theta)$ at the current parameters after the E-step, any improvement of the ELBO in the M-step is an improvement (or at worst no change) of the true log-likelihood too.

E-step (Expectation). Fix the current parameters $\theta^{(t)}$. For every data point $i$, compute the posterior responsibility of each latent value $z$:

For a Gaussian mixture, this is the familiar responsibility

No parameters are changed in this step; only beliefs about the hidden labels are updated. This closes the KL gap exactly (setting it to zero), so the ELBO now equals the true log-likelihood at $\theta^{(t)}$.

M-step (Maximization). With $q^{(t)}$ fixed, maximize the ELBO over $\theta$. Since $q$ does not depend on $\theta$ in this step, this reduces to maximizing the expected complete-data log-likelihood:

For a Gaussian mixture this maximization is closed-form: it is weighted maximum likelihood, where each point $x_i$ contributes to cluster $k$ with weight $\gamma_{ik}$:

These two steps are repeated until the log-likelihood stops improving (by more than a tolerance) or a maximum number of iterations is reached. ``` ALGORITHM EM-for-GMM(X, K, tol, max_iter) Input: data X = {x_1, ..., x_N}, number of components K Output: parameters theta = {pi_k, mu_k, Sigma_k}_{k=1..K} Initialize pi_k, mu_k, Sigma_k for k = 1..K // e.g. via k-means++ seeding prev_ll <- -infinity repeat // ---- E-step: compute responsibilities ---- for i = 1..N: for k = 1..K: num[i][k] <- pi_k Gaussian(x_i | mu_k, Sigma_k) denom <- sum_k num[i][k] for k = 1..K: gamma[i][k] <- num[i][k] / denom // posterior p(z_i=k | x_i, theta) // ---- M-step: re-estimate parameters ---- for k = 1..K: N_k <- sum_i gamma[i][k] pi_k <- N_k / N mu_k <- (1 / N_k) sum_i gamma[i][k] x_i Sigma_k <- (1 / N_k) sum_i gamma[i][k] (x_i - mu_k)(x_i - mu_k)^T // ---- check convergence via observed log-likelihood ---- ll <- sum_i log( sum_k pi_k Gaussian(x_i | mu_k, Sigma_k) ) if |ll - prev_ll| < tol: break prev_ll <- ll until max_iter reached return theta ``` Two edge cases the pseudocode must handle in practice: (1) if $N_k \to 0$ for some component (an empty cluster), $\Sigma_k$ is undefined — guard by re-seeding that component or adding a small floor to $N_k$; (2) if a component collapses onto a single point, $\Sigma_k \to 0$ and the Gaussian density diverges — regularize by adding $\epsilon I$ to every covariance estimate.

Time complexity per iteration (Gaussian mixture, $d$-dimensional data, $K$ components, $N$ points). The E-step evaluates $N \cdot K$ Gaussian densities; each density evaluation costs $O(d^2)$ for the Mahalanobis term when using a precomputed Cholesky factor of $\Sigma_k$ (or $O(d^3)$ once per component per iteration to factor $\Sigma_k$, amortized over the $N$ points). So the E-step costs

The M-step computes weighted sums for $\mu_k$ ($O(NKd)$) and weighted outer products for $\Sigma_k$ ($O(NKd^2)$), giving

Total per iteration: $O(NKd^2 + Kd^3)$, dominated by $O(NKd^2)$ when $N \gg d$, which is the typical regime. Space complexity is $O(NK)$ for the responsibility matrix $\gamma$ (or $O(K)$ if streamed one point at a time and parameters are accumulated online) plus $O(Kd^2)$ for the covariances. The number of iterations to convergence is not bounded by the problem size alone — it depends on how close the initialization is to a stationary point — but each individual iteration is cheap and the iteration count is usually small (tens to low hundreds) in practice.

Correctness — the monotonicity theorem. The central claim is:

The proof follows directly from the ELBO decomposition above. Write, for any $q$,

At the start of iteration $t$, the E-step sets $q^{(t)} = p(Z \mid X, \theta^{(t)})$, which makes the KL term exactly zero, so

The M-step then chooses $\theta^{(t+1)} = \arg\max_\theta \mathrm{ELBO}(q^{(t)}, \theta)$, so in particular

Finally, because the ELBO is a lower bound for any $q$ (Jensen's inequality applied to the concave $\log$), in particular for $q^{(t)}$ evaluated at the new parameters,

Chaining the two inequalities gives $\ell(\theta^{(t+1)}) \ge \ell(\theta^{(t)})$, proving monotone non-decrease. Since $\ell$ is bounded above for well-posed likelihood models (under regularity/regularization to avoid degenerate variance collapse), a monotone bounded sequence converges, so $\ell(\theta^{(t)})$ converges to some limit.

What this does not prove. Convergence of the likelihood sequence does not imply convergence to a global maximum, nor even that the parameter sequence $\theta^{(t)}$ converges (though under mild conditions it does, to a stationary point of $\ell$). EM is exactly coordinate ascent on the joint objective $F(q,\theta) = \mathrm{ELBO}(q,\theta)$: the E-step is exact coordinate maximization over $q$ (closed form, since the optimal $q$ is just the posterior), and the M-step is exact or partial coordinate maximization over $\theta$. Coordinate ascent on a function bounded above converges to a point where neither coordinate block can be locally improved — a stationary point, which for non-convex $\ell$ may be a local maximum or saddle. This is the fundamental tradeoff of EM: guaranteed monotonic progress, no guarantee of global optimality. In practice this is mitigated by multiple random restarts, sensible initialization (e.g. k-means++), and model-selection criteria like BIC to choose $K$.

Consider $N=6$ one-dimensional points, $K=2$ Gaussian components, and a deliberately asymmetric dataset so no accidental symmetry hides an error:

Initialize (deliberately imperfectly): $\pi_1=\pi_2=0.5$, $\mu_1 = 2$, $\mu_2 = 8$, $\sigma_1^2 = \sigma_2^2 = 4$.

Iteration 1 — E-step. Using $\mathcal N(x\mid\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\!\big(-\frac{(x-\mu)^2}{2\sigma^2}\big)$, compute unnormalized responsibilities $\pi_k \mathcal N(x_i\mid\mu_k,\sigma_k^2)$ and normalize:

(*The far outlier $x=25$ is roughly $10^{18}$ times more likely under component 2's density than component 1's at this scale, so it is assigned almost entirely — but not literally 100% — to component 2; the exact ratio still depends on the tails.)

Iteration 1 — M-step. Compute effective counts and weighted moments:

Notice how the single outlier $x=25$, softly assigned to component 2, pulls $\mu_2$ up sharply from the initial 8 to about 14.65 — a preview of EM's known sensitivity to outliers in Gaussian mixtures. New variances are computed the same weighted way, and new mixing weights are $\pi_1 = N_1/6 \approx 0.666$, $\pi_2 \approx 0.334$.

Subsequent iterations. As iterations continue, component 1 tightens around the tight cluster $\{1,2,3\}$, and component 2 stretches to cover $\{9,10,25\}$ with an inflated variance to accommodate the outlier — since a single Gaussian must explain both 9–10 and 25 simultaneously. The observed log-likelihood $\ell(\theta^{(t)}) = \sum_i \log\hat p(x_i)$ increases at every step (monotonicity theorem above); after a handful of iterations the parameters stabilize near $\mu_1\approx 2$, $\sigma_1^2$ small, $\mu_2 \approx 14.7$, $\sigma_2^2$ large, and $\pi_1 \approx 0.5$, matching the visually obvious two-cluster-plus-outlier structure of the data — while also exposing that this 2-component model is not a great fit for a true three-group structure (two tight clusters and one outlier), a diagnostic used heavily in Anomaly Detection Fundamentals.

The 1977 paper by Dempster, Laird, and Rubin unified a scattered family of ad hoc "fill in the missing data, re-estimate, repeat" procedures — already used piecemeal in genetics, factor analysis, and hidden Markov model estimation — into a single, formally justified algorithm with a general monotonicity proof. This transformed EM from a collection of tricks into a principled tool applicable to any exponential-family latent variable model.

EM (and its variational generalizations) is embedded in an enormous range of production systems. Speech recognition and bioinformatics pipelines historically relied on the Baum–Welch algorithm — EM specialized to hidden Markov models — to train phoneme and gene-sequence models. Gaussian mixture models fit with EM are used for speaker verification, background subtraction in video, and soft clustering in customer segmentation, directly building on Gaussian Mixture Models and Density Estimation Fundamentals. Recommender systems and survey-data pipelines use EM-style imputation to handle missing entries under a probabilistic model rather than discarding incomplete rows. In probabilistic topic modeling (e.g., variants of latent Dirichlet allocation), variational EM extends the same E-step/M-step alternation to approximate rather than exact posteriors when the exact posterior is intractable. Fitted mixture densities from EM feed directly into likelihood-based outlier scoring, the starting point for Anomaly Detection Fundamentals: points with very low $\hat p(x)$ under the learned mixture are flagged as anomalies, making EM not just a clustering tool but a density-estimation engine for downstream detection systems.

Ask the tutor: Why does adding more mixture components always increase training likelihood, and what does that imply about how you should choose K?