Gaussian Mixture Models

Chapter 1: The GMM

A Gaussian Mixture Model is a weighted sum of K Gaussian densities. Each component k has its own mean μ_k, covariance Σ_k, and mixing weight π_k. The overall density at a point x is:

where N(x | μ, Σ) is the multivariate Gaussian density:

The mixing weights π_k satisfy two constraints: each weight is non-negative (π_k ≥ 0), and they sum to one (∑ π_k = 1). Think of π_k as the probability that a randomly drawn data point came from component k. If π₁ = 0.7, then 70% of the data was generated by the first Gaussian.

In 1D, each Gaussian is a bell curve. In 2D, each is an ellipse (the shape determined by its covariance matrix). A GMM in 2D looks like a landscape of overlapping hills, each centered on a different μ_k, each with its own tilt and spread from Σ_k.

Chapter 2: Log-Likelihood

Given N data points x₁, …, x_N, we want to find the parameters θ = {π_k, μ_k, Σ_k}_k=1..K that best explain the data. The standard approach: maximum likelihood estimation. We maximize the probability of observing our data under the model.

Assuming data points are i.i.d., the likelihood is a product over data points:

Taking the log (to turn products into sums) gives the log-likelihood:

Here is the problem. In standard maximum likelihood for a single Gaussian, the log passes through to the exponent and everything simplifies to a clean closed-form solution. But here the log sits outside a sum over k. The log of a sum does not simplify. Setting derivatives to zero gives equations where every parameter depends on every other parameter in a tangled mess.

There is a second subtle issue: the likelihood surface has singularities. If a component's mean μ_k collapses exactly onto a data point and its variance shrinks to zero, the density at that point goes to infinity. The likelihood becomes unbounded. In practice, we avoid this with regularization (adding a small diagonal to Σ_k) or by detecting degenerate components and resetting them.

Chapter 3: Responsibilities

At the heart of EM is a simple question: for each data point x_n, how much does each component k "claim" it? This is the responsibility r_nk — the posterior probability that component k generated point n.

We compute it with Bayes' theorem. The prior probability that point n came from component k is π_k. The likelihood of x_n given component k is N(x_n | μ_k, Σ_k). By Bayes' rule:

Each r_nk is between 0 and 1, and for each point n, the responsibilities sum to 1 across components: ∑_k r_nk = 1. A point near the center of cluster k gets r_nk close to 1. A point equidistant from two clusters might have r_nk ≈ 0.5 for each.

We also define the effective number of points assigned to component k:

This is a soft count. If all 100 points have r_n1 = 0.6, then N₁ = 60. It tells us how much data each component "owns." We will use N_k to normalize the M-step updates.

Chapter 4: Updating Means

Once we have the responsibilities (the E-step), we can update the parameters to better fit the data — this is the M-step. Let's start with the means. The new mean for component k is the responsibility-weighted average of all data points:

This makes perfect intuitive sense. Each point contributes to the mean of component k in proportion to how much component k claims it. A point with r_nk = 0.9 pulls the mean strongly. A point with r_nk = 0.01 barely matters. The normalizer N_k = ∑ r_nk ensures the weights sum to 1.

Why does this update increase the likelihood? Setting the derivative of the expected complete-data log-likelihood with respect to μ_k to zero gives exactly this formula. It is the MLE for a single Gaussian, but with each data point weighted by r_nk. If we knew the true assignments, fitting each cluster's mean would be trivial. The responsibilities act as a surrogate for those unknown assignments.

Chapter 5: Updating Covariances

The mean tells us where a cluster is; the covariance tells us its shape. Is it a tight circle? A wide ellipse tilted 45 degrees? The M-step updates each covariance Σ_k as a responsibility-weighted outer product of the residuals:

Read this carefully. For each data point, we compute its deviation from the new mean, form the outer product (which gives a matrix capturing direction and magnitude of deviation), and weight it by the responsibility. Sum them up, normalize by N_k, and you have the new covariance.

There is a practical concern: if too few points are assigned to a component (N_k very small), the covariance estimate becomes unreliable and can collapse. A common fix is to add a small regularization term: Σ_k^new ← Σ_k^new + ε I, where I is the identity matrix and ε is tiny (e.g. 10⁻⁶). This prevents singularities and keeps the matrix invertible.

python
def update_covariance(X, mu_k, r_k, N_k, eps=1e-6):
    # X: (N, D), mu_k: (D,), r_k: (N,), N_k: scalar
    diff = X - mu_k                       # (N, D)
    weighted = r_k[:, None] * diff        # (N, D) weighted residuals
    Sigma_k = (weighted.T @ diff) / N_k   # (D, D)
    Sigma_k += eps * np.eye(X.shape[1])  # regularize
    return Sigma_k

Chapter 6: Updating Weights

The mixing weights π_k represent how prevalent each cluster is. If cluster 1 claims most of the data (high N₁), its weight should be large. The M-step update is the simplest of the three:

That's it. The new weight for component k is just the fraction of the total (soft) data count assigned to it. If N = 100 data points and N_k = 30 (counting soft assignments), then π_k = 0.30. The constraint ∑ π_k = 1 is automatically satisfied because ∑ N_k = N (every point's responsibilities sum to 1 across components).

Here is the complete set of M-step updates, side by side:

Each formula is a soft-weighted version of the maximum likelihood estimator for a single Gaussian. If the responsibilities were hard (0 or 1), these would be exactly the MLE formulas for each component independently. The soft assignments gracefully interpolate.

Chapter 7: The EM Algorithm

Now we put it all together. The Expectation-Maximization algorithm for GMMs alternates between two steps until convergence:

Initialization matters. EM finds a local maximum, not necessarily the global one. Different starting points lead to different solutions. Common strategies:

When to stop: Track the log-likelihood after each iteration. When the increase drops below a threshold (e.g. |Δ ln L| < 10⁻⁴), declare convergence. Alternatively, stop after a fixed number of iterations (e.g. 100).

python
def em_gmm(X, K, max_iter=100, tol=1e-4):
    N, D = X.shape
    # Initialize
    mu = X[np.random.choice(N, K, replace=False)]
    Sigma = [np.eye(D) for _ in range(K)]
    pi = np.ones(K) / K

    for it in range(max_iter):
        # E-step: compute responsibilities
        r = np.zeros((N, K))
        for k in range(K):
            r[:, k] = pi[k] * multivariate_normal(mu[k], Sigma[k]).pdf(X)
        r /= r.sum(axis=1, keepdims=True)

        # M-step: update parameters
        Nk = r.sum(axis=0)
        for k in range(K):
            mu[k] = (r[:, k] @ X) / Nk[k]
            diff = X - mu[k]
            Sigma[k] = (r[:, k][:, None] * diff).T @ diff / Nk[k]
            Sigma[k] += 1e-6 * np.eye(D)
        pi = Nk / N

        # Check convergence (log-likelihood)
        ll = np.sum(np.log(r.sum(axis=1)))  # simplified
        if it > 0 and abs(ll - prev_ll) < tol: break
        prev_ll = ll
    return mu, Sigma, pi

Chapter 8: EM in Action

Time to see it work. Below is a 2D dataset with K = 3 true clusters. The algorithm starts with random initial means (the large circles) and spherical covariances. Click through the E-step and M-step to watch EM converge, or hit Auto-Run to animate the full process.

Try resetting and watching again. Different random initializations lead to different convergence paths — sometimes a component "steals" points from another, sometimes two components start near the same cluster and one drifts away. This is the local-optima issue in practice.

Chapter 9: Latent Variables

There is a deeper way to think about GMMs that connects to a broad class of models. Introduce a latent variable z_n for each data point. z_n is a one-hot vector of length K — it indicates which component generated x_n. If z_n = [0, 1, 0], then point n came from component 2.

The generative story becomes explicit:

The joint distribution of x and z is:

When z_k = 1 (point assigned to component k), only the k-th term survives. Marginalizing out z (summing over all possible assignments) recovers the GMM density: p(x) = ∑_k π_k N(x | μ_k, Σ_k). This is exactly where the mixture formula comes from.

This perspective reveals EM as a general algorithm. Any model where some variables are hidden fits the pattern: define the joint p(x, z | θ), compute the expected complete-data log-likelihood (E-step), and maximize it (M-step). The latent variables could be cluster assignments (GMMs), hidden states (HMMs), missing data (imputation), or topic assignments (LDA).

Chapter 10: Connections

Gaussian Mixture Models sit at a crossroads of several important ideas in machine learning. Let's map the connections.

Looking ahead: Chapter 12 takes a completely different approach to the classification problem. Instead of modelling the full density p(x), Support Vector Machines find the decision boundary directly. Where GMMs are generative (model how data is produced), SVMs are discriminative (model the boundary between classes). Both perspectives have their place.