Unsupervised learning has no labels, which means there is no score to optimise and no test set to check against. That makes it far easier to produce a confident, meaningless result - so the discipline is different, and mostly consists of validating the output against something outside the algorithm.
k-means
k-means partitions the data into k groups by minimising the distance from each point to its group's centre.
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(3)
# Three genuine groups.
groups = [
rng.normal([0, 0], 0.8, (300, 2)),
rng.normal([5, 5], 0.8, (300, 2)),
rng.normal([0, 6], 0.8, (300, 2)),
]
X = np.vstack(groups)
truth = np.repeat([0, 1, 2], 300)
model = KMeans(n_clusters=3, n_init=10, random_state=0).fit(X)
print('cluster centres:')
print(np.round(model.cluster_centers_, 2))
print()
print('cluster sizes:', np.bincount(model.labels_))
print(f'inertia (total squared distance to centres): {model.inertia_:.0f}')
print()
# The labels are arbitrary, so compare the partition rather than the ids.
crosstab = pd.crosstab(truth, model.labels_)
print(crosstab.to_string())
print()
print(f'each true group maps to one cluster: '
f'{bool((crosstab > 0).sum(axis=1).max() == 1)}')
Open this example in the editor and change it. Nothing you type is saved.
Three things k-means assumes, and each one is a way it goes wrong:
Scale matters. Distance is computed across all features, so a feature in pounds dominates one in years.
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(5)
n = 600
frame = pd.DataFrame({
'tenure_years': np.r_[rng.normal(2, 0.6, n // 2), rng.normal(8, 0.6, n // 2)],
'annual_spend': np.r_[rng.normal(50_000, 9_000, n // 2),
rng.normal(52_000, 9_000, n // 2)],
})
truth = np.repeat([0, 1], n // 2)
for label, data in [('raw', frame.to_numpy()),
('standardised', StandardScaler().fit_transform(frame))]:
labels = KMeans(n_clusters=2, n_init=10, random_state=0).fit_predict(data)
agreement = max((labels == truth).mean(), (labels != truth).mean())
print(f'{label:<14} recovers the true split {agreement:.1%} of the time')
print()
print('Unstandardised, the clustering is entirely about spend, which')
print('varies by thousands, and ignores tenure, which varies by years.')
Open this example in the editor and change it. Nothing you type is saved.
Clusters are assumed round and similar in size. k-means draws straight boundaries at equal distance between centres, so it cannot find an elongated or a curved group.
import numpy as np
from sklearn.cluster import KMeans, DBSCAN
from sklearn.datasets import make_moons
X, truth = make_moons(n_samples=600, noise=0.06, random_state=0)
for name, labels in [
('k-means', KMeans(n_clusters=2, n_init=10, random_state=0).fit_predict(X)),
('DBSCAN', DBSCAN(eps=0.2, min_samples=5).fit_predict(X)),
]:
agreement = max((labels == truth).mean(), (labels != truth).mean())
print(f'{name:<10} recovers the two crescents {agreement:.1%} of the time')
print()
print('k-means cuts a straight line through both crescents. DBSCAN')
print('follows density instead, and gets them exactly.')
Open this example in the editor and change it. Nothing you type is saved.
k must be chosen. The algorithm will happily produce whatever number you ask for.
Choosing k, and being honest
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
rng = np.random.default_rng(7)
X = np.vstack([rng.normal([0, 0], 0.8, (250, 2)),
rng.normal([5, 5], 0.8, (250, 2)),
rng.normal([0, 6], 0.8, (250, 2))])
rows = []
for k in range(2, 9):
model = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
rows.append({'k': k, 'inertia': round(model.inertia_),
'silhouette': round(silhouette_score(X, model.labels_), 3)})
table = pd.DataFrame(rows)
table['inertia_drop'] = (-table['inertia'].diff()).fillna(0).astype(int)
print(table.to_string(index=False))
print()
print(f'the true k is 3')
print(f'silhouette picks k={int(table.loc[table["silhouette"].idxmax(), "k"])}')
print('inertia falls monotonically and always will - it cannot choose k')
Open this example in the editor and change it. Nothing you type is saved.
Inertia always falls as k rises, because more centres means shorter distances; at k = n it reaches zero. The "elbow" is a visual judgement and often ambiguous.
Silhouette measures how much closer each point is to its own cluster than to the next nearest, and it does have a maximum. It is the better default.
But the honest test is external:
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
rng = np.random.default_rng(11)
# NO clusters at all: one uniform blob.
X = rng.normal(0, 1, (900, 4))
rows = []
for k in range(2, 8):
labels = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(X)
rows.append({'k': k, 'silhouette': round(silhouette_score(X, labels), 3),
'sizes': np.bincount(labels).tolist()})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('There is no structure here, and k-means still returns k tidy')
print('groups with a positive silhouette for every k. The algorithm')
print('cannot tell you that there is nothing to find.')
Open this example in the editor and change it. Nothing you type is saved.
So the validation has to come from outside the clustering:
- Do the clusters differ on a variable that was not used to build them? If a "high-value" segment does not actually spend more on a held-out measure, it is not a segment.
- Are they stable? Re-cluster on a random half of the data and see whether the same groups appear.
- Can somebody name them? A cluster nobody can describe in a sentence will not be used.
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(13)
n = 1200
# Three real segments, each with its own behaviour.
segment = rng.choice([0, 1, 2], n, p=[.5, .3, .2])
frame = pd.DataFrame({
'sessions': rng.poisson([4, 14, 28][0] if False else np.choose(segment, [4, 14, 28])),
'lessons': rng.poisson(np.choose(segment, [1, 6, 18])),
'days_active': rng.poisson(np.choose(segment, [2, 9, 22])),
})
# A variable NOT used for clustering, held back for validation.
frame['paid'] = rng.random(n) < np.choose(segment, [0.03, 0.18, 0.52])
CLUSTER_ON = ['sessions', 'lessons', 'days_active']
labels = KMeans(n_clusters=3, n_init=10, random_state=0).fit_predict(
StandardScaler().fit_transform(frame[CLUSTER_ON]))
frame['cluster'] = labels
profile = frame.groupby('cluster').agg(
n=('paid', 'size'), sessions=('sessions', 'mean'),
lessons=('lessons', 'mean'), days=('days_active', 'mean'),
paid_rate=('paid', 'mean')).round(2).sort_values('lessons')
print(profile.to_string())
print()
spread = profile['paid_rate'].max() - profile['paid_rate'].min()
print(f'the held-out conversion rate ranges {spread:.1%} across clusters')
print('That is the evidence the segments are real: they differ on a')
print('variable the clustering never saw.')
Open this example in the editor and change it. Nothing you type is saved.
Stability
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score
rng = np.random.default_rng(17)
real = np.vstack([rng.normal([0, 0], 0.7, (300, 2)),
rng.normal([4, 4], 0.7, (300, 2)),
rng.normal([0, 5], 0.7, (300, 2))])
noise = rng.normal(0, 1, (900, 2))
for name, X in [('three real clusters', real), ('no structure', noise)]:
agreements = []
for seed in range(10):
half = rng.permutation(len(X))[:len(X) // 2]
a = KMeans(3, n_init=10, random_state=0).fit(X[half])
b = KMeans(3, n_init=10, random_state=seed + 1).fit(X)
agreements.append(adjusted_rand_score(a.predict(X), b.labels_))
print(f'{name:<22} stability {np.mean(agreements):.3f} '
f'(1.0 = identical partitions)')
print()
print('Real structure reproduces across resamples. Noise does not.')
print('This is the single most useful check in clustering.')
Open this example in the editor and change it. Nothing you type is saved.
PCA
Principal component analysis rotates the data so that the first axis carries the most variance, the second the most of what remains, and so on. Keeping the first few is lossy compression.
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(19)
n = 1000
# Twelve columns generated from three underlying factors.
factors = rng.normal(0, 1, (n, 3))
loadings = rng.normal(0, 1, (3, 12))
X = factors @ loadings + rng.normal(0, 0.25, (n, 12))
scaled = StandardScaler().fit_transform(X)
pca = PCA().fit(scaled)
explained = pd.DataFrame({
'component': range(1, 13),
'variance': pca.explained_variance_ratio_.round(3),
'cumulative': pca.explained_variance_ratio_.cumsum().round(3),
})
print(explained.head(6).to_string(index=False))
print()
print(f'three components carry '
f'{pca.explained_variance_ratio_[:3].sum():.1%} of the variance')
print('- which is right, because three factors generated the data')
print()
for target in [0.80, 0.90, 0.95, 0.99]:
needed = int(np.searchsorted(pca.explained_variance_ratio_.cumsum(), target) + 1)
print(f'{target:.0%} of the variance needs {needed} of 12 components')
Open this example in the editor and change it. Nothing you type is saved.
Three rules for PCA:
Standardise first, always. PCA maximises variance, so an unscaled feature with a large range becomes the first component by arithmetic rather than meaning.
Components are not features. The first component is a weighted mixture of all the originals; it rarely has a name, and interpreting it as "customer value" is a story rather than a finding.
Fit on the training set only. PCA is a transformation that is learned, so it belongs inside the pipeline.
import numpy as np
from sklearn.decomposition import PCA
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(23)
n = 2000
factors = rng.normal(0, 1, (n, 3))
X = factors @ rng.normal(0, 1, (3, 40)) + rng.normal(0, 0.4, (n, 40))
y = ((factors[:, 0] + 0.6 * factors[:, 1] + rng.normal(0, 0.6, n)) > 0).astype(int)
folds = StratifiedKFold(5, shuffle=True, random_state=0)
plain = cross_val_score(HistGradientBoostingClassifier(random_state=0), X, y,
cv=folds, scoring='roc_auc').mean()
reduced = cross_val_score(
Pipeline([('scale', StandardScaler()), ('pca', PCA(n_components=4)),
('model', HistGradientBoostingClassifier(random_state=0))]),
X, y, cv=folds, scoring='roc_auc').mean()
print(f'40 raw features: AUC {plain:.3f}')
print(f'4 principal components: AUC {reduced:.3f}')
print()
print('Ten times fewer features, no loss - because the 40 columns only')
print('ever carried three factors of information.')
Open this example in the editor and change it. Nothing you type is saved.
The warning
import pandas as pd
print(pd.DataFrame([
{'temptation': 'the algorithm found 5 segments',
'reality': 'you asked for 5; it would have found 9 just as happily'},
{'temptation': 'this cluster is our premium customers',
'reality': 'name it only if a held-out variable agrees'},
{'temptation': 'the segments are stable',
'reality': 'measure it: re-cluster a random half and compare'},
{'temptation': 'PC1 represents engagement',
'reality': 'PC1 is a weighted sum of everything; it means nothing on its own'},
{'temptation': 'we will act on these segments',
'reality': 'decide the action first, then see whether the split supports it'},
]).to_string(index=False))
Open this example in the editor and change it. Nothing you type is saved.
The honest framing: clustering is a hypothesis generator, not a result. It suggests groups; something else has to confirm that the groups are real and that they matter.
When it goes wrong
| Symptom | Cause |
|---|---|
| One feature dominates the clustering | Unstandardised; scale first |
| Clusters found in random data | Expected; k-means always returns k groups |
| The clusters change every run | No real structure; measure stability |
| k-means missed an obvious shape | It assumes round clusters; try DBSCAN |
| The elbow is ambiguous | Inertia always falls; use silhouette |
| A segment nobody can describe | It will not be used; validate externally |
| PCA components are uninterpretable | They are mixtures; that is normal |
| PCA hurt the model | Fitted on all the data, or too few components |
A check you can run
Cluster a random half of your data, cluster the whole of it, and compute adjusted_rand_score between the two labellings on the overlapping rows.
Above about 0.7 the structure is real. Near zero it is not, and the segmentation deck built on it is describing noise. That is five lines, and it is the check that almost nobody runs.