Regularization Intuition

The previous two lessons restricted capacity with a crude knob: the depth of the tree is an integer, there is no intermediate value between 3 and 5, and the same goes for the leaf lower bound. Yet the complexity of a model can be penalized instead of being cut off. Next to the error measure, a term that measures how assertive the model is is added, and the model is forced into a compromise between lowering the error and lowering its assertiveness. The penalty coefficient is a continuous number; between two integers there are as many stops as one likes.

This lesson raises that coefficient from 0.001 to 1000 and prints three things at every step: the magnitude of the weights, the training error and the held-out error. The expected behavior is not one-directional — the training error worsens without interruption, the held-out error first improves and then worsens — and from which set the best point is to be read is the last question of the lesson.

• GE23. The model is linear and has 29 columns: five region indicators, household size, satisfaction, number of periods, their squares, their pairwise products and the product of each with the region.
• GE24. The penalty is a quadratic term: $\lambda$ times the sum of squares of the weights is added to the error measure. Instead of a constant term, the training mean is subtracted from the target, so that all the weights are penalized and, as the penalty grows, the model converges to the baseline.
• GE25. The columns are scaled by the mean and the scatter of the training set. The penalty is sensitive to scale and from which set the scale is computed is a decision.
• GE26. The training set is deliberately cut down to 120 rows. 120 rows for twenty-nine columns is the ratio that makes overfitting measurable.
• GE27. The validation and test sets are the 252-row sets of the previous lessons; the seed and the split are the same.
• GE28. The measure is the root mean squared error (m³); the decomposition is made in squared units.
• GE29. The baseline is the model in which all the weights are zero: telling every subscriber the mean of the training set.
• GE30. The test set is looked at once in this lesson, after the selection is finished.

The Penalty Moves the Two Numbers in Opposite Directions

The setup generates the set with the same seed, builds the linear model and raises the penalty coefficient in eleven steps. The solution is in closed form: the weights that minimize the squared error are the solution of a system of equations with the penalty added to its diagonal.

# regularization.py — MODEL. The same regional measurement network (a fictional setup), the same seed. The target is
# the subscriber's mean consumption; the model is linear, the penalty quadratic.
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 = []
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})
DATA = []
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 v:
        DATA.append({"id": k["id"], "region": k["region"][0], "household": k["household"],
                     "satisfaction": k["satisfaction"], "period": len(v),
                     "target": round(sum(v) / len(v), 2)})


def columns(x):
    t = [float(x["region"] == b[0]) for b in REGION]
    h, m, d = float(x["household"]), float(x["satisfaction"]), float(x["period"])
    return (t + [h, m, d, h * h, m * m, d * d, h * m, h * d, m * d]
            + [g * h for g in t] + [g * m for g in t] + [g * d for g in t])


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[:120], K[756:1008], K[1008:]
P = len(columns(TRAIN[0]))
RAW_TRAIN = [columns(x) for x in TRAIN]
MEANS = [sum(s[j] for s in RAW_TRAIN) / len(RAW_TRAIN) for j in range(P)]
SPREADS = [max(1e-9, math.sqrt(sum((s[j] - MEANS[j]) ** 2 for s in RAW_TRAIN) / len(RAW_TRAIN)))
           for j in range(P)]
MEAN = sum(x["target"] for x in TRAIN) / len(TRAIN)


def scale(x):
    s = columns(x)
    return [(s[j] - MEANS[j]) / SPREADS[j] for j in range(P)]


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


def fit(rows, penalty):
    X = [scale(x) for x in rows]
    y = [x["target"] - MEAN for x in rows]
    A = [[sum(X[i][a] * X[i][b] for i in range(len(X))) + (penalty if a == b 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 rmse(w, s):
    return math.sqrt(sum((MEAN + sum(a * b for a, b in zip(w, scale(x))) - x["target"]) ** 2
                         for x in s) / len(s))


print(f"train {len(TRAIN)}, validation {len(VAL)}, test {len(TEST)}, columns {P}")
print(f"baseline: train {rmse([0.0] * P, TRAIN):.3f}  validation {rmse([0.0] * P, VAL):.3f}  "
      f"test {rmse([0.0] * P, TEST):.3f}")
print(f"\n{'penalty':>8} {'weight':>9} {'train':>9} {'validation':>10} {'test':>8}")
for c in (0.001, 0.01, 0.1, 0.3, 1, 3, 10, 30, 100, 300, 1000):
    w = fit(TRAIN, c)
    print(f"{c:>8} {math.sqrt(sum(v * v for v in w)):>9.3f} {rmse(w, TRAIN):>9.3f} "
          f"{rmse(w, VAL):>10.3f} {rmse(w, TEST):>8.3f}")

train 120, validation 252, test 252, columns 29
baseline: train 6.121  validation 6.439  test 6.556

 penalty    weight     train validation     test
   0.001    10.040     4.603      5.924    5.841
    0.01     9.963     4.603      5.921    5.839
     0.1     9.334     4.604      5.894    5.815
     0.3     8.396     4.606      5.850    5.777
       1     6.687     4.620      5.770    5.704
       3     4.891     4.653      5.693    5.634
      10     3.291     4.716      5.614    5.571
      30     2.222     4.804      5.556    5.534
     100     1.534     4.920      5.553    5.559
     300     1.099     5.089      5.657    5.694
    1000     0.601     5.453      5.932    6.008

The weight column shows directly what the penalty does: the magnitude of the weights falls from 10.040 to 0.601, that is, to one seventeenth. The penalty does nothing but shrink the weights of the model; it does not touch the columns, it does not drop rows, it does not delete rules.

The third column is one-directional and it is worsening: the training error rises from 4.603 to 5.453 m³. This is expected, because the unpenalized solution is already the solution that makes the squared error on the training set smallest; every penalty put on top of it can only spoil that number. The worsening of the training error is not a side effect of regularization but the proof that it is working.

The fourth column is two-directional: it falls from 5.924 to 5.553, then returns to 5.932. The lowest point is at penalty 100 and gives a gain of 0.371 m³ over the unpenalized model. The distance between the fall and the rise is not narrow: every value of the penalty between 3 and 300 is better than the unpenalized model, but at 1000 the model is worse than even the unpenalized model.

The direction of the last row is the intuition of the lesson. As the penalty grows the weights approach zero, and when all the weights are zero the model tells every subscriber the training mean — that is, it turns into the baseline itself. Regularization is a force that pulls the model toward the baseline with a continuous knob. At penalty 1000 the weight, 0.601, is still not zero, and the validation error, 5.932, is below the baseline’s 6.439; that knob has not yet been turned all the way.

Looking at the weights themselves shows what the penalty drives out.

NAMES = ([b[0] for b in REGION] + ["household", "satisfaction", "period", "household^2",
         "satisfaction^2", "period^2", "household*satisfaction", "household*period",
         "satisfaction*period"] + [b[0] + "*household" for b in REGION]
         + [b[0] + "*satisfaction" for b in REGION] + [b[0] + "*period" for b in REGION])
for c in (0.001, 100):
    w = fit(TRAIN, c)
    top = sorted(range(P), key=lambda j: -abs(w[j]))[:4]
    print(f"penalty {c}: four largest weights")
    for j in top:
        print(f"    {NAMES[j]:<20} {w[j]:>8.3f}")
    t = sum(abs(w[j]) for j in range(P))
    print(f"    share of the top four weights in the total weight: "
          f"{sum(abs(w[j]) for j in top) / t:.3f}")

penalty 0.001: four largest weights
    south                  -4.138
    period                  2.981
    household^2             2.965
    household              -2.962
    share of the top four weights in the total weight: 0.303
penalty 100: four largest weights
    west*period            -0.574
    east*period             0.553
    west                   -0.487
    east                    0.479
    share of the top four weights in the total weight: 0.302

The weights of household and household^2 in the unpenalized solution are -2.962 and +2.965: they almost cancel each other out. Because the two columns are two forms of the same information, the model can give both large and opposite weights; the predictions do not change, the weights swell. The penalty forbids exactly this. At penalty 100 the four largest weights are the west and east regions and the terms of these tied to the number of periods — this is the real structure in the fiction, because the region determines the base consumption.

The last line says how the penalty works. The share of the first four weights in the total weight is the same at both penalties: 0.303 and 0.302. The quadratic penalty does not eliminate columns, it does not zero out any of them, it does not change the shape of the distribution; it shrinks them all together. This is the difference from restricting capacity — tree pruning deletes a rule, the penalty softens all the rules by the same proportion.

Which Component of the Error the Penalty Shrinks

The decomposition of the previous lesson can be applied as is. If the penalty is a knob that restricts capacity, it must lower the variance component of the error and raise the bias component.

B = 30
SAMPLES = []
for b in range(B):
    r = rng(SEED + 61000 + b * 13)
    SAMPLES.append([TRAIN[int(r() * len(TRAIN))] for _ in range(len(TRAIN))])

print(f"{B} samples drawn with replacement, {len(TRAIN)} rows each")
print(f"\n{'penalty':>8} {'total':>9} {'variance':>9} {'bias':>9}")
for c in (0.001, 1, 30, 100, 1000):
    PRED = []
    for s in SAMPLES:
        w = fit(s, c)
        PRED.append([MEAN + sum(a * b for a, b in zip(w, scale(x))) for x in VAL])
    total = variance = bias = 0.0
    for i, x in enumerate(VAL):
        p = [PRED[b][i] for b in range(B)]
        m = sum(p) / B
        total += sum((q - x["target"]) ** 2 for q in p) / B
        variance += sum((q - m) ** 2 for q in p) / B
        bias += (m - x["target"]) ** 2
    n = len(VAL)
    print(f"{c:>8} {total / n:>9.3f} {variance / n:>9.3f} {bias / n:>9.3f}")

30 samples drawn with replacement, 120 rows each

 penalty     total  variance      bias
   0.001    47.027    12.988    34.039
       1    38.330     5.133    33.196
      30    32.689     1.622    31.066
     100    32.075     1.098    30.976
    1000    35.476     0.181    35.295

The variance column is the real target of the penalty: it falls from 12.988 to 0.181, that is, to one seventy-second. As the penalty coefficient grows, which 120 rows are drawn affects the answer the model gives less. At penalty 1000 the thirty models say almost the same thing, because all of them stand close to the baseline.

The bias column splits into two regions. The 34.039 at the unpenalized end is higher than expected: a system of equations with 29 columns on 120 rows is nearly singular and the solution swings so much that even the mean of the thirty models is corrupted. Up to penalty 100 this column falls to 30.976. Then it jumps to 35.295: there the penalty erases not only the noise but also the real region structure in the set.

The point where the total column is smallest is penalty 100, and this is the same as the point where the validation column in the previous table was smallest. Two independent measurements — one the root error on a single split, the other the component decomposition on thirty samples — point to the same coefficient.

From Which Set the Best Point Is Read

Penalty 100 was read from the table, but from which table it was read is the question of the first lesson of the course.

PENALTY = [0.001, 0.01, 0.1, 0.3, 1, 3, 10, 30, 100, 300, 1000]
SCORES = [(c, rmse(fit(TRAIN, c), VAL), rmse(fit(TRAIN, c), TEST)) for c in PENALTY]
by_val = min(SCORES, key=lambda t: t[1])
by_test = min(SCORES, key=lambda t: t[2])
print(f"best penalty if chosen on validation {by_val[0]:>3}: validation {by_val[1]:.3f}, "
      f"test {by_val[2]:.3f}")
print(f"best penalty if chosen on test       {by_test[0]:>3}: validation {by_test[1]:.3f}, "
      f"test {by_test[2]:.3f}")
print(f"if chosen on test, reported number {by_test[2]:.3f}, "
      f"honest {by_val[2]:.3f}, difference {by_val[2] - by_test[2]:.3f} m3")
print(f"\non the test set: baseline {rmse([0.0] * P, TEST):.3f}, unpenalized "
      f"{rmse(fit(TRAIN, 0.001), TEST):.3f}, chosen {by_val[2]:.3f}")
print(f"gain from the penalty {rmse(fit(TRAIN, 0.001), TEST) - by_val[2]:.3f} m3, "
      f"gain of the model over the baseline {rmse([0.0] * P, TEST) - by_val[2]:.3f} m3")

best penalty if chosen on validation 100: validation 5.553, test 5.559
best penalty if chosen on test        30: validation 5.556, test 5.534
if chosen on test, reported number 5.534, honest 5.559, difference 0.025 m3

on the test set: baseline 6.556, unpenalized 5.841, chosen 5.559
gain from the penalty 0.282 m3, gain of the model over the baseline 0.997 m3

The two sets choose two different coefficients: validation 100, test 30. That both are defensible is the heart of the matter — the difference between them is not a real difference in quality but the own scatter of the two 252-row sets. But if the selection had been made on the test set, the reported number would have been 5.534; the honest number is 5.559. The difference is small at 0.025 m³ and its sign is fixed: the number on the set where the selection is made can never come out worse. Here there were eleven candidates; in the first lesson 48 candidates had raised the difference to 2 points. The share of selection grows with the number of candidates, but it is not zero even with zero candidates.

The last two lines close the lesson’s account. On the test set the baseline is 6.556 m³, the unpenalized linear model 5.841 m³, the regularized model 5.559 m³. The gain of the model over the baseline is 0.997 m³, that is, 15.2%. Of this, 0.282 m³, that is, 28.3%, comes from the penalty coefficient alone. Regularization did not add a new feature or teach a new rule; it chose how much of what had been learned would be kept, and took more than a quarter of the gain from there.

Summary

• When the penalty coefficient is raised from 0.001 to 1000, the weight magnitude falls from 10.040 to 0.601; the training error worsens without interruption from 4.603 to 5.453 m³, and this is the proof that the penalty is working.
• The held-out error falls from 5.924 to 5.553 and returns to 5.932; the lowest point is at penalty 100 and, had the penalty been taken to infinity, the model would have turned into the baseline itself.
• The decomposition shows the target of the penalty: the variance falls from 12.988 to 0.181, while the bias hits a floor at 30.976 and jumps to 35.295 at penalty 1000; there the penalty erases the real region structure.
• The best penalty read from the validation set is 100, and it would have been 30 if read from the test set; choosing from the test would have lowered the reported number from 5.559 to 5.534, and the sign of this difference never turns around.
• On the test set the baseline is 6.556, the unpenalized model 5.841, the regularized model 5.559 m³; the gain of the model over the baseline is 0.997 m³ and 28.3% of it comes from the penalty coefficient alone.

Next Step

In four lessons the separation has been established and the baseline has been written next to every number. In the setup of this lesson there is a decision that was silently made correctly: the scale of the columns was computed from the training set only. The decision can be reversed, and when it is reversed the held-out number improves — and the model has learned nothing. The next lesson applies four such corruptions one by one and prints how much each of them inflates the held-out number. The close of the course is there too.