Lesson 19: Clustering — k-Means, GMMs, and DBSCAN

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).

In [1]:
import numpy as np
import cv2
import matplotlib.pyplot as plt

k-means on real color data

Lesson 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 Lab*, and cluster the raw pixel colors.

In [2]:
SIZE = 48

def make_region_image(size=SIZE, noise_std=8, seed=3):
    rng = np.random.default_rng(seed)
    img = np.zeros((size, size, 3), dtype=np.uint8)
    region_map = np.zeros((size, size), dtype=np.int64)
    colors_bgr = [(180, 40, 30), (40, 160, 60), (30, 80, 190), (200, 200, 60)]
    centers = [(12, 12), (12, 36), (36, 12), (36, 36)]
    yy, xx = np.mgrid[0:size, 0:size]
    for k, ((cy, cx), col) in enumerate(zip(centers, colors_bgr)):
        mask = (xx - cx) ** 2 + (yy - cy) ** 2 <= 20 ** 2
        img[mask] = col
        region_map[mask] = k + 1
    noise = rng.normal(0, noise_std, img.shape)
    img = np.clip(img.astype(np.float64) + noise, 0, 255).astype(np.uint8)
    return img, region_map

img, region_map = make_region_image()
lab = cv2.cvtColor(img, cv2.COLOR_BGR2LAB).astype(np.float64)
X_pixels = lab.reshape(-1, 3)

def kmeans(X, k, n_iter=50, seed=0):
    rng = np.random.default_rng(seed)
    centers = X[rng.choice(len(X), k, replace=False)].copy()
    for _ in range(n_iter):
        dists = ((X[:, None, :] - centers[None, :, :]) ** 2).sum(-1)
        labels = dists.argmin(1)
        new_centers = np.array([X[labels == j].mean(0) if (labels == j).any() else centers[j]
                                 for j in range(k)])
        if np.allclose(new_centers, centers):
            break
        centers = new_centers
    return labels, centers

def purity(labels, true_flat):
    correct = 0
    for u in np.unique(labels):
        mask = labels == u
        majority = np.bincount(true_flat[mask]).argmax()
        correct += (true_flat[mask] == majority).sum()
    return correct / len(true_flat)

km_labels, _ = kmeans(X_pixels, k=4, seed=3)
km_seg = km_labels.reshape(SIZE, SIZE)
print(f'k-means segmentation purity vs. true regions: {purity(km_labels, region_map.ravel()):.3f}')

fig, axes = plt.subplots(1, 3, figsize=(9, 3.2))
axes[0].imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB)); axes[0].set_title('input (noisy)', fontsize=9); axes[0].axis('off')
axes[1].imshow(region_map, cmap='tab10'); axes[1].set_title('true regions', fontsize=9); axes[1].axis('off')
axes[2].imshow(km_seg, cmap='tab10'); axes[2].set_title('k-means segmentation', fontsize=9); axes[2].axis('off')
plt.show()
k-means segmentation purity vs. true regions: 1.000
No description has been provided for this image

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 Lab* 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.

Where k-means breaks: non-convex clusters

K-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.

In [3]:
rng = np.random.default_rng(1)

def ring(n, r, noise=0.05):
    theta = rng.uniform(0, 2 * np.pi, n)
    rad = r + rng.normal(0, noise, n)
    return np.stack([rad * np.cos(theta), rad * np.sin(theta)], axis=1)

X_rings = np.vstack([ring(150, 1.0), ring(150, 2.5)])
true_rings = np.array([0] * 150 + [1] * 150)

def best_binary_acc(labels, true):
    return max((labels == true).mean(), (labels == (1 - true)).mean())

km_ring_labels, km_ring_centers = kmeans(X_rings, k=2, seed=0)
print(f'k-means accuracy on concentric rings: {best_binary_acc(km_ring_labels, true_rings):.1%}  (2-class chance = 50%)')

fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
axes[0].scatter(*X_rings.T, c=true_rings, cmap='coolwarm', s=8)
axes[0].set_title('true rings', fontsize=9); axes[0].set_aspect('equal')
axes[1].scatter(*X_rings.T, c=km_ring_labels, cmap='coolwarm', s=8)
axes[1].scatter(*km_ring_centers.T, c='black', marker='x', s=100, label='centers')
axes[1].set_title('k-means result', fontsize=9); axes[1].set_aspect('equal'); axes[1].legend(fontsize=7)
plt.show()
k-means accuracy on concentric rings: 54.0%  (2-class chance = 50%)
No description has been provided for this image

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.

Gaussian mixture models: soft, elliptical clusters

A 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.

In [4]:
def gmm_em(X, k, n_iter=100, seed=0):
    rng = np.random.default_rng(seed)
    n, d = X.shape
    means = X[rng.choice(n, k, replace=False)].copy()
    covs = np.array([np.cov(X.T) + 1e-3 * np.eye(d) for _ in range(k)])
    weights = np.full(k, 1 / k)
    for _ in range(n_iter):
        resp = np.zeros((n, k))
        for j in range(k):
            diff = X - means[j]
            inv = np.linalg.inv(covs[j])
            expo = -0.5 * np.sum(diff @ inv * diff, axis=1)
            norm = 1.0 / np.sqrt((2 * np.pi) ** d * np.linalg.det(covs[j]))
            resp[:, j] = weights[j] * norm * np.exp(expo)
        resp /= resp.sum(1, keepdims=True) + 1e-12
        Nk = resp.sum(0)
        weights = Nk / n
        means = (resp.T @ X) / Nk[:, None]
        for j in range(k):
            diff = X - means[j]
            covs[j] = (resp[:, j:j + 1] * diff).T @ diff / Nk[j] + 1e-6 * np.eye(d)
    return resp.argmax(1), means, covs

gmm_ring_labels, gmm_means, gmm_covs = gmm_em(X_rings, k=2, seed=0)
print(f'GMM accuracy on concentric rings: {best_binary_acc(gmm_ring_labels, true_rings):.1%}')
GMM accuracy on concentric rings: 55.7%

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.

DBSCAN: clusters defined by density, not shape

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.

In [5]:
def dbscan(X, eps, min_samples):
    n = len(X)
    labels = np.full(n, -1)
    visited = np.zeros(n, dtype=bool)
    dist = np.sqrt(((X[:, None, :] - X[None, :, :]) ** 2).sum(-1))
    cluster_id = 0
    for i in range(n):
        if visited[i]:
            continue
        visited[i] = True
        neighbors = list(np.where(dist[i] <= eps)[0])
        if len(neighbors) < min_samples:
            continue  # stays noise (-1)
        labels[i] = cluster_id
        j = 0
        while j < len(neighbors):
            q = neighbors[j]
            if not visited[q]:
                visited[q] = True
                q_neighbors = np.where(dist[q] <= eps)[0]
                if len(q_neighbors) >= min_samples:
                    neighbors.extend([x for x in q_neighbors if x not in neighbors])
            if labels[q] == -1:
                labels[q] = cluster_id
            j += 1
        cluster_id += 1
    return labels

db_ring_labels = dbscan(X_rings, eps=0.6, min_samples=4)
n_clusters = len(set(db_ring_labels.tolist()) - {-1})
n_noise = (db_ring_labels == -1).sum()
print(f'DBSCAN found {n_clusters} clusters, {n_noise} noise points')
print(f'DBSCAN accuracy on concentric rings: {best_binary_acc(db_ring_labels, true_rings):.1%}')

fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))
for ax, labels, title in zip(axes, [true_rings, km_ring_labels, db_ring_labels],
                              ['true rings', 'k-means (fails)', 'DBSCAN (succeeds)']):
    ax.scatter(*X_rings.T, c=labels, cmap='coolwarm', s=8)
    ax.set_title(title, fontsize=9); ax.set_aspect('equal')
plt.show()
DBSCAN found 2 clusters, 0 noise points
DBSCAN accuracy on concentric rings: 100.0%
No description has been provided for this image

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.

Noise robustness: a real advantage of "no forced assignment"

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.

In [6]:
noise_pts = rng.uniform(-3, 3, (60, 2))
X_cluttered = np.vstack([X_rings, noise_pts])
true_cluttered = np.concatenate([true_rings, np.full(60, -1)])  # -1 = "not really a ring"

km_clut_labels, _ = kmeans(X_cluttered, k=2, seed=0)
db_clut_labels = dbscan(X_cluttered, eps=0.6, min_samples=4)

km_noise_correctly_flagged = 0  # k-means has no noise concept at all
db_noise_correctly_flagged = (db_clut_labels[300:] == -1).sum()

print(f'{"":>10} {"noise points correctly flagged":>32} {"out of":>8}')
print(f'{"k-means":>10} {km_noise_correctly_flagged:>32} {60:>8}  (no noise concept — every point forced into a cluster)')
print(f'{"DBSCAN":>10} {db_noise_correctly_flagged:>32} {60:>8}')
             noise points correctly flagged   out of
   k-means                                0       60  (no noise concept — every point forced into a cluster)
    DBSCAN                                8       60

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.

Where this goes next

This 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.

Exercises

  1. 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?
  2. 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?
  3. 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?
Take the Lesson 19 Quiz →