--- title: 'Unsupervised Learning' source: 'https://academia.sh/en/courses/unsupervised-learning' language: en updated: '2026-10-21T05:03:31+00:00' license: 'CC BY-SA 4.0' lessons: 10 outcomes: - 'Chooses clustering methods according to the structure of the data and builds the baseline of each result itself' - 'Measures dimensionality reduction by reconstruction error and neighborhood preservation and writes down what is lost' - 'Evaluates unsupervised results without labels and shows that the choice of metric determines the answer' --- # Unsupervised Learning How the baseline must also be built on a set without labels, and measuring how far the result departs from a random partition: the definition of similarity producing three entirely different partitions from the same set, the result that initialization moves at the same cluster count, the number of points the linkage criterion separates, the number of noise points in the density-based approach, the uncertainty that soft assignment carries and that vanishes when switching to hard assignment, three internal validity metrics choosing three different cluster counts, the reconstruction error of a projection and the geometry it preserves, the trade-off between rank and error, the neighborhood preservation rate of an embedding, how many of the flagged anomalies are real, and confidence alone being misleading. ## Clustering Counting how far a result of unlabeled grouping of points departs from a random partition: three entirely different partitions that three different similarity definitions draw from the same dataset, the different results that different starting points produce at the same number of clusters, and the uncertainty of the elbow, the number of points that single, complete and average linkage separate and the cut height being a decision, the number of clusters and noise points that the neighborhood radius and the minimum point count move, the uncertainty carried by soft assignment disappearing in hard assignment, three internal validity measures choosing three different numbers of clusters, and finding clusters even in a dataset with no structure. ### The Clustering Problem Source: https://academia.sh/en/courses/unsupervised-learning/clustering-problem Building the measure itself once the label is removed: the same 756 subscribers are split into four clusters with three different similarity definitions and three entirely different partitions come out — the partitions disagree on a share of point pairs between 0.3750 and 0.4801. Three baselines are defined; under the standard score definition a single cluster gives 6.0000, random assignment 5.9790, slicing by a single column 4.9044 and the method 3.8910, so the method's difference is 2.0880 or 1.0134 depending on which baseline is written. The numbers of the three definitions cannot be compared with each other; the proportional gains are 0.525, 0.349 and 0.295. The Supervised Learning course tried seven model families in twenty lessons on the same fictional subscriber table and wrote its budget next to every number. All of those measurements shared one foundation that was never questioned — the **label**. The baseline predicted the most frequent label, the measure compared the predicted label with the true one, the tuning search read the validation labels. Remove the label column and all three mechanisms stop working at once. This lesson starts there. There is a table without labels and the goal is to gather similar subscribers together. The trouble is not finding a method — producing a partition is easy, and random assignment produces one too. The trouble is **what** the partition is measured against. Before the measure is built, there is an even more basic decision: the definition of "similar". - **CL1.** A **cluster**, in this course, is a subgroup of points; the database cluster of the Databases curriculum was a group of servers, while the **dataset** is all the rows. - **CL2.** The dataset is the **fictional** subscriber table from the Introduction to Machine Learning, Data Preparation and Feature Engineering, and Supervised Learning courses: seed 20260218, 1,260 subscribers, split 756/252/252. The procedure works on the 756 subscribers of the training set. - **CL3.** **The suspicious column is not generated in this course.** No method sees the label; the label will be used only once, in the last lesson of the topic, when external validity is measured. - **CL4.** The columns are six numeric fields: mean consumption, volatility, per-person consumption, household size, satisfaction and period count. Region indicators are left out, because the "distance" between two regions is always the same constant in binary columns. - **CL5.** A **partition** is an arrangement that assigns every point to exactly one cluster. Its measure is the **within-cluster sum**: the sum of every point's distance to the **centroid** of its own cluster, divided by the number of points; smaller is better. The numbers printed in the application areas lesson of the Introduction to Machine Learning course are not repeated; they came from a different set of columns. - **CL6.** Three baselines are used: **random assignment**, **single cluster** and **slicing by a single column**. The number for random assignment is the mean of five seeds; a single seed hides its fluctuation. - **CL7.** The number of clusters is held at **4**, and so is the procedure: four starting points are chosen with farthest-point seeding and two refinement rounds are run. The number of clusters as a hyperparameter and the starting point moving the result belong to the next lesson; here only the **similarity definition** changes. - **CL8.** The numbers of different definitions **cannot be compared with each other**, because each measures its own distance; only the difference against its own baseline is comparable. #### Similarity Definition Is a Choice A clustering procedure does not learn what "close" means from the data; it takes it from outside. The first definition leaves the columns on the **raw scale** and sums squared differences: because the range of mean consumption is about eight and a half times as wide as the range of household size, the sum is quietly a weighting. The second turns every column into a **standard score**, so the six columns get equal weight. The third turns every value into the **percentage of its rank within the column** and sums absolute differences; extreme values are squeezed and only the ordering remains. The definition determines not only the distance but also **what the center is**: the center that minimizes the squared distance is the mean vector, and the one that minimizes the absolute difference is the median vector. The three definitions do not pick the same nearest neighbor even for a single subscriber. ```python # clustering.py — MODEL: the FICTIONAL subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering and Supervised Learning courses is rebuilt with the same seed; # the suspicious column is not generated in this course, no procedure sees the label. 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] N = len(TRAIN) 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)] STD_VALUES = [[(v[j] - MEANS[j]) / STD_DEVS[j] for j in range(len(FIELDS))] for v in RAW_VALUES] def rank_percent(column): # value turns into the percentage of its mean rank in the column s, place = sorted(range(len(column)), key=lambda i: column[i]), {} for p, i in enumerate(s): place.setdefault(column[i], []).append(p) mid = {v: statistics.fmean(p) / (len(s) - 1) for v, p in place.items()} return [mid[x] for x in column] RANK_COLUMNS = [rank_percent(c) for c in zip(*RAW_VALUES)] RANK_VALUES = [[RANK_COLUMNS[j][i] for j in range(len(FIELDS))] for i in range(N)] squared = lambda a, b: sum((p - q) * (p - q) for p, q in zip(a, b)) absolute = lambda a, b: sum(abs(p - q) for p, q in zip(a, b)) mean_centroid = lambda g: [statistics.fmean(c) for c in zip(*g)] median_centroid = lambda g: [statistics.median(c) for c in zip(*g)] DEFINITIONS = [("raw scale, squared distance", RAW_VALUES, squared, mean_centroid), ("standard score, squared distance", STD_VALUES, squared, mean_centroid), ("rank percentage, absolute difference", RANK_VALUES, absolute, median_centroid)] print(f"clustered subscribers {N}, columns {len(FIELDS)}, no label column") print(f"{'':<38}{'neighbor':>10}" + "".join(f"{a:>18}" for a in FIELDS)) print(f"{'subscriber 48 (itself)':<38}{'--':>10}" + "".join(f"{TRAIN[48][a]:>18.3f}" for a in FIELDS)) for name, V, dist, _ in DEFINITIONS: j = min((x for x in range(N) if x != 48), key=lambda x: (dist(V[48], V[x]), x)) print(f"{name:<38}{j:>10}" + "".join(f"{TRAIN[j][a]:>18.3f}" for a in FIELDS)) ``` ``` clustered subscribers 756, columns 6, no label column neighbor mean_consumption volatility per_person household satisfaction period subscriber 48 (itself) -- 16.090 10.560 3.218 5.000 4.000 3.000 raw scale, squared distance 711 16.280 10.600 4.070 4.000 4.000 2.000 standard score, squared distance 52 13.540 11.990 2.708 5.000 4.000 3.000 rank percentage, absolute difference 241 17.400 9.830 2.900 6.000 4.000 3.000 ``` The three rows give three different neighbors for the same subscriber, and the reason is visible in each row. Subscriber 711, chosen on the raw scale, is only 0.19 away from subscriber 48 in mean consumption; the household and period columns do not match, but their differences add at most 2 to the sum. Subscriber 52, chosen with the standard score, is 2.55 away in consumption, but the household, satisfaction and period columns all match: a perfect match in one column offsets a difference of one standard deviation. None of the three neighbors is wrong. All three fit "nearest"; only for three different definitions. #### Three Baselines The within-cluster sum cannot be interpreted alone; whether the number is large or small cannot be known without a comparison point from the same measure. Three baselines answer three questions. **Random assignment** produces the same number of clusters without any structure: "how much better is this partition than chance". **A single cluster** keeps the points together: "what did splitting gain over not splitting at all". **Slicing by a single column** splits the points into four equal slices by mean consumption: "what did looking at all the columns gain over a solution that looks at a single column". The third is the hardest baseline: it uses the data even though it does not learn. ```python CLUSTERS = 4 def within_sum(V, dist, centroid, assign): # within-cluster sum, per point t = 0.0 for c in set(assign): g = [V[i] for i in range(len(V)) if assign[i] == c] t += sum(dist(v, centroid(g)) for v in g) return t / len(V) def random_assignment(n, k, seed): # BASELINE 1: uses no structure r = rng(seed) return [int(r() * k) for _ in range(n)] def by_column(V, j, k): # BASELINE 3: splits a single column into k equal slices s, a = sorted(range(len(V)), key=lambda i: (V[i][j], i)), [0] * len(V) for p, i in enumerate(s): a[i] = min(k - 1, p * k // len(V)) return a def partition(V, dist, centroid, k, rounds=2): # MODEL: farthest-point seeding + two rounds starts = [min(range(len(V)), key=lambda i: (dist(V[i], centroid(V)), i))] while len(starts) < k: starts.append(max(range(len(V)), key=lambda i: (min(dist(V[i], V[c]) for c in starts), -i))) M = [V[c] for c in starts] for _ in range(rounds + 1): a = [min(range(k), key=lambda c: dist(v, M[c])) for v in V] M = [centroid([V[i] for i in range(len(V)) if a[i] == c]) for c in range(k)] return a _, V2, D2, C2 = DEFINITIONS[1] RND = [within_sum(V2, D2, C2, random_assignment(N, CLUSTERS, SEED + 400 + t)) for t in range(5)] A2 = partition(V2, D2, C2, CLUSTERS) single, rnd = within_sum(V2, D2, C2, [0] * N), statistics.fmean(RND) sliced = within_sum(V2, D2, C2, by_column(V2, 0, CLUSTERS)) method = within_sum(V2, D2, C2, A2) print(f"with the standard score, squared distance definition, k = 4\n" f" baseline - single cluster {single:8.4f}\n" f" baseline - random assignment {rnd:8.4f} (five seeds: {min(RND):.4f} - {max(RND):.4f})\n" f" baseline - mean_consumption slice {sliced:8.4f}\n" f" method {method:8.4f}" f" cluster sizes {sorted(A2.count(c) for c in range(CLUSTERS))}\n" f" difference: against random {rnd - method:+.4f}, against slice {sliced - method:+.4f}") ``` ``` with the standard score, squared distance definition, k = 4 baseline - single cluster 6.0000 baseline - random assignment 5.9790 (five seeds: 5.9763 - 5.9819) baseline - mean_consumption slice 4.9044 method 3.8910 cluster sizes [61, 67, 117, 511] difference: against random +2.0880, against slice +1.0134 ``` The number for the single cluster is **exactly 6.0000**, that is, the number of columns. This is no coincidence: a standard score carries unit variability, so the per-point sum of the distance to a single center equals the number of columns. The measure's upper bound is known by calculation, and every partition is read by how far down it pulls this 6.0000. The second row is the most important number of this course. Random assignment gives 5.9790, only 0.0210 below the single cluster. Four random groups gain almost nothing over not splitting, because the center of a random cluster falls on the overall center of the data. The fluctuation across the five seeds is 0.0056, negligible beside the measured differences. The third row changes the table. Slicing into four equal slices by mean consumption alone gives 4.9044. The method's number is 3.8910 and the difference is **2.0880 or 1.0134, depending on which baseline is written**. Had the slice not been written, the method would have credited itself with the 1.0134 that looking at all six columns gained. The choice of baseline is a measurement decision; when it is not written with its reasoning, the reported difference is not a finding but a choice. #### The Same Cluster, Three Partitions The third block measures what happens when the procedure and the number of clusters are held fixed and only the similarity definition changes. The disagreement of two partitions is counted over point **pairs**, without matching the clusters: 756 subscribers have 285,390 pairs, and for every pair each partition says either "in the same cluster" or "in different clusters". The disagreement is the ratio of the pairs on which they answer differently. ```python PARTITIONS = [] print(f"{'definition':<38}{'single cluster':>16}{'random':>16}{'slice':>16}" f"{'method':>16}{'share over random':>21}") for name, V, dist, centroid in DEFINITIONS: a = partition(V, dist, centroid, CLUSTERS) PARTITIONS.append((name, a)) single = within_sum(V, dist, centroid, [0] * N) rnd = statistics.fmean(within_sum(V, dist, centroid, random_assignment(N, CLUSTERS, SEED + 400 + t)) for t in range(5)) cut = within_sum(V, dist, centroid, by_column(V, 0, CLUSTERS)) y = within_sum(V, dist, centroid, a) print(f"{name:<38}{single:>16.4f}{rnd:>16.4f}{cut:>16.4f}{y:>16.4f}{(rnd - y) / rnd:>21.3f}") def disagreement(a, b): # on how many point pairs the two partitions differ differing = sum((a[i] == a[j]) != (b[i] == b[j]) for i in range(len(a)) for j in range(i + 1, len(a))) return 2 * differing / (len(a) * (len(a) - 1)) print("\npair disagreement between partitions (ratio of differing point pairs)") for i in range(3): for j in range(i + 1, 3): print(f" {PARTITIONS[i][0]:<38} - {PARTITIONS[j][0]:<38} {disagreement(PARTITIONS[i][1], PARTITIONS[j][1]):.4f}") print("\nmeans of the largest cluster") for name, a in PARTITIONS: g = [TRAIN[i] for i in range(N) if a[i] == max(a, key=a.count)] print(f" {name:<38} n {len(g):>4}" + "".join( f" {x}={statistics.fmean(k[x] for k in g):.2f}" for x in ("mean_consumption", "volatility"))) ``` ``` definition single cluster random slice method share over random raw scale, squared distance 130.6679 130.2948 92.4562 61.8604 0.525 standard score, squared distance 6.0000 5.9790 4.9044 3.8910 0.349 rank percentage, absolute difference 1.3983 1.3907 1.1475 0.9809 0.295 pair disagreement between partitions (ratio of differing point pairs) raw scale, squared distance - standard score, squared distance 0.4192 raw scale, squared distance - rank percentage, absolute difference 0.4801 standard score, squared distance - rank percentage, absolute difference 0.3750 means of the largest cluster raw scale, squared distance n 509 mean_consumption=17.40 volatility=6.86 standard score, squared distance n 511 mean_consumption=19.18 volatility=11.97 rank percentage, absolute difference n 325 mean_consumption=21.50 volatility=12.66 ``` One rule applies when reading the first table: **no comparison across rows.** 61.8604 and 0.9809 are not two values of the same measure; one is on the order of cubic meters squared, the other a sum of percentage differences. The only comparable column is the share on the right: each definition's gain over its own random baseline, 0.525, 0.349 and 0.295 in order. The share is highest on the raw scale, but that does not make the definition better — there almost all of the total scatter comes from a single column, so splitting is the same as cutting that column. The share measures not the success of the definition but **how one-sidedly it looks at the data**. The second table is the real result of the lesson. The three partitions part on a share of point pairs between 0.3750 and 0.4801. The 0.4801 between the raw scale and the rank percentage means the two partitions give opposite answers on 137,002 of the 285,390 pairs. These are three runs on the same data, with the same procedure and number of clusters; only the definition of distance changes. The third table shows where the difference sits. On the raw scale the largest cluster has 509 subscribers and a mean volatility of 6.86 — mid-consumption, low-volatility subscribers together. With the standard score the cluster is of similar size but its mean volatility is 11.97: the volatility column now has as much say as consumption. With the rank percentage the largest cluster drops to 325, because extreme values, once turned into ranks, do not push points apart. No measure can say that one of these three partitions is the correct one. Each definition finds its own partition best on its own measure, and no higher court can be set up between definitions, because setting one up means choosing a distance definition once more. In the supervised course the label closed this question; here it stays open. #### Summary - When the label is removed, all three components of the measure must be built: the baseline, the method's number and the difference. In clustering the baseline is not given ready-made; it is chosen and justified. - Under the standard score definition a single cluster gives 6.0000 and random assignment 5.9790; a partition that does not learn moves the measure by only 0.0210 per point, and its fluctuation across five seeds is 0.0056. - The method's number is 3.8910; its difference is 2.0880 against the random baseline and 1.0134 against the slice baseline. Which baseline is written moves the reported gain nearly twofold. - When three similarity definitions split the same 756 subscribers with the same procedure, three partitions come out that disagree on a share of point pairs between 0.3750 and 0.4801; the raw scale and the rank percentage give opposite answers on 137,002 of 285,390 pairs. - The numbers of the definitions cannot be compared; the only thing that can be read is the proportional gain over their own baseline, and that share measures not the correctness of the definition but how one-sidedly it looks at the data. #### Next Step In this lesson the procedure was deliberately frozen: the starting points came from a single rule, two refinement rounds were run, the number of clusters was held at 4. All three are decisions. The next lesson resolves two of them: how many different partitions come out at the same number of clusters when the iteration continues until it stops and the starting points are chosen at random, and how large the gap is between the best and the worst run. After that the number of clusters is treated as a hyperparameter; because the within-cluster sum falls in only one direction as the number of clusters grows, it cannot point to the best value, and the "elbow" formed by the successive differences is counted as a decision point. ### k-Means Source: https://academia.sh/en/courses/unsupervised-learning/k-means Measuring iteration and the number of clusters separately: thirty random starts at the same k produce 26 distinct partitions, the numbers stay between 3.5649 and 3.7763 but the best and the worst run disagree on 0.2844 of the point pairs — even two runs that give the same number to four decimals disagree on 0.0037 of them. The number of clusters is a hyperparameter; in a sweep of 120 runs the random baseline falls from 6.0000 to 5.9181 while the method falls from 6.0000 to 2.2570, and three separate readings of the elbow choose three separate values of k: 3, 4 and 5. The previous lesson froze the procedure deliberately: the starting points were chosen with a single rule, two refinement rounds were run, the number of clusters was held at 4. The only thing that changed was the similarity definition. This lesson fixes the definition — standard score with squared distance — and frees two of the three frozen constants: the iteration continues until it stops, and the starting points are chosen at random. The resulting procedure is short: every point is assigned to the nearest centroid, every centroid moves to the mean of its own cluster, and this repeats until nothing changes. The lesson measures that this produces a decision in two separate places. The first is the start: how many distinct partitions come out when the same number of clusters is run from different starting points. The second is the number of clusters itself; k is a **hyperparameter** given from outside, and choosing it is a search. - **CL9.** The dataset, columns, split and measure are as in the previous lesson: the fictional subscriber table, seed 20260218, the 756 subscribers of the training set, six numeric columns, standard score, squared distance. The label is again not generated. - **CL10.** The measure is again the **within-cluster sum**, per point, smaller is better. - **CL11.** The baseline is **random assignment** and is computed separately for every k; because the baseline changes with the number of clusters, a single number cannot be written for the whole table. - **CL12.** The starting points are chosen from **distinct points** of the dataset; if a cluster ends up empty, its centroid is left where it was. - **CL13.** The iteration stops when the assignment does not change for two rounds in a row; the upper limit is 60 rounds and no run reaches it. A floating-point number is never compared with `==` — the stopping criterion is equality of the assignment lists, that is, an integer comparison. - **CL14.** The difference between two partitions is measured with the **differing point pair ratio** built in the previous lesson. - **CL15.** In the sweep over the number of clusters, 12 starts are tried for every k and the best is taken; **the number of candidates tried is written**. The reported number is also a function of the search budget. #### Where the Iteration Stops The procedure is the alternating repetition of two steps, and each step lowers the measure from a different direction. The assignment step holds the centroids fixed and gives every point to its nearest one. The update step holds the assignment fixed and moves every centroid to the mean of its own cluster; because the point that minimizes the sum of squared distances is the mean, this step too lowers the measure or leaves it the same. Since the measure never rises and is bounded below, the iteration must stop. About **where** it stops, this argument says nothing. ```python # k_means.py — MODEL. The FICTIONAL set and the standard score definition of lesson 01 # are rebuilt with the same seed; the label is still not generated. 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] N = len(V) 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)] def within_sum(assign, k): # within-cluster sum, per point t = 0.0 for c in range(k): g = [V[i] for i in range(N) if assign[i] == c] if g: t += sum(squared(v, centroid(g)) for v in g) return t / N def starting_points(r, k): # pick k distinct points as initial centroids s = [] while len(s) < k: j = int(r() * N) if j not in s: s.append(j) return [V[j] for j in s] def k_means(M, limit=60): # MODEL: assign, update centroids, repeat k, previous = len(M), None for it in range(1, limit + 1): a = [min(range(k), key=lambda c: squared(v, M[c])) for v in V] if a == previous: return a, it previous = a M = [centroid([V[i] for i in range(N) if a[i] == c]) or M[c] for c in range(k)] return a, limit r = rng(SEED + 5100) M0 = starting_points(r, 4) print("convergence trace from a single start (k = 4)") for it in range(1, 9): a = [min(range(4), key=lambda c: squared(v, M0[c])) for v in V] print(f" round {it} within-cluster sum {within_sum(a, 4):.4f}" f" cluster sizes {sorted(a.count(c) for c in range(4))}") M0 = [centroid([V[i] for i in range(N) if a[i] == c]) or M0[c] for c in range(4)] ``` ``` convergence trace from a single start (k = 4) round 1 within-cluster sum 4.2484 cluster sizes [109, 159, 233, 255] round 2 within-cluster sum 3.9366 cluster sizes [148, 166, 213, 229] round 3 within-cluster sum 3.7729 cluster sizes [153, 180, 198, 225] round 4 within-cluster sum 3.6519 cluster sizes [145, 175, 201, 235] round 5 within-cluster sum 3.6125 cluster sizes [140, 170, 215, 231] round 6 within-cluster sum 3.6058 cluster sizes [137, 172, 221, 226] round 7 within-cluster sum 3.6050 cluster sizes [138, 171, 222, 225] round 8 within-cluster sum 3.6049 cluster sizes [138, 170, 223, 225] ``` The number that was 4.2484 in the first round falls to 3.6049 by the eighth, and the fall slows quickly: from round one to round three 0.4755 is gained, from round five to round eight only 0.0076. The cluster sizes settle too, the 109–255 range of the first round narrowing to 138–225. The two-round procedure of the previous lesson had been cut off at 3.8910; here even the third round goes below it. What this trace does not say matters more. The number falling and stopping does not show that where it stopped is **the lowest number attainable**. Every round tries partitions neighboring the previous one; no round tries a placement far away. The procedure therefore stops at a solution near where it started. #### The Same k, a Different Start How decisive the start is gets measured by running thirty different sets of starting points at the same number of clusters. The identity of a partition, independent of cluster numbers, is which points stand together. ```python def identity(a, k): # the label-independent identity of a partition return frozenset(frozenset(i for i in range(N) if a[i] == c) for c in range(k)) def disagreement(a, b): # on how many point pairs the two partitions differ differing = sum((a[i] == a[j]) != (b[i] == b[j]) for i in range(N) for j in range(i + 1, N)) return 2 * differing / (N * (N - 1)) r, RUNS = rng(SEED + 7700), 30 RESULTS = [] for s in range(RUNS): a, it = k_means(starting_points(r, 4)) RESULTS.append((within_sum(a, 4), it, a)) RESULTS.sort(key=lambda x: x[0]) DISTINCT = {identity(a, 4) for _, _, a in RESULTS} print(f"{RUNS} random starts, k = 4 -> distinct partitions {len(DISTINCT)}") print(" sorted within-cluster sums") for b in range(0, RUNS, 10): print(" " + " ".join(f"{u:.4f}" for u, _, _ in RESULTS[b:b + 10])) print(f" best {RESULTS[0][0]:.4f}, median {RESULTS[RUNS // 2][0]:.4f}, " f"worst {RESULTS[-1][0]:.4f}, range {RESULTS[-1][0] - RESULTS[0][0]:.4f}") print(f" convergence rounds at least {min(t for _, t, _ in RESULTS)}, " f"at most {max(t for _, t, _ in RESULTS)}") print(f" lesson 01's single-pass procedure gave 3.8910; improvement " f"{3.8910 - RESULTS[0][0]:.4f}") print(f" best - second best: number difference {RESULTS[1][0] - RESULTS[0][0]:.6f}, " f"differing point pairs {disagreement(RESULTS[0][2], RESULTS[1][2]):.4f}") print(f" best - worst: number difference {RESULTS[-1][0] - RESULTS[0][0]:.4f}, " f"differing point pairs {disagreement(RESULTS[0][2], RESULTS[-1][2]):.4f}") ``` ``` 30 random starts, k = 4 -> distinct partitions 26 sorted within-cluster sums 3.5649 3.5649 3.5656 3.5656 3.5656 3.5658 3.5711 3.5711 3.5711 3.5711 3.5711 3.5712 3.5712 3.5727 3.5903 3.5906 3.5936 3.5936 3.5941 3.5955 3.5956 3.5957 3.5981 3.5995 3.6026 3.6077 3.6165 3.7467 3.7744 3.7763 best 3.5649, median 3.5906, worst 3.7763, range 0.2115 convergence rounds at least 11, at most 52 lesson 01's single-pass procedure gave 3.8910; improvement 0.3261 best - second best: number difference 0.000028, differing point pairs 0.0037 best - worst: number difference 0.2115, differing point pairs 0.2844 ``` Thirty runs gave **26 distinct partitions**. The same data, the same procedure, the same number of clusters; only the starting points changed. The numbers sit in a narrow band — between 3.5649 and 3.7763 — and twenty-seven of the thirty runs are below 3.62. Since the random assignment baseline is 5.9730 at this number of clusters, **every run beats the baseline by a wide margin**; even the worst is 2.1967 below it. The bottom two lines are the real finding of the lesson. The difference between the numbers of the best two runs is **0.000028**, the same number to four decimals; yet the two partitions disagree on **0.0037** of the point pairs. Between the best and the worst run the difference in numbers is 0.2115 and the disagreement is **0.2844** — more than a quarter of the pairs. The within-cluster sum reduces a partition to a single number, and what is preserved when that number stays the same is not who stands in the same cluster as whom. The convergence rounds vary between 11 and 52, and no run reaches the 60-round limit. The iteration really does stop, but where it stops depends on where the starting points fall. #### The Number of Clusters Is a Hyperparameter k is not something the model learns from the data but a number given from outside — a **hyperparameter** in the terminology of the Supervised Learning course. Choosing it by looking at the measure does not work directly: as the number of clusters grows, the within-cluster sum **falls in only one direction**, and the smallest value comes from the partition in which every point is its own cluster. The reading is therefore made not from the number itself but from the **successive differences**; the place where the differences shrink sharply is called the **elbow**. ```python def random_baseline(k, seed): r = rng(seed) return within_sum([int(r() * k) for _ in range(N)], k) STARTS, CANDIDATES, WITHIN = 12, 0, {} print(f"{'k':>3}{'baseline':>9}{'method':>9}{'difference':>11}{'successive difference':>23}{'drop ratio':>13}") for k in range(1, 11): r = rng(SEED + 8800 + k) WITHIN[k] = min(within_sum(k_means(starting_points(r, k))[0], k) for _ in range(STARTS)) CANDIDATES += STARTS baseline = statistics.fmean(random_baseline(k, SEED + 9900 + 50 * k + t) for t in range(5)) d = "" if k == 1 else f"{WITHIN[k - 1] - WITHIN[k]:23.4f}" o = "" if k == 1 else f"{(WITHIN[k - 1] - WITHIN[k]) / WITHIN[k - 1]:13.3f}" print(f"{k:>3}{baseline:>9.4f}{WITHIN[k]:>9.4f}{baseline - WITHIN[k]:>11.4f}{d}{o}") print(f"candidates tried {CANDIDATES} runs ({STARTS} starts x 10 cluster counts)") STEP = {k: WITHIN[k - 1] - WITHIN[k] for k in range(2, 11)} biggest_drop = max(range(2, 10), key=lambda k: STEP[k] - STEP[k + 1]) threshold_k = min(k for k in range(2, 11) if STEP[k] / WITHIN[k - 1] < 0.10) slope = (WITHIN[10] - WITHIN[1]) / 9 deviation = {k: (WITHIN[1] + slope * (k - 1)) - WITHIN[k] for k in range(2, 10)} chord_k = max(deviation, key=deviation.get) runner_up = sorted(deviation, key=deviation.get)[-2] print("\nthree separate readings of the elbow") print(f" k where the successive difference falls the most {biggest_drop}") print(f" k where the drop ratio first falls below 0.10 {threshold_k}") print(f" k farthest from the chord joining the endpoints {chord_k}" f" (deviation {deviation[chord_k]:.4f}; runner-up k = {runner_up}, deviation {deviation[runner_up]:.4f})") print(f" successive differences are not one-directional: k = 9 gives {STEP[9]:.4f}, k = 10 gives {STEP[10]:.4f}") ``` ``` k baseline method difference successive difference drop ratio 1 6.0000 6.0000 0.0000 2 5.9928 4.8521 1.1407 1.1479 0.191 3 5.9894 3.9996 1.9897 0.8524 0.176 4 5.9730 3.5655 2.4075 0.4342 0.109 5 5.9721 3.2151 2.7570 0.3504 0.098 6 5.9635 2.8969 3.0666 0.3182 0.099 7 5.9486 2.6732 3.2754 0.2237 0.077 8 5.9384 2.5020 3.4365 0.1712 0.064 9 5.9443 2.3874 3.5568 0.1146 0.046 10 5.9181 2.2570 3.6611 0.1304 0.055 candidates tried 120 runs (12 starts x 10 cluster counts) three separate readings of the elbow k where the successive difference falls the most 3 k where the drop ratio first falls below 0.10 5 k farthest from the chord joining the endpoints 4 (deviation 1.1869; runner-up k = 3, deviation 1.1686) successive differences are not one-directional: k = 9 gives 0.1146, k = 10 gives 0.1304 ``` The baseline column alone is a warning. The number for random assignment falls too as k grows: from 6.0000 to 5.9181. The fall is small but real, because every new cluster estimates its centroid from its own points and the estimation noise pulls the measure down. If even the number of a procedure that does not learn improves as k grows, the whole of the method's fall cannot be attributed to structure. That is why the difference column is needed. The elbow, on the other hand, does not give a single answer from the table. Three common readings were tried on the same ten rows and **the three chose three separate numbers of clusters**: the place where the successive difference falls the most is k = 3, the place where the drop ratio first falls below 0.10 is k = 5, and the point farthest from the chord joining the endpoints is k = 4. In the third, the winner's deviation is 1.1869 and the runner-up's 1.1686 — the gap is under 2%. The last line shuts the table altogether: the successive differences are not even one-directional, rising from 0.1146 at k = 9 to 0.1304 at k = 10. The number of candidates tried is **120 runs**, and this number stands beside every reported cell. The 3.5655 in the k = 4 row is the best of twelve starts; the search with thirty starts had found 3.5649 for the same k. The 0.0006 between them comes not from the data but from the **search budget**. #### Summary - The procedure is the alternating repetition of two steps and each step lowers the measure or leaves it the same; the iteration must therefore stop, but there is no guarantee that where it stops is the best. - Thirty random starts gave **26 distinct partitions** at k = 4; the numbers varied between 3.5649 and 3.7763, and convergence between 11 and 52 rounds. - A difference in numbers and a difference in partitions are not the same thing: two runs that give the same number to four decimals disagree on 0.0037 of the point pairs, the best and the worst run on 0.2844. - The random assignment baseline also falls with the number of clusters, from 6.0000 to 5.9181; once this baseline is subtracted, the method's difference rises from 2.4075 to 3.6611. - Three separate readings of the elbow chose three separate numbers of clusters — 3, 5 and 4 — and the successive differences are not even one-directional; the chosen number is the result of a search budget of 120 runs. #### Next Step The two constraints of this procedure are tied together: the number of clusters must be given **before starting**, and every run produces a single partition. Three separate searches were made for k = 3, k = 4 and k = 5 and the three results came out unrelated — the clusters of k = 3 are not unions of the clusters of k = 5. The next lesson takes a family that reverses the order: every point starts as its own cluster, the nearest two clusters are merged, and a **nested family of clusters** is produced up to a single cluster. The number of clusters is here not an input but the answer to the question of which height the tree is cut at. In return a new decision appears: "the nearest two clusters" cannot be derived from point distance, and the distance between clusters must be defined separately. The lesson counts on how many points three **linkage criteria** disagree on the same data and prints the dendrogram as text. ### Hierarchical Clustering Source: https://academia.sh/en/courses/unsupervised-learning/hierarchical-clustering A family that turns the number of clusters from an input into a cut decision, and the price of that decision: on the same 200 subscribers, single, complete and average linkage start from the same first merge (0.079) and end at 2.245, 8.917 and 4.994; cut at a height of 2.0 they give 3, 56 and 37 clusters respectively. At four clusters, single linkage beats the random assignment baseline by only 0.1492 with 6.0849 and produces [1, 1, 1, 197], while complete linkage beats it by 1.6922 with 4.5419; the two linkages place 86 of the 200 subscribers in different clusters. 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. ```python # 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. ```python 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. ```python 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. ### Probabilistic Clustering Source: https://academia.sh/en/courses/unsupervised-learning/probabilistic-clustering Measuring a mixture model in which a point is attached to every cluster with a responsibility: the fit score of the random assignment baseline is -4.2424, that of the model found by iterative estimation is -3.7755, a difference of +0.4669. The same difference comes out as +0.3358 on the dataset whose columns are shuffled within themselves, so only +0.1311 of the difference comes from the data. For 74 of the 756 points the largest share is below 0.60, and for these points the mean largest share is 0.5418. On moving to hard assignment the mean uncertainty falls from 0.3797 to zero and the cluster sizes deviate by up to 28.8 points from the soft mass. At a second start, 43.24% of the undecided points and 7.20% of the certain points change cluster. The procedure of the previous lesson wrote a point into a cluster, placed it on the border, or left it outside. All three cases were definite, and only the result stood in the ledger of decisions. A subscriber who falls right between two dense regions and a subscriber who stands at the heart of a cluster were recorded in the same way; how hard the decision was to make was written nowhere. The same held for centroid-based and hierarchical partitioning: the output was a cluster number, not a degree of confidence. This lesson turns the assignment into a number. Every point is attached to every cluster with a **responsibility**, the shares of a point sum to one, and the closeness of the shares directly measures the undecidedness of that point. For this, the clusters are no longer written as lists of centroids but as **components** with their own center and their own variability; the dataset as a whole is a **mixture model** of these components. There are two things to measure: how much the model beats a random assignment, and how much of the information carried by the soft assignment is lost when one steps down to hard assignment. - **CL37.** 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. The label column stays aside in this lesson as well. - **CL38.** Three continuous columns are used: mean consumption, volatility and consumption per person. Columns with integer values are left out, because a component can collapse onto a single integer value and reduce the variability in that column to zero. - **CL39.** Every component is defined by a weight, one center per column and one variability per column. Dependence between columns is not modeled; this is a **model decision** and it requires the components to be axis-aligned. - **CL40.** A **lower bound of 0.02** is placed on the variability. Without a lower bound a component can collapse onto a single point, the density at that point grows without limit, and the fit score loses its meaning. The bound is one fiftieth of unit variability on the scaled columns. - **CL41.** The **fit score** is the per-point mean of the natural logarithm of the density that the model gives to a point. Higher is better, and the score is comparable only on the same set of columns. - **CL42.** The baseline is **random assignment**: the points are distributed randomly into three groups and the parameters are computed from that distribution. This baseline is chosen because it uses **the same model form** and the same number of parameters; the difference between them comes only from the assignment, not from the model family. - **CL43.** The number of clusters is three and is a **hyperparameter**; the justification for the choice and how arbitrary the choice is are measured in the closing lesson of the topic. - **CL44.** The structureless reference is built, as in the previous lesson, by **within-column shuffling**; the tool belongs to the Statistical Analysis and Exploratory Data Analysis course and the testing theory is not repeated. #### A Mixture of Components The model is this: every point in the dataset is taken to have come from one of three components, but which one it came from is not recorded. Every component $c$ carries a weight $a_c$ (the share of the points that come from that component) and keeps, for each column, a center $m_{cd}$ and a variability $v_{cd}$. The density that the component gives to a point fades quickly with distance from the center: $$f_c(x) = \prod_{d} \frac{1}{\sqrt{2\pi v_{cd}}}\, \exp\left(-\frac{(x_d - m_{cd})^2}{2\,v_{cd}}\right)$$ The share by which a point is attached to a component is the ratio of the contribution of that component to the total contribution: $$r_{ic} = \frac{a_c\, f_c(x_i)}{\sum_j a_j\, f_j(x_i)}$$ Since the parameters are unknown the shares cannot be computed, and since the shares are unknown the parameters cannot be computed. **Iterative estimation** breaks this loop: it temporarily fixes one side, computes the other, then swaps. Every round raises the fit score or leaves it where it is; where it arrives depends on the start. ```python # mixture.py — MODEL. The FICTIONAL subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering and Supervised Learning courses is rebuilt with the same seed; # the label column stays aside in this course as well. 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 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[1], household, rating)) for sub_id, base, 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( base * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100) r() if v: mean_value = round(sum(v) / len(v), 2) DATA.append({"mean_consumption": mean_value, "volatility": round(max(v) - min(v), 2), "per_person": round(mean_value / household, 3)}) CHOSEN = ["mean_consumption", "volatility", "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 CHOSEN} P = [[(x[a] - SCALE[a][0]) / SCALE[a][1] for a in CHOSEN] for x in TRAIN] K, FLOOR = 3, 0.02 # number of clusters and variability floor def parameters(P, member, k): # MODEL: weight, center, variability A, M, V = [], [], [] for c in range(k): w = [u[c] for u in member] t = sum(w) + 1e-9 m = [sum(w[i] * P[i][d] for i in range(len(P))) / t for d in range(len(P[0]))] V.append([max(sum(w[i] * (P[i][d] - m[d]) ** 2 for i in range(len(P))) / t, FLOOR) for d in range(len(P[0]))]) A.append(t / len(P)) M.append(m) return A, M, V def responsibility(P, A, M, V): # soft assignment and fit score R, total = [], 0.0 for x in P: g = [math.log(A[c]) + sum(-0.5 * (math.log(2 * math.pi * V[c][d]) + (x[d] - M[c][d]) ** 2 / V[c][d]) for d in range(len(x))) for c in range(len(A))] b = max(g) s = sum(math.exp(u - b) for u in g) total += b + math.log(s) R.append([math.exp(u - b) / s for u in g]) return R, total / len(P) def random_membership(n, k, seed): # BASELINE: random assignment r, U = rng(seed), [] for _ in range(n): c = int(r() * k) U.append([1.0 if j == c else 0.0 for j in range(k)]) return U def estimate(P, k, seed, rounds=80): # MODEL: iterative estimation U = random_membership(len(P), k, seed) for _ in range(rounds): U, score = responsibility(P, *parameters(P, U, k)) return U, score def column_shuffle(P, seed): Y = [list(x) for x in P] for s in range(len(P[0])): for i, j in enumerate(shuffle(list(range(len(P))), seed + s * 17)): Y[i][s] = P[j][s] return Y def baseline_score(P, k, seed): return responsibility(P, *parameters(P, random_membership(len(P), k, seed), k))[1] U, SCORE = estimate(P, K, SEED + 120000) BASELINE = baseline_score(P, K, SEED + 120000) SHUFFLED = column_shuffle(P, SEED + 610000) # FICTIONAL: within-column reshuffling _, SHUF_SCORE = estimate(SHUFFLED, K, SEED + 120000) SHUF_BASELINE = baseline_score(SHUFFLED, K, SEED + 120000) print(f"points {len(P)}, columns {len(CHOSEN)}, components {K}") print(f"{'dataset':<16}{'baseline':>15}{'mixture model':>16}{'difference':>12}") print(f"{'data':<16}{BASELINE:>15.4f}{SCORE:>16.4f}{SCORE - BASELINE:>+12.4f}") print(f"{'shuffled':<16}{SHUF_BASELINE:>15.4f}{SHUF_SCORE:>16.4f}{SHUF_SCORE - SHUF_BASELINE:>+12.4f}") print(f"share from the data {(SCORE - BASELINE) - (SHUF_SCORE - SHUF_BASELINE):+.4f}") ``` ``` points 756, columns 3, components 3 dataset baseline mixture model difference data -4.2424 -3.7755 +0.4669 shuffled -4.2443 -3.9085 +0.3358 share from the data +0.1311 ``` The random assignment baseline is **-4.2424**, the mixture model **-3.7755**, the difference **+0.4669**. The baseline has been chosen deliberately hard here: random assignment also builds three components, computes the same number of parameters and uses the same density form. Its only lack is that it **does not learn** which point belongs to which component. In this way the difference comes not from the model family but only from the assignment. The second row says how much of this difference comes from the data. On the dataset whose columns are shuffled within themselves the same procedure gains **+0.3358**. That is, about three quarters of the gain stands even when no relationship between columns is left. The reason is clear: three components wrap a skewed, long-tailed distribution axis by axis better than a single component could; this says nothing about the cluster structure in the data, it says something about the flexibility of the model. The share that falls to the structure between columns is **+0.1311**, about 28% of the total difference. The question of the course finds its answer exactly here: the answer to whether the difference comes from the measure or from the data is often "from both", and without separating the shares the number cannot be read. #### The Uncertainty That Soft Assignment Carries The responsibility keeps the information that hard assignment has to throw away. If the shares of a point are $0.98 / 0.01 / 0.01$ the assignment is definite; if they are $0.40 / 0.38 / 0.22$ the assignment is a choice, and how fine the choice is stands written in the shares. To reduce the uncertainty to a single number, the **sum of the shares weighted by their natural logarithm** is used; this number is zero when all the weight is given to a single share and reaches its largest value when the shares are distributed equally. ```python def uncertainty(r): return -sum(p * math.log(p) for p in r if p > 1e-12) HARD = [max(range(K), key=lambda c: r[c]) for r in U] uncertain = [i for i, r in enumerate(U) if max(r) < 0.60] A, _, _ = parameters(P, U, K) print(f"points with largest responsibility below 0.60: {len(uncertain)}; " f"for these points the mean largest share is {sum(max(U[i]) for i in uncertain) / len(uncertain):.4f}") print(f"points with the top two shares closer than 0.10: {sum(1 for r in U if sorted(r)[-1] - sorted(r)[-2] < 0.10)}" f"; largest share 0.99 and above: {sum(1 for r in U if max(r) >= 0.99)}") print(f"mean uncertainty: soft assignment " f"{sum(uncertainty(r) for r in U) / len(U):.4f}, hard assignment 0.0000, " f"upper bound {math.log(K):.4f}") print(f"{'cluster':<8}{'soft mass':>15}{'hard count':>12}{'difference':>12}") for c in range(K): print(f"{c:<8}{A[c] * len(P):>15.1f}{HARD.count(c):>12}" f"{HARD.count(c) - A[c] * len(P):>+12.1f}") ``` ``` points with largest responsibility below 0.60: 74; for these points the mean largest share is 0.5418 points with the top two shares closer than 0.10: 28; largest share 0.99 and above: 94 mean uncertainty: soft assignment 0.3797, hard assignment 0.0000, upper bound 1.0986 cluster soft mass hard count difference 0 217.8 189 -28.8 1 200.3 213 +12.7 2 337.9 354 +16.1 ``` The threshold of 0.60 is chosen because an equal distribution over three components comes to 0.3333, and a largest share below 0.60 means that **the majority of the weight of the point is not gathered in a single cluster**. **74** of the 756 points are in this situation, and for these points the largest share is **0.5418** on average. For twenty-eight points the top two shares differ from each other by less than 0.10 — for these points the cluster chosen by hard assignment and the cluster in second place are almost equally justified. At the other end, 94 points are attached to a single cluster with a share of 0.99 and above. The mean uncertainty is **0.3797** in soft assignment, against an upper bound of **1.0986**. On moving to hard assignment this number falls to **zero** by definition: every point is written to a single cluster with full weight, and the information that 74 points are undecided does not remain in the table. The concrete counterpart of the loss is in the cluster sizes. In soft assignment the mass of the first cluster comes to 217.8 points, while hard assignment writes 189 points to that cluster; the deviation of **28.8 points** comes from all the points on the border being rounded to one side. Every decision that rests on the relative size of the clusters inherits this rounding. #### Testing the Shares Up to here it has been assumed that the responsibility is meaningful. This is a testable claim: if the share really measures how hard the decision is, undecided points should change cluster when the conditions change a little, and certain points should not. The condition to change is ready — iterative estimation is sensitive to the start, and the second lesson of the topic measured this for centroid-based partitioning. The model run from a second start is compared with the first after the clusters are matched through their centers. ```python U2, SCORE2 = estimate(P, K, SEED + 555000) _, M, _ = parameters(P, U, K) _, M2, _ = parameters(P, U2, K) MATCH = [min(range(K), key=lambda c: math.dist(M2[j], M[c])) for j in range(K)] HARD2 = [MATCH[max(range(K), key=lambda c: r[c])] for r in U2] changed = [i for i in range(len(P)) if HARD[i] != HARD2[i]] certain = [i for i, r in enumerate(U) if max(r) >= 0.90] changed_uncertain = sum(1 for i in changed if max(U[i]) < 0.60) changed_certain = sum(1 for i in changed if max(U[i]) >= 0.90) print(f"fit score of the second start {SCORE2:.4f} (first {SCORE:.4f}); " f"cluster matching {MATCH}") print(f"points that changed cluster {len(changed)} / {len(P)}") print(f"change among uncertain points {changed_uncertain}/{len(uncertain)} = {changed_uncertain / len(uncertain):.4f}") print(f"change among certain points (share >= 0.90) {changed_certain}/{len(certain)} = {changed_certain / len(certain):.4f}") ``` ``` fit score of the second start -3.7827 (first -3.7755); cluster matching [1, 0, 2] points that changed cluster 125 / 756 change among uncertain points 32/74 = 0.4324 change among certain points (share >= 0.90) 29/403 = 0.0720 ``` The two runs stand in different places: the fit scores are **-3.7755** and **-3.7827**. The difference between them is small, but the difference between the partitions is not small — **125** of the 756 points change cluster. This shows that looking only at the fit score and saying "the model converged" is not enough. The real number is in the disagreement. **32** of the 74 undecided points, that is **43.24%**, change cluster. Of the 403 points with a share of 0.90 and above, only **29**, that is **7.20%**. The ratio is more than sixfold. The responsibility correctly estimates something outside the model it was computed from: how fragile the assignment of a point is. Hard assignment does not carry this information — a cluster number shows a number that will change with a 43% probability and a number that will change with a 7% probability in the same way. #### Summary - The mixture model writes the clusters as components defined by a center and a variability; every point is attached to every component with a responsibility, and the shares sum to one. - The random assignment baseline is -4.2424, the fit score of the model is -3.7755, the difference is +0.4669; because the baseline uses the same model form and the same number of parameters, the difference comes only from the assignment. - The same difference comes out as +0.3358 on the dataset shuffled within its columns; the share that falls to the structure between columns is +0.1311, that is, about 28% of the total difference. The rest is the flexibility of the model. - For 74 of the 756 points the largest share is below 0.60 (mean 0.5418), and for 28 points the top two shares are closer than 0.10. On moving to hard assignment the mean uncertainty falls from 0.3797 to zero, and the count of one cluster deviates by 28.8 points from its soft mass. - The responsibility also measures something outside the model: at a second start, 43.24% of the undecided points and 7.20% of the certain points change cluster. #### Next Step In this lesson the number of clusters was taken as three and this choice was not justified. Similar choices were made in the previous three lessons as well: how many clusters, which cut height, which neighborhood radius. Each time it was shown that the choice moved the result, but how to make the choice **itself** was postponed. The closing lesson of the topic takes up this question and shows with numbers that there is no answer: on the same data three different internal validity measures choose three different numbers of clusters, none of them more correct than another, and on a structureless dataset generated entirely at random the same three measures still choose a number of clusters. After that, the label is opened, for the first and only time in this course, for **external validity**, and it is counted how far the number of clusters that the label points to disagrees with those chosen by the internal measures. ### Clustering Evaluation Source: https://academia.sh/en/courses/unsupervised-learning/clustering-evaluation Measuring the measures that choose the number of clusters: on the same 756 subscribers the silhouette coefficient chooses 3 clusters, the scatter ratio 6 and the gap measure 10, and none of them is more correct than another. On a random cloud that carries no structure the same three measures still choose 9, 7 and 9; the highest silhouette is 0.2142 on the real data and 0.1255 on the structureless cloud. The label is opened for the first and only time in the course, and the pair agreement rate comes out as 0.5235 for the best partition: the single cluster baseline, at 0.6539, beats all nine partitions. Twenty-seven comparisons were made, and none of the internal measures chose the number of clusters that the label points to. 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. ```python # 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. ```python 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. ```python 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. ## Dimensionality Reduction and Anomalies Measuring what is kept and what is lost when the number of columns is reduced, and how a point that stands alone is flagged: the reconstruction error of the projection that keeps the variability, how much of the pairwise distances survive, and the scaling decision changing the projection; the trade-off between the error that grows as the rank falls and the storage cost; how many of the nearest neighbors remain neighbors in a neighborhood-preserving embedding and the distance between two points being uninterpretable; how many of the flagged anomalies are real and the multivariate anomaly the univariate procedure cannot catch; high-confidence rules coming only from a frequent item and lift correcting this. ### Principal Component Analysis Source: https://academia.sh/en/courses/unsupervised-learning/principal-component-analysis Measuring what is kept and what is lost when a twelve-column subscriber table is reduced to fewer numbers: reconstruction with column means gives 3.3411 while reconstruction with two components gives 2.7762, and the difference is only 0.169 of the baseline; the same two components carry 0.341 of the component share, yet the mean of the pairwise distances falls to 0.4965 times the original and only 0.0057 of the 31,626 pairs survive within a 10% margin. In the unscaled set the load of the first component binds to a single column with +0.96 and the share of the first two components rises to 0.863, but the reconstruction error measured in a common unit, 2.8473, stays behind the scaled projection. Six lessons grouped points. The cluster count was chosen, the linkage criterion was changed, density thresholds were swept, and the results were read against a random partition. Underneath all of these operations there was a single operation: computing the distance between two subscribers. The distance was found each time by looking at all the columns at once, and how many columns there were was never questioned. The question of this topic is this: how many numbers does it really take to describe a point. In a table of twelve columns each subscriber stands on twelve numbers, but the columns are not independent — consumption per person is derived from mean consumption and household size, and the region indicators exclude each other. Instead of dropping columns, mixing them and reducing them to fewer numbers is an option, and what this lesson measures is the **cost** of that reduction: when a point is reconstructed from a few numbers, how far from its original place does it land. - **DR1.** The fictional subscriber table is the same one used in the Introduction to Machine Learning, Data Preparation and Feature Engineering, and Supervised Learning courses, and it **is fictional**: seed 20260218, 1,260 rows, split 756/252/252. The label column is **never generated** in this course. - **DR2.** Twelve columns are used: six numeric measurements, the zero reading indicator and five region indicators. - **DR3.** The scale and the component directions are learned **from the training set only** and applied to the test set with the same numbers. If they were learned from the whole set, the measurement would be inflated. - **DR4.** The **component share** computation is not repeated here. The Statistical Analysis and Exploratory Data Analysis course had built the share as an exploration tool, and the feature selection topic of the Data Preparation and Feature Engineering course had built it as a feature measure; in this lesson the share is only an **input** and is used to be read. - **DR5.** The **reconstruction error** is the per-point mean of the distance between a point's original place in the scaled space and its place reconstructed from a few numbers. Its unit is the standard deviation. - **DR6.** The baseline is **reconstruction with column means**: no component is used, and every point is set equal to the mean point of the training set. This is the plainest reconstruction that does not learn. - **DR7.** A second baseline is read as well: **k random orthogonal directions**. Only this separates whether the difference of the components comes from the number of directions or from their choice. Three draws are averaged. - **DR8.** The directions are found by iteration within the lesson and are the **model**: the covariance matrix is built, the largest direction is searched for in 200 iterations, and each new direction is made orthogonal to the previous ones. No library is called. #### Reconstruction Error and Its Baseline The procedure has three steps. The columns are brought to a common scale, the covariance matrix of the training set is built, and the direction that carries the largest variability of this matrix is found by iteration and appended to the list. Each new direction is forced to be orthogonal to those before it. When a point is projected onto k directions and reconstructed, the remaining distance is everything those k directions do not carry. ```python # COMMON — MODEL: the FICTIONAL subscriber table of the Introduction to Machine Learning course is rebuilt with the same seed. # The label column is never generated in this course. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 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 SUBSCRIBERS, DATA = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] 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({"id": 10001 + i, "region": b, "household": household, "satisfaction": rating}) for k in SUBSCRIBERS: r, v = rng(SEED + 7000 + k["id"]), [] for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): v.append(0.0 if r() < 0.038 else math.floor( k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - n * 0.05) * 100 + 0.5) / 100) r() if not v: continue mean = sum(v) / len(v) x = {"id": k["id"], "household": k["household"], "satisfaction": k["satisfaction"], "period": len(v), "mean_consumption": round(mean, 2), "volatility": round(max(v) - min(v), 2), "zero_readings": int(min(v) == 0.0), "per_person": round(mean / k["household"], 3)} for b in REGION: x["area_" + b[0]] = int(k["region"][0] == b[0]) DATA.append(x) FIELDS = ["mean_consumption", "volatility", "per_person", "household", "satisfaction", "period", "zero_readings", "area_north", "area_south", "area_east", "area_west", "area_central"] 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] K = shuffle(DATA, SEED + 90000) TRAIN, VAL, TEST = K[:756], K[756:1008], K[1008:] ``` ```python n = len(FIELDS) MEAN = {a: sum(x[a] for x in TRAIN) / len(TRAIN) for a in FIELDS} STD_DEV = {a: max(1e-9, (sum((x[a] - MEAN[a]) ** 2 for x in TRAIN) / len(TRAIN)) ** 0.5) for a in FIELDS} def multiply(m, v): return [sum(m[i][j] * v[j] for j in range(len(v))) for i in range(len(m))] def orthogonalize(v, B): for u in B: c = sum(v[i] * u[i] for i in range(len(v))) v = [v[i] - c * u[i] for i in range(len(v))] return v def find_directions(Z): C = [[sum(z[i] * z[j] for z in Z) / len(Z) for j in range(n)] for i in range(n)] B, share = [], [] for _ in range(n): v, lam = [1.0 / (i + 1) for i in range(n)], 0.0 for _ in range(200): w = orthogonalize(multiply(C, v), B) u = sum(q * q for q in w) ** 0.5 if u < 1e-9: break v, lam = [q / u for q in w], u if lam < 1e-9: break B.append(v) share.append(lam) return B, share def reconstruct(z, B): c = [sum(z[j] * b[j] for j in range(n)) for b in B] return [sum(c[i] * B[i][j] for i in range(len(B))) for j in range(n)] def error(Z, B): return sum(sum((z[j] - reconstruct(z, B)[j]) ** 2 for j in range(n)) ** 0.5 for z in Z) / len(Z) def random_directions(k, seed): r, B = rng(seed), [] while len(B) < k: v = orthogonalize([r() + r() + r() - 1.5 for _ in range(n)], B) u = sum(q * q for q in v) ** 0.5 if u > 1e-6: B.append([q / u for q in v]) return B Z_TRAIN = [[(x[a] - MEAN[a]) / STD_DEV[a] for a in FIELDS] for x in TRAIN] Z_TEST = [[(x[a] - MEAN[a]) / STD_DEV[a] for a in FIELDS] for x in TEST] COMPONENTS, SHARE = find_directions(Z_TRAIN) TOTAL, BASELINE = sum(SHARE), error(Z_TEST, []) print(f"columns {n}, directions found {len(COMPONENTS)}, baseline {BASELINE:.4f}") print(f"{'k':>3} {'share':>7} {'train':>8} {'held-out':>9} {'random':>9} {'difference':>10}") for k in (1, 2, 3, 4, 6, 8, 10, 11): B, s = COMPONENTS[:k], 0.0 for t in range(3): s += error(Z_TEST, random_directions(k, SEED + 900 * (t + 1))) print(f"{k:>3} {sum(SHARE[:k]) / TOTAL:>7.3f} {error(Z_TRAIN, B):>8.4f} " f"{error(Z_TEST, B):>9.4f} {s / 3:>9.4f} {error(Z_TEST, B) - BASELINE:>+10.4f}") ``` ``` columns 12, directions found 11, baseline 3.3411 k share train held-out random difference 1 0.190 3.0359 3.0295 3.1971 -0.3116 2 0.341 2.7663 2.7762 3.0531 -0.5650 3 0.458 2.5071 2.5419 2.8281 -0.7993 4 0.571 2.1700 2.2211 2.6875 -1.1200 6 0.771 1.5551 1.5474 2.2723 -1.7937 8 0.935 0.7878 0.8316 1.8246 -2.5096 10 0.989 0.2625 0.2668 1.1637 -3.0743 11 1.000 0.0000 0.0000 0.8436 -3.3411 ``` Three numbers are in the table. The baseline is 3.3411: when no component is used, a point is this far from the mean point. Reconstruction with two components gives 2.7762; the difference is 0.5650, that is, 0.169 of the baseline. At the same k the mean of three random orthogonal directions is 3.0531 — the difference of the components from random directions is 0.2769. The rest of the drop, 0.2880, comes solely from **using two directions**; even if it were not known which two directions they were, the baseline falls that much. The contribution of the components is only the part that remains above this portion, and it is close to half of the difference from the baseline. The second reading of the table is that the share column and the error column **do not go at the same speed**. At two components the share is 0.341 — it sounds like more than a third — but the error has dropped by only 0.169. The share is read over squared quantities, the reconstruction error over distance; half of the share is not half of the distance. Reading these two columns in place of each other is the most frequent mistake in dimensionality reduction. The third reading is in the last row. There are twelve columns but the procedure could find only **eleven directions**, and at the eleventh direction the error fell to exactly zero. Because the five region indicators sum to one in every row, one of them can be read from the others; there is no such thing as a twelfth direction. The table holds twelve numbers, but it carries no more than eleven. #### How Much of the Distances Survive The reconstruction error measures the **place** of a point. Clustering, however, uses not the place of a point but the **distance between points**. A projection does not have to preserve these two to the same extent, because a projection can only shorten distances — it can never lengthen them. ```python def distance_preservation(Z, B): P = [[sum(z[j] * b[j] for j in range(n)) for b in B] for z in Z] total, count, preserved = 0.0, 0, 0 for a in range(len(Z)): for b in range(a + 1, len(Z)): d = sum((Z[a][j] - Z[b][j]) ** 2 for j in range(n)) ** 0.5 if d < 1e-9: continue dp = sum((P[a][i] - P[b][i]) ** 2 for i in range(len(B))) ** 0.5 total, count = total + dp / d, count + 1 preserved += dp / d >= 0.9 return total / count, preserved / count, count print(f"{'k':>3} {'share':>7} {'mean ratio':>10} {'preserved':>9} {'random ratio':>14}") for k in (1, 2, 4, 6, 8, 10): o, g, c = distance_preservation(Z_TEST, COMPONENTS[:k]) r = distance_preservation(Z_TEST, random_directions(k, SEED + 900))[0] print(f"{k:>3} {sum(SHARE[:k]) / TOTAL:>7.3f} {o:>10.4f} {g:>9.4f} {r:>14.4f}") print(f"pairs measured {c}") ``` ``` k share mean ratio preserved random ratio 1 0.190 0.3398 0.0004 0.2341 2 0.341 0.4965 0.0057 0.3747 4 0.571 0.7029 0.0542 0.5358 6 0.771 0.8535 0.4299 0.6945 8 0.935 0.9540 0.8910 0.7924 10 0.989 0.9943 0.9994 0.9049 pairs measured 31626 ``` The second table is harsher than the first. At two components the pairwise distances fall on average to 0.4965 times the original — half of the distances are gone — and only 0.0057 of the 31,626 pairs survive within a 10% margin. At six components the mean ratio rises to 0.8535 and 0.4299 of the pairs are preserved. The first row where most of the distances survive is eight components: 0.9540 and 0.8910. The conclusion drawn from this is clear. Clustering on a two-component projection means running a distance-based procedure after throwing away half of the distances. Random directions giving 0.3747 at the same k shows that the components gain something, but the gain lies between 0.4965 and 0.3747, and both numbers are far from full preservation. The component share does not foretell any of these rows: when the share is 0.341, the preserved-pair rate is 0.0057, and the two numbers are not even of the same order of magnitude. #### The Scaling Decision Changes the Projection Components look for the directions that carry the most variability. Because variability is read from the unit of the column, when the columns are in different units the column with the largest unit pulls the largest direction to itself. This is not a flaw but the definition of the procedure — but if the scaling decision is made without measuring where it turns the projection, it is not a decision but an accident. ```python for scaled in (True, False): divisor = (lambda a: STD_DEV[a]) if scaled else (lambda a: 1.0) H_TRAIN = [[(x[a] - MEAN[a]) / divisor(a) for a in FIELDS] for x in TRAIN] H_TEST = [[(x[a] - MEAN[a]) / divisor(a) for a in FIELDS] for x in TEST] B, P = find_directions(H_TRAIN) order = sorted((-abs(B[0][j]), FIELDS[j], B[0][j]) for j in range(n)) print(f"{'scaled' if scaled else 'raw':<8} first two share " f"{(P[0] + P[1]) / sum(P):.3f} c0 loadings: " + ", ".join(f"{a} {w:+.2f}" for _, a, w in order[:4])) remainder, total = [0.0] * n, 0.0 for z in H_TEST: g = reconstruct(z, B[:2]) d = [(z[j] - g[j]) / (1.0 if scaled else STD_DEV[FIELDS[j]]) for j in range(n)] total += sum(q * q for q in d) ** 0.5 for j in range(n): remainder[j] += d[j] * d[j] / len(H_TEST) print(f"{'':8} k=2 error in common unit {total / len(H_TEST):.4f} column remainders: " + ", ".join(f"{FIELDS[j]} {remainder[j]:.2f}" for j in (0, 1, 4, 11))) ``` ``` scaled first two share 0.341 c0 loadings: per_person +0.55, mean_consumption +0.54, area_west -0.34, household -0.31 k=2 error in common unit 2.7762 column remainders: mean_consumption 0.36, volatility 0.19, satisfaction 1.00, area_central 1.07 raw first two share 0.863 c0 loadings: volatility +0.96, mean_consumption +0.24, per_person +0.12, period +0.04 k=2 error in common unit 2.8473 column remainders: mean_consumption 0.10, volatility 0.00, satisfaction 1.00, area_central 1.21 ``` The two rows take the same data to two different places. In the unscaled set the first component is almost entirely a single column: the `volatility` loading is +0.96 and none of the remaining eleven columns exceeds 0.25 in that component. The share of the first two components rises to 0.863. In the scaled set the same component carries `per_person` +0.55 and `mean_consumption` +0.54, `volatility` does not even make the top four, and the share of the first two components stays at 0.341. Someone looking only at the share column would choose the unscaled projection: 0.863 is two and a half times 0.341. The reconstruction error measured in the common unit says the opposite — 2.8473 unscaled, 2.7762 scaled. The unscaled projection collects 0.863 of the share because the raw unit of a single column swallows the whole share; it reconstructs that column perfectly (`volatility` remainder 0.00) and in return leaves the `satisfaction` column with a remainder of 1.00, that is, in a place no better than the baseline. The remainder in the `area_central` column is 1.21 — the projection reconstructs that column **worse** than the column mean, because the direction it picks pushes that column the wrong way. The scaled projection is not perfect either: `satisfaction` stands at a remainder of 1.00 there too, and `area_central` at 1.07. The difference lies in where the loss is distributed. The scaling decision is the choice of which column will be saved, and this choice was made in the scaling step, long before the components were computed. #### Summary - The baseline of the reconstruction error is reconstruction with column means and it is 3.3411; reconstruction with two components gives 2.7762, and the difference is 0.169 of the baseline. - The second baseline is random orthogonal directions: 3.0531 at the same k. The contribution of the components is 0.2769, and the remaining 0.2880 of the drop comes solely from using two directions. - The component share and the reconstruction error do not go at the same speed: with a share of 0.341 the error drops by only 0.169, because the share is read over squared quantities and the error over distance. - The preserved geometry is less than the share: at two components the pairwise distances fall to 0.4965 times the original, and only 0.0057 of the 31,626 pairs survive within a 10% margin. - Scaling changes the projection: in the unscaled set the share of the first two components rises to 0.863, but the error in the common unit, 2.8473, stays behind the scaled one, and some columns are reconstructed worse than the baseline. #### Next Step In this lesson the components were found iteratively through the covariance matrix, and every point was projected onto the directions found and reconstructed. The procedure worked, but it never asked one question: how many numbers are needed to store the reconstructed table. A 252-row table of twelve columns is 3,024 numbers; with how many numbers does the two-component representation hold the same table, and by what curve does the error rise as that count falls. The next lesson writes the projection as a matrix decomposition, brings the notion of rank into the measure, and prints the trade-off between the storage cost and the reconstruction error together with the same baseline. ### Singular Value Decomposition Application Source: https://academia.sh/en/courses/unsupervised-learning/singular-value-decomposition-application Measuring the cost of storing a matrix with a low-rank approximation, in error and in storage cost (in numbers): in a hand-built 6x5 matrix the total error of the rank 2 approximation is 0.6746 and this number equals exactly the root of the squares of the discarded singular values, but at rank 3 the storage cost rises to 36 and exceeds the original 30 numbers. In the 252x12 subscriber table reconstruction with column means gives 3.3411 while the rank 2 approximation gives 2.7464 and stores 554 numbers instead of 3,024; storing two random raw columns with the same budget gives 3.0617, and the advantage of the low-rank approximation is 0.3153 at rank 2 and 0.9663 at rank 8. The previous lesson found the component directions through the covariance matrix and projected each point onto k directions and reconstructed it. The measure was the mean distance per point, and reconstruction with column means gave 3.3411. The question the procedure did not ask was this: how many numbers does it take to store this representation. A table of 252 rows with twelve columns is 3,024 numbers; with how many does the two-component representation hold the same table. The moment the question is asked, the projection turns into a **matrix operation**. The table is a matrix whose rows are subscribers and whose columns are measurements; the projection is rewriting this matrix as the product of two small matrices. The **singular value decomposition** is the decomposition that does this: it splits the matrix into three — left singular vectors, right singular vectors and **singular values** that give a weight to each pair of directions. Keeping the first r singular values and discarding the rest produces a **low-rank approximation**. What this lesson measures is the two sides of that approximation: how much the error rises and how much the storage cost falls. - **DR9.** The matrix measured is the test set: 252 rows, 12 columns, 3,024 numbers. The columns are scaled with the scale of the previous lesson, that is, with the mean and standard deviation read **from the training set only**. - **DR10.** The decomposition is written within the lesson and is the **model**: the column product matrix of the matrix is built, the largest direction is searched for in 300 iterations, each new direction is made orthogonal to the previous ones, and the singular value is taken as the root of the magnitude found. No library is called and the scale is kept small. The decomposition is built **on the measured matrix itself**; the component directions of the previous lesson were learned from the training set and applied to the test set, so the reconstruction errors of the two lessons at the same number of dimensions do not match exactly. - **DR11.** Singular values are non-negative and sorted from large to small. The **rank** is the number of nonzero singular values; the threshold is taken as $10^{-10}$, because the iterative search produces zero not exactly but approximately. - **DR12.** The **storage cost (in numbers)** is counted explicitly: for rank r the left singular vectors are $252r$, the right singular vectors $12r$, the singular values $r$, and to return the table to its original unit 12 means and 12 standard deviations. The total is $265r + 24$. - **DR13.** The definition of the **reconstruction error** carries over from the previous lesson: the mean distance per point, in the scaled space, with the standard deviation as the unit. - **DR14.** The baseline is again **reconstruction with column means** and it is the rank 0 row. - **DR15.** The small matrix is built by hand, does not come from the subscriber table, and exists only to make visible what the decomposition does. - **DR16.** The second comparison is **storing r raw columns with the same budget**: r random columns are kept, the remaining columns are reconstructed with the mean, and three draws are averaged. #### The Decomposition on a Small Matrix A matrix of six rows and five columns is built by hand. Most of the rows are the weighted sum of two basic patterns: one high at the ends and low in the middle, the other the exact opposite. The first cell of the third row has been deliberately shifted by one unit, so the matrix does not fit the two patterns exactly. ```python def multiply(m, v): return [sum(m[i][j] * v[j] for j in range(len(v))) for i in range(len(m))] def orthogonalize(v, B): for u in B: c = sum(v[i] * u[i] for i in range(len(v))) v = [v[i] - c * u[i] for i in range(len(v))] return v def decompose(A): # MODEL: the singular value decomposition is written within the lesson, using only the standard library. m = len(A[0]) G = [[sum(s[i] * s[j] for s in A) for j in range(m)] for i in range(m)] right, singular = [], [] for _ in range(m): v, lam = [1.0 / (i + 1) for i in range(m)], 0.0 for _ in range(300): w = orthogonalize(multiply(G, v), right) u = sum(q * q for q in w) ** 0.5 if u < 1e-12: break v, lam = [q / u for q in w], u if lam < 1e-10: break right.append(v) singular.append(lam ** 0.5) left = [[sum(s[j] * v[j] for j in range(m)) / singular[i] for s in A] for i, v in enumerate(right)] return left, singular, right def approximate(left, singular, right, r, a, b): return sum(left[i][a] * singular[i] * right[i][b] for i in range(r)) SMALL = [[6, 3, 0, 3, 6], [4, 3, 2, 3, 4], [3, 3, 4, 3, 2], [0, 3, 6, 3, 0], [4, 4, 4, 4, 4], [2, 2, 2, 2, 2]] LEFT, SINGULAR, RIGHT = decompose([[float(q) for q in s] for s in SMALL]) print("singular values: " + ", ".join(f"{s:.4f}" for s in SINGULAR)) print(f"{'rank':>4} {'total error':>12} {'root of discarded':>17} {'stored':>9}") for r in range(len(SINGULAR) + 1): e = sum((SMALL[a][b] - approximate(LEFT, SINGULAR, RIGHT, r, a, b)) ** 2 for a in range(6) for b in range(5)) ** 0.5 print(f"{r:>4} {e:>12.4f} {sum(s * s for s in SINGULAR[r:]) ** 0.5:>17.4f} " f"{r * (6 + 5 + 1):>9}") print("original matrix 30 numbers. rank 2 approximation:") for a in range(6): print(" " + " ".join(f"{approximate(LEFT, SINGULAR, RIGHT, 2, a, b):6.2f}" for b in range(5))) ``` ``` singular values: 17.0637, 7.3058, 0.6746 rank total error root of discarded stored 0 18.5742 18.5742 0 1 7.3369 7.3369 12 2 0.6746 0.6746 24 3 0.0000 0.0000 36 original matrix 30 numbers. rank 2 approximation: 6.03 2.99 0.01 2.99 5.98 4.07 2.98 2.02 2.98 3.95 2.54 3.12 3.90 3.12 2.34 0.16 2.96 6.04 2.96 -0.12 4.13 3.97 4.03 3.97 3.90 2.06 1.98 2.01 1.98 1.95 ``` Two columns are exactly equal to each other. The total error of the rank r approximation is **the root of the sum of the squares of the discarded singular values** — at rank 2 the error is 0.6746 and the single discarded singular value is also 0.6746. This is the most useful thing the decomposition gives: the error is known before the approximation is built, by looking only at the list of singular values. There is no need to run the procedure and measure. The third column reverses the reading. The rank 2 approximation stores 24 numbers, the original matrix 30; the gain is 6 numbers and the error is 0.6746. The rank 3 approximation is error-free but stores 36 numbers — 6 more than the original matrix. **A full decomposition is not compression.** The gain comes only from the discarded rank, and if the matrix itself is small it may not come at all. The output of the approximation can be read too. The first, second, fifth and sixth rows come back almost exactly, because they really were sums of two patterns. In place of `3, 3, 4, 3, 2` the shifted third row returns `2.54, 3.12, 3.90, 3.12, 2.34`. The whole of the error has collected in this single row. The low-rank approximation does not distribute the error evenly; it **corrects the row that does not fit the pattern** and discards the information about that row. #### Rank and Stored Count in the Subscriber Table The same decomposition is applied to the subscriber table of 252 rows and 12 columns. The columns are scaled with the scale of the previous lesson, so that the baseline does not change. ```python import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 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 SUBSCRIBERS, DATA = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] 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({"id": 10001 + i, "region": b, "household": household, "satisfaction": rating}) for k in SUBSCRIBERS: r, v = rng(SEED + 7000 + k["id"]), [] for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): v.append(0.0 if r() < 0.038 else math.floor( k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - n * 0.05) * 100 + 0.5) / 100) r() if not v: continue mean = sum(v) / len(v) x = {"household": k["household"], "satisfaction": k["satisfaction"], "period": len(v), "mean_consumption": round(mean, 2), "volatility": round(max(v) - min(v), 2), "zero_readings": int(min(v) == 0.0), "per_person": round(mean / k["household"], 3)} for b in REGION: x["area_" + b[0]] = int(k["region"][0] == b[0]) DATA.append(x) FIELDS = ["mean_consumption", "volatility", "per_person", "household", "satisfaction", "period", "zero_readings", "area_north", "area_south", "area_east", "area_west", "area_central"] r, s = rng(SEED + 90000), 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] K = [DATA[i] for i in s] TRAIN, TEST = K[:756], K[1008:] n = len(FIELDS) MEAN = {a: sum(x[a] for x in TRAIN) / len(TRAIN) for a in FIELDS} STD_DEV = {a: max(1e-9, (sum((x[a] - MEAN[a]) ** 2 for x in TRAIN) / len(TRAIN)) ** 0.5) for a in FIELDS} Z = [[(x[a] - MEAN[a]) / STD_DEV[a] for a in FIELDS] for x in TEST] N = len(Z) LEFT, SINGULAR, RIGHT = decompose(Z) def error(r): return sum(sum((Z[a][b] - approximate(LEFT, SINGULAR, RIGHT, r, a, b)) ** 2 for b in range(n)) ** 0.5 for a in range(N)) / N print(f"original table {N} x {n} = {N * n} numbers, rank found {len(SINGULAR)}") print("singular values: " + ", ".join(f"{s:.2f}" for s in SINGULAR)) print(f"{'rank':>4} {'error':>8} {'difference':>10} {'stored':>9} {'ratio':>6}") for r in (0, 1, 2, 3, 4, 6, 8, 10, 11): print(f"{r:>4} {error(r):>8.4f} {error(r) - error(0):>+10.4f} " f"{24 + r * (N + n + 1):>9} {(24 + r * (N + n + 1)) / (N * n):>6.3f}") ``` ``` original table 252 x 12 = 3024 numbers, rank found 11 singular values: 24.19, 21.44, 19.46, 18.48, 18.14, 16.64, 15.66, 14.15, 10.50, 7.14, 5.69 rank error difference stored ratio 0 3.3411 +0.0000 24 0.008 1 3.0129 -0.3282 289 0.096 2 2.7464 -0.5947 554 0.183 3 2.4284 -0.9127 819 0.271 4 2.1363 -1.2048 1084 0.358 6 1.5054 -1.8357 1614 0.534 8 0.7925 -2.5486 2144 0.709 10 0.2568 -3.0843 2674 0.884 11 0.0000 -3.3411 2939 0.972 ``` The singular values in this table never drop steeply: the largest is 24.19, the smallest 5.69, and there is no cliff in between. In the small matrix two singular values were 17.06 and 7.31 while the third fell to 0.67; there rank 2 was a natural cut-off point. Here there is no such point. This means the cut comes not from the data but from **the decision**, and the cost of the decision can be read in the table. The rank 2 approximation stores 554 numbers instead of 3,024 — 0.183 of the original — and pulls the error 0.5947 below the baseline. When the rank is raised to 8 the error falls to 0.7925, but the storage cost rises to 2,144, that is, to 0.709 of the original. At rank 11 the error is zero and the storage is 2,939 numbers; only 85 numbers fewer than the original 3,024. Error-free compression practically does not exist in this table, because **the rank is eleven, and the only difference between it and twelve columns is the sum of the region indicators** — the error of the previous lesson that fell to zero at the eleventh direction finds its name here as the rank. #### Storing Raw Columns with the Same Budget The storage cost is a budget, and there is another way to spend it: keeping some of the raw columns as they are and reconstructing the rest with the mean. This resembles the decision measured in the feature selection topic of the Data Preparation and Feature Engineering course, but the criterion here is different — the question is not which column discriminates but which column reconstructs the table. So that the comparison is meaningful, the columns to keep are chosen at random and three draws are averaged. ```python def column_baseline(r, seed): rr, idx = rng(seed), list(range(n)) for i in range(n - 1, 0, -1): j = int(rr() * (i + 1)) idx[i], idx[j] = idx[j], idx[i] keep = set(idx[:r]) return sum(sum(z[j] * z[j] for j in range(n) if j not in keep) ** 0.5 for z in Z) / N print(f"{'r':>3} {'low rank':>11} {'stored':>9} {'r raw columns':>13} " f"{'stored':>9} {'difference':>10}") for r in (1, 2, 4, 6, 8, 11): col = sum(column_baseline(r, SEED + 500 * (t + 1)) for t in range(3)) / 3 print(f"{r:>3} {error(r):>11.4f} {24 + r * (N + n + 1):>9} {col:>13.4f} " f"{r * N + r:>9} {col - error(r):>+10.4f}") ``` ``` r low rank stored r raw columns stored difference 1 3.0129 289 3.2008 253 +0.1879 2 2.7464 554 3.0617 506 +0.3153 4 2.1363 1084 2.6787 1012 +0.5423 6 1.5054 1614 2.2730 1518 +0.7676 8 0.7925 2144 1.7588 2024 +0.9663 11 0.0000 2939 0.7857 2783 +0.7857 ``` The two paths spend almost the same budget: 554 numbers for rank 2, 506 for two raw columns. The difference is 48 numbers, less than ten percent of the budget. The error difference received in return is 0.3153 at rank 2. The low-rank approximation wins, but the gain is small — and this is an important finding for this table, because the advantage of dimensionality reduction over storing raw columns is assumed to be **automatic**. The advantage widens as the rank grows: 0.5423 at 4, 0.7676 at 6, 0.9663 at 8. The reason sits in the list of singular values. At each rank the low-rank approximation takes the **best** part of the matrix that can be carried at that rank; storing raw columns takes the columns as they are and never uses the shared structure among the columns. As the rank grows this shared structure accumulates. At rank 11 the difference begins to reverse (0.7857), because rank 11 is already error-free and the other side carries the error of dropping a single column. The decision is read from these two tables. Using a low-rank approximation for compression means being able to know the error from the singular values in advance and getting the most reconstruction for every number. But when the table is narrow and the rank is high, the gain is large neither in storage nor in error — the fall of 3,024 numbers to 554 is paid for with the error moving away from the baseline by only 0.5947. #### Summary - The singular value decomposition writes the matrix as left singular vectors, right singular vectors and singular values; keeping the first r terms is the low-rank approximation. - The reconstruction error is read from the discarded singular values: in the small matrix the total error of rank 2 is 0.6746 and the single discarded singular value is also 0.6746. - A full decomposition is not compression: in the small matrix the rank 3 approximation is error-free but stores 36 numbers, while the original matrix is 30 numbers. - In the subscriber table the baseline is 3.3411; the rank 2 approximation gives 2.7464 and stores 554 numbers instead of 3,024. Rank 11 is error-free and saves only 85 numbers. - Storing random raw columns with the same budget gives 3.0617 at rank 2; the advantage of the low-rank approximation is 0.3153 there and rises only to 0.9663 at rank 8. #### Next Step In the two lessons the projection was always of the same form: every new number was a **weighted sum** of the old columns. This constraint made the error computable and let it be read from the singular values, but it has a cost — lines go only to lines. If the points stand not on a plane but on a curved surface, no projection built with a weighted sum can follow that curve, and neighboring points land far from each other in the projection. The next lesson takes up approaches that remove this constraint: embeddings that aim to preserve not distances but **neighborhoods**. The measure changes too, because the reconstruction error is undefined there — the question to ask is how many of the nearest neighbors in the original space remain neighbors in the embedding. ### Nonlinear Embedding Source: https://academia.sh/en/courses/unsupervised-learning/nonlinear-embedding Measuring two-dimensional layouts that aim to preserve neighborhoods against the linear projection: a random layout preserves 0.0390 of the ten nearest neighbors, the linear projection 0.2639, and the neighborhood-targeted embedding 0.5032. The same embedding preserves only 0.0071 of the ten farthest pairs while the linear projection preserves 0.6012; the two methods preserve the opposite of each other. After the mean pairwise distance is equalized, the original distance of pairs standing at the same distance in the embedding varies between 3.47 and 9.46, and in a second run with the same setting the embedded distance of the same pair wanders by a mean relative difference of 0.2906. In the two lessons the projection was always of the same form: every new number was a weighted sum of the old columns. This constraint made the error computable and let it be read from the singular values. It had a cost too — a weighted sum projects onto a flat subspace, and if the points do not stand on a flat subspace, that plane cannot follow them. Two points standing close to each other on a curved surface may land far from each other in a projection that presses the surface onto a plane. The approaches this lesson covers remove that constraint. What they share is this: **they target neighborhoods, not distances**. Whom each point is close to is determined once in the original space, and then a layout is sought in a low-dimensional space that rebuilds that neighborhood relation as well as it can. There are dozens of distinct procedures in the family and they differ in the details; here the common core is written and measured. The measure has to change too, because the reconstruction error is undefined here: there is no path from the embedding back to the original space. The question to ask is this — **of the k nearest neighbors in the original space, how many remain neighbors in the embedding.** - **DR17.** The set measured is the test set: 252 points, 12 scaled columns. The scale is again read from the training set only. - **DR18.** The **neighborhood preservation rate** is the average, over all points, of how many of each point's $k$ nearest neighbors in the original space are also among its $k$ nearest neighbors in the embedding. $k = 10$ is taken. - **DR19.** The baseline is **the random layout**: 252 points are scattered randomly on the plane. Its expected value is $k / (N-1)$ and three draws are averaged. This baseline was chosen because it directly gives the neighborhood preservation of a layout that does not learn. - **DR20.** The second comparison is **the linear projection**: the first two and the first three of the component directions of the first lesson. Two dimensions is the basis, so that a comparison can be made at the same number of dimensions. - **DR21.** The layout is written within the lesson and is the **model**: it is started randomly, at every step neighbors are pulled toward each other and randomly chosen non-neighbor points are pushed away, and the step size falls linearly to zero. No library is called. - **DR22.** For the pushing, 10 random samples are drawn per point. If all pairs were scanned, the cost would grow with the square of the number of points. - **DR23.** The tuning budget is **six candidates**: two step counts with three step sizes. All the candidates tried are printed. - **DR24.** The selection is made **with the measured number itself** and no separate set is used. This means the number takes part of its value from the criterion itself; the range of the six candidates is printed for this reason. - **DR25.** The embedding has no unit. In the distance comparisons the embedding is scaled so that its mean pairwise distance equals that of the original space; otherwise the comparison would measure the choice of scale. - **DR26.** Two runs are made with the same setting from two separate starting points. The sensitivity to the starting point is part of the measure. #### The Neighborhood Preservation Rate and Its Baseline The measure is the intersection of two sets. For a point, the ten nearest neighbors in the original space are one set and the ten nearest neighbors in the embedding are a second; the size of the intersection is divided by ten and averaged over all points. In a random layout the expected value is as much as ten neighbors being chosen at random among 251 candidates. ```python # COMMON — MODEL: the FICTIONAL subscriber table of the Introduction to Machine Learning course is rebuilt with the same seed. # The label column is never generated in this course. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 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 SUBSCRIBERS, DATA = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] 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({"id": 10001 + i, "region": b, "household": household, "satisfaction": rating}) for k in SUBSCRIBERS: r, v = rng(SEED + 7000 + k["id"]), [] for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): v.append(0.0 if r() < 0.038 else math.floor( k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - n * 0.05) * 100 + 0.5) / 100) r() if not v: continue mean = sum(v) / len(v) x = {"household": k["household"], "satisfaction": k["satisfaction"], "period": len(v), "mean_consumption": round(mean, 2), "volatility": round(max(v) - min(v), 2), "zero_readings": int(min(v) == 0.0), "per_person": round(mean / k["household"], 3)} for b in REGION: x["area_" + b[0]] = int(k["region"][0] == b[0]) DATA.append(x) FIELDS = ["mean_consumption", "volatility", "per_person", "household", "satisfaction", "period", "zero_readings", "area_north", "area_south", "area_east", "area_west", "area_central"] r, s = rng(SEED + 90000), 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] K = [DATA[i] for i in s] TRAIN, TEST = K[:756], K[1008:] n = len(FIELDS) MEAN = {a: sum(x[a] for x in TRAIN) / len(TRAIN) for a in FIELDS} STD_DEV = {a: max(1e-9, (sum((x[a] - MEAN[a]) ** 2 for x in TRAIN) / len(TRAIN)) ** 0.5) for a in FIELDS} Z_TRAIN = [[(x[a] - MEAN[a]) / STD_DEV[a] for a in FIELDS] for x in TRAIN] Z = [[(x[a] - MEAN[a]) / STD_DEV[a] for a in FIELDS] for x in TEST] N = len(Z) ``` ```python def multiply(m, v): return [sum(m[i][j] * v[j] for j in range(len(v))) for i in range(len(m))] def orthogonalize(v, B): for u in B: c = sum(v[i] * u[i] for i in range(len(v))) v = [v[i] - c * u[i] for i in range(len(v))] return v def find_directions(Y): C = [[sum(z[i] * z[j] for z in Y) / len(Y) for j in range(n)] for i in range(n)] B = [] for _ in range(n): v, lam = [1.0 / (i + 1) for i in range(n)], 0.0 for _ in range(200): w = orthogonalize(multiply(C, v), B) u = sum(q * q for q in w) ** 0.5 if u < 1e-9: break v, lam = [q / u for q in w], u if lam < 1e-9: break B.append(v) return B COMPONENTS, NEIGHBOR_COUNT = find_directions(Z_TRAIN), 10 def distance(P, a, b): return sum((P[a][i] - P[b][i]) ** 2 for i in range(len(P[a]))) ** 0.5 def nearest(P, k): return [{b for _, b in sorted((distance(P, a, b), b) for b in range(len(P)) if b != a)[:k]} for a in range(len(P))] NEIGHBORS = nearest(Z, NEIGHBOR_COUNT) def preservation(P, k=NEIGHBOR_COUNT): q = nearest(P, k) return sum(len(NEIGHBORS[a] & q[a]) for a in range(N)) / (N * k) def project(d): return [[sum(z[j] * COMPONENTS[i][j] for j in range(n)) for i in range(d)] for z in Z] def random_layout(seed, d=2): r = rng(seed) return [[r() * 2 - 1 for _ in range(d)] for _ in range(N)] rs = sum(preservation(random_layout(SEED + 7 * (t + 1))) for t in range(3)) / 3 print(f"points {N}, columns {n}, neighbor count {NEIGHBOR_COUNT}") print(f"baseline: random layout {rs:.4f} (expected {NEIGHBOR_COUNT / (N - 1):.4f})") print(f"linear projection 2 dimensions {preservation(project(2)):.4f}") print(f"linear projection 3 dimensions {preservation(project(3)):.4f}") ``` ``` points 252, columns 12, neighbor count 10 baseline: random layout 0.0390 (expected 0.0398) linear projection 2 dimensions 0.2639 linear projection 3 dimensions 0.3849 ``` The baseline is 0.0390 and it sits on the theoretical expectation, that is, 0.0398. The linear projection gives 0.2639 in two dimensions: 0.2249 above the baseline, but only about a quarter of the ten nearest neighbors. In the distance table of the first lesson, two components had brought the mean of the pairwise distances down to 0.4965 times; this row is the same loss in terms of neighborhoods. Distances do not shorten proportionally, and their ordering is disturbed too. #### A Layout That Targets Neighborhoods The procedure has three parts. For each point the ten nearest neighbors in the original space are determined once. A random layout is started on the plane. Then at each step two forces are applied: neighbors are pulled toward each other, and randomly chosen non-neighbor points are pushed away. The pushing force falls with distance, so that pairs that already stand far apart do not push each other to infinity. The step size falls linearly to zero. ```python def embed(seed, steps, eta, k=NEIGHBOR_COUNT, neg=10): # MODEL: the neighborhood-preserving layout is written at small scale, using only the standard library. r = rng(seed) Y = [[(r() - 0.5) * 0.1 for _ in range(2)] for _ in range(N)] nl = [sorted(NEIGHBORS[a]) for a in range(N)] for it in range(steps): h = eta * (1 - it / steps) for a in range(N): gx = gy = 0.0 for b in nl[a]: gx += Y[a][0] - Y[b][0] gy += Y[a][1] - Y[b][1] for _ in range(neg): c = int(r() * N) if c == a or c in NEIGHBORS[a]: continue dx, dy = Y[a][0] - Y[c][0], Y[a][1] - Y[c][1] w = 1.0 / (1.0 + dx * dx + dy * dy) gx -= w * dx gy -= w * dy Y[a][0] -= h * gx / k Y[a][1] -= h * gy / k return Y print(f"{'steps':>5} {'step size':>10} {'preserved':>9}") best = None for steps in (200, 400): for eta in (0.05, 0.10, 0.20): d = preservation(embed(SEED + 11, steps, eta)) print(f"{steps:>5} {eta:>10.2f} {d:>9.4f}") if best is None or d > best[0]: best = (d, steps, eta) print(f"candidates tried 6, chosen steps {best[1]}, step size {best[2]:.2f}") EMBEDDING = embed(SEED + 11, best[1], best[2]) SECOND = embed(SEED + 22, best[1], best[2]) print(f"embedding {preservation(EMBEDDING):.4f}, second run with the same setting {preservation(SECOND):.4f}") ``` ``` steps step size preserved 200 0.05 0.4099 200 0.10 0.4869 200 0.20 0.4286 400 0.05 0.5032 400 0.10 0.4377 400 0.20 0.3940 candidates tried 6, chosen steps 400, step size 0.05 embedding 0.5032, second run with the same setting 0.4806 ``` Three numbers stand side by side. The baseline 0.0390, the linear projection 0.2639, the embedding 0.5032. The embedding's difference from the linear projection is 0.2393, and this is about the same size as the difference from the baseline to the linear projection (0.2249). In the same two dimensions, targeting neighborhoods almost doubles the neighborhood preservation. Where the difference comes from is read in the tuning table. The best of the six candidates is 0.5032, the worst 0.3940. Even the worst candidate exceeds the linear projection by 0.1301, so the advantage does not come from the tuning choice. But because the chosen number was chosen with the measured number itself, 0.5032 is an upper bound; the honest range is 0.3940 to 0.5032. The sensitivity to the starting point points the same way: the second run with the same setting gives 0.4806. As with clustering, an embedding produces not a single number but a range of numbers. #### What Distance in the Embedding Says A high neighborhood preservation rate does not mean the embedding preserves distances. They are separate questions, and the embedding answers only the first. There are two ways to measure this: does the embedding preserve the **farthest** pairs too, and how scattered are the original distances of pairs that stand at the same distance in the embedding. ```python def farthest(P, k): return [{b for _, b in sorted((-distance(P, a, b), b) for b in range(len(P)) if b != a)[:k]} for a in range(len(P))] FARTHEST = farthest(Z, NEIGHBOR_COUNT) def mean_distance(P): return sum(distance(P, a, b) for a in range(N) for b in range(a + 1, N)) / (N * (N - 1) / 2) ORIG = mean_distance(Z) print(f"original mean pairwise distance {ORIG:.3f}") for name, P in (("linear 2", project(2)), ("embedding", EMBEDDING), ("second run", SECOND)): q = farthest(P, NEIGHBOR_COUNT) far = sum(len(FARTHEST[a] & q[a]) for a in range(N)) / (N * NEIGHBOR_COUNT) scale = ORIG / mean_distance(P) Q = [[y * scale for y in p] for p in P] band = sorted(distance(Z, a, c) for a in range(N) for c in range(a + 1, N) if 0.95 * ORIG <= distance(Q, a, c) <= 1.05 * ORIG) print(f"{name:<13} far preservation {far:.4f} band {len(band):>4} pairs original distance " f"{band[0]:.2f}-{band[-1]:.2f} percentile 10-90 {band[len(band) // 10]:.2f}-" f"{band[len(band) * 9 // 10]:.2f}") f1, f2 = ORIG / mean_distance(EMBEDDING), ORIG / mean_distance(SECOND) Q1 = [[y * f1 for y in p] for p in EMBEDDING] Q2 = [[y * f2 for y in p] for p in SECOND] s, c, large = 0.0, 0, 0 for a in range(N): for b in range(a + 1, N): d1, d2 = distance(Q1, a, b), distance(Q2, a, b) o = abs(d1 - d2) / max(d1, d2, 1e-9) s, c, large = s + o, c + 1, large + (o > 0.5) print(f"relative difference of the same pair between the two runs: mean {s / c:.4f}, " f"fraction of pairs deviating by more than half {large / c:.4f}") ``` ``` original mean pairwise distance 4.695 linear 2 far preservation 0.6012 band 1899 pairs original distance 2.87-7.20 percentile 10-90 3.93-5.74 embedding far preservation 0.0071 band 2085 pairs original distance 3.47-9.46 percentile 10-90 4.04-6.79 second run far preservation 0.0000 band 5439 pairs original distance 3.39-9.21 percentile 10-90 3.97-6.01 relative difference of the same pair between the two runs: mean 0.2906, fraction of pairs deviating by more than half 0.2088 ``` The first column is unexpectedly sharp. The linear projection preserves 0.6012 of the ten farthest pairs; the embedding 0.0071, and the second run none. **The two approaches preserve the opposite of each other.** The linear projection carries large distances because it has chosen the directions that carry the most variability, and that is exactly what separates far pairs; it crushes close neighbors. The embedding rescues close neighborhoods and in return throws away the large scale of distance entirely. If two points stand far apart in the embedding, that does not mean they are far apart in the original space — it only means they are not neighbors. The second measurement says the same thing with a single range. After the mean pairwise distance is equalized, the original distance of the 2,085 pairs that stand within five percent of 4.695 in the embedding varies between 3.47 and 9.46; even the range between the 10th and 90th percentiles is 4.04–6.79. In the same band the range of the linear projection is 2.87–7.20 and its 10th–90th percentile range is 3.93–5.74, that is, markedly narrower. To look at the distance between two points in an embedding and say "they differ this much" is to read a random number from a range whose largest is 2.7 times its smallest. The last line closes the decision. Changing only the starting point with the same setting moves the embedded distance of the same pair by a mean relative difference of 0.2906, and in 0.2088 of the pairs the difference is larger than half. The coordinates of an embedding, and the distances between the coordinates, depend on the run. The only number that can be reported for an embedding is the neighborhood preservation rate; axes, directions and distances cannot be reported. #### Summary - The baseline of the neighborhood preservation rate is the random layout and it is 0.0390; with the theoretical expectation 0.0398 it sits on it. - The linear projection gives 0.2639 in two dimensions and the embedding that targets neighborhoods 0.5032; the difference is 0.2393, and since even the worst of the six candidates, 0.3940, exceeds the linear projection, it does not come from the tuning choice. - The number is a range: the tuning range is 0.3940–0.5032, and the second run with the same setting is 0.4806. - The two approaches preserve the opposite of each other: the linear projection preserves 0.6012 of the farthest pairs, the embedding 0.0071. - The embedding makes no claim of preserving distances: after the mean distance is equalized, the original distance of pairs standing in the same band in the embedding varies between 3.47 and 9.46, and in the second run the distance of the same pair wanders by a mean relative difference of 0.2906. #### Next Step In three lessons the points were reduced to fewer numbers, and each time what was lost was measured: the reconstruction error, the storage cost, the neighborhood preservation rate. The three measures had one thing in common — all of them were averaged **over all points**. A projection was counted good if it was good on the average of 252 points, and an embedding was counted as preserving half if it preserved half of the neighborhoods. The single points that stay below the average — a subscriber, for example, whose reconstruction error is ten times that of the others — were never looked at. Yet such a point may itself be a finding: a record that does not fit the table and stands alone. The next lesson turns this question into a measure — how a point that stands alone is flagged, and how many of those flagged really are so. ### Anomaly Detection Source: https://academia.sh/en/courses/unsupervised-learning/anomaly-detection Measuring flagging on 63 anomalies deliberately planted into the 1,260 rows of the fictional subscriber table: at a 5% rate random flagging gives a precision of 0.0488, flagging by a single column 0.3968, the procedure that looks at each column on its own 0.6190, and the reconstruction error with four components 0.4921. The precision ranking is misleading — the three detectors with the highest precision find none of the 24 multivariate anomalies that go to the extreme in no column, while the one with the lowest precision finds fifteen of them, and when the contamination rate rises to 10% it finds all twenty-four. At a rate of 2% three detectors give a precision of 1.0000 and at the same time catch only 25 of the 63 anomalies. The first three lessons of this topic did the same job in three different ways: the points were reduced to fewer numbers and each time what was lost was measured — the reconstruction error, the cost of rank, the neighborhood preservation rate. All three measures were read by taking **all points together**. What a point that stands alone is, or even whether such a point exists, was not asked. Anomaly detection asks this question, and asks it as an unlabeled problem: no column says that a row is an anomaly, so the definition of "anomaly" comes from the method itself. In such a setup the only way to say which flagging is good is to be able to count **how many of the flagged are real**. This lesson does that by **deliberately planting** anomalies into the fictional dataset. - **DR31.** The dataset is the **fictional** subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering, and Supervised Learning courses, seed 20260218, 1,260 rows and five numeric columns. The label column is not read in this course and is not computed here either. - **DR32.** The anomalies are **deliberately planted** into the table. Which row was planted is known **only for measurement**; no detector sees this list. This is the place of the fiction: in a real dataset this list does not exist and the numbers of this lesson cannot be computed there. - **DR33.** `anomaly` and `outlier` are kept apart. Outlier is the term of the Data Collection and Preparation course and describes **a single value in one column** standing at the extreme. Anomaly here is **the row as a whole** not fitting the dataset; a row with no value at the extreme may be an anomaly. - **DR34.** The baseline is the interquartile range, fence and standard score procedure of the Data Collection and Preparation course. The procedure was built there and is **not repeated** here; because a fixed number of rows must be flagged, instead of the fence the continuous form of the same procedure, the magnitude of the standard score, is called as the ranking criterion. - **DR35.** The **contamination rate** is a hyperparameter: how many rows will be counted as anomalies does not come out of the data. It is swept and the number of each rate is written separately. #### The Planted Anomaly ```python # COMMON — FICTIONAL. The subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering and Supervised Learning courses is rebuilt with the same seed; # the label column is never read in this course. The anomalies are PLANTED into the table ON PURPOSE. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 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 SUBSCRIBERS, DATA = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] 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({"id": 10001 + i, "region": b, "household": household, "satisfaction": rating}) for k in SUBSCRIBERS: r, v = rng(SEED + 7000 + k["id"]), [] for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): v.append(0.0 if r() < 0.038 else math.floor( k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - n * 0.05) * 100 + 0.5) / 100) r() if v: DATA.append({"id": k["id"], "household": k["household"], "satisfaction": k["satisfaction"], "period": len(v), "mean_consumption": round(sum(v) / len(v), 2), "volatility": round(max(v) - min(v), 2)}) # Planting. Three types: A and B are moved to the extreme in a single column, C breaks the # relationship between columns without going to the extreme in any column (three readings, with almost no difference among them). r, taken, PLANTED = rng(SEED + 31000), set(), {} def choose(cond, count): s = [i for i in range(len(DATA)) if cond(DATA[i]) and i not in taken] for i in range(len(s) - 1, 0, -1): j = int(r() * (i + 1)) s[i], s[j] = s[j], s[i] taken.update(s[:count]) return s[:count] for i in choose(lambda x: 12 <= x["mean_consumption"] <= 30 and x["period"] >= 2, 24): DATA[i]["mean_consumption"], PLANTED[i] = round(DATA[i]["mean_consumption"] * 5.0, 2), "A" for i in choose(lambda x: x["period"] == 3 and 12 <= x["volatility"] <= 20, 15): DATA[i]["volatility"], PLANTED[i] = round(DATA[i]["volatility"] * 5.0, 2), "B" for i in choose(lambda x: x["period"] == 3 and 22 <= x["mean_consumption"] <= 33, 24): DATA[i]["volatility"], PLANTED[i] = round(DATA[i]["mean_consumption"] * 0.02, 2), "C" FIELDS = ["mean_consumption", "volatility", "household", "satisfaction", "period"] n, p = len(DATA), len(FIELDS) MEAN = [sum(x[a] for x in DATA) / n for a in FIELDS] STD_DEV = [math.sqrt(sum((x[a] - MEAN[j]) ** 2 for x in DATA) / n) for j, a in enumerate(FIELDS)] P = [[(x[a] - MEAN[j]) / STD_DEV[j] for j, a in enumerate(FIELDS)] for x in DATA] print(f"fictional table {n} rows, {p} columns; planted anomalies {len(PLANTED)} " f"({len(PLANTED) / n:.4f})") for t, name in (("A", "scale shift"), ("B", "volatility spike"), ("C", "steady high reading")): g = [i for i in PLANTED if PLANTED[i] == t] print(f" type {t} {name:<21} {len(g):>3} rows " f"mean_consumption {min(DATA[i]['mean_consumption'] for i in g):>6.2f}-" f"{max(DATA[i]['mean_consumption'] for i in g):<6.2f} " f"volatility {min(DATA[i]['volatility'] for i in g):>5.2f}-" f"{max(DATA[i]['volatility'] for i in g):.2f}") print(f" non-planted rows have mean_consumption at most " f"{max(DATA[i]['mean_consumption'] for i in range(n) if i not in PLANTED):.2f}, " f"volatility at most " f"{max(DATA[i]['volatility'] for i in range(n) if i not in PLANTED):.2f}") ``` ``` fictional table 1260 rows, 5 columns; planted anomalies 63 (0.0500) type A scale shift 24 rows mean_consumption 62.45-144.40 volatility 1.71-23.22 type B volatility spike 15 rows mean_consumption 9.50-35.29 volatility 61.10-98.55 type C steady high reading 24 rows mean_consumption 22.09-31.68 volatility 0.44-0.63 non-planted rows have mean_consumption at most 43.48, volatility at most 51.15 ``` The difference among the three types stands in the last line. The consumption of type A is between 62.45 and 144.40, while the largest among the non-planted is 43.48; the volatility of type B is between 61.10 and 98.55, while the largest among the non-planted is 51.15. Both stand **at the extreme in one column** and, in the term of the Data Collection and Preparation course, each carries an outlier in that column. Type C goes to the extreme in no column. Its consumption is between 22.09 and 31.68, in the upper quarter of the dataset but well inside the limit. Its volatility is between 0.44 and 0.63, and there are 187 rows in the dataset with volatility below 1, because in a subscriber with a single reading the volatility is by definition zero. The anomaly is **not in the values but in the togetherness of the values**: a subscriber read over three periods, with high consumption and almost no difference among the three readings. #### Five Detectors, Three Baselines ```python # Five detectors. The first three are baselines, the last two are multivariate. def components(k, seed): # Covariance matrix and the first k eigenvectors: power iteration and deflation. # The component itself was built in the first two lessons of this topic; here it is an input. C = [[sum(P[i][a] * P[i][b] for i in range(n)) / n for b in range(p)] for a in range(p)] rr, found = rng(seed), [] for _ in range(k): v = [rr() - 0.5 for _ in range(p)] for _ in range(300): w = [sum(C[a][b] * v[b] for b in range(p)) for a in range(p)] u = math.sqrt(sum(x * x for x in w)) or 1.0 v = [x / u for x in w] lam = sum(v[a] * sum(C[a][b] * v[b] for b in range(p)) for a in range(p)) found.append(v) C = [[C[a][b] - lam * v[a] * v[b] for b in range(p)] for a in range(p)] return found def reconstruction_error(k): V = components(k, SEED + 31500) s = [] for i in range(n): pr = [sum(P[i][a] * v[a] for a in range(p)) for v in V] rec = [sum(pr[j] * V[j][a] for j in range(k)) for a in range(p)] s.append(math.sqrt(sum((P[i][a] - rec[a]) ** 2 for a in range(p)))) return s FAR = [] for i in range(n): d = sorted(sum((P[i][a] - P[j][a]) ** 2 for a in range(p)) for j in range(n) if j != i) FAR.append(math.sqrt(d[9])) # distance of the 10th neighbor SCORES = [("random", "random flagging", None), ("single column", "by single column: mean_consumption", [abs(P[i][0]) for i in range(n)]), ("column by column", "each column on its own: largest standard score", [max(abs(P[i][a]) for a in range(p)) for i in range(n)]), ("neighbor", "10th neighbor distance", FAR), ("reconstruct 2", "reconstruction error with two components", reconstruction_error(2)), ("reconstruct 4", "reconstruction error with four components", reconstruction_error(4))] def flag(score, m, seed=0): if score is None: rr, s = rng(seed), list(range(n)) for i in range(n - 1, 0, -1): j = int(rr() * (i + 1)) s[i], s[j] = s[j], s[i] return s[:m] return sorted(range(n), key=lambda i: (-score[i], i))[:m] M = 63 print(f"contamination rate 0.0500 -> {M} rows are flagged; planted 63\n") print(f"{'detector':<46} {'real':>7} {'precision':>9} " f"{'A':>4} {'B':>4} {'C':>4}") for _, name, score in SCORES: if score is None: d = [sum(1 for i in flag(None, M, SEED + 32000 + t) if i in PLANTED) for t in range(200)] g = sum(d) / len(d) print(f"{name:<46} {g:>7.2f} {g / M:>9.4f} {'':>4} {'':>4} {'':>4}") continue s = flag(score, M) counts = {t: sum(1 for i in s if PLANTED.get(i) == t) for t in "ABC"} print(f"{name:<46} {sum(counts.values()):>7} {sum(counts.values()) / M:>9.4f} " f"{counts['A']:>4} {counts['B']:>4} {counts['C']:>4}") ``` ``` contamination rate 0.0500 -> 63 rows are flagged; planted 63 detector real precision A B C random flagging 3.08 0.0488 by single column: mean_consumption 25 0.3968 24 1 0 each column on its own: largest standard score 39 0.6190 24 15 0 10th neighbor distance 39 0.6190 24 15 0 reconstruction error with two components 40 0.6349 23 15 2 reconstruction error with four components 31 0.4921 1 15 15 ``` The baseline is in the first row: when 63 rows are flagged at random, the mean of 200 draws is 3.08 real anomalies, a precision of **0.0488** — the contamination rate itself. Next to this number no detector counts as having found anything. The second baseline looks at a single column and gives 0.3968 — eight times the random one. What it found is in the third column: it caught **all 24** of type A, only one of type B, none of type C. The third baseline applies the same procedure to the five columns one by one and ranks by the largest standard score; the precision rises to 0.6190 and the missing type B is completed too. Type C is still zero. The last three rows are multivariate. The distance of the tenth neighbor gives the **same** number as the column-by-column baseline: 39 real, the same distribution. The reconstruction error with two components is slightly better, 0.6349. The reconstruction error with four components, however, falls behind all of them in precision, **0.4921** — and by itself finds fifteen of type C. The precision ranking is misleading in this table. The three detectors with the highest precision see **none** of the 24 anomalies that go to the extreme in no column; the one with the lowest precision sees fifteen. The difference comes not from the sensitivity of the method but from **what it looks at**. Type A has gone far along the strongest component of the dataset, and four components reconstruct it almost exactly; its residual is close to zero. Type C, on the other hand, is far along no component; it stands **in the direction left outside the components**, because what it breaks is not a value but the relationship between the number of periods and the volatility and consumption. For the same reason the reconstruction with four components misses type A: the residual it looks for is not a residual A ever produces. #### When the Contamination Rate Is Swept ```python # The contamination rate is a hyperparameter: how many rows will be counted as anomalies # comes not from the data but from a decision. Five rates are swept. RATES = [0.02, 0.05, 0.10, 0.15, 0.20] HEADER = f"{'rate':>5} {'flagged':>7} " + " ".join(f"{a:>16}" for a, _, _ in SCORES) print("precision: how many of the flagged are planted anomalies\n" + HEADER) for o in RATES: m, row = int(round(o * n)), [] for _, _, score in SCORES: if score is None: d = [sum(1 for i in flag(None, m, SEED + 32000 + t) if i in PLANTED) for t in range(200)] row.append(sum(d) / len(d) / m) else: row.append(sum(1 for i in flag(score, m) if i in PLANTED) / m) print(f"{o:>5.2f} {m:>7} " + " ".join(f"{x:>16.4f}" for x in row)) print("\ncaught: how many of the 63 anomalies / how many of the 24 of type C\n" + HEADER) for o in RATES: m, row = int(round(o * n)), [] for _, _, score in SCORES: s = flag(score, m, SEED + 32000) if score is None else flag(score, m) row.append(f"{sum(1 for i in s if i in PLANTED)}/" f"{sum(1 for i in s if PLANTED.get(i) == 'C')}") print(f"{o:>5.2f} {m:>7} " + " ".join(f"{x:>16}" for x in row)) ``` ``` precision: how many of the flagged are planted anomalies rate flagged random single column column by column neighbor reconstruct 2 reconstruct 4 0.02 25 0.0506 0.9600 1.0000 1.0000 1.0000 0.6000 0.05 63 0.0488 0.3968 0.6190 0.6190 0.6349 0.4921 0.10 126 0.0496 0.2143 0.3175 0.3175 0.3333 0.3333 0.15 189 0.0501 0.1534 0.2116 0.2169 0.2275 0.2328 0.20 252 0.0501 0.1230 0.1587 0.1706 0.1706 0.1865 caught: how many of the 63 anomalies / how many of the 24 of type C rate flagged random single column column by column neighbor reconstruct 2 reconstruct 4 0.02 25 1/1 24/0 25/0 25/0 25/0 15/0 0.05 63 2/1 25/0 39/0 39/0 40/2 31/15 0.10 126 3/1 27/0 40/1 40/1 42/3 42/24 0.15 189 6/3 29/1 40/1 41/2 43/4 44/24 0.20 252 7/3 31/2 40/1 43/4 43/4 47/24 ``` The first row shows how little it takes to raise the precision. When the rate is pulled to 2%, three detectors give **1.0000**: all 25 of the 25 flagged rows are real. In the second table the same row says the opposite — only 25 of the 63 anomalies have been found, none of them of type C. Perfect precision is the product of a narrow threshold and has been bought by never seeing 38 anomalies. As the rate grows, the precision of every detector falls, because the number of flagged rows exceeds the number of planted anomalies and the excess is necessarily wrong. The only thing that stays constant is **random flagging**: 0.0506, 0.0488, 0.0496, 0.0501, 0.0501. The baseline is independent of the rate, and precision can be read only next to this number. The second table counts the limit of the univariate procedures. Flagging by a single column carries the caught anomalies only from 24 to 31 even at a rate of 20%; the column-by-column procedure stops at 40 and does not progress at all after 10%. In these two procedures raising the rate does not show what is unseen, it only multiplies the wrongly flagged rows. The reconstruction error with four components, however, has found **all 24** of type C at 10%. At a rate of 20% the best precision is 0.1865: 252 rows flagged, 47 of them real. The remaining 205 rows have been set aside for inspection and none of them holds a planted anomaly. Choosing the flag rate is accepting the cost of these 205 rows; the criterion does not choose the rate, the rate determines the number of the criterion. #### Summary - Of the 63 anomalies deliberately planted into the fictional table, 24 go to the extreme in no column: the anomaly is not in the values but in the togetherness of the values. - At a contamination rate of 5% the baseline is 0.0488; flagging by a single column gives 0.3968, the procedure that looks at each column on its own 0.6190, the distance of the tenth neighbor 0.6190, and the reconstruction error with four components 0.4921. - The precision ranking is misleading: the three detectors with the highest precision find none of the 24 multivariate anomalies, the one with the lowest precision finds fifteen; at a rate of 10% it finds all twenty-four. - When the rate is pulled to 2%, three detectors give a precision of 1.0000 and at the same time catch only 25 of the 63 anomalies; perfect precision is the product of a narrow threshold. - The precision of random flagging stays between 0.0488 and 0.0506 at all rates; the baseline is independent of the rate and the other numbers are read only next to it. #### Next Step In this lesson the row was handled with five numeric columns and whether a point stands alone was measured. Subscriber records have one more side that is not a number: which events occur **together**. Is the joint occurrence of two events in the same subscriber a rule, or can it already be explained by how frequent the two are separately. The next lesson builds association rules, counts how many of the high-confidence rules come only from a frequent item, and reads the baseline from the support that random matching produces. ### Association Rules Source: https://academia.sh/en/courses/unsupervised-learning/association-rules Measuring the 166 rules that pass the support floor in the 1,329 baskets of the fictional subscriber table: all six of the top six rules by confidence name the same frequent item as the consequent and their lifts are between 1.0139 and 1.0489; of the 97 rules with confidence above 0.50, 75 have a lift closer to 1 than 0.05. The baseline is the support that random matching produces — when the item columns are shuffled separately, the support of 0.1279 for delay and installment plan does not exceed 0.0467 in any of sixty shuffles, but the support of a rule with a confidence of 0.9246 stays inside the random range. The course wrap-up gathers the baseline, the method's number and the difference of eleven lessons in a single table. The previous lesson asked whether a row stands alone and looked for the answer in five numeric columns. Subscriber records have one more side that is not a number: **which events occur together** in that subscriber within the period. The question keeps its form: is the joint occurrence of two events a finding, or can it already be explained by how frequent the two are separately. Association rules answer this question with three measures, and one of them misleads when taken alone. The number of the lesson is **how many of the rules chosen by confidence come only from a frequent item**; its baseline is the support that comes out when the items are **matched to each other at random**. - **DR36.** The set is the **fictional** subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering, and Supervised Learning courses, seed 20260218, 1,329 subscribers. Each subscriber is read as an **itemset**: the service records that occur in that subscriber within the period. - **DR37.** The dependencies are **deliberately** placed in the generator, and which of them is real is known only for measurement; no measure sees this information. - **DR38.** Three measures: **support**, the share of the baskets that contain the itemset; **confidence**, the share of the baskets containing the antecedent that also contain the consequent; **lift**, the confidence divided by the support of the consequent. Lift is read as a **departure from 1**; 1 is the absence of association. - **DR39.** The baseline is **the support that random matching produces**: the column of each item is shuffled separately across the baskets, the frequencies are kept, and the bond of occurring together is broken. - **DR40.** This lesson builds only **the measures of extracting rules**; which rule will be shown to whom is left to the Reinforcement Learning and Recommender Systems course. #### Baskets and Item Supports ```python # COMMON — FICTIONAL. The subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering, and Supervised Learning courses is rebuilt with the same seed; this time the # service records of each subscriber within the period are kept as an ITEMSET. # The dependencies are PLACED in the generator ON PURPOSE and are known only for measurement. SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF 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 SUBSCRIBERS = [] for i in range(RAW): r = rng(SEED + i) r(), r(), r() # region, household and satisfaction are not read in this lesson if r() >= 0.046: SUBSCRIBERS.append(10001 + i) # (item, condition item, probability when the condition holds, probability when it does not) TABLE = [("mobile_notification", None, 0.91, 0.91), # independent, very frequent ("outage_notification", None, 0.55, 0.55), # independent, very frequent ("delay", None, 0.22, 0.22), ("installment_plan", "delay", 0.58, 0.04), # placed dependency ("auto_payment", "delay", 0.13, 0.47), # placed inverse dependency ("objection", None, 0.16, 0.16), ("call_center", "objection", 0.90, 0.58), # weak dependency, frequent item ("address_change", None, 0.08, 0.08), ("meter_replacement", "address_change", 0.36, 0.09), ("field_visit", "meter_replacement", 0.86, 0.19), ("tariff_change", None, 0.12, 0.12), ("illicit_use_report", None, 0.05, 0.05)] ITEMS = [t[0] for t in TABLE] BASKETS = [] for id in SUBSCRIBERS: r, s = rng(SEED + 37000 + id), set() for item, cond, p1, p0 in TABLE: if r() < (p1 if cond is None or cond in s else p0): s.add(item) BASKETS.append(s) N = len(BASKETS) D1 = {o: sum(1 for s in BASKETS if o in s) / N for o in ITEMS} print(f"fictional: {N} subscriber baskets, {len(ITEMS)} items; " f"mean {sum(len(s) for s in BASKETS) / N:.2f} items per basket") for o in sorted(ITEMS, key=lambda x: -D1[x]): print(f" {o:<20} support {D1[o]:.4f}") ``` ``` fictional: 1329 subscriber baskets, 12 items; mean 3.68 items per basket mobile_notification support 0.9120 call_center support 0.6388 outage_notification support 0.5403 auto_payment support 0.3973 field_visit support 0.2671 delay support 0.2220 installment_plan support 0.1670 objection support 0.1497 meter_replacement support 0.1332 tariff_change support 0.1151 address_change support 0.0835 illicit_use_report support 0.0512 ``` The supports are spread over a wide range: the most frequent item is in 0.9120 of the baskets, the rarest in 0.0512. This distribution sets the trap of the next section. #### Confidence Alone ```python # Extracting rules: one or two items in the antecedent, one item in the consequent. The support floor # is a decision and is fixed before measuring. from itertools import combinations SUPPORT_FLOOR, CONFIDENCE_FLOOR = 0.05, 0.50 def support(itemset): return sum(1 for s in BASKETS if itemset <= s) / N RULES = [] antecedents = [frozenset([o]) for o in ITEMS] + \ [frozenset(c) for c in combinations(ITEMS, 2)] for X in antecedents: dx = support(X) if dx < SUPPORT_FLOOR: continue for y in ITEMS: if y in X: continue d = support(X | {y}) if d < SUPPORT_FLOOR: continue g = d / dx RULES.append({"X": X, "y": y, "support": d, "confidence": g, "lift": g / D1[y]}) SINGLE = sum(1 for k in RULES if len(k["X"]) == 1) near = [k for k in RULES if abs(k["lift"] - 1) <= 0.05] print(f"support floor {SUPPORT_FLOOR:.2f} and above: {len(RULES)} rules " f"({SINGLE} single-item antecedents, {len(RULES) - SINGLE} two-item)") print(f"rules with lift within 0.05 of 1: {len(near)} " f"({len(near) / len(RULES):.4f})\n") def show(title, ranked, count=8): print(title) print(f" {'rule':<46} {'support':>7} {'confidence':>10} {'lift':>9}") for k in ranked[:count]: name = " + ".join(sorted(k["X"])) + " -> " + k["y"] print(f" {name:<46} {k['support']:>7.4f} {k['confidence']:>10.4f} " f"{k['lift']:>9.4f}") by_confidence = sorted(RULES, key=lambda k: (-k["confidence"], sorted(k["X"]), k["y"])) show("top six rules by confidence", by_confidence, 6) top = [k for k in by_confidence if k["confidence"] >= CONFIDENCE_FLOOR] frequent = sum(1 for k in top if abs(k["lift"] - 1) <= 0.05) print(f"\nthere are {len(top)} rules with confidence {CONFIDENCE_FLOOR:.2f} and above; " f"{frequent} of them have a lift within 0.05 of 1 ({frequent / len(top):.4f})") last = {} for k in top: last[k["y"]] = last.get(k["y"], 0) + 1 print(" by consequent item: " + ", ".join( f"{o} {last[o]} (support {D1[o]:.4f})" for o in sorted(last, key=lambda x: -last[x]))) single = [k for k in RULES if len(k["X"]) == 1] show("\ntop four rules by lift (single-item antecedent)", sorted(single, key=lambda k: (-k["lift"], sorted(k["X"]), k["y"])), 4) lowest = min(single, key=lambda k: (k["lift"], k["confidence"], sorted(k["X"]), k["y"])) print(f" lowest-lift rule: {' + '.join(sorted(lowest['X']))} -> " f"{lowest['y']}, support {lowest['support']:.4f}, " f"confidence {lowest['confidence']:.4f}, lift {lowest['lift']:.4f}") ``` ``` support floor 0.05 and above: 166 rules (64 single-item antecedents, 102 two-item) rules with lift within 0.05 of 1: 115 (0.6928) top six rules by confidence rule support confidence lift outage_notification + tariff_change -> mobile_notification 0.0662 0.9565 1.0489 tariff_change -> mobile_notification 0.1076 0.9346 1.0249 auto_payment + objection -> mobile_notification 0.0602 0.9302 1.0200 call_center + tariff_change -> mobile_notification 0.0700 0.9300 1.0198 address_change -> mobile_notification 0.0775 0.9279 1.0175 objection -> mobile_notification 0.1384 0.9246 1.0139 there are 97 rules with confidence 0.50 and above; 75 of them have a lift within 0.05 of 1 (0.7732) by consequent item: mobile_notification 32 (support 0.9120), call_center 28 (support 0.6388), outage_notification 25 (support 0.5403), field_visit 4 (support 0.2671), delay 4 (support 0.2220), installment_plan 4 (support 0.1670) top four rules by lift (single-item antecedent) rule support confidence lift installment_plan -> delay 0.1279 0.7658 3.4498 delay -> installment_plan 0.1279 0.5763 3.4498 field_visit -> meter_replacement 0.1159 0.4338 3.2572 meter_replacement -> field_visit 0.1159 0.8701 3.2572 lowest-lift rule: outage_notification -> field_visit, support 0.1362, confidence 0.2521, lift 0.9437 ``` All six of the top six rules by confidence have **the same item as the consequent**: the record that occurs in 0.9120 of the baskets. Their confidences are 0.9246–0.9565 and their antecedents have nothing to do with each other. **Whichever antecedent is chosen**, the confidence settles around 0.91; the lift column says this, between 1.0139 and 1.0489. There are 97 rules with confidence above 0.50 and **75 of them** have a lift closer to 1 than 0.05: 0.7732. By consequent item, 32, 28 and 25 rules point to the three most frequent items — 85 of the 97. Raising the threshold does not fix this; as the threshold rises, **only frequent items** can pass it. Sorting by lift gives another list: two pairs in the top four, each pair in both directions — installment plan and delay 3.4498, meter replacement and field visit 3.2572. In the same rows confidence is read as **asymmetric** and lift as **symmetric**: the confidence from delay to installment plan is 0.5763, in the opposite direction 0.7658, and the lift is the same in both. The lowest lift is 0.9437; a lift that falls below 1 says that occurring together is **rarer** than chance. #### The Random Matching Baseline ```python # BASELINE: the support that random matching produces. The column of each item is # shuffled separately across the baskets; item frequencies are kept exactly, # the bond of occurring together is broken. Lift is the closed form of this baseline. BIT = {o: 1 << b for b, o in enumerate(ITEMS)} COLUMNS = [[1 if o in s else 0 for s in BASKETS] for o in ITEMS] REPEATS = 60 def superset_sums(masks): # For each pattern, the number of baskets that cover it: a single 12-bit pass. h = [0] * (1 << len(ITEMS)) for m in masks: h[m] += 1 for b in range(len(ITEMS)): for m in range(len(h)): if not m & (1 << b): h[m] += h[m | (1 << b)] return h def shuffled(seed): rr, new = rng(seed), [0] * N for b, col in enumerate(COLUMNS): c = list(col) for i in range(N - 1, 0, -1): j = int(rr() * (i + 1)) c[i], c[j] = c[j], c[i] for i in range(N): if c[i]: new[i] |= 1 << b return new RANDOM_RUNS = [superset_sums(shuffled(SEED + 41000 + t)) for t in range(REPEATS)] EXAMPLES = (("delay", "installment_plan"), ("meter_replacement", "field_visit"), ("tariff_change", "mobile_notification"), ("objection", "mobile_notification")) pattern = lambda k: sum(BIT[o] for o in k["X"]) | BIT[k["y"]] random_support = lambda k: [R[pattern(k)] / N for R in RANDOM_RUNS] print(f"{REPEATS} shuffles: the support random matching produces\n") print(f" {'rule':<38} {'real':>7} {'random':>9} {'highest':>7} {'lift':>9}") for k in sorted((k for k in RULES if len(k["X"]) == 1 and (sorted(k["X"])[0], k["y"]) in EXAMPLES), key=lambda k: -k["lift"]): r = random_support(k) print(f" {sorted(k['X'])[0] + ' -> ' + k['y']:<38} {k['support']:>7.4f} " f"{sum(r) / len(r):>9.4f} {max(r):>7.4f} {k['lift']:>9.4f}") OUTSIDE = {i: k["support"] > max(random_support(k)) or k["support"] < min(random_support(k)) for i, k in enumerate(RULES)} print(f"\n{len(RULES)} rules separated by two criteria at once") print(f" {'':<22} {'inside range':>14} {'outside range':>15}") for name, cond in (("lift near 1", True), ("lift far from 1", False)): g = [i for i, k in enumerate(RULES) if (abs(k["lift"] - 1) <= 0.05) == cond] inside = sum(1 for i in g if not OUTSIDE[i]) print(f" {name:<22} {inside:>14} {len(g) - inside:>15}") dt = support(frozenset(["delay", "auto_payment"])) print(f"\na dependency that stays under the support floor: delay + auto_payment " f"support {dt:.4f}, lift {dt / (D1['delay'] * D1['auto_payment']):.4f}; " f"since the floor is {SUPPORT_FLOOR:.2f} it never enters the rule list") ``` ``` 60 shuffles: the support random matching produces rule real random highest lift delay -> installment_plan 0.1279 0.0362 0.0467 3.4498 meter_replacement -> field_visit 0.1159 0.0354 0.0451 3.2572 tariff_change -> mobile_notification 0.1076 0.1043 0.1099 1.0249 objection -> mobile_notification 0.1384 0.1363 0.1407 1.0139 166 rules separated by two criteria at once inside range outside range lift near 1 107 8 lift far from 1 26 25 a dependency that stays under the support floor: delay + auto_payment support 0.0316, lift 0.3584; since the floor is 0.05 it never enters the rule list ``` The baseline is in the first two rows. Delay and installment plan occur together in 0.1279 of the baskets; when the columns are shuffled, the same pair falls to a mean of 0.0362 and the **highest of sixty shuffles is 0.0467**. The difference is 0.0917 and it comes from the data. The last two rows show the opposite. The rule with a confidence of 0.9246 has a support of 0.1384; random matching carries the same pair to a mean of 0.1363 and a highest of 0.1407. **The real support is inside the random range**; the rule takes its high confidence entirely from the frequency of the consequent itself. Lift is the closed form of this baseline: dividing by the product of the supports is dividing by what shuffling produces on average. When the two criteria are crossed, the 166 rules are spread over four cells: of the 115 rules with a lift near 1, 107 are inside the random range, and of the 51 rules with a lift far from 1, 26 are as well. In total 133 rules cannot be told apart from what random matching produces. The 8 rows in the opposite direction are meaningful too: as the support grows the range narrows, and even a lift of 1.04 falls outside. A departure of the lift is not enough; **how often the departure is produced by chance** is needed as well. The last row opens another face. The fictional dependency between delay and auto payment has a lift of 0.3584, farther from 1 than anything in the list; but since the two occur together in only 0.0316 of the baskets, it stays under the support floor and **never enters** the list. If the floor were 0.03, this rule would be at the top of the list. #### Summary - All six of the top six rules by confidence have as the consequent the same item that occurs in 0.9120 of the baskets; their lifts are 1.0139–1.0489. - Of the 97 rules with confidence above 0.50, 75 have a lift closer to 1 than 0.05, and 85 have one of the three most frequent items as the consequent. - The baseline is the support that random matching produces: the support of 0.1279 for delay and installment plan does not exceed 0.0467 in any of sixty shuffles, but the support of the rule with a confidence of 0.9246 stays inside the random range. - 133 of the 166 rules cannot be told apart from what random matching produces. - The support floor is a decision: a real dependency with a lift of 0.3584 stays under the floor of 0.05 because its support is 0.0316, and never enters the list. #### Course Wrap-Up The account of eleven lessons is in a single table, and each number was read from its own lesson. The triple is the same in every row, but **the rows do not carry a common unit**: when there are no labels, each lesson chooses its own metric. The last column writes this — in the distance and error rows the smaller is better, in the rate and score rows the larger. The difference is the number of the method minus the baseline. | Lesson | Baseline | Method's number | Difference | Metric | |---|---|---|---|---| | Clustering 01 | random assignment 5.9790 | 3.8910 | −2.0880 | within-cluster distance | | Clustering 02 | random assignment 5.9181 | 2.2570 | −3.6611 | within-cluster distance | | Clustering 03 | random assignment 6.2340 | 4.5419 | −1.6922 | within-cluster distance | | Clustering 04 | random flagging 1.0200 | 1.7551 | +0.7351 | neighbor distance | | Clustering 05 | random assignment −4.2424 | −3.7755 | +0.4669 | fit score | | Clustering 06 | single cluster 0.6539 | 0.5235 | −0.1304 | agreement index | | Dimensionality reduction 01 | column mean 3.3411 | 2.7762 | −0.5650 | reconstruction error | | Dimensionality reduction 02 | column mean 3.3411 | 2.7464 | −0.5947 | reconstruction error | | Dimensionality reduction 03 | random layout 0.0390 | 0.5032 | +0.4642 | neighborhood preservation rate | | Dimensionality reduction 04 | random flagging 0.0488 | 0.6190 | +0.5702 | precision | | Dimensionality reduction 05 | random matching 0.0362 | 0.1279 | +0.0917 | support of the rule | The second column was produced in each lesson by a separate procedure: random assignment, single cluster, random flagging, reconstruction with the column mean, random layout, random matching. In the supervised course the baseline was ready — saying the most frequent label — here each lesson built its baseline itself. **This is the rule of the course: an unlabeled result exists only as much as it departs from a random result; a clustering whose baseline is not built counts as unmeasured.** Three rows cannot be read without a footnote: in Clustering 06 the 0.5235 is below its baseline, the +0.7351 of Clustering 04 comes not from the data but from the metric, the 0.5032 of Dimensionality reduction 03 is the best of six candidates and the honest range is 0.3940–0.5032. The second claim of the course is read from the same table: in an unsupervised result there is no right answer, **the choice of metric determines the answer**. The number of lessons in which the answer changed when the metric was changed is **ten** — three definitions of similarity gave three partitions, three linkage criteria chose 3, 37 and 56 clusters at the same height, three internal validity metrics chose 3, 6 and 10; in clustering 04 two metrics gave opposite rankings, in dimensionality reduction 03 what one preserved the other threw away, and here confidence and lift produced two separate lists. In twenty-one lessons — eleven of them unlabeled — every method worked on the assumption that the features were **given**. Scaling, encoding, generating and reducing were always decided **by hand**, and the cost of each decision was measured separately. What a procedure that also learns these decisions from the data itself gains, and what it makes invisible, was never asked. The next course, Deep Learning, opens with this question.