The EM Algorithm and Gaussian-Mixture / Latent-Variable Models

Intuition Beginner

Suppose you have a pile of measurements that came from two different sources mixed together — say, the heights of fish from two species netted in one haul — but the labels are gone. You want the average height of each species, yet to compute those averages you would need to know which fish belongs to which species, and to guess the species you would need the averages. Each piece of knowledge depends on the other. This chicken-and-egg loop is the situation the EM algorithm is built to break.

EM breaks it by guessing, then refining. Start with a rough guess for the two averages. Now, for every fish, ask: given these guessed averages, how likely is it that this fish came from species one versus species two? A short fish near the first average gets assigned mostly to species one; a fish sitting between the two averages gets split, say sixty percent to one and forty percent to the other. These soft, fractional memberships are called responsibilities. That is the first half of the cycle.

The second half: now that every fish carries a fractional vote for each species, recompute each average as a weighted average, where each fish counts according to its vote. The new averages are better than the guesses you started with. Feed them back in, recompute the votes, recompute the averages, and repeat. Each pass improves the fit, and the numbers settle down. You alternate between guessing the hidden labels and re-estimating the parameters, and the loop converges.

Visual Beginner

Picture a number line with a cloud of dots — the fish heights. Two bell curves sit on the line, one centred on each guessed average. For a given dot, the two bell curves report how tall they are at that spot; the taller curve at that point wins more of the dot's vote. A dot directly under the left bell goes almost entirely to the left species; a dot in the overlap zone between the bells splits its vote.

step	what you do	what comes out
E-step	with the bells fixed, score each dot's vote for each bell	a soft membership (responsibility) per dot
M-step	with the votes fixed, slide and reshape each bell	updated averages, spreads, and mixing fractions

The takeaway: EM is two alternating moves. One move fills in the missing labels softly using the current parameters; the other re-estimates the parameters using those soft labels. Repeating the pair drives the fit uphill until it stops improving.

Worked example Beginner

Take four measurements on a line: $1,2,8,9$. You believe they come from two clusters and you want each cluster's centre. Start with guessed centres $c_1=0$ and $c_2=10$.

Step 1, score the votes (E-step). For each point, give more vote to the nearer centre. A simple scoring is to make the vote for a centre larger when the point is closer. Point $1$ is near $c_1=0$ (distance $1$) and far from $c_2=10$ (distance $9$), so almost all of its vote goes to cluster one; say $0.9$ to cluster one and $0.1$ to cluster two. By symmetry point $2$ also votes about $0.9$ for cluster one. Points $8$ and $9$ are near $c_2=10$, so each votes about $0.9$ for cluster two and $0.1$ for cluster one.

Step 2, recompute the centres (M-step). Each new centre is the vote-weighted average of all four points. For cluster one, the weights are $0.9,0.9,0.1,0.1$ on the points $1,2,8,9$:

$$c_1=\frac{0.9\cdot 1+0.9\cdot 2+0.1\cdot 8+0.1\cdot 9}{0.9+0.9+0.1+0.1}=\frac{0.9+1.8+0.8+0.9}{2.0}=\frac{4.4}{2.0}=\mathrm{2.2.}$$

For cluster two the weights flip to $0.1,0.1,0.9,0.9$, giving $c_2=(0.1\cdot 1+0.1\cdot 2+0.9\cdot 8+0.9\cdot 9)/2.0=(0.1+0.2+7.2+8.1)/2.0=15.6/2.0=7.8$.

Step 3, repeat. The centres moved from $0$ and $10$ to $2.2$ and $7.8$ — much closer to the two natural clumps near $1.5$ and $8.5$. Feeding $2.2$ and $7.8$ back in sharpens the votes (points $1,2$ now vote even more strongly for cluster one), and the next M-step pulls the centres to roughly $1.5$ and $8.5$, where they stop moving.

What this tells us: starting from a crude guess, one E-step plus one M-step already produced sensible centres, and a second pass nearly finishes the job. The fractional votes are what let the two unknowns be estimated together instead of waiting on each other.

Check your understanding Beginner

Formal definition Intermediate+

Let $Y$ denote the observed data with density $p(y\mid \theta )$ for $\theta \in \Theta \subseteq {\mathbb{R}}^d$, and suppose this density arises by marginalising a latent (unobserved) variable $Z$ from a complete-data density $p(y,z\mid \theta )$:

$$p(y\mid \theta )=\int p(y,z\mid \theta )dz\text{(or }\sum _z\text{ for discrete }Z\text{)}.$$

The incomplete-data log-likelihood is $\ell (\theta )=\log p(y\mid \theta )$ and the complete-data log-likelihood is ${\ell }_c(\theta )=\log p(y,z\mid \theta )$; the latter cannot be evaluated because $z$ is not observed. Write $p(z\mid y,\theta )=p(y,z\mid \theta )/p(y\mid \theta )$ for the posterior of the latent variable.

Definition (Q-function). For a current parameter value ${\theta }^{\prime }$, the expected complete-data log-likelihood, or Q-function, is the conditional expectation 37.04.01

$$Q(\theta \mid {\theta }^{\prime })={\mathbb{E}}_{Z\mid Y,{\theta }^{\prime }}[\log p(Y,Z\mid \theta )|Y=y]=\int \log p(y,z\mid \theta )p(z\mid y,{\theta }^{\prime })dz.$$

The averaging is over the latent posterior at the current parameter ${\theta }^{\prime }$, while the argument $\theta$ inside the logarithm is the free variable to be optimised.

Definition (the EM iteration). Given an initial ${\theta }^{(0)}$, the EM algorithm produces a sequence ${\theta }^{(0)},{\theta }^{(1)},\dots$ by alternating:

• E-step. Using the current ${\theta }^{(k)}$, form the function $\theta \mapsto Q(\theta \mid {\theta }^{(k)})$ (in practice: compute the latent posterior $p(z\mid y,{\theta }^{(k)})$ and the expectations it induces).
• M-step. Maximise the surrogate, ${\theta }^{(k+1)}\in \arg {\max }_{\theta \in \Theta }Q(\theta \mid {\theta }^{(k)})$.

A generalised EM (GEM) algorithm relaxes the M-step to any update with $Q({\theta }^{(k+1)}\mid {\theta }^{(k)})\ge Q({\theta }^{(k)}\mid {\theta }^{(k)})$.

Definition (the evidence lower bound). For any probability density $q(z)$ over the latent space, the evidence lower bound (ELBO) is

$$\mathcal{L}(q,\theta )={\mathbb{E}}_q\left[\log \frac{p(y,Z\mid \theta )}{q(Z)}\right]=\log p(y\mid \theta )-D_{\mathrm{KL}}(q\Vert p(\cdot \mid y,\theta )),$$

where $D_{\mathrm{KL}}(q\Vert r)={\mathbb{E}}_q[\log (q(Z)/r(Z))]\ge 0$ is the Kullback-Leibler divergence. Because the KL term is non-negative and vanishes exactly when $q=p(\cdot \mid y,\theta )$, the ELBO is a lower bound on $\ell (\theta )$ that is tight at the latent posterior.

Definition (Gaussian mixture model). A $K$-component Gaussian mixture in ${\mathbb{R}}^p$ has density $p(y\mid \theta )={\sum }_{k=1}^K{\pi }_k\varphi (y;{\mu }_k,{\Sigma }_k)$, with mixing weights ${\pi }_k\ge 0$, ${\sum }_k{\pi }_k=1$, component means ${\mu }_k$, covariances ${\Sigma }_k$, and $\varphi (\cdot ;\mu ,\Sigma )$ the multivariate normal density. The latent $Z\in \{1,\dots ,K\}$ is the component label, $p(Z=k\mid \theta )={\pi }_k$ and $p(y\mid Z=k,\theta )=\varphi (y;{\mu }_k,{\Sigma }_k)$. The responsibility of component $k$ for observation $y_i$ at parameter ${\theta }^{\prime }$ is the posterior class probability

$${\gamma }_{ik}=p(Z_i=k\mid y_i,{\theta }^{\prime })=\frac{{\pi }_k^{\prime }\varphi (y_i;{\mu }_k^{\prime },{\Sigma }_k^{\prime })}{{\sum }_{j=1}^K{\pi }_j^{\prime }\varphi (y_i;{\mu }_j^{\prime },{\Sigma }_j^{\prime })}.$$

The new symbols are recorded in _meta/NOTATION.md: $Q(\theta \mid {\theta }^{\prime })$ (Q-function), $\ell (\theta )$ (incomplete-data log-likelihood), ${\ell }_c(\theta )$ (complete-data log-likelihood), $\mathcal{L}(q,\theta )$ (ELBO), $D_{\mathrm{KL}}$ (Kullback-Leibler divergence), ${\gamma }_{ik}$ (responsibility).

Counterexamples to common slips Intermediate+

• "EM maximises the Q-function, so it maximises the likelihood." The M-step maximises $Q(\cdot \mid {\theta }^{(k)})$, a surrogate tangent to $\ell$ at ${\theta }^{(k)}$, not $\ell$ itself. Each step increases $\ell$, but a single M-step lands at the surrogate's peak, which is generally not $\ell$'s peak; the guarantee is monotone ascent, not one-shot maximisation.

• "EM converges to the global maximum likelihood estimate." The mixture log-likelihood is non-concave and multimodal. EM ascends to a stationary point — a local maximum or a saddle — determined by the starting value. Different initialisations reach different solutions, which is why multiple random restarts (or k-means initialisation) are standard.

• "The responsibilities are hard 0/1 labels." They are posterior probabilities in $[0,1]$ summing to one across components. A point exactly between two equal-variance components splits its responsibility evenly; hard assignment is the degenerate limit recovered only as the component variances shrink to zero.

• "The mixture likelihood has a well-defined global maximum to find." For Gaussian mixtures with free covariances the likelihood is unbounded: placing one component's mean on a single data point and shrinking its variance to zero sends the likelihood to $+\infty$. The meaningful target is a sensible interior local maximum, not the (degenerate) global supremum.

Key theorem with proof Intermediate+

The signature result is the monotonicity theorem: each EM iteration never decreases the observed-data log-likelihood. The proof exhibits $Q$ as the variable part of a lower bound on $\ell$ and isolates the gap as a Kullback-Leibler divergence, which the conditional form of Jensen's inequality 37.04.01 forces to move the right way ^{[Dempster, A. P., Laird, N. M. & Rubin, D. B. — Maximum Likelihood from Incomplete Data via the EM Algorithm]}.

Theorem (EM ascent / monotonicity). Let ${\theta }^{(k+1)}$ satisfy $Q({\theta }^{(k+1)}\mid {\theta }^{(k)})\ge Q({\theta }^{(k)}\mid {\theta }^{(k)})$ (in particular if ${\theta }^{(k+1)}\in \arg {\max }_{\theta }Q(\theta \mid {\theta }^{(k)})$). Then the observed-data log-likelihood does not decrease,

$$\ell ({\theta }^{(k+1)})\ge \ell ({\theta }^{(k)}),$$

with equality only if $Q({\theta }^{(k+1)}\mid {\theta }^{(k)})=Q({\theta }^{(k)}\mid {\theta }^{(k)})$ and $p(z\mid y,{\theta }^{(k+1)})=p(z\mid y,{\theta }^{(k)})$ almost everywhere.

Proof. Take logarithms of the identity $p(y,z\mid \theta )=p(z\mid y,\theta )p(y\mid \theta )$ and rearrange:

$$\log p(y\mid \theta )=\log p(y,z\mid \theta )-\log p(z\mid y,\theta ).$$

The left side is free of $z$. Take the expectation of both sides over $z\sim p(\cdot \mid y,{\theta }^{(k)})$; the left side passes through unchanged, giving the decomposition

$$\ell (\theta )=Q(\theta \mid {\theta }^{(k)})-H(\theta \mid {\theta }^{(k)}),H(\theta \mid {\theta }^{(k)})={\mathbb{E}}_{Z\mid y,{\theta }^{(k)}}[\log p(Z\mid y,\theta )].$$

Evaluate this at $\theta ={\theta }^{(k+1)}$ and at $\theta ={\theta }^{(k)}$ and subtract:

$$\ell ({\theta }^{(k+1)})-\ell ({\theta }^{(k)})=\underset{\ge 0\text{ by the M-step}}{\underbrace{[Q({\theta }^{(k+1)}\mid {\theta }^{(k)})-Q({\theta }^{(k)}\mid {\theta }^{(k)})]}}-\underset{=:\Delta H}{\underbrace{[H({\theta }^{(k+1)}\mid {\theta }^{(k)})-H({\theta }^{(k)}\mid {\theta }^{(k)})]}}.$$

It remains to show $\Delta H\le 0$, i.e. $H$ is maximised in its first argument at ${\theta }^{(k)}$. Write $r(z)=p(z\mid y,{\theta }^{(k)})$ and $s(z)=p(z\mid y,{\theta }^{(k+1)})$. Then

$$H({\theta }^{(k+1)}\mid {\theta }^{(k)})-H({\theta }^{(k)}\mid {\theta }^{(k)})={\mathbb{E}}_{Z\sim r}\left[\log \frac{s(Z)}{r(Z)}\right]=-D_{\mathrm{KL}}(r\Vert s)\le 0,$$

where the inequality is the Gibbs inequality, itself the conditional Jensen inequality 37.04.01 applied to the concave logarithm: ${\mathbb{E}}_r[\log (s/r)]\le \log {\mathbb{E}}_r[s/r]=\log \int s(z)dz=\log 1=0$. Therefore both bracketed terms in the difference push $\ell$ upward, so $\ell ({\theta }^{(k+1)})-\ell ({\theta }^{(k)})\ge 0$. Equality forces both the $Q$-difference to vanish and $D_{\mathrm{KL}}(r\Vert s)=0$, i.e. $r=s$ almost everywhere. $\square$

Bridge. This monotonicity is the foundational reason EM is a usable algorithm despite the latent-variable likelihood having no closed-form maximiser, and it builds toward the variational reformulation: the decomposition $\ell (\theta )=Q(\theta \mid {\theta }^{\prime })-H(\theta \mid {\theta }^{\prime })$ is exactly the ELBO identity $\ell (\theta )=\mathcal{L}(q,\theta )+D_{\mathrm{KL}}(q\Vert p(\cdot \mid y,\theta ))$ read with $q=p(\cdot \mid y,{\theta }^{\prime })$, so EM is coordinate ascent on $\mathcal{L}(q,\theta )$ — the E-step maximises over $q$ (setting the KL gap to zero) and the M-step maximises over $\theta$. This is dual to the Newton/Fisher-scoring route of direct maximum likelihood 44.03.03, which climbs $\ell$ directly using its gradient and curvature; EM trades that curvature information for a sequence of easy complete-data maximisations. The central insight is that the awkward marginal likelihood is replaced by a surrogate that touches it from below and is maximised separately, and this is exactly the minorize-maximize principle, of which EM is the canonical probabilistic instance. Putting these together, the same identity reappears in the Gaussian-mixture M-step, where maximising $Q$ becomes responsibility-weighted moment matching, the latent-data echo of the exponential-family likelihood equation ${\mathbb{E}}_{\hat{\theta }}[T]=\bar{T}$ from 45.01.04; the convergence behaviour this guarantees — ascent to a stationary point, not the global optimum — appears again in the analysis of the algorithm's fixed points and rate in the Advanced results.

Exercises Intermediate+

Advanced results Master

The results below organise the EM algorithm around the variational identity of the Key theorem, the fixed-point and convergence theory that the monotonicity guarantee only opens, the rate of convergence governed by missing information, and the canonical instances — Gaussian mixtures, the soft/hard k-means limit, and the hidden-Markov Baum-Welch recursion. Direct gradient and Newton maximisation of $\ell$ remain the contrast class throughout 44.03.03, 45.01.04.

Theorem 1 (EM as coordinate ascent on the ELBO). Define $F(q,\theta )=\mathcal{L}(q,\theta )=\ell (\theta )-D_{\mathrm{KL}}(q\Vert p(\cdot \mid y,\theta ))$. The E-step is the exact maximisation $q^{(k+1)}=\arg {\max }_{q}F(q,{\theta }^{(k)})=p(\cdot \mid y,{\theta }^{(k)})$, attaining $F(q^{(k+1)},{\theta }^{(k)})=\ell ({\theta }^{(k)})$; the M-step is ${\theta }^{(k+1)}=\arg {\max }_{\theta }F(q^{(k+1)},\theta )$, which (dropping the $\theta$-free entropy of $q^{(k+1)}$) equals $\arg {\max }_{\theta }Q(\theta \mid {\theta }^{(k)})$. EM is therefore alternating maximisation of a single objective over two coordinate blocks, and monotone ascent of $\ell$ is the immediate corollary $\ell ({\theta }^{(k+1)})\ge F(q^{(k+1)},{\theta }^{(k+1)})\ge F(q^{(k+1)},{\theta }^{(k)})=\ell ({\theta }^{(k)})$. This is the foundation of the variational generalisations in which the E-step is restricted to a tractable family $\mathcal{Q}$, giving variational EM with an inexact but still monotone bound ^{[Hastie, T., Tibshirani, R. & Friedman, J. — The Elements of Statistical Learning (2nd ed.)]}.

Theorem 2 (Wu's convergence theorem). Suppose $\ell$ is continuous on $\Theta$ and the Q-function $Q(\theta \mid {\theta }^{\prime })$ is continuous in both arguments. If the EM map $M$ is closed and the level set $\{\theta :\ell (\theta )\ge \ell ({\theta }^{(0)})\}$ is compact, then all limit points of the EM sequence $\{{\theta }^{(k)}\}$ are stationary points of $\ell$, and $\ell ({\theta }^{(k)})$ increases monotonically to a value ${\ell }^{\ast }=\ell ({\theta }^{\ast })$ for some stationary point ${\theta }^{\ast }$. Convergence of the parameter sequence itself (not merely of its limit-point set or of $\ell$) requires the additional condition that the stationary points achieving ${\ell }^{\ast }$ are isolated; without it the iterates can in principle drift along a ridge of equal likelihood. Dempster, Laird, and Rubin's original convergence assertion overstated the conclusion; Wu's analysis supplied the correct hypotheses ^{[McLachlan, G. J. & Krishnan, T. — The EM Algorithm and Extensions (2nd ed.)]}.

Theorem 3 (linear rate and the fraction of missing information). Near a fixed point ${\theta }^{\ast }$ the EM map $M$ has Jacobian $DM({\theta }^{\ast })=I_{\mathrm{com}}({\theta }^{\ast }{)}^{-1}{I}_{\mathrm{mis}}({\theta }^{\ast })=I-I_{\mathrm{com}}^{-1}{I}_{\mathrm{obs}}$, where $I_{\mathrm{com}}$ is the expected complete-data information, $I_{\mathrm{obs}}$ the observed-data information, and $I_{\mathrm{mis}}=I_{\mathrm{com}}-I_{\mathrm{obs}}$ the missing information. Consequently ${\theta }^{(k+1)}-{\theta }^{\ast }\approx DM({\theta }^{\ast })({\theta }^{(k)}-{\theta }^{\ast })$, so EM converges linearly (geometrically) with rate equal to the largest eigenvalue of $DM({\theta }^{\ast })$ — the largest fraction of missing information. When little information is missing the rate is fast; when the latent structure is highly informative (heavily overlapping mixture components) EM crawls, which motivates acceleration via Aitken extrapolation, ECME, or switching to a Newton step on $\ell$ near the optimum ^{[McLachlan, G. J. & Krishnan, T. — The EM Algorithm and Extensions (2nd ed.)]}.

Theorem 4 (Gaussian mixtures: closed-form M-step and the k-means limit). For the $K$-component Gaussian mixture the M-step has the closed form ${\hat{\pi }}_k=N_k/n$, ${\hat{\mu }}_k=N_k^{-1}{\sum }_i{\gamma }_{ik}{y}_i$, ${\hat{\Sigma }}_k=N_k^{-1}{\sum }_i{\gamma }_{ik}(y_i-{\hat{\mu }}_k)(y_i-{\hat{\mu }}_k{)}^{\top }$ with $N_k={\sum }_i{\gamma }_{ik}$, i.e. responsibility-weighted means and covariances, the soft analogue of the cluster statistics in 45.08.08. Constrain all covariances to ${\Sigma }_k={\sigma }^{2}I$ and let ${\sigma }^2\to 0$: the responsibilities ${\gamma }_{ik}={\pi }_k\exp (-\Vert y_i-{\mu }_k{\Vert }^2/2{\sigma }^2)/{\sum }_j{\pi }_j\exp (-\Vert y_i-{\mu }_j{\Vert }^2/2{\sigma }^2)$ converge to the hard indicator $1\{k=\arg {\min }_j\Vert y_i-{\mu }_j{\Vert }^2\}$, the M-step means become unweighted cluster averages, and EM degenerates exactly into Lloyd's k-means algorithm. Soft k-means is EM for the spherical-Gaussian mixture at finite ${\sigma }^2$; hard k-means is its zero-temperature limit ^{[Hastie, T., Tibshirani, R. & Friedman, J. — The Elements of Statistical Learning (2nd ed.)]}.

Theorem 5 (Baum-Welch as EM for hidden Markov models). For a hidden Markov model with latent state path $Z_{1:T}$, observations $Y_{1:T}$, initial distribution $\rho$, transition matrix $A$, and emission densities $b_k$, the complete-data log-likelihood is additive over time, $\log p(y_{1:T},z_{1:T})=\log {\rho }_{z_1}+{\sum }_t\log A_{z_{t-1}{z}_t}+{\sum }_t\log b_{z_t}(y_t)$. The E-step requires only the single-state posteriors ${\gamma }_t(k)=p(Z_t=k\mid y_{1:T})$ and pairwise posteriors ${\xi }_t(j,k)=p(Z_{t-1}=j,Z_t=k\mid y_{1:T})$, computed in $O(T{K}^2)$ by the forward-backward recursion; the M-step then updates $\rho$, $A$, and the emission parameters as posterior-expected-count frequencies. The resulting algorithm is the Baum-Welch algorithm, which predates and is a special case of the general EM formulation ^{[McLachlan, G. J. & Krishnan, T. — The EM Algorithm and Extensions (2nd ed.)]}.

Synthesis. The EM algorithm is one identity worked in several settings: the decomposition $\ell (\theta )=Q(\theta \mid {\theta }^{\prime })-H(\theta \mid {\theta }^{\prime })$, equivalently $\ell (\theta )=\mathcal{L}(q,\theta )+D_{\mathrm{KL}}(q\Vert p(\cdot \mid y,\theta ))$, is the foundational reason monotone ascent holds and the reason EM is exactly coordinate ascent on the ELBO. The central insight is that an intractable marginal likelihood is replaced by a surrogate that touches it from below at the current iterate and is maximised separately — the minorize-maximize principle — so the latent integral is handled once, in the E-step, and the M-step is an ordinary complete-data fit. This is exactly the latent-data generalisation of 45.01.04: for exponential-family complete data the M-step is moment matching with the conditional expectation of the sufficient statistic in place of the unobserved statistic, so EM and maximum likelihood share one stationarity condition. The construction is dual to direct gradient/Newton maximisation of $\ell$ 44.03.03: EM buys stable, constraint-respecting, derivative-free ascent at the price of only linear convergence, the rate set by the fraction of missing information of Theorem 3, while Newton converges quadratically but must contend with the non-concave, possibly unbounded mixture surface. Putting these together, the Gaussian-mixture responsibilities specialise to hard k-means 45.08.08 in the vanishing-variance limit, and the forward-backward E-step specialises EM to Baum-Welch; the same surrogate-ascent skeleton recurs each time, and the bridge forward is that restricting the E-step to a tractable family is precisely the variational-inference programme that scales EM to models where the latent posterior is itself intractable.

Full proof set Master

Proposition 1 (the H-function inequality). For any $\theta$ and fixed ${\theta }^{\prime }$, $H(\theta \mid {\theta }^{\prime })\le H({\theta }^{\prime }\mid {\theta }^{\prime })$, with equality iff $p(\cdot \mid y,\theta )=p(\cdot \mid y,{\theta }^{\prime })$ almost everywhere, where $H(\theta \mid {\theta }^{\prime })={\mathbb{E}}_{Z\mid y,{\theta }^{\prime }}[\log p(Z\mid y,\theta )]$.

Proof. Write $r(z)=p(z\mid y,{\theta }^{\prime })$ and $t(z)=p(z\mid y,\theta )$, both probability densities in $z$. Then

$$H(\theta \mid {\theta }^{\prime })-H({\theta }^{\prime }\mid {\theta }^{\prime })={\mathbb{E}}_{Z\sim r}\left[\log \frac{t(Z)}{r(Z)}\right]=-D_{\mathrm{KL}}(r\Vert t)\le 0,$$

the inequality being the Gibbs inequality. It is the conditional Jensen inequality 37.04.01 for the strictly concave $\log$: ${\mathbb{E}}_r[\log (t/r)]\le \log {\mathbb{E}}_r[t/r]=\log \int t(z)dz=0$, with equality iff $t/r$ is constant $r$-almost-surely, i.e. $t=r$ a.e. since both integrate to one. $\square$

Proposition 2 (monotonicity / GEM ascent). If $Q({\theta }^{(k+1)}\mid {\theta }^{(k)})\ge Q({\theta }^{(k)}\mid {\theta }^{(k)})$ then $\ell ({\theta }^{(k+1)})\ge \ell ({\theta }^{(k)})$.

Proof. The identity $\ell (\theta )=Q(\theta \mid {\theta }^{(k)})-H(\theta \mid {\theta }^{(k)})$ holds for every $\theta$ (take the conditional expectation under $p(\cdot \mid y,{\theta }^{(k)})$ of $\log p(y\mid \theta )=\log p(y,z\mid \theta )-\log p(z\mid y,\theta )$, whose left side is $z$-free). Subtracting the instances at ${\theta }^{(k+1)}$ and ${\theta }^{(k)}$,

$$\ell ({\theta }^{(k+1)})-\ell ({\theta }^{(k)})=[Q({\theta }^{(k+1)}\mid {\theta }^{(k)})-Q({\theta }^{(k)}\mid {\theta }^{(k)})]+[H({\theta }^{(k)}\mid {\theta }^{(k)})-H({\theta }^{(k+1)}\mid {\theta }^{(k)})].$$

The first bracket is $\ge 0$ by hypothesis; the second is $\ge 0$ by Proposition 1. $\square$

Proposition 3 (Gaussian-mixture M-step). Maximising $Q(\theta \mid {\theta }^{\prime })={\sum }_{i=1}^n{\sum }_{k=1}^K{\gamma }_{ik}[\log {\pi }_k+\log \varphi (y_i;{\mu }_k,{\Sigma }_k)]$ over $\{{\pi }_k\ge 0,{\sum }_k{\pi }_k=1\}$ and $\{{\mu }_k,{\Sigma }_k\succ 0\}$ yields ${\hat{\pi }}_k=N_k/n$, ${\hat{\mu }}_k=N_k^{-1}{\sum }_i{\gamma }_{ik}{y}_i$, ${\hat{\Sigma }}_k=N_k^{-1}{\sum }_i{\gamma }_{ik}(y_i-{\hat{\mu }}_k)(y_i-{\hat{\mu }}_k{)}^{\top }$, $N_k={\sum }_i{\gamma }_{ik}$.

Proof. The objective separates: the $\pi$-part is ${\sum }_k{N}_k\log {\pi }_k$ and the $(\mu ,\Sigma )$-part splits across $k$. For $\pi$, form the Lagrangian ${\sum }_k{N}_k\log {\pi }_k-\lambda ({\sum }_k{\pi }_k-1)$; stationarity gives $N_k/{\pi }_k=\lambda$, and ${\sum }_k{\pi }_k=1$ with ${\sum }_k{N}_k=n$ forces $\lambda =n$, so ${\hat{\pi }}_k=N_k/n$. For component $k$, the relevant terms are ${\sum }_i{\gamma }_{ik}[-\frac{1}{2}\log |{\Sigma }_k|-\frac{1}{2}(y_i-{\mu }_k{)}^{\top }{\Sigma }_k^{-1}(y_i-{\mu }_k)]$. Differentiating in ${\mu }_k$: ${\Sigma }_k^{-1}{\sum }_i{\gamma }_{ik}(y_i-{\mu }_k)=0$ gives ${\hat{\mu }}_k=N_k^{-1}{\sum }_i{\gamma }_{ik}{y}_i$. Differentiating in ${\Sigma }_k^{-1}$ using $\partial \log |{\Sigma }_k|/\partial {\Sigma }_k^{-1}=-{\Sigma }_k$ and $\partial \mathrm{tr}({\Sigma }_k^{-1}{S}_i)/\partial {\Sigma }_k^{-1}=S_i$ with $S_i=(y_i-{\mu }_k)(y_i-{\mu }_k{)}^{\top }$: ${\sum }_i{\gamma }_{ik}({\Sigma }_k-S_i)=0$, so ${\hat{\Sigma }}_k=N_k^{-1}{\sum }_i{\gamma }_{ik}{S}_i$ at ${\mu }_k={\hat{\mu }}_k$. Each is the weighted Gaussian maximum-likelihood estimate. $\square$

Proposition 4 (hard k-means as the ${\sigma }^2\to 0$ limit). For the spherical mixture ${\Sigma }_k={\sigma }^{2}I$ with ${\pi }_k=1/K$, as ${\sigma }^2\to 0^{+}$ the responsibility ${\gamma }_{ik}\to 1\{k=k_i^{\ast }\}$ where $k_i^{\ast }=\arg {\min }_j\Vert y_i-{\mu }_j{\Vert }^2$ (assuming a unique nearest centre), and the M-step mean updates converge to the unweighted cluster means of k-means.

Proof. With equal weights and ${\Sigma }_k={\sigma }^{2}I$,

$${\gamma }_{ik}=\frac{\exp (-\Vert y_i-{\mu }_k{\Vert }^2/2{\sigma }^2)}{{\sum }_j\exp (-\Vert y_i-{\mu }_j{\Vert }^2/2{\sigma }^2)}=\frac{1}{{\sum }_j\exp ((\Vert y_i-{\mu }_k{\Vert }^2-\Vert y_i-{\mu }_j{\Vert }^2)/2{\sigma }^2)}.$$

Let $d_{ij}=\Vert y_i-{\mu }_j{\Vert }^2$ and $k_i^{\ast }=\arg {\min }_j{d}_{ij}$, unique by assumption. For $k=k_i^{\ast }$ every exponent $(d_{ik}-d_{ij})/2{\sigma }^2\le 0$ and the $j=k$ term is $1$, while every $j\ne k_i^{\ast }$ term has $d_{ij}>d_{ik}$, so its exponent $\to -\infty$ and the term $\to 0$; hence ${\gamma }_{i{k}_i^{\ast }}\to 1$. For $k\ne k_i^{\ast }$ the term with $j=k_i^{\ast }$ contributes $\exp ((d_{ik}-d_{i{k}_i^{\ast }})/2{\sigma }^2)\to +\infty$ in the denominator, so ${\gamma }_{ik}\to 0$. The responsibilities become hard indicators; the M-step mean $N_k^{-1}{\sum }_i{\gamma }_{ik}{y}_i$ becomes the average of the points assigned to centre $k$, which is the Lloyd update. $\square$

Proposition 5 (linear convergence rate of EM). Let $M$ be the EM map with fixed point ${\theta }^{\ast }$, $M({\theta }^{\ast })={\theta }^{\ast }$, and suppose $M$ is differentiable near ${\theta }^{\ast }$. Then $DM({\theta }^{\ast })=I_{\mathrm{com}}({\theta }^{\ast }{)}^{-1}{I}_{\mathrm{mis}}({\theta }^{\ast })$, and the iteration converges linearly with rate $\rho ={\lambda }_{\max }(DM({\theta }^{\ast }))\in [0,1)$ when ${\theta }^{\ast }$ is an isolated local maximum.

Proof (structure). Differentiating the self-consistency relation that defines the M-step, ${\nabla }_{10}Q(M(\theta )\mid \theta )=0$ (the gradient of $Q$ in its first argument vanishes at the maximiser), with respect to $\theta$ gives ${\nabla }_{20}QDM(\theta )+{\nabla }_{11}Q=0$ where the subscripts denote the order of differentiation in the first and second arguments. At ${\theta }^{\ast }$ one identifies $-{\nabla }_{20}Q({\theta }^{\ast }\mid {\theta }^{\ast })=I_{\mathrm{com}}({\theta }^{\ast })$ (the complete-data information) and ${\nabla }_{11}Q({\theta }^{\ast }\mid {\theta }^{\ast })=I_{\mathrm{mis}}({\theta }^{\ast })$ (the missing information, from differentiating the H-function relation $\ell =Q-H$ and using ${\nabla }^2\ell ({\theta }^{\ast })=-I_{\mathrm{obs}}=-(I_{\mathrm{com}}-I_{\mathrm{mis}})$). Hence $DM({\theta }^{\ast })=I_{\mathrm{com}}^{-1}{I}_{\mathrm{mis}}=I-I_{\mathrm{com}}^{-1}{I}_{\mathrm{obs}}$. Because $0⪯I_{\mathrm{mis}}\prec I_{\mathrm{com}}$ at an isolated local maximum, the eigenvalues of $DM({\theta }^{\ast })$ lie in $[0,1)$, and the linearisation ${\theta }^{(k+1)}-{\theta }^{\ast }\approx DM({\theta }^{\ast })({\theta }^{(k)}-{\theta }^{\ast })$ gives geometric convergence at rate ${\lambda }_{\max }(DM({\theta }^{\ast }))$. $\square$

Connections Master

• Maximum likelihood and the method of moments 45.01.04 is the parent estimation theory: EM is the device that computes the maximum likelihood estimate when the likelihood marginalises a latent variable and the score equations have no closed form. The exponential-family M-step is exactly that unit's likelihood equation ${\mathbb{E}}_{\theta }[T]=\bar{T}$ with the unobserved sufficient statistic replaced by its conditional expectation, so EM inherits the moment-matching characterisation and the strict-concavity-of-the-complete-data-likelihood structure that makes each M-step well posed, while the observed-data likelihood it ascends is in general non-concave.

• Newton's method for optimization 44.03.03 is the contrast and the accelerator: direct maximisation of the observed-data log-likelihood by Newton or Fisher scoring converges quadratically but must handle the non-concave, possibly unbounded mixture surface and supply the Hessian; EM trades that for derivative-free, monotone, linear-rate ascent. The two are routinely hybridised — EM from a crude start to enter a good basin, then a Newton step near the optimum where the fraction-of-missing-information rate of EM (Theorem 3) becomes the bottleneck.

• k-means and prototype clustering 45.08.08 is the hard-assignment specialisation: k-means is EM for an equal-weight spherical-Gaussian mixture in the limit of vanishing component variance, where soft responsibilities collapse to nearest-centre indicators and the responsibility-weighted M-step means become plain cluster averages. The soft/hard duality explains both why k-means is fast (no probabilities to track) and why it is brittle to non-spherical or unequal-size clusters that the full mixture handles.

• Conditional expectation and martingales 37.04.01 supplies the analytic engine: the Q-function is a conditional expectation of the complete-data log-likelihood, and the monotonicity proof rests on conditional Jensen / the non-negativity of Kullback-Leibler divergence applied to latent posteriors. The same conditional-expectation calculus reappears in the forward-backward E-step for hidden Markov models and in the variational-inference generalisations where the E-step is projected onto a tractable family.

Historical & philosophical context Master

Special cases of the alternate-and-fill-in idea long predate the general formulation: Hartley treated the discrete-count missing-data case in 1958 ^{[Hartley 1958]}, Baum and colleagues developed the forward-backward re-estimation for probabilistic functions of Markov chains through the late 1960s and early 1970s ^{[Baum 1970]}, and the same recursion arose independently in genetics and in factor analysis. The unifying recognition that these were one algorithm — and the name "EM" — is due to Arthur Dempster, Nan Laird, and Donald Rubin in their 1977 Journal of the Royal Statistical Society B paper, which gave the general incomplete-data formulation, the E- and M-steps, the GEM relaxation, and the monotone-ascent argument via the $H$-function and Jensen's inequality ^{[Dempster, A. P., Laird, N. M. & Rubin, D. B. — Maximum Likelihood from Incomplete Data via the EM Algorithm]}.

The 1977 convergence claim was incomplete: it asserted convergence to a maximiser where only convergence of the likelihood values and stationarity of limit points hold under regularity. C. F. Jeff Wu corrected the theory in 1983, separating convergence of $\ell ({\theta }^{(k)})$ from convergence of the parameter sequence and giving the conditions (closedness of the EM map, compact level sets, isolated stationary points) under which each holds ^{[Wu 1983]}. The book-length synthesis, including the missing-information rate, standard-error estimation by Louis's method and the supplemented-EM algorithm, and the catalogue of extensions (ECM, ECME, stochastic and Monte Carlo EM), is McLachlan and Krishnan's monograph ^{[McLachlan, G. J. & Krishnan, T. — The EM Algorithm and Extensions (2nd ed.)]}. The variational reading of EM as coordinate ascent on the evidence lower bound, which became central to scalable Bayesian computation, is articulated by Neal and Hinton in 1998 ^{[Neal-Hinton 1998]}.

Bibliography Master

@article{dempsterlairdrubin1977,
  author  = {Dempster, Arthur P. and Laird, Nan M. and Rubin, Donald B.},
  title   = {Maximum Likelihood from Incomplete Data via the {EM} Algorithm},
  journal = {Journal of the Royal Statistical Society, Series B},
  volume  = {39},
  number  = {1},
  year    = {1977},
  pages   = {1--38}
}

@book{mclachlankrishnan2008,
  author    = {McLachlan, Geoffrey J. and Krishnan, Thriyambakam},
  title     = {The {EM} Algorithm and Extensions},
  edition   = {2},
  series    = {Wiley Series in Probability and Statistics},
  publisher = {Wiley},
  year      = {2008}
}

@book{hastietibshiranifriedman2009,
  author    = {Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome},
  title     = {The Elements of Statistical Learning: Data Mining, Inference, and Prediction},
  edition   = {2},
  series    = {Springer Series in Statistics},
  publisher = {Springer},
  year      = {2009}
}

@article{wu1983,
  author  = {Wu, C. F. Jeff},
  title   = {On the Convergence Properties of the {EM} Algorithm},
  journal = {The Annals of Statistics},
  volume  = {11},
  number  = {1},
  year    = {1983},
  pages   = {95--103}
}

@article{baum1970,
  author  = {Baum, Leonard E. and Petrie, Ted and Soules, George and Weiss, Norman},
  title   = {A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of {Markov} Chains},
  journal = {The Annals of Mathematical Statistics},
  volume  = {41},
  number  = {1},
  year    = {1970},
  pages   = {164--171}
}

@article{hartley1958,
  author  = {Hartley, Herman O.},
  title   = {Maximum Likelihood Estimation from Incomplete Data},
  journal = {Biometrics},
  volume  = {14},
  number  = {2},
  year    = {1958},
  pages   = {174--194}
}

@incollection{nealhinton1998,
  author    = {Neal, Radford M. and Hinton, Geoffrey E.},
  title     = {A View of the {EM} Algorithm that Justifies Incremental, Sparse, and Other Variants},
  booktitle = {Learning in Graphical Models},
  editor    = {Jordan, Michael I.},
  publisher = {Kluwer Academic Publishers},
  year      = {1998},
  pages     = {355--368}
}

@book{bishop2006,
  author    = {Bishop, Christopher M.},
  title     = {Pattern Recognition and Machine Learning},
  series    = {Information Science and Statistics},
  publisher = {Springer},
  year      = {2006}
}