Expectation–Maximization

The expectation–maximization (EM) algorithm computes maximum likelihood estimates for models in which some of the variables are never observed: the component that generated each pixel colour in a Gaussian mixture, the hidden state sequence of a hidden Markov model, the motion layer a pixel belongs to, or simply the values of pixels that are missing. It alternates two steps. The E step uses the current parameters to compute a distribution over the hidden variables; the M step refits the parameters as if that distribution were the truth, which is usually as easy as fitting the model to fully labelled data. Each pair of steps increases the likelihood of the observed data, or leaves it unchanged at a stationary point. Dempster, Laird and Rubin (1977) presented the algorithm at a general level, with theory for the monotone behaviour of the likelihood and many examples (their published summary); special cases had been in use for several decades before (Neal & Hinton, 1998, §1).

The Gaussian Mixture Model article covers the model the examples fit, its singularities, covariance types and the choice of the number of components; the Evidence Lower Bound article derives the bound in the setting of deep generative models.

Problem

Let $X$ be the observed data and $Z$ the unobserved variables, with a joint model $p(X, Z \mid \theta)$. Maximum likelihood asks for the $\theta$ that maximizes the probability of what was observed, the incomplete-data log-likelihood

$$\ell(\theta) = \log p(X \mid \theta) = \log \sum_{Z} p(X, Z \mid \theta),$$

with the sum replaced by an integral for continuous $Z$ (Bishop, 2006, §9.3, Eq. 9.29). The sum inside the logarithm is the difficulty. Even when the joint distribution is in the exponential family, so that $\log p(X, Z \mid \theta)$ is simple and the complete-data problem (with $Z$ known) is much easier, the marginal $p(X \mid \theta)$ usually is not, and for a mixture, setting the gradient of $\ell$ to zero gives coupled equations without a closed-form solution (Bishop, 2006, §9.2.2 and §9.3). For a Gaussian mixture fitted to the pixels of an image, knowing which region every pixel came from would reduce the fit to one mean and covariance per region; not knowing it couples all components through the logarithm of a sum.

Inputs and Outputs

• Inputs: the observed data, here a data set $X \in \mathbb{R}^{n \times d}$ with one sample $\mathbf{x}_i^\top$ per row as scikit-learn stores it (pixel colours, feature vectors, or partially observed pixels); a model $p(\mathbf{x}, \mathbf{z} \mid \theta)$ whose posterior $p(\mathbf{z} \mid \mathbf{x}, \theta)$ can be computed and whose complete-data likelihood can be maximized; an initial value $\theta^{(0)}$; and a stopping rule.
• Outputs: parameters $\hat\theta$ at a stationary point of $\ell$, usually a local maximum; the posterior over the hidden variables at $\hat\theta$, which for a mixture is the matrix of responsibilities (soft labels for every pixel); and the sequence of log-likelihood values, which must never decrease.

Intuition

If the hidden labels were known, fitting would be easy; if the parameters were known, inferring the labels would be easy. EM breaks the circularity by alternating. With the current parameters, compute for every pixel how probable each label is. Then treat those probabilities as fractional labels: a pixel that is 70% “background” contributes 0.7 of a sample to the background’s mean and covariance and 0.3 to the other components. Refit, and repeat.

Why this should improve the likelihood is less obvious, and the explanation is a bound. Any guess $q$ about the hidden labels defines a function of $\theta$ that lies below the log-likelihood everywhere. The E step chooses the guess that makes the bound touch the log-likelihood at the current parameters; the M step climbs to the top of that bound. Since the bound touches at the old parameters and lies below the log-likelihood everywhere, the log-likelihood at the new parameters is at least the bound’s maximum, which is at least the old value. Murphy (2022, §8.7.1) presents EM as one case of this pattern, bound optimization, also called minorize–maximize (MM).

Algorithm

The general algorithm (Bishop, 2006, §9.3):

1. Choose initial parameters $\theta^{(0)}$; set $t = 0$.
2. E step. Compute the posterior over the hidden variables, $q_t(Z) = p(Z \mid X, \theta^{(t)})$, or the expected sufficient statistics it implies.
3. M step. Set $\theta^{(t+1)} = \arg\max_\theta Q(\theta, \theta^{(t)})$, where $Q(\theta, \theta^{(t)}) = \sum_Z q_t(Z) \log p(X, Z \mid \theta)$ is the expected complete-data log-likelihood.
4. Evaluate $\ell(\theta^{(t+1)})$. Stop if the increase in $\ell$, or the change in the parameters, is below a threshold; otherwise increase $t$ and return to step 2.

For a Gaussian mixture the E step computes the responsibilities $\gamma_{ik}$ and the M step sets each weight, mean and covariance to a responsibility-weighted average (Bishop, 2006, §9.2.2, Eqs. 9.23–9.28); the updates are derived below. Since the log-likelihood is a by-product of the E step for a mixture (it is the log of the normalizer of the responsibilities), step 4 costs nothing extra.

Mathematical Formulation

The notation follows the site’s convention: bold lowercase column vectors, $p(\cdot)$ for densities, $P(\cdot)$ for probabilities of events, and natural logarithms.

The lower bound

For any distribution $q(Z)$ over the hidden variables, the log-likelihood splits into two terms (Bishop, 2006, §9.4, Eqs. 9.70–9.72):

$$\log p(X \mid \theta) = \mathcal{L}(q, \theta) + \operatorname{KL}\big(q \,\Vert\, p(Z \mid X, \theta)\big),$$ $$\mathcal{L}(q, \theta) = \sum_{Z} q(Z) \log \frac{p(X, Z \mid \theta)}{q(Z)}, \qquad \operatorname{KL}\big(q \,\Vert\, p(Z \mid X, \theta)\big) = -\sum_{Z} q(Z) \log \frac{p(Z \mid X, \theta)}{q(Z)}.$$

The identity follows by writing $p(X, Z \mid \theta) = p(Z \mid X, \theta)\, p(X \mid \theta)$ inside $\mathcal{L}$ and using $\sum_Z q(Z) = 1$. The Kullback–Leibler divergence is non-negative and zero only when $q$ equals the posterior, so $\mathcal{L}(q, \theta)$ is a lower bound on the log-likelihood for every $q$, and the bound is tight exactly at $q = p(Z \mid X, \theta)$. In variational inference and in variational autoencoders the same $\mathcal{L}$ is called the evidence lower bound (ELBO) (Murphy, 2022, §8.7.2.1); the Evidence Lower Bound article gives both of its derivations. Neal and Hinton (1998, §2, Eqs. 2–3) wrote the same function as $F(\tilde P, \theta)$, an expected complete-data log-likelihood plus the entropy of $\tilde P$, the distribution over the unobserved variables, and noted its analogy with a negative free energy in statistical physics.

E step and M step as coordinate ascent

EM maximizes $\mathcal{L}(q, \theta)$ alternately over its two arguments (Bishop, 2006, §9.4; Neal & Hinton, 1998, §2, Eq. 5):

• E step: with $\theta^{(t)}$ fixed, $\log p(X \mid \theta^{(t)})$ does not depend on $q$, so maximizing $\mathcal{L}$ over $q$ is minimizing the KL term. The maximizer is the posterior, $q_t(Z) = p(Z \mid X, \theta^{(t)})$, and the bound then equals the log-likelihood.
• M step: with $q_t$ fixed, maximize $\mathcal{L}(q_t, \theta)$ over $\theta$. Substituting $q_t$,

$$\mathcal{L}(q_t, \theta) = \underbrace{\sum_{Z} p(Z \mid X, \theta^{(t)}) \log p(X, Z \mid \theta)}_{Q(\theta,\, \theta^{(t)})} \; \underbrace{- \sum_{Z} q_t(Z) \log q_t(Z)}_{\text{entropy of } q_t,\ \text{independent of } \theta},$$

so the M step maximizes the Q function, the expected complete-data log-likelihood (Bishop, 2006, Eq. 9.74). The logarithm now acts on the joint distribution, not on a sum. When the joint is in the exponential family, $Q$ depends on the data only through the expected sufficient statistics, and the M step is the complete-data maximum likelihood estimate with those statistics in place of the observed ones (Murphy, 2022, §8.7.2.3).

For independent samples the posterior factorizes over them, $p(Z \mid X, \theta) = \prod_i p(\mathbf{z}_i \mid \mathbf{x}_i, \theta)$, so the E step is a separate posterior computation for every pixel (Bishop, 2006, Eq. 9.75).

Monotone increase

Chaining the two steps gives the central property:

$$\ell(\theta^{(t+1)}) = \mathcal{L}(q_t, \theta^{(t+1)}) + \operatorname{KL}\big(q_t \,\Vert\, p(Z \mid X, \theta^{(t+1)})\big) \;\ge\; \mathcal{L}(q_t, \theta^{(t+1)}) \;\ge\; \mathcal{L}(q_t, \theta^{(t)}) = \ell(\theta^{(t)}).$$

The first inequality holds because the KL divergence is non-negative, the second because the M step maximizes over $\theta$, and the final equality because the E step made the bound tight. The log-likelihood therefore rises by at least as much as the bound, and by more whenever the old posterior $q_t$ differs from the new one (Bishop, 2006, §9.4, Figs. 9.11–9.14). Murphy (2022, §8.7.1) notes that if the objective does not increase monotonically there is an error in the mathematics or the code, which makes the check a debugging tool. At the point of contact the bound and the log-likelihood also have the same gradient with respect to $\theta$ (Bishop, 2006, §9.4).

Neal and Hinton (1998, §2, Theorem 2) proved the result that justifies partial steps: if $F$ has a local maximum at $(\tilde P^\ast, \theta^\ast)$, then $\ell$ has a local maximum at $\theta^\ast$, and the same holds for global maxima (assuming, for the local case, that the joint probability is positive and continuous in $\theta$). A procedure that only increases $\mathcal{L}$ in one argument at a time, or optimizes both at once, is therefore justified whenever it reaches a maximum of $\mathcal{L}$; this is the basis of the variants below.

EM for a Gaussian mixture

A mixture of $K$ Gaussians with weights $\pi_k$, means $\boldsymbol{\mu}_k$ and covariances $\Sigma_k$ introduces a hidden label $z_i \in \{1, \dots, K\}$ for every pixel, with $P(z_i = k) = \pi_k$ and $p(\mathbf{x}_i \mid z_i = k) = \mathcal{N}(\mathbf{x}_i; \boldsymbol{\mu}_k, \Sigma_k)$; the Gaussian Mixture Model article defines the model. With the labels written as indicators $z_{ik} = [z_i = k]$, the complete-data log-likelihood is (Bishop, 2006, §9.3.1, Eq. 9.36)

$$\log p(X, Z \mid \theta) = \sum_{i=1}^{n} \sum_{k=1}^{K} z_{ik} \left\{ \log \pi_k + \log \mathcal{N}(\mathbf{x}_i; \boldsymbol{\mu}_k, \Sigma_k) \right\}.$$

E step. The posterior factorizes over pixels, and the expected indicator is the responsibility (Bishop, 2006, Eqs. 9.38–9.39):

$$\gamma_{ik} = \mathbb{E}[z_{ik}] = P(z_i = k \mid \mathbf{x}_i, \theta^{(t)}) = \frac{\pi_k \, \mathcal{N}(\mathbf{x}_i; \boldsymbol{\mu}_k, \Sigma_k)}{\sum_{j=1}^{K} \pi_j \, \mathcal{N}(\mathbf{x}_i; \boldsymbol{\mu}_j, \Sigma_j)},$$

evaluated at the current parameters. Since $\log p(X, Z \mid \theta)$ is linear in the indicators, the Q function replaces each $z_{ik}$ by $\gamma_{ik}$ (Bishop, 2006, Eq. 9.40; Murphy, 2022, Eq. 8.164):

$$Q(\theta, \theta^{(t)}) = \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} \left\{ \log \pi_k - \tfrac{1}{2} \log \det \Sigma_k - \tfrac{1}{2} (\mathbf{x}_i - \boldsymbol{\mu}_k)^\top \Sigma_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k) \right\} + \text{const}.$$

M step. The function separates into one weighted single-Gaussian fit per component plus a term in the weights (Bishop, 2006, §9.3.1). Write $N_k = \sum_i \gamma_{ik}$, the effective number of pixels assigned to component $k$.

Means. The gradient with respect to $\boldsymbol{\mu}_k$ is $\sum_i \gamma_{ik} \Sigma_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k)$. Setting it to zero and multiplying by $\Sigma_k$ gives

$$\boldsymbol{\mu}_k^{\text{new}} = \frac{1}{N_k} \sum_{i=1}^{n} \gamma_{ik}\, \mathbf{x}_i .$$

Covariances. Writing the $k$-th term in the precision $\Lambda_k = \Sigma_k^{-1}$ gives $\frac{1}{2} N_k \log \det \Lambda_k - \frac{1}{2} \operatorname{tr}\big(\Lambda_k \sum_i \gamma_{ik} (\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^\top\big)$. With $\partial \log \det \Lambda / \partial \Lambda = \Lambda^{-1}$ for symmetric $\Lambda$ and $\partial \operatorname{tr}(\Lambda B) / \partial \Lambda = B^\top$, the stationarity condition is $N_k \Sigma_k = \sum_i \gamma_{ik} (\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^\top$, and with the new mean substituted (our derivation of Bishop’s Eq. 9.19):

$$\Sigma_k^{\text{new}} = \frac{1}{N_k} \sum_{i=1}^{n} \gamma_{ik}\, (\mathbf{x}_i - \boldsymbol{\mu}_k^{\text{new}})(\mathbf{x}_i - \boldsymbol{\mu}_k^{\text{new}})^\top .$$

The denominator is $N_k$, the weighted analogue of the maximum likelihood covariance with denominator $n$, not $n - 1$.

Weights. Maximize $\sum_k N_k \log \pi_k$ subject to $\sum_k \pi_k = 1$ with a Lagrange multiplier $\lambda$: $N_k / \pi_k + \lambda = 0$; multiplying by $\pi_k$ and summing over $k$ gives $\lambda = -n$, and (Bishop, 2006, §9.3.1, and the same steps on the incomplete-data likelihood in Eqs. 9.20–9.22)

$$\pi_k^{\text{new}} = \frac{N_k}{n}.$$

These are the updates of Bishop (2006, Eqs. 9.24–9.27) and Murphy (2022, Eqs. 8.165–8.167). Bishop first obtains them in §9.2.2 by setting the derivatives of the incomplete-data log-likelihood to zero: the same weighted averages appear as stationarity conditions, but there the responsibilities depend on the parameters, so the conditions are not a solution; EM turns them into an iteration by freezing the responsibilities in each M step. Because the means enter the covariance update, the new means are computed first (Bishop, 2006, §9.2.2). By the same derivation, a mixture with other component densities, such as the Bernoulli components for binary images of Bishop (2006, §9.3.3), keeps the E step as Bayes’ rule with those densities, and only the M step changes.

Missing values

The same machinery fits a model to data with entries missing, such as pixels whose red channel was lost: the missing entries are the hidden variables, and the likelihood of the observed entries is obtained by marginalizing them out (Bishop, 2006, §9.3). For a single Gaussian over RGB colours, split each pixel into observed entries $o$ and missing entries $m$. The E step needs the conditional distribution of the missing part, which for a Gaussian is again Gaussian (Bishop, 2006, §2.3.1, Eqs. 2.81–2.82):

$$\mathbb{E}[\mathbf{x}_m \mid \mathbf{x}_o] = \boldsymbol{\mu}_m + \Sigma_{mo} \Sigma_{oo}^{-1} (\mathbf{x}_o - \boldsymbol{\mu}_o), \qquad \operatorname{Cov}[\mathbf{x}_m \mid \mathbf{x}_o] = \Sigma_{mm} - \Sigma_{mo} \Sigma_{oo}^{-1} \Sigma_{om}.$$

The sufficient statistics of a Gaussian are $\sum_i \mathbf{x}_i$ and $\sum_i \mathbf{x}_i \mathbf{x}_i^\top$, so the M step needs $\mathbb{E}[\mathbf{x}_i]$ and $\mathbb{E}[\mathbf{x}_i \mathbf{x}_i^\top]$; the second is the outer product of the filled-in vector plus the conditional covariance in the missing block (our derivation, from the expected sufficient statistics of Murphy, 2022, §8.7.2.3). Dropping that covariance, which is what filling in the conditional mean and then computing an ordinary covariance does, underestimates the spread. Bishop (2006, §9.3) adds the condition for the whole approach to be valid: values must be missing at random, meaning that the mechanism that loses them does not depend on the unobserved values, which fails, for instance, for a sensor that returns nothing when its measurement exceeds a threshold. Murphy (2022, §1.5.5) distinguishes missing completely at random (MCAR, missingness independent of all values), missing at random (MAR, dependent only on observed values) and not missing at random (NMAR); in the first two cases the mechanism can be ignored. The third example below shows all three on pixel colours.

Parameters

• Initialization. EM converges to a stationary point near where it starts, so the start matters. For mixtures, a common choice is a k-means solution, with covariances set to the sample covariances of the clusters and weights to their fractions of the data, since EM needs many more iterations than k-means and each is costlier (Bishop, 2006, §9.2.2); Murphy (2022, §8.7.3.5) suggests k-means++ seeding. Random starts at data points are the cheap alternative.
• Restarts. Because of local maxima, running EM from several starts and keeping the solution with the highest likelihood is standard; Wu (1983, §3) recommended starting points representative of the parameter space. scikit-learn’s n_init keeps the run with the largest lower bound.
• Stopping rule. Stop when the gain in log-likelihood, or the change in the parameters, falls below a threshold (Bishop, 2006, §9.2.2), with a cap on the number of iterations. A small gain does not always mean convergence: EM can creep across a plateau of small gains and then speed up again (first example).
• Regularization. The likelihood of a mixture is unbounded, so implementations add a small constant to the covariance diagonals or use a prior (MAP-EM, under Variants); Maximum Likelihood Estimation and Gaussian Mixture Model explain the collapse and its remedies.

scikit-learn’s GaussianMixture implements EM with the defaults listed in Gaussian Mixture Model, among them tol=1e-3 and max_iter=100; its documentation for version 1.7.2 says that iterations stop when the average gain of the lower bound falls below tol. Three details of the 1.7.2 source matter when reading its output. The lower bound for a Gaussian mixture is the mean log-likelihood per sample, so the tolerance applies per pixel: on a 100,000-pixel image the total log-likelihood may still be rising by up to 100 per iteration when the fit stops. The value is computed in the E step, before the M step of the same iteration, so lower_bound_ and the entries of lower_bounds_ belong to the parameters one M step earlier than the ones returned; by monotonicity, score on the training data is at least lower_bound_. And the default k-means initialization runs KMeans with n_init=1, so restarts come only from n_init; after the best run is chosen, a final E step makes fit_predict(X) agree with fit(X).predict(X).

Complexity

By our count, one iteration for a $K$-component mixture with full covariances in $d$ dimensions costs $O(K d^3)$ to factor the covariances, $O(n K d^2)$ to evaluate the $nK$ Mahalanobis distances in the E step, and $O(n K d^2)$ to accumulate the weighted means and covariances in the M step, with $O(nK)$ memory for the responsibilities. For pixel colours ($d = 3$) the cost is dominated by the $nK$ terms and grows linearly in the number of pixels. The number of iterations is the larger uncertainty: in the second example it ranges from 6 to 326 depending on how much the components overlap and where EM starts. In practice the E step is computed in the log domain with the log-sum-exp trick, as in the examples, because component densities far from a pixel underflow in floating point.

Implementation

Each example renders a synthetic image with known regions and uses its pixel colours as data, one pixel per row. The code runs with NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2.

EM for a colour mixture, with the bound made visible

The first script implements the E and M steps above in NumPy, starts from a deliberately poor guess (three random pixels as means, the covariance of the whole image for every component, equal weights), and prints at each iteration the log-likelihood, the bound right after the E step, the bound after the M step with the E-step responsibilities still in place, and the new log-likelihood. The last column is the KL divergence between the old and the new posteriors. It then runs scikit-learn from the same starting values with reg_covar=0 and compares.

import os
os.environ["OMP_NUM_THREADS"] = "1"
import warnings
import numpy as np
from scipy.special import logsumexp
from sklearn.mixture import GaussianMixture

# A 64 x 64 RGB image with three regions whose colours come from known Gaussians:
# a background, a disk and a bar. One pixel per row of X, as scikit-learn expects.
rng = np.random.default_rng(0)
H = W = 64
yy, xx = np.mgrid[0:H, 0:W]
label = np.zeros((H, W), int)
label[(xx - 22) ** 2 + (yy - 24) ** 2 < 14**2] = 1
label[(np.abs(xx - 46) < 10) & (np.abs(yy - 44) < 16)] = 2
true_mu = np.array([[70.0, 90.0, 120.0], [180.0, 70.0, 60.0], [120.0, 150.0, 80.0]])
true_S = []
for m, s_shade in zip(true_mu, [10.0, 15.0, 12.0]):
    a = m / np.linalg.norm(m)                   # shading moves along the colour direction
    true_S.append(s_shade**2 * np.outer(a, a) + 5.0**2 * np.eye(3))
X = np.empty((H * W, 3))
for k in range(3):
    idx = np.flatnonzero(label.ravel() == k)
    X[idx] = rng.multivariate_normal(true_mu[k], true_S[k], size=idx.size)
n, d = X.shape


def log_joint(X, pi, mu, S):
    """log pi_k + log N(x_i; mu_k, S_k) for every pixel i and component k."""
    out = np.empty((len(X), len(pi)))
    for k in range(len(pi)):
        L = np.linalg.cholesky(S[k])
        z = np.linalg.solve(L, (X - mu[k]).T)  # whitened residuals
        out[:, k] = (np.log(pi[k]) - 0.5 * (z**2).sum(0)
                     - np.log(np.diag(L)).sum() - 0.5 * X.shape[1] * np.log(2 * np.pi))
    return out


def e_step(X, pi, mu, S):
    lj = log_joint(X, pi, mu, S)
    lp = logsumexp(lj, axis=1)                  # log p(x_i | theta)
    return np.exp(lj - lp[:, None]), lp.sum()   # responsibilities, log-likelihood


def m_step(X, R):
    Nk = R.sum(0)
    mu = (R.T @ X) / Nk[:, None]
    S = np.array([((R[:, k, None] * (X - mu[k])).T @ (X - mu[k])) / Nk[k] for k in range(len(Nk))])
    return Nk / len(X), mu, S


def bound(X, R, pi, mu, S):
    """L(q, theta) = sum_i sum_k q_ik [log p(x_i, z_i = k | theta) - log q_ik]."""
    lj = log_joint(X, pi, mu, S)
    return np.sum(R * (lj - np.log(np.clip(R, 1e-300, None))))


# A poor start: three random pixels as means, the global covariance, equal weights.
pi = np.full(3, 1 / 3)
start = rng.choice(n, 3, replace=False)
mu = X[start].copy()
S = np.array([np.cov(X.T, bias=True)] * 3)
theta0 = (pi.copy(), mu.copy(), S.copy())
print("regions of the three starting pixels (0 = background):", label.ravel()[start].tolist())

print(" it   loglik(t)   bound after E   bound after M   loglik(t+1)   KL gap")
lls, touch = [], []
for t in range(40):
    R, ll = e_step(X, pi, mu, S)                # E step: q = posterior at theta_t
    lls.append(ll)
    bE = bound(X, R, pi, mu, S)                 # equals ll: the bound touches
    touch.append(abs(bE - ll) / n)
    pi, mu, S = m_step(X, R)                    # M step: maximize the bound over theta
    bM = bound(X, R, pi, mu, S)
    _, ll_next = e_step(X, pi, mu, S)
    if t < 8 or t in (10, 13, 20, 39):
        print(f"{t:3d} {ll / n:11.6f} {bE / n:15.6f} {bM / n:15.6f} {ll_next / n:13.6f} {(ll_next - bM) / n:9.2e}")
print(f"largest |bound after E - loglik| per pixel: {max(touch):.1e}")
gain = np.diff(lls) / n
print("log-likelihood never decreased:", bool(np.all(gain >= -1e-12)),
      "| first iteration with gain below 1e-8 per pixel:", int(np.argmax(gain < 1e-8)) + 1)

# Cross-check: scikit-learn 1.7.2 from the same start, no covariance regularization.
with warnings.catch_warnings():
    warnings.simplefilter("ignore")              # tol=0 never "converges"
    gm = GaussianMixture(3, covariance_type="full", reg_covar=0.0, tol=0.0, max_iter=40,
                         weights_init=theta0[0], means_init=theta0[1],
                         precisions_init=np.linalg.inv(theta0[2])).fit(X)
# sklearn's lower_bounds_[t] is the mean log-likelihood at theta_t, computed in the E step.
print("max |sklearn lower_bounds_ - ours| per pixel:", f"{np.max(np.abs(gm.lower_bounds_ - np.array(lls) / n)):.1e}")
print("max |sklearn means_ - ours|:", f"{np.max(np.abs(gm.means_ - mu)):.1e}")
order = [int(np.argmin(np.linalg.norm(mu - m, axis=1))) for m in true_mu]
print("fitted means (rounded):", np.round(mu[order]).astype(int).tolist())
print("true weights", np.round(np.bincount(label.ravel()) / n, 3), "fitted", np.round(pi[order], 3))

Output:

regions of the three starting pixels (0 = background): [0, 0, 0]
 it   loglik(t)   bound after E   bound after M   loglik(t+1)   KL gap
  0  -14.159684      -14.159684      -13.653078    -13.621449  3.16e-02
  1  -13.621449      -13.621449      -13.517675    -13.337153  1.81e-01
  2  -13.337153      -13.337153      -12.786919    -12.561861  2.25e-01
  3  -12.561861      -12.561861      -12.361377    -12.337773  2.36e-02
  4  -12.337773      -12.337773      -12.324562    -12.312438  1.21e-02
  5  -12.312438      -12.312438      -12.299470    -12.283736  1.57e-02
  6  -12.283736      -12.283736      -12.258297    -12.217567  4.07e-02
  7  -12.217567      -12.217567      -12.044640    -11.887409  1.57e-01
 10  -11.088614      -11.088614      -11.022413    -10.990292  3.21e-02
 13  -10.828349      -10.828349      -10.760638    -10.760569  6.86e-05
 20  -10.760569      -10.760569      -10.760569    -10.760569  0.00e+00
 39  -10.760569      -10.760569      -10.760569    -10.760569  0.00e+00
largest |bound after E - loglik| per pixel: 1.8e-15
log-likelihood never decreased: True | first iteration with gain below 1e-8 per pixel: 16
max |sklearn lower_bounds_ - ours| per pixel: 4.3e-14
max |sklearn means_ - ours|: 0.0e+00
fitted means (rounded): [[70, 90, 120], [180, 70, 60], [120, 150, 80]]
true weights [0.708 0.149 0.144] fitted [0.708 0.149 0.144]

The three columns that matter tell the story of the derivation. Right after every E step the bound equals the log-likelihood to within $2 \times 10^{-15}$ per pixel: the posterior makes it tight. The M step then raises the bound, and the new log-likelihood exceeds the raised bound by the KL gap, which is positive while the parameters are moving and vanishes at convergence; the log-likelihood never decreased. The gains are far from uniform. After four gains of 0.2 to 0.8 per pixel, iterations 4 and 5 gain only about 0.03, and then the gains grow again (0.33 at iteration 7). All three starting pixels lay in the background, so in our reading the components first had to spread out within the background before two of them could move to the disk and the bar. A stopping rule with a threshold of a few hundredths per pixel would have quit on that plateau, about 1.5 per pixel below the final log-likelihood. The gain first falls below $10^{-8}$ per pixel at iteration 16, and the fit recovers the three regions’ means and their area fractions. scikit-learn, started from the same parameters, follows the same trajectory to rounding error, and its lower_bounds_ are our log-likelihoods at the start of each iteration, as the 1.7.2 source implies.

Overlap, slow convergence and local maxima

The second script renders two versions of a scene with a background, a red disk and a shadow on the background: in the first the shadow is at 60% of the background’s brightness, in the second at 90%, so that its colour population overlaps the background’s. It runs the NumPy EM from 20 random starts on each and records the number of iterations, the rate at which the late gains shrink, and the final log-likelihoods. It then compares scikit-learn’s default stopping rule with running the same k-means starts to a tolerance of $10^{-12}$.

import os
os.environ["OMP_NUM_THREADS"] = "1"
import warnings
import numpy as np
from scipy.special import logsumexp
from sklearn.mixture import GaussianMixture


def render(shadow_level, seed=0, H=48, W=48):
    """RGB image: a background, a red disk, and a shadow on the background at
    `shadow_level` times its brightness. Isotropic noise of sd 6 grey levels."""
    rng = np.random.default_rng(seed)
    yy, xx = np.mgrid[0:H, 0:W]
    label = np.zeros((H, W), int)
    label[(xx - 16) ** 2 + (yy - 16) ** 2 < 10**2] = 1
    label[(np.abs(xx - 33) < 12) & (np.abs(yy - 33) < 12) & (label == 0)] = 2
    bg = np.array([110.0, 120.0, 135.0])
    mu = np.array([bg, [190.0, 60.0, 55.0], shadow_level * bg])
    X = mu[label.ravel()] + rng.normal(0, 6.0, (H * W, 3))
    return X, label


def em(X, K, seed, tol=1e-12, max_iter=5000):
    """EM for a full-covariance GMM from a random start (means at K random pixels,
    global covariance, equal weights). Returns per-pixel log-likelihoods per iteration."""
    rng = np.random.default_rng(seed)
    n, d = X.shape
    pi, mu = np.full(K, 1 / K), X[rng.choice(n, K, replace=False)].copy()
    S = np.array([np.cov(X.T, bias=True)] * K)
    lls = []
    for _ in range(max_iter):
        lj = np.empty((n, K))
        for k in range(K):                       # E step
            L = np.linalg.cholesky(S[k])
            z = np.linalg.solve(L, (X - mu[k]).T)
            lj[:, k] = np.log(pi[k]) - 0.5 * (z**2).sum(0) - np.log(np.diag(L)).sum() - 0.5 * d * np.log(2 * np.pi)
        lp = logsumexp(lj, axis=1)
        lls.append(lp.mean())
        if len(lls) > 1 and lls[-1] - lls[-2] < tol:
            break
        R = np.exp(lj - lp[:, None])
        Nk = R.sum(0)                            # M step
        pi, mu = Nk / n, (R.T @ X) / Nk[:, None]
        S = np.array([((R[:, k, None] * (X - mu[k])).T @ (X - mu[k])) / Nk[k] + 1e-6 * np.eye(d)
                      for k in range(K)])
    return np.array(lls), mu


for name, level in [("separated (shadow at 0.60)", 0.60), ("overlapping (shadow at 0.90)", 0.90)]:
    X, label = render(level)
    runs = [em(X, 3, s) for s in range(20)]
    iters = np.array([len(ll) for ll, _ in runs])
    final = np.array([ll[-1] for ll, _ in runs])
    # Linear convergence: late gains shrink by a constant factor per iteration.
    rates = []
    for ll, _ in runs:
        g = np.diff(ll)
        ok = np.flatnonzero((g > 1e-12) & (g < 1e-6))
        if len(ok) > 2:
            rates.append(np.median(g[ok[1:]] / g[ok[1:] - 1]))
    best = final.max()
    print(f"{name}, random starts, seeds 0-19:")
    print(f"  iterations until the gain is below 1e-12 per pixel: min {iters.min()}, "
          f"median {int(np.median(iters))}, max {iters.max()}")
    print(f"  runs with a slow linear tail: {len(rates)}"
          + (f", gain ratio per iteration {np.min(rates):.2f} to {np.max(rates):.2f}" if rates else ""))
    print(f"  final log-likelihood per pixel: best {best:.4f}, worst {final.min():.4f}, "
          f"distinct optima {len(np.unique(np.round(final, 4)))}, runs at the best {np.sum(final > best - 1e-4)}")
    worst = runs[int(np.argmin(final))][1]
    print("  means of the worst run:", np.round(worst[np.argsort(worst[:, 0])]).astype(int).tolist())

# scikit-learn 1.7.2 defaults (k-means start, tol=1e-3 on the mean log-likelihood)
# against the same k-means starts run to tol=1e-12, on the overlapping image.
X, label = render(0.90)
print("region means:", np.round([X[label.ravel() == k].mean(0) for k in range(3)]).astype(int).tolist())
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    fits = [(GaussianMixture(3, random_state=s).fit(X),
             GaussianMixture(3, random_state=s, tol=1e-12, max_iter=5000).fit(X)) for s in range(10)]
it_a = [a.n_iter_ for a, b in fits]
it_b = [b.n_iter_ for a, b in fits]
dll = max(b.score(X) - a.score(X) for a, b in fits)
dmu = max(np.abs(np.sort(a.means_, 0) - np.sort(b.means_, 0)).max() for a, b in fits)
print(f"k-means starts, seeds 0-9: default tol stops after {min(it_a)} to {max(it_a)} iterations, "
      f"tol 1e-12 after {min(it_b)} to {max(it_b)}; best fit {max(b.score(X) for a, b in fits):.4f}")
best = max(b.score(X) for a, b in fits)
print(f"  runs at the best fit: {sum(b.score(X) > best - 1e-4 for a, b in fits)} of 10")
print(f"  largest log-likelihood gain from continuing {dll:.1e} per pixel; largest mean shift {dmu:.2f} grey levels")

Output:

separated (shadow at 0.60), random starts, seeds 0-19:
  iterations until the gain is below 1e-12 per pixel: min 6, median 13, max 36
  runs with a slow linear tail: 0
  final log-likelihood per pixel: best -10.5091, worst -10.5091, distinct optima 1, runs at the best 20
  means of the worst run: [[66, 72, 81], [110, 120, 135], [189, 60, 55]]
overlapping (shadow at 0.90), random starts, seeds 0-19:
  iterations until the gain is below 1e-12 per pixel: min 83, median 124, max 326
  runs with a slow linear tail: 20, gain ratio per iteration 0.77 to 0.97
  final log-likelihood per pixel: best -10.4280, worst -10.5238, distinct optima 4, runs at the best 17
  means of the worst run: [[107, 117, 132], [188, 61, 50], [190, 59, 58]]
region means: [[110, 120, 135], [189, 60, 55], [99, 108, 122]]
k-means starts, seeds 0-9: default tol stops after 3 to 3 iterations, tol 1e-12 after 72 to 74; best fit -10.4280
  runs at the best fit: 10 of 10
  largest log-likelihood gain from continuing 2.2e-04 per pixel; largest mean shift 0.30 grey levels

When the colour populations are well separated, every start ends at the same maximum within 6 to 36 iterations, and no run has a slow tail of three or more iterations with gains between $10^{-12}$ and $10^{-6}$ per pixel. By our reading, the responsibilities near the solution are almost exactly 0 or 1, so one M step nearly lands on the answer, close to the situation that Xu and Jordan’s Theorem 2 below assumes. When the shadow overlaps the background, the same algorithm needs 83 to 326 iterations, and late in every run each gain is a fixed fraction, 0.77 to 0.97 depending on the start, of the one before: convergence is linear, and the closer that fraction is to one, the slower it is. Four distinct optima appear among the 20 starts. Three runs ended below the best; the worst one split the red disk between two components and merged the background with the shadow into the third, a stationary point where EM stopped improving, not an unfinished run. All ten k-means starts reached the best fit. With that start, scikit-learn’s default tolerance stopped after 3 iterations, while running on to $10^{-12}$ took 72 to 74 iterations and gained only $2.2 \times 10^{-4}$ per pixel and 0.3 grey levels in the means. Here the default rule was harmless because the k-means start was already close to the maximum; from a poor start on overlapping data, the first example’s plateau shows how it could stop early.

Missing colour values

The third script takes the colours of a shaded surface, whose channels are strongly correlated, and deletes the red value of some pixels under three mechanisms: at random (MCAR), where the green value is high (MAR, since green is observed), and where the red value itself is high, as when a channel clips (NMAR). It estimates the red mean and standard deviation by dropping the incomplete pixels, by filling in the mean of the observed red values, and by EM for a Gaussian with the E and M steps derived above.

import os
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np

# The 4,096 pixel colours of a 64 x 64 image of one shaded surface, drawn from a
# known Gaussian elongated along the colour direction: channels rise and fall together.
rng = np.random.default_rng(3)
H = W = 64
mu_true = np.array([150.0, 100.0, 70.0])
a = mu_true / np.linalg.norm(mu_true)
S_true = 25.0**2 * np.outer(a, a) + 5.0**2 * np.eye(3)
X = rng.multivariate_normal(mu_true, S_true, size=H * W)    # one pixel per row
full_mu, full_var = X.mean(0), X.var(0)                      # what complete data would give


def em_missing(X, miss, iters=200):
    """ML estimate of a Gaussian's mean and covariance when X[i, j] is missing
    where miss[i, j] is True. E step: conditional mean and covariance of the
    missing entries given the observed ones; M step: moments of the completed data."""
    n, d = X.shape
    mu = np.nanmean(np.where(miss, np.nan, X), axis=0)
    S = np.diag(np.nanvar(np.where(miss, np.nan, X), axis=0))
    for _ in range(iters):
        Xh, C = np.where(miss, 0.0, X), np.zeros((d, d))
        for pattern in np.unique(miss, axis=0):              # rows sharing a pattern
            rows = np.all(miss == pattern, axis=1)
            m, o = pattern, ~pattern
            if not m.any():
                continue
            K = S[np.ix_(m, o)] @ np.linalg.inv(S[np.ix_(o, o)])
            Xh[np.ix_(rows, m)] = mu[m] + (X[np.ix_(rows, o)] - mu[o]) @ K.T
            C[np.ix_(m, m)] += rows.sum() * (S[np.ix_(m, m)] - K @ S[np.ix_(o, m)])
        mu = Xh.mean(0)
        S = (Xh - mu).T @ (Xh - mu) / n + C / n              # E[x x^T] needs the conditional covariance
    return mu, S


red, green = X[:, 0], X[:, 1]
mechanisms = {
    "MCAR: red lost at random (40%)": rng.random(len(X)) < 0.4,
    "MAR: red lost where green > 105": green > 105,
    "NMAR: red lost where red > 160 (clipping)": red > 160,
}
print(f"complete data: red mean {full_mu[0]:.1f}, red sd {np.sqrt(full_var[0]):.1f}")
for name, lost in mechanisms.items():
    miss = np.zeros_like(X, bool)
    miss[:, 0] = lost
    kept = ~lost
    cc = (red[kept].mean(), red[kept].std())                       # complete-case analysis
    imp = np.where(lost, red[kept].mean(), red)                    # mean imputation
    mu_em, S_em = em_missing(X, miss)
    print(f"{name}  [{lost.mean():.0%} missing]")
    print(f"  red mean / sd: drop pixels {cc[0]:6.1f} / {cc[1]:4.1f};  mean-fill {imp.mean():6.1f} / {imp.std():4.1f};"
          f"  EM {mu_em[0]:6.1f} / {np.sqrt(S_em[0, 0]):4.1f}")

Output:

complete data: red mean 149.9, red sd 20.2
MCAR: red lost at random (40%)  [40% missing]
  red mean / sd: drop pixels  149.6 / 20.1;  mean-fill  149.6 / 15.5;  EM  149.8 / 20.2
MAR: red lost where green > 105  [36% missing]
  red mean / sd: drop pixels  139.1 / 15.2;  mean-fill  139.1 / 12.1;  EM  150.0 / 20.4
NMAR: red lost where red > 160 (clipping)  [32% missing]
  red mean / sd: drop pixels  139.3 / 14.2;  mean-fill  139.3 / 11.7;  EM  147.2 / 17.8

With values lost at random, every method gets the mean right, but filling in a constant shrinks the standard deviation from 20.2 to 15.5: a quarter of the spread disappears because 40% of the pixels were given the same value. When the loss depends on green, the observed pixels are a biased sample: dropping or mean-filling puts the red mean 11 grey levels too low. EM recovers it within 0.1 and the spread within 0.2, because the E step predicts each missing red value from that pixel’s observed green and blue through the fitted covariance, and the conditional covariance term restores the variance. When the red channel is lost precisely where it is bright, the missingness is informative about the missing value, and EM is biased too, though less than the alternatives (147.2 against 149.9) because the correlated channels still carry part of the information. In our view, the model of the data cannot repair a wrong model of why data are missing: a clipped sensor calls for a model of the clipping, not EM under the missing-at-random assumption.

Properties and Behavior

Convergence

The monotone increase of the likelihood is not by itself a proof that the parameters converge, or that they converge to a maximum. Wu (1983, §1) pointed out that Dempster, Laird and Rubin’s proof that EM sequences converge contains an error, an incorrect use of the triangle inequality in their Theorem 2, while their results on the monotone likelihood and on the rate of convergence remain valid. Under standing assumptions that include compact sets $\{\theta : \ell(\theta) \ge \ell(\theta^{(0)})\}$ and a likelihood continuous on the parameter space and differentiable inside it, Wu proved (§2):

• if $Q(\theta, \theta')$ is continuous in both arguments, as it is for a curved exponential family with a compact parameter space, every limit point of an EM sequence is a stationary point of the likelihood, and the likelihood values converge monotonically to its value at a stationary point (Theorem 2);
• limit points are local maxima under a further condition that excludes EM stalling at a stationary point that is not a maximum, which he notes is hard to verify (Theorem 3);
• the parameters themselves converge to a point under additional conditions, for example when no two stationary points share the same likelihood value (Theorem 4), or when the likelihood is unimodal with a single stationary point and a differentiability condition holds (Corollary 1).

Stationary points need not be maxima. In an example of Murray with missing values from a bivariate normal distribution, an EM sequence started at zero correlation keeps zero correlation and converges to a saddle point of the likelihood (Wu, 1983, §2.1). Neal and Hinton (1998, §2) state that for most models EM converges to a local maximum, with exceptions, and set aside convergence to a saddle point. Even where it does converge to a maximum, EM can spend many iterations on a plateau of small gains before reaching it, as in the first example.

Local maxima and initialization

The likelihood of a mixture generally has several local maxima, and EM is not guaranteed to find the largest (Bishop, 2006, §9.2.2). Besides the $K!$ equivalent relabellings, which are equally good, there can be an exponential number of modes, and finding a global optimum is NP-hard according to the results Murphy cites, so in practice one settles for a good local optimum and seeds the search well, for example with k-means++ (Murphy, 2022, §8.7.3.5). In the second example, three of 20 random starts on the overlapping image ended at worse fits, one of which split one region and merged two others. Restarts and a careful initialization are the standard remedies; the likelihoods of the runs are directly comparable, so keeping the best one is well defined.

Rate of convergence

EM is a first-order method. Writing one iteration as a map $\theta^{(t+1)} = M(\theta^{(t)})$ with fixed point $\theta^\ast$, near the fixed point $\theta^{(t+1)} - \theta^\ast \approx \frac{\partial M}{\partial \theta}(\theta^\ast)\,(\theta^{(t)} - \theta^\ast)$, so the error shrinks by a roughly constant factor per iteration; Xu and Jordan (1995, §1) give this as one reason for caution about EM, and recall that Redner and Walker (1984) reported slow convergence on a Gaussian mixture whose components were not well separated and recommended superlinear and second-order methods instead. Murphy (2022, §8.7.2) states that the speed of convergence depends on the amount of missing data, which affects how tight the bound is. The second example shows both regimes on images. Overlapping colour populations leave much uncertainty about each pixel’s label, and the late gains shrink only by a factor between 0.77 and 0.97 per iteration; separated populations leave almost none, and the iteration finishes in a few steps. In our reading, the uncertain labels are the missing information: the more of it, the less each E step changes the bound and the smaller each M step’s move.

Xu and Jordan (1995) gave the other half of the picture for Gaussian mixtures. The EM step equals the gradient of the log-likelihood multiplied by a matrix $P(\theta)$ that they give explicitly and that is positive definite (for the means and covariances, with probability one when the sample is large enough) (their Theorem 1), so EM is a variable-metric gradient ascent with no step size to tune, and its search direction always has a positive projection on the gradient (§3). In their experiments $P$ reduced the condition number of the effective Hessian by large factors, and for the means of well-separated components, with the weights and covariances known and the responsibilities nearly 0 or 1, they proved that the effective Hessian’s condition number is smaller than for gradient ascent and that its largest eigenvalue tends to one, as for Newton’s method (their Theorem 2). They argued that EM’s slowness for mixtures had been overstated: where it is slow, the Hessian is poorly conditioned and other gradient methods also struggle, and convergence in likelihood is usually fast even when the parameters converge slowly (§6). Mean shift with a Gaussian kernel is itself an EM algorithm, and Carreira-Perpiñán (2007) used this to show that its convergence is in general linear, with a rate often close to one at intermediate bandwidths and sublinear at bandwidths where modes merge; as for EM in general, this does not by itself prove that the iterates converge to a point. Mean Shift covers both questions.

What EM does well

Xu and Jordan (1995, §6) list EM’s advantages for mixtures as the automatic satisfaction of the probabilistic constraints, monotone convergence without a learning rate, and low computational overhead. The constraints are satisfied because the updated weights are average responsibilities and the updated covariances weighted averages of outer products, hence positive semidefinite (our reasoning). Each step is usually cheap and in closed form, and the responsibilities it produces are themselves useful output: soft segmentations, ownership probabilities, posterior probabilities of hidden states.

Limitations

• Local, not global. The answer depends on the start; restarts multiply the cost (second example).
• Slow when the hidden variables are uncertain. Overlapping components, much missing data or weakly identified parameters give linear convergence with a rate close to one, and a stopping rule on the likelihood gain can be fooled by a plateau.
• Needs tractable steps. The E step requires the posterior over the hidden variables, or at least its expected sufficient statistics, and the M step a maximization of the expected complete-data log-likelihood. When either is intractable, EM must be replaced by one of its approximate variants.
• Inherits the likelihood’s problems. EM climbs the likelihood it is given, including the unbounded likelihood of a mixture, so it needs the regularization described under Parameters (Bishop, 2006, §9.2.2).
• Assumes the missingness can be ignored. For missing data EM is valid when values are missing at random; when the reason a value is missing depends on the value itself, as for clipped pixels, the estimates are biased (Bishop, 2006, §9.3; third example).
• Gives a point estimate. EM returns a maximum likelihood or MAP estimate, not a posterior distribution over the parameters; uncertainty about the parameters needs other methods.

Variants

• Generalized EM (GEM). The M step only increases the bound instead of maximizing it, for example with a nonlinear optimizer such as conjugate gradients; each cycle still increases the likelihood unless it is already at a local maximum (Bishop, 2006, §9.4). The name comes from Dempster, Laird and Rubin (Neal & Hinton, 1998, §1; Wu, 1983, §1).
• Expectation conditional maximization (ECM). A form of GEM in which the M step is split into several constrained maximizations, each over one group of parameters with the others fixed, due to Meng and Rubin (Bishop, 2006, §9.4; Neal & Hinton, 1998, §1).
• Incremental and online EM. Neal and Hinton (1998, §3) used their bound to justify an E step that updates the posterior of one data item at a time, followed by an M step from incrementally updated sufficient statistics, so both steps take constant time; Bishop (2006, §9.4, Eqs. 9.78–9.79) gives the incremental update of a Gaussian mixture’s means. On 1,000 points from a two-component one-dimensional mixture, the incremental version reached any given likelihood in about half as many passes through the data as standard EM, but each pass took about twice as long; recomputing ten points per E step kept almost the same rate with a pass only about 10% slower than standard EM (§4). A variant due to Nowlan keeps exponentially decaying statistics instead of exact ones; with a decay factor of 0.99 it reached a good but not optimal point faster than the incremental algorithm, and with 0.95 a rather poor one (§3–4).
• Sparse EM. For many components, most of which have negligible responsibility for a given pixel, freeze those small responsibilities and update only the plausible ones at most iterations, with an occasional full E step; the bound still increases at every iteration, although the likelihood may fall at some of them (Neal & Hinton, 1998, §5).
• Hard (winner-take-all) EM. Restrict $q$ to put all its probability on one value of each hidden variable, the most probable one. This maximizes the bound over a restricted family, so it cannot in general find the maximum of the likelihood, and the likelihood itself is not guaranteed to increase at each iteration; it may still save computation in the early iterations, before switching to a variant that can find the true maximum (Neal & Hinton, 1998, §6). k-means is the incremental winner-take-all version of EM for a Gaussian mixture with fixed variances and weights (§6), and Bishop (2006, §9.3.2) derives it as the limit of EM as a shared spherical variance goes to zero; a hard-assignment mixture with general covariances has been called elliptical k-means (Sung and Poggio, as cited by Bishop).
• Stochastic and Monte Carlo EM. When the expectation in the E step has no closed form, approximate $Q$ by an average over samples of the hidden variables drawn from their posterior (Monte Carlo EM); stochastic EM draws a single sample per E step, which for a mixture makes a random hard assignment of every point (Bishop, 2006, §11.1.6).
• MAP-EM. With a prior $p(\theta)$, the E step is unchanged and the M step maximizes $Q(\theta, \theta^{(t)}) + \log p(\theta)$ (Bishop, 2006, §9.3 and Eq. 9.77). A suitable prior removes the singularities of a mixture’s likelihood (Bishop, 2006, §9.3); Gaussian Mixture Model covers that case.
• Variational EM. When the exact posterior cannot be computed, use a tractable approximation to it as $q$; the bound is then no longer tight after the E step, and the method is called variational EM (Murphy, 2022, §8.7.2.2, citing Neal & Hinton). Bishop (2006, §9.4) presents the general EM derivation as the basis of variational inference.
• Amortized inference. A variational autoencoder maximizes the same bound for a latent-variable model whose decoder is a neural network, but instead of computing a separate $q$ for each image by an E step it trains an inference network that maps any image to its approximate posterior, and updates both networks together by gradient ascent (Murphy, 2022, §20.3.5). Murphy notes that for a fixed inference network, increasing the ELBO should increase the log-likelihood, similar to EM.

Where It Is Used

• Mixture models. Fitting Gaussian mixtures to pixel colours, local descriptors and feature vectors, as in colour models for segmentation and the vocabulary of Fisher vectors; see Gaussian Mixture Model.
• Hidden Markov models. Maximum likelihood for an HMM uses EM, with the E step computed by the forward–backward algorithm, also called the Baum–Welch algorithm, because the joint distribution of an HMM, unlike that of a mixture, does not factorize over time steps (Bishop, 2006, §13.2.1–13.2.2). Xu and Jordan (1995, §3) relate the mixing-proportion update to Baum and Sell’s growth transformations of 1968.
• Missing and incomplete data. Missing values, grouped, censored or truncated data, variance components and factor analysis are among the examples listed in Dempster, Laird and Rubin’s summary; in images the missing values may be masked, dropped or occluded pixels, under the missing-at-random condition.
• Background subtraction and interactive segmentation. Stauffer and Grimson’s per-pixel background mixtures replace exact EM by an on-line approximation (Stauffer & Grimson, 1999, §2.1), and GrabCut re-estimates its colour mixtures from hard assignments of pixels to components between graph cuts (Rother, Kolmogorov & Blake, 2004, §3.1); in our reading these are an online and a hard-EM scheme in the sense of the variants above. Gaussian Mixture Model describes both.
• Mixtures of motions. Within an image patch that straddles an occlusion boundary or contains transparency, a single optical flow vector is wrong. Jepson and Black (1993) modelled the motion constraints in a patch as a mixture of a few smooth motion layers plus an outlier process and fitted it by a simple modification of EM; the responsibilities, which they call ownership probabilities, assign each constraint to a motion, and in a $32 \times 32$ patch on an occlusion boundary the method recovered the velocities on both sides. They compare the approach with robust regression, which also estimates the dominant motion robustly; the mixture in addition recovers several coherent motions when they are present and rejects outliers.