Dimensionality Reduction as a Feature Tool

The previous lesson measured one way of reducing the column count: keeping or dropping the columns one by one. The embedded method kept fourteen of the twenty-five candidates and gave 0.8254 on the held-out set; the thirteen-column set built by hand-measuring over two lessons gave the same number. Every one of the kept columns has a name, and what they mean is known.

The second way of reducing the column count is mixing the columns instead of dropping them. Thirteen columns can be turned into a few new columns, each a weighted sum of the thirteen; these are called components, and the components are ordered starting from the directions that carry the most variability in the set. The Statistical Analysis and Exploratory Data Analysis course measured the component share there, and the theory is not repeated here. The only question this lesson asks is this: what do components do as a feature — at how many components is the held-out number of the raw features caught up with, and what is lost in return.

• FS22. The raw representation is thirteen columns: the eleven columns of the Introduction to Machine Learning course, per_person and kw_seal. Its number on the held-out set is 0.8254 and it is the line before the decision of this lesson.
• FS23. The scale and the component directions are learned from the training set only and applied to all three sets with the same numbers. If they were learned from the whole set, the measurement would inflate.
• FS24. The components are written within the lesson: the covariance matrix of the standardized training columns is built, its largest direction is found by iteration, the direction found is subtracted from the matrix, and the operation is repeated thirteen times. No library is called.
• FS25. The number of iterations is 200 for every component; this number was not chosen by looking at the result. If the remaining magnitude drops below $10^{-12}$, the iteration stops early.
• FS26. The component share is the ratio of the variability the component carries to the total variability.
• FS27. The components do not see the label. Their ranking is by variability, not by discriminating power.
• FS28. The model does not change: a decision tree of depth 6. The baseline is 0.7579.
• FS29. The interpretability measure is which column the splits of the tree are read from; each split is weighted by the number of rows it splits.

Thirteen Columns, Thirteen Components

# COMMON — MODEL: the FICTIONAL subscriber table of the Introduction to Machine Learning course is rebuilt with the same seed;
# the model is a decision tree of depth 6. The only new column left from the previous lesson is per_person.
import math

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [("north", 0.28, 21, 0.00), ("south", 0.22, 17, -0.10),
          ("east", 0.18, 26, 0.30), ("west", 0.14, 14, -0.05),
          ("central", 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)
    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, vol, r = sum(v) / len(v), max(v) - min(v), rng(SEED + 51000 + k["id"])
    z = (0.052 * (mean - 20) + 0.85 * (min(v) == 0.0) + 0.026 * vol
         + 0.24 * (k["household"] >= 5) - 0.20 * (k["satisfaction"] >= 4) + k["region"][3]
         + 0.9 * (k["id"] - 10001) / 1399 + (r() + r() + r() - 1.5) * 1.30)
    x = {"id": k["id"], "household": k["household"], "satisfaction": k["satisfaction"],
         "period": len(v), "mean_consumption": round(mean, 2), "volatility": round(vol, 2),
         "zero_readings": int(min(v) == 0.0), "suspicious": int(z > 1.35),
         "region_name": k["region"][0]}
    for b in REGION:
        x["area_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)

BASE = ["mean_consumption", "volatility", "household", "satisfaction", "period",
        "zero_readings", "area_north", "area_south", "area_east", "area_west", "area_central"]
BASE_FIELDS = BASE + ["per_person"]
for x in DATA:
    x["per_person"] = round(x["mean_consumption"] / x["household"], 3)


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


def gini(rows):
    p = sum(x["suspicious"] for x in rows) / len(rows) if rows else 0.0
    return 2 * p * (1 - p)


def tree(rows, depth, fields, min_size=2):
    p, best = sum(x["suspicious"] for x in rows) / len(rows), None
    if depth and len(rows) >= 2 * min_size and 0.0 < p < 1.0:
        for f in fields:
            values = sorted({x[f] for x in rows})
            for v in (values[1:] if len(values) < 10 else
                      [values[int(i * len(values) / 10)] for i in range(1, 10)]):
                left = [x for x in rows if x[f] < v]
                right = [x for x in rows if x[f] >= v]
                if min(len(left), len(right)) < min_size:
                    continue
                k = gini(rows) - (len(left) * gini(left) + len(right) * gini(right)) / len(rows)
                if best is None or k > best[0]:
                    best = (k, f, v, left, right)
    if best is None or best[0] <= 1e-9:
        return {"prediction": int(p > 0.5)}
    return {"field": best[1], "threshold": best[2],
            "left": tree(best[3], depth - 1, fields, min_size),
            "right": tree(best[4], depth - 1, fields, min_size)}


def predict(node, x):
    while "prediction" not in node:
        node = node["left"] if x[node["field"]] < node["threshold"] else node["right"]
    return node["prediction"]


def accuracy(node, rows):
    return sum(predict(node, x) == x["suspicious"] for x in rows) / len(rows)


def baseline(train, rows):
    c = int(sum(x["suspicious"] for x in train) * 2 > len(train))
    return sum(x["suspicious"] == c for x in rows) / len(rows)


def measure(fields, depth=6):
    m = tree(TRAIN, depth, fields)
    return accuracy(m, TRAIN), accuracy(m, TEST)

Three steps are needed to find the components: bringing the columns to a common scale, building the covariance matrix, and finding the direction that carries the largest variability and subtracting it from the matrix. When the third step is repeated thirteen times, thirteen components come out. The scale and all the directions are read from the training set only.

for x in DATA:
    r = rng(SEED + 61000 + x["id"])
    r(), r()
    x["kw_seal"] = int(r() < (0.42 if x["suspicious"] else 0.10))
RAW_FIELDS = BASE_FIELDS + ["kw_seal"]

MEAN = {f: sum(x[f] for x in TRAIN) / len(TRAIN) for f in RAW_FIELDS}
SPREAD = {f: max(1e-9, (sum((x[f] - MEAN[f]) ** 2 for x in TRAIN) / len(TRAIN)) ** 0.5)
          for f in RAW_FIELDS}
for x in DATA:
    x["_z"] = [(x[f] - MEAN[f]) / SPREAD[f] for f in RAW_FIELDS]

n = len(RAW_FIELDS)
C = [[sum(x["_z"][i] * x["_z"][j] for x in TRAIN) / len(TRAIN) for j in range(n)]
     for i in range(n)]


def multiply(m, v):
    return [sum(m[i][j] * v[j] for j in range(len(v))) for i in range(len(m))]


COMPONENTS, SHARE = [], []
A = [row[:] for row in C]
for _ in range(n):
    v = [1.0 / (i + 1) for i in range(n)]
    for _ in range(200):
        w = multiply(A, v)
        u = sum(q * q for q in w) ** 0.5
        if u < 1e-12:
            break
        v = [q / u for q in w]
    lam = sum(v[i] * multiply(A, v)[i] for i in range(n))
    COMPONENTS.append(v)
    SHARE.append(max(0.0, lam))
    A = [[A[i][j] - lam * v[i] * v[j] for j in range(n)] for i in range(n)]

TOTAL = sum(SHARE)
for x in DATA:
    for i, v in enumerate(COMPONENTS):
        x["c" + str(i)] = round(sum(x["_z"][j] * v[j] for j in range(n)), 4)

RAW_RESULT = measure(RAW_FIELDS)
print(f"{'representation':<22} {'columns':>7} {'share':>7} {'train':>9} {'held_out':>9} {'diff':>8}")
print(f"{'raw features':<22} {n:>7} {'':>7} {RAW_RESULT[0]:>9.4f} {RAW_RESULT[1]:>9.4f}")
for k in (1, 2, 3, 5, 8, 11, 13):
    e, s = measure(["c" + str(i) for i in range(k)])
    print(f"{'component ' + str(k):<22} {k:>7} {sum(SHARE[:k]) / TOTAL:>7.3f} {e:>9.4f} "
          f"{s:>9.4f} {s - RAW_RESULT[1]:>+8.4f}")
print(f"\nbaseline {baseline(TRAIN, TEST):.4f}")

representation         columns   share     train  held_out     diff
raw features                13            0.8876    0.8254
component 1                  1   0.177    0.8003    0.7381  -0.0873
component 2                  2   0.323    0.8651    0.8175  -0.0079
component 3                  3   0.432    0.8677    0.8056  -0.0198
component 5                  5   0.630    0.9074    0.7778  -0.0476
component 8                  8   0.877    0.9074    0.7460  -0.0794
component 11                11   0.990    0.8995    0.7817  -0.0437
component 13                13   1.000    0.9048    0.7897  -0.0357

baseline 0.7579

The answer to the first half of the question is in the table: no number of components catches up with the raw features. The closest row is two components; with 0.8175 it is 0.0079 below the raw representation, that is, two subscribers. Thirteen columns have been reduced to two columns and two subscribers have been lost. This is a defensible balance as a compression decision.

The real finding of the table comes after this. As the number of components grows, the number does not improve: 0.8056 with three components, 0.7778 with five, 0.7460 with eight. The component share, on the other hand, grows steadily — from 0.323 to 0.877. The two columns do not go in the same direction. When all thirteen components are used there is no loss of information, the representation is a rotated version of the raw representation, and still the number stays at 0.7897; the difference is 0.0357. A decision tree makes its splits perpendicular to the axes, and building the same separation on a rotated axis takes more splits. The component representation is not independent of the model.

Share Ranking Is Not Separation Ranking

The components never see the label; their ranking looks only at which direction has more variability. Whether the label sits in those directions is a separate question and has to be measured.

def usage(node, rows):
    counts = {}

    def walk(node, rows):
        if "prediction" in node:
            return
        counts[node["field"]] = counts.get(node["field"], 0) + len(rows)
        walk(node["left"], [x for x in rows if x[node["field"]] < node["threshold"]])
        walk(node["right"], [x for x in rows if x[node["field"]] >= node["threshold"]])
    walk(node, rows)
    total = sum(counts.values())
    return ", ".join(f"{f} {counts[f] / total:.2f}"
                     for _, f in sorted((-counts[f], f) for f in counts)[:4])


print("most split in raw features:", usage(tree(TRAIN, 6, RAW_FIELDS), TRAIN))
print("in component features     :",
      usage(tree(TRAIN, 6, ["c" + str(i) for i in range(13)]), TRAIN))
for i in (0, 1):
    y = sorted(((-abs(COMPONENTS[i][j]), RAW_FIELDS[j], COMPONENTS[i][j]) for j in range(n)))
    print(f"c{i} loadings (share {SHARE[i] / TOTAL:.3f}): " +
          ", ".join(f"{f} {w:+.2f}" for _, f, w in y[:4]) +
          f", ... total weight of the remaining {n - 4} columns "
          f"{sum(abs(w) for _, _, w in y[4:]):.2f}")
print()
print(f"{'single component':<16} {'share':>7} {'train':>9} {'held_out':>9}")
for i in range(5):
    e, s = measure(["c" + str(i)])
    print(f"{'c' + str(i):<16} {SHARE[i] / TOTAL:>7.3f} {e:>9.4f} {s:>9.4f}")

most split in raw features: mean_consumption 0.35, volatility 0.19, kw_seal 0.17, per_person 0.14
in component features     : c0 0.29, c1 0.25, c8 0.12, c2 0.11
c0 loadings (share 0.177): per_person +0.54, mean_consumption +0.53, area_west -0.34, area_east +0.30, ... total weight of the remaining 9 columns 1.16
c1 loadings (share 0.146): volatility +0.62, zero_readings +0.51, period +0.42, kw_seal +0.29, ... total weight of the remaining 9 columns 0.86

single component   share     train  held_out
c0                 0.177    0.8003    0.7381
c1                 0.146    0.8267    0.7421
c2                 0.109    0.8188    0.7579
c3                 0.105    0.7937    0.7579
c4                 0.093    0.8148    0.7659

The last table shows with numbers that the share ranking is not the separation ranking. The component with the largest share, c0, gives 0.7381 on its own — below the baseline by 0.0198. c4, with about half of its share, gives 0.7659 and stays above the baseline. “The direction that carries the most variability” and “the direction that separates the most” are two different things, and since the components are ordered only by the first, taking the first k components is an arbitrary cut as far as separation is concerned.

It is also understood from here why two components together give 0.8175: c0 carries the consumption magnitude, while c1 carries the volatility and the zero readings. Two separate causes of the label have fallen onto two separate components, and the model can combine the two.

What Is Lost from Interpretability

The second half of the loss is not in the number but in the names. In the raw representation, 0.35 of the splits are read from mean_consumption, 0.19 from volatility and 0.17 from the seal indicator. These three rows can be given to a field crew: high consumption, volatile readings and a seal record. In the component representation the same list comes out as c0 0.29, c1 0.25, c8 0.12 and says nothing.

What a component is can be understood only by reading the loadings, and the loadings do not point to a single column. c0 carries mainly consumption per person and mean consumption, but the total weight of the remaining nine columns is 1.16 — that is, more than half of the component comes from columns outside the four large loadings. The answer to “why was this subscriber flagged as suspicious” is a column name and a threshold in the raw representation; in the component representation it is a weighted sum of thirteen numbers.

The decision is a comparison of two losses. Going down to two components buys 0.0079 of number and a saving of eleven columns of storage; in return the ability to write the justification of every decision the model gives into a table is lost. In a field such as illicit-use suspicion the second loss is usually more expensive, because the justification of an action taken against a subscriber can be asked for. The feature role of the component representation is therefore a narrow role: it is valuable if the column count is truly too large to manage or if the columns are copies of each other, and not when there are thirteen columns measured one by one.

Summary

• The thirteen raw columns give 0.8254 on the held-out set; no number of components catches up with this, the closest row being two components with 0.8175 (difference 0.0079, two subscribers).
• While the component share grows, the held-out number does not grow: 0.8056 with three components, 0.7778 with five, 0.7460 with eight; with thirteen components, with no loss of information, the number still stays at 0.7897.
• The component with the largest share falls below the baseline on its own with 0.7381, while the fifth component with a share of 0.093 gives 0.7659; share ranking is not separation ranking.
• Interpretability is a measurable loss: in the raw representation 0.35 of the splits are read from the mean consumption, in the component representation the same share falls on a mixture with no name.
• The scale and the component directions are learned from the training set only; if they were learned from the whole set, the measurement would inflate.

Next Step

In four lessons the number of preparation steps grew, and in each of them the same sentence was repeated: this step is learned from the training set only. The scale mean, the region median, the selected column set, the component directions — all four are trained objects, and all four have to be applied by hand, separately, in the right order. When this chain is applied by hand, one has to count how many steps were learned from the wrong set. The next lesson combines the preparation and the model in a single object and measures how many steps leak without a pipeline.