Clustering Evaluation

Over five lessons every method asked for a number, always given from the outside: how many clusters, which cut height, which neighborhood radius, how many components. Each time it was measured that the choice moved the result, but how to make the choice itself was postponed. In the Supervised Learning course this question was easy: one looked at the labels in the validation set and took the candidate with the best number. Here there is no label, so the measure must be built from the data.

This lesson writes three internal validity measures and applies them to the same data. All three are reasonable, all three turn a common intuition into a number, and the three choose three different numbers of clusters. Then the measures are run on a cloud with no structure; they choose a number of clusters there as well. At the very end the label is opened, for the first and only time in this course, for external validity, and compared with the internal measures’ choices.

• CL45. The dataset is the training portion of 756 rows of the fictional subscriber table from the Introduction to Machine Learning, Data Preparation and Feature Engineering, and Supervised Learning courses, seed 20260218, six scaled numeric columns. The label is generated but enters no computation until the third section.
• CL46. The partition is the centroid-based procedure of the second lesson of the topic and is run for every number of clusters with the same seed; the partitioning procedure is held fixed so that the measures are comparable.
• CL47. Silhouette coefficient: for a point, take the mean distance $a$ to the other points in its own cluster and the mean distance $b$ to the points in the nearest other cluster, and compute $(b-a)/\max(a,b)$. The mean over the points is taken, and the larger is chosen.
• CL48. Scatter ratio: the mean distance of a cluster to its centroid is counted as the radius of that cluster; for a pair of clusters the sum of the two radii is divided by the distance between the centroids, the worst partner is taken for each cluster, and the mean is computed. The smaller is chosen.
• CL49. Gap measure: the natural logarithm of the total distance of the points to their centroids is subtracted from the mean of the same quantity over three reference clouds that carry no structure. The larger is chosen.
• CL50. The tried numbers of clusters are nine candidates from 2 to 10; with three measures, twenty-seven comparisons are made. The Statistical Analysis and Exploratory Data Analysis course measured that when many questions are asked of the same data the best result is a choice; those numbers are not repeated, and the rule stands.
• CL51. The structureless cloud is 756 points whose columns are generated independently of each other. Its difference from the within-column shuffling of the earlier lessons is that there the real distributions are preserved.
• CL52. The agreement index is built from its own definition: all pairs of points are visited, and the share of pairs for which the “same cluster?” answer of the partition agrees with the “same class?” answer of the label is taken. 756 points make 285,390 pairs.
• CL53. Two baselines are used in external validity: a random partition into the same number of clusters, and a single cluster that puts everyone into one cluster. The second is necessary, because when the label is imbalanced a single cluster produces a high agreement rate.

Three Measures, Three Different Numbers of Clusters

All three measures are forms of the same intuition: if a partition is good, points in the same cluster should be close to each other and points in different clusters far apart. They part ways in how they balance “close” against “far” and how they penalize growth in the number of clusters.

# validity.py — MODEL. The FICTIONAL subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering and Supervised Learning courses; the label column
# is generated but is NEVER used in this block or the next.
import math
import statistics

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


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


SUBSCRIBERS, DATA = [], []
for i in range(RAW):
    r = rng(SEED + i)
    region = REGION[pick(r(), [x[0] for x in REGION])]
    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:
        SUBSCRIBERS.append((10001 + i, region, household, rating))
for sub_id, region, household, rating in SUBSCRIBERS:
    r, v = rng(SEED + 7000 + sub_id), []
    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(
            region[1] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - d * 0.05) * 100 + 0.5) / 100)
        r()
    if not v:
        continue
    avg, spread, q = round(sum(v) / len(v), 2), round(max(v) - min(v), 2), rng(SEED + 51000 + sub_id)
    z = (0.052 * (sum(v) / len(v) - 20) + 0.85 * (min(v) == 0.0) + 0.026 * spread
         + 0.24 * (household >= 5) - 0.20 * (rating >= 4) + region[2]
         + 0.9 * (sub_id - 10001) / 1399 + (q() + q() + q() - 1.5) * 1.30)
    DATA.append({"mean_consumption": avg, "volatility": spread, "household": household,
                 "satisfaction": rating, "period": len(v), "per_person": round(avg / household, 3), "suspicious": int(z > 1.35)})

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


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]
SCALE = {a: (statistics.fmean(x[a] for x in TRAIN),
             statistics.pstdev([x[a] for x in TRAIN])) for a in NUMERIC}
N = [[(x[a] - SCALE[a][0]) / SCALE[a][1] for a in NUMERIC] for x in TRAIN]
LABEL = [x["suspicious"] for x in TRAIN]          # not opened until the third section
K_RANGE = range(2, 11)


def distance_matrix(P):
    D = [[0.0] * len(P) for _ in P]
    for i in range(len(P)):
        for j in range(i + 1, len(P)):
            D[i][j] = D[j][i] = math.dist(P[i], P[j])
    return D


def partition(P, k, seed, rounds=40):             # MODEL: centroid-based partition
    r = rng(seed)
    M, e = [list(P[int(r() * len(P))]) for _ in range(k)], [0] * len(P)
    for _ in range(rounds):
        for i, x in enumerate(P):
            e[i] = min(range(k), key=lambda c: math.dist(x, M[c]))
        for c in range(k):
            g = [P[i] for i in range(len(P)) if e[i] == c]
            M[c] = [statistics.fmean(x[d] for x in g) for d in range(len(P[0]))] if g else M[c]
    return e, M


def silhouette(D, e, k):                          # MEASURE 1: larger is chosen
    G = [[i for i in range(len(e)) if e[i] == c] for c in range(k)]
    t = 0.0
    for i in range(len(e)):
        o = e[i]
        if len(G[o]) < 2:
            continue
        a = sum(D[i][j] for j in G[o]) / (len(G[o]) - 1)
        b = min(sum(D[i][j] for j in G[c]) / len(G[c]) for c in range(k) if c != o and G[c])
        t += (b - a) / max(a, b)
    return t / len(e)


def scatter(P, e, M, k):                          # MEASURE 2: smaller is chosen
    R = [statistics.fmean([math.dist(P[i], M[c]) for i in range(len(P)) if e[i] == c] or [0.0])
         for c in range(k)]
    return statistics.fmean(
        max((R[c] + R[j]) / max(math.dist(M[c], M[j]), 1e-9) for j in range(k) if j != c)
        for c in range(k))


def within_total(P, e, M):
    return sum(math.dist(P[i], M[e[i]]) for i in range(len(P)))


def structureless(n, d, seed):                    # FICTIONAL: a cloud that carries no structure
    r = rng(seed)
    return [[(r() + r() + r() - 1.5) * 1.4142 for _ in range(d)] for _ in range(n)]


REFERENCE = [structureless(756, len(NUMERIC), SEED + 300000 + 1000 * b) for b in range(3)]


def gap(P, k, W):                                 # MEASURE 3: larger is chosen
    return statistics.fmean(math.log(within_total(R, *partition(R, k, SEED + 200000 + k)))
                            for R in REFERENCE) - math.log(W)


def sweep(P, D):
    out = {}
    for k in K_RANGE:
        e, M = partition(P, k, SEED + 200000 + k)
        out[k] = (silhouette(D, e, k), scatter(P, e, M, k), gap(P, k, within_total(P, e, M)))
    return out


def choice(t):
    return (max(t, key=lambda k: t[k][0]), min(t, key=lambda k: t[k][1]),
            max(t, key=lambda k: t[k][2]))


D = distance_matrix(N)
VT = sweep(N, D)
print(f"points {len(N)}, columns {len(NUMERIC)}, tried cluster counts {len(K_RANGE)}")
print(f"{'k':>4}{'silhouette':>12}{'scatter':>10}{'gap':>10}")
for k in K_RANGE:
    print(f"{k:>4}{VT[k][0]:>12.4f}{VT[k][1]:>10.4f}{VT[k][2]:>10.4f}")
print(f"chosen k -> silhouette {choice(VT)[0]}, scatter {choice(VT)[1]}, gap {choice(VT)[2]}")
print(f"compared candidates {len(K_RANGE) * 3}")

points 756, columns 6, tried cluster counts 9
   k  silhouette   scatter       gap
   2      0.1954    1.8924   -0.2894
   3      0.2142    1.5466   -0.2481
   4      0.1760    1.6043   -0.2288
   5      0.1799    1.4596   -0.2148
   6      0.1928    1.3328   -0.1882
   7      0.1838    1.3481   -0.1865
   8      0.1810    1.3773   -0.1856
   9      0.1781    1.3388   -0.1646
  10      0.1833    1.3947   -0.1520
chosen k -> silhouette 3, scatter 6, gap 10
compared candidates 27

Three measures, three answers: 3, 6 and 10. The same data, the same partitioning procedure, the same seed. Only the definition of a “good partition” changes, and the answer changes with it.

The size of the numbers must be read too. The highest value of the silhouette is 0.2142; the upper bound of the measure is 1, and 0.2142 means that the clusters largely overlap. The ranking between 3 and 6 rests on a difference of 0.0214, so the measure shows no peak, only a slight slope. The gap measure does not peak at all: it rises in a single direction from 2 to 10, and the value it chooses is the last candidate of the tried range. Its answer comes not from the data but from where the sweep was cut off.

On top of that, the best of twenty-seven comparisons is taken. The Statistical Analysis and Exploratory Data Analysis course measured that when many questions are asked of the same data the highest result is a choice, not a finding; the rule applies here to the number of clusters as it is.

The Measures on a Structureless Dataset

If a measure produces only a number, one needs to know when that number says “there is no structure”. All three measures always find a best candidate; none has an output saying “there is no cluster here”. The way to see this is to run the measures on a cloud with no structure.

FLAT = structureless(756, len(NUMERIC), SEED + 777000)
YT = sweep(FLAT, distance_matrix(FLAT))
print(f"{'k':>4}{'silhouette':>12}{'scatter':>10}{'gap':>10}")
for k in K_RANGE:
    print(f"{k:>4}{YT[k][0]:>12.4f}{YT[k][1]:>10.4f}{YT[k][2]:>10.4f}")
print(f"k chosen on the structureless cloud -> silhouette {choice(YT)[0]}, scatter {choice(YT)[1]}, "
      f"gap {choice(YT)[2]}")
print(f"highest silhouette: data {VT[choice(VT)[0]][0]:.4f}, structureless {YT[choice(YT)[0]][0]:.4f}, "
      f"difference {VT[choice(VT)[0]][0] - YT[choice(YT)[0]][0]:+.4f}")

   k  silhouette   scatter       gap
   2      0.1169    2.6101   -0.0336
   3      0.1093    2.3145   -0.0405
   4      0.1136    2.0117   -0.0374
   5      0.1151    1.8629   -0.0408
   6      0.1142    1.8107   -0.0388
   7      0.1252    1.6415   -0.0360
   8      0.1135    1.7448   -0.0417
   9      0.1255    1.6649   -0.0301
  10      0.1201    1.6989   -0.0331
k chosen on the structureless cloud -> silhouette 9, scatter 7, gap 9
highest silhouette: data 0.2142, structureless 0.1255, difference +0.0887

There is no cluster in this cloud. The six columns were generated independently of each other, with no region of concentration. Still the three measures give three answers: 9, 7 and 9. None says “no structure”, because none has such an output — every measure has to point to the best among the candidates it is given.

Comparing the values is even more uncomfortable. The highest silhouette on the structureless cloud is 0.1255 and on the real data 0.2142; the difference is 0.0887, a ratio of about 1.7 times. This shows a nonzero structure, but also that a 0.21 read on its own does not mean “distinct clusters”. That the gap measure stays between -0.0301 and -0.0417 on the structureless cloud is expected: the structureless data does not separate from the reference clouds. That the same measure also stays negative on the real data is a separate warning; the centroid-based partition leaves a longer total distance on this dataset than on the reference clouds.

Opening the Label

So far in this course no method has seen the label. Here, once, it is opened. The aim is not to train a model with it but to ask whether the number of clusters chosen by the internal measures agrees with the one the label points to. The measure is the agreement index: all pairs of points are visited, and the pairs for which the “same cluster” answer of the partition agrees with the “same class” answer of the label are counted.

def agreement(e, y):                              # EXTERNAL VALIDITY: pair agreement rate
    n, same = len(y), 0
    for i in range(n):
        for j in range(i + 1, n):
            same += (e[i] == e[j]) == (y[i] == y[j])
    return same / (n * (n - 1) // 2)


def random_assignment(n, k, seed):                # BASELINE: random partition
    r = rng(seed)
    return [int(r() * k) for _ in range(n)]


print(f"label: {sum(LABEL)} of {len(LABEL)} subscribers are suspicious "
      f"({sum(LABEL) / len(LABEL):.4f}); number of clusters the label points to 2")
print(f"{'k':>4}{'agreement':>10}{'random':>10}{'difference':>12}")
AG = {}
for k in K_RANGE:
    e, _ = partition(N, k, SEED + 200000 + k)
    AG[k] = (agreement(e, LABEL),
             agreement(random_assignment(len(LABEL), k, SEED + 900000 + k), LABEL))
    print(f"{k:>4}{AG[k][0]:>10.4f}{AG[k][1]:>10.4f}{AG[k][0] - AG[k][1]:>+12.4f}")
SINGLE = agreement([0] * len(LABEL), LABEL)
best = max(AG, key=lambda k: AG[k][0])
print(f"single cluster baseline {SINGLE:.4f}; partitions that exceed it "
      f"{sum(1 for k in K_RANGE if AG[k][0] > SINGLE)}/{len(K_RANGE)}")
print(f"partition with the highest agreement k={best} ({AG[best][0]:.4f}), "
      f"{AG[best][0] - SINGLE:+.4f} below the single cluster baseline")
print(f"k chosen by the internal measures {list(choice(VT))}; k the label points to 2")

label: 168 of 756 subscribers are suspicious (0.2222); number of clusters the label points to 2
   k agreement    random  difference
   2    0.5235    0.4996     +0.0239
   3    0.4690    0.4484     +0.0206
   4    0.4290    0.4225     +0.0065
   5    0.4195    0.4073     +0.0121
   6    0.4209    0.3967     +0.0241
   7    0.4051    0.3898     +0.0153
   8    0.4044    0.3840     +0.0204
   9    0.3936    0.3804     +0.0132
  10    0.3843    0.3778     +0.0065
single cluster baseline 0.6539; partitions that exceed it 0/9
partition with the highest agreement k=2 (0.5235), -0.1304 below the single cluster baseline
k chosen by the internal measures [3, 6, 10]; k the label points to 2

The label has two classes: 168 of 756 subscribers are suspicious, a rate of 0.2222. The label therefore points to 2 clusters. The numbers chosen by the internal measures are 3, 6 and 10; none of them said 2.

Against the random partition baseline the difference is positive at every number of clusters but small: the largest is 0.0241. The partition agrees with the label slightly more than random does. The real number is in the second baseline. The procedure that puts everyone in a single cluster gets 0.6539, and all nine partitions stay below it; the best one is 0.1304 behind. The reason is the imbalance of the label: at a suspicious rate of 0.2222, two randomly chosen subscribers are already likely to be in the same class, and saying “all the same” gets all of those pairs right.

The conclusion is not that clustering does not work but that it does not find the label — and there was no reason for it to, since no method saw the label. The clusters separated by consumption level, volatility and household size; the suspicion label overlaps with none of these axes. Measuring an unsupervised result against an external label often measures not the clustering’s success but whether the label is visible in the data.

Summary

• Three internal validity measures choose three different numbers of clusters on the same data: the silhouette coefficient 3, the scatter ratio 6, the gap measure 10. The only thing that changes is the definition of a “good partition”.
• The gap measure does not peak, it chooses the last candidate of the tried range; its answer comes not from the data but from where the sweep was cut off. The best of twenty-seven comparisons was taken.
• On a cloud that carries no structure the same three measures still choose a number of clusters (9, 7 and 9); the highest silhouette there is 0.1255, on the real data 0.2142, a difference of 0.0887.
• The label was opened for the first and only time in the course: the agreement index is 0.5235 for the best partition, 0.0239 above the random partition baseline but 0.1304 below the single cluster baseline (0.6539); all nine partitions stay below this baseline.
• The numbers of clusters chosen by the internal measures (3, 6, 10) do not agree with the one the label points to (2); measuring an unsupervised result with an external label is a measure not of the clustering but of the visibility of the label in the data.

Next Step

In six lessons points were grouped, and every group was read beside a baseline. The methods differed greatly from each other — centroid, linkage, density, mixture — but they shared one thing: all of them looked at all the columns of a point at once and read the distance between two points from the sum of these columns. As the number of columns grows, what distance itself is was never asked; it was assumed that six columns and three hundred columns produce the same operation in the same sense. Yet as the number of columns grows, distances approach one another and the word “nearest” loses its meaning. The next topic opens with this question: is it possible to show the points with fewer numbers, and what is lost in this representation. The measure is no longer the partition but the reconstruction error — how far the points rebuilt from a few numbers stay from their original places.