Unsupervised learning

Sometimes there are no labels: no column saying which customers churned, no file saying which documents matter. This section covers what can still be learned from the shape of the data alone, grouping similar things together, squeezing many dimensions down to a few you can plot, and spotting the examples unlike everything else. It also covers the honest limit, which is that an algorithm will always return groups, whether or not any exist.

Clustering#

K-means partitions data into k groups by repeatedly assigning points to the nearest centre and moving the centres. It is fast and assumes clusters are roughly spherical and similarly sized, which is often not true.

The algorithm is four steps and worth knowing exactly, because every one of its failure modes is visible in them. Place k centres at random points. Assign each data point to the nearest centre. Move each centre to the mean of the points assigned to it. Repeat the last two until nothing changes. Note what is missing: nothing anywhere in that loop asks whether the groups it is producing are real. It will terminate, and it will return k groups, on any data you give it.

from sklearn.cluster import KMeans

km = KMeans(n_clusters=5, n_init=10).fit(X)
km.labels_           # which cluster each point belongs to
km.inertia_          # total within-cluster distance, lower is tighter

Watch it invent clusters that are not there#

That note is the most important sentence in this section and it should not be taken on trust. Generate 500 points spread uniformly at random over a square — data with no groups in it by construction, the flattest possible null case — and run the whole standard procedure on it.

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

rng = np.random.default_rng(0)
noise = rng.uniform(0, 1, (500, 2))                     # no structure at all
blobs, _ = make_blobs(n_samples=500, centers=4, cluster_std=0.6, random_state=0)

for name, data in (("uniform", noise), ("blobs", blobs)):
    sil = []
    for k in range(2, 9):
        labels = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(data)
        sil.append(round(silhouette_score(data, labels), 3))
    print(name, sil)

# uniform [0.359, 0.384, 0.406, 0.385, 0.392, 0.389, 0.392]
# blobs   [0.556, 0.597, 0.684, 0.601, 0.498, 0.43,  0.313]

print(np.bincount(KMeans(n_clusters=4, n_init=10, random_state=0).fit_predict(noise)))
# [141 109 136 114]   four tidy, balanced groups, out of nothing

The elbow behaves no better. Inertia on the uniform data falls by factors of 1.62, 1.59 and 1.54 for the first three splits and then by only 1.22 — that bend is exactly what people point at and call an elbow at k equal to 4. On the real four-cluster data the same ratios are 2.41, 2.20 and 2.59, then 1.11. The difference is not whether there is a bend, it is how sharp the bend is, and reading sharpness off a plot by eye is how a whole team ends up with four customer segments that do not exist.

The advice usually offered next is to rerun with a different seed and check the clusters are stable. Check that too, because it also fails. Comparing runs from five different seeds with the adjusted Rand index gives 0.93 to 1.00 on the uniform data and exactly 1.00 on the real clusters. The noise partitions stably, because the partition is a property of the square's geometry rather than of any structure in the points.

What does work is comparison against a null. Generate random data with the same shape and range as yours, cluster it the same way, and compare the scores: 0.406 against 0.684 here. If your data's silhouette is not clearly above what pure noise of the same shape achieves, you have found nothing. That comparison, formalised, is the gap statistic, and it is the only one of these methods that asks the right question.

DBSCAN is worth knowing as the alternative: it finds clusters of arbitrary shape, decides the number itself, and labels sparse points as noise rather than forcing them into a group.

What the silhouette number actually is#

It is worth knowing what you are reading rather than treating it as a score out of one. For each point, take a, the average distance to the other points in its own cluster, and b, the average distance to the points in the nearest other cluster. The point's silhouette is b minus a, divided by whichever of the two is larger. It is 1 when the point's own cluster is far closer than any other, 0 when the two are equally far, and negative when the point is actually nearer to another cluster than its own. The reported score is the average over all points.

from sklearn.metrics import silhouette_score

points = np.array([[0.0, 0.0], [0.0, 1.0], [5.0, 0.0], [5.0, 1.0]])
pair_labels = np.array([0, 0, 1, 1])

# By hand, for the first point:
#   a = distance to its one neighbour [0, 1]                    = 1.0
#   b = mean distance to [5, 0] and [5, 1] = (5 + sqrt(26)) / 2 = 5.0495
#   silhouette = (b - a) / max(a, b) = 4.0495 / 5.0495          = 0.802
print(silhouette_score(points, pair_labels))       # 0.8019609728144304

That last point matters more than the number. Cluster the two interleaved crescents of the classic two-moons dataset correctly and the silhouette is low, because the crescents are long and the ends of one are closer to the other crescent than to their own far end. A method scoring well on silhouette is a method producing round blobs, which is a statement about k-means rather than about your data.

Dimensionality reduction#

• PCA: finds the directions of greatest variance and projects onto them. Linear, fast, and the projection is interpretable as a combination of features.
• t-SNE and UMAP: non-linear, for visualising high-dimensional data in two dimensions. Good at showing local structure.

Anomaly detection#

The third use of unlabelled data is finding what does not fit. It is the same machinery pointed the other way: instead of describing the groups, you describe what normal looks like and flag whatever sits far from it. The useful framing is that an anomaly is rarely an extreme value in one column. It is an unusual combination — an ordinary amount at an unusual hour with an unusual number of items — which is why a rule on a single column misses almost everything worth catching.

Two cautions apply from the start. Unusual is not the same as wrong: a flag is a request for someone to look, never a verdict, and presenting it as one is how these systems lose the trust of the people who have to act on them. And because nobody can review a thousand flags, the number to measure is how many of the top twenty were worth a human's time, which is precision at k rather than any overall rate.

Where this comes back#

These unsupervised techniques reappear constantly in LLM work under different names. Clustering embeddings to discover topics in a corpus. Reducing dimensions to visualise what your retrieval index looks like. Anomaly detection to find eval cases that are unlike everything else.

classical name            what it is called in LLM work
clustering                topic discovery over a month of user questions
clustering                grouping eval failures by kind so you fix a cause, not a case
dimensionality reduction  the two-dimensional plot of a retrieval index
near-duplicate detection  deduplicating pretraining data, and spotting eval leakage
anomaly detection         finding the odd inputs your eval set never anticipated
anomaly detection         drift monitoring: today's inputs against the validated ones

Some of what you need to find is too rare and too varied to collect examples of: fraud, faults, corrupted records, the eval case unlike anything else in the set. You cannot train a classifier on them, so instead you learn what normal looks like and flag whatever departs from it. This section covers how that is done, then turns the same unsupervised tools on embeddings, where clustering a corpus and reading a handful of items from each cluster is the quickest honest answer to what is in all this text.

Finding the unusual#

Anomaly detection flags data unlike the rest: fraud, faults, corrupted records, eval cases that do not resemble anything else. Because anomalies are rare and varied, you usually cannot train a classifier on them and must learn what normal looks like instead.

• Statistical: flag points far from the mean in standard deviations. Simple, assumes roughly normal data.
• Isolation forest: anomalies are easier to isolate with random splits, so they end up in short paths. Fast and effective on tabular data.
• Distance-based: points far from their nearest neighbours are unusual. Works well on embeddings.

Those three sentences hide an enormous difference in how well they work, so run them against each other on a problem where the right answer is known. A thousand ordinary transactions with a skewed amount, an hour between 8 and 20, and a handful of items; then twenty planted anomalies whose amounts are perfectly ordinary but which happen at three in the morning with thirty items in the basket. The only honest measure without labels in production is precision at k, so score the top twenty flagged by each method.

import numpy as np
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(0)
n = 1000
normal = np.c_[rng.lognormal(4.0, 1.0, n),                   # amount: skewed
               rng.integers(8, 20, n).astype(float),         # hour: daytime
               rng.poisson(3, n).astype(float) + 1]          # items
planted = np.c_[rng.lognormal(4.0, 1.0, 20),                 # ordinary amounts
                rng.integers(2, 5, 20).astype(float),        # the middle of the night
                rng.poisson(30, 20).astype(float) + 1]       # far too many items

X = np.vstack([normal, planted])
y = np.r_[np.zeros(n), np.ones(20)]

def precision_at_20(scores):
    top = np.argsort(-scores)[:20]
    return y[top].sum() / 20

amount_z = np.abs((X[:, 0] - X[:, 0].mean()) / X[:, 0].std())
all_z = np.abs((X - X.mean(0)) / X.std(0)).max(axis=1)
median = np.median(X, axis=0)
mad = np.median(np.abs(X - median), axis=0)
robust_z = np.abs((X - median) / (1.4826 * mad)).max(axis=1)

logged = np.c_[np.log(X[:, 0]), X[:, 1], X[:, 2]]
forest = -IsolationForest(random_state=0).fit(logged).score_samples(logged)
scaled = StandardScaler().fit_transform(logged)
distances, _ = NearestNeighbors(n_neighbors=11).fit(scaled).kneighbors(scaled)
knn = distances[:, 1:].mean(axis=1)

for name, s in (("z on amount", amount_z), ("max z, all columns", all_z),
                ("robust z (median/MAD)", robust_z), ("isolation forest", forest),
                ("10-NN distance", knn)):
    print(name, precision_at_20(s))

# z on amount            0.0
# max z, all columns     0.75
# robust z (median/MAD)  0.8
# isolation forest       0.95
# 10-NN distance         0.65

Three lessons come out of that table and they transfer to any anomaly problem. A z-score assumes a symmetric distribution, so on skewed data such as money, latency or counts it reports the fat tail as anomalous and ignores everything else — take logs first or use the median and the median absolute deviation instead of the mean and the standard deviation. An anomaly is usually a combination rather than an extreme value in one column, which is why the multivariate methods win here: nothing about a 3am transaction or a 30-item basket is individually impossible. And the scale of the input still matters: log the skewed column and the isolation forest goes from 0.85 to 0.95.

Topics from embeddings#

Embed a corpus, cluster the embeddings, and read a sample from each cluster to name it. It is a fast way to understand what a large collection of text is about, and a practical use of these methods in LLM work.

from sklearn.cluster import KMeans

vectors = embed(documents)
labels = KMeans(n_clusters=12, n_init=10).fit_predict(vectors)

for k in range(12):
    sample = [d for d, l in zip(documents, labels) if l == k][:5]
    print(k, sample)          # read these and name the cluster yourself

One detail decides whether that works. k-means measures straight-line distance, while embedding similarity is cosine, which ignores length. The two agree exactly once the vectors are normalised to length 1: for unit vectors the squared distance is 2 minus twice the cosine, so ranking by one is ranking by the other. If your embeddings are not normalised, normalise them before clustering, or you are grouping by how long the documents are.

a = np.array([3.0, 4.0])
b = np.array([1.0, 0.0])

ua, ub = a / np.linalg.norm(a), b / np.linalg.norm(b)
cosine = ua @ ub
print(cosine)                                      # 0.6
print(((ua - ub) ** 2).sum(), 2 - 2 * cosine)      # 0.8000000000000002 0.8

The other decision is how many clusters, and here the honest answer is different from the one in the tabular case. You are not trying to recover a true number of topics; there is no such number. You are choosing a resolution for reading. Twelve clusters over ten thousand support tickets is a summary a person can read in an afternoon; two hundred is a second dataset. Pick the number by what you can act on, say so, and check the result by reading rather than by score.

There is one way to get a real number, and it is worth the hour it costs. Label a sample by hand — two hundred items is usually enough — and use it as a test set for whatever the unsupervised method was supposed to do. For an anomaly detector, that is precision at k on the flags a person actually reviewed. For clustering, take pairs of items and record whether a person says they belong together, then measure how often the clustering agrees. The point is not to label the whole corpus; it is to stop reporting a method's confidence in itself as though it were evidence.

The alternative number people report instead is the silhouette or the anomaly score, and both are internal: they say the algorithm did what it was designed to do on your data, which it always does. An internal score falling sharply when you change the input is useful information. An internal score on its own is not a result.

Clustering groups records that resemble each other when nobody has labelled anything. It is genuinely useful and easy to misuse, because the algorithm will hand you groups whether or not any exist. This section covers running k-means, why scaling the features first is not optional, choosing the number of groups with evidence rather than by eye, a density-based method for data that does not come in tidy blobs, and judging the result by reading what landed in each group.

k-means#

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

X_scaled = StandardScaler().fit_transform(X)       # otherwise the largest-valued feature dominates
for k in range(2, 11):
    labels = KMeans(n_clusters=k, n_init="auto", random_state=0).fit_predict(X_scaled)
    print(k, round(silhouette_score(X_scaled, labels), 3))

k-means places k centres and assigns each point to the nearest, then moves each centre to the middle of its points, repeating until stable. It assumes roughly round clusters of similar size. The silhouette score compares how close points are to their own cluster against the next nearest; higher is better, and it helps choose k.

The scaling line in that snippet is doing more work than it looks. Cluster 300 people on income in euros and age in years without scaling and the distance between two people is, for practical purposes, the difference in their incomes: a gap of 20,000 euros swamps a gap of 30 years. The result is three income bands whose average ages are 39.8, 39.9 and 39.9 — the age column had no effect whatsoever on a result presented as a segmentation by income and age.

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

rng = np.random.default_rng(0)
income = rng.normal(50_000, 15_000, 300)
age = rng.normal(40, 12, 300)
people = np.c_[income, age]

for name, data in (("unscaled", people), ("scaled", StandardScaler().fit_transform(people))):
    labels = KMeans(n_clusters=3, n_init=10, random_state=0).fit_predict(data)
    summary = [(round(income[labels == c].mean()), round(float(age[labels == c].mean()), 1))
               for c in range(3)]
    print(name, summary)

# unscaled [(67168, 39.8), (47117, 39.9), (29254, 39.9)]
# scaled   [(66167, 41.3), (40840, 49.1), (42458, 28.2)]

The n_init argument in that snippet is not decoration either. k-means starts from random centres and only finds a local optimum, so a single unlucky start can leave two centres inside one group and none in another. Run it ten times with different starts and keep the best, which is what n_init does. On six well-separated blobs, ten single-start runs gave inertias of 272, 300, 300, 300, 300, 272, 300, 2649, 3356 and 2350: three of the ten landed somewhere catastrophic, and four more settled ten per cent worse than the best.

from sklearn.datasets import make_blobs

blobs, _ = make_blobs(n_samples=600, centers=6, cluster_std=0.5, random_state=3)

singles = [round(KMeans(n_clusters=6, n_init=1, init="random", random_state=s)
                 .fit(blobs).inertia_, 1) for s in range(10)]
print(singles)          # [299.8, 300.2, 272.3, 2649.1, 299.8, 300.0, 271.9, 3355.7, 299.6, 2350.3]
print(round(KMeans(n_clusters=6, n_init=10, random_state=0).fit(blobs).inertia_, 1))
# 271.9   the best of the ten, which is what n_init gives you

When clusters are not round#

from sklearn.cluster import HDBSCAN

labels = HDBSCAN(min_cluster_size=20).fit_predict(X_scaled)
# -1 means 'noise': points that belong to no dense region

Density-based methods find clusters of any shape and leave outliers unassigned instead of forcing them into a group. You choose a minimum cluster size rather than the number of clusters, which is often easier to reason about.

The two interleaved crescents of the two-moons dataset are the standard case where this matters, and running both methods on it also shows why you cannot pick between them with a score. k-means gets it badly wrong, with an adjusted Rand index against the true grouping of 0.25, and HDBSCAN gets it exactly right, at 1.00. The silhouette scores go the other way: 0.489 for the wrong answer and 0.330 for the right one.

from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score, silhouette_score

moons, truth = make_moons(n_samples=500, noise=0.06, random_state=0)

km = KMeans(n_clusters=2, n_init=10, random_state=0).fit_predict(moons)
hd = HDBSCAN(min_cluster_size=20).fit_predict(moons)

print(round(adjusted_rand_score(truth, km), 3), round(silhouette_score(moons, km), 3))
# 0.253 0.489      wrong grouping, better silhouette
print(round(adjusted_rand_score(truth, hd), 3), round(silhouette_score(moons, hd), 3))
# 1.0 0.33         correct grouping, worse silhouette

The practical rule that falls out of this: use an internal score to compare runs of the same method at different settings, never to choose between methods that make different shape assumptions. Choosing between methods is done by looking at the result.

what your data looks like                   reach for            the setting you choose
round groups, similar sizes, you know k     k-means              k
round groups, you do not know k             k-means + gap stat   k, by comparison with noise
elliptical or overlapping groups            Gaussian mixture     number of components
odd shapes, varying density, real outliers  HDBSCAN              minimum cluster size
you want a hierarchy to cut at any level    agglomerative        linkage, then where to cut
millions of rows                            MiniBatchKMeans      k, and the batch size

What to skip here: the dozens of other clustering algorithms, the internal indices beyond silhouette and inertia, and any attempt to prove a clustering correct. Two methods with different assumptions, an honest comparison against noise, and the discipline of naming every cluster is the whole of what this subject gives you in practice.

Judge by reading#

Scores say clusters are separated, not that they mean anything. For each cluster, look at typical members and the features that distinguish it from the rest, and try to name it. Clusters you cannot name are often artefacts of scaling or of the method.

Make that a procedure rather than an intention, because it is the step everyone skips. For each cluster print its size, the five features whose mean differs most from the overall mean in standard deviations, and the three members closest to its centre. Then write a name for it in a sentence, in a document. A cluster you cannot describe in a sentence is not a finding, and the number of clusters you can name is a better answer to the question of how many clusters there are than any score on this page.

def describe(data, labels, feature_names, X_original):
    overall = data.mean(axis=0)
    spread = data.std(axis=0)
    for c in sorted(set(labels) - {-1}):
        members = labels == c
        gap = (data[members].mean(axis=0) - overall) / spread     # in standard deviations
        top = np.argsort(-np.abs(gap))[:5]
        print(f"cluster {c}: {members.sum()} members")
        for i in top:
            print(f"   {feature_names[i]:20} {gap[i]:+.2f} sd")
        centre = data[members].mean(axis=0)
        closest = np.argsort(((data[members] - centre) ** 2).sum(axis=1))[:3]
        print("   typical members:", X_original[members][closest])

These methods earn their place in language work on two jobs in particular. One is finding out what people are actually asking your product, by embedding a month of real questions and grouping them, which tends to contradict what the team assumed. The other is cleaning data: spotting near-duplicates at scale, and surfacing the unusual inputs that belong in an eval set precisely because nobody thought of them.

What are users asking?#

Embed a month of user questions, cluster the embeddings, and count cluster sizes. You get a map of real demand: which topics are common, which are growing, which your product never anticipated. It is one of the most useful analyses an AI engineer can run on a live product.

emb = embed(questions)                                   # [n, d], normalised
labels = HDBSCAN(min_cluster_size=25).fit_predict(emb)

for c in sorted(set(labels) - {-1}, key=lambda c: -(labels == c).sum())[:15]:
    sample = [q for q, l in zip(questions, labels) if l == c][:20]
    name = ask(f"Give a short name for what these questions have in common:\n{sample}")
    print((labels == c).sum(), name)

Having a model name clusters from samples is fast and usually good. Read a few clusters yourself to check the names are honest.

Two habits make the difference between an interesting picture and something a team acts on. Rank clusters by size and by growth rate between two time windows, not by size alone: a cluster that is fifth largest and has tripled this month is the one to take to a product meeting. And keep the noise label rather than forcing every question into a cluster, because the questions that fit nowhere are where the new use cases and the abuse both live.

Near-duplicates#

from datasketch import MinHash, MinHashLSH

def signature(text):
    m = MinHash(num_perm=128)
    for shingle in {text[i:i + 5] for i in range(len(text) - 4)}:
        m.update(shingle.encode())
    return m

lsh = MinHashLSH(threshold=0.8, num_perm=128)
for doc_id, text in docs.items():
    sig = signature(text)
    if lsh.query(sig):
        continue                     # near-duplicate of something already kept
    lsh.insert(doc_id, sig)

That code is worth understanding rather than pasting, because it rests on two ideas that are each one line of arithmetic. The first is the measure: Jaccard similarity is the size of the overlap between two sets divided by the size of their union. Cut each document into overlapping five-character pieces, called shingles, and two documents differing by one word share almost all their shingles.

The second is the trick that makes it cheap. Comparing every pair of a million documents is half a trillion comparisons. Instead, hash every shingle of a document with 128 different hash functions and keep the smallest result from each — that list of 128 numbers is the signature. The probability that two documents give the same minimum under one hash function is exactly their Jaccard similarity, so the fraction of the 128 positions where their signatures agree estimates it. Two short lists compared instead of two full documents.

import hashlib
import random

def shingles(text, k=5):
    return {text[i:i + k] for i in range(len(text) - k + 1)}

def jaccard(a, b):
    return len(a & b) / len(a | b)

def signature(shingle_set, permutations=128, seed=0):
    rng = random.Random(seed)
    salts = [rng.randrange(1 << 30) for _ in range(permutations)]
    return [min(int(hashlib.md5((str(salt) + s).encode()).hexdigest()[:8], 16)
                for s in shingle_set)
            for salt in salts]

def estimate(sig_a, sig_b):
    return sum(x == y for x, y in zip(sig_a, sig_b)) / len(sig_a)

d1 = "the quick brown fox jumps over the lazy dog"
d2 = "the quick brown fox jumped over the lazy dog"      # one word changed
d3 = "completely different text about machine learning models"

s1, s2, s3 = shingles(d1), shingles(d2), shingles(d3)
print(round(jaccard(s1, s2), 4), round(jaccard(s1, s3), 4))        # 0.7556 0.0
print(estimate(signature(s1), signature(s2)),
      estimate(signature(s1), signature(s3)))                       # 0.6953125 0.0

So choose the threshold knowing it is fuzzy. At 0.8 with 128 permutations you will keep some documents that are really 0.85 similar and drop some that are really 0.75. For deduplicating pretraining data that is fine, because the cost of either error is small and the alternative is not deduplicating at all. For deciding whether two customer records are the same person it is not fine, and you should be using the estimate to shortlist candidates and something exact to decide.

Unusual inputs make good test cases#

Points that fall in no cluster, or far from every cluster centre, are the odd requests: other languages, pasted logs, attempts to misuse the system. They are exactly what an eval set usually lacks.

Selecting them is four lines once the embeddings exist, and it is worth doing deliberately rather than taking whatever HDBSCAN happened to label as noise. Score every item by its distance to the nearest cluster centre, take the top few hundred, then sample across the range rather than only the extreme: the very furthest items are often broken input rather than interesting input.

import numpy as np

def outliers(emb, labels, per_bucket=10, buckets=5):
    """Rank by distance to the nearest cluster centre, then sample across the range."""
    centres = np.array([emb[labels == c].mean(axis=0) for c in sorted(set(labels) - {-1})])
    distance = np.linalg.norm(emb[:, None, :] - centres[None, :, :], axis=2).min(axis=1)
    order = np.argsort(-distance)[:per_bucket * buckets * 4]     # the odd tail
    picked = []
    for chunk in np.array_split(order, buckets):
        picked.extend(chunk[:per_bucket])
    return picked, distance

What this is worth, honestly#

None of these three techniques produces a number you can report. Clustering gives you a map to read, deduplication gives you a smaller corpus, outlier selection gives you a list of cases to look at. That is the right expectation to hold: unsupervised methods in language work are instruments for looking at data, and their output is always a thing a person reads next. The failure mode is treating a cluster id as a label and building a pipeline on top of it, at which point an arbitrary partition of your data becomes a fact about your product.