Gaussian Mixture Models and the EM Lower Bound

A mixture model treats the data as a sample from a weighted sum of component densities. For k Gaussian components in dimension d, the density of a point x is p(x) = sum_k pi_k N(x | mu_k, Sigma_k), with weights pi_k that are positive and sum to one. Fitting means maximizing the total log-likelihood

$$\log L = \sum_{i=1}^{n} \log \left( \sum_{j=1}^{k} \pi_j , N(x_i \mid \mu_j, \Sigma_j) \right)$$

and the log outside the sum blocks the closed-form solution that a single Gaussian has. Each component’s mean depends on which points belong to it, and which points belong to it depends on all the other components, so the equations are coupled. The parameters number k d + k d(d+1)/2 + (k - 1) for full covariances: k d means, k d(d+1)/2 covariance entries, and k - 1 free weights. With d = 20 and k = 5, that is 1,195 parameters.

A latent variable turns the coupling into two tractable steps

Introduce a one-hot latent indicator z_i naming the component that generated point i. The complete-data log-likelihood, which pretends the indicators are known, is

$$\sum_i \sum_j z_{ij} \left[ \log \pi_j + \log N(x_i \mid \mu_j, \Sigma_j) \right]$$

and maximizing it is trivial: the weights become sample proportions and each component’s mean and covariance are the ordinary estimates over its assigned points. The indicators are not known, so fix any distribution q(z_i) over them and apply Jensen’s inequality to each term of the log-likelihood:

$$\log \sum_z q(z) \frac{p(x, z)}{q(z)} \ge \sum_z q(z) \log \frac{p(x, z)}{q(z)}$$

The right-hand side is a lower bound on log L for any q, and the slack is exactly the Kullback-Leibler divergence between q and the true posterior p(z | x). Two facts follow. The bound is tight when q equals the posterior, which happens to be computable in closed form for a mixture because

$$p(z_{ij} = 1 \mid x_i) = \frac{\pi_j N(x_i \mid \mu_j, \Sigma_j)}{\sum_l \pi_l N(x_i \mid \mu_l, \Sigma_l)}$$

and second, given the bound, the parameters can be improved by maximizing the expected complete-data log-likelihood under that posterior. The expectation step computes the responsibilities gamma_ij above; the maximization step solves the weighted estimation problem, where the weighting enters the normal equations.

The maximization step is a weighted mean and a weighted covariance

With N_j = sum_i gamma_ij, the effective number of points in component j, the updates are

• pi_j = N_j / n
• mu_j = (1 / N_j) sum_i gamma_ij x_i
• Sigma_j = (1 / N_j) sum_i gamma_ij (x_i - mu_j)(x_i - mu_j)^T

The weight N_j shows how much the component is used. A component whose responsibilities decay toward zero contributes nothing, and its parameters drift toward the last values they held, so implementations either keep it and report it as collapsed or drop it. Nothing enforces a minimum component size, unlike the hard partition of a centroid method.

Since the expectation step maximizes the bound over q and the maximization step maximizes it over the parameters, the bound is non-decreasing, and because it is tight after each expectation step, the log-likelihood is non-decreasing as well. The iteration therefore converges to a stationary point, but not necessarily to a global maximum: the objective has many stationary points, and different initializations reach different ones. Convergence is linear in general, and slow when components overlap, because the responsibilities stay ambiguous and each step moves the parameters a little.

k-means is the hard limit of the same iteration

Fix every covariance to epsilon I with equal weights and let epsilon tend to zero. The responsibilities become numerically dominated by the smallest Mahalanobis distance and collapse onto a one-hot vector, so the expectation step becomes a nearest centroid assignment, and the maximization step becomes the arithmetic mean. That is exactly the assignment and update pair from Lloyd’s Algorithm and the Local Minima of k-Means, which is why the two methods share their local-minimum behaviour and their dependence on initialization.

The intermediate case, sometimes called hard EM, assigns each point to its most likely component and then re-estimates, keeping some locality but discarding the soft weights. It is cheaper per iteration and is the basis of the classification variant of the algorithm. The full soft version retains the information a point carries about several components at once, which is what allows overlapping groups to be described by shared membership rather than by a boundary.

Degenerate solutions and the arithmetic that avoids them

The likelihood is unbounded above. Shrink one component’s covariance toward zero while placing its mean exactly on a data point: the density at that point grows without bound as the determinant of the covariance shrinks, while all other points contribute negligibly. So no finite parameter set maximizes the likelihood, and an unguarded iteration can drive a component onto a single point, reporting a very large likelihood and a useless model. The guard used in practice is a floor on the covariance diagonal, a small constant added to Sigma_j before inversion, typically on the order of 1e-6 in the units of the data; setting it too small leaves the degeneracy reachable, and setting it too large washes out genuine structure. A conjugate prior on the covariances achieves the same thing more principledly, turning the maximization into a maximum a posteriori step.

The arithmetic itself needs care. Gaussian log-densities involve -d/2 log(2 pi) - 1/2 log det(Sigma) - 1/2 (x-mu)^T Sigma^{-1} (x-mu), and for moderate d a single density can be smaller than the smallest positive double, around 1e-308. Computing the mixture density as a sum of raw densities then underflows to zero and the log-likelihood becomes negative infinity. The stable form computes the unnormalized log-weights log pi_j + log N(x | mu_j, Sigma_j), subtracts their maximum, exponentiates, and adds back the maximum at the end. Each covariance is best factored once per maximization step with a Cholesky decomposition, giving O(k d^3) work for the factors plus O(n k d^2) per expectation step to evaluate the quadratic forms.

import numpy as np

def e_step(X, pi, mu, Sigma):
    # Sigma: (k, d, d); returns responsibilities (n, k)
    logw = np.empty((len(X), len(pi)))
    for j in range(len(pi)):
        L = np.linalg.cholesky(Sigma[j])
        diff = X - mu[j]
        # solve L y = diff^T  ->  squared Mahalanobis distance is sum(y*y)
        y = np.linalg.solve(L, diff.T)
        logdet = 2.0 * np.log(np.diag(L)).sum()
        logw[:, j] = (np.log(pi[j]) - 0.5 * (y * y).sum(0)
                      - 0.5 * logdet - 0.5 * X.shape[1] * np.log(2 * np.pi))
    m = logw.max(1, keepdims=True)                 # log-sum-exp for stability
    w = np.exp(logw - m)
    return w / w.sum(1, keepdims=True)

Choosing a mixture structure and comparing fits

Covariance structure is the main modelling lever. A shared, tied covariance for all components assumes every group has the same orientation and elongation; a diagonal covariance assumes features are independent within a group; a spherical one assumes a single variance per group. The parameter counts differ by orders of magnitude: k d(d+1)/2 for full, d(d+1)/2 for tied, k d for diagonal, k for spherical. A tied full covariance is equivalent to whitening the data by that covariance and then fitting spherical components, the transform described in Principal Component Analysis Through the Singular Value Decomposition.

Comparing mixtures with different k usually borrows a penalized likelihood. The Bayesian information criterion subtracts a complexity penalty from twice the log-likelihood: BIC = -2 log L + p log n, with p the parameter count above. Its justification assumes the maximum likelihood estimate is an interior point, which the degeneracy above violates, so a regularized covariance is a precondition rather than an optional refinement. Where the likelihood is a poor guide, the geometric criteria used for hard partitions in Cluster Validity, Stability, and Choosing the Number of Clusters can be applied to the responsibilities collapsed to a hard assignment, at the cost of discarding the uncertainty the model provides.

One structural property has no counterpart in centroid methods: the likelihood is invariant under relabelling of the components. Any permutation of the k triples (pi_j, mu_j, Sigma_j) yields the same density and the same likelihood, so k! parameter vectors describe every distinct fit. Comparisons of density or responsibility are unaffected; comparisons of parameters between runs are not, and the index of a component carries no identity across restarts. Restarts are therefore the usual remedy for the local-minimum problem, and initializing the components from a cheap hard partition, such as the outcome of a few centroid iterations, is the standard practical choice.

Finally, the mixture is a parametric density and can be used as one: evaluate p(x) to score new points, or compare it with the non-parametric alternative in Kernel Density Estimation and the Bandwidth Trade-off, which places a kernel on every observation instead of fitting k components. The mixture is the better choice when the components are real and few, and the kernel estimator is the better choice when the shape of the density matters more than any group structure.