Clustering (K-Means, DBSCAN, GMM)

TL;DR

Intuition

Three families partition data in three different ways:

• Partitional (K-Means, GMM). Pick $k$ representatives (centroids or Gaussians). Assign each point to the closest one. Re-fit the representatives. Repeat. Output: a flat partition of the data.
• Density-based (DBSCAN, OPTICS). A cluster is a maximal region of high point density, separated from other clusters by low-density gaps. Output: clusters of arbitrary shape, plus a “noise” label for points in low-density regions.
• Hierarchical (agglomerative, divisive). Start with each point in its own cluster, repeatedly merge the closest pair (agglomerative); or start with everything together and recursively split (divisive). Output: a tree of nested partitions (dendrogram), readable at any granularity.

Math

DBSCAN: density-based clustering

DBSCAN [4] defines clusters by density rather than centroids, which lets it find arbitrary shapes and label outliers as noise.

Key concepts

Fix two hyperparameters: a radius $\epsilon$ and a minimum count $m = \mathrm{minPts}$.

• Core point. A point $x$ is core if its $\epsilon$-neighborhood contains at least $m$ points (including itself): $|N_\epsilon(x)| \ge m$.
• Directly density-reachable. $y$ is directly density-reachable from $x$ if $x$ is a core point and $y \in N_\epsilon(x)$.
• Density-reachable. $y$ is density-reachable from $x$ if there is a chain of core points $x = p_0, p_1, \dots, p_\ell = y$ with each $p_{i+1}$ directly density-reachable from $p_i$.
• Density-connected. $x$ and $y$ are density-connected if there exists a core point $z$ from which both are density-reachable.

A cluster is a maximal set of density-connected points. A noise point is a non-core point that is not density-reachable from any core point.

Algorithm sketch

1. For each unvisited point $x$: count its $\epsilon$-neighbors. If at most $m - 1$, label noise (provisionally) and continue.
2. Otherwise mark $x$ as the next cluster id and BFS-expand the neighborhood: any neighbor that is itself a core point contributes its own neighbors.
3. Border points (in some core point’s neighborhood but not core themselves) get the cluster’s id; truly isolated points keep the noise label.

Why it succeeds where K-Means fails

K-Means assumes spherical, equal-variance clusters of comparable size. DBSCAN makes no shape assumption: it finds two interlocking moons, a ring around a ball, or any density-connected shape, and tags the low-density tail as noise rather than forcing it into the nearest cluster.

Practical knobs

• $\epsilon$: the most sensitive parameter. The k-distance plot (sort the distances to the $m$-th nearest neighbor and look for the elbow) is the standard heuristic.
• minPts: typically $\ge d + 1$; many practitioners use $m = 2 d$ as a starting point.
• Limitation: variable density. A single $\epsilon$ struggles when clusters have very different densities. OPTICS (1999) and HDBSCAN (2013) are the standard fixes.

Hierarchical clustering

Agglomerative hierarchical clustering builds a binary tree (dendrogram) [10]:

1. Start with each point as its own cluster.
2. Merge the two clusters with the smallest linkage distance.
3. Repeat until one cluster remains.

Cutting the dendrogram at any height yields a flat partition. Linkage choices:

Hierarchical clustering is $O(n^2)$ memory (the full distance matrix) and $O(n^2 \log n)$ to $O(n^3)$ time. It does not scale beyond a few tens of thousands of points, but it is interpretable: the dendrogram is the model.

Spectral clustering

Spectral clustering [6] embeds points via the eigenvectors of the graph Laplacian and runs K-Means in that space. The recipe:

1. Build a similarity graph (often a $k$-nearest-neighbor or RBF-weighted graph) with affinity matrix $W$.
2. Compute the symmetric normalized Laplacian $L_{\mathrm{sym}} = I - D^{-1/2} W D^{-1/2}$, where $D$ is the diagonal degree matrix.
3. Stack the bottom $k$ eigenvectors of $L_{\mathrm{sym}}$ as columns of $U \in \mathbb{R}^{n \times k}$.
4. Normalize each row of $U$ to unit length.
5. Run K-Means on the rows of $U$.

The geometry is non-obvious. The Laplacian’s bottom eigenvectors encode the graph’s connected-component structure; on a graph that is almost disconnected (a few weak inter-cluster edges), the bottom eigenvectors are nearly piecewise-constant on the components, so a simple K-Means in spectral space recovers them. Spectral clustering shines on non-convex shapes that defeat K-Means but does not scale to millions of points without approximations.

Cluster validity: choosing $k$

Without labels, “the right $k$” is ill-defined. The three standard scores:

In practice run all three over a range of $k$ and look for agreement. Silhouette is the cheapest and most popular; the gap statistic has the best statistical justification [7][8].

How the families compare

Implementations

The numpy version implements Lloyd’s algorithm from scratch in about 30 lines. The sklearn version is the production default. The MiniBatchKMeans version trades a sliver of accuracy for the ability to fit millions of rows.

import numpy as np

def kmeans(X, k, n_iters=100, seed=0):
  """Lloyd's algorithm with K-Means++ seeding."""
  rng = np.random.default_rng(seed)
  n, d = X.shape

  # K-Means++ seeding
  centers = np.empty((k, d))
  centers[0] = X[rng.integers(n)]
  for j in range(1, k):
      d2 = np.min(((X[:, None, :] - centers[:j][None, :, :]) ** 2).sum(-1), axis=1)
      probs = d2 / d2.sum()
      centers[j] = X[rng.choice(n, p=probs)]

  # Lloyd iterations
  for _ in range(n_iters):
      # Assignment
      dists = ((X[:, None, :] - centers[None, :, :]) ** 2).sum(-1)
      assign = dists.argmin(axis=1)
      # Update
      new_centers = np.array([
          X[assign == c].mean(axis=0) if (assign == c).any() else centers[c]
          for c in range(k)
      ])
      if np.allclose(new_centers, centers):
          break
      centers = new_centers

  inertia = float(((X - centers[assign]) ** 2).sum())
  return centers, assign, inertia

rng = np.random.default_rng(0)
X = np.vstack([rng.normal([-2, 0], 0.5, (40, 2)),
             rng.normal([+2, 0], 0.5, (40, 2)),
             rng.normal([0, 3], 0.5, (40, 2))])
centers, assign, inertia = kmeans(X, k=3)
print("centers:\n", np.round(centers, 2))
print("inertia:", round(inertia, 2))

Run sklearn’s K-Means live on a synthetic 3-blob dataset and inspect the cluster centers and inertia:

(no output yet)

Interactive demo

Pick a method (K-Means, DBSCAN, GMM) and a tunable hyperparameter to see how each algorithm partitions three overlapping Gaussian blobs. Note how DBSCAN labels gray points as noise when $\epsilon$ is too small, and how GMM with $k$ greater than 3 wastes a component on the densest blob.

Real-world uses

Pitfalls

Pop quiz

Interview Q&A