What it actually does
Start with every point as its own cluster, find the two closest clusters, merge them, and repeat until one cluster remains. The record of those merges — which two joined, and at what distance — is the dendrogram, and cutting it at any height or at any K gives a flat clustering for free. Everything interesting hides in one word, closest: the distance between two clusters is not defined until you choose a linkage. Single linkage uses the nearest pair of points, so clusters join along thin bridges and chain into long streaks. Complete linkage uses the farthest pair, so clusters stay compact. Average linkage splits the difference. Ward merges whichever pair increases the within-cluster sum of squares least, which is k-means’ objective applied greedily, and it is the only one of the four that needs Euclidean geometry. The others need only a dissimilarity — no mean, no coordinates, no metric axioms — so this is the clustering you reach for when your data is word sets compared by Jaccard, embeddings compared by cosine, or time series compared by dynamic time warping.
- A dissimilarity exists that reflects the similarity that matters; it is the model, and no linkage repairs a poor one — for Ward that dissimilarity must be Euclidean
- Clusters are nested or at least that a hierarchy is a useful description, which is a weaker assumption than being blobs but not an empty one
- The linkage’s notion of cluster shape matches the data: compact for complete and Ward, connected for single, somewhere between for average
- The dataset fits an n-by-n matrix, or can be reduced to something that does, which is the binding constraint in practice
- Early merges are reliable, because they are irrevocable — noisy or heavily overlapping data violates this and the tree inherits the error at every level above it
Loss function
At every step, join the pair of clusters whose union raises the total squared scatter by the least. That is a loss, and for Ward it is exactly k-means’ loss — but it is minimized greedily, one irrevocable merge at a time, rather than globally. A K-cut of the Ward tree is therefore generally not the K-means optimum on the same data; it is a good starting point for k-means rather than a substitute for it. Be honest about the other linkages: single, complete and average have no global objective being minimized. Each is a rule for what cluster-to-cluster distance means, applied greedily. Two narrow exceptions are worth stating. Single linkage cut at K is exactly the maximum-spacing clustering — it is Kruskal’s minimum spanning tree with the K-1 longest edges removed — so there is an optimality statement there. Complete linkage is a greedy heuristic for small cluster diameter with no such guarantee. For average linkage there is no objective at all.
- two current clusters; the algorithm considers every pair of them at every step
- Ward’s merge cost: how much the total within-cluster sum of squares rises if A and B are joined
- sum of squared Euclidean distances from a cluster’s points to its own mean — the same quantity k-means minimizes
- the cluster mean, which is why Ward exists only in a space where means and squared distances make sense
- the size factor that makes merging two large clusters expensive even when their centres are close — it is what pushes Ward toward balanced, compact clusters
How it is minimized
Greedy bottom-up merging, with cluster-to-cluster distances updated by the Lance-Williams recurrence
- the distance from any other cluster k to the cluster just formed by merging i and j — computed from old distances only, so the original points are never revisited
- weights on the two old distances: one half each for single and complete, size-proportional for average, and (n_i + n_k) over (n_i + n_j + n_k) for Ward
- weight on the distance between the two merged clusters: zero for single, complete and average, and minus n_k over (n_i + n_j + n_k) for Ward
- weight on the absolute difference: minus one half collapses the formula to min (single), plus one half to max (complete), and zero otherwise
- the term that lets one formula express both min and max, since min is the average minus half the gap and max is the average plus half the gap
The recurrence is the whole reason this is practical. Without it, every merge would recompute cluster-to-cluster distances from the original points, which for Ward means revisiting every pair; with it, each merge is one O(n) row update of the dissimilarity matrix and the data is never touched again. The linkage is nothing more than the choice of coefficients. A second decision decides the running time: the naive loop scans the whole matrix for the closest pair on every merge, which is O(n^3). A heap of nearest neighbours gets O(n^2 log n), and the nearest-neighbour chain algorithm gets O(n^2) outright — follow a chain of nearest neighbours until two clusters point at each other, merge them, and continue. That is only valid for reducible linkages, meaning merging two clusters can never bring the result closer to a third than either parent was. Single, complete, average and Ward are reducible. Centroid and median linkage are not, which is why they can produce inversions — a later merge at a smaller height than an earlier one — and why they should be avoided.
linkage
The definition of distance between clusters, and therefore the geometry of the result. Single chains along thin bridges. Complete and Ward favour compact clusters. Average sits between them. The choice is a statement about what a cluster looks like (Ward for Euclidean features; average for a general dissimilarity; single only for deliberately connectivity-driven or outlier-hunting work)
dissimilarity (metric)
What “close” means between two points. Euclidean is required by Ward. Cosine suits text and embedding vectors where direction carries meaning and magnitude is noise; Jaccard suits sets and binary indicators, such as shared devices or tokens; correlation distance suits series compared by shape; DTW or edit distance suit sequences and arrive as a precomputed matrix. Swapping the metric changes what the tree represents, not just a number (chosen from the data representation, never from the algorithm)
cut (n_clusters or distance_threshold)
Where to read the tree. A fixed K gives a flat partition; a height gives a partition whose number of clusters depends on the data, which is usually what you want for deduplication and entity resolution. This is the decision that k-means forces before fitting and this method defers until after
connectivity
An optional sparse graph restricting which pairs may merge — a spatial adjacency, a kNN graph. It injects structure the dissimilarity does not carry, and it lowers cost well below quadratic on sparse graphs (kNN graph with 10 to 20 neighbours, or a domain adjacency)
scaling
Standardize features before any Euclidean-based linkage, or the largest-variance feature decides the tree
There is nothing to converge: the procedure is deterministic, has no initialization and no restarts, and always terminates after exactly n-1 merges, which is a real operational advantage over k-means. Ties in the distance matrix are the one source of arbitrariness, and they make the tree depend on input order. What the method does not guarantee is that the tree is a good description of the data, and its failure modes are specific. Chaining: single linkage merges clusters through a single bridging pair, so a thin trail of points joins two blobs into one. Size bias: complete linkage breaks large clusters apart to keep diameters small, and Ward prefers balanced, spherical clusters. Irrevocability: every merge is final, so an early mistake caused by noise is carried into every level above it. Outlier sensitivity: under Ward an outlier inflates the sum of squares and distorts the partition at levels where it finally joins. And it always returns a tree whatever the data is — uniform noise produces a perfectly well-formed dendrogram. Cophenetic correlation, the correlation between the original dissimilarities and the dendrogram heights at which each pair first joins, is the standard diagnostic for whether the tree distorts the geometry.
Memory is O(n^2) for the dissimilarity matrix, which is the binding constraint: n(n-1)/2 float64 values is roughly 400 MB at 10,000 points and 10 GB at 50,000. Time depends on the algorithm. Naive scan-and-merge is O(n^3). A priority-queue implementation is O(n^2 log n). The nearest-neighbour chain algorithm is O(n^2) for any reducible linkage — single, complete, average, Ward. SLINK computes single linkage in O(n^2) time and O(n) extra memory by building the pointer representation one point at a time, never materializing the matrix, provided distances can be computed on the fly. With a sparse connectivity graph the cost drops well below quadratic. For anything past a few tens of thousands of points the honest answer is a two-stage scheme: reduce with mini-batch k-means, then cluster the centroids hierarchically.
How, where, and why it is used
Time-Series Forecasting
AdaptedCluster the series by shape using a distance chosen for that purpose, then use the tree two ways. Cut it to form groups over which one pooled model is fitted, which helps when no single series has enough history. Or use it as the aggregation structure for hierarchical forecasting when no natural hierarchy exists — forecast every node, then reconcile so that forecasts at different levels add up. The distance is the whole design: correlation distance, or Euclidean distance on z-normalized series (the two are equivalent, since squared Euclidean distance between z-normalized series of length T is 2T times one minus the correlation), or dynamic time warping when series are the same shape at different phases.
- Deriving a data-driven aggregation tree for demand series that have no clean product or region hierarchy, as input to forecast reconciliation
- Grouping SKUs, meters or stores by demand profile so one global model is fitted per group
- Load-profile segmentation in energy, where the dendrogram lets an analyst pick a coarser or finer segmentation after the fact
- Sensor and asset grouping where series are similar in shape but offset in time
It accepts a dissimilarity you can compute but not average, and that is the common situation with time series: DTW has no mean, so k-means cannot use it, and a DTW distance matrix goes straight into average or complete linkage. The dendrogram also answers the question k-means cannot, which is what granularity the groups should have — you cut once for pooling and again for reconciliation levels, from one fit. The costs bind quickly. DTW is itself quadratic in series length, so an n-by-n DTW matrix is quadratic in both, and DTW is not a metric, which rules out Ward entirely. There is no way to place a new series without redoing the work, so this is an offline step.
- Z-normalize each series before computing correlation or Euclidean distance, or the tree recovers magnitude rather than shape
- Use correlation distance or z-normalized Euclidean for Ward; use DTW only with single, complete or average linkage, where a non-metric dissimilarity is legitimate
- Constrain DTW with a warping band — it both bounds the cost and stops it matching unrelated shapes
- Compute distances on the training window only; a matrix built over the full history leaks future co-movement into the grouping
Rolling-origin backtest of the pooled or reconciled forecasts against the ungrouped baseline, reported at every level of the hierarchy. Cophenetic correlation is a diagnostic for whether the tree distorts the distances, not a measure of whether the grouping helps forecasts.
- Euclidean distance on unaligned, unnormalized series, where level and phase dominate any difference in shape
- Pairing Ward with a DTW matrix, which is not a Euclidean space, so the sum-of-squares criterion has no meaning
- Reading the tree’s top split as meaningful when it is only the largest merge height of a tree that always has one
- Forgetting that a refit can reshuffle the groups, silently changing which pooled model scores a series
Anomaly Detection
ViableBuild the tree on the data and score each point by how late it stops being isolated: the height at which it first joins a cluster of at least m points. Normal points join dense clusters early; an outlier remains a singleton until the tree is nearly complete. Equivalently, cut at a height and flag singletons and tiny clusters. Under single linkage a point’s first merge height is exactly its nearest-neighbour distance, so the score is a nearest-neighbour outlier score with the dendrogram supplying the context.
- Retrospective quality control on a bounded batch, such as sample or sensor screening before analysis
- Reviewing accounts or entities whose behaviour profile never merges with any established group
- Finding stray documents in a corpus clustered by cosine similarity — the ones that attach to nothing until the very top
- Anomaly review where an analyst needs to see why a point is isolated, which the tree shows directly
Its virtue is the dissimilarity: it works for data where no density or mean exists, such as Jaccard over set-valued records or Gower distance over mixed types, and it produces an explanation rather than a bare score. Against that, it is offline, quadratic, and has no way to score a new point without rebuilding. For streaming or large-scale detection an isolation forest or a density model is the right tool; this is a defensible choice for a bounded dataset reviewed once. Be specific about linkage: single linkage isolates outliers cleanly but chains them to clusters through any bridging point, while Ward and complete resist bridging but let an outlier distort the partition when it finally merges.
- Pick the dissimilarity for the data type — cosine for embeddings, Jaccard for sets — since the isolation it finds is isolation in that metric and no other
- Standardize numeric features before any Euclidean distance, or one wide feature decides what counts as isolated
- Choose the minimum cluster size m from the smallest group you consider a legitimate segment, not from the data
- Cluster a confirmed-clean window where one exists, or colluding anomalies form their own tidy cluster
Precision@k against confirmed anomalies, and PR-AUC rather than ROC-AUC. Compare against an isolation forest or LOF on the same data — if they agree, the tree’s extra explanation is free; if they disagree, the linkage choice is usually what differs.
- A group of anomalies larger than m, which forms its own cluster and is scored as normal — the same masking failure every cluster-based detector has
- A single bridging point under single linkage, which attaches an outlier to a dense cluster at a low height and hides it
- Comparing merge heights across different linkages or metrics, which are on different scales and mean different things
- Treating the merge height as a calibrated probability when it is a distance in an arbitrary metric
Optimization
Not applicableIt is a greedy procedure that produces a tree, not a method for solving a stated optimization problem: there is no feasible set to search, no fitness to rank candidates, and the only objective in the picture, Ward’s sum of squares, is minimized greedily and incidentally rather than offered as a solver for anyone’s problem.
Natural Language
ViableCluster documents, queries or support tickets by cosine distance on TF-IDF or sentence-embedding vectors, and use the dendrogram as a topic hierarchy a person can browse and cut at the granularity they need. The metric matters: cosine ignores vector length, which is document length or embedding norm and carries no topical information. If the vectors are L2-normalized first, squared Euclidean distance equals two times one minus the cosine, so Ward becomes legitimate on them as well.
- Building a topic taxonomy over support tickets or search queries, where the hierarchy is the deliverable
- Near-duplicate document grouping, cut at a distance threshold so the number of groups is whatever the data holds
- Exploratory analysis of an embedding space, where the dendrogram shows which topics are siblings
The number of topics is unknown, a fixed K would be arbitrary, and a browsable hierarchy is exactly the artefact an analyst wants. It scales poorly: a few tens of thousands of documents is the quadratic ceiling, so the practical pattern on a real corpus is two-stage — mini-batch k-means to a few thousand centroids, then hierarchical clustering on those. It is also the wrong choice when topics overlap heavily, because a hard tree cannot express a document that belongs to two branches, and a topic model or soft clustering fits that case better.
- L2-normalize embedding or TF-IDF rows before clustering, which makes cosine and Euclidean agree and permits Ward
- Remove boilerplate and near-duplicates first, or the tree’s lowest levels are dominated by templates
- Reduce very high-dimensional sparse vectors, or distances concentrate and every merge height looks alike
Human review of sampled clusters at the chosen cut, plus downstream task quality such as routing accuracy. Cophenetic correlation shows whether the tree preserves cosine distances; it does not show whether the topics are sensible.
- Single linkage chaining unrelated topics together through a generic boilerplate document
- Running the full quadratic method on a corpus that needed a first-stage reduction
- Choosing the cut height by eye on one sample and not revalidating it when the corpus grows
Risk & Fraud
ViableEntity resolution and account linking. Represent each account or record as a set of attributes — devices, addresses, phone numbers, name tokens — measure Jaccard distance between sets, and cut the tree at a distance threshold. Accounts that land in the same cluster plausibly belong to one actor. The linkage choice is a modelling decision here, and it is consequential: single linkage is transitive closure over the threshold graph, so any shared attribute links accounts; average or complete linkage requires the group to be mutually similar.
- Linking duplicate or synthetic identities across applications before underwriting
- Grouping accounts that share devices and contact details into candidate fraud rings for investigation
- Merging customer records across systems where no reliable key exists
Jaccard over set-valued attributes has no mean, so k-means is not even defined, and a distance threshold is the natural stopping rule because the number of real-world entities is unknown. The characteristic failure is a shared hub attribute — a public IP address, a common last name — which under single linkage bridges unrelated entities into one enormous cluster. That is why complete or average linkage is usually the better default for linking. The cost is the usual one: the quadratic matrix means blocking, comparing only within candidate buckets, is mandatory beyond modest volumes.
- Remove or down-weight high-frequency attributes such as shared IPs and common names, which carry no identifying information
- Block the data on a cheap key before computing pairwise distances, so the quadratic matrix is built per block rather than globally
- Use Jaccard or a weighted set similarity for set-valued attributes, never Euclidean distance on a one-hot encoding of them
Pairwise precision and recall of linked entities against confirmed matches, reported at several cut heights. Cluster-size distribution at the chosen threshold is a useful sanity check: one giant cluster means chaining.
- Single linkage turning one shared hub attribute into a single enormous false cluster
- A global distance threshold applied across blocks with very different attribute density
- Treating a linked cluster as a confirmed ring when it is a candidate for review
What it costs to run
Dominated by the n-by-n dissimilarity matrix, not the merging. Illustrative arithmetic: a full float64 matrix is about 800 MB at 10,000 points and 20 GB at 50,000, so memory is the ceiling long before time is. With nearest-neighbour chain the merge phase is O(n^2), and building the matrix is O(n^2 d) — or far more for an expensive dissimilarity such as DTW, where computing the matrix, not clustering it, is the real bill.
There is no inference. The tree describes the points it was built from and has no predict method. Placing a new point means recomputing, assigning it to the nearest cluster by the chosen linkage as an approximation, or — the usual pattern — fitting a classifier or storing centroids on the labels from a cut and serving that instead.
Batch, offline, and as rarely as the application allows. The procedure is deterministic, so a refit on the same data reproduces the same tree, which is better than k-means; but a refit on changed data can change the tree substantially, especially under single linkage, where one new bridging point can merge two clusters and rewrite everything above them. Anything downstream keyed to a cluster must be remapped after a refit.
- Track cophenetic correlation across refits; a fall means the tree distorts the dissimilarities more than it used to
- Watch the sorted top merge heights, since a collapse of the gap between them means the structure being cut no longer exists
- Monitor cluster-size distribution at the fixed cut height — one cluster absorbing most points is chaining, and it appears before downstream metrics move
- Compare partitions across refits with a set-overlap measure rather than by cluster id, which is arbitrary between runs
- Ward assumes Euclidean geometry and libraries cannot check it: handed a precomputed matrix of cosine or Jaccard distances, the tree comes back without complaint and the sum-of-squares criterion it claims to minimize does not exist
- The merge height has a convention that differs by library: the Lance-Williams recurrence on squared distances maintains twice the SSE increase, and scipy reports its square root, so heights from different tools are not comparable
- A tree built on a matrix that no longer fits in memory is not a configuration problem — the O(n^2) matrix is the algorithm, and the fix is a first-stage reduction or a sparse connectivity graph
- Chaining produces a dendrogram of depth close to n, and recursive traversal or plotting code that assumes balanced trees hits the recursion limit
- There is no predict: teams build a pipeline around the fit and discover afterward that a new record cannot be assigned without redoing the work
Pros and cons, in context
Needs only a dissimilarity, not a mean or coordinates
Single, complete and average linkage accept cosine, Jaccard, DTW, edit distance, Gower, or any precomputed matrix. That opens problems where centroid methods are not merely worse but undefined, and it is the main reason to choose this over k-means.
The number of clusters is chosen after fitting
One fit gives every flat clustering from one cluster to n, and the cut can be a K or a height. That removes the guess k-means requires and lets different consumers read the same tree at different granularities.
Deterministic, with no initialization or restarts
The same dissimilarity matrix always yields the same tree, up to tie-breaking. That makes results reproducible and debuggable in a way a stochastic method is not, and removes a whole class of tuning.
The dendrogram is an explanation, not just a label
It shows which groups are siblings and how far apart they are, so an analyst can see why two items were grouped. Useful for taxonomies, exploration and review; irrelevant when only a flat label feeds a downstream system.
Quadratic memory is a hard ceiling
The n-by-n matrix limits the method to roughly tens of thousands of points, and a more expensive dissimilarity such as DTW lowers that further. This decides most hierarchical-versus-k-means arguments before cluster quality is discussed.
Merges are greedy and irrevocable
An early mistake caused by noise is carried into every level above it, and nothing in the procedure revisits it. Ward is a greedy approximation to the k-means objective, so a K-cut of the tree is generally worse in SSE than k-means on the same data.
The linkage choice changes the answer fundamentally
Single chains, complete splits large clusters, Ward prefers balanced spherical ones, and on the same data they can disagree completely. There is no neutral default, so the choice has to be argued from the shape of the clusters expected.
No way to place a new point
The tree exists only for the points it was fitted on. Teams routinely discover this after building a pipeline around it, and the standard workaround is to train a classifier or store centroids from a cut, which is a different model.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""Agglomerative clustering - merge the closest pair, update, repeat. Transcribed.
Every point starts as its own cluster. Find the two closest clusters, merge
them, work out the distance from the new cluster to every other cluster with
the Lance-Williams formula, and record the merge. After n - 1 merges one
cluster remains, and the list of merges IS the dendrogram.
The algorithm takes a DISSIMILARITY MATRIX, not points. That is the main
advantage over k-means: single, complete and average linkage accept any
dissimilarity - cosine, Jaccard, a DTW matrix - because no cluster mean is ever
formed. Ward is the exception: it needs squared Euclidean distances.
"""
import math
def euclidean_matrix(X, squared):
"""Pairwise Euclidean distances. Ward needs them SQUARED; the rest do not."""
n = len(X)
D = [[0.0] * n for _ in range(n)]
for i in range(n):
for j in range(i + 1, n):
total = 0.0
for k in range(len(X[i])):
difference = X[i][k] - X[j][k]
total += difference * difference
D[i][j] = total if squared else math.sqrt(total)
D[j][i] = D[i][j]
return D
def lance_williams(linkage, d_ki, d_kj, d_ij, n_i, n_j, n_k):
"""Distance from cluster k to the merge of i and j, from old distances only.
d(k, i+j) = a_i d(k,i) + a_j d(k,j) + b d(i,j) + g |d(k,i) - d(k,j)|
One formula, four sets of coefficients. The coefficients ARE the linkage.
"""
if linkage == "single": # a = 1/2, b = 0, g = -1/2: this is min
return 0.5 * d_ki + 0.5 * d_kj - 0.5 * abs(d_ki - d_kj)
if linkage == "complete": # g = +1/2: this is max
return 0.5 * d_ki + 0.5 * d_kj + 0.5 * abs(d_ki - d_kj)
if linkage == "average": # size-weighted mean of the two distances
return (n_i * d_ki + n_j * d_kj) / (n_i + n_j)
if linkage == "ward": # valid ONLY on squared Euclidean distances
total = n_i + n_j + n_k
return ((n_i + n_k) * d_ki + (n_j + n_k) * d_kj - n_k * d_ij) / total
raise ValueError(f"unknown linkage: {linkage}")
def fit(D, linkage="ward"):
"""Returns merges as (left_id, right_id, height, size), in merge order.
Ids below n are original points; the cluster formed at step s gets id n + s.
"""
n = len(D)
D = [list(row) for row in D] # work on a copy
size = [1] * n
label = list(range(n))
alive = [True] * n
merges = []
for step in range(n - 1):
# ---- find the closest pair of live clusters: the O(n^2) scan ------
best = None
for i in range(n):
if not alive[i]:
continue
for j in range(i + 1, n):
if not alive[j]:
continue
if best is None or D[i][j] < best[0]:
best = (D[i][j], i, j)
height, i, j = best
# ---- update every other cluster's distance to the merged one -----
for k in range(n):
if not alive[k] or k == i or k == j:
continue
updated = lance_williams(
linkage, D[k][i], D[k][j], D[i][j], size[i], size[j], size[k]
)
D[k][i] = updated
D[i][k] = updated
merges.append((label[i], label[j], height, size[i] + size[j]))
size[i] += size[j]
alive[j] = False
label[i] = n + step
return merges
def cut(merges, n, n_clusters):
"""Flat labels from the first n - n_clusters merges."""
members = {point: [point] for point in range(n)}
for step in range(n - n_clusters):
left, right, _, _ = merges[step]
members[n + step] = members.pop(left) + members.pop(right)
labels = [0] * n
for cluster, points in enumerate(members.values()):
for point in points:
labels[point] = cluster
return labelsO(n^3) time — every one of the n-1 merges rescans an O(n^2) matrix for the minimum — and O(n^2) memory, over lists of lists.
How this language actually wants it written.
"""Hierarchical clustering - typed, nearest-neighbour chain, Ward kept honest."""
from dataclasses import dataclass
from enum import Enum
import numpy as np
from numpy.typing import NDArray
Matrix = NDArray[np.float64]
Labels = NDArray[np.int64]
class Linkage(str, Enum):
SINGLE = "single"
COMPLETE = "complete"
AVERAGE = "average"
WARD = "ward"
@dataclass(frozen=True)
class Dendrogram:
"""n - 1 merges as rows of [left_id, right_id, height, size].
Ids below n_points are original points; id n_points + s is the cluster
formed at step s. Rows are sorted by height, so a cut at a height is a
prefix of the merges.
"""
merges: Matrix
n_points: int
def cut_into(self, n_clusters: int) -> Labels:
if not 1 <= n_clusters <= self.n_points:
raise ValueError(f"n_clusters must lie in [1, {self.n_points}], got {n_clusters}")
return self._replay(self.n_points - n_clusters)
def cut_at(self, height: float) -> Labels:
if not np.isfinite(height) or height < 0.0:
raise ValueError(f"height must be finite and non-negative, got {height}")
steps = int(np.searchsorted(self.merges[:, 2], height, side="right"))
return self._replay(steps)
def cophenetic_correlation(self, dissimilarity: Matrix) -> float:
"""Correlation between the original dissimilarities and the heights at
which each pair first joins - how faithfully the tree keeps the
geometry. Compare against the SAME dissimilarity the tree was built from."""
n = self.n_points
cophenetic = np.zeros((n, n))
members = {point: np.array([point]) for point in range(n)}
for step, (left, right, height, _) in enumerate(self.merges):
first = members.pop(int(left))
second = members.pop(int(right))
cophenetic[np.ix_(first, second)] = height
cophenetic[np.ix_(second, first)] = height
members[n + step] = np.concatenate((first, second))
upper = np.triu_indices(n, k=1)
return float(np.corrcoef(dissimilarity[upper], cophenetic[upper])[0, 1])
def _replay(self, steps: int) -> Labels:
parent = np.arange(self.n_points + len(self.merges))
for step in range(steps):
parent[int(self.merges[step, 0])] = self.n_points + step
parent[int(self.merges[step, 1])] = self.n_points + step
roots = np.empty(self.n_points, dtype=np.int64)
for point in range(self.n_points):
node = point
while parent[node] != node:
node = parent[node]
roots[point] = node
_, labels = np.unique(roots, return_inverse=True)
return labels.astype(np.int64)
def _lance_williams(
linkage: Linkage,
d_ki: Matrix,
d_kj: Matrix,
d_ij: float,
n_i: float,
n_j: float,
n_k: Matrix,
) -> Matrix:
"""Vectorized over every other cluster k at once.
With the coefficients substituted, single and complete collapse to min and
max, so the absolute-value term never has to be evaluated.
"""
if linkage is Linkage.SINGLE:
return np.minimum(d_ki, d_kj)
if linkage is Linkage.COMPLETE:
return np.maximum(d_ki, d_kj)
if linkage is Linkage.AVERAGE:
return (n_i * d_ki + n_j * d_kj) / (n_i + n_j)
total = n_i + n_j + n_k
return ((n_i + n_k) * d_ki + (n_j + n_k) * d_kj - n_k * d_ij) / total
def _nearest_neighbour_chain(
D: Matrix, linkage: Linkage
) -> list[tuple[int, int, float]]:
"""O(n^2) agglomeration, valid because every supported linkage is reducible:
merging two clusters never brings the result closer to a third than either
parent was, so mutual nearest neighbours can be merged in ANY order and the
merges sorted by height afterwards."""
n = D.shape[0]
D = D.copy() # the caller's matrix is not ours to destroy
np.fill_diagonal(D, np.inf)
sizes = np.ones(n)
active = np.ones(n, dtype=bool)
chain: list[int] = []
raw: list[tuple[int, int, float]] = []
while len(raw) < n - 1:
if not chain:
chain.append(int(np.flatnonzero(active)[0]))
while True:
tail = chain[-1]
nearest = int(np.argmin(D[tail]))
if len(chain) > 1:
previous = chain[-2]
if D[tail, previous] <= D[tail, nearest]: # ties go backward: terminates
nearest = previous
if nearest == previous:
break
chain.append(nearest)
b = chain.pop()
a = chain.pop()
height = float(D[a, b])
others = np.flatnonzero(active)
others = others[(others != a) & (others != b)]
updated = _lance_williams(
linkage, D[others, a], D[others, b], height, sizes[a], sizes[b], sizes[others]
)
D[others, a] = updated
D[a, others] = updated
D[b, :] = np.inf
D[:, b] = np.inf
sizes[a] += sizes[b]
active[b] = False
raw.append((a, b, height))
return raw
def _relabel(raw: list[tuple[int, int, float]], n: int, take_root: bool) -> Matrix:
"""Sort merges by height and rewrite slot indices as dendrogram ids."""
order = sorted(range(len(raw)), key=lambda step: raw[step][2]) # stable
parent = list(range(n))
cluster_id = list(range(n))
size = [1] * n
merges = np.empty((n - 1, 4))
def find(slot: int) -> int:
while parent[slot] != slot:
parent[slot] = parent[parent[slot]]
slot = parent[slot]
return slot
for step, index in enumerate(order):
slot_a, slot_b, height = raw[index]
root_a, root_b = find(slot_a), find(slot_b)
left, right = sorted((cluster_id[root_a], cluster_id[root_b]))
merges[step] = (left, right, np.sqrt(height) if take_root else height,
size[root_a] + size[root_b])
parent[root_b] = root_a
size[root_a] += size[root_b]
cluster_id[root_a] = n + step
return merges
def linkage_from_points(X: Matrix, linkage: Linkage = Linkage.WARD) -> Dendrogram:
"""Euclidean geometry is available, so every linkage is legitimate."""
if X.ndim != 2 or X.shape[0] < 2:
raise ValueError(f"X must be 2-D with at least 2 rows, got shape {X.shape}")
difference = X[:, None, :] - X[None, :, :]
squared = (difference * difference).sum(axis=2)
if linkage is Linkage.WARD:
# Lance-Williams for Ward is exact on SQUARED distances, where it
# maintains twice the increase in SSE; the root is the usual height.
raw = _nearest_neighbour_chain(squared, linkage)
return Dendrogram(_relabel(raw, X.shape[0], take_root=True), X.shape[0])
raw = _nearest_neighbour_chain(np.sqrt(squared), linkage)
return Dendrogram(_relabel(raw, X.shape[0], take_root=False), X.shape[0])
def linkage_from_dissimilarity(D: Matrix, linkage: Linkage = Linkage.AVERAGE) -> Dendrogram:
"""Any dissimilarity: cosine, Jaccard, DTW, edit distance. NOT Ward."""
if linkage is Linkage.WARD:
raise ValueError(
"ward minimises a sum of squares around cluster means, which exists "
"only in Euclidean space; use linkage_from_points, or choose single, "
"complete or average for a general dissimilarity"
)
if D.ndim != 2 or D.shape[0] != D.shape[1] or D.shape[0] < 2:
raise ValueError(f"D must be square with at least 2 rows, got shape {D.shape}")
if not np.isfinite(D).all() or (D < 0.0).any():
raise ValueError("D must be finite and non-negative")
if not np.allclose(D, D.T):
raise ValueError("D must be symmetric")
raw = _nearest_neighbour_chain(D, linkage)
return Dendrogram(_relabel(raw, D.shape[0], take_root=False), D.shape[0])What changed: Three changes, and the first is the one that changes the complexity class. The full-matrix scan on every merge, which made the previous stage O(n^3), becomes the nearest-neighbour chain: follow nearest neighbours until two clusters point at each other, merge them, and carry on. That is valid only because single, complete, average and Ward are reducible, so merges can be found out of order and sorted by height afterwards. Second, the Lance-Williams update is vectorized across all other clusters, with single and complete collapsing to min and max. Third, Ward is moved behind an explicit boundary: a precomputed dissimilarity cannot carry Euclidean geometry, so linkage_from_dissimilarity raises on Ward and linkage_from_points computes the squared distances itself. The result is a frozen dataclass that can cut by K or by height and report its own cophenetic correlation.
O(n^2) time in the chain, but the broadcast difference materializes an (n, n, d) intermediate — which is the memory ceiling, not the arithmetic.
Faster, and the justification for each change.
"""Hierarchical clustering - GEMM distances, in-place nearest-neighbour chain."""
import numpy as np
from numpy.typing import NDArray
Matrix = NDArray[np.float64]
def pairwise_distances(X: Matrix, squared: bool) -> Matrix:
"""Distances written straight into ONE preallocated n-by-n buffer.
||a - b||^2 = ||a||^2 - 2 a.b + ||b||^2: only the cross term touches both
operands, and it is a GEMM. The broadcast form it replaces materializes an
(n, n, d) array - at 10,000 points and 50 features that is 40 GB before the
algorithm has started.
"""
design = np.ascontiguousarray(X, dtype=np.float64)
n = design.shape[0]
norms = np.einsum("ij,ij->i", design, design)
D = np.empty((n, n), dtype=np.float64)
np.dot(design, design.T, out=D) # GEMM directly into D
D *= -2.0
D += norms[:, None]
D += norms[None, :]
np.maximum(D, 0.0, out=D) # cancellation can go slightly negative
if not squared:
np.sqrt(D, out=D)
np.fill_diagonal(D, np.inf)
return D
def nearest_neighbour_chain(D: Matrix, linkage: str) -> Matrix:
"""Consumes D in place. Returns raw merges as rows of [slot_a, slot_b, height].
Retired clusters are marked by +inf rows and columns rather than removed, so
D never changes shape, the row update runs over the full row with no index
arrays, and an argmin can never pick a dead cluster.
"""
n = D.shape[0]
sizes = np.ones(n)
active = np.ones(n, dtype=bool)
scratch = np.empty(n)
scratch_b = np.empty(n)
raw = np.empty((n - 1, 3))
chain = np.empty(n, dtype=np.int64)
depth = 0
first_active = 0
for step in range(n - 1):
if depth == 0:
while not active[first_active]:
first_active += 1
chain[0] = first_active
depth = 1
while True:
tail = chain[depth - 1]
nearest = int(D[tail].argmin())
if depth > 1:
previous = int(chain[depth - 2])
if D[tail, previous] <= D[tail, nearest]:
nearest = previous
if nearest == previous:
break
chain[depth] = nearest
depth += 1
b = int(chain[depth - 1])
a = int(chain[depth - 2])
depth -= 2
n_a, n_b, d_ab = sizes[a], sizes[b], float(D[a, b])
row_a, row_b = D[a], D[b]
if linkage == "single":
np.minimum(row_a, row_b, out=row_a)
elif linkage == "complete":
np.maximum(row_a, row_b, out=row_a)
elif linkage == "average":
row_a *= n_a
np.multiply(row_b, n_b, out=scratch)
row_a += scratch
row_a /= n_a + n_b
else:
# Ward: ((n_a + n_k) d_ak + (n_b + n_k) d_bk - n_k d_ab) / (n_a + n_b + n_k)
np.add(sizes, n_a, out=scratch)
scratch *= row_a
np.add(sizes, n_b, out=scratch_b)
scratch_b *= row_b
scratch += scratch_b
np.multiply(sizes, d_ab, out=scratch_b)
scratch -= scratch_b
np.add(sizes, n_a + n_b, out=scratch_b)
np.divide(scratch, scratch_b, out=row_a)
D[:, a] = row_a # keep the matrix symmetric
D[b, :] = np.inf
D[:, b] = np.inf
D[a, a] = np.inf
sizes[a] = n_a + n_b
active[b] = False
raw[step] = (a, b, d_ab)
# Sorting by height and relabelling slots as dendrogram ids is unchanged from
# the previous stage: O(n log n), nowhere near the hot path.
return rawWhat changed: The previous stage’s broadcast difference materializes an (n, n, d) array, which is the real memory ceiling, so distances become a GEMM written directly into the one n-by-n buffer the algorithm needs anyway. Inside the chain, the Lance-Williams update is performed in place on row views with two preallocated scratch vectors, and retired clusters are marked with infinity instead of being indexed out, which removes the per-merge index arrays and temporaries. The algorithm and its O(n^2) memory are unchanged — this stage buys a constant factor and a smaller peak, not a different complexity class. Past a few tens of thousands of points the right move is not faster dense code but a first-stage reduction or a sparse connectivity graph.
Vectorize to NumPy/BLAS (@, np.dot, einsum)
The cross term of every pairwise distance is one GEMM that writes directly into the output, and the row norms are einsum reductions computed once rather than per pair.
Trade-off: The squared-norm expansion loses precision through cancellation, so two nearly-equidistant pairs can be misranked — which changes the merge order, and under single linkage which edges end up in the spanning tree, not merely a distance value.
Pre-allocate output arrays and use in-place operations
The distance matrix, two scratch vectors and the chain are allocated once for the whole fit, and every Lance-Williams update runs in place on a row view instead of allocating temporaries per merge.
Trade-off: The scratch buffers and the matrix are shared across merges, so the function is not reentrant and it destroys its input; a reordered in-place line corrupts the tree silently rather than failing.
Fuse adjacent operations so intermediates are never materialized
Single and complete linkage reduce to one np.minimum or np.maximum with an output argument, and Ward’s update is assembled in two scratch buffers rather than five temporaries.
Trade-off: The Ward expression is far less readable than the formula it implements, and the order of the in-place steps matters — each one overwrites a buffer a later one needs.
O(n^2 d) distances in BLAS and an O(n^2) in-place chain, with one n-by-n buffer and no (n, n, d) intermediate. Illustrative, not a measured benchmark.