Hierarchical Clustering

The procedure of the previous lesson asked for the number of clusters before starting, and every run produced a single partition. When three different searches were made for three different values of k, the three results were unrelated to each other: the clusters of the three-cluster partition were not a union of the clusters of the five-cluster partition. Since three readings of the elbow chose three different values of k, this meant three unrelated answers at a point where no decision could be made.

This lesson reverses the order. Every subscriber starts as its own cluster, at every step the two nearest clusters are merged, and a nested family is produced from 200 subscribers down to a single cluster. The number of clusters is no longer an input; it is the answer to the question of which height the tree is cut at. In return a new decision appears, and it is what the lesson measures: the phrase “two nearest clusters” cannot be derived from the point distance, because the distance between two clusters must be defined separately.

• CL16. The dataset, columns and similarity definition are as in the previous lesson; the scale is again learned from the training set of 756 subscribers. Clustering, however, is done on the first 200 subscribers: the cost of the merging procedure grows with the cube of the number of points. This is a scale decision, and the numbers here are for these 200 subscribers, not for 756.
• CL17. The point distance is the square root of the squared distance in standard scores; the square root is taken so that the merge heights can be read in the unit of the columns.
• CL18. The linkage criterion is the definition of the distance between two clusters. Three are tried: single linkage takes the nearest pair of points, complete linkage the farthest pair of points, average linkage the average over all pairs of points.
• CL19. The merge height is the linkage distance between the two clusters merged at that step; it is specific to the criterion and cannot be compared across criteria.
• CL20. The procedure is deterministic; there is no choice of starting points, and if two candidates come out at equal distance, the pair with the smaller numbers is chosen.
• CL21. The measure is again the within-cluster sum, the baselines are random assignment and single cluster. To read k clusters from the tree, one starts at the root and opens the highest branch until k branches remain.
• CL22. The disagreement of two partitions is given by two measures: the ratio of disagreeing point pairs and the number of points that land apart under the best cluster matching; for the second, the 24 matchings of four clusters are tried.
• CL23. The measure is the number that the procedure of the previous lesson tries to minimize directly; hierarchical procedures do not target it. The comparison is therefore in favor of the measure.

The Tree Built by Merging

The skeleton of the procedure is short: the distances between clusters are kept, the smallest is found, two clusters are merged, and the distance of the new cluster to the others is computed from the old ones by the rule of the linkage criterion. This last step makes the procedure cheaper — after a merge, no pair of points is looked at again. The three criteria differ only in that single line.

# hierarchical.py — MODEL. The same FICTIONAL subscriber table, the same seed, the same standard
# score definition; the label is still not generated.
import itertools
import math
import statistics

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [(0.28, 21), (0.22, 17), (0.18, 26), (0.14, 14), (0.18, 23)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


DATA = []
for i in range(RAW):
    r = rng(SEED + i)
    base = REGION[pick(r(), [x[0] for x in REGION])][1]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() < 0.046:
        continue
    r, v = rng(SEED + 17001 + i), []
    for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            base * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - d * 0.05) * 100 + 0.5) / 100)
        r()
    if v:
        x = {"mean_consumption": round(sum(v) / len(v), 2), "volatility": round(max(v) - min(v), 2),
             "household": household, "satisfaction": rating, "period": len(v)}
        x["per_person"] = round(x["mean_consumption"] / household, 3)
        DATA.append(x)

FIELDS = ["mean_consumption", "volatility", "per_person", "household", "satisfaction", "period"]


def shuffle(data, seed):
    r, s = rng(seed), list(range(len(data)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return [data[i] for i in s]


TRAIN = shuffle(DATA, SEED + 90000)[:756]
RAW_VALUES = [[x[a] for a in FIELDS] for x in TRAIN]
MEANS = [statistics.fmean(c) for c in zip(*RAW_VALUES)]
STD_DEVS = [statistics.pstdev(c) for c in zip(*RAW_VALUES)]
V = [[(v[j] - MEANS[j]) / STD_DEVS[j] for j in range(len(FIELDS))] for v in RAW_VALUES]

SUBSET = 200
W = V[:SUBSET]
squared = lambda a, b: sum((p - q) * (p - q) for p, q in zip(a, b))
centroid = lambda g: [statistics.fmean(c) for c in zip(*g)]
dist = lambda a, b: math.sqrt(squared(a, b))

LINKAGE = {"single": lambda da, db, na, nb: min(da, db),
           "complete": lambda da, db, na, nb: max(da, db),
           "average": lambda da, db, na, nb: (na * da + nb * db) / (na + nb)}


def merge(rule):                            # MODEL: merge the two nearest clusters
    D = {(i, j): dist(W[i], W[j]) for i in range(SUBSET) for j in range(i + 1, SUBSET)}
    pair = lambda x, y: (x, y) if x < y else (y, x)
    active, size, children, height, new = list(range(SUBSET)), {i: 1 for i in range(SUBSET)}, {}, {}, SUBSET
    while len(active) > 1:
        a, b = min(((x, y) for i, x in enumerate(active) for y in active[i + 1:]),
                   key=lambda p: (D[pair(*p)], p))
        children[new], height[new] = (a, b), D[pair(a, b)]
        for c in active:
            if c not in (a, b):
                D[pair(c, new)] = rule(D[pair(a, c)], D[pair(b, c)], size[a], size[b])
        size[new] = size[a] + size[b]
        active = [c for c in active if c not in (a, b)] + [new]
        new += 1
    return children, height, size, new - 1


TREE = {name: merge(rule) for name, rule in LINKAGE.items()}
print(f"subset {SUBSET} subscribers; every linkage criterion made {SUBSET - 1} merges")
for name, (_, height, _, _) in TREE.items():
    h = sorted(height.values())
    print(f"  {name:<9} smallest merge {h[0]:.3f}, median {statistics.median(h):.3f}, "
          f"last three {h[-3]:.3f} {h[-2]:.3f} {h[-1]:.3f}")

subset 200 subscribers; every linkage criterion made 199 merges
  single    smallest merge 0.079, median 0.896, last three 1.977 2.173 2.245
  complete  smallest merge 0.079, median 1.211, last three 6.370 6.791 8.917
  average   smallest merge 0.079, median 1.086, last three 4.536 4.715 4.994

All three criteria start from the same place. The first merge is 0.079 in all three; between two clusters that each contain a single point, the nearest, the farthest and the average pair are the same pair, so the criteria cannot differ on single-element clusters. The difference begins as the clusters grow, and there is a gulf at the top: the last merge is 2.245 in single linkage, 4.994 in average linkage, 8.917 in complete linkage. Height is the unit of the criterion, not of the data, and the heights of two trees cannot be read side by side.

The Dendrogram and the Cut Height

The tree is printed as text: each line gives a branch, the indentation gives the depth of the branch, h gives the merge height, n gives the number of subscribers under the branch. It is limited to three levels; deeper branches are summarized in a single line.

def dendrogram(name, top=3, min_size=15):
    children, height, size, root = TREE[name]
    print(f"dendrogram — {name} linkage, top {top} levels")
    stack = [(root, 0)]
    while stack:
        node, depth = stack.pop()
        if node in children and size[node] >= min_size and depth < top:
            print(f"  {'  ' * depth}h={height[node]:.3f}  n={size[node]}")
            stack.extend((c, depth + 1) for c in sorted(children[node], key=lambda c: size[c]))
        else:
            print(f"  {'  ' * depth}branch    n={size[node]}")


dendrogram("average")
dendrogram("single")


def cluster_count(name, h):                 # how many merges remain above height h
    _, height, _, _ = TREE[name]
    return sum(1 for y in height.values() if y > h) + 1


print(f"\n{'cut height':>18}" + "".join(f"{name:>12}" for name in LINKAGE))
for h in (1.0, 1.5, 2.0, 3.0, 4.0, 5.0):
    print(f"{h:>18.1f}" + "".join(f"{cluster_count(name, h):>12}" for name in LINKAGE))

dendrogram — average linkage, top 3 levels
  h=4.994  n=200
    h=4.715  n=198
      h=3.704  n=181
        branch    n=140
        branch    n=41
      h=4.536  n=17
        branch    n=15
        branch    n=2
    branch    n=2
dendrogram — single linkage, top 3 levels
  h=2.245  n=200
    h=2.173  n=199
      h=1.977  n=198
        branch    n=197
        branch    n=1
      branch    n=1
    branch    n=1

        cut height      single    complete     average
               1.0          62         123         114
               1.5          16          80          59
               2.0           3          56          37
               3.0           1          26          10
               4.0           1          13           4
               5.0           1           9           1

The two diagrams show two different trees. In average linkage the root splits 200 into 198 and 2, then 181 and 17, then 140 and 41; at every step there is a split with both sides populated. In single linkage, on the other hand, a single subscriber breaks off at every step: from 200 to 199, from that to 198, from that to 197. The reason is in the definition of the criterion — because the nearest pair is what counts, every point that lies close to a large cluster is added to that cluster. The tree turns into a chain, and when it is cut from the top, what it yields is not a partition but one cluster plus a few single points.

The lower table shows what the cut height is. Cutting at a height of 2.0 gives 3 clusters in single linkage, 37 in average linkage, 56 in complete linkage; a single number corresponds to a number of clusters that varies eighteenfold across three criteria. At 5.0, single and average linkage give 1 cluster while complete linkage still counts 9. The cut height is therefore not an observation but a decision; the claim of reading the number of clusters from the tree does nothing but move the decision to another place.

The Partition That the Linkage Criterion Moves

The three trees are opened to four clusters and compared with the same measure and the same baseline. Next to them, the centroid-based procedure of the previous lesson is placed on the same 200 subscribers.

def cut_k(name, k):                         # starting at the root, open the highest branch until k
    children, height, _, root = TREE[name]
    branches = [root]
    while len(branches) < k:
        d = max((x for x in branches if x in children), key=lambda x: height[x])
        branches = [x for x in branches if x != d] + list(children[d])
    a = [0] * SUBSET
    for c, d in enumerate(branches):
        stack = [d]
        while stack:
            x = stack.pop()
            if x in children:
                stack.extend(children[x])
            else:
                a[x] = c
    return a


def within_sum(a, k):                       # within-cluster sum, per point
    t = 0.0
    for c in range(k):
        g = [W[i] for i in range(SUBSET) if a[i] == c]
        if g:
            t += sum(squared(v, centroid(g)) for v in g)
    return t / SUBSET


def disagreement(a, b):
    ay = sum((a[i] == a[j]) != (b[i] == b[j]) for i in range(SUBSET) for j in range(i + 1, SUBSET))
    return 2 * ay / (SUBSET * (SUBSET - 1))


def mismatched_points(a, b, k):             # points that disagree under the best cluster matching
    return SUBSET - max(sum(p[a[i]] == b[i] for i in range(SUBSET))
                        for p in itertools.permutations(range(k)))


def k_means(k, seed, runs=12):              # the procedure of lesson 02, on the same subset
    r, best = rng(seed), None
    for _ in range(runs):
        s = []
        while len(s) < k:
            j = int(r() * SUBSET)
            if j not in s:
                s.append(j)
        M, a = [W[j] for j in s], None
        for _round in range(60):
            y = [min(range(k), key=lambda c: squared(v, M[c])) for v in W]
            if y == a:
                break
            a = y
            M = [centroid([W[i] for i in range(SUBSET) if a[i] == c]) or M[c] for c in range(k)]
        u = within_sum(a, k)
        best = u if best is None else min(best, u)
    return best


def random_assignment(k, seed):             # BASELINE: uses no structure at all
    r = rng(seed)
    return [int(r() * k) for _ in range(SUBSET)]


K = 4
P = {name: cut_k(name, K) for name in LINKAGE}
baseline = statistics.fmean(within_sum(random_assignment(K, SEED + 33300 + t), K) for t in range(5))
print(f"k = {K}, {SUBSET} subscribers")
print(f"  {'baseline - single cluster':<28} {within_sum([0] * SUBSET, K):.4f}")
print(f"  {'baseline - random assignment':<28} {baseline:.4f}")
for name in LINKAGE:
    u = within_sum(P[name], K)
    print(f"  {name + ' linkage':<28} {u:.4f}   difference {baseline - u:+.4f}"
          f"   cluster sizes {sorted(P[name].count(c) for c in range(K))}")
print(f"  {'k-means (12 runs)':<28} {k_means(K, SEED + 44400):.4f}")
print("\ndisagreement between linkage criteria")
A = list(LINKAGE)
for i in range(3):
    for j in range(i + 1, 3):
        print(f"  {A[i]:<9} - {A[j]:<9} point pairs {disagreement(P[A[i]], P[A[j]]):.4f}, "
              f"points apart {mismatched_points(P[A[i]], P[A[j]], K):>3} / {SUBSET}")

k = 4, 200 subscribers
  baseline - single cluster    6.3473
  baseline - random assignment 6.2340
  single linkage               6.0849   difference +0.1492   cluster sizes [1, 1, 1, 197]
  complete linkage             4.5419   difference +1.6922   cluster sizes [6, 16, 65, 113]
  average linkage              5.0080   difference +1.2260   cluster sizes [2, 2, 15, 181]
  k-means (12 runs)            3.4245

disagreement between linkage criteria
  single    - complete  point pairs 0.5451, points apart  86 / 200
  single    - average   point pairs 0.1477, points apart  17 / 200
  complete  - average   point pairs 0.4673, points apart  82 / 200

The number for single linkage is 6.0849, and it beats the random assignment baseline by only 0.1492. The cluster sizes give the reason: [1, 1, 1, 197]. Four clusters were asked for, three of them are a single subscriber; the result is too close to a procedure that does not learn to be told apart from it. Complete linkage beats the baseline by 1.6922 with 4.5419 and its sizes are [6, 16, 65, 113] — the most balanced partition is its, because a criterion that looks at the farthest pair penalizes wide clusters. Average linkage is between the two: 5.0080, difference 1.2260. The choice of linkage criterion moves the difference over the baseline between 0.1492 and 1.6922, that is, by more than elevenfold.

The lower table counts the same thing at the partition level. Single and complete linkage disagree on 0.5451 of the point pairs, and even under the best matching 86 subscribers land in different places — 43% of the 200 subscribers. The disagreement between single and average linkage is much smaller, 17 subscribers; average linkage also leans toward one large cluster, so it stays related to single linkage. Complete and average linkage, in turn, disagree on 82 subscribers.

The 3.4245 in the last row must be read carefully. The centroid-based procedure of the previous lesson beats all three by a clear margin at this measure, but the measure is the number that it tries to minimize directly; hierarchical procedures do not even compute this sum at any step. The comparison shows not the worth of the procedures but which side the measure stands on. Had we chosen the consistency of the linkage distances or resistance to chaining as the measure, the ranking would have reversed.

Summary

• The merging procedure turns the number of clusters from an input into a nested family of clusters, and the decision moves to the cut height.
• The three linkage criteria start from the same first merge (0.079) but their roots are 2.245, 4.994 and 8.917; height is the unit of the criterion and cannot be compared across criteria.
• Cutting at a height of 2.0 gives 3 clusters in single linkage, 37 in average linkage, 56 in complete linkage.
• At four clusters, single linkage beats the random assignment baseline by only 0.1492 with 6.0849 and produces [1, 1, 1, 197]; complete linkage beats it by 1.6922 with 4.5419. The choice of criterion moves the difference by more than elevenfold.
• Single and complete linkage place 86 of the 200 subscribers in different places even under the best matching; the 3.4245 of the centroid-based procedure is read together with the fact that the measure is its target.

Next Step

In three lessons, four procedures were tried — single-pass centroid selection, iterative centroid update and three linkage criteria — and all of them shared a single assumption: every point must be assigned to a cluster. The [1, 1, 1, 197] partition of single linkage shows the price of this assumption. Three subscribers stand on their own, but the procedure had to declare them “clusters”, because it has no other box. In the same way, the points that stand at the edge of the 113-subscriber cluster of complete linkage and are not really close to any centroid were also written into a cluster. There was no such thing as a point that belongs to no cluster. The next lesson takes up a family that removes this assumption — a procedure that looks at the density of the points, leaves the points in sparse regions outside as noise, and never asks for the number of clusters — and the number of points left outside becomes a measure in itself.