EM Algorithm

EM Algorithm

The Expectation-Maximization (EM) algorithm is a deterministic iterative procedure for computing maximum-likelihood (MLE) or maximum-a-posteriori (MAP) estimates in models with latent variables $z$, whose marginal likelihood $p(y\mid \theta )=\int p(y,z\mid \theta )dz$ is intractable but whose complete-data likelihood $p(y,z\mid \theta )$ has a tractable form. EM alternates two steps:

1. E-step (Expectation). Given the current parameter estimate ${\theta }^{(t)}$, compute the responsibilities — the posterior of the latent variables under the current parameters:
   $q^{(t)}(z):=p(z\mid y,{\theta }^{(t)}).$
   Then form the Q-function (expected complete-data log-likelihood):
   $Q(\theta \mid {\theta }^{(t)}):={\mathbb{E}}_{z\sim q^{(t)}}[\log p(y,z\mid \theta )].$
2. M-step (Maximization). Update the parameters by maximizing $Q$:
   ${\theta }^{(t+1)}:=\arg {\max }_{\theta }Q(\theta \mid {\theta }^{(t)}).$

EM guarantees monotonic improvement of the observed-data log-likelihood: $\log p(y\mid {\theta }^{(t+1)})\ge \log p(y\mid {\theta }^{(t)})$ at every step. Under mild regularity conditions the iterates converge to a stationary point (local max, saddle, or local min), so multiple restarts are standard practice. EM is not a Bayesian posterior sampler — it produces a point estimate, not samples — but it sits adjacent to Bayesian methods both historically and algorithmically.

Relation to Gibbs sampling. EM is the deterministic analog of Gibbs for the same latent-variable augmentation. Gibbs draws $z^{(t)}\sim p(z\mid y,{\theta }^{(t)})$ and ${\theta }^{(t+1)}\sim p(\theta \mid y,z^{(t)})$ and produces samples from the joint posterior; EM replaces the Gibbs draw of $z$ with its expectation (the E-step) and replaces the Gibbs draw of $\theta$ with an argmax (the M-step), producing the MAP estimate instead of posterior samples. The same latent-variable augmentation and the same full conditionals power both algorithms — they differ only in “sample” vs. “replace with conditional expectation / argmax.”

Canonical application: Gaussian mixture models. For $y_i\stackrel{\text{iid}}{\sim }{\sum }_{j=1}^k{w}_j\mathcal{N}({\mu }_j,{\sigma }_j^2)$, the latent $z_i$ is the component assignment. The E-step computes responsibility $r_{ij}^{(t)}=\mathbb{P}[z_i=j\mid y_i,{\theta }^{(t)}]$; the M-step updates $w_j,{\mu }_j,{\sigma }_j^2$ to their responsibility-weighted MLEs. The algorithm reduces to the soft-assignment version of Lloyd’s $k$-means.

Prerequisites

• 3.11.2 Gibbs Sampling — Gibbs samples from the same latent-variable augmentation; EM replaces the draws with expectations/argmaxes.
• 3.11 MCMC — parent context; EM is not MCMC but the latent-variable augmentation idea is shared.
• 3.1 Bayesian Framework — prior, likelihood, posterior; EM computes the MAP rather than sampling the posterior.

Intuition

For many real models, the marginal likelihood $p(y\mid \theta )={\sum }_{z}p(y,z\mid \theta )$ is a sum over a combinatorially large latent space, and maximizing it directly by gradient ascent is numerically unstable or impossible. The EM trick is to introduce the latent $z$ explicitly and alternate between (E) inferring $z$ given current $\theta$, and (M) re-estimating $\theta$ given current inferred $z$. Each step is an easy sub-problem: the E-step is a posterior calculation on $z$ with $\theta$ held fixed, and the M-step is a weighted MLE calculation on $\theta$ with $z$ summarized by its responsibilities. Repeating drives the observed-data log-likelihood uphill.

The central tension is between the tractable complete data log-likelihood $\log p(y,z\mid \theta )$ (which EM uses through its expectation) and the intractable marginal log-likelihood $\log p(y\mid \theta )=\log {\sum }_{z}p(y,z\mid \theta )$ (what we actually want to maximize). EM exploits a variational identity (Jensen’s inequality applied to the marginal) to turn the intractable log-sum into a tractable sum-of-logs that can be optimized in closed form at each M-step. The cost is that EM finds only a local optimum — the landscape of $\log p(y\mid \theta )$ typically has many modes in latent-variable models — so multiple restarts from different initializations are the practical norm.

Properties

(Monotonic Ascent) Observed-data log-likelihood never decreases across an iteration.

The EM update satisfies $\log p(y\mid {\theta }^{(t+1)})\ge \log p(y\mid {\theta }^{(t)})$ for every $t$, with equality only at stationary points. This is proved via the evidence lower bound (ELBO): for any distribution $q(z)$, $\log p(y\mid \theta )\ge {\mathbb{E}}_q[\log p(y,z\mid \theta )]-{\mathbb{E}}_q[\log q(z)]=Q(\theta \mid q)+H(q)$. Setting $q=p(z\mid y,{\theta }^{(t)})$ makes the bound tight at ${\theta }^{(t)}$; maximizing $Q$ over $\theta$ (the M-step) can only increase the bound and hence cannot decrease $\log p(y\mid \theta )$. See Derivation.

(E-step as Posterior Computation) Responsibilities are the posterior $p(z\mid y,{\theta }^{(t)})$.

For discrete $z$, the responsibility $r_{ij}^{(t)}:=\mathbb{P}[z_i=j\mid y_i,{\theta }^{(t)}]$ is computed by Bayes’ rule using the current parameters: $r_{ij}^{(t)}\propto w_j^{(t)}p(y_i\mid z_i=j,{\theta }^{(t)})$, normalized over $j$. For continuous $z$ (linear factor models, Kalman smoothers) the E-step is a conditional Gaussian computation. In every case, the E-step is a posterior inference problem on the latent variables with parameters frozen; it is the step where the latent structure of the model is exploited.

(M-step as Weighted MLE) $\arg {\max }_{\theta }Q(\theta \mid {\theta }^{(t)})$ is usually a weighted complete-data MLE.

When $\log p(y,z\mid \theta )$ is a sum over data points with separable $\theta$-terms, $Q(\theta \mid {\theta }^{(t)})$ decomposes into the same sum, now weighted by responsibilities. For a Gaussian mixture,
$Q(\theta \mid {\theta }^{(t)})={\sum }_{i=1}^n{\sum }_{j=1}^k{r}_{ij}^{(t)}[\log w_j+\log \mathcal{N}(y_i;{\mu }_j,{\sigma }_j^2)],$
and maximizing over $(w,\mu ,\sigma )$ gives the responsibility-weighted sample estimators:
$w_j^{(t+1)}=\frac{{\sum }_i{r}_{ij}^{(t)}}{n},{\mu }_j^{(t+1)}=\frac{{\sum }_i{r}_{ij}^{(t)}{y}_i}{{\sum }_i{r}_{ij}^{(t)}},({\sigma }_j^2{)}^{(t+1)}=\frac{{\sum }_i{r}_{ij}^{(t)}(y_i-{\mu }_j^{(t+1)}{)}^2}{{\sum }_i{r}_{ij}^{(t)}}.$
These are the MLE formulas with “soft counts” $r_{ij}^{(t)}$ replacing hard $\{0,1\}$ indicators.

(Relation to Gibbs Sampling) EM replaces Gibbs draws with expectations and argmaxes.

For the same latent-variable augmentation $(y,z,\theta )$, Gibbs alternates $z^{(t)}\sim p(z\mid y,{\theta }^{(t)})$ and ${\theta }^{(t+1)}\sim p(\theta \mid y,z^{(t)})$ — producing samples from the joint posterior. EM replaces $z^{(t)}$ with the full responsibility distribution $q^{(t)}(z)=p(z\mid y,{\theta }^{(t)})$ (E-step) and replaces ${\theta }^{(t+1)}$ with its argmax $\arg {\max }_{\theta }Q(\theta \mid {\theta }^{(t)})$ (M-step). Gibbs gives the full posterior distribution but is slower; EM gives the MAP point estimate much faster but loses uncertainty quantification. The two algorithms sit at opposite ends of a compute-vs-uncertainty tradeoff on the same model.

Explicit comparison for Gaussian mixtures: Gibbs draws $z_i^{(t)}\sim \text{Categorical}(r_{i1}^{(t)},\dots ,r_{ik}^{(t)})$ (hard stochastic assignment at iteration $t$), then updates $w,\mu ,\sigma$ by drawing from their conjugate posteriors given $z^{(t)}$. EM uses the responsibilities $r_{ij}^{(t)}$ as soft probabilistic assignments and sets $w,\mu ,\sigma$ to their MLEs given the weighted soft assignments. As the prior becomes uninformative (Dirichlet $a_j\to 0$, Gaussian $s^2\to \infty$, Inverse-Gamma $\alpha ,\beta \to 0$), the Gibbs update for $w,\mu ,{\sigma }^2$ shrinks toward its MLE, and EM emerges as the $\sigma \to 0$ limit of the Gibbs update.

(MAP Extension) EM computes the MAP instead of the MLE by adding a prior to the M-step.

Replace the M-step objective $Q(\theta \mid {\theta }^{(t)})$ with the posterior-aware quantity $Q(\theta \mid {\theta }^{(t)})+\log p(\theta )$. The monotonic-improvement property generalizes immediately: EM now climbs $\log p(\theta \mid y)$ rather than $\log p(y\mid \theta )$. Conjugate priors turn the M-step into a closed-form posterior-mode computation — for Gaussian mixtures with Dirichlet / Normal / Inverse-Gamma priors, the MAP M-step is
$w_j^{(t+1)}=\frac{a_j-1+{\sum }_i{r}_{ij}^{(t)}}{{\sum }_{j^{\prime }}(a_{j^{\prime }}-1)+n},{\mu }_j^{(t+1)}=\frac{m/s^2+{\sum }_i{r}_{ij}^{(t)}{y}_i/({\sigma }_j^2{)}^{(t)}}{1/s^2+{\sum }_i{r}_{ij}^{(t)}/({\sigma }_j^2{)}^{(t)}},$
which is recognizably the posterior mode (mean) of the conjugate Gibbs update rather than a pure sample.

(Convergence to Stationary Points, Not Necessarily Global) EM can get trapped in local maxima or saddles.

The monotonic-ascent guarantee is qualitative: the sequence $\log p(y\mid {\theta }^{(t)})$ is nondecreasing and bounded above, so it converges. Under regularity (Wu 1983), ${\theta }^{(t)}$ converges to a stationary point of $\log p(y\mid \theta )$, but which stationary point depends on the initialization. In mixture models, label-switching and empty clusters are classical local-optimum pathologies. Running EM from multiple random restarts and keeping the best result is the standard workaround; in critical applications, combining EM with simulated annealing or running a short Gibbs warmup before EM is also common.

(Convergence Speed) Linear near the optimum, with rate determined by the fraction of missing information.

Near a local optimum ${\theta }^{\ast }$, EM converges linearly: $\Vert {\theta }^{(t+1)}-{\theta }^{\ast }\Vert \approx \rho \Vert {\theta }^{(t)}-{\theta }^{\ast }\Vert$ with $\rho \in [0,1)$. The rate $\rho$ is the largest eigenvalue of the information ratio $I_{\text{miss}}/I_{\text{comp}}$, where $I_{\text{comp}}$ and $I_{\text{miss}}$ are the complete-data and missing-data Fisher informations at ${\theta }^{\ast }$. When $z$ carries little information beyond $y$ (responsibilities near $0$ or $1$), EM converges fast; when $z$ is highly informative about $\theta$ and the responsibilities are near $1/k$ for all components (highly ambiguous cluster membership), EM converges slowly. Direct gradient-based or Newton-based maximization of $\log p(y\mid \theta )$ can be faster near the optimum but requires evaluating the observed-data log-likelihood, which is the reason EM exists in the first place.

(Generalizations) ECM, GEM, and stochastic EM handle intractable M-steps and scale to large data.

When the M-step does not have a closed form, the Generalized EM (GEM) algorithm only requires $Q({\theta }^{(t+1)}\mid {\theta }^{(t)})\ge Q({\theta }^{(t)}\mid {\theta }^{(t)})$ (a partial ascent) rather than the exact argmax — monotonic ascent of $\log p(y\mid \theta )$ still holds. Expectation Conditional Maximization (ECM) replaces the joint M-step with coordinate-wise maximizations. Stochastic EM replaces the full responsibility computation at each iteration with a Monte Carlo estimate (useful for large $n$ or intractable posteriors), trading monotonic ascent for scalability. Each variant preserves some version of the ELBO inequality.

Examples

Gaussian mixture EM vs Gibbs on the same model

Model. $y_i\stackrel{\text{iid}}{\sim }{\sum }_{j=1}^k{w}_j\mathcal{N}({\mu }_j,{\sigma }_j^2)$ with priors $w\sim \text{Dirichlet}(a)$, ${\mu }_j\stackrel{\text{iid}}{\sim }\mathcal{N}(m,s^2)$, ${\sigma }_j^2\stackrel{\text{iid}}{\sim }\text{IG}(\alpha ,\beta )$.

Gibbs sampler (from 3.11.2 Gibbs Sampling):

1. Draw $z_i^{(t)}\sim \text{Categorical}(r_{i1}^{(t)},\dots ,r_{ik}^{(t)})$ with $r_{ij}^{(t)}\propto w_j^{(t)}\varphi (y_i;{\mu }_j^{(t)},({\sigma }_j^2{)}^{(t)})$.
2. Let $n_j^{(t)}={\sum }_{i}1\{z_i^{(t)}=j\}$, the integer component count.
3. Draw $w^{(t+1)}\sim \text{Dirichlet}(a_1+n_1^{(t)},\dots ,a_k+n_k^{(t)})$.
4. Draw ${\mu }_j^{(t+1)}\sim \mathcal{N}\left(\frac{m/s^2+{\sum }_{i:z_i^{(t)}=j}{y}_i/({\sigma }_j^2{)}^{(t)}}{1/s^2+n_j^{(t)}/({\sigma }_j^2{)}^{(t)}},\frac{1}{1/s^2+n_j^{(t)}/({\sigma }_j^2{)}^{(t)}}\right)$.
5. Draw $({\sigma }_j^2{)}^{(t+1)}\sim \text{IG}(\alpha +n_j^{(t)}/2,\beta +\frac{1}{2}{\sum }_{i:z_i^{(t)}=j}(y_i-{\mu }_j^{(t+1)}{)}^2)$.

EM algorithm (MLE variant — take flat priors):

1. E-step. Compute responsibilities $r_{ij}^{(t)}=w_j^{(t)}\varphi (y_i;{\mu }_j^{(t)},({\sigma }_j^2{)}^{(t)})/{\sum }_l{w}_l^{(t)}\varphi (y_i;{\mu }_l^{(t)},({\sigma }_l^2{)}^{(t)})$.
2. Compute effective (soft) counts $n_j^{(t)}={\sum }_i{r}_{ij}^{(t)}$ (now a real number, not an integer).
3. M-step — weights. $w_j^{(t+1)}=n_j^{(t)}/n$.
4. M-step — means. ${\mu }_j^{(t+1)}={\sum }_i{r}_{ij}^{(t)}{y}_i/n_j^{(t)}$ (responsibility-weighted sample mean).
5. M-step — variances. $({\sigma }_j^2{)}^{(t+1)}={\sum }_i{r}_{ij}^{(t)}(y_i-{\mu }_j^{(t+1)}{)}^2/n_j^{(t)}$ (responsibility-weighted sample variance).

What changes. The E-step replaces the stochastic Categorical draw of $z_i^{(t)}$ with the deterministic expectation $r_{ij}^{(t)}=\mathbb{E}[1\{z_i=j\}\mid y_i,{\theta }^{(t)}]$; the effective count $n_j^{(t)}={\sum }_i{r}_{ij}^{(t)}$ replaces the integer count; the Dirichlet / Normal / Inverse-Gamma draws of $w,\mu ,\sigma$ are replaced with their MLEs (posterior modes with uninformative priors). The algorithmic skeleton is identical; only the randomness is removed.

Missing-data regression (hidden labels)

Suppose we want to fit a linear regression $y_i=x_i^{\top }\beta +{\epsilon }_i$ with ${\epsilon }_i\sim \mathcal{N}(0,{\sigma }^2)$ but some observations are missing their response $y_i$. Introduce latent $z_i=y_i$ for missing entries; the complete-data log-likelihood is
$\log p(y,z\mid \beta ,{\sigma }^2)=-\frac{n}{2}\log (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}{\sum }_{i=1}^n({\tilde{y}}_i-x_i^{\top }\beta {)}^2,$
where ${\tilde{y}}_i=y_i$ if observed, $z_i$ if not.

E-step. Given $({\beta }^{(t)},({\sigma }^2{)}^{(t)})$, the conditional distribution of missing $z_i$ given $y_i$ observed and parameters is $\mathcal{N}(x_i^{\top }{\beta }^{(t)},({\sigma }^2{)}^{(t)})$. The Q-function replaces ${\tilde{y}}_i^2$ with $(x_i^{\top }{\beta }^{(t)}{)}^2+({\sigma }^2{)}^{(t)}$ for missing $i$.

M-step. Maximize $Q$ in $\beta ,{\sigma }^2$: this is ordinary least squares on the imputed data set where missing $y_i$ is replaced by $x_i^{\top }{\beta }^{(t)}$, with a variance correction. Iterating until convergence gives the MLE of $\beta ,{\sigma }^2$ under missing-at-random. This is the classical data augmentation example that gave EM its name.

EM for the Bernoulli mixture (coin-flipping)

Two biased coins $A,B$ with success probabilities $p_A,p_B$. We observe $n$ trials, each consisting of $m$ flips; the coin identity $z_i\in \{A,B\}$ is hidden. Let $y_i\in \{0,1,\dots ,m\}$ be the number of heads in trial $i$.

E-step. $\mathbb{P}(z_i=A\mid y_i,p_A,p_B)=\frac{\left(\genfrac{}{}{0ex}{}{m}{y_i}\right)p_A^{y_i}(1-p_A{)}^{m-y_i}}{\left(\genfrac{}{}{0ex}{}{m}{y_i}\right)p_A^{y_i}(1-p_A{)}^{m-y_i}+(\genfrac{}{}{0ex}{}{m}{y_i})p_B^{y_i}(1-p_B{)}^{m-y_i}}=\frac{p_A^{y_i}(1-p_A{)}^{m-y_i}}{p_A^{y_i}(1-p_A{)}^{m-y_i}+p_B^{y_i}(1-p_B{)}^{m-y_i}}$.

M-step. With responsibilities $r_i:=\mathbb{P}(z_i=A\mid y_i,p_A^{(t)},p_B^{(t)})$,
$p_A^{(t+1)}=\frac{{\sum }_i{r}_i{y}_i}{m{\sum }_i{r}_i},p_B^{(t+1)}=\frac{{\sum }_i(1-r_i)y_i}{m{\sum }_i(1-r_i)}.$
Iterate until the responsibilities stabilize. The MLE of $(p_A,p_B)$ is not identifiable up to swapping labels $A\leftrightarrow B$ — the two modes of the likelihood are reflections — so the particular solution EM returns depends on initialization.

Derivation

Goal. Prove the EM monotonic ascent property $\log p(y\mid {\theta }^{(t+1)})\ge \log p(y\mid {\theta }^{(t)})$, and show how the E-step and M-step together guarantee it.

Step 1 — ELBO decomposition. For any parameter $\theta$ and any distribution $q(z)$ over the latent variable $z$,
$\log p(y\mid \theta )=\underset{\text{ELBO}(\theta ,q)}{\underbrace{{\mathbb{E}}_q[\log p(y,z\mid \theta )]-{\mathbb{E}}_q[\log q(z)]}}+\underset{\ge 0}{\underbrace{\mathrm{KL}(q\parallel p(z\mid y,\theta ))}}.$
This identity follows from writing $\log p(y\mid \theta )=\log p(y,z\mid \theta )-\log p(z\mid y,\theta )$, taking expectations under $q$, and rearranging. Since $\mathrm{KL}\ge 0$,
$\log p(y\mid \theta )\ge \mathrm{ELBO}(\theta ,q)$
for every $q$, with equality when $q(z)=p(z\mid y,\theta )$.

Step 2 — E-step makes the bound tight. At iteration $t$, set $q^{(t)}(z):=p(z\mid y,{\theta }^{(t)})$. Then $\mathrm{KL}(q^{(t)}\parallel p(z\mid y,{\theta }^{(t)}))=0$, so
$\log p(y\mid {\theta }^{(t)})=\mathrm{ELBO}({\theta }^{(t)},q^{(t)})=Q({\theta }^{(t)}\mid {\theta }^{(t)})+H(q^{(t)}),$
where $H(q)=-{\mathbb{E}}_q\log q$ is the entropy of $q$. The Q-function $Q(\theta \mid {\theta }^{(t)}):={\mathbb{E}}_{q^{(t)}}[\log p(y,z\mid \theta )]$ is the part of ELBO that depends on $\theta$; the entropy $H(q^{(t)})$ does not.

Step 3 — M-step increases ELBO at $\theta$. With $q=q^{(t)}$ held fixed, maximize the ELBO over $\theta$:
${\theta }^{(t+1)}=\arg {\max }_{\theta }\mathrm{ELBO}(\theta ,q^{(t)})=\arg {\max }_{\theta }Q(\theta \mid {\theta }^{(t)}),$
since the entropy $H(q^{(t)})$ does not depend on $\theta$. By construction,
$\mathrm{ELBO}({\theta }^{(t+1)},q^{(t)})\ge \mathrm{ELBO}({\theta }^{(t)},q^{(t)})=\log p(y\mid {\theta }^{(t)}).$

Step 4 — Conclude monotonic ascent. Apply the ELBO inequality at the new parameter ${\theta }^{(t+1)}$ and any distribution $q$:
$\log p(y\mid {\theta }^{(t+1)})\ge \mathrm{ELBO}({\theta }^{(t+1)},q^{(t)})\ge \log p(y\mid {\theta }^{(t)}).$
The first inequality is ELBO $\le$ log-likelihood for any $q$; the second is the M-step improvement. Combined,
$\log p(y\mid {\theta }^{(t+1)})\ge \log p(y\mid {\theta }^{(t)}),$
which is the promised monotonic ascent. The sequence $\log p(y\mid {\theta }^{(t)})$ is nondecreasing and bounded above (by the global maximum of the log-likelihood), so it converges. Under mild regularity conditions (continuity of $Q$ in $\theta$), ${\theta }^{(t)}$ itself converges to a stationary point of $\log p(y\mid \theta )$ — local max, saddle, or local min — but not necessarily the global max.

Why this works algorithmically. The E-step replaces an intractable $\log \int p(y,z\mid \theta )dz$ with a tractable expectation ${\mathbb{E}}_q\log p(y,z\mid \theta )$, using the best possible $q$ (the posterior under ${\theta }^{(t)}$) so that nothing is thrown away. The M-step then optimizes the tractable surrogate. Because the ELBO touches $\log p(y\mid \theta )$ at ${\theta }^{(t)}$ and lies beneath it everywhere, improving the ELBO is guaranteed to improve $\log p(y\mid \theta )$.

See Also

• 3.11.2 Gibbs Sampling — stochastic version of EM on the same latent-variable augmentation; gives posterior samples rather than point estimate.
• 3.11 MCMC — parent index; EM is adjacent to rather than an instance of MCMC.
• 3.1 Bayesian Framework — prior, likelihood, posterior; EM finds MAP (or MLE) rather than sampling the full posterior.
• 3.9 Laplace Approximation — related tool: approximates the posterior around the MAP (which EM can find) using the Hessian.
• 3.7 Marginal Likelihood — the ELBO is a lower bound on the log marginal likelihood, linking EM to variational inference.