Polynomial Regression

The previous lesson built the model as a plane and used the columns as given: household size one column, satisfaction one column, period count one column. What that model says is that when household size rises by one, consumption rises by the same amount everywhere. There is no reason for this to be true. The consumption that grows in going from a two-person household to a four-person one may not grow at the same rate in going from four to six; the effect of the region may change with household size.

Building this relation does not require changing the model family. A linear model is linear in its coefficients, not in its columns. When the square of household size is added as a column, the model is still solved by least squares, but now it carries curvature. Polynomial features is the name of this expansion, and it has a single hyperparameter: the degree. This lesson sweeps the degree from 1 to 8, prints the train and held-out numbers of each degree side by side, and shows which set the choice was read from as two separate numbers.

• LM10. The set, the split and the target are as in the previous lesson: the same fictional subscriber table, seed 20260218, 756/252/252, the target is mean consumption (m³), the baseline is the train mean and a root error of 6.5324 m³ on the test set.
• LM11. The polynomial is expanded over three numeric columns only: household size, satisfaction, period count. For degree $d$, all the monomials that satisfy $1 \le i+j+k \le d$ are added.
• LM12. The four region indicators stay as they are at every degree; they do not enter the polynomial. North is the reference as in the previous lesson.
• LM13. The columns are scaled with the mean and spread learned from the train set. The eighth power of household size exceeds one million; without scaling the elimination collapses numerically.
• LM14. At high degrees the columns turn into linear combinations of one another, because household size takes only six discrete values, satisfaction five and period count three. To make the system of equations solvable, $10^{-6}$ is added to the diagonal. This is not a regularization decision but a numerical threshold set against a singular system; the penalty decision is the subject of the next lesson.
• LM15. The number of candidates tried is 8: eight degrees, no other tuning.
• LM16. The choice is read from the validation set. The test set numbers appear in the table but do not enter the choice; the reason they appear is to measure the difference between the two choices.
• LM17. The ceiling model is the procedure that gives every discrete combination of the same four features the train set mean. No model built with these features can be better than it on the train set.

As the Degree Grows, Two Numbers Diverge

The setup is the setup of the previous lesson; polynomial features and the degree sweep are added on top of it.

# polynomial.py — MODEL: the same FICTIONAL subscriber table, the same seed and the same split.
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
    x = {"region_name": k["region"][0], "household": k["household"],
         "satisfaction": k["satisfaction"], "period": len(v),
         "target": round(sum(v) / len(v), 2)}
    for b in REGION:
        x["area_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)


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:]
MEAN = sum(x["target"] for x in TRAIN) / len(TRAIN)


def rmse(predict, s):
    return math.sqrt(sum((predict(x) - x["target"]) ** 2 for x in s) / len(s))


print(f"train {len(TRAIN)}, validation {len(VAL)}, test {len(TEST)}, "
      f"baseline on test {rmse(lambda x: MEAN, TEST):.4f} m3")

train 756, validation 252, test 252, baseline on test 6.5324 m3

Polynomial features, scaling and the solution stand together. The scaling step follows the rule of the previous course and has two faces: fitting only from the train set, applying to every set.

import itertools

NUMERIC = ["household", "satisfaction", "period"]
INDICATORS = ["area_south", "area_east", "area_west", "area_central"]


def powers(degree):                         # all monomials with 1 <= i+j+k <= degree
    return [u for u in itertools.product(range(degree + 1), repeat=len(NUMERIC))
            if 1 <= sum(u) <= degree]


def raw(x, pw):
    return ([float(x[a]) for a in INDICATORS]
            + [math.prod(x[a] ** k for a, k in zip(NUMERIC, u)) for u in pw])


def fit_scale(train, pw):                   # LEARN: from the train set only
    H = [raw(x, pw) for x in train]
    p = len(H[0])
    o = [sum(r[j] for r in H) / len(H) for j in range(p)]
    s = [max(1e-9, math.sqrt(sum((r[j] - o[j]) ** 2 for r in H) / len(H)))
         for j in range(p)]
    return o, s


def apply(x, pw, o, s):                     # APPLY
    h = raw(x, pw)
    return [1.0] + [(h[j] - o[j]) / s[j] for j in range(len(h))]


def solve(A, b):
    n = len(b)
    for i in range(n):
        p = max(range(i, n), key=lambda k: abs(A[k][i]))
        A[i], A[p], b[i], b[p] = A[p], A[i], b[p], b[i]
        for k in range(i + 1, n):
            f = A[k][i] / A[i][i]
            for j in range(i, n):
                A[k][j] -= f * A[i][j]
            b[k] -= f * b[i]
    w = [0.0] * n
    for i in range(n - 1, -1, -1):
        w[i] = (b[i] - sum(A[i][j] * w[j] for j in range(i + 1, n))) / A[i][i]
    return w


EPSILON = 1e-6                              # numerical threshold that makes the singular system solvable


def polynomial(train, pw, o, s):            # MODEL: polynomial regression
    X = [apply(x, pw, o, s) for x in train]
    y = [x["target"] for x in train]
    p = len(X[0])
    A = [[sum(X[i][a] * X[i][b] for i in range(len(X)))
          + (EPSILON if a == b and a else 0.0) for b in range(p)] for a in range(p)]
    return solve(A, [sum(X[i][a] * y[i] for i in range(len(X))) for a in range(p)])


def predictor(w, pw, o, s):
    return lambda x: sum(a * b for a, b in zip(w, apply(x, pw, o, s)))


RESULTS = []
print(f"{'degree':>7} {'columns':>7} {'train':>9} {'validation':>10} {'test':>8} "
      f"{'weight':>10}")
for degree in range(1, 9):
    pw = powers(degree)
    o, s = fit_scale(TRAIN, pw)
    w = polynomial(TRAIN, pw, o, s)
    g = predictor(w, pw, o, s)
    RESULTS.append((degree, len(pw) + len(INDICATORS), rmse(g, TRAIN), rmse(g, VAL),
                    rmse(g, TEST)))
    print(f"{degree:>7} {RESULTS[-1][1]:>7} {RESULTS[-1][2]:>9.4f} {RESULTS[-1][3]:>10.4f} "
          f"{RESULTS[-1][4]:>8.4f} {math.sqrt(sum(v * v for v in w[1:])):>10.2f}")
print(f"\nbaseline on test {rmse(lambda x: MEAN, TEST):.4f} m3, candidates tried 8")

 degree columns     train validation     test     weight
      1       7    5.0757     5.3225   5.3096       3.26
      2      13    5.0531     5.3219   5.2648       5.12
      3      23    5.0355     5.3184   5.3183      17.22
      4      38    5.0013     5.2569   5.3775     212.08
      5      59    4.9228     5.2635   5.2009    1520.80
      6      87    4.8933     5.3610   5.3224    3487.91
      7     123    4.8733     5.4923   5.5168    6959.24
      8     168    4.8437     5.6877   5.9085    9344.45

baseline on test 6.5324 m3, candidates tried 8

The train column is one-directional and improves without interruption: from 5.0757 to 4.8437. This is not a finding but an identity. The columns of degree $d$ contain the columns of degree $d-1$, so the solution of the smaller degree is also reachable at the larger degree; least squares can only improve on it. It is impossible for the train number to get worse along a growing column set, and for that reason the train number can never be used in choosing the degree.

The validation column is two-directional: it falls from 5.3225 to 5.2569, then climbs to 5.6877. The weight column shows the mechanism of the climb. The coefficient magnitude rises from 3.26 to 9,344.45, that is, about three thousand times. High-degree monomials are nearly linear combinations of one another — the sixth power of a column that takes six discrete values is a combination of its lower powers — and the model gives these nearly identical columns huge coefficients that cancel one another. On the train set this is free; on the held-out set those coefficients are not cancelled exactly and the difference turns into error.

Which Set the Choice Was Read From

One degree has to be chosen from the table of eight rows, and this choice can be made from two different sets.

best_val = min(RESULTS, key=lambda t: t[3])
best_test = min(RESULTS, key=lambda t: t[4])
first = RESULTS[0]
print(f"degree chosen on validation {best_val[0]}: validation "
      f"{best_val[3]:.4f}, test {best_val[4]:.4f}")
print(f"degree if chosen on test {best_test[0]}: validation "
      f"{best_test[3]:.4f}, test {best_test[4]:.4f}")
print(f"difference between the two choices {best_val[4] - best_test[4]:.4f} m3")
print(f"degree 1 on test {first[4]:.4f}; degree chosen on validation "
      f"{best_val[4] - first[4]:+.4f} m3 behind it\n")
print(f"{'degree':>7} {'divergence (validation - train)':>31}")
for d, cols, e, v, t in RESULTS:
    print(f"{d:>7} {v - e:>31.4f}")

degree chosen on validation 4: validation 5.2569, test 5.3775
degree if chosen on test 5: validation 5.2635, test 5.2009
difference between the two choices 0.1766 m3
degree 1 on test 5.3096; degree chosen on validation +0.0679 m3 behind it

 degree divergence (validation - train)
      1                          0.2467
      2                          0.2687
      3                          0.2828
      4                          0.2556
      5                          0.3407
      6                          0.4677
      7                          0.6189
      8                          0.8440

The first three lines show in a single example why the rule of the course is necessary. The validation set chooses degree 4; that degree gives 5.3775 on test. Had the choice been made from the test set, degree 5 would have come and 5.2009 would have been reported. The 0.1766 m³ between the two numbers is the price of the lesson’s eight-candidate budget and comes from nothing the model learned.

The fourth line is more unsettling. The degree 4 chosen from validation is 0.0679 m³ behind degree 1 on test. That is, trying eight candidates and choosing by looking at the validation set produced a worse result than not trying and staying with the simplest model. This is not a calculation error; it is the consequence of the validation set also being 252 rows; the smallest number on that set includes that set’s own scatter. The second claim of the course is measured for the first time on this line: a stronger family did not beat the weaker one.

The divergence column reduces the table to a single number. The difference between validation and train stays between 0.25 and 0.28 up to degree 4, rises to 0.3407 at degree 5 and becomes 0.8440 at degree 8. The place where the divergence begins is degree 5; one step after the place where the validation number bottoms out.

The Ceiling of These Features

At degree 8 there are 168 columns. All of these columns are functions of four features, and those four features produce only a finite number of combinations.

cell = {}
for x in TRAIN:
    cell.setdefault((x["region_name"], x["household"], x["satisfaction"], x["period"]),
                    []).append(x["target"])
CELL_MEANS = {k: sum(v) / len(v) for k, v in cell.items()}


def cell_prediction(x):                     # every cell gets the train mean: the ceiling model
    return CELL_MEANS.get((x["region_name"], x["household"], x["satisfaction"], x["period"]), MEAN)


print(f"discrete feature combinations 5*6*5*3 = {5 * 6 * 5 * 3}, "
      f"filled cells in train {len(cell)}")
print(f"cell mean model: train {rmse(cell_prediction, TRAIN):.4f}, "
      f"test {rmse(cell_prediction, TEST):.4f}")
print(f"cells seen in test but absent from train "
      f"{sum(1 for x in TEST if (x['region_name'], x['household'], x['satisfaction'], x['period']) not in CELL_MEANS)}")

discrete feature combinations 5*6*5*3 = 450, filled cells in train 268
cell mean model: train 3.8418, test 6.3887
cells seen in test but absent from train 33

The cell mean model is the most flexible procedure that can be built with these four features: it gives every combination its own train mean and carries no shape assumption. On the train set it gives 3.8418, that is, it comfortably beats the 4.8437 of the degree 8 polynomial. On the test set it gives 6.3887 and comes closer to the baseline’s 6.5324 by only 0.1437 m³.

This chasm between the two numbers says where the degree sweep is heading. 756 rows fall into 268 filled cells, that is, three subscribers per cell on average; the mean of three subscribers is not a prediction, it is those three subscribers themselves. As the polynomial degree grows, the model approaches this procedure and the test number goes toward its number. The last line shows the same limit from the other direction: of the 252 subscribers in the test set, 33 fall into a cell never seen in the train set, and the only thing that can be said to them is the baseline.

Summary

• Polynomial features do not change the model family; a linear model is linear in its coefficients, not in its columns. Going from degree 1 to 8 the number of columns rises from 7 to 168.
• The train root error falls without interruption from 5.0757 to 4.8437 m³ and this is an identity: the columns of the larger degree contain those of the smaller. The train number cannot be used in choosing the degree.
• The validation root error bottoms out at 5.2569 and climbs to 5.6877; the weight magnitude rises over the same range from 3.26 to 9,344.45. The mechanism of the climb is the huge coefficients that cancel one another, given to nearly identical columns.
• The degree read from validation, 4, gives 5.3775 on test; had it been read from test, degree 5 would have reported 5.2009 and the difference is 0.1766 m³. The same choice stays 0.0679 m³ behind degree 1 on test.
• Four features produce 450 discrete combinations and 268 of them are filled in train; the ceiling model that gives every cell its own mean yields 3.8418 on train and 6.3887 on test. The 33 subscribers in the test set fall into a cell never seen in training.

Next Step

The degree sweep turned a knob in whole-number steps: there is no intermediate stop between degree 4 and 5, and moving to degree 5 adds twenty-one columns at once. Between adding columns and deleting them there is a third way — keeping them all and penalizing their coefficients. The Introduction to Machine Learning course had built this intuition with a single penalty type and had measured that the two numbers go in opposite directions as the penalty coefficient grows. The next lesson changes the form of the penalty: the ridge, lasso and elastic net penalties are applied separately to the same column set and the measure is no longer the error but sparsity — how many coefficients each penalty sets exactly to zero, which one shrinks them all and leaves them, and what the zeroed columns are worth on the held-out set.