Principal Component Analysis

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.

# 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:]

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.

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.

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.