Hierarchical clustering

What this chapter covers

k-means made you pick k before it would tell you anything. Hand it three clusters when the data has five and it will still hand you three, splitting real groups to hit the number you demanded. This chapter takes that demand away. We build a clustering that never fixes k — it builds the whole tree of groupings, from every point alone at the bottom to one big group at the top, and lets you cut it wherever the structure tells you to.

The method is agglomerative hierarchical clustering, and the idea is almost insultingly simple: start with every point as its own cluster, find the two closest clusters, fuse them, and repeat until one cluster is left. What "closest" means between two groups of points is the one real choice you make — the linkage rule — and it changes everything, so we implement three of them and make it selectable. We build the whole thing in NumPy, record every merge and the distance it happened at, and use those to draw a dendrogram: the tree of merges that is the real output here, not a set of labels. Then we cut the tree at k and check that the flat clusters we get match scikit-learn's exactly.

The data is thirty 2-D points in three tidy blobs — small on purpose, because a dendrogram of thirty leaves is readable and a dendrogram of ten thousand is a black smear. That size limit is not incidental. It's the whole trade this chapter is about.

A bit of history

Hierarchical clustering grew up in biology, not computer science. The people who needed it were taxonomists trying to build the tree of life from measurements — if you have forty beetles and a table of their traits, which two are most alike, and which groups of beetles form a genus? In 1963 Robert Sokal and Peter Sneath published Principles of Numerical Taxonomy, the book that argued you could classify organisms by algorithm instead of expert judgment, computing similarities from raw measurements and letting the groupings fall out. The dendrograms in this chapter are their diagrams. The whole framing — merge the most similar things, draw the tree, read the branches — comes from that program of putting numbers under classification.

The same year, 1963, Joe Ward wrote down a different linkage rule from a different angle: instead of merging the closest clusters by distance, merge the pair whose fusion increases the within-cluster variance the least. Ward's method is the one you'll most often see as a default in practice, because it tends to produce compact, balanced clusters, and it ties the whole procedure back to the same sum-of-squares objective k-means minimizes. The single-link rule we lead with is older still in spirit, and reads like the taxonomist's instinct: two groups are as close as their two nearest members. Sixty years on, the algorithm hasn't changed. What changed is that we stopped drawing the trees by hand.

The intuition

Forget centroids. There are none here, and no k to guess. Put every point in a cluster by itself — thirty points, thirty clusters — and then do the one thing the algorithm knows how to do: find the two clusters that are closest together and glue them into one. Now you have twenty-nine clusters. Do it again: twenty-eight. Keep going and the little clusters fuse into bigger ones, the bigger ones fuse into blobs, and the blobs finally fuse into a single cluster holding everything. Twenty-nine merges for thirty points, and you're done.

Nothing about that asked how many groups there are. The algorithm doesn't decide — it records. Every merge happened at some distance, and if you write those distances down in order you get a history: these two points were basically on top of each other, this pair of clumps was a little farther apart, and right at the end two big well-separated blobs got dragged together across a huge gap because nothing else was left to merge. That history is the dendrogram, and the gaps in it are where the real structure lives. A merge across a tiny distance is joining things that belong together; a merge across a big jump is fusing two groups that didn't really want to be one. You pick your number of clusters by cutting the tree below the big jumps, after you've seen them, instead of guessing before you start.

Here's the raw data — thirty points, no labels, gray, three blobs your eye finds without help:

The math

The only quantity the algorithm needs is a distance between two clusters, and every cluster is just a set of the original points it contains. Write the point distance as $\lVert x_a - x_b \rVert$, the ordinary Euclidean distance between points $a$ and $b$. For two clusters $A$ and $B$, the linkage rule says how to turn all the point-to-point distances across them into one number.

Single linkage takes the closest pair — the smallest distance between any point of $A$ and any point of $B$:

Complete linkage takes the farthest pair instead, so two clusters count as close only if even their most distant members are near:

Average linkage splits the difference and takes the mean over all cross-cluster pairs, where $|A|$ and $|B|$ are the number of points in each cluster:

Those three formulas are the whole algorithm's personality. Single linkage follows chains of near neighbors and can string separate groups together through a bridge of points; complete linkage insists on tight balls and can chop a long cluster in half; average sits between them and is usually the safe first try, which is why it's the one we headline. Once you've picked a rule, the procedure is fixed: at every step merge the pair $(A, B)$ with the smallest $d(A, B)$, record the distance as the merge height, and repeat. No objective is being minimized here the way k-means minimizes inertia — the tree is just the transcript of a greedy sequence of nearest-merges.

What it's good at, what it isn't

The appeal is that it answers a question k-means can't even hear: how many clusters are there? You don't tell it — you build the full tree once and read the answer off the merge heights, cutting below whatever gap looks real. You get the whole nested structure for free, so you can see that these two subgroups sit inside that bigger group, which matters when the thing you're clustering genuinely has a hierarchy: taxonomies, document topics, gene families. It's deterministic — same data, same linkage, same tree, every run, with no seed to babysit and no random restart to hope on. And with single linkage it can trace out long, thin, non-round clusters that k-means would carve straight through, because it only ever cares about near neighbors, not distance to a center.

The catch is that it does not scale, full stop. You need the distances between all pairs of points, which is an $n \times n$ matrix — quadratic memory before you've merged anything — and the straightforward algorithm scans the live pairs at every one of the n−1 steps, which lands you around $O(n^3)$ time. On thousands of points that's already slow; on millions it's a non-starter, and k-means, linear in the number of points, wins by default. It's also greedy and committed: a merge, once made, is never revisited, so an early wrong fusion propagates all the way up the tree. And the linkage choice isn't a detail you can skip — single, complete, and average can hand you three genuinely different trees on the same data, so "hierarchical clustering says these are the groups" is never a complete sentence without naming the rule.

The data

Three Gaussian blobs from scikit-learn's make_blobs, thirty points in 2-D, ten per blob, with a fixed random seed so the snapshot never moves. Small on purpose: this is the one algorithm in the course where the visualization degrades with size, because a dendrogram draws one leaf per point and thirty leaves is exactly as many as you can read. Two features so the whole thing fits on a page; three blobs spaced far enough apart that the right cut is obvious to a human, which is what you want when you're checking whether the tree agrees. The generator also hands back a true blob id per point. We commit those and keep them sealed — the coordinates build the tree, the ids come out only at the end to score how well the flat cut recovered the structure that was really there.

Build it, one function at a time

Four functions. One computes distances, one turns a linkage rule into a cluster-to-cluster distance, one grows the whole tree, and one cuts it. They stack in that order.

Everything rests on the pairwise distances between points, computed once:

def pairwise(X):
    """Euclidean distance from every point to every other point.

    X is (N, D); the result is an (N, N) matrix where entry (i, j) is
    ||x_i - x_j||. This is computed once, up front — every linkage rule below
    is just a way of reducing a block of this matrix down to one number.
    """
    diff = X[:, None, :] - X[None, :, :]      # (N, N, D)
    return np.sqrt((diff ** 2).sum(axis=2))   # (N, N)

That (N, N) matrix is the raw material. Every linkage rule is a way of reducing a rectangular block of it — the distances between the points of one cluster and the points of another — down to a single number:

def cluster_distance(members_a, members_b, D, linkage):
    """Distance between two clusters, reduced from the point distance matrix D.

    A cluster is just a list of the original point indices it contains. Pull
    the block of D between the two member sets and collapse it to a single
    number the way the chosen linkage says to:

      single   = the closest pair of points across the two clusters (min)
      complete = the farthest pair (max)
      average  = the mean over all cross-cluster pairs

    That one choice is the whole personality of the algorithm.
    """
    block = D[np.ix_(members_a, members_b)]
    if linkage == "single":
        return float(block.min())
    if linkage == "complete":
        return float(block.max())
    if linkage == "average":
        return float(block.mean())
    raise ValueError(f"unknown linkage: {linkage}")

Single takes the block's minimum, complete its maximum, average its mean. Three lines, three completely different clusterings; the rest of the code doesn't care which one you pick. Now the loop that grows the tree. Start with every point in its own cluster, then n−1 times over, scan all pairs of live clusters, find the closest under the linkage rule, and fuse them:

def agglomerate(X, linkage="average"):
    """Grow the merge tree: fuse the two nearest clusters, N-1 times.

    Clusters are held in a dict id -> list of member point indices. Every point
    starts as its own cluster with id 0..N-1. At each step we scan all pairs of
    live clusters, find the closest under the linkage rule, merge them into a
    new cluster (given the next fresh id, scipy-style: N, N+1, ...), and record
    the merge.

    Returns a scipy-compatible linkage matrix Z with one row per merge —
    [id_a, id_b, height, size] — plus a `history` list, one entry per merge,
    that the chapter replays frame by frame to animate the agglomeration.
    """
    D = pairwise(X)
    n = len(X)
    clusters = {i: [i] for i in range(n)}     # id -> member point indices
    next_id = n
    Z = []
    history = []

    for _ in range(n - 1):
        # find the closest pair of live clusters under the linkage rule
        ids = list(clusters)
        best = None
        for a_pos in range(len(ids)):
            for b_pos in range(a_pos + 1, len(ids)):
                ia, ib = ids[a_pos], ids[b_pos]
                d = cluster_distance(clusters[ia], clusters[ib], D, linkage)
                if best is None or d < best[0]:
                    best = (d, ia, ib)

        height, ia, ib = best
        merged = clusters[ia] + clusters[ib]
        history.append({
            "id_a": ia, "id_b": ib, "new_id": next_id,
            "members_a": list(clusters[ia]), "members_b": list(clusters[ib]),
            "height": height, "size": len(merged),
            "labels": _labels_now(clusters, n),   # cluster id per point, pre-merge
        })
        Z.append([ia, ib, height, len(merged)])

        del clusters[ia]
        del clusters[ib]
        clusters[next_id] = merged
        next_id += 1

    return np.array(Z, dtype=float), history

Every merge is recorded twice: once in the linkage matrix Z in scipy's own format — [id_a, id_b, height, size], with each new cluster getting the next fresh id — and once in a fuller history entry that also snapshots which cluster every point belonged to just before the merge, which is what the animation replays. The merged cluster inherits both children's members and a new id, so the ids climb from n upward and the last merge's id is the root of the tree. This is the greedy, readable version: O(n²) memory for the distance matrix and a full pair scan at every step. It is not how you'd cluster a million points, and that's the point of the chapter.

The tree by itself isn't a clustering — it's every clustering at once, one for each place you could cut it. To get flat labels at a specific k, replay the tree from the bottom and stop early:

def cut_tree(Z, n, k):
    """Cut the merge tree to get k flat clusters.

    Every row of Z is a merge; the first N-k of them are the merges that happen
    below the cut line. Replay exactly those with a union-find, and whatever
    connected components remain are the k flat clusters. Labels are remapped to
    a dense 0..k-1 so they line up with what a library returns.
    """
    parent = list(range(2 * n - 1))

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x

    for row, (a, b, _h, _s) in enumerate(Z[:n - k]):
        new = n + row
        parent[int(a)] = new
        parent[int(b)] = new

    roots = [find(i) for i in range(n)]
    order = {}
    labels = np.empty(n, dtype=int)
    for i, r in enumerate(roots):
        if r not in order:
            order[r] = len(order)
        labels[i] = order[r]
    return labels

The first n−k merges are the ones below the cut line; play exactly those with a union-find and whatever connected components survive are your k clusters. Cut at k=1 and you get the whole dataset; cut at k=n and you get every point alone; cut anywhere between and you read a clustering off the same tree without ever rebuilding it.

Watch it work

This is why thirty points and not thirty thousand. Below is a real run of the agglomerate function above on the three blobs, one merge per frame, average linkage. The top panel is the points, colored by their current cluster — gray while a point is still a singleton, a solid color once it's joined something. The red dashed line snaps between the two clusters about to fuse this step. The bottom panel is the dendrogram building itself bar by bar: every merge adds a bracket at its own height, the current one in red, so you watch the tree grow from the leaves up while the points condense on the plane above it.

Press play and watch the order the merges happen in. The first fusions are tiny — points that are nearly on top of each other, joining at distance 0.19, 0.29, 0.33 — and they all happen inside the blobs, low on the dendrogram. The gray dust condenses into three colored clumps well before anything crosses between blobs. Then the character changes. Once each blob is a single cluster, the only merges left are between blobs, and they happen across a huge gap: the last within-blob merge is at distance 1.90, and then the tree jumps to 5.94 to join two blobs and 8.97 to swallow the third. Those two tall brackets at the top, stranded far above everything below them, are the algorithm telling you there are three real groups here. You didn't pick three. You watched the merge distances announce it.

Reset and run it again — it's deterministic, so it draws the exact same tree every time. There's no seed here, no random start, nothing to get lucky or unlucky on. Same data, same linkage, same transcript.

The full implementation

The whole file, no library, top to bottom. This is exactly what the animation ran:

"""Agglomerative hierarchical clustering, built from scratch.

No k up front, no centroids. Start with every point as its own cluster, then
repeatedly fuse the two closest clusters — where "closest" is decided by a
linkage rule — until a single cluster holds everything. Record every merge and
its height (the distance the two clusters merged at) so we can draw a dendrogram
and cut the tree at any k we like. Pure NumPy; no ML library in this file.

Every function below appears in the chapter one step at a time (the
`# region:` markers are what the book's include directives pull in).
"""

import numpy as np
import pandas as pd


# region: pairwise
def pairwise(X):
    """Euclidean distance from every point to every other point.

    X is (N, D); the result is an (N, N) matrix where entry (i, j) is
    ||x_i - x_j||. This is computed once, up front — every linkage rule below
    is just a way of reducing a block of this matrix down to one number.
    """
    diff = X[:, None, :] - X[None, :, :]      # (N, N, D)
    return np.sqrt((diff ** 2).sum(axis=2))   # (N, N)
# endregion


# region: linkage_rules
def cluster_distance(members_a, members_b, D, linkage):
    """Distance between two clusters, reduced from the point distance matrix D.

    A cluster is just a list of the original point indices it contains. Pull
    the block of D between the two member sets and collapse it to a single
    number the way the chosen linkage says to:

      single   = the closest pair of points across the two clusters (min)
      complete = the farthest pair (max)
      average  = the mean over all cross-cluster pairs

    That one choice is the whole personality of the algorithm.
    """
    block = D[np.ix_(members_a, members_b)]
    if linkage == "single":
        return float(block.min())
    if linkage == "complete":
        return float(block.max())
    if linkage == "average":
        return float(block.mean())
    raise ValueError(f"unknown linkage: {linkage}")
# endregion


# region: agglomerate
def agglomerate(X, linkage="average"):
    """Grow the merge tree: fuse the two nearest clusters, N-1 times.

    Clusters are held in a dict id -> list of member point indices. Every point
    starts as its own cluster with id 0..N-1. At each step we scan all pairs of
    live clusters, find the closest under the linkage rule, merge them into a
    new cluster (given the next fresh id, scipy-style: N, N+1, ...), and record
    the merge.

    Returns a scipy-compatible linkage matrix Z with one row per merge —
    [id_a, id_b, height, size] — plus a `history` list, one entry per merge,
    that the chapter replays frame by frame to animate the agglomeration.
    """
    D = pairwise(X)
    n = len(X)
    clusters = {i: [i] for i in range(n)}     # id -> member point indices
    next_id = n
    Z = []
    history = []

    for _ in range(n - 1):
        # find the closest pair of live clusters under the linkage rule
        ids = list(clusters)
        best = None
        for a_pos in range(len(ids)):
            for b_pos in range(a_pos + 1, len(ids)):
                ia, ib = ids[a_pos], ids[b_pos]
                d = cluster_distance(clusters[ia], clusters[ib], D, linkage)
                if best is None or d < best[0]:
                    best = (d, ia, ib)

        height, ia, ib = best
        merged = clusters[ia] + clusters[ib]
        history.append({
            "id_a": ia, "id_b": ib, "new_id": next_id,
            "members_a": list(clusters[ia]), "members_b": list(clusters[ib]),
            "height": height, "size": len(merged),
            "labels": _labels_now(clusters, n),   # cluster id per point, pre-merge
        })
        Z.append([ia, ib, height, len(merged)])

        del clusters[ia]
        del clusters[ib]
        clusters[next_id] = merged
        next_id += 1

    return np.array(Z, dtype=float), history
# endregion


def _labels_now(clusters, n):
    """A cluster id per original point given the current set of live clusters."""
    lab = np.empty(n, dtype=int)
    for cid, members in clusters.items():
        for m in members:
            lab[m] = cid
    return lab.tolist()


# region: cut_tree
def cut_tree(Z, n, k):
    """Cut the merge tree to get k flat clusters.

    Every row of Z is a merge; the first N-k of them are the merges that happen
    below the cut line. Replay exactly those with a union-find, and whatever
    connected components remain are the k flat clusters. Labels are remapped to
    a dense 0..k-1 so they line up with what a library returns.
    """
    parent = list(range(2 * n - 1))

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x

    for row, (a, b, _h, _s) in enumerate(Z[:n - k]):
        new = n + row
        parent[int(a)] = new
        parent[int(b)] = new

    roots = [find(i) for i in range(n)]
    order = {}
    labels = np.empty(n, dtype=int)
    for i, r in enumerate(roots):
        if r not in order:
            order[r] = len(order)
        labels[i] = order[r]
    return labels
# endregion


def load_data(path="../data/blobs.csv"):
    """make_blobs snapshot: 30 2-D points (x, y) plus the true blob id.

    The blob column is ground truth we NEVER cluster on — agglomerate sees only
    the coordinates. It exists so we can score the flat clustering afterwards.
    """
    return pd.read_csv(path)

The library version

Nobody hand-rolls this in production, and once you've built it you don't need to. scikit-learn's AgglomerativeClustering is the same bottom-up merging with a C inner loop and the same choice of linkage, and it takes k directly as n_clusters:

def sklearn_agglomerative(X, k, linkage="average"):
    """Cluster X into k groups with scikit-learn's agglomerative clustering.

    Same linkage rule as our from-scratch build, Euclidean metric. Returns the
    flat labels at the k-cut — the one thing the library gives you directly.
    """
    model = AgglomerativeClustering(n_clusters=k, linkage=linkage, metric="euclidean")
    return model.fit_predict(X)

Notice what it returns: flat labels, and nothing else. That's the library making a choice for you — it runs the merges internally and hands back only the cut at the k you asked for, throwing the tree away unless you go out of your way to keep it. Our version returns the whole tree because the tree is the interesting part; the library optimizes for the common case where you already know k and just want the groups. To validate our merge heights against a trusted reference we borrow scipy's linkage, which produces a linkage matrix in the same format ours does:

def scipy_reference(X, linkage="average"):
    """scipy's linkage matrix, used only to check our recorded merge heights.

    Same [id_a, id_b, height, size] rows our agglomerate() produces, computed by
    a battle-tested implementation. We compare the sorted merge heights against
    these to prove the from-scratch tree is the real tree.
    """
    return scipy_linkage(X, method=linkage, metric="euclidean")

We don't cluster with scipy — we only use it to prove our recorded heights are the real heights, by sorting both and comparing. They agree to 2.2e-16, which is machine epsilon: our from-scratch tree is byte-for-byte the same tree scipy builds.

Scratch versus library

Same data, same average linkage, cut at k=3 — our from-scratch flat labels against scikit-learn's. Since there's no inertia to compare, the honest metric is agreement: the adjusted Rand index, which scores how well two labelings match after correcting for chance, where 1.0 is a perfect match and 0 is random.

All three bars are pinned at 1.0. Our clustering and sklearn's agree perfectly, and both recover the true blobs perfectly — every point landed in the group it was really drawn from, without either implementation ever seeing a label. On three clean, well-separated blobs there's only one sensible way to cut the tree at k=3, both trees are identical, and the cut lands in the same place. The value of the library here isn't a better answer; it's the C loop that finds the same answer on data far bigger than thirty points.

And this is the payoff of building the tree instead of just the labels: we don't have to trust the cut, we can see it. Here's the full dendrogram, colored by the three clusters below the cut with the between-blob merges in gray above it, and the dashed line marking where k=3 falls — at height 3.92, dropped into the empty gap between the last within-blob merge at 1.90 and the first between-blob merge at 5.94:

The tall gray brackets stranded at the top are the whole story: everything below 1.90 is a merge inside a blob, everything above 5.94 is a merge between blobs, and the red line cuts through the empty band between them. Any cut in that band gives three clusters. That gap is what k-means never showed you — it would have taken your k and gone to work, but it would never have told you that three was sitting there in the data, obvious, the moment you drew the tree.

Takeaways

Reach for hierarchical clustering when you don't know k and don't want to pretend you do, or when the thing you're clustering is genuinely nested and the tree is the answer you actually want — a taxonomy, a topic hierarchy, a family of related documents. Build it once, read the merge heights, cut below the biggest gap. On small-to-medium data where you can afford the $n^2$ distance matrix, it's the honest first move for exploration, because it shows you the structure before it makes you commit to a number.

Reach past it the moment size becomes the issue. The quadratic memory and cubic time are hard walls, not constants you can tune away, and past a few thousand points you either sample down to a size the tree can handle or you go back to k-means, which stays linear and doesn't care how big the data gets. And whatever you do, name your linkage — single, complete, and average are three different algorithms wearing the same name, and single linkage in particular will happily chain two real clusters together through a bridge of stray points, so a clustering that looks wrong is often a linkage that was wrong for the data. Between this chapter and the last you now have the two clustering instincts that cover most of what you'll meet: k-means when you know how many groups and the groups are round, hierarchical when you don't and they might not be. Almost everything fancier is a variation on one of these two — which is the best reason to have built both by hand.