Gaussian mixtures and EM

k-means answers “which cluster?” with a hard label and describes every cluster as a point. Both simplifications cost information: boundary points get confident labels they don’t deserve, and elongated clusters get described by a center that captures nothing of their shape. The fix is to do what this curriculum always does when a procedure feels ad hoc — write down the probabilistic model the procedure is secretly assuming, then do maximum likelihood properly. The model is the Gaussian mixture; the algorithm that maximizes it is EM, the most important algorithm in classical unsupervised learning and the conceptual ancestor of variational autoencoders.

The model: data with a latent label

Assume each point was generated by a two-stage story: first draw a hidden component $z_i \in \{1..k\}$ with probabilities $\pi_1, \ldots, \pi_k$; then draw the point from that component’s Gaussian:

Clustering becomes inference of the latent $z$, and Bayes’ theorem (lesson 4 of Probability) gives the exact answer — the responsibility of component $j$ for point $i$:

Soft labels, with calibrated uncertainty at the boundaries — and full covariances $\Sigma_j$, so clusters can be elongated, tilted, and different sizes. Everything k-means couldn’t say.

EM: maximum likelihood when a variable is missing

Direct MLE on $\log p(x) = \log \sum_j (\cdots)$ is awkward — the sum sits inside the log. EM’s move: lower-bound the log-likelihood using Jensen’s inequality with any distribution $q$ over the latent labels,

then alternate: the E-step sets $q(z) = p(z \mid x, \theta)$ — the responsibilities — which makes the bound tight; the M-step maximizes the bound over $\theta$, which has closed form because the log now sits inside the expectation:

with $N_j = \sum_i r_{ij}$ the “effective count”. These are exactly the ordinary Gaussian MLEs from the Probability module, with each point counted fractionally by responsibility. Tight bound + maximized bound ⟹ the log-likelihood never decreases — EM’s celebrated monotonicity, and the number the widget audits every step:

k-means is a frozen mixture

Fix every $\Sigma_j = \sigma^2 I$ and take $\sigma \to 0$: responsibilities collapse to one-hot (nearest center wins), the M-step becomes the plain mean, and EM becomes Lloyd’s algorithm — k-means is the zero-temperature limit of a spherical Gaussian mixture. This is why k-means implicitly assumes round, equal clusters: it is doing MLE in a model that asserts them. The pattern — “popular heuristic = probabilistic model + degenerate limit” — recurs across ML, and knowing the model tells you exactly which dial you froze.

The fine print, all of it load-bearing

• Singularities. The likelihood is unbounded: park a component on a single point and shrink its variance, and $\log p \to \infty$. Real implementations add a covariance floor or a prior (MAP-EM) — the widget adds a small ridge, the Bayesian lesson’s medicine again.
• Choosing k, again — but now likelihood-based tools exist: BIC ($-2\log\hat L + p\log n$) and its kin penalize the extra parameters each component costs. Better grounded than the elbow, still not oracle.
• Covariance budgets. Full $\Sigma_j$ costs $O(kd^2)$ parameters — in high dimensions you tie covariances, make them diagonal, or go spherical: a bias–variance dial identical in spirit to regularization.
• Identifiability. Components can be permuted freely (“label switching”), which is harmless for clustering, treacherous for interpreting “component 2” across runs.
• EM beyond mixtures. The same E/M scaffolding fits hidden Markov models (the Forecasting track’s state-space cousins), factor analysis, and — with $q$ restricted to a tractable family instead of made exact — becomes variational inference, the engine of VAEs in the Deep Learning track. The bound $\mathcal{L}(q,\theta)$ you met here is the ELBO you will meet there.

import numpy as np
from sklearn.mixture import GaussianMixture
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=400, centers=3, random_state=2)
X = X @ np.array([[1.4, 0.7], [0.0, 0.4]])          # tilt: break sphericity

gm = GaussianMixture(n_components=3, covariance_type="full", random_state=0).fit(X)
r = gm.predict_proba(X)                              # responsibilities
print(gm.lower_bound_)                               # avg log-lik at convergence
print(r.max(axis=1).min())                           # least-confident point
print(gm.bic(X), GaussianMixture(6, random_state=0).fit(X).bic(X))  # k by BIC

Exercises