Singular Value Decomposition Application

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.

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.

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.

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.