08 - K-Means Clustering and Hierarchical Clustering

You have a cloud of unlabeled points and you want to know whether it contains groups. Clustering is the task of dividing the points into groups whose members resemble one another more than they resemble the others, without anyone telling you in advance which groups exist. K-Means does it with two repeated moves: assign every point to the nearest center, then move every center to the mean of its points.

It always converges, but to a result that depends on where the centers start, it wants you to decide $k$ in advance, and it fails systematically when the groups are not compact. This article, the eighth in the ML Fundamentals series, writes the algorithm from scratch, measures what changes with the initialization, compares four ways of choosing $k$ and tries K-Means, hierarchical clustering, DBSCAN and Gaussian mixtures on data where you know the right answer.

To know whether a clustering is right you need a reference answer, and in practice there isn't one. Here there is, because the data are synthetic (generated with a fixed seed) and the true labels are known: to compare them with the ones found I use the Adjusted Rand Index (ARI), which is 1 if the two partitions coincide up to the names of the groups and about 0 if the agreement is what chance would give.

You do not need to have read anything else in the series, but the SVM of the previous article already showed you why units of measurement matter: here they still do.

Lloyd's algorithm, written from scratch

K-Means minimizes the inertia, that is the sum of the squared distances between every point and the center (centroid) of its group: $J = \sum_i \|x_i - \mu_{c(i)}\|^2$. Finding the absolute minimum is a hard problem; Lloyd's algorithm finds a local minimum by alternating two steps.

In the assignment step every point goes to the nearest centroid, with the centers fixed; in the update step every centroid becomes the mean of the assigned points, with the assignments fixed. Neither step can make $J$ grow, and since the possible partitions are finite, the algorithm stops: Arthur and Vassilvitskii's paper says so too, noting that the total error decreases monotonically.

Try it on six numbers on the line, 1, 2, 4, 7, 8 and 11, with two badly chosen starting centers, 1 and 2. In the first round point 1 goes to center 1 and all the others to center 2, which moves to 6.4; in the second round point 2 also moves to the first group, and so on.

The listing does the same computations, prints $J$ at every iteration and then compares your lloyd function with scikit-learn's KMeans on 600 points, starting from the same centers (init with an array and n_init=1), because only that way do the two algorithms follow exactly the same path.

# Lloyd's algorithm written from scratch, tried on 6 numbers and then compared with KMeans
import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

def lloyd(X, centers, max_iter=100, verbose=False):
    """Alternate assignment to the nearest centers and recomputation of the means."""
    C = np.array(centers, dtype=float)
    prev_lab = None
    for it in range(1, max_iter + 1):
        D = ((X[:, None, :] - C[None, :, :]) ** 2).sum(axis=2)   # squared point-center distances
        lab = D.argmin(axis=1)
        J = D[np.arange(len(X)), lab].sum()                      # inertia with the CURRENT centers
        if verbose:
            print(f"iter {it}: labels={lab.tolist()} centers={np.round(C.ravel(), 3).tolist()} J={J:.3f}")
        if prev_lab is not None and (lab == prev_lab).all():
            return lab, C, J, it
        C = np.array([X[lab == k].mean(axis=0) if (lab == k).any() else C[k] for k in range(len(C))])
        prev_lab = lab
    return lab, C, J, max_iter

# 1) Six numbers on a line, two badly chosen initial centers: 1 and 2
x = np.array([[1.0], [2.0], [4.0], [7.0], [8.0], [11.0]])
lab, C, J, it = lloyd(x, [[1.0], [2.0]], verbose=True)
print("final centers:", C.ravel().tolist(), " final J:", J, " iterations:", it)

# 2) Same initial centers on 600 2D points: our Lloyd against sklearn
X, _ = make_blobs(n_samples=600, centers=4, cluster_std=1.1, random_state=3)
C0 = X[:4]
my_lab, my_C, my_J, my_it = lloyd(X, C0)
km = KMeans(n_clusters=4, init=C0, n_init=1, algorithm="lloyd", tol=0.0, max_iter=100).fit(X)
print("same labels:", bool((my_lab == km.labels_).all()))
print(f"my inertia = {my_J:.6f}  sklearn = {km.inertia_:.6f}")
print("max |my centers - sklearn| =", f"{np.abs(my_C - km.cluster_centers_).max():.2e}")

# Output:
# iter 1: labels=[0, 1, 1, 1, 1, 1] centers=[1.0, 2.0] J=146.000
# iter 2: labels=[0, 0, 1, 1, 1, 1] centers=[1.0, 6.4] J=30.840
# iter 3: labels=[0, 0, 0, 1, 1, 1] centers=[1.5, 7.5] J=19.500
# iter 4: labels=[0, 0, 0, 1, 1, 1] centers=[2.333, 8.667] J=13.333
# final centers: [2.3333333333333335, 8.666666666666666]  final J: 13.333333333333332  iterations: 4
# same labels: True
# my inertia = 1454.239467  sklearn = 1454.239467
# max |my centers - sklearn| = 3.55e-15

The inertia drops at every round, 146, then 30.84, then 19.5, then 13.333, and the algorithm stops when the labels no longer change; the final centers 2.333 and 8.667 are the means of {1, 2, 4} and {7, 8, 11}.

On the 600 points the labels coincide with those of scikit-learn, the inertia is equal to the sixth decimal digit (1454.239467) and the centers differ by 3.55e-15: KMeans with algorithm="lloyd" does exactly this, with speed optimizations. The other side of guaranteed convergence is that the arrival point depends on the starting point, as the next section shows.

Initialization matters: k-means++ measured

The classic method picks $k$ points at random as initial centers. The paper by Arthur and Vassilvitskii (SODA 2007) proposes to pick the first at random and each of the following with probability proportional to the squared distance from the nearest center already chosen: the initial centers thus tend to be far from one another.

The authors prove that the expected inertia is at most $O(\log k)$ times the optimal one and report experiments in which k-means++ beats classic K-Means, often by a substantial margin. The KMeans documentation adds an operational detail: with init='k-means++' the default n_init='auto' means a single start, with init='random' ten.

To see how big the margin is, I ran 300 times, with different seeds, a single-start K-Means (n_init=1) with each of the two methods on 1,500 points in two dimensions that come from 10 true groups. Then I compared the distribution of the final inertias with the best found in the 600 starts, used as a reference for «almost optimal».

# k-means++ against random initialization: 300 starts per method, one at a time (n_init=1)
import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=1500, centers=10, cluster_std=0.9, center_box=(-12, 12), random_state=7)

def trial(init, n=300):
    J, it = [], []
    for seed in range(n):
        km = KMeans(n_clusters=10, init=init, n_init=1, random_state=seed).fit(X)
        J.append(km.inertia_); it.append(km.n_iter_)
    return np.array(J), np.array(it)

J_r, it_r = trial("random")
J_p, it_p = trial("k-means++")
best = min(J_r.min(), J_p.min())
print(f"best inertia found in 600 starts: {best:.2f}")
print(f"{'init':10s} {'mean':>9s} {'median':>9s} {'max':>9s} {'within 1% of best':>22s} {'mean iter':>11s}")
for name, J, it in [("random", J_r, it_r), ("k-means++", J_p, it_p)]:
    frac = (J <= best * 1.01).mean()
    print(f"{name:10s} {J.mean():9.1f} {np.median(J):9.1f} {J.max():9.1f} {frac:21.1%} {it.mean():11.1f}")
print(f"mean inertia of k-means++ relative to random: {J_p.mean() / J_r.mean() - 1:+.1%}")

# Output:
# best inertia found in 600 starts: 2224.20
# init            mean    median       max      within 1% of best   mean iter
# random        5530.3    4833.1   19243.9                  6.0%         8.4
# k-means++     2390.9    2435.6    4778.1                 45.7%         5.1
# mean inertia of k-means++ relative to random: -56.8%

The mean inertia with random initialization is 5530 and the median 4833, against 2391 and 2436 for k-means++: on average 57% less, on this dataset. The worst case of the random start reaches 19,244, more than eight times the best (2224). It also matters how many starts end close to the best: 6% for random initialization, 45.7% for k-means++, which therefore does not guarantee the optimum, and does not even find it in one start out of two.

The 57% figure must not be generalized: it depends on the separation of the groups and on their number. The advice that holds after this measurement is the documentation's: k-means++ as the default, and n_init=10 instead of 1, keeping the start with the minimum inertia.

How to choose k

The number of groups does not come out of the algorithm: you decide it. Inertia alone is not enough, because it always decreases when $k$ grows (with $k = n$ it is zero). Four tools use different criteria. The elbow looks at where the descent of the inertia slows down. Rousseeuw's silhouette gives every point $s = (b - a)/\max(a, b)$, with $a$ the mean distance from the points of its own group and $b$ the mean distance from the points of the nearest group; it goes from -1 to 1 and is averaged over the dataset.

The gap statistic by Tibshirani, Walther and Hastie (2001) compares the logarithm of the observed inertia with the one expected on reference data without structure, generated uniformly inside the bounding box of the data.

The paper's rule is precise: you compute $\mathrm{Gap}(k) = \frac{1}{B}\sum_b \log W^*_{kb} - \log W_k$ on $B$ reference sets, with $s_k = \mathrm{sd}_k\sqrt{1 + 1/B}$, and you choose the smallest $k$ such that $\mathrm{Gap}(k) \geq \mathrm{Gap}(k+1) - s_{k+1}$. The listing implements it with $B = 20$ and sets it beside inertia, silhouette and the Davies-Bouldin index (lower is better) on four true groups, two of which are close. First it checks the silhouette formula by hand on eight points.

# How to choose k: inertia, silhouette, Davies-Bouldin and gap statistic (Tibshirani, Walther, Hastie 2001)
import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import davies_bouldin_score, silhouette_samples, silhouette_score

X, _ = make_blobs(n_samples=500, centers=4, cluster_std=1.3, random_state=11)   # 4 true groups, slightly overlapping

# --- Hand check of the silhouette on 8 points ---
P = np.array([[0, 0], [0, 1], [1, 0], [5, 5], [5, 6], [6, 5], [3, 3], [9, 9]], float)
lab = np.array([0, 0, 0, 1, 1, 1, 0, 1])
D = np.sqrt(((P[:, None] - P[None]) ** 2).sum(-1))
s_mine = []
for i in range(len(P)):
    same = (lab == lab[i]) & (np.arange(len(P)) != i)
    a = D[i, same].mean()                        # mean distance within its own cluster
    b = D[i, lab != lab[i]].mean()               # mean distance from the other cluster
    s_mine.append((b - a) / max(a, b))
print("hand silhouette = sklearn:", np.allclose(s_mine, silhouette_samples(P, lab)))

# --- Table for k = 1..8 ---
rng = np.random.default_rng(0)
lo, hi, B = X.min(axis=0), X.max(axis=0), 20

def W(M, k):
    return KMeans(n_clusters=k, n_init=5, random_state=0).fit(M).inertia_

references = [rng.uniform(lo, hi, size=X.shape) for _ in range(B)]   # method (a) of the paper
ks = list(range(1, 9))
gap, s = [], []
for k in ks:
    log_ref = np.log([W(R, k) for R in references])
    gap.append(log_ref.mean() - np.log(W(X, k)))
    s.append(log_ref.std() * np.sqrt(1 + 1 / B))
gap, s = np.array(gap), np.array(s)

print(f"{'k':>2} {'inertia':>9} {'silhouette':>11} {'Davies-B':>9} {'gap':>7} {'gap(k)-gap(k+1)+s(k+1)':>24}")
for i, k in enumerate(ks):
    inert = W(X, k)
    if k > 1:
        km = KMeans(n_clusters=k, n_init=5, random_state=0).fit(X)
        sil, db = silhouette_score(X, km.labels_), davies_bouldin_score(X, km.labels_)
        ss, dd = f"{sil:11.3f}", f"{db:9.3f}"
    else:
        ss, dd = f"{'-':>11}", f"{'-':>9}"
    crit = f"{gap[i] - gap[i + 1] + s[i + 1]:24.3f}" if k < 8 else f"{'':>24}"
    print(f"{k:2d} {inert:9.1f} {ss} {dd} {gap[i]:7.3f} {crit}")
k_gap = next(k for i, k in enumerate(ks[:-1]) if gap[i] >= gap[i + 1] - s[i + 1])
print("k chosen by the gap statistic (the smallest with gap(k) >= gap(k+1) - s(k+1)):", k_gap)

# Output:
# hand silhouette = sklearn: True
#  k   inertia  silhouette  Davies-B     gap   gap(k)-gap(k+1)+s(k+1)
#  1   21602.0           -         -   0.150                   -0.154
#  2    9250.7       0.540     0.796   0.328                   -0.749
#  3    3010.8       0.675     0.439   1.100                   -0.213
#  4    1649.0       0.636     0.524   1.338                    0.164
#  5    1463.7       0.579     0.664   1.204                    0.128
#  6    1310.4       0.490     0.894   1.102                    0.057
#  7    1174.1       0.450     0.959   1.070                    0.036
#  8    1034.9       0.332     1.123   1.059
# k chosen by the gap statistic (the smallest with gap(k) >= gap(k+1) - s(k+1)): 4

The hand-computed silhouette coincides with scikit-learn's. In the table the inertia has an elbow between $k=3$ and $k=4$ (from 3011 to 1649, then only 1464), the silhouette is maximal at $k=3$ (0.675) and the Davies-Bouldin is minimal at 3 (0.439), while the gap statistic chooses 4, which is the true number.

The two closest centers are 4.87 apart with a standard deviation of 1.3: two of the four groups touch, and the silhouette, which rewards well-separated groups, prefers to merge them. No indicator always wins: if the criteria disagree, the extra information you have is the starting question (how many groups can you use?) and a plot of the points.

Where K-Means goes wrong

The documentation states that the inertia assumes convex and isotropic groups (same extent in every direction) and that it reacts badly to elongated groups or irregular shapes. To measure how badly, I generated five datasets with true labels: three spherical blobs; three horizontal stripes with staggered centers; three groups of very different sizes (500, 50 and 50 points); two moons; two concentric circles.

Each algorithm receives the right number of groups where needed ($k$ for K-Means, Ward and single linkage), and for DBSCAN I took eps without looking at the labels, as the 95th percentile of the distance from the fifth nearest point (the point itself included, as in min_samples). All the data are standardized beforehand. The listing also has a second experiment: a feature multiplied by 1000, like an amount in euros next to a count.

# Where K-Means goes wrong: five synthetic datasets with known true labels, measured with the Adjusted Rand Index
import numpy as np
from sklearn.cluster import DBSCAN, HDBSCAN, AgglomerativeClustering, KMeans
from sklearn.datasets import make_blobs, make_circles, make_moons
from sklearn.metrics import adjusted_rand_score as ari
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler

rs = 0
data = {}
X, y = make_blobs(n_samples=600, centers=3, cluster_std=0.8, random_state=rs); data["spherical blobs"] = (X, y, 3)
def stripes(seed=0, dx=2.0):
    """Three horizontally elongated groups (sigma 3 against 0.39), with staggered centers."""
    rng = np.random.default_rng(seed)
    cov = np.diag([9.0, 0.15])
    means = [(0, 0), (dx, 1.6), (2 * dx, 3.2)]
    X = np.vstack([rng.multivariate_normal(m, cov, 200) for m in means])
    return X, np.repeat([0, 1, 2], 200)
X, y = stripes(); data["three stripes"] = (X, y, 3)
X, y = make_blobs(n_samples=[500, 50, 50], centers=[[0, 0], [4, 0], [2, 3.5]], cluster_std=[1.3, 0.4, 0.4], random_state=rs)
data["different sizes"] = (X, y, 3)
X, y = make_moons(n_samples=500, noise=0.07, random_state=rs); data["two moons"] = (X, y, 2)
X, y = make_circles(n_samples=500, noise=0.05, factor=0.4, random_state=rs); data["two circles"] = (X, y, 2)

print(f"{'dataset':20s} {'KMeans':>7s} {'Ward':>7s} {'single':>7s} {'DBSCAN':>7s} {'HDBSCAN':>8s}")
for name, (X, y, k) in data.items():
    Z = StandardScaler().fit_transform(X)
    # eps without looking at the labels: the 95th percentile of the distance to the 5th neighbor
    eps = np.quantile(NearestNeighbors(n_neighbors=5).fit(Z).kneighbors(Z)[0][:, -1], 0.95)
    r = [ari(y, KMeans(k, n_init=10, random_state=0).fit_predict(Z)),
         ari(y, AgglomerativeClustering(k, linkage="ward").fit_predict(Z)),
         ari(y, AgglomerativeClustering(k, linkage="single").fit_predict(Z)),
         ari(y, DBSCAN(eps=eps, min_samples=5).fit_predict(Z)),
         ari(y, HDBSCAN(min_cluster_size=20, copy=True).fit_predict(Z))]
    print(f"{name:20s} " + " ".join(f"{v:7.2f}" for v in r))

# Units of measurement: a feature in "euros" (x1000) squashes the other in the Euclidean distance
X, y = make_blobs(n_samples=600, centers=3, cluster_std=0.8, random_state=rs)
Xu = X.copy(); Xu[:, 1] *= 1000
print("\nblobs with the second feature multiplied by 1000:")
print(f"  KMeans without scaling: ARI = {ari(y, KMeans(3, n_init=10, random_state=0).fit_predict(Xu)):.2f}")
print(f"  KMeans with StandardScaler: ARI = {ari(y, KMeans(3, n_init=10, random_state=0).fit_predict(StandardScaler().fit_transform(Xu))):.2f}")

# What single linkage and DBSCAN do on the three touching blobs
X, y = make_blobs(n_samples=600, centers=3, cluster_std=0.8, random_state=rs)
Z = StandardScaler().fit_transform(X)
print("\nsingle linkage, sizes of the 3 groups:", np.bincount(AgglomerativeClustering(3, linkage="single").fit_predict(Z)).tolist())
eps = np.quantile(NearestNeighbors(n_neighbors=5).fit(Z).kneighbors(Z)[0][:, -1], 0.95)
lab = DBSCAN(eps=eps, min_samples=5).fit_predict(Z)
print(f"DBSCAN with eps = {eps:.3f}: clusters found = {len(set(lab)) - (1 if -1 in lab else 0)}, noise points = {(lab == -1).sum()}")

# Output:
# dataset               KMeans    Ward  single  DBSCAN  HDBSCAN
# spherical blobs         0.93    0.89    0.00    0.00    0.46
# three stripes           0.54    0.61   -0.00   -0.00    0.42
# different sizes         0.29    0.26   -0.01   -0.05    0.62
# two moons               0.48    0.59    1.00    0.98    0.98
# two circles            -0.00    0.09    1.00    0.99    1.00
#
# blobs with the second feature multiplied by 1000:
#   KMeans without scaling: ARI = 0.53
#   KMeans with StandardScaler: ARI = 0.93
#
# single linkage, sizes of the 3 groups: [597, 1, 2]
# DBSCAN with eps = 0.323: clusters found = 1, noise points = 13

No algorithm wins on every row. K-Means and Ward do well on the spherical blobs (ARI 0.93 and 0.89) and badly elsewhere: 0.54 on the stripes, 0.29 on the groups of different sizes, 0.48 on the moons, -0.00 on the circles, where K-Means cuts the plane into two half-planes and cannot isolate the inner ring. Single linkage and DBSCAN, which follow the connection between nearby points, solve moons and circles (ARI 1.00 and 0.98-0.99) but on the spherical blobs they score zero: the three blobs touch, single linkage detaches only three isolated points (the last lines of the listing: groups of 597, 1 and 2 points) and DBSCAN merges everything into a single cluster and discards 13 points as noise. With different densities, on the groups of different sizes, only HDBSCAN exceeds 0.6 (0.62).

The last experiment concerns the unit of measurement. Multiplying one of the two coordinates by 1000 makes the ARI of K-Means drop from 0.93 to 0.53: the Euclidean distance is dominated by the large coordinate and the other no longer counts. With StandardScaler the ARI goes back to 0.93. Before a distance-based clustering, the question about units of measurement must always be asked.

Hierarchical clustering

Agglomerative clustering starts with one group per point and at every step merges the two closest groups, until only one is left. The sequence of merges, with the distance at which they happen, is a tree called a dendrogram; cutting it at a certain height gives a partition, so you do not have to fix $k$ before seeing the tree. What changes is the meaning of «distance between two groups», the linkage:

• single: the minimum distance between a point of one and a point of the other, which chains nearby points and follows elongated shapes;
• complete: the maximum distance, which produces compact groups of limited diameter;
• average: the mean of all the pairwise distances;
• ward: the merge that makes the total inertia grow the least, therefore the closest in spirit to K-Means.

The four items are not cosmetic variants: on the same set of points they produce different trees, and the choice changes the shapes you manage to find (single) or to avoid (complete). The scikit-learn documentation presents Ward as similar to the K-Means objective function, and indeed in the previous table Ward and K-Means obtain similar results on the blobs.

The listing computes the merge heights of the four linkages on five numbers, 1, 2, 6, 7 and 15, which can be worked out by hand, and checks that scikit-learn and SciPy give the same partition.

# Hierarchical clustering: merges computed by hand on five numbers and compared with SciPy
import numpy as np
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import pdist
from sklearn.cluster import AgglomerativeClustering
from sklearn.datasets import make_blobs
from sklearn.metrics import adjusted_rand_score

x = np.array([[1.0], [2.0], [6.0], [7.0], [15.0]])
for method in ["single", "complete", "average", "ward"]:
    Z = linkage(x, method=method)           # rows: [cluster a, cluster b, height, size]
    print(f"{method:8s} heights of the 4 merges:", np.round(Z[:, 2], 3).tolist())

# Hand check of the heights:
#   single:   (1,2)=1  (6,7)=1  {1,2}+{6,7}: min = 6-2 = 4   +{15}: min = 15-7 = 8
#   complete: (1,2)=1  (6,7)=1  {1,2}+{6,7}: max = 7-1 = 6   +{15}: max = 15-1 = 14
#   average:  {1,2}+{6,7}: (5+6+4+5)/4 = 5   +{15}: (14+13+9+8)/4 = 11
#   ward:     height = sqrt(2 nA nB / (nA + nB)) * |centroid A - centroid B|
for nA, nB, cA, cB in [(2, 2, 1.5, 6.5), (4, 1, 4.0, 15.0)]:
    print("ward by hand:", round(float(np.sqrt(2 * nA * nB / (nA + nB)) * abs(cA - cB)), 3))

# Cutting the tree into 2 and 3 groups
Z = linkage(x, method="single")
for k in (2, 3):
    print(f"cut into {k} groups:", fcluster(Z, t=k, criterion="maxclust").tolist())

# sklearn and SciPy do the same thing: same partitions on 300 points
X, _ = make_blobs(n_samples=300, centers=4, cluster_std=1.0, random_state=2)
a = AgglomerativeClustering(n_clusters=4, linkage="ward").fit_predict(X)
b = fcluster(linkage(X, method="ward"), t=4, criterion="maxclust")
print("ARI between Agglomerative(ward) and linkage+fcluster:", adjusted_rand_score(a, b))

# The memory cost: the distance matrix has n(n-1)/2 entries of 8 bytes
for n in (1_000, 10_000, 100_000):
    entries = n * (n - 1) // 2
    print(f"n = {n:>6}: {entries:>11} distances, {entries * 8 / 1e6:>7.0f} MB")
print("pdist check on 5 points:", len(pdist(x)), "entries (5*4/2 = 10)")

# Output:
# single   heights of the 4 merges: [1.0, 1.0, 4.0, 8.0]
# complete heights of the 4 merges: [1.0, 1.0, 6.0, 14.0]
# average  heights of the 4 merges: [1.0, 1.0, 5.0, 11.0]
# ward     heights of the 4 merges: [1.0, 1.0, 7.071, 13.914]
# ward by hand: 7.071
# ward by hand: 13.914
# cut into 2 groups: [1, 1, 1, 1, 2]
# cut into 3 groups: [1, 1, 2, 2, 3]
# ARI between Agglomerative(ward) and linkage+fcluster: 1.0
# n =   1000:      499500 distances,       4 MB
# n =  10000:    49995000 distances,     400 MB
# n = 100000:  4999950000 distances,   40000 MB
# pdist check on 5 points: 10 entries (5*4/2 = 10)

The heights coincide with the hand computation written in the comments. With single linkage the merges are (1,2) and (6,7) at distance 1, then the two groups at distance 4 (from 2 to 6) and finally 15 at distance 8 (from 7 to 15); with complete they become 6 and 14, because they take the maximum distance.

For Ward the formula $\sqrt{2 n_A n_B / (n_A + n_B)}\,\|c_A - c_B\|$ holds, and the heights 7.071 and 13.914 computed by hand in the listing are equal to SciPy's. The cost is in the last part of the output: the pairwise distance matrix has $n(n-1)/2$ entries of 8 bytes, 400 MB at ten thousand points and 40,000 MB at a hundred thousand. That is why classic hierarchical clustering remains a tool for a few thousand points.

DBSCAN and HDBSCAN: density

DBSCAN (Ester, Kriegel, Sander and Xu, KDD 1996) does not look for centers, it looks for dense regions. A point is a core point if there are at least min_samples points within the radius eps; a cluster is the set of core points connected to one another, plus the border points close to a core; the rest is noise (label -1). You do not need $k$ and isolated points are declared noise instead of being forced into a group, but everything depends on eps.

HDBSCAN, available in scikit-learn as sklearn.cluster.HDBSCAN, extends the idea to several density scales and does not require eps. The listing varies eps on two moons of 500 points, with min_samples=5, and prints how many clusters and how much noise result; it also computes the quantiles of the distance from the fifth nearest point, which is the information of the «k-distance plot» without the plot.

# DBSCAN and the eps radius: on the two moons the result changes drastically with a single parameter
import numpy as np
from sklearn.cluster import DBSCAN, HDBSCAN
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score as ari

X, y = make_moons(n_samples=500, noise=0.07, random_state=0)
print(f"{'eps':>5s} {'clusters':>8s} {'noise points':>12s} {'ARI':>6s}")
for eps in [0.02, 0.05, 0.1, 0.2, 0.3, 0.5, 1.0]:
    lab = DBSCAN(eps=eps, min_samples=5).fit_predict(X)
    n_cl = len(set(lab)) - (1 if -1 in lab else 0)
    print(f"{eps:5.2f} {n_cl:8d} {(lab == -1).sum():12d} {ari(y, lab):6.2f}")

# The k-distance plot, without the plot: the distance to the 5th neighbor, sorted
from sklearn.neighbors import NearestNeighbors
d = NearestNeighbors(n_neighbors=5).fit(X).kneighbors(X)[0][:, -1]   # d[:, 0] is the point itself
q = np.quantile(d, [0.5, 0.9, 0.95, 0.99, 1.0])
print("distance to the 5th neighbor, quantiles 50/90/95/99/100%:", np.round(q, 3).tolist())

h = HDBSCAN(min_cluster_size=20, copy=True).fit(X)
print(f"HDBSCAN without eps: {h.labels_.max() + 1} clusters, {(h.labels_ == -1).sum()} noise, ARI = {ari(y, h.labels_):.2f}")

# Output:
#   eps clusters noise points    ARI
#  0.02        0          500   0.00
#  0.05       30          253   0.02
#  0.10        3           18   0.71
#  0.20        2            0   1.00
#  0.30        2            0   1.00
#  0.50        1            0   0.00
#  1.00        1            0   0.00
# distance to the 5th neighbor, quantiles 50/90/95/99/100%: [0.062, 0.1, 0.117, 0.148, 0.194]
# HDBSCAN without eps: 2 clusters, 5 noise, ARI = 0.98

With eps=0.02 no point has enough neighbors and all 500 are noise; with 0.05 DBSCAN breaks the moons into 30 fragments; with 0.10 it finds 3 clusters and 18 noise points (ARI 0.71); with 0.20 and 0.30 it finds the two moons (ARI 1.00); from 0.50 up the two moons merge into a single cluster.

The good interval is wide, but it has sharp edges. The quantiles help place it: 95% of the points have their fifth nearest point (the point itself included) within 0.117, so an eps slightly above that value is a good starting candidate, not a guaranteed value. HDBSCAN, without eps, finds 2 clusters, 5 noise points and ARI 0.98. When the densities are similar DBSCAN is convenient; when they are very different, or you do not want to search for eps, HDBSCAN.

Gaussian mixtures: elliptical groups and soft membership

A Gaussian mixture (Gaussian Mixture Model, GMM) describes the data as a sum of $k$ normal distributions, each with its own mean, covariance and weight. It is trained with the EM algorithm, which has the same two-step structure as Lloyd's but assigns every point a probability of belonging to each group instead of a single label.

With 'full' covariance the groups can be tilted ellipses, and the number of components is chosen with the BIC (Bayesian Information Criterion, lower is better), which penalizes the added parameters. The listing takes the three stripes of the previous comparison, standardized, and compares K-Means and GMM; then it chooses $k$ with the BIC and shows the probabilities of the most ambiguous point.

# Gaussian Mixture: elliptical clusters and probabilistic memberships, k chosen with the BIC
import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score as ari
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(0)                          # three horizontal stripes with staggered centers
cov = np.diag([9.0, 0.15])
X = np.vstack([rng.multivariate_normal(m, cov, 200) for m in [(0, 0), (2, 1.6), (4, 3.2)]])
y = np.repeat([0, 1, 2], 200)
X = StandardScaler().fit_transform(X)                    # as in the previous comparison

km = KMeans(3, n_init=10, random_state=0).fit(X)
gm = GaussianMixture(3, covariance_type="full", n_init=5, random_state=0).fit(X)
print(f"ARI K-Means = {ari(y, km.labels_):.2f}   ARI GMM full = {ari(y, gm.predict(X)):.2f}")

bic = [GaussianMixture(k, covariance_type="full", n_init=5, random_state=0).fit(X).bic(X) for k in range(1, 7)]
print("BIC for k = 1..6:", np.round(bic, 1).tolist(), "-> minimum at k =", int(np.argmin(bic)) + 1)

# The most ambiguous point: membership probabilities (soft) against a single label (hard)
P = gm.predict_proba(X)
i = int(np.argsort(P.max(axis=1))[0])
print("most uncertain point:", np.round(X[i], 2).tolist(), " probabilities:", np.round(P[i], 3).tolist(), " hard label:", int(gm.predict(X[[i]])[0]))
print("component weights:", np.round(gm.weights_, 3).tolist())

# Output:
# ARI K-Means = 0.54   ARI GMM full = 0.93
# BIC for k = 1..6: [3281.8, 3224.3, 3051.9, 3089.8, 3115.4, 3154.3] -> minimum at k = 3
# most uncertain point: [-0.25, -0.58]  probabilities: [0.0, 0.479, 0.521]  hard label: 2
# component weights: [0.329, 0.337, 0.334]

On the three stripes K-Means stays at 0.54 and the mixture rises to 0.93, because the full covariance describes the elongated shape. The BIC has its minimum at $k=3$ (3051.9, against 3224.3 for $k=2$ and 3089.8 for $k=4$), that is the true number.

The most uncertain point has probabilities 0.479 and 0.521 for two components: a hard label assigns it to the second, without saying that the margin is four percentage points. The mixture costs more, and with few points per component the full covariance can become unstable: in that case covariance_type='diag' or 'tied' reduces the parameters.

A summary of when to use each method, drawn from the previous experiments and not from a general rule:

• K-Means or Ward, if the groups are compact, similar in size, and the features are scaled;
• Gaussian mixture, if the groups are elliptical or you need membership probabilities;
• DBSCAN or HDBSCAN, if the groups have any shape and you want to discard isolated points as noise;
• hierarchical clustering, if you want to see the tree of merges and you have at most a few thousand points.

No item is a verdict: the ARI table shows that with different data each method falls to zero. The criterion that counts is whether the assumptions of the method (compactness, shape, uniform density) are reasonable for your data, and to check that you have to look at them, for example projected onto two dimensions with the PCA of the next article.

Worked exercise

Run K-Means by hand with $k=2$ on the numbers 1, 3, 4, 10, 12, with initial centers 1 and 3. For every iteration write the labels, the new centers and the value of $J$. After how many iterations do you stop and what is the final inertia? Then check with the listing.

# Exercise: K-Means (k=2) on 1, 3, 4, 10, 12 with initial centers 1 and 3
import numpy as np

x = np.array([1.0, 3.0, 4.0, 10.0, 12.0])
c = np.array([1.0, 3.0])
prev = None
for it in range(1, 10):
    lab = np.abs(x[:, None] - c[None, :]).argmin(axis=1)
    J = ((x - c[lab]) ** 2).sum()
    print(f"iter {it}: centers {np.round(c, 3).tolist()}  labels {lab.tolist()}  J = {J:.3f}")
    if prev is not None and (lab == prev).all():
        break
    c = np.array([x[lab == k].mean() for k in range(2)])
    prev = lab

# Output:
# iter 1: centers [1.0, 3.0]  labels [0, 1, 1, 1, 1]  J = 131.000
# iter 2: centers [1.0, 7.25]  labels [0, 0, 0, 1, 1]  J = 43.125
# iter 3: centers [2.667, 11.0]  labels [0, 0, 0, 1, 1]  J = 6.667

What you can do now

After this article you can do five things, all with a numerical check at hand:

• write Lloyd's algorithm and show that the inertia never grows;
• choose the initialization and n_init knowing how much they weigh (57% less inertia on a measured dataset);
• choose $k$ by crossing elbow, silhouette, Davies-Bouldin and gap statistic, and explain why they can disagree;
• recognize the data on which K-Means fails and try hierarchical, DBSCAN or Gaussian mixture;
• check the units of measurement before any distance-based method.

The five abilities revolve around one idea: clustering methods make assumptions about the shape of the groups and do not tell you, so it is up to you to test them on data where the answer is known. In the next article you need the same discipline to reduce dimensions: PCA. To evaluate a supervised model with the right metrics there is instead the article on metrics; another take on clustering is the one on K-Means and DBSCAN without labels, in the Machine Learning Engineer series. Neural networks too are trained by improving a loss function step by step, but with the gradient in place of Lloyd's alternation; and a clustering in use must be monitored like any model, as the article on deployment shows.