Unsupervised Learning — Clustering

Clustering in NumPy from scratch: K-means, Gaussian Mixture Models, the Expectation-Maximization algorithm, BIC model selection, and the sklearn/scipy APIs that replace them.

---

Learning Objectives

---

Task-by-Task Reference

Each task entry captures only what is new relative to all previous tasks.

---

Task 0 — Initialize K-means (0-initialize.py)

Challenge: Initialize (k, d) centroids without loops, with exactly one numpy.random.uniform call, each dimension bounded by the data's range.

Approach: Take X.min(axis=0) and X.max(axis=0) (per-dimension ranges, shape (d,)), then let NumPy broadcast them into a single np.random.uniform(mins, maxs, size=(k, d)) draw.

New techniques introduced:

Key takeaway: Broadcasting turns a k-loop of per-dimension random draws into a single array expression.

---

Task 1 — K-means (1-kmeans.py)

Challenge: Assign every point to its nearest centroid without per-point loops; re-seed empty clusters; stop when centroids stop moving.

Approach: Compute all pairwise distances with broadcasting (X[:, None, :] - C), take argmin(axis=1) for assignments, and recompute each centroid as the boolean-mask mean of its points. A cluster with no points gets re-seeded with np.random.uniform(mins, maxs). Loop breaks when np.all(C == C_old).

New techniques introduced:

Key takeaway: The assignment step is a vectorized nearest-neighbour query; the loop budget goes to iterations and clusters, never to points.

---

Task 2 — Variance (2-variance.py)

Challenge: Total intra-cluster variance with no loops.

Approach: Chain reductions in one expression: (diffs ** 2).sum(axis=2).min(axis=1).sum() — squared distance to every centroid, keep the smallest per point, sum over points.

New techniques introduced:

Key takeaway: Variance measures how tightly points hug their assigned centroid — the score that K-means minimizes.

---

Task 3 — Optimize K (3-optimum.py)

Challenge: Find the "right" k by measuring how variance drops as clusters grow; kmax=None defaults to the number of points; two loops max.

Approach: Sweep k from kmin to kmax, run K-means and variance for each, then express every variance as the drop from the kmin baseline (d_vars). The elbow — where extra clusters stop buying much variance — is the answer.

New techniques introduced:

Key takeaway: More clusters always reduce variance; the elbow is where the reduction stops being worth it.

---

Task 4 — Initialize GMM (4-initialize.py)

Challenge: Seed a Gaussian Mixture Model — priors, means, covariances — with no loops, reusing K-means for the means.

Approach: pi = np.full(k, 1/k) (neutral priors), m from kmeans(X, k)[0], and covariances as identity matrices stacked to (k, d, d) via np.identity(d)[np.newaxis].repeat(k, axis=0).

New techniques introduced:

Key takeaway: GMM initialization is K-means for the means plus neutral priors and identity covariances — EM does the rest.

---

Task 5 — Multivariate Gaussian PDF (5-pdf.py)

Challenge: Vectorize the quadratic form (x−μ)ᵀΣ⁻¹(x−μ) for all points and avoid numeric underflow in the exponent.

Approach: Compute the normalization constant from (2π)^(−d/2) · det(S)^(−1/2), evaluate the quadratic form with the matmul-multiply trick ((X − m) @ np.linalg.inv(S)) * (X − m) summed over axis=1, and floor the result with np.maximum(P, 1e-300).

New techniques introduced:

Key takeaway: The quadratic form is a row-wise dot product — express it as a matmul followed by an element-wise multiply to avoid a loop over points.

---

Task 6 — Expectation step (6-expectation.py)

Challenge: Compute soft assignments g[k, n] = π_k · pdf_k / Σ_j π_j · pdf_j and the model log-likelihood, with one loop.

Approach: For each component, g[i] = pi[i] * pdf(X, m[i], S[i]); accumulate the mixture evidence in denominator; normalize; then ll = np.sum(np.log(denominator)).

New techniques introduced:

Key takeaway: ll scores how well the mixture explains the data (Σ log P(x)); normalized g would give log 1 = 0, which is why the log is taken on the denominator.

---

Task 7 — Maximization step (7-maximization.py)

Challenge: Update (π, μ, Σ) from the soft assignments — the M-step of EM.

Approach: For each component: nk = g[i].sum() (soft count N_k), pi[i] = nk / n, m[i] = g[i] @ X / nk, and the weighted covariance (g[i][:, None] * diff).T @ diff / nk with diff = X − m[i].

New techniques introduced:

Key takeaway: Every M-step formula is the familiar statistic with hard assignments replaced by soft responsibilities.

---

Task 8 — Expectation Maximization (8-EM.py)

Challenge: Orchestrate E and M steps until the log-likelihood stops moving (|Δll| ≤ tol) or the budget runs out, with a verbose reporting protocol.

Approach: Initial E step, then a single loop: print every 10 iterations, run M, run E, check abs(ll − prev_ll) <= tol and break. The final reported g/ll come from the E step run after the last M step, so the return is self-consistent.

New techniques introduced:

Key takeaway: EM climbs the likelihood hill; when the altimeter stops moving by more than tol, you have arrived.

---

Task 9 — Bayesian Information Criterion (9-BIC.py)

Challenge: Pick the best k by balancing fit against complexity: more clusters always raise ll, so charge a parameter tax.

Approach: Sweep k, run EM for each, count free parameters p = (k−1) + k·d + k·d(d+1)/2 (priors minus one, means, symmetric covariances), score b = p·ln(n) − 2·ll, and take argmin.

New techniques introduced:

Key takeaway: The tax p·ln(n) grows with k, so the best k is where the fit gain stops covering the complexity cost.

---

Task 10 — K-means, sklearn (10-kmeans.py)

Challenge: Replace the hand-rolled K-means with scikit-learn's estimator.

Approach: sklearn.cluster.KMeans(k).fit(X), then read cluster_centers_ and labels_ — the same C and clss as before, via fit-then-read.

New techniques introduced:

Key takeaway: sklearn exposes results as fitted attributes, not return values — fit first, read afterwards.

---

Task 11 — GMM, sklearn (11-gmm.py)

Challenge: Same pipeline for a Gaussian Mixture Model, plus a BIC score straight from the library.

Approach: sklearn.mixture.GaussianMixture(k).fit(X), read weights_, means_, covariances_, assign with predict(X), score with bic(X).

New techniques introduced:

Key takeaway: bic is a bound method: without (X) you get the function, not the score.

---

Task 12 — Agglomerative (12-agglomerative.py)

Challenge: Hierarchical clustering with imports restricted to scipy.cluster.hierarchy and matplotlib.pyplot — no sklearn, no own algorithm.

Approach: Build the tree with linkage(X, method='ward'), cut it at the cophenetic threshold with fcluster(Z, dist, criterion='distance'), and draw it with dendrogram(Z, color_threshold=dist). Cluster indices are 1-based by scipy convention.

New techniques introduced:

Key takeaway: The whole hierarchy is one linkage call; fcluster is simply "cut the tree at this height".

---

Technique Inventory

---

Resources

• scipy.cluster.hierarchy — Hierarchical clustering
• sklearn.mixture.GaussianMixture
• sklearn.cluster.KMeans
• NumPy Array Broadcasting
• Bayesian Information Criterion (Wikipedia)
• Multivariate normal distribution (Wikipedia)