{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lesson 19: Clustering — k-Means, GMMs, and DBSCAN\n",
"\n",
"Every prior lesson in this part assumed some *target* — a threshold, a template, a known transform. **Clustering** asks a different kind of question, with no target at all: given a pile of points (or pixels), group them into a small number of sensible categories, using only how similar they are to each other. This is **unsupervised learning**'s classical starting point — no labels, just a distance function and an assumption about what \"a cluster\" should look like. This lesson builds three clustering algorithms from scratch, each built on a different assumption about cluster shape, and shows concretely where each one wins and fails: **k-means** (round, similarly-sized clusters), **Gaussian mixture models** (elliptical clusters, soft assignment), and **DBSCAN** (arbitrary shape, defined by density rather than a center)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c23ca5c1",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import cv2\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"id": "3f61f144",
"metadata": {},
"source": "## k-means on real color data\n\nLesson 18 built a perceptually meaningful color distance (L\\*a\\*b\\* Euclidean distance). Clustering is the natural next step: given a pile of pixel colors and no labels at all, group them into regions automatically. Build a synthetic 4-region color image (with noise, so it's not trivially separable by exact color match), convert to L*a*b*, and cluster the raw pixel colors."
},
{
"cell_type": "code",
"execution_count": null,
"id": "e28b4f52",
"metadata": {},
"outputs": [],
"source": [
"SIZE = 48\n",
"\n",
"def make_region_image(size=SIZE, noise_std=8, seed=3):\n",
" rng = np.random.default_rng(seed)\n",
" img = np.zeros((size, size, 3), dtype=np.uint8)\n",
" region_map = np.zeros((size, size), dtype=np.int64)\n",
" colors_bgr = [(180, 40, 30), (40, 160, 60), (30, 80, 190), (200, 200, 60)]\n",
" centers = [(12, 12), (12, 36), (36, 12), (36, 36)]\n",
" yy, xx = np.mgrid[0:size, 0:size]\n",
" for k, ((cy, cx), col) in enumerate(zip(centers, colors_bgr)):\n",
" mask = (xx - cx) ** 2 + (yy - cy) ** 2 <= 20 ** 2\n",
" img[mask] = col\n",
" region_map[mask] = k + 1\n",
" noise = rng.normal(0, noise_std, img.shape)\n",
" img = np.clip(img.astype(np.float64) + noise, 0, 255).astype(np.uint8)\n",
" return img, region_map\n",
"\n",
"img, region_map = make_region_image()\n",
"lab = cv2.cvtColor(img, cv2.COLOR_BGR2LAB).astype(np.float64)\n",
"X_pixels = lab.reshape(-1, 3)\n",
"\n",
"def kmeans(X, k, n_iter=50, seed=0):\n",
" rng = np.random.default_rng(seed)\n",
" centers = X[rng.choice(len(X), k, replace=False)].copy()\n",
" for _ in range(n_iter):\n",
" dists = ((X[:, None, :] - centers[None, :, :]) ** 2).sum(-1)\n",
" labels = dists.argmin(1)\n",
" new_centers = np.array([X[labels == j].mean(0) if (labels == j).any() else centers[j]\n",
" for j in range(k)])\n",
" if np.allclose(new_centers, centers):\n",
" break\n",
" centers = new_centers\n",
" return labels, centers\n",
"\n",
"def purity(labels, true_flat):\n",
" correct = 0\n",
" for u in np.unique(labels):\n",
" mask = labels == u\n",
" majority = np.bincount(true_flat[mask]).argmax()\n",
" correct += (true_flat[mask] == majority).sum()\n",
" return correct / len(true_flat)\n",
"\n",
"km_labels, _ = kmeans(X_pixels, k=4, seed=3)\n",
"km_seg = km_labels.reshape(SIZE, SIZE)\n",
"print(f'k-means segmentation purity vs. true regions: {purity(km_labels, region_map.ravel()):.3f}')\n",
"\n",
"fig, axes = plt.subplots(1, 3, figsize=(9, 3.2))\n",
"axes[0].imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB)); axes[0].set_title('input (noisy)', fontsize=9); axes[0].axis('off')\n",
"axes[1].imshow(region_map, cmap='tab10'); axes[1].set_title('true regions', fontsize=9); axes[1].axis('off')\n",
"axes[2].imshow(km_seg, cmap='tab10'); axes[2].set_title('k-means segmentation', fontsize=9); axes[2].axis('off')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "fe03e92f",
"metadata": {},
"source": "K-means, also known as Lloyd's algorithm, is here unchanged from its textbook form: pick `k` random points as initial centers, assign every point to its nearest center, recompute each center as the mean of its assigned points, repeat until nothing moves. With compact, roughly equal-sized, roughly round clusters in L*a*b* space (exactly what this synthetic image has), it recovers the true regions almost perfectly. Two real caveats worth flagging even in this success case: `k` (the number of clusters) has to be chosen in advance — nothing in the algorithm discovers it — and different random initializations can converge to different, sometimes noticeably worse, local optima (try changing `seed=3` to `seed=0` above and rerun); **k-means++** (Arthur & Vassilvitskii, 2007★) fixes this by seeding centers to be spread apart rather than purely at random."
},
{
"cell_type": "markdown",
"id": "ca8a3a20",
"metadata": {},
"source": "## Where k-means breaks: non-convex clusters\n\nK-means' update rule — assign to the nearest *center* — can only carve space into convex regions. Build a dataset that isn't convex at all: two concentric rings. No placement of two centers can separate them correctly, no matter how many times the algorithm iterates."
},
{
"cell_type": "code",
"execution_count": null,
"id": "074ce5cd",
"metadata": {},
"outputs": [],
"source": [
"rng = np.random.default_rng(1)\n",
"\n",
"def ring(n, r, noise=0.05):\n",
" theta = rng.uniform(0, 2 * np.pi, n)\n",
" rad = r + rng.normal(0, noise, n)\n",
" return np.stack([rad * np.cos(theta), rad * np.sin(theta)], axis=1)\n",
"\n",
"X_rings = np.vstack([ring(150, 1.0), ring(150, 2.5)])\n",
"true_rings = np.array([0] * 150 + [1] * 150)\n",
"\n",
"def best_binary_acc(labels, true):\n",
" return max((labels == true).mean(), (labels == (1 - true)).mean())\n",
"\n",
"km_ring_labels, km_ring_centers = kmeans(X_rings, k=2, seed=0)\n",
"print(f'k-means accuracy on concentric rings: {best_binary_acc(km_ring_labels, true_rings):.1%} (2-class chance = 50%)')\n",
"\n",
"fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))\n",
"axes[0].scatter(*X_rings.T, c=true_rings, cmap='coolwarm', s=8)\n",
"axes[0].set_title('true rings', fontsize=9); axes[0].set_aspect('equal')\n",
"axes[1].scatter(*X_rings.T, c=km_ring_labels, cmap='coolwarm', s=8)\n",
"axes[1].scatter(*km_ring_centers.T, c='black', marker='x', s=100, label='centers')\n",
"axes[1].set_title('k-means result', fontsize=9); axes[1].set_aspect('equal'); axes[1].legend(fontsize=7)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "0c4646e1",
"metadata": {},
"source": [
"K-means lands right at chance — it splits the two rings radially, half of each ring going to each cluster, because that's the only kind of boundary (a straight line, equidistant between two centers) the algorithm is capable of drawing. The clusters here aren't ambiguous to a human eye at all; the *algorithm's assumption* (round, center-based clusters) is simply the wrong tool for this shape."
]
},
{
"cell_type": "markdown",
"id": "f580d6d3",
"metadata": {},
"source": "## Gaussian mixture models: soft, elliptical clusters\n\nA **Gaussian mixture model (GMM)** generalizes k-means two ways: each cluster gets its own full covariance (so clusters can be elongated ellipses, not just circles), and assignment is *soft* — every point gets a probability of belonging to each cluster, not a hard label. Fit by **expectation-maximization (EM)** (Dempster, Laird, & Rubin, 1977★): the E-step computes each point's responsibility (probability) under the current Gaussians; the M-step re-fits each Gaussian's mean, covariance, and weight to its responsibility-weighted points."
},
{
"cell_type": "code",
"execution_count": null,
"id": "f8bef68e",
"metadata": {},
"outputs": [],
"source": [
"def gmm_em(X, k, n_iter=100, seed=0):\n",
" rng = np.random.default_rng(seed)\n",
" n, d = X.shape\n",
" means = X[rng.choice(n, k, replace=False)].copy()\n",
" covs = np.array([np.cov(X.T) + 1e-3 * np.eye(d) for _ in range(k)])\n",
" weights = np.full(k, 1 / k)\n",
" for _ in range(n_iter):\n",
" resp = np.zeros((n, k))\n",
" for j in range(k):\n",
" diff = X - means[j]\n",
" inv = np.linalg.inv(covs[j])\n",
" expo = -0.5 * np.sum(diff @ inv * diff, axis=1)\n",
" norm = 1.0 / np.sqrt((2 * np.pi) ** d * np.linalg.det(covs[j]))\n",
" resp[:, j] = weights[j] * norm * np.exp(expo)\n",
" resp /= resp.sum(1, keepdims=True) + 1e-12\n",
" Nk = resp.sum(0)\n",
" weights = Nk / n\n",
" means = (resp.T @ X) / Nk[:, None]\n",
" for j in range(k):\n",
" diff = X - means[j]\n",
" covs[j] = (resp[:, j:j + 1] * diff).T @ diff / Nk[j] + 1e-6 * np.eye(d)\n",
" return resp.argmax(1), means, covs\n",
"\n",
"gmm_ring_labels, gmm_means, gmm_covs = gmm_em(X_rings, k=2, seed=0)\n",
"print(f'GMM accuracy on concentric rings: {best_binary_acc(gmm_ring_labels, true_rings):.1%}')"
]
},
{
"cell_type": "markdown",
"id": "cfeb07c0",
"metadata": {},
"source": [
"GMM does no better than k-means here — also near chance. That's an important, easy-to-miss point: GMM fixes k-means' *round-only* assumption, but it's still a **unimodal-per-cluster** model — each cluster is described by one Gaussian bump. A ring isn't elliptical any more than it's circular; no single Gaussian, of any shape, fits an annulus well. GMM and k-means fail this dataset for the same underlying reason."
]
},
{
"cell_type": "markdown",
"id": "f8687440",
"metadata": {},
"source": "## DBSCAN: clusters defined by density, not shape\n\n**DBSCAN** (density-based spatial clustering, Ester et al., 1996★) never fits a parametric shape at all. A point is a **core point** if at least `min_samples` other points lie within distance `eps` of it. Clusters are formed by chaining together core points that are within `eps` of each other (and their neighbors), so a cluster can be any shape — including a ring — as long as it's a *connected, sufficiently dense* region of space. Points that end up in no core point's neighborhood are labeled **noise**, not forced into the nearest cluster."
},
{
"cell_type": "code",
"execution_count": null,
"id": "9d19677c",
"metadata": {},
"outputs": [],
"source": [
"def dbscan(X, eps, min_samples):\n",
" n = len(X)\n",
" labels = np.full(n, -1)\n",
" visited = np.zeros(n, dtype=bool)\n",
" dist = np.sqrt(((X[:, None, :] - X[None, :, :]) ** 2).sum(-1))\n",
" cluster_id = 0\n",
" for i in range(n):\n",
" if visited[i]:\n",
" continue\n",
" visited[i] = True\n",
" neighbors = list(np.where(dist[i] <= eps)[0])\n",
" if len(neighbors) < min_samples:\n",
" continue # stays noise (-1)\n",
" labels[i] = cluster_id\n",
" j = 0\n",
" while j < len(neighbors):\n",
" q = neighbors[j]\n",
" if not visited[q]:\n",
" visited[q] = True\n",
" q_neighbors = np.where(dist[q] <= eps)[0]\n",
" if len(q_neighbors) >= min_samples:\n",
" neighbors.extend([x for x in q_neighbors if x not in neighbors])\n",
" if labels[q] == -1:\n",
" labels[q] = cluster_id\n",
" j += 1\n",
" cluster_id += 1\n",
" return labels\n",
"\n",
"db_ring_labels = dbscan(X_rings, eps=0.6, min_samples=4)\n",
"n_clusters = len(set(db_ring_labels.tolist()) - {-1})\n",
"n_noise = (db_ring_labels == -1).sum()\n",
"print(f'DBSCAN found {n_clusters} clusters, {n_noise} noise points')\n",
"print(f'DBSCAN accuracy on concentric rings: {best_binary_acc(db_ring_labels, true_rings):.1%}')\n",
"\n",
"fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))\n",
"for ax, labels, title in zip(axes, [true_rings, km_ring_labels, db_ring_labels],\n",
" ['true rings', 'k-means (fails)', 'DBSCAN (succeeds)']):\n",
" ax.scatter(*X_rings.T, c=labels, cmap='coolwarm', s=8)\n",
" ax.set_title(title, fontsize=9); ax.set_aspect('equal')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "38457105",
"metadata": {},
"source": [
"DBSCAN separates the two rings perfectly — it never assumed a center or a shape, only that points *within* a ring are densely connected to their neighbors, while the gap between rings is not. This comes at a real cost: `eps` and `min_samples` have to be chosen by hand (a poor choice fragments each ring into dozens of tiny arcs — worth trying `eps=0.15` above to see it happen), and DBSCAN struggles when different true clusters have very different densities, since one global `eps` can't be simultaneously right for a sparse cluster and a dense one."
]
},
{
"cell_type": "markdown",
"id": "22ef3473",
"metadata": {},
"source": [
"## Noise robustness: a real advantage of \"no forced assignment\"\n",
"\n",
"Add clutter — points scattered uniformly across the whole region, belonging to neither ring — and compare how each algorithm handles them. k-means and GMM have no concept of \"doesn't belong to any cluster\": every point gets assigned somewhere, dragging cluster centers toward the clutter. DBSCAN can simply call clutter what it is."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "adae561f",
"metadata": {},
"outputs": [],
"source": [
"noise_pts = rng.uniform(-3, 3, (60, 2))\n",
"X_cluttered = np.vstack([X_rings, noise_pts])\n",
"true_cluttered = np.concatenate([true_rings, np.full(60, -1)]) # -1 = \"not really a ring\"\n",
"\n",
"km_clut_labels, _ = kmeans(X_cluttered, k=2, seed=0)\n",
"db_clut_labels = dbscan(X_cluttered, eps=0.6, min_samples=4)\n",
"\n",
"km_noise_correctly_flagged = 0 # k-means has no noise concept at all\n",
"db_noise_correctly_flagged = (db_clut_labels[300:] == -1).sum()\n",
"\n",
"print(f'{\"\":>10} {\"noise points correctly flagged\":>32} {\"out of\":>8}')\n",
"print(f'{\"k-means\":>10} {km_noise_correctly_flagged:>32} {60:>8} (no noise concept — every point forced into a cluster)')\n",
"print(f'{\"DBSCAN\":>10} {db_noise_correctly_flagged:>32} {60:>8}')"
]
},
{
"cell_type": "markdown",
"id": "dcbfdd7d",
"metadata": {},
"source": [
"k-means, by construction, flags exactly zero of the 60 clutter points as anything other than \"definitely part of a ring\" — the algorithm has no vocabulary for \"doesn't belong anywhere.\" DBSCAN correctly isolates a minority of them; the rest happen to land within `eps` of a real ring point (or of each other, forming a small cluster of their own), which is a fair outcome, not a bug — density-based clustering only calls a point noise if it's genuinely isolated, and uniform random scatter still produces occasional clumps. The qualitative point stands regardless of the exact count: DBSCAN has a real, usable category for \"this doesn't fit anywhere,\" whereas k-means structurally cannot."
]
},
{
"cell_type": "markdown",
"id": "2428ca35",
"metadata": {},
"source": "## Where this goes next\n\nThis is the end of Part 1. Many of the algorithms encountered so far operate directly on **raw or hand-picked features** — pixel color, `(x, y)` coordinates, texture — and rely on specific rules, such as the human-chosen distance functions used here. Future lessons revisit many of these problems using deep learning to automatically discover the features and the rules for combining them. For example, Lesson 43's instance segmentation uses a similar clustering idea, and Lesson 48's autoencoder learns the feature space itself, purely from reconstruction error, before any grouping is asked for at all. Where this lesson's methods assume a shape and fit it, the deep-learning methods later in this course learn a representation and let structure emerge from it.\n\n### Exercises\n\n1. Rerun the color-segmentation k-means with `k=3` and `k=6` instead of the true `4`. Since nothing in the algorithm can \"know\" the right number of clusters, how does the segmentation degrade in each direction, and can you tell purely from the output which `k` is likely wrong?\n2. The GMM above was evaluated only on the ring dataset, where it fails for the same reason k-means does. Construct a two-cluster dataset where the clusters are genuinely elliptical (e.g. `np.random.multivariate_normal` with a non-diagonal covariance) and compare GMM against k-means there. Does GMM's extra flexibility (fitting the shape, not just the center) show a real accuracy advantage this time?\n3. In the noisy-clutter DBSCAN demo, shrink `eps` from `0.6` to `0.3`. Does DBSCAN start mislabeling real ring points as noise, real clutter points as belonging to a ring, or both — and does that match the earlier warning about `eps` being a hand-tuned, dataset-specific choice rather than something the algorithm discovers?"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.10.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}