38 - Clustering: K-Means and DBSCAN without labels

You have a box of a thousand photographs with no captions and you have to sort them into piles. Nobody has told you which piles exist: “sea”, “mountains”, “parties” emerge by looking at which photos resemble each other. This is clustering: grouping similar examples without anyone having supplied the labels. It is unsupervised learning, and it brings with it a difficulty that classification did not have: there is no correct answer to measure yourself against. The result depends on the idea of “similarity” you choose, and an algorithm will always give you groups, even if there are none in the data.

Let us look at two algorithms with opposite philosophies. K-Means looks for K compact groups around centres. DBSCAN looks for dense regions separated by empty regions, and knows how to recognize noise.

K-Means: the objective

The number of groups K is fixed. Each group C_k has a centre μ_k (centroid). The algorithm looks for assignments and centres that minimize the inertia, the sum of the squared distances of each point from the centre of its own group:

J = Σ_{k=1..K} Σ_{x ∈ C_k} ‖x − μ_k‖²

The loop is the scheme to keep in mind: assign, update, repeat. A centroid is the “midpoint” of a group, the arithmetic mean of the coordinates of its points; it is not necessarily one of the data points.

Lloyd's algorithm and why it converges

1. Choose K initial centres.
2. Assignment: each point goes to the nearest centre.
3. Update: each centre becomes the mean of the points assigned to it.
4. Repeat 2-3 until the assignments no longer change.

Steps 2 and 3 are an alternation, not two independent phases: the first fixes the centres and chooses the assignments, the second fixes the assignments and chooses the centres. Each minimizes J with respect to only one of the two unknowns, and this is what makes possible the convergence proof that follows. Step 1 is the only non-deterministic one, and as you will see it matters more than it seems.

Why the mean? With the assignments fixed, we want the μ that minimizes Σ_x ‖x − μ‖² in a group. Differentiating with respect to μ and setting to zero:

∂/∂μ Σ_x ‖x − μ‖² = −2 Σ_x (x − μ) = 0      →      μ = (1/|C|) Σ_x x

Each step cannot increase J: the assignment chooses, for each point, the centre that minimizes its distance (so J does not grow with the centres fixed), and the update chooses the centre that minimizes J with the assignments fixed. J never goes up and can take a finite number of values, so the algorithm terminates. But it terminates in a local minimum, not necessarily the global one: it depends on the initial centres.

Initialization: k-means++

Starting with randomly chosen centres can give terrible solutions (two centres in the same true group, a true group with no centre). k-means++ chooses the first centre at random, then each following one among the points with probability proportional to D(x)², where D(x) is the distance from the nearest centre already chosen: the centres move away from one another.

Implementation by hand

import numpy as np
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score, silhouette_score

X, y_vero = make_blobs(n_samples=600, centers=4, cluster_std=[1.0, 1.5, 0.5, 1.0], random_state=7)

def dist2(X, C):
    return ((X[:, None, :] - C[None, :, :]) ** 2).sum(axis=2)      # (n, K)

def init_pp(X, k, rng):
    centri = [X[rng.integers(len(X))]]
    for _ in range(k - 1):
        d = dist2(X, np.array(centri)).min(axis=1)                   # distance to the nearest centre already chosen
        centri.append(X[rng.choice(len(X), p=d / d.sum())])          # probability proportional to D(x)^2
    return np.array(centri)

def kmeans(X, k, init="pp", seed=0, max_iter=100):
    rng = np.random.default_rng(seed)
    C = init_pp(X, k, rng) if init == "pp" else X[rng.choice(len(X), k, replace=False)]
    for _ in range(max_iter):
        etichette = dist2(X, C).argmin(axis=1)                        # assignment step
        nuovi = np.array([X[etichette == j].mean(axis=0) if (etichette == j).any() else C[j]
                          for j in range(k)])                         # update step
        if np.allclose(nuovi, C):
            break
        C = nuovi
    inerzia = dist2(X, C).min(axis=1).sum()
    return C, etichette, inerzia

C, lab, J = kmeans(X, 4)
sk = KMeans(n_clusters=4, n_init=10, random_state=0).fit(X)
print("inertia by hand: %.1f    sklearn: %.1f" % (J, sk.inertia_))
print("agreement with the true groups (ARI) by hand: %.3f   sklearn: %.3f" % (
    adjusted_rand_score(y_vero, lab), adjusted_rand_score(y_vero, sk.labels_)))

inertia by hand: 1295.9    sklearn: 1295.9
agreement with the true groups (ARI) by hand: 0.996   sklearn: 0.996

In the listing, dist2 computes the squared distances between all the points and all the centres with broadcasting, returning an (n, K) matrix. init_pp implements k-means++: it chooses the first centre at random, then for each following centre it computes for every point the distance to the nearest centre already chosen and picks a point with probability proportional to that value (rng.choice(..., p=d/d.sum())). kmeans alternates the two steps: argmin for the assignment, the group mean for the update (if a group stays empty its centre does not move); it stops when the centres do not change (np.allclose). It also returns the final inertia. The data are 4 Gaussian groups with different deviations (0.5-1.5), and y_vero holds the groups we generated.

The ARI (Adjusted Rand Index) measures the agreement between two partitions: 1 if they coincide, about 0 if they are independent. We can use it here only because we generated the data ourselves and know the true groups: in reality you do not know them, and that is exactly the problem of clustering.

A single start is not enough: local minima

for init in ("random", "pp"):
    inerzie = [kmeans(X, 4, init="pp" if init == "pp" else "random", seed=s)[2] for s in range(30)]
    print("init %-8s  30 starts: best inertia %.1f   worst %.1f   starts in the best minimum: %d/30" % (
        init, min(inerzie), max(inerzie), sum(abs(v - min(inerzie)) < 1.0 for v in inerzie)))

init random    30 starts: best inertia 1295.9   worst 14265.5   starts in the best minimum: 19/30
init pp        30 starts: best inertia 1295.9   worst 1295.9   starts in the best minimum: 30/30

The listing repeats the algorithm with 30 different seeds, with random initialization (X[rng.choice(len(X), k, replace=False)]: k points taken at random as centres) or with k-means++, and counts how many times the best inertia is reached (with a tolerance of 1.0).

With random initial centres, out of 30 starts only 19 reach the best inertia (1295.9) and the worst ends at 14265.5, eleven times worse. With k-means++ all 30 starts reach the best minimum on this dataset. Different starts therefore end in different solutions. The practical rule is the one scikit-learn applies: run the algorithm several times (n_init) with k-means++ and keep the solution with the lowest inertia.

centre        : (0,0)   (1,0)   (0,1)
university    : (10,0)  (11,0)  (10,1)
railway station: (20,0)  (21,0)  (20,1)

How do you choose K?

The inertia always goes down as K increases (with K = n it is zero), so it cannot be minimized: you look for the “elbow”, the point beyond which the gains become marginal. The elbow is often ambiguous. A more informative criterion is the silhouette coefficient. For a point i, let a(i) be the mean distance from the other points of its own group and b(i) the mean distance from the nearest of the other groups:

s(i) = ( b(i) − a(i) ) / max( a(i), b(i) )         values in [−1, +1]

Close to +1: the point is well inside its group and far from the others. Close to 0: it is on the border. Negative: it is probably in the wrong group. The silhouette of the clustering is the mean over all points. Below we compute it by hand and use it, together with the inertia, to choose K.

def silhouette_a_mano(X, lab):
    D = np.sqrt(dist2(X, X))
    s = np.zeros(len(X))
    for i in range(len(X)):
        stesso = (lab == lab[i]); stesso[i] = False
        a = D[i, stesso].mean()
        b = min(D[i, lab == j].mean() for j in np.unique(lab) if j != lab[i])
        s[i] = (b - a) / max(a, b)
    return s.mean()

print("silhouette by hand (K=4): %.4f   sklearn: %.4f" % (silhouette_a_mano(X, lab), silhouette_score(X, lab)))
for k in range(2, 9):
    m = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    print("K=%d  inertia %8.1f   silhouette %.3f" % (k, m.inertia_, silhouette_score(X, m.labels_)))

silhouette by hand (K=4): 0.7982   sklearn: 0.7982
K=2  inertia  18575.6   silhouette 0.580
K=3  inertia   5430.6   silhouette 0.761
K=4  inertia   1295.9   silhouette 0.798
K=5  inertia   1078.1   silhouette 0.705
K=6  inertia    921.1   silhouette 0.712
K=7  inertia    811.3   silhouette 0.619
K=8  inertia    710.7   silhouette 0.488

The function silhouette_a_mano computes the full distance matrix (square root of dist2), then for each point i it computes a as the mean of the distances from the points of the same group (excluding itself with stesso[i] = False) and b as the minimum, among the other groups, of the mean distance from that group. The final loop retrains K-Means for K from 2 to 8 with 10 starts and prints inertia and silhouette.

The silhouette peaks at K = 4 (0.798), the value with which we generated the data. The inertia always goes down (18575, 5430, 1296, 1078, ...); here the elbow at K = 4 is visible because the groups are well separated, but on real data it is rarely so sharp. The silhouette by hand coincides with the scikit-learn one (0.7982). Neither metric is an oracle: they are tools for formulating a hypothesis, which you then check with domain knowledge (do the groups make sense to someone who knows the data?).

Where K-Means gets it wrong: non-spherical shapes and different scales

K-Means implicitly assumes spherical groups, of similar size and separated by their centres: it assigns to the nearest centre, so the boundaries between groups are straight lines. Two interleaved moons do not have this shape. In the same way, if one feature has a much larger scale than the others, the Euclidean distance is dominated by it (as in KNN).

from sklearn.datasets import make_moons

Xm, ym = make_moons(n_samples=500, noise=0.07, random_state=0)
km = KMeans(n_clusters=2, n_init=10, random_state=0).fit(Xm)
print("K-Means on the two moons: ARI = %.3f   (1.0 = true groups recovered)" % adjusted_rand_score(ym, km.labels_))

K-Means on the two moons: ARI = 0.245   (1.0 = true groups recovered)

The listing generates 500 points on two semicircles with little noise (0.07) and runs K-Means with K = 2; it compares the result with the true groups through the ARI.

The ARI is 0.245, very far from 1: K-Means cuts the moons in half with a straight line. We need an algorithm that reasons about density, not about distance from a centre.

DBSCAN: density-based clustering

DBSCAN (Density-Based Spatial Clustering of Applications with Noise) has two parameters: a radius ε (eps) and a minimum number of points min_samples. It classifies each point in one of three ways:

• Core point: within radius ε around itself (itself included) there are at least min_samples points. It is in a dense region.
• Border point: it is not core, but lies within radius ε of a core point.
• Noise: neither core nor border. It receives the label −1.

The three categories have a logical order: first you decide who is core by looking at the counts alone, then you assign the border points by looking at the cores, and what is left is noise. According to the scikit-learn guide a cluster is a set of core points close to one another plus the non-core points close to a core. In our code min_samples also counts the point itself, as the parameter of DBSCAN does.

A cluster is the set of core points connected by density (reachable from one another through chains of nearby core points), plus the attached border points. The algorithm starts from an unvisited core point, explores the neighbours of the cores in a chain (a breadth-first visit), and when the chain is exhausted it closes the cluster:

from sklearn.cluster import DBSCAN

def dbscan(X, eps, min_samples):
    n = len(X)
    D = np.sqrt(dist2(X, X))
    vicini = [np.where(D[i] <= eps)[0] for i in range(n)]
    core = np.array([len(v) >= min_samples for v in vicini])
    lab = np.full(n, -1)                       # -1 = noise, until it is attached
    c = 0
    for i in range(n):
        if lab[i] != -1 or not core[i]:
            continue
        lab[i] = c
        coda = [i]
        while coda:
            p = coda.pop()
            if not core[p]:
                continue                       # border points get attached but do not expand
            for q in vicini[p]:
                if lab[q] == -1:
                    lab[q] = c
                    coda.append(q)
        c += 1
    return lab

lab_db = dbscan(Xm, eps=0.2, min_samples=5)
lib_db = DBSCAN(eps=0.2, min_samples=5).fit(Xm).labels_
print("clusters found (by hand): %d, noise points: %d" % (len(set(lab_db)) - (1 if -1 in lab_db else 0), (lab_db == -1).sum()))
print("ARI DBSCAN by hand on the moons:  %.3f" % adjusted_rand_score(ym, lab_db))
print("ARI DBSCAN sklearn:               %.3f" % adjusted_rand_score(ym, lib_db))
print("agreement between the two implementations (ARI): %.3f" % adjusted_rand_score(lab_db, lib_db))

clusters found (by hand): 2, noise points: 0
ARI DBSCAN by hand on the moons:  1.000
ARI DBSCAN sklearn:               1.000
agreement between the two implementations (ARI): 1.000

In the listing, vicini[i] holds the indices of the points within radius ε of i (itself included); core is a boolean vector. The outer loop scans the points and, for each core point not yet assigned, opens a new cluster c and uses a stack (coda) to explore: it pops a point, and if it is core it adds to the cluster all its neighbours not yet assigned. Border points are labelled but, not being core, do not extend the search. Those never reached stay at −1 (noise). The final comparison with scikit-learn's DBSCAN gives ARI 1.000: same partitions.

DBSCAN follows the shape of the moons because it reasons by connection, not by closeness to a centre. There is no need to fix the number of groups: it emerges from the data.

The DBSCAN trap: choosing ε

The results depend enormously on ε. In the experiment below, on the moons, the 99th percentile of the k-distance is 0.170: a radius of 0.2 covers almost all the points. With ε = 0.05 we find 30 fragments and 253 noise points, with 0.10 three clusters, with 0.20 the two true ones (ARI 1.000), with 0.40 everything merges into a single cluster (ARI 0).

In general, an ε that is too small fragments everything or downgrades it to noise, an ε that is too large merges everything into a single cluster. A practical guide is the k-distance plot: for each point compute the distance from its min_samples-th neighbour, sort the values in increasing order and look for the “elbow”, where the points go from dense zones to isolated ones.

kd = np.sort(np.sort(np.sqrt(dist2(Xm, Xm)), axis=1)[:, 5])          # distance to the 5th neighbour (index 0 = itself)
print("percentiles of the k-distance (k=5):", {p: round(float(np.percentile(kd, p)), 3) for p in (50, 90, 95, 99)})
for eps in (0.05, 0.1, 0.2, 0.4, 1.0):
    l = dbscan(Xm, eps=eps, min_samples=5)
    print("eps=%.2f  clusters %d  noise %3d   ARI %.3f" % (eps, len(set(l)) - (1 if -1 in l else 0), (l == -1).sum(), adjusted_rand_score(ym, l)))

percentiles of the k-distance (k=5): {50: 0.071, 90: 0.11, 95: 0.128, 99: 0.17}
eps=0.05  clusters 30  noise 253   ARI 0.020
eps=0.10  clusters 3  noise  18   ARI 0.706
eps=0.20  clusters 2  noise   0   ARI 1.000
eps=0.40  clusters 1  noise   0   ARI 0.000
eps=1.00  clusters 1  noise   0   ARI 0.000

The k-distance line sorts the distance matrix row by row and takes column 5 (the distance from the fifth neighbour, column 0 being the point itself), then sorts the values again; some percentiles are printed. In the loop DBSCAN is tried with five values of ε, comparing the number of clusters, the noise points and the ARI.

Other pitfalls: ε is a single radius for the whole dataset, so DBSCAN struggles when the groups have different densities; it is distance-based, so it requires scaled features; and in high dimension distances concentrate (the curse of dimensionality described in KNN) and the concept of density weakens.

Trade-off between the two

• K-Means: fast, scales to millions of points, interpretable centres (average profiles), but requires K, assumes spherical groups, does not handle noise (every point goes to a group), depends on the initialization.
• DBSCAN: arbitrary shapes, recognizes noise, does not require K, but is sensitive to ε and min_samples, does not handle different densities, costs more on large datasets (the neighbour search; with spatial indexes it is acceptable in low dimension).

The two lists are mirror images: K-Means asks for the number of groups and imposes their shape, DBSCAN asks for the density and leaves the shape free. Neither wins in general. Choose K-Means if the groups are compact and you need an average profile per group, DBSCAN if the shape is irregular or you care about knowing which points belong to nothing.

Real-world traps

• Taking the clusters for truth: the algorithm finds groups even in uniform noise. Always ask whether the groups make sense and are stable (rerun with other seeds and on subsamples).
• Forgetting scaling: age in years and income in euros in the same K-Means = income decides everything.
• Categorical features encoded as numbers: K-Means on means makes no sense on categories. Variants (K-Modes) or other representations are needed.
• Using ARI or accuracy in production: they need true labels, which by definition you do not have. Use silhouette, stability and validation on the domain.
• Leakage in the preparation phase: if you cluster to create a feature to give to a supervised model, the clustering must be trained on the training set only (inside the pipeline).

The first three traps are preparation errors, the last two evaluation errors. In all of them the clustering runs and returns groups, because an unsupervised algorithm has no way of telling you that the answer is senseless. The section that follows measures two of these errors with numbers.

Feature scale and useless features: two ways to ruin a clustering

In clustering, where there are no labels to check, data-preparation errors are even more insidious than in the supervised case: no test metric warns you that something is wrong. We show two experiments. In the first, a column of pure noise with standard deviation 50 is added to the three well-separated groups from before (that is, a “large” quantity with no information, like a numeric identifier or an unnormalized amount). In the second, only the unit of measure of one of the two true features is changed (from feet to centimetres, a factor of 30.48), without adding anything.

from sklearn.metrics import silhouette_score
from sklearn.preprocessing import StandardScaler

rng2 = np.random.default_rng(0)
X_true, y_true = make_blobs(n_samples=600, centers=3, cluster_std=1.0, random_state=4)
X_noisy = np.hstack([X_true, rng2.normal(0, 50, (600, 1))])                      # one noise column on a very large scale

def evaluate(name, A):
    km = KMeans(3, n_init=10, random_state=0).fit(A)
    print("%-38s  ARI %.3f    silhouette %.3f" % (name, adjusted_rand_score(y_true, km.labels_), silhouette_score(A, km.labels_)))

evaluate("only the two true features", X_true)
evaluate("plus a noise column (std 50)", X_noisy)
evaluate("... after standardization", StandardScaler().fit_transform(X_noisy))

print()
evaluate("true features, original units", X_true)
evaluate("first feature multiplied by 30.48", np.hstack([X_true[:, :1] * 30.48, X_true[:, 1:]]))

only the two true features              ARI 0.932    silhouette 0.613
plus a noise column (std 50)            ARI 0.001    silhouette 0.508
... after standardization               ARI 0.498    silhouette 0.439

true features, original units           ARI 0.932    silhouette 0.613
first feature multiplied by 30.48       ARI 0.492    silhouette 0.626

The function evaluate runs K-Means with K = 3 and reports two numbers: the ARI with the true groups (which we know here because the data are synthetic) and the silhouette (which in practice is the only thing you would have). With the true features alone the ARI is 0.932 and the silhouette 0.613. Adding the noise column the ARI collapses to 0.001: the clustering is random, because the Euclidean distance is dominated by the column with deviation 50 and K-Means separates the points on the basis of the noise. Standardizing helps, but is not enough: the ARI rises only to 0.498, because the noise column still has the same weight as the two informative columns (a third each) and keeps interfering.

Note the silhouette: 0.508 for the clustering ruined by noise, a value that would seem “good”, because the groups are well separated along the noise dimension. The silhouette measures the compactness and geometric separation of the groups found, not their meaning: it does not tell you whether the groups correspond to anything real.

The second experiment shows that even just changing the unit of measure of a feature, without adding anything, moves the ARI from 0.932 to 0.492 (the silhouette, meanwhile, rises to 0.626): a different scale weighs the feature more in the distance and changes the groups. Clustering is not invariant with respect to units, unlike many supervised models. The practical consequences are three: always standardize (or choose the scales deliberately); select the features (remove those without information and the identifiers); and do not trust an internal metric such as the silhouette to say that a clustering “is right”.

Gaussian Mixture and the EM algorithm: K-Means with uncertainty and shapes

K-Means assigns each point to a single group (“hard assignment”) and assumes spherical groups of similar size. A Gaussian mixture model (GMM) generalizes: the data are generated by K Gaussians, each with its own mean μ_k, variance (or covariance matrix) σ²_k and weight π_k (the fraction of points that belong to it). Each point has a probability of belonging to each component (soft assignment), and the groups can be elliptical and of different sizes. It is trained with the EM algorithm (Expectation-Maximization), which has the same two-step structure as K-Means:

E (expectation): compute the responsibility of each component for each point
    r_ik = π_k · N(x_i | μ_k, σ²_k) / Σ_j π_j · N(x_i | μ_j, σ²_j)

M (maximization): update the parameters as means weighted by the responsibilities
    N_k = Σ_i r_ik        μ_k = Σ_i r_ik x_i / N_k        σ²_k = Σ_i r_ik (x_i − μ_k)² / N_k        π_k = N_k / n

The E step is Bayes' theorem (the same as in Naive Bayes) with the components as hidden classes: the “posterior probability” that point i comes from component k. The M step is K-Means with weighted means instead of “hard” ones. A theoretical result guarantees that the likelihood of the data never decreases at each iteration: we verify it by hand on one-dimensional data generated from two Gaussians (300 points from N(−2, 1) and 200 from N(3, 1.5)).

from sklearn.mixture import GaussianMixture

data = np.concatenate([rng2.normal(-2, 1.0, 300), rng2.normal(3, 1.5, 200)])
mu, sd, pi = np.array([-1.0, 1.0]), np.array([1.0, 1.0]), np.array([0.5, 0.5])      # rough starting point
normal_pdf = lambda x, m, s: np.exp(-0.5 * ((x - m) / s) ** 2) / (s * np.sqrt(2 * np.pi))

log_likelihood = []
for iteration in range(40):
    weighted = np.stack([pi[k] * normal_pdf(data, mu[k], sd[k]) for k in range(2)], axis=1)     # (n, 2): pi_k N(x_i | mu_k, sd_k)
    log_likelihood.append(np.log(weighted.sum(axis=1)).sum())
    r = weighted / weighted.sum(axis=1, keepdims=True)                              # E: responsibilities
    Nk = r.sum(axis=0)                                                              # M: new parameters
    mu = (r * data[:, None]).sum(axis=0) / Nk
    sd = np.sqrt((r * (data[:, None] - mu) ** 2).sum(axis=0) / Nk)
    pi = Nk / len(data)

print("log-likelihood at the first 5 iterations:", np.round(log_likelihood[:5], 2), "...  final %.2f" % log_likelihood[-1])
print("never decreases: %s" % np.all(np.diff(log_likelihood) > -1e-9))
print("EM by hand:     means %s   deviations %s   weights %s" % (np.round(mu, 3), np.round(sd, 3), np.round(pi, 3)))
gmm = GaussianMixture(2, random_state=0).fit(data[:, None])
print("scikit-learn:   means %s   deviations %s   weights %s" % (np.round(gmm.means_.ravel(), 3), np.round(np.sqrt(gmm.covariances_.ravel()), 3), np.round(gmm.weights_, 3)))
print("true values:    means [-2, 3]   deviations [1, 1.5]   weights [0.6, 0.4]")

log-likelihood at the first 5 iterations: [-1619.43 -1078.2  -1075.2  -1074.98 -1074.96] ...  final -1074.95
never decreases: True
EM by hand:     means [-2.064  2.75 ]   deviations [0.888 1.576]   weights [0.591 0.409]
scikit-learn:   means [-2.05   2.792]   deviations [0.9  1.54]   weights [0.596 0.404]
true values:    means [-2, 3]   deviations [1, 1.5]   weights [0.6, 0.4]

The listing starts from rough estimates (means −1 and 1, deviations 1, weights 0.5) and for 40 iterations repeats: compute the weighted probabilities of each point for each component (with the Gaussian density written by hand), sum them to obtain the log-likelihood, normalize to get the responsibilities (E step), and update means, deviations and weights as weighted means (M step). The log-likelihood goes from −1,619.43 at the first iteration to −1,078.2 at the second and converges to −1,074.95: it rises at every step and never falls, as guaranteed by the theory (the check np.all(np.diff(...) > -1e-9) gives true).

The parameters found (means −2.064 and 2.750, deviations 0.888 and 1.576, weights 0.591 and 0.409) are close to the values that generated the data (−2 and 3; 1 and 1.5; 0.6 and 0.4), and coincide with those of scikit-learn up to small convergence differences (−2.050, 2.792; 0.900, 1.540; 0.596, 0.404). The gap from the true values is the statistical error of a sample of 500 points, not a defect of the algorithm.

Elliptical groups: where K-Means gets it wrong and the mixture does not

T = np.array([[2.5, 0.0], [1.8, 0.4]])                                              # matrix that stretches and tilts the groups
X_ell = np.vstack([rng2.normal(size=(300, 2)) @ T + shift for shift in ([0, 0], [0, 6], [5, 3])])
y_ell = np.repeat([0, 1, 2], 300)

km = KMeans(3, n_init=10, random_state=0).fit(X_ell)
gmm = GaussianMixture(3, covariance_type="full", n_init=5, random_state=0).fit(X_ell)
print("ARI with the true groups:   K-Means %.3f    Gaussian mixture %.3f" % (adjusted_rand_score(y_ell, km.labels_), adjusted_rand_score(y_ell, gmm.predict(X_ell))))

print("\nchoosing the number of components with the BIC criterion (lower = better):")
for k in range(1, 7):
    print("   K=%d   BIC %.1f" % (k, GaussianMixture(k, n_init=3, random_state=0).fit(X_ell).bic(X_ell)))

ARI with the true groups:   K-Means 0.555    Gaussian mixture 1.000

choosing the number of components with the BIC criterion (lower = better):
   K=1   BIC 9165.7
   K=2   BIC 8041.7
   K=3   BIC 7196.1
   K=4   BIC 7238.3
   K=5   BIC 7266.5
   K=6   BIC 7305.1

The listing generates three elongated, tilted Gaussian groups (by applying a matrix T to spherical points). K-Means, which assumes spherical groups, cuts the groups with straight boundaries and obtains an ARI of 0.555; the mixture with full covariances, which models the elliptical shape of each component, recovers them perfectly (ARI 1.000).

The second part shows how to choose K with the BIC criterion (Bayesian Information Criterion): the log-likelihood penalized by the number of parameters, which unlike the likelihood itself does not always grow with K. The BIC falls from 9,165.7 (one group) to 8,041.7 (two) and reaches its minimum, 7,196.1, at K = 3, the true value; for K = 4, 5, 6 it rises again (7,238.3, 7,266.5, 7,305.1). Unlike the K-Means inertia, which always decreases, the BIC has a true minimum, and that is why it is a more principled selection criterion (under the assumption, to be checked, that the groups really are Gaussian).

Stability and hierarchical clustering

A “right” clustering should be stable: if you resample the data and redo the analysis, you should get the same groups. It is measured with the bootstrap: samples are drawn with replacement, K-Means is trained on each, the original points are assigned to the clustering thus obtained, and the ARI between this and the clustering on the whole dataset is computed. Values close to 1 indicate robust groups. Then we look at an alternative that does not require fixing K: agglomerative hierarchical clustering, which starts with every point as a group of its own and repeatedly merges the two closest groups, according to a linkage criterion: single (distance between the two closest points of the groups), complete (distance between the farthest), average (mean), ward (minimizes the increase in internal variance).

from sklearn.cluster import AgglomerativeClustering

X_four, _ = make_blobs(n_samples=600, centers=4, cluster_std=0.9, random_state=3)

def stability(X, k, repetitions=20):
    reference = KMeans(k, n_init=5, random_state=0).fit(X)
    aris = []
    for b in range(repetitions):
        sample = np.random.default_rng(b).integers(0, len(X), len(X))                # bootstrap sample
        km = KMeans(k, n_init=5, random_state=b).fit(X[sample])
        aris.append(adjusted_rand_score(reference.labels_, km.predict(X)))            # comparison on the whole dataset
    return np.mean(aris), np.std(aris)

print("K    stability (mean ARI between repeated clusterings)")
for k in (2, 3, 4, 5, 6):
    mean_ari, std_ari = stability(X_four, k)
    print("%-4d  %.3f +- %.3f" % (k, mean_ari, std_ari))

X_moons, y_moons = make_moons(500, noise=0.07, random_state=0)
print("\nhierarchical clustering on the two moons, ARI with the true groups:")
for linkage in ("ward", "complete", "average", "single"):
    print("   %-9s %.3f" % (linkage, adjusted_rand_score(y_moons, AgglomerativeClustering(2, linkage=linkage).fit_predict(X_moons))))

K    stability (mean ARI between repeated clusterings)
2     1.000 +- 0.000
3     0.899 +- 0.171
4     1.000 +- 0.000
5     0.834 +- 0.047
6     0.752 +- 0.080

hierarchical clustering on the two moons, ARI with the true groups:
   ward      0.445
   complete  0.523
   average   0.429
   single    1.000

The function stability compares, with the ARI, the reference clustering with those obtained from 20 resamplings. On four true groups, stability is perfect (1.000) with K = 2 and K = 4, and lower with K = 3 (0.899 ± 0.171: sometimes two true groups are merged one way, sometimes the other), K = 5 (0.834) and K = 6 (0.752). Two lessons. First, stability signals the “natural” K (here 4), but K = 2 also turns out perfectly stable, because two large pairs of groups are a robust structure: stability does not identify a single K, it identifies the scales at which structure exists. Second, K that are too large are unstable because K-Means arbitrarily splits compact groups, and this is a good warning sign of over-clustering.

The second experiment is a reminder of the importance of the linkage criterion: on the two moons, single linkage obtains a perfect ARI (1.000), because it follows chains of nearby points (like DBSCAN), while Ward (0.445), complete (0.523) and average (0.429) look for compact groups and cut the moons in half. Here too the choice of algorithm encodes a hypothesis about what a group is.

Three exercises to try

Exercise 1. In the noise experiment, compute the ARI after deliberately removing the noise column with an automatic criterion (for example by measuring, for each feature, the silhouette of the clustering obtained with that feature alone). Exercise 2. With the Gaussian mixture on the elliptical groups, compare covariance_type equal to “spherical”, “diag” and “full”: which one reconstructs the groups best, and which has the lowest BIC? Exercise 3. Compute the stability of DBSCAN as ε varies on the two moons: which values of ε are stable?

To practice on the domain of bike-sharing station usage profiles there is a Python exercise skeleton with three parts to complete: Lloyd's algorithm written from scratch with the check that the SSE does not grow, the initialization that ends in a local minimum 38 times worse, and DBSCAN labelled by hand. As long as the points to complete are empty, it prints which exercises are missing.

What you can do now

You can write K-Means and DBSCAN from scratch and explain why K-Means converges but to a local minimum, you can choose K with inertia, silhouette, BIC and stability knowing that none is an oracle, and you can recognize when the shape of the groups or the scale of the features makes the result meaningless. You can use a Gaussian mixture when the groups are elliptical. In the next lesson, dimensionality reduction, we reduce the number of dimensions while keeping the information that matters.

Sources

Previous lesson in the course: Boosting: AdaBoost, Gradient Boosting and a note on XGBoost. Other lessons in the series that link to this one: Decision trees: choosing a question with Gini and entropy.