Maximum Likelihood from Incomplete Data via the EM Algorithm

“Maximum Likelihood from Incomplete Data via the EM Algorithm” is the paper by Arthur Dempster, Nan Laird and Donald Rubin of Harvard University and the Educational Testing Service that was read before the Royal Statistical Society on 8 December 1976 and published, with the discussion it provoked, in the Society’s Journal, Series B, in 1977. It took an iteration that statisticians had rediscovered for decades in special problems, gave it one general definition and a name, the EM algorithm for its expectation and maximization steps, and proved what every such iteration has in common: the likelihood never decreases, and near a solution the error shrinks by a factor set by how much information the missing data would have carried. Its second half is a survey showing that missing values, grouped and censored data, finite mixtures, variance components, empirical Bayes, robust regression and factor analysis are all instances. The method as it is understood and used today, including its derivation as coordinate ascent on a lower bound and the mixture updates used for images, is the subject of Expectation–Maximization; this article is about the paper.

Section, equation, theorem and table numbers below are those of the paper as printed in the journal (pages 1–22), read in the JSTOR scan of the paper with its discussion; the discussion and the authors’ reply (pages 22–38) have their own DOI and are cited separately. The paper writes parameters and sufficient statistics as row vectors; here they are columns.

Problem

The paper’s examples show how common incomplete data were in the statistics of the time: plots lost from a designed experiment, variables observed for some individuals and not others, data censored at known limits or reported only in groups. For many of these problems the maximum likelihood estimate had no closed form, while the same model fitted to complete data was easy. The paper sets the problem up abstractly (§1). The complete data $\mathbf{x}$ live in a sample space $\mathcal{X}$ and are seen only through a many-to-one map $\mathbf{y} = \mathbf{y}(\mathbf{x})$; a family of complete-data densities $f(\mathbf{x} \mid \boldsymbol{\phi})$ induces the density of what is observed,

$$g(\mathbf{y} \mid \boldsymbol{\phi}) = \int_{\mathcal{X}(\mathbf{y})} f(\mathbf{x} \mid \boldsymbol{\phi}) \, d\mathbf{x}, \qquad (1.1)$$

and the task is to maximize $g$ over $\boldsymbol{\phi}$. The authors stress that the complete data are a modelling choice: one $g$ is generated by many $f$, and $\mathbf{x}$ may include quantities usually called parameters.

Earlier special cases, as the paper presents them

The paper is candid that the algorithm itself was not new (§1). It credits Hartley (1958) with three multinomial examples similar to its own, and in §4.1.1 states that Hartley introduced the procedure explicitly for samples in which some observations are assigned only to groups of cells. It finds parts of the theory in Orchard and Woodbury (1972), who called the idea a “missing information principle”, in Sundberg (1974, 1976), and, “buried” in a special case, in Baum, Petrie, Soules and Weiss (1970), whose argument already contained the monotonicity proof. The key identity for exponential families, that the score is a conditional minus an unconditional expectation of the sufficient statistics, was given in general form by Sundberg, who attributed it to unpublished 1966 lecture notes of Martin-Löf (§2).

The survey in §4 adds many more. For missing values it names Healy and Westmacott (1956) for filling in missing plots of an analysis of variance, and Orchard and Woodbury and Beale and Little (1975) for the multivariate normal. For grouped, censored and truncated data it names Grundy (1952), Blight (1970), whom it calls a key paper for recognizing the two-step interpretation in exponential families, Efron (1967) and Turnbull (1974, 1976), and notes that its fuller treatment of truncation was used for truncated zero counts by Hartley (1958) and Irwin (1959, 1963), Irwin ascribing the idea to McKendrick (1926). For mixtures it names Hasselblad (1966, 1969), Day (1969) and Wolfe (1970), who, it says, apparently worked independently apart from Wolfe’s reference to Hasselblad (1966), and all suggested EM, and Ceppellini and co-authors (1955) in genetics. It also names the mixtures of Baum and Eagon (1967), Baum et al. (1970) and Baum (1972), whose hidden states follow a Markov chain; Hartley and Rao (1967) for mixed models; and Good (1956), who used the iteration in a Bayesian setting without connecting it to maximum likelihood. Hasselblad had reported that the likelihood always increased in practice, but, the paper notes, none of these authors proved it.

Contribution

• One general definition. The E and M steps, defined first for regular exponential families, then for curved ones, then for any complete-data model through the function $Q$ (§2), and the name “EM”.
• Generalized EM. A class of iterations that only increase $Q$ instead of maximizing it, with the same guarantee (§3).
• Monotone likelihood. A short proof that every EM or GEM step leaves the likelihood unchanged or increases it (Lemma 1, Theorem 1).
• Convergence and rate. Conditions under which the iterates converge (Theorems 2 and 3), the conditions satisfied at a limit, and the rate near it as the fraction of missing information (Theorem 4).
• Posterior modes. The same algorithm and theory with a log-prior added in the M step (§2, §3).
• A survey of applications. Seven families of examples (§4), with two that the authors present as new: a general complete-data model for truncated samples (§4.2) and EM for maximum likelihood factor analysis (§4.7).

Method

The genetics example and Table 1

The paper introduces the algorithm with data from Rao’s textbook: 197 animals distributed over four categories as $\mathbf{y} = (125, 18, 20, 34)$, with cell probabilities $\big(\tfrac12 + \tfrac14\pi, \tfrac14(1-\pi), \tfrac14(1-\pi), \tfrac14\pi\big)$ under a genetic model with $0 \le \pi \le 1$ (Eq. 1.2). Splitting the first cell into two with probabilities $\tfrac12$ and $\tfrac14\pi$ makes the complete data a five-cell multinomial $\mathbf{x} = (x_1, \dots, x_5)$ in which $\pi$ has a closed-form estimate (Eq. 1.3). The E step divides the 125 animals of the first cell in proportion to the two sub-cell probabilities at the current $\pi^{(p)}$, and the M step is the binomial-type estimate from the filled-in counts:

$$x_2^{(p)} = 125 \, \frac{\tfrac14 \pi^{(p)}}{\tfrac12 + \tfrac14 \pi^{(p)}}, \qquad \pi^{(p+1)} = \frac{x_2^{(p)} + 34}{x_2^{(p)} + 34 + 18 + 20}. \qquad (1.4, 1.5)$$

Table 1 runs this from $\pi^{(0)} = 0.5$ for eight steps to $\pi = 0.626821484$, against the exact maximum $\pi^* = (15 + \sqrt{53809})/394 \approx 0.6268214980$ obtained by solving the fixed-point equation as a quadratic. Its last column, the ratio of successive deviations from $\pi^*$, settles at 0.1328 from $p = 3$, and §3 later computes the theoretical rate, $DM(\pi^*) \approx 0.132778$, which the table confirms. In our recomputation from Eqs. (1.4) and (1.5), the iterates rise towards $\pi^*$ from below, so the column headed $\pi^{(p)} - \pi^*$ prints the magnitudes of negative deviations; the iterates, deviations and ratios agree with the table to the last printed digit except for one-unit differences in the ninth decimal of the iterate and of the deviation at $p = 2$ and $p = 7$.

Three levels of definition

Section 2 defines the algorithm three times, each more general and each, the authors say, less simple to use. For a regular exponential family

$$f(\mathbf{x} \mid \boldsymbol{\phi}) = b(\mathbf{x}) \exp\big(\boldsymbol{\phi}^\top \mathbf{t}(\mathbf{x})\big) / a(\boldsymbol{\phi}), \qquad (2.1)$$

with natural parameters $\boldsymbol{\phi}$ in a convex set $\Omega$ and complete-data sufficient statistics $\mathbf{t}(\mathbf{x})$, the steps are:

• E step: estimate the sufficient statistics, $\mathbf{t}^{(p)} = \mathbb{E}\big[\mathbf{t}(\mathbf{x}) \mid \mathbf{y}, \boldsymbol{\phi}^{(p)}\big]$ (2.2);
• M step: solve the complete-data likelihood equations $\mathbb{E}\big[\mathbf{t}(\mathbf{x}) \mid \boldsymbol{\phi}\big] = \mathbf{t}^{(p)}$ for $\boldsymbol{\phi}^{(p+1)}$ (2.3).

The filled-in statistics need not correspond to any possible $\mathbf{x}$; in the genetics example they are fractional counts. Why the fixed point is the maximum likelihood estimate follows from writing the conditional density of the complete data, $k(\mathbf{x} \mid \mathbf{y}, \boldsymbol{\phi}) = f / g$ (2.5), as another exponential family with the same natural parameters: differentiating its normalizer gives

$$DL(\boldsymbol{\phi}) = -\mathbb{E}(\mathbf{t} \mid \boldsymbol{\phi}) + \mathbb{E}(\mathbf{t} \mid \mathbf{y}, \boldsymbol{\phi}), \qquad (2.13)$$

where $L(\boldsymbol{\phi}) = \log g(\mathbf{y} \mid \boldsymbol{\phi})$ (2.4) and $D$ denotes differentiation, so a fixed point of (2.2)–(2.3) is a stationary point of $L$. The paper calls this a “striking representation” noticed in special cases by many authors, and records Sundberg’s extension to all orders: derivatives of the log-likelihood are differences of conditional and unconditional cumulants, and the Hessian is a difference of two covariance matrices.

For a curved exponential family, with $\boldsymbol{\phi}$ restricted to a curved subset $\Omega_0$, the M step becomes maximization of $-\log a(\boldsymbol{\phi}) + \boldsymbol{\phi}^\top \mathbf{t}^{(p)}$ over $\Omega_0$; the same form is used for regular families whose likelihood equations have no solution inside $\Omega$. The final level drops exponential families entirely. With

$$Q(\boldsymbol{\phi}' \mid \boldsymbol{\phi}) = \mathbb{E}\big[\log f(\mathbf{x} \mid \boldsymbol{\phi}') \mid \mathbf{y}, \boldsymbol{\phi}\big], \qquad (2.17)$$

the E step computes $Q(\boldsymbol{\phi} \mid \boldsymbol{\phi}^{(p)})$ as a function of $\boldsymbol{\phi}$, and the M step chooses $\boldsymbol{\phi}^{(p+1)}$ to maximize it. The heuristic the paper gives is to maximize the complete-data log-likelihood’s current expectation, since the log-likelihood itself is unknown. A footnote records a referee’s objection that “algorithm” is a misnomer when the computing steps inside each E and M step are left unspecified; the authors concede that implementations vary widely in complexity and feasibility. For a posterior mode, the M step maximizes $Q(\boldsymbol{\phi} \mid \boldsymbol{\phi}^{(p)}) + G(\boldsymbol{\phi})$, where $G$ is the log of the prior density, and with a conjugate prior this often has the same form as $Q$ alone.

Generalized EM and monotonicity

Section 3 writes an iterative algorithm as a map, $\boldsymbol{\phi}^{(p+1)} = M(\boldsymbol{\phi}^{(p)})$ (3.4), and defines a generalized EM (GEM) algorithm as any such map with

$$Q\big(M(\boldsymbol{\phi}) \mid \boldsymbol{\phi}\big) \ge Q(\boldsymbol{\phi} \mid \boldsymbol{\phi}) \quad \text{for every } \boldsymbol{\phi} \in \Omega, \qquad (3.5)$$

of which EM, which maximizes $Q$, is a special case (3.6). With $H(\boldsymbol{\phi}' \mid \boldsymbol{\phi}) = \mathbb{E}\big[\log k(\mathbf{x} \mid \mathbf{y}, \boldsymbol{\phi}') \mid \mathbf{y}, \boldsymbol{\phi}\big]$ (3.1), $Q$ splits as $Q = L + H$ (3.2), and Lemma 1, a consequence of Jensen’s inequality, states $H(\boldsymbol{\phi}' \mid \boldsymbol{\phi}) \le H(\boldsymbol{\phi} \mid \boldsymbol{\phi})$. Theorem 1 follows in one line:

$$L(M(\boldsymbol{\phi})) - L(\boldsymbol{\phi}) = \big\{Q(M(\boldsymbol{\phi}) \mid \boldsymbol{\phi}) - Q(\boldsymbol{\phi} \mid \boldsymbol{\phi})\big\} + \big\{H(\boldsymbol{\phi} \mid \boldsymbol{\phi}) - H(M(\boldsymbol{\phi}) \mid \boldsymbol{\phi})\big\} \ge 0, \qquad (3.10)$$

with equality only if both braces vanish. The paper reads its Corollaries 1 and 2 as showing that a maximum likelihood estimate is a fixed point of a GEM algorithm (Corollary 2, for a unique maximum, gives $M(\boldsymbol{\phi}^*) = \boldsymbol{\phi}^*$). The paper is explicit that these results “stop short” of convergence to a maximum likelihood estimate.

Convergence theorems

Theorem 2 states that a GEM sequence converges to some $\boldsymbol{\phi}^*$ in the closure of $\Omega$ if (1) the sequence $L(\boldsymbol{\phi}^{(p)})$ is bounded and (2) each step increases $Q$ by at least a multiple of its squared length,

$$Q(\boldsymbol{\phi}^{(p+1)} \mid \boldsymbol{\phi}^{(p)}) - Q(\boldsymbol{\phi}^{(p)} \mid \boldsymbol{\phi}^{(p)}) \ge \lambda \, \lVert \boldsymbol{\phi}^{(p+1)} - \boldsymbol{\phi}^{(p)} \rVert^2 \quad \text{for some } \lambda > 0. \qquad (2)$$

The proof uses the convergence of the bounded, non-decreasing likelihood to make the total gain from step $p$ onwards smaller than any $\varepsilon$, and then condition (2) to bound the sum of the squared step lengths, $\varepsilon > \lambda \sum_{j=1}^{r} \lVert \boldsymbol{\phi}^{(p+j)} - \boldsymbol{\phi}^{(p+j-1)} \rVert^2$ (3.13), from which it concludes $\varepsilon > \lambda \lVert \boldsymbol{\phi}^{(p+r)} - \boldsymbol{\phi}^{(p)} \rVert^2$ (3.14), a Cauchy condition. That last step is the error Wu found (see below).

After introducing differentiability, Lemma 2 gives the familiar score identities for the conditional density, (3.15) and (3.16). Theorem 3 shows that when each step solves $D^{10} Q(\boldsymbol{\phi}^{(p+1)} \mid \boldsymbol{\phi}^{(p)}) = \mathbf{0}$ (the superscripts count derivatives in the first and second arguments), a Taylor expansion (3.17) turns negative definiteness of $D^{20} Q$, with eigenvalues bounded away from zero along the path, into condition (2) of Theorem 2, so with a bounded likelihood the sequence converges.

The rate of convergence

Theorem 4 assumes that the sequence converges, that each step solves $D^{10} Q(\boldsymbol{\phi}^{(p+1)} \mid \boldsymbol{\phi}^{(p)}) = \mathbf{0}$, and that $D^{20} Q(\boldsymbol{\phi}^{(p+1)} \mid \boldsymbol{\phi}^{(p)})$ is negative definite with eigenvalues bounded away from zero. Then the limit is a stationary point, $DL(\boldsymbol{\phi}^*) = \mathbf{0}$ (3.18), and the derivative of the EM map there is

$$DM(\boldsymbol{\phi}^*) = D^{20} H(\boldsymbol{\phi}^* \mid \boldsymbol{\phi}^*) \, \big[D^{20} Q(\boldsymbol{\phi}^* \mid \boldsymbol{\phi}^*)\big]^{-1}, \qquad (3.19)$$

which for a regular exponential family in its natural parameters becomes a ratio of covariance matrices of the complete-data sufficient statistics,

$$DM(\boldsymbol{\phi}^*) = V(\mathbf{t} \mid \mathbf{y}, \boldsymbol{\phi}^*) \, V(\mathbf{t} \mid \boldsymbol{\phi}^*)^{-1}. \qquad (3.26)$$

(These are in the paper’s row-vector convention; with column vectors the Jacobian is the transpose, with the same eigenvalues.) The paper’s interpretation (p. 10): the rate depends on the relative sizes of $D^{2}L(\boldsymbol{\phi}^*)$ and $D^{20}H(\boldsymbol{\phi}^* \mid \boldsymbol{\phi}^*)$, where $-D^{2}L$ measures the information in the observed data $\mathbf{y}$ about $\boldsymbol{\phi}$ and $-D^{20}H$ is the expected information in the unobserved part of $\mathbf{x}$; when all eigenvalues of $DM(\boldsymbol{\phi}^*)$ are below one, the largest gives the rate. If little information is lost to incompleteness, EM converges quickly; and because the loss can differ across components of $\boldsymbol{\phi}$, some components may settle fast while others need many iterations. Since $D^{20}L(\boldsymbol{\phi}^*) = \big(I - DM(\boldsymbol{\phi}^*)\big) D^{20}Q(\boldsymbol{\phi}^* \mid \boldsymbol{\phi}^*)$, an eigenvalue above one signals a saddle point, from which a small perturbation would drive EM away, and an eigenvalue equal to one throughout a neighbourhood signals a ridge of the likelihood, as in factor analysis, where the incomplete data determine the factors only up to an orthogonal transformation. The paper adds that Theorem 2 shows EM converges to a particular point on such a ridge. For posterior modes, $D^{20}Q$ is replaced by $D^{20}Q + D^2 G$ (the text refers here to “equation (3·15)”, which is the score identity of Lemma 2; by our reading the rate formula (3.19) is meant), so an informative prior reduces the missing information and speeds convergence.

The paper lists three ways to speed things up (p. 10): a Newton step with second derivatives estimated from the iterates (Aitken’s acceleration for one parameter) or from (3.19) or (3.26); a partial M step that only increases $Q$, which is GEM; and, for missing-data patterns whose likelihood factors into pieces with distinct parameters (Rubin, 1974), EM on each factor, which converges faster because each factor has less missing data.

Examples in Section 4

Section 4 states for each class of problems the model and enough detail to derive the E and M steps, but, by the authors’ own description, neither studies the algorithms in detail nor their rates. The seven subsections:

• 4.1 Missing data, under two stated conditions: the parameters are a priori independent of those of the missing-data process, and the data are missing at random in Rubin’s (1976) sense, which roughly excludes values missing because of what they would have been. For multinomial data the E step assigns fractional counts to the cells a unit might belong to; for the normal linear model it fills in missing responses with their expectations, the Healy–Westmacott procedure; for the multivariate normal it cannot simply fill in values, because some sufficient statistics are quadratic: the missing parts of the sums of squares and products are replaced by their conditional expectations, which include the conditional covariances of the missing values.
• 4.2 Grouping, censoring and truncation. For truncated data the paper proposes a further completion: the unknown number $m$ of truncated draws is given a negative-binomial distribution (4.2.3), chosen so that the observed-data likelihood is unchanged and the complete data become a plain random sample of size $n + m$; the expected complete sample size is then $n / P(\boldsymbol{\phi}^{(p)})$ (4.2.8). The probit model of quantal assay is given as an example of censoring at varying points.
• 4.3 Finite mixtures. Unobserved indicator vectors $\mathbf{z}_i$ of the $R$ states make the complete-data log-likelihood linear in the indicators, so the E step replaces them by the current posterior state probabilities and the M step fits each component with every observation weighted by them; the hidden-Markov-chain case of Baum and co-authors replaces independent indicators by a Markov chain.
• 4.4 Variance components. For random-effects models the E step computes conditional expectations of $\mathbf{x}_i^\top \mathbf{x}_i$ and the M step is $\sigma_i^{2(p+1)} = t_i^{(p)} / r_i$ (4.4.8); treating the fixed effects as having infinite variance gives estimates identical to restricted maximum likelihood (REML).
• 4.5 Hyperparameter estimation, taking $\mathbf{x} = (\mathbf{y}, \boldsymbol{\theta})$ so that EM gives maximum likelihood estimates of the hyperparameters of a prior, as used in empirical Bayes.
• 4.6 Iteratively reweighted least squares. When $(y_i - \mu)\sqrt{q_i}/\sigma$ is standard normal given an unobserved $q_i$, as for the $t$ distribution, EM becomes weighted least squares with weights $w_i = \mathbb{E}(q_i \mid y_i, \mu^{(p)}, \sigma^{(p)})$ recomputed at each step (4.6.2).
• 4.7 Factor analysis, treating the factor scores as missing, so that EM alternates least-squares regressions on the cross-products matrix of variables and factors; the authors write that the method has apparently not been proposed before.

The paper has no experiments beyond Table 1. Its evidence is the theory and the unification: the same two steps, the same monotonicity proof and the same rate theory apply to every subsection.

The rate theory on a clipped image

Section 4.2’s censoring and Theorem 4’s rate can be checked together on an image. The script renders a uniformly lit grey wall at three exposures with sensor noise of standard deviation 30, before quantization, and clips every value at 255, so a saturated pixel reports only that its intensity is at least 255, a right-censored observation. EM’s E step replaces each clipped pixel’s $z$ and $z^2$ by their conditional expectations above 255 under the current Gaussian, and the M step is the complete-data estimate (variance with denominator $n$). The script compares the result with direct maximization of the censored likelihood, then compares the observed ratio of successive errors in the mean with the eigenvalues of $DM$ computed from Eq. (3.26), with $\mathbf{t} = (\sum z, \sum z^2)$. It uses NumPy 2.3.5 and SciPy 1.16.3.

import numpy as np
from scipy import stats, optimize

# A 128 x 128 patch of a uniformly lit grey wall, photographed at three exposures.
# True intensity mu with sensor noise sigma = 30 (values before quantization);
# the sensor clips at c = 255, so a bright pixel reports only "at least 255".
rng = np.random.default_rng(0)
c, sigma_true = 255.0, 30.0
noise = rng.standard_normal(128 * 128)

def trunc_moments(mu, sd, k):
    # E[z^j | z > c] for j = 1..k under N(mu, sd^2), via scipy's truncated normal
    a = (c - mu) / sd
    return [stats.truncnorm.moment(j, a, np.inf, loc=mu, scale=sd) for j in range(1, k + 1)]

for mu_true in (220.0, 240.0, 260.0):
    z = mu_true + sigma_true * noise
    obs = z[z < c]                      # fully reported pixels
    n, m = z.size, int((z >= c).sum())  # m pixels clipped at 255

    # EM (DLR Section 4.2): the E step replaces each clipped pixel's z and z^2 by
    # their conditional expectations given z > c; the M step is the complete-data MLE.
    mu, sd = obs.mean(), obs.std()      # start from the unclipped pixels
    path = [mu]
    for it in range(500):
        e1, e2 = trunc_moments(mu, sd, 2)
        mu_new = (obs.sum() + m * e1) / n
        sd = np.sqrt((np.sum(obs**2) + m * e2) / n - mu_new**2)
        path.append(mu_new)
        if abs(mu_new - mu) < 1e-12:
            break
        mu = mu_new
    mu_s, sd_s = mu_new, sd

    # Cross-check: maximize the censored-data log-likelihood directly.
    nll = lambda p: -(stats.norm.logpdf(obs, p[0], p[1]).sum() + m * stats.norm.logsf(c, p[0], p[1]))
    ref = optimize.minimize(nll, [obs.mean(), obs.std()], method="Nelder-Mead",
                            options={"xatol": 1e-10, "fatol": 1e-12, "maxiter": 5000}).x

    # Predicted rate, DLR Eq. (3.26): DM = V(t | y, phi*) V(t | phi*)^{-1}, t = (sum z, sum z^2).
    def cov_t(moms):                    # covariance of (z, z^2) from raw moments 1..4
        m1, m2, m3, m4 = moms
        return np.array([[m2 - m1**2, m3 - m1 * m2], [m3 - m1 * m2, m4 - m2**2]])
    full = [stats.norm.moment(j, loc=mu_s, scale=sd_s) for j in range(1, 5)]
    V_t = n * cov_t(full)               # unconditional: every pixel unknown
    V_ty = m * cov_t(trunc_moments(mu_s, sd_s, 4))  # given y: only clipped pixels vary
    lam = np.sort(np.linalg.eigvals(V_ty @ np.linalg.inv(V_t)).real)[::-1]

    dev = np.abs(np.array(path) - mu_s)
    ratios = dev[1:] / dev[:-1]
    print(f"true mean {mu_true:.0f}: {m / n:5.1%} of pixels clipped, {len(path) - 1} EM iterations")
    print(f"  EM estimate  mean {mu_s:7.3f}, sd {sd_s:6.3f};  direct MLE {ref[0]:7.3f}, {ref[1]:6.3f}")
    print(f"  observed ratio of successive errors in the mean (iterations 8-11): "
          + " ".join(f"{r:.3f}" for r in ratios[7:11]))
    print(f"  eigenvalues of DM from Eq. (3.26): {lam[0]:.3f}, {lam[1]:.3f}")

Output:

true mean 220: 12.1% of pixels clipped, 19 EM iterations
  EM estimate  mean 220.213, sd 29.970;  direct MLE 220.213, 29.970
  observed ratio of successive errors in the mean (iterations 8-11): 0.193 0.193 0.193 0.193
  eigenvalues of DM from Eq. (3.26): 0.193, 0.001
true mean 240: 31.6% of pixels clipped, 38 EM iterations
  EM estimate  mean 240.378, sd 30.197;  direct MLE 240.378, 30.197
  observed ratio of successive errors in the mean (iterations 8-11): 0.456 0.457 0.457 0.457
  eigenvalues of DM from Eq. (3.26): 0.457, 0.004
true mean 260: 56.5% of pixels clipped, 92 EM iterations
  EM estimate  mean 259.946, sd 29.749;  direct MLE 259.946, 29.749
  observed ratio of successive errors in the mean (iterations 8-11): 0.724 0.725 0.726 0.727
  eigenvalues of DM from Eq. (3.26): 0.728, 0.026

Because the model knows that each clipped value lies above 255, the censoring is modelled rather than ignored as values missing at random, the case in which Expectation–Maximization shows clipped pixels biasing the estimates. At all three exposures EM lands on the same mean and standard deviation as direct optimization of the censored likelihood, even when more than half the wall is saturated and the true mean lies above the clipping level. The ratio of successive errors settles on the largest eigenvalue that Eq. (3.26) predicts from two covariance matrices at the solution, as Table 1 settles on 0.1328: 0.193, 0.457 and, still approaching it after eleven iterations, 0.728. The rate grows with the exposure and so does the iteration count, from 19 to 92 to reach a change below $10^{-12}$. The rate is not the fraction of pixels lost: it is above it at every exposure (0.193 against 12.1%), while the second eigenvalue is close to zero. By our reading, the clipped pixels are the brightest ones, which carry more than their share of the information about the combination of mean and spread along the slow direction, and almost none is lost in the other direction. That is the paper’s point about components converging at different speeds, and Dr Little’s in the discussion: what sets the speed is the fraction of missing information, which depends on how the complete data are chosen, not simply on how much data is missing.

The Discussion

The paper was read at a meeting of the Society’s Research Section chaired by S. D. Silvey, and the printed discussion, with contributions received in writing afterwards and the authors’ reply, fills pages 22–38. Its main themes:

• Terminology and generality. E. M. L. Beale, opening, praised the technical content but was less enthusiastic about the new terminology, and judged the separation into E and M steps in the general form after Eq. (2.17) impractical; J. A. Nelder, seconding, made a similar point, since the general E step requires a whole function of continuous variables. B. Torsney suggested that in Scotland the technique would be called “a wee GEM”.
• Speed. Nelder, D. A. Preece, S. C. Pearce and M. J. R. Healy described accelerated “stretching” versions of the missing-plot procedure, and Healy noted that EM on designed experiments is numerically the Jacobi method, slower than Gauss–Seidel. R. Sundberg reported cases from mixtures and convolutions with a rate above 0.98, too slow in his judgement for practical use, and advocated Newton–Raphson or scoring there. R. Thompson gave an approximate rate for a one-way random-effects model that tends to one as the between-group variance vanishes. S. J. Haberman wrote Newton–Raphson and scoring in the same moment notation and found EM often painfully slow in latent-structure models.
• Standard errors. R. J. A. Little contrasted EM with scoring, which converges quadratically and gives an asymptotic covariance matrix; C. A. B. Smith showed that for the genetics example a binomial variance, computed with the effective count of about 101.83 animals and divided by $1 - 0.1328$, gives a standard error of 0.051, and W. H. Carter and H. O. Hartley pointed to Hartley’s own methods for obtaining variances from the iterates.
• Choosing the complete data. Little showed that adding fewer imaginary observations to an unbalanced two-way table reduces the missing information and hence speeds convergence.
• Multiple maxima. G. D. Murray gave twelve observations from a bivariate normal with four values of each variable missing whose likelihood has a saddle point at zero correlation and two maxima at correlation $\pm\tfrac12$, and reported that about 5 per cent of simulated data sets converged to different local maxima from two starting points. S. E. Fienberg, citing Haberman (1974), noted that an incomplete-data likelihood need not be concave even when the complete-data one is, and can have its maximum on the boundary of the parameter space, as in the Heywood cases of factor analysis.
• Priority. L. E. Baum wrote that Theorem 2.1 of Baum et al. (1970), used with an EM transformation, contains the paper’s Lemma 1 and Theorem 1. B. Efron pointed out that Fisher (1925) had already used the score identity behind Eq. (2.13). Hartley described himself as a minstrel who had sung the song for 18 years, and argued that his 1958 method covered any grouped discrete distribution, that his treatment of truncation fell outside the paper’s framework, and that his 1971 convergence conditions with Hocking were weaker than Theorems 2–4.

In their reply the authors agreed that strict separation of the E and M steps is numerically feasible only when the M step is easily computed from the expectations of a finite number of sufficient statistics, and said they kept the general E step for the statistical principle it expresses. They accepted that Newton–Raphson converges quadratically near a maximum but noted that it does not always increase the likelihood, can move towards a local minimum, and needs a second-derivative matrix whose repeated evaluation or storage can be infeasible. On priority they acknowledged Hartley as the originator of many of the ideas reviewed, while rejecting his claim about truncation by pointing to §4.2, and said they had not meant to question the mathematical generality of Baum and co-authors, only their recognition of its range of statistical application. On ridges they maintained that their convergence results hold where Hartley and Hocking’s do not; Wu’s later analysis bears on exactly this claim.

Impact

Wu (1983, §1) summarized the paper’s three contributions as recognizing the E and M steps in their general forms, giving theoretical properties of the algorithm, and recognizing a wide range of applications in statistics; he also observed that EM had become very popular despite slow numerical convergence, because the steps are easy to implement and need little storage. Neal and Hinton (1998, §1) date special cases of the algorithm back several decades and say its use grew further once Dempster, Laird and Rubin set out its generality. The paper is widely cited, and later work refers partial-M-step variants to it under the name GEM (Neal & Hinton, 1998, §1). In vision, EM fits the colour mixtures, background models and motion layers described in Expectation–Maximization and Gaussian Mixture Model; in our reading, those methods use the paper’s §4.3 formulation, with each pixel’s component indicator as the missing data.

Limitations

The paper and its discussants named most of the limitations that later work studied:

• No implementation. The general E step is a function over the whole parameter space, and the paper specifies neither its computation nor efficient mixing of the two steps; the referee’s footnote and the opening discussants made this point, and the authors conceded it.
• Linear convergence. Theorem 4 makes the rate the largest fraction of missing information, which can be close to one; the acceleration methods suggested in §3 were not developed, and the reply admits that routine acceleration for many parameters was unclear.
• No standard errors. EM does not produce the information matrix that Newton–Raphson and scoring compute along the way.
• Local maxima and saddle points. Theorem 1 guarantees only a non-decreasing likelihood, and the paper’s own discussion of saddle points is qualitative; Murray’s example showed EM stopping at one.
• Examples sketched, not tested. Section 4 gives no computations, and S. R. Searle questioned several statements in the variance-components treatment and saw no advantage for EM over methods designed for those models.

The Error in Theorem 2

Wu (1983, §1) found that the convergence proof contains an error: the step from (3.13) to (3.14) in Theorem 2 fails, through an incorrect use of the triangle inequality. Inequality (3.13) bounds the sum of the squared lengths of the steps; (3.14) needs a bound on the squared length of their sum. The triangle inequality bounds the length of the sum by the sum of the lengths, not of the squared lengths, and the two differ: $r$ equal steps of length $s$ in one direction have squared lengths summing to $r s^2$ but a total displacement whose square is $r^2 s^2$ (our illustration). Small squared steps can therefore add up to an unbounded path. Because Theorem 3 reaches its convergence conclusion by verifying condition (2) of Theorem 2, Wu concluded that the convergence of EM sequences as proved in Theorems 2 and 3 is cast in doubt, while Theorems 1 and 4, on the monotone likelihood and the rate, remain valid.

Wu’s §2.2 adds that the argument would work if condition (2) bounded the gain in $Q$ by $\lambda \lVert \boldsymbol{\phi}^{(p+1)} - \boldsymbol{\phi}^{(p)} \rVert$ instead of its square, but then the gradient at any limit would have norm at least $\lambda$, so the limit could not be a stationary point; he saw no simple way to repair the proof. He reports that Boyles (1983) gave a GEM sequence, though not an EM one, as a counterexample to Theorem 2, and calls the claim on p. 10 that EM converges to a particular point of a ridge rather than moving along it indefinitely unwarranted: condition (2) makes the squared step lengths summable, but the step lengths themselves can still sum to infinity. He had no counterexample for EM, but expected that convergence on a ridge needs stronger conditions than the paper assumed.

Instead of patching the proof, Wu treated EM as an optimization algorithm and used Zangwill’s global convergence theorem. Under assumptions that include compact level sets of the likelihood and continuity of $Q$ in both arguments, all limit points of an EM sequence are stationary points and the likelihood converges to a stationary value (his Theorem 2); a further condition makes them local maxima (Theorem 3); and the parameters themselves converge when only one stationary point has the limiting likelihood value, when the steps shrink to zero and such points are isolated, or when the likelihood is unimodal with a single stationary point (Theorems 4–6, Corollary 1). Murray had reported in the discussion that EM started at zero correlation converges to the saddle point of his example; Wu explained this as a failure, at that point, of the extra condition his Theorem 3 needs for local maxima (§2.1). Wu’s results, and what they mean for EM in practice, are set out in Expectation–Maximization.

What Came After

Wu (1983) replaced Theorems 2 and 3 with results proved under conditions that can be checked, and recommended running EM from several starting points representative of the parameter space (§3). Neal and Hinton (1998, §1) list the expectation conditional maximization (ECM) algorithms of Meng and Rubin (1992), which split the M step into several conditional maximizations, as a widely applicable subclass of GEM. Their own paper recast EM as maximization of one function, resembling negative free energy, over the parameters in the M step and over the distribution of the unobserved variables in the E step; that view justified incremental and sparse E steps, which the paper’s definitions did not cover, and is the lower-bound derivation that Expectation–Maximization uses and that the evidence lower bound generalizes. k-means clustering is the hard-assignment limit of EM for a Gaussian mixture, and the paper’s factor-analysis model, with its diagonal residual covariance restricted to a single variance, becomes the probabilistic form of principal component analysis, as those articles explain.