k-means, Gaussian Mixtures, and EM

1. \(J\) is non-increasing in \(K\). Any clustering with \(K\) centroids can be reproduced with \(K+1\) by adding a redundant centroid at an existing location (which changes nothing) — so the \(K+1\) optimum is at least as good. In practice \(J\) strictly decreases, since an extra centroid can always split some cluster.

2. Because a criterion that always prefers larger \(K\) does not select anything. It is exactly the \(R^2\) problem from Lecture 02: a measure of fit that cannot get worse when you add capacity is useless for model selection. You need a penalty for complexity, or a held-out estimate.

3. \(J=0\). Each point is its own centroid, so every distance is zero. Perfect fit, zero information — the clustering analogue of the lookup table from Lecture 01.

4. Two defensible options.

• A penalised criterion. Fit a Gaussian mixture instead and use BIC \(= -2\ell + p\log n\), which explicitly charges for parameters and therefore has an interior minimum. This is principled but leans on the mixture model being roughly right.
• The downstream task. If the clusters feed something you can measure — a classifier, a recommendation, a segmentation someone will act on — tune \(K\) by that metric with proper validation. Same argument as choosing PCA's \(k\) in Lecture 08: an external criterion beats an internal one whenever you have access to it.

The elbow method and the silhouette score are the common alternatives. Both are reasonable heuristics; neither is a principle, and both reward the kind of structure the algorithm was built to find.

The M-step. Let \(g(\vm)=\sum_{i\in C_k}\|\vx_i-\vm\|^2\). Using Phase A Lecture 12's rules,

The Hessian is \(2|C_k|\mI\succ0\), so \(g\) is strictly convex and this is the unique global minimum (Phase A Lecture 7). ∎ Note it is the same computation as "the mean minimises squared deviation" (Phase A Lecture 13), lifted to vectors.

The pairwise identity. Write \(m=|C_k|\) and \(\vm_k=\frac1m\sum_i\vx_i\). Expand the right-hand side:

And the left-hand side:

Dividing the first display by \(2m\) gives exactly the second. ∎

Why it is useful — three reasons.

It removes the centroid. The objective depends only on distances between data points, never on any point in the ambient space. So the method extends immediately to settings where averaging is undefined — strings, graphs, permutations — as long as you can compute pairwise distances. That is kernel k-means, and it is why the kernel trick of Lecture 12 applies here too.

It gives an interpretation. Inertia is (up to the \(2|C_k|\) factor) the total pairwise dissimilarity within each cluster. "Minimise distance to the centre" and "minimise average distance between members" are the same instruction.

It explains the spherical bias. Because only Euclidean distances enter, the objective is invariant to rotation but not to anisotropic scaling. A cluster that is long and thin has large internal pairwise distances and is therefore penalised, which is precisely why k-means prefers to cut it in half. The bias is visible directly in the objective.

1. With \(\pi_1=\pi_2=\tfrac12\) the prefactors and weights cancel:

Expand the exponent: \((x-\mu_1)^2-(x-\mu_2)^2 = -2x(\mu_1-\mu_2)+\mu_1^2-\mu_2^2\), which is linear in \(x\). So

— a logistic sigmoid in \(x\), with slope \((\mu_1-\mu_2)/\sigma^2\) and its 0.5 crossing at the midpoint \(\frac{\mu_1+\mu_2}{2}\). ∎

This is the generative-versus-discriminative connection of Lecture 10, appearing early: a Gaussian mixture with shared variance induces exactly a logistic posterior, which is why Gaussian naive Bayes and logistic regression are so closely related.

2. The slope is \((\mu_1-\mu_2)/\sigma^2\).

\(\sigma\to0\): the slope diverges, the sigmoid becomes a step function at the midpoint, and every point is assigned with probability 1 to its nearer mean. Hard assignment — the k-means E-step.

\(\sigma\to\infty\): the slope goes to zero and \(\gamma_{i1}\to\tfrac12\) for every point. The components are so diffuse that the data carries no information about which generated what, and the posterior reverts to the prior. EM then makes no progress: both means converge to the global mean.

3. Write \(\gamma_{ik}=\operatorname{softmax}_k(-d_{ik}^2/T)\) with \(T=2\sigma^2\). Then \(\sigma^2\) is a temperature in the statistical-mechanics sense. High temperature gives a soft, near-uniform distribution; low temperature concentrates all mass on the minimum. \(T\to0\) turns softmax into argmax — the same limit that appears in Gumbel-softmax, in simulated annealing, and in the temperature parameter of a language model's sampler.

So the k-means/GMM relationship is one instance of a general pattern: a hard combinatorial algorithm is often the zero-temperature limit of a soft probabilistic one, and the soft version is differentiable while the hard one is not. That observation is the engine behind a great deal of modern method design.

import numpy as np
rng = np.random.default_rng(0)

def kmeans(X, K, seed=0, T=100):
    r = np.random.default_rng(seed)
    M = X[r.choice(len(X), K, replace=False)].copy()
    for _ in range(T):
        D = ((X[:,None,:]-M[None,:,:])**2).sum(-1); a = D.argmin(1)
        newM = np.array([X[a==k].mean(0) if (a==k).any() else M[k] for k in range(K)])
        if np.allclose(newM, M): break
        M = newM
    D = ((X[:,None,:]-M[None,:,:])**2).sum(-1); a = D.argmin(1)
    return float(D[np.arange(len(X)), a].sum()), a

def silhouette(X, a):
    D = np.sqrt(((X[:,None,:]-X[None,:,:])**2).sum(-1)); s = []
    for i in range(len(X)):
        own = a == a[i]; own[i] = False
        if own.sum() == 0: s.append(0.0); continue
        ai = D[i, own].mean()
        bi = min(D[i, a==k].mean() for k in set(a) - {a[i]})
        s.append((bi-ai)/max(ai, bi))
    return float(np.mean(s))

def gap(X, Kmax=10, B=20):
    lo, hi = X.min(0), X.max(0); out = []
    for K in range(1, Kmax+1):
        Jd = min(kmeans(X, K, seed=s)[0] for s in range(5))
        ref = [np.log(min(kmeans(rng.uniform(lo, hi, X.shape), K, seed=s)[0]
                          for s in range(3))) for _ in range(B)]
        out.append((K, np.log(Jd), np.mean(ref), np.mean(ref)-np.log(Jd), np.std(ref)))
    return out

def clustered(n=500):
    c = np.array([[0,0],[4,0],[2,3.5]])
    return np.vstack([c[k] + 0.5*rng.normal(size=(n//3+1,2)) for k in range(3)])[:n]

for name, X in (("uniform noise", rng.uniform(0,1,(500,2))), ("three clusters", clustered())):
    print(f"\\n=== {name} ===")
    print(f"{'K':>3} {'inertia':>10} {'gap':>9} {'sd_ref':>8} {'silhouette':>11}")
    for K, lJ, lref, g, sd in gap(X):
        sil = silhouette(X, kmeans(X, K, seed=0)[1]) if K > 1 else float('nan')
        print(f"{K:3d} {np.exp(lJ):10.3f} {g:9.4f} {sd:8.4f} {sil:11.4f}")

1. On uniform noise, inertia falls smoothly and steadily — roughly like \(K^{-1}\) in two dimensions, since \(K\) cells each shrink in area. There is no elbow, but the curve looks perfectly plausible, and a determined eye will find a bend somewhere around \(K=3\) or \(4\). This is the central weakness of the elbow method: it produces an answer whether or not there is one.

2. The gap statistic compares \(\log J_K\) against what you would get on structureless data of the same extent. On uniform noise the gap stays near zero at every \(K\) (well within one reference standard deviation) — correctly reporting no structure, and the only diagnostic here that does. Tibshirani's rule is to choose the smallest \(K\) with \(\text{gap}(K)\ge\text{gap}(K+1)-s_{K+1}\); on noise no \(K\) qualifies beyond 1.

3. On the three-cluster data the inertia curve has a sharp, unmistakable elbow at \(K=3\); the gap statistic peaks at \(K=3\) and is several reference standard deviations above zero; and the silhouette score peaks at \(K=3\) at around 0.7–0.8. All three agree, which is what agreement is for.

4. Which to trust, and the honest conclusion. The silhouette on uniform noise typically comes out around 0.3–0.4 — a value many people would read as "moderate structure". It is not zero, because silhouette measures whether the found partition is internally coherent, and \(k\)-means will always produce a reasonably coherent partition of anything. Silhouette can validate a clustering; it cannot tell you whether clusters exist.

The gap statistic is the only one of the three with a built-in null hypothesis, and that is exactly why it can detect the absence of structure. The general principle is worth taking away: a diagnostic without a reference distribution cannot distinguish signal from its absence — which is the same lesson as "compare against the majority-class baseline" in Lecture 07, and "compare against random subspaces" in Lecture 08.

Practical advice: always run your clustering pipeline on shuffled or resampled data as a null, and report what it produces alongside your real result.