Expectation-maximization algorithm

The expectation-maximization algorithm fits a latent-variable model by alternating a conditional expectation over the latent variables with a maximization over the statistical parameters. Its target is the observed-data likelihood function, although its updates use the simpler complete-data log-likelihood.^[ref]

Write $y$ for the observations, $z$ for the latent variables, and $p_{\theta }(y,z)$ for their joint probability density function or probability mass function. The observed-data log-likelihood iswhere the integral becomes a sum for discrete $z$. Directly maximizing a logarithm of this sum can be awkward. The E-step instead computes a conditional distribution at the current statistical parameters ${\theta }^{(t)}$; the M-step uses it to construct a new fit. The attractive idea is to retain uncertainty about the hidden labels while recovering simple weighted means and other complete-data updates. The Known-variance Gaussian mixture EM example shows the mechanism without estimating a variance or mixture weight.

Arthur Pentland Dempster, Nan M. Laird, and Donald B. Rubin gave the general framework its name in their 1977 paper, which also discusses earlier special-case algorithms.^[ref]

Hold the old statistical parameters fixed when taking the conditional expectation, while keeping the candidate statistical parameters $\theta$ variable:For independent observations in a finite mixture model, $Z_i$ is a hidden component label. Bayes' theorem gives its mixture responsibilities:These posterior probabilities are the conditional expectations of the mixture allocation indicators; they turn hidden counts into fractional counts.^[ref] They condition on a current fixed fit, rather than averaging over a posterior distribution of the statistical parameters. Taking the conditional expectation of the complete-data log-likelihood is essential: substituting a single guessed label, or merely substituting the conditional mean into a nonlinear complete-data expression, generally gives a different algorithm.

Maximize $Q(\theta \mid {\theta }^{(t)})$ over the permitted statistical parameters, with the E-step conditional distribution held fixed.^[ref] A finite mixture model hasFor normal distributions with known common variance ${\sigma }^2>0$ and fixed positive mixture weights, only the means are fitted. Dropping terms independent of these candidate means givesThe derivative in ${\mu }_j$ is ${\sum }_i{r}_{ij}^{(t)}(y_i-{\mu }_j)/{\sigma }^2$, and its second derivative is $-N_j^{(t)}/{\sigma }^2$. Thus a positive weighted count gives a unique maximizing weighted mean. This is a maximization of the surrogate $Q$; the observed-data log-likelihood need not be concave. Recompute the mixture responsibilities after updating the means. The broader EM for Gaussian mixtures with a common variance also estimates mixture weights and the shared variance.

Choose six illustrative values $y=(-2.2,-1.8,-1,1,1.8,2.2)$, fixed mixture weights ${\pi }_1={\pi }_2=1/2$, known standard deviation $\sigma =1$, and initial means $({\mu }_1^{(0)},{\mu }_2^{(0)})=(-1,1)$. These are constructed data for a Gaussian mixture distribution, not measurements. With $\varphi$ the standard normal distribution probability density function, the fitted log-likelihood isAt the initial means, the component density ratio is $\varphi (y_i-1)/\varphi (y_i+1)=e^{2y_i}$. The first E-step therefore gives

The paired values $y$ and $-y$ have complementary mixture responsibilities, so $N_1^{(0)}=N_2^{(0)}=3$ exactly. The first M-step givesThe hyperbolic tangent expression follows by pairing the contributions at $y$ and $-y$. Assigning every negative value wholly to the first mixture component and every positive value wholly to the second would instead give means $\pm 5/3$. Those hard assignments discard the overlap: the expectation-maximization algorithm is not the same update as K-means clustering.

The observed-data log-likelihood rises from $-11.420338905$ initially to $-10.368745063$ after one update and $-10.340655142$ after two. Further updates approach means $\pm 1.638023833$, with log-likelihood approximately $-10.340509617$. These values describe this initialization; EM likelihood monotonicity supplies no certificate of a global maximum.

For stable numerical evaluation, put $a_{ij}=\log (1/2)-\frac{1}{2}\log (2\pi )-\frac{1}{2}(y_i-{\mu }_j{)}^2$ and $b_i={\max }_j{a}_{ij}$. The log-sum-exp function givesSubtracting $b_i$ keeps at least one exponential equal to one, avoiding the all-zero sum that direct evaluation of tiny component densities can produce.

In the Known-variance Gaussian mixture EM example, initialize both means at the same value $m$. The mixture responsibilities are then $1/2$ for every observation. Each M-step sets both means to the sample mean, zero here; the next E-step again gives $1/2$. Thus $(0,0)$ is a fixed point with log-likelihood $-14.593631199$, despite the better fit found from $(-1,1)$.

One can check the failure without trusting an iteration plot. Separate the means as $(-a,a)$ and use the hyperbolic cosine identity:The Taylor expansion uses ${\sum }_i{y}_i^2=18.16$. Small nonzero separation therefore improves the log-likelihood, whereas translating both means to $m$ changes it by $-3m^2$. This stationary point has both improving and worsening directions. Identical initialization cannot discover the separation because every observation keeps identical label probabilities. Trying several distinct initial fits and comparing their final log-likelihoods is useful, but does not prove global optimality.

Keeping the variance fixed also removes a different pathology. With known $\sigma =1$, every component probability density function is at most $(2\pi {)}^{-1/2}$, so $\ell \le -3\log (2\pi )$. If separate component variances are allowed to shrink without restriction, one mixture component can concentrate at an observation while another retains positive density at the other observations, making the likelihood function unbounded. The fixed positive variance in this example prevents that collapse; it does not remove initialization traps.

For $N$ observations from each of two labeled exponential distributions, an E-step may produce positive completed time totals $U,V$. Its expected complete-data log-likelihood is $Q=N\log {\lambda }_1+N\log {\lambda }_2-U{\lambda }_1-V{\lambda }_2$. This is strictly concave. On the closed ordering constraint ${\lambda }_1\ge {\lambda }_2>0$, its maximum is $(N/U,N/V)$ when $U<V$; otherwise it lies at ${\lambda }_1={\lambda }_2=2N/(U+V)$. To see the boundary case, first maximize along the equality line. At that maximum, moving into either feasible rate-separating direction has directional derivative $(V-U)/2\le 0$, so concavity proves global optimality. For the open constraint ${\lambda }_1>{\lambda }_2$, this boundary value is a supremum rather than an attained M-step.

// Target: probability-and-statistics.bigb

For a multinomial distribution with two unobserved cells merged into one observed count, the conditional allocation is a binomial distribution with probabilities proportional to the two cell probabilities. The expectation-maximization algorithm replaces the missing counts by their conditional expectations in the parameter-dependent part of the complete-data log-likelihood. The expectation of a log-factorial normalization generally differs from its value at the expected count, but this difference is independent of the candidate parameter and does not affect the maximization step.

For independent normal coordinates, an unobserved first coordinate has conditional law equal to its current marginal normal distribution, even when the second coordinate of that observation is known. Thus its E-step first moment is ${\mu }_1^{(t)}$ and second moment is $({\mu }_1^{(t)}{)}^2+v_1^{(t)}$. With three observed first coordinates and one missing value, the mean update is $({\sum }_{i=1}^3{x}_{i1}+{\mu }_1^{(t)})/4$ and the variance update is $[{\sum }_{i=1}^3(x_{i1}-{\mu }_1^{(t+1)}{)}^2+v_1^{(t)}+({\mu }_1^{(t)}-{\mu }_1^{(t+1)}{)}^2]/4$. The conditional variance term is essential: filling in only the conditional mean is not the expectation-maximization algorithm.

For a common unconstrained Markov kernel on the edges of a rooted tree, let $N_{ab}$ be the conditional expectation of the number of transitions $a\to b$, given observed leaves and the old Markov kernel. The expectation-maximization algorithm maximizes ${\sum }_{a,b}{N}_{ab}\log K_{ab}$ over row probability distributions, giving $K_{ab}^{\mathrm{new}}=N_{ab}/{\sum }_c{N}_{ac}$. A row with zero total count is unrestricted by this objective. Belief propagation computes the counts exactly on a tree.

In the expectation-maximization algorithm for a stationary autoregressive process of order one, a single missing interior value has a normal distribution conditional on its neighbors. The E-step uses both its conditional mean and its conditional second moment, not just mean imputation. The M-step maximizes the expected stationary complete-data log-likelihood, including the initial-observation density and the conditional-variance contribution to the two adjacent innovation squares.

For $Q(\theta \mid {\theta }_0)={\mathbb{E}}_{{\theta }_0}[\log p_{\theta }(Y,Z)\mid Y]$, the Jensen inequality givesThus any expectation-maximization algorithm M-step that increases $Q$ cannot decrease the observed-data likelihood function.^[ref]

For a short proof, assume common support and finite expressions, and keep the E-step conditional distribution $q_t(z)=p_{{\theta }^{(t)}}(z\mid y)$ fixed. The concavity of the logarithm and Jensen's inequality giveThis evidence lower bound equals $Q(\theta \mid {\theta }^{(t)})-{\mathbb{E}}_{q_t}\log q_t(Z)$. At $\theta ={\theta }^{(t)}$, the ratio inside the expectation is the constant $p_{{\theta }^{(t)}}(y)$, so equality holds. During the M-step, the second term is constant; increasing $Q$ raises a lower bound that touches the current log-likelihood. HenceIncreasing likelihood values alone do not establish convergence of the statistical parameters or attainment of a global maximum. In a differentiable finite mixture model, the gradient of $\ell$ equals the candidate-parameter gradient of $Q$ evaluated at the current fit. Thus an interior fixed point of an exact M-step is a stationary point, which need not be a local maximum; the Identical-mean EM trap is an explicit counterexample. Small changes in log-likelihood can also occur before the means have settled, so inspecting both is useful.