The EM Algorithm for Gaussian Mixtures

Why Gaussian Mixtures Need an Iterative Algorithm

A Gaussian mixture describes a profile as arising from one of several Gaussian components, while leaving its component membership unobserved. The mixture likelihood sums over possible memberships for every profile. That sum makes direct maximization difficult: the component memberships that would make the Gaussian parameter estimates straightforward are not known, while their probabilities depend on those same parameters.

The expectation-maximization (EM) algorithm handles this dependence by alternating between estimating membership probabilities and updating the component parameters. These are called the E step and the M step. The steps can be derived for a Gaussian mixture and guarantee that the observed-data likelihood does not decrease under appropriate conditions. They do not, however, guarantee that the algorithm finds the best solution among all possible parameter values.

Model and Notation

Suppose \(x_1,\ldots,x_n\) are independent profile vectors in \(p\) dimensions. A \(K\)-component Gaussian mixture has component weights \(\pi_k\), means \(\mu_k\), and covariance matrices \(\Sigma_k\), for \(k=1,\ldots,K\). The weights satisfy \(\pi_k>0\) and \(\sum_{k=1}^K\pi_k=1\); each \(\Sigma_k\) is positive definite. Conditional on an unobserved component label \(z_i=k\), the model assumes \(x_i\sim N(\mu_k,\Sigma_k)\). The marginal density of a profile is therefore a weighted sum of component densities.

Here \(\phi_p\) is the \(p\)-dimensional Gaussian density and \(\theta\) denotes all mixture parameters. Introduce an indicator \(z_{ik}\), equal to 1 if profile \(i\) belongs to component \(k\), and 0 otherwise. If these indicators were observed, the complete-data log likelihood would be

The EM algorithm replaces each unknown indicator by its conditional expectation under the current parameters. Write \(\theta^{(t)}\) for the parameters at iteration \(t\), and let \(r_{ik}^{(t)}=P(z_{ik}=1\mid x_i,\theta^{(t)})\). These conditional probabilities, also called responsibilities, are fractional rather than hard cluster assignments.

The M Step: Updating Mixture Parameters

With responsibilities held fixed, the expected complete-data log likelihood is

The M step maximizes this expression over the mixture parameters. Define the effective membership count \(N_k^{(t)}=\sum_i r_{ik}^{(t)}\). For unrestricted component means and positive-definite covariance matrices, the updates are as follows when each weighted scatter matrix is positive definite; if a weighted scatter matrix is singular, the unrestricted M-step may have no finite maximizer, so suitable covariance constraints or regularization are needed:

These formulas use fractional counts and weighted means and covariances. The covariance denominator is \(N_k^{(t)}\), not \(N_k^{(t)}-1\): this is a maximum-likelihood update, not an unbiased sample covariance estimate. If a component has zero effective count, these formulas are undefined, and the implementation must handle that component, for example by reinitializing it. Constrained covariance structures require maximizing \(Q\) subject to those constraints; the unrestricted formulas need not apply unchanged.

Why the Observed-Data Likelihood Does Not Decrease

The monotonicity guarantee follows by expressing the observed-data log likelihood as a lower bound plus a nonnegative divergence. Let \(q_{ik}\) be any probability distribution over the component labels for profile \(i\), so \(q_{ik}\geq0\) and \(\sum_k q_{ik}=1\). For positive joint terms \(a_{ik}=\pi_k\phi_p(x_i;\mu_k,\Sigma_k)\), Jensen's inequality gives

The inequality applies as written when \(a_{ik}>0\) wherever \(q_{ik}>0\). In a Gaussian mixture with positive weights and positive-definite covariance matrices, all component densities are positive, so this condition holds. The right-hand side, summed across profiles, is a lower bound on the log likelihood. Equality holds when \(q_{ik}\) is the posterior probability of component \(k\) for profile \(i\) under the parameters being evaluated.

The E step sets \(q_{ik}\) to those posterior probabilities at the current parameters, making the bound tight there. In the M step, the entropy terms involving only \(q\) stay fixed, while maximizing \(Q\) increases or preserves the parameter-dependent part of the bound. The bound is therefore at least as large at the updated parameters as it was at the old parameters. Since the old bound equaled the old log likelihood and the new log likelihood is at least its new bound, the new log likelihood cannot be smaller.

This justifies a useful diagnostic: for an exact E step and exact M step under the stated model, the observed-data log likelihood should not fall from one iteration to the next. A substantial decrease suggests a coding error, an inconsistent likelihood calculation, or an update that is not actually maximizing the stated \(Q\).

Convergence, Local Optima, and Practical Failure Modes

Nondecreasing likelihood is not the same as guaranteed convergence to the global maximum. EM is a local iterative procedure: its updates depend on the starting weights, means, and covariance matrices. Under regularity conditions, limit points of exact EM iterations are stationary points of the likelihood, but a stationary point can be a local maximum or a saddle point. Different initializations can therefore produce different solutions. Mixture labels are also interchangeable: permuting component labels does not change the fitted mixture or its likelihood.

A further issue is covariance collapse. With unrestricted Gaussian covariances, the likelihood of a mixture can be unbounded. For example, a component can center its mean on one observed profile while its covariance shrinks toward zero. In one dimension, the density of an observation at the component mean is \(1/(\sqrt{2\pi}\sigma)\), which grows without bound as \(\sigma\) approaches zero. A mixture component with positive weight can consequently make the likelihood diverge, rather than converge to a finite global maximum. This can appear as a covariance becoming extremely small, a component effective count approaching one, or rapidly increasing likelihood with unstable parameter estimates.

In practice, monitor the observed-data log likelihood, parameter changes, effective component counts, and covariance eigenvalues. A small change in likelihood is a reasonable stopping rule, but it only indicates that the algorithm has nearly stopped improving from its current trajectory. It does not establish global optimality. Multiple starts can expose sensitivity to initialization; covariance constraints or a structured covariance model can reduce degeneracy, though they alter the optimization problem and may require different M-step updates.

Compared with hard-assignment procedures such as Lloyd's algorithm for K-means, Gaussian-mixture EM retains uncertainty in membership and estimates component covariance as well as location and weight. Hard assignment is simpler and can be faster, but it does not optimize the same likelihood and cannot represent overlapping component membership in the same way. Direct numerical optimization of the observed likelihood is another alternative; it avoids the E-M decomposition but still faces nonconvexity and can be difficult to optimize. EM is attractive because its steps have a clear conditional-probability interpretation and a monotonicity guarantee, not because it removes the problem of local solutions.