Clustering

Problem

You have customer records, documents or embedding vectors and no labels. You want to know whether they fall into natural groups, how many, and which point belongs where, so you can segment customers, deduplicate near-copies, compress a dataset into a few representatives, or build the coarse partitions that a vector index searches.

"Natural groups" is not one well-defined thing, which is why there are several algorithms. Should a group be a compact ball around a centre, a dense region of any shape, or a bump in a probability density? Should every point belong somewhere, or can some be noise? Can a point belong partly to two groups? And with no labels to score against, how do you even know you got it right? Each algorithm below answers these questions differently.

Intuition

k-means is the "post offices" picture: place k offices so that the total squared walking distance from every house to its nearest office is as small as possible. Alternate two easy steps: every house goes to its nearest office, then every office moves to the middle of its houses. Repeat until nothing moves. It always finds compact, roughly round groups, even when the data has none.

DBSCAN is the "crowds" picture: a point standing in a crowd (enough neighbours within a fixed radius) is a core point, crowds that touch merge into one cluster, people at the edge of a crowd join it, and loners are noise. It finds groups of any shape and refuses to force outliers into a cluster, but needs one density scale for the whole dataset.

A Gaussian mixture is the "overlapping spotlights" picture: the data was drawn from k bell-shaped blobs, each with its own centre, shape and size. A point between two spotlights is lit by both, so it belongs, say, 70% to one and 30% to the other. That is a soft version of k-means that also learns each cluster's spread and orientation.

Both runs are fixed points of Lloyd's algorithm: no point wants to switch centre and no centre wants to move. Only one is the good answer, which is why the starting centres matter.

Mechanics

k-means and Lloyd's algorithm

Given points x₁ … x_n and a number k, k-means looks for centres μ₁ … μ_k minimizing the inertia, the within-cluster sum of squares:

inertia = Σ_i min_j ‖x_i − μ_j‖²

Finding the exact minimum is NP-hard, so everyone uses Lloyd's algorithm, a local search with two steps:

1. Assign: give each point the label of its nearest centre (a nearest-neighbour query against k points).
2. Update: move each centre to the mean of its points, the position that minimizes that cluster's squared distances.

Each step can only lower the inertia or leave it unchanged, and there are finitely many ways to label the points, so the loop stops at a fixed point. That is a local minimum: the picture above shows one with inertia seven times the best. Each iteration costs O(n·k·d). scikit-learn stops when the centres' total squared shift falls below tol times the data's mean per-feature variance, or after max_iter = 300 rounds; an empty cluster gets re-seeded with a far-away point. Its algorithm="elkan" option gives the same answer faster on some data by using the triangle inequality to skip distance computations.

k-means++ seeding

Arthur and Vassilvitskii pick starting centres one at a time: the first uniformly at random, then each next centre with probability proportional to D(x)², the squared distance from x to its nearest already-chosen centre:

P(choose x) = D(x)² / Σ_x′ D(x′)²

Points far from every current centre, likely in an uncovered cluster, are favoured, while outliers are only occasionally picked because they are few. The seeding alone is guaranteed to be within O(log k) of the optimal inertia in expectation, and Lloyd's iterations then improve it. scikit-learn uses a greedy variant: at each step it draws 2 + ⌊ln k⌋ candidates and keeps the one that lowers the inertia most. Because restarts are still cheap insurance, KMeans runs n_init seedings and keeps the lowest inertia; with k-means++ the default n_init="auto" means one run, with random seeding ten.

Choosing k

Inertia always falls as k grows (k = n gives zero), so it cannot choose k by itself. Three common tools:

• Elbow: plot inertia against k and look for the point where adding a cluster stops paying off much. Simple, but the elbow is often vague.
• Silhouette: for each point, a = mean distance to the other points in its own cluster, b = mean distance to the points of the nearest other cluster, and s = (b − a) / max(a, b). s is near 1 for a well-placed point, near 0 on a boundary, negative if it is probably in the wrong cluster; a point alone in its cluster gets 0. Pick the k with the highest mean s.
• BIC for mixtures: BIC = −2 · log-likelihood + p · ln n, where p counts free parameters. More components raise the likelihood but pay ln n per parameter; pick the lowest BIC.

Domain knowledge beats all three. If you need 5 customer segments for 5 sales teams, k is 5. The worked example below shows the three criteria agreeing, which happens when the clusters really are separate blobs; on real data they often disagree.

DBSCAN

Ester, Kriegel, Sander and Xu define clusters by density with two parameters: a radius ε (eps) and a count min_samples.

• A core point has at least min_samples points, itself included, within distance ε.
• A border point is not core but lies within ε of a core point.
• Noise is everything else, labelled −1.

A cluster is a maximal set of core points connected through ε-neighbourhoods, plus their border points. scikit-learn computes every point's ε-neighbourhood with a radius query, marks the core points, then runs a depth-first search from each unlabelled core point, expanding only through core points. There is no k: the number of clusters falls out of the density. Clusters can have any shape. The catch is the single ε: clusters of very different densities cannot all be captured at once, and memory can grow as O(n²) when ε is large. HDBSCAN (in sklearn.cluster) varies the density threshold to fix that. A border point within ε of two clusters' core points goes to whichever cluster the search reaches first.

Gaussian mixtures and EM

A Gaussian mixture models the data's density as a weighted sum of k Gaussians with weights π_j, means μ_j and covariances Σ_j:

p(x) = Σ_j π_j · N(x; μ_j, Σ_j)

Maximizing the likelihood directly is awkward because each point's cluster is unknown. Expectation–maximization (Dempster, Laird and Rubin) alternates two steps, the soft analogue of Lloyd's:

1. E-step: compute each point's responsibilities, the posterior probability of each component, r_ij = π_j N(x_i; μ_j, Σ_j) / Σ_l π_l N(x_i; μ_l, Σ_l), in log space for stability.
2. M-step: re-estimate each component from responsibility-weighted points: N_j = Σ_i r_ij, π_j = N_j/n, μ_j = Σ_i r_ij x_i / N_j, Σ_j = Σ_i r_ij(x_i − μ_j)(x_i − μ_j)ᵀ / N_j, plus a small reg_covar on the diagonal so a collapsing component cannot reach zero variance.

Each round never lowers the log-likelihood, and like Lloyd's it stops at a local optimum. scikit-learn initializes the responsibilities from a k-means run by default and stops when the mean log-likelihood changes by less than tol = 10⁻³. With a full covariance per component in d dimensions, the parameter count is k·d(d+1)/2 + k·d + (k − 1); covariance_type can be "diag", "tied" or "spherical" to cut it. k-means is the limit of a mixture with equal weights and shared spherical covariances shrinking to zero, where responsibilities harden to 0 or 1.

Evaluating a clustering

With no ground truth, use internal scores that reward tight, well-separated clusters: the silhouette above, the Davies–Bouldin index (lower is better) or the Calinski–Harabasz score (higher is better). All of them favour round, convex clusters, so they will rate DBSCAN's correct answer on curved data poorly. When labels exist for a sample, use external scores. The adjusted Rand index (ARI) counts pairs of points that both clusterings put together or apart, corrected so a random labelling scores about 0 and an identical one scores 1, whatever the cluster names; normalized mutual information is the information-theoretic alternative. Finally, look at the clusters: plot them, read samples from each, and check whether they are stable across seeds and subsamples.

Worked example

DBSCAN never had to be told how many clusters there were, and it left the outliers out instead of dragging a centre toward them.

Implementation

import numpy as np
from sklearn.cluster import KMeans, DBSCAN
from sklearn.mixture import GaussianMixture
from sklearn.metrics import silhouette_score, adjusted_rand_score
from sklearn.preprocessing import StandardScaler

X = StandardScaler().fit_transform(X_raw)         # distances need comparable feature scales

for k in range(2, 8):
    km = KMeans(n_clusters=k, n_init="auto", random_state=0).fit(X)
    gm = GaussianMixture(n_components=k, covariance_type="full", random_state=0).fit(X)
    print(k, km.inertia_, silhouette_score(X, km.labels_), gm.bic(X))

km = KMeans(n_clusters=3, n_init="auto", random_state=0).fit(X)
db = DBSCAN(eps=0.3, min_samples=5).fit(X)        # label -1 = noise
gm = GaussianMixture(n_components=3, random_state=0).fit(X)
soft = gm.predict_proba(X)                        # (n, 3) responsibilities

if y_sample is not None:                          # labels for a sample, if you have them
    print(adjusted_rand_score(y_sample, km.labels_[sample_idx]))

def lloyd(X, centers, iters=100):                 # the algorithm itself, in NumPy
    for _ in range(iters):
        d2 = ((X[:, None, :] - centers[None]) ** 2).sum(-1)    # (n, k)
        labels = d2.argmin(1)
        new = np.array([X[labels == j].mean(0) if (labels == j).any() else centers[j]
                        for j in range(len(centers))])
        if np.allclose(new, centers):
            break
        centers = new
    return labels, centers, d2.min(1).sum()        # inertia

Scale the features first: k-means and DBSCAN use raw Euclidean distance, so a feature measured in thousands drowns out one measured in fractions (see feature preprocessing). A common way to set DBSCAN's ε is to sort every point's distance to its min_samples-th nearest neighbour and pick the distance where that curve bends upward. For millions of points, MiniBatchKMeans updates centres from small random batches at some cost in inertia.

Tradeoffs

• k-means is the default for speed and simplicity and for "compress into k representatives" jobs such as colour quantization and the coarse quantizer of a vector index, where convex cells are exactly what you want.
• DBSCAN suits spatial data and anomaly detection, where noise matters and the shapes are irregular. In high dimensions distances concentrate, so a good ε is hard to find; reduce dimensions first (with PCA, say).
• Gaussian mixtures give probabilities and a density you can sample from or use for anomaly scores, but full covariances need many points per component and can collapse onto a few points without reg_covar.
• All of them find clusters in random noise. Check stability across seeds and subsamples, and prefer a clustering that serves a decision over one that merely scores well.

Connections

Nearest neighbours supplies k-means' assign step and DBSCAN's radius queries; vector indexes use k-means to partition embeddings for IVF and to build product-quantization codebooks; PCA is the usual step before clustering high-dimensional data.

References

• Lloyd. Least Squares Quantization in PCM. IEEE Transactions on Information Theory 1982. doi:10.1109/TIT.1982.1056489
• Arthur and Vassilvitskii. k-means++: The Advantages of Careful Seeding. SODA 2007.
• Ester, Kriegel, Sander and Xu. A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. KDD 1996.
• Dempster, Laird and Rubin. Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society B 1977.
• Rousseeuw. Silhouettes: a Graphical Aid to the Interpretation and Validation of Cluster Analysis. Journal of Computational and Applied Mathematics 1987. doi:10.1016/0377-0427(87)90125-7
• Hubert and Arabie. Comparing Partitions. Journal of Classification 1985. doi:10.1007/BF01908075
• Schwarz. Estimating the Dimension of a Model. Annals of Statistics 1978. doi:10.1214/aos/1176344136
• Campello, Moulavi and Sander. Density-Based Clustering Based on Hierarchical Density Estimates. PAKDD 2013.
• scikit-learn 1.9.1: KMeans, DBSCAN, GaussianMixture, silhouette_score.