Scikit-Learn Clustering

Clustering groups rows without labels: customer segments, document topics, anomalous events. Every example here ran with scikit-learn 1.9.1 on generated data, so the true groups are known and each algorithm’s answer can be checked.

K-Means and Choosing K

K-Means needs the number of clusters up front. The silhouette score compares candidates.

from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import silhouette_score

X, _ = make_blobs(n_samples=600, centers=4, cluster_std=1.0, random_state=7)

for k in range(2, 8):
    km = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    print(f"K={k}  silhouette {silhouette_score(X, km.labels_):.3f}"
          f"  inertia {km.inertia_:8.0f}")

$ python choose_k.py
K=2  silhouette 0.588  inertia    18363
K=3  silhouette 0.762  inertia     5241
K=4  silhouette 0.811  inertia     1127
K=5  silhouette 0.680  inertia     1020
K=6  silhouette 0.577  inertia      921
K=7  silhouette 0.446  inertia      814

Silhouette peaks at K=4, the true number of blobs. Inertia, the total squared distance to the nearest centroid, always falls as K grows, so it can’t pick K by itself; the useful signal is the elbow, where the big drop from 5241 to 1127 flattens into small steps after K=4.

Scaling Changes the Clusters

K-Means measures Euclidean distance, so a feature in large units takes over. Here monthly spend separates three real groups, while income is unrelated noise measured in the tens of thousands.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(0)
y_true = np.repeat([0, 1, 2], 200)
spend = np.r_[rng.normal(20, 3, 200), rng.normal(50, 3, 200), rng.normal(80, 3, 200)]
income = rng.normal(60_000, 15_000, 600)              # no relation to the groups
X = np.column_stack([spend, income])

for name, data in [("raw", X), ("scaled", StandardScaler().fit_transform(X))]:
    labels = KMeans(n_clusters=3, n_init=10, random_state=0).fit_predict(data)
    print(f"{name:6s} agreement with true groups (ARI): "
          f"{adjusted_rand_score(y_true, labels):.3f}")

$ python kmeans_scale.py
raw    agreement with true groups (ARI): 0.002
scaled agreement with true groups (ARI): 0.476

On raw data the clusters follow income, whose spread of about 15,000 dwarfs spend’s gaps of 30, so the result has almost no relation to the real groups (ARI near 0; 1.0 is a perfect match). Standardising lets spend count, but the recovery is only partial: income still gets an equal vote even though it is pure noise. Scaling fixes units; it can’t tell the algorithm which features matter.

When K-Means Fails: Non-Round Shapes

Two interleaved half-moons are obvious to the eye but not round. DBSCAN groups points by density instead of distance to a centre.

from sklearn.cluster import DBSCAN, KMeans
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score
from sklearn.preprocessing import StandardScaler

X, y_true = make_moons(n_samples=500, noise=0.07, random_state=0)
X = StandardScaler().fit_transform(X)

km = KMeans(n_clusters=2, n_init=10, random_state=0).fit_predict(X)
db = DBSCAN(eps=0.3, min_samples=5).fit_predict(X)

print(f"K-Means ARI: {adjusted_rand_score(y_true, km):.3f}")
print(f"DBSCAN  ARI: {adjusted_rand_score(y_true, db):.3f}")
print(f"K-Means accuracy: {max((km == y_true).mean(), (km != y_true).mean()):.3f}")
print(f"DBSCAN found {len(set(db) - {-1})} clusters, {(db == -1).sum()} noise points")

$ python moons.py
K-Means ARI: 0.478
DBSCAN  ARI: 1.000
K-Means accuracy: 0.846
DBSCAN found 2 clusters, 0 noise points

K-Means cuts the moons with a straight boundary between two centroids, so the tips of each moon land in the wrong cluster: about 1 point in 7 is misassigned, and the ARI, which penalises such mixing heavily, falls below 0.5. DBSCAN follows chains of nearby points, so each curved moon becomes one cluster and the match is perfect. Points too isolated to belong anywhere would get the label -1; here none did.

DBSCAN’s eps Is the Whole Game

DBSCAN doesn’t ask for K, but its neighbourhood radius eps decides everything.

from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons
from sklearn.preprocessing import StandardScaler

X, _ = make_moons(n_samples=500, noise=0.07, random_state=0)
X = StandardScaler().fit_transform(X)

for eps in [0.05, 0.1, 0.2, 0.3, 0.5, 1.0]:
    labels = DBSCAN(eps=eps, min_samples=5).fit_predict(X)
    print(f"eps={eps:<4}  clusters {len(set(labels) - {-1}):3d}"
          f"  noise {(labels == -1).sum():3d}")

$ python eps_sweep.py
eps=0.05  clusters   8  noise 454
eps=0.1   clusters  25  noise 123
eps=0.2   clusters   2  noise   3
eps=0.3   clusters   2  noise   0
eps=0.5   clusters   2  noise   0
eps=1.0   clusters   1  noise   0

Too small an eps and most points have too few neighbours, so they become noise and the rest shatter into many fragments (8 at 0.05, 25 at 0.1). Too large and the two moons merge into one. Only a middle band, here 0.2 to 0.5, recovers the true structure, and it has to be found by experiment for each dataset.

Agglomerative Clustering and Linkage

Agglomerative clustering starts with every point alone and repeatedly merges the closest pair of clusters. The linkage rule defines “closest”.

from sklearn.cluster import AgglomerativeClustering
from sklearn.datasets import make_blobs, make_moons
from sklearn.metrics import adjusted_rand_score
from sklearn.preprocessing import StandardScaler

blobs, yb = make_blobs(n_samples=500, centers=3, cluster_std=1.2, random_state=1)
moons, ym = make_moons(n_samples=500, noise=0.05, random_state=0)

for linkage in ["ward", "complete", "average", "single"]:
    scores = []
    for X, y in [(blobs, yb), (moons, ym)]:
        X = StandardScaler().fit_transform(X)
        labels = AgglomerativeClustering(n_clusters=len(set(y)), linkage=linkage).fit_predict(X)
        scores.append(adjusted_rand_score(y, labels))
    print(f"{linkage:8s} blobs ARI {scores[0]:.3f}   moons ARI {scores[1]:.3f}")

$ python linkage.py
ward     blobs ARI 0.936   moons ARI 0.467
complete blobs ARI 0.930   moons ARI 0.261
average  blobs ARI 0.882   moons ARI 0.583
single   blobs ARI 0.570   moons ARI 1.000

No linkage wins everywhere. Ward, which merges the pair that least increases within-cluster variance, is best on round blobs and poor on moons. Single linkage, which merges by the closest pair of points, follows the moons perfectly but chains blobs together. Choose the linkage to match the shape you expect.

K-Means Always Finds Clusters

K-Means returns K groups no matter what. Compare its silhouette on real clusters with its silhouette on uniform noise.

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

rng = np.random.default_rng(0)
clustered, _ = make_blobs(n_samples=600, centers=4, random_state=7)
noise = rng.uniform(clustered.min(), clustered.max(), size=clustered.shape)

for name, X in [("real clusters", clustered), ("uniform noise", noise)]:
    labels = KMeans(n_clusters=4, n_init=10, random_state=0).fit_predict(X)
    sizes = np.bincount(labels).tolist()
    print(f"{name:14s} sizes {sizes}  silhouette {silhouette_score(X, labels):.3f}")

$ python noise_clusters.py
real clusters  sizes [150, 150, 150, 150]  silhouette 0.811
uniform noise  sizes [174, 142, 129, 155]  silhouette 0.410

On uniform noise K-Means still returns four tidy groups of similar size, and nothing in the labels warns you. The silhouette does: 0.41 against 0.81 for real clusters. Compare against a noise baseline like this before presenting segments as real.

Practice

import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(1)
df = pd.DataFrame({
    "spend": np.r_[rng.normal(2000, 300, 100), rng.normal(300, 80, 300)],
    "orders": np.r_[rng.poisson(20, 100), rng.poisson(3, 300)],
    "days_inactive": np.r_[rng.exponential(10, 100), rng.exponential(60, 300)],
})
X = StandardScaler().fit_transform(df)
df["segment"] = KMeans(n_clusters=2, n_init=10, random_state=0).fit_predict(X)
print(df.groupby("segment").mean().round(1).assign(size=df["segment"].value_counts()))

$ python p1_profile.py
          spend  orders  days_inactive  size
segment                                     
0         293.2     3.1           57.1   300
1        1977.9    19.9            8.0   100

import numpy as np
from sklearn.cluster import DBSCAN

rng = np.random.default_rng(0)
normal = rng.normal(0, 1, size=(300, 2))
outliers = np.array([[6.0, 6.0], [-7.0, 5.0], [8.0, -6.0]])
X = np.vstack([normal, outliers])

labels = DBSCAN(eps=1.0, min_samples=5).fit_predict(X)
print("rows labelled noise:", np.where(labels == -1)[0].tolist())

$ python p2_outliers.py
rows labelled noise: [151, 239, 300, 301, 302]

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score

rng = np.random.default_rng(0)
y_true = np.repeat([0, 1, 2], 200)
spend = np.r_[rng.normal(20, 3, 200), rng.normal(50, 3, 200), rng.normal(80, 3, 200)]

labels = KMeans(n_clusters=3, n_init=10, random_state=0).fit_predict(spend.reshape(-1, 1))
print(f"spend only ARI: {adjusted_rand_score(y_true, labels):.3f}")

$ python p3_drop_noise.py
spend only ARI: 1.000