--- title: 'Singular Value Decomposition Application' source: 'https://academia.sh/en/courses/unsupervised-learning/singular-value-decomposition-application' course: 'Unsupervised Learning' language: en updated: '2026-10-21T05:03:32+00:00' license: 'CC BY-SA 4.0' objectives: - 'Writes the singular value decomposition on a small matrix and verifies that the reconstruction error can be read from the discarded singular values' - 'Prints the trade-off between rank and the storage cost (in numbers): as the rank falls the error rises, the storage cost falls, and past some point the compression gain ends' - 'Measures the low-rank approximation against the baseline of storing raw columns with the same budget' --- # 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.