--- title: 'Principal Component Analysis' source: 'https://academia.sh/en/courses/unsupervised-learning/principal-component-analysis' course: 'Unsupervised Learning' language: en updated: '2026-10-21T05:03:32+00:00' license: 'CC BY-SA 4.0' objectives: - 'Measures the reconstruction error against the baseline of reconstruction with column means and against random orthogonal directions' - 'Shows that the component share and the preserved geometry do not say the same thing: with a share of 0.341 the mean of the pairwise distances falls to about half' - 'Compares how the scaling decision changes the component directions and the reconstruction error measured in a common unit' --- # 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.