Regularization Penalties

The previous lesson turned a knob in whole-number steps. There was no intermediate stop between degree 4 and 5, and moving to degree 5 added twenty-one columns at once. Besides adding or deleting columns there is a third way: keeping the columns and penalizing their coefficients. The Generalization topic of the Introduction to Machine Learning course had built this intuition with a single penalty type, had measured that as the penalty coefficient grows the train number worsens without interruption while the held-out number first improves and then worsens, and had shown that the best point is read from the validation set. That sweep is not repeated here.

The question of this lesson is not the size of the penalty but its form. Penalizing the sum of squares of the coefficients and penalizing the sum of absolute values produce two different procedures, and the difference shows not in the error but in how many coefficients are exactly zero. This is the measure: sparsity. How many columns each penalty drops from the model, which one shrinks them all and leaves them, and what the dropped columns are worth on the held-out set.

• LM18. The set, the split, the target and the baseline are as in the previous two lessons: the fictional subscriber table, seed 20260218, 756/252/252, a root error of 6.5324 m³ on the test set.
• LM19. The column set is the degree 6 polynomial: the four region indicators and all the monomials of the three numeric columns whose total degree does not exceed six, 87 columns in all. The degree does not change in this lesson; the only thing that changes is the penalty.
• LM20. The columns are scaled with the mean and spread learned from the train set. The penalty is sensitive to scale: on unscaled columns the same penalty coefficient is applied to different columns with different force.
• LM21. The intercept is not penalized. Had it been penalized, the model would say zero when all the coefficients were zeroed; when it is not penalized, it returns to the baseline, that is, to the train mean.
• LM22. Three penalty forms: the ridge penalty penalizes the sum of squares of the coefficients, the lasso penalty the sum of absolute values, the elastic net penalty an even mixture of the two. No product or figurative name is used; the penalty is referred to by its form.
• LM23. The solution is coordinate descent: each coefficient is updated in turn, with the others held fixed. Because the derivative of the absolute penalty is undefined at zero, there is no closed form; the descent sets a coefficient that stays below the threshold exactly to zero.
• LM24. The number of rounds is 120 and is fixed in all settings. The number of rounds is not swept as a hyperparameter; had it been swept, the number of candidates would have tripled.
• LM25. The “exactly zero” comparison is not made with ==; the threshold is $10^{-9}$, and since the descent assigns exactly 0.0 to the coefficient it zeroes, this threshold is only a safeguard.
• LM26. The number of candidates tried is 18: three penalty forms, six penalty coefficients.

Three Penalty Forms

The setup is the setup of the previous two lessons.

# penalty.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

Coordinate descent gathers the three penalty forms in a single procedure. For each coefficient the best value is computed with the others held fixed; the share of the lasso penalty pushes this value toward zero, and a coefficient that cannot pass the threshold lands exactly on zero.

import itertools

NUMERIC = ["household", "satisfaction", "period"]
INDICATORS = ["area_south", "area_east", "area_west", "area_central"]
DEGREE = 6
POWERS = [u for u in itertools.product(range(DEGREE + 1), repeat=3) if 1 <= sum(u) <= DEGREE]
# h household, s satisfaction, p period; h^2s^1 = household squared times satisfaction
NAMES = INDICATORS + ["".join(f"{a}^{k}" for a, k in zip("hsp", u) if k) for u in POWERS]


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


H = [raw(x) for x in TRAIN]
P = len(NAMES)
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)]


def apply(x):
    h = raw(x)
    return [(h[j] - O[j]) / S[j] for j in range(P)]


COLS = [[(H[i][j] - O[j]) / S[j] for i in range(len(H))] for j in range(P)]
SQ = [sum(v * v for v in COLS[j]) for j in range(P)]
Y = [x["target"] for x in TRAIN]


def descent(penalty, absolute_share, rounds=120):        # MODEL: penalized linear regression
    w, intercept = [0.0] * P, sum(Y) / len(Y)
    resid = [y - intercept for y in Y]
    for _ in range(rounds):
        for j in range(P):
            if SQ[j] < 1e-12:
                continue
            r = sum(COLS[j][i] * resid[i] for i in range(len(resid))) + SQ[j] * w[j]
            c1, c2 = penalty * absolute_share, penalty * (1 - absolute_share)
            if r > c1:
                updated = (r - c1) / (SQ[j] + c2)
            elif r < -c1:
                updated = (r + c1) / (SQ[j] + c2)
            else:
                updated = 0.0
            delta = updated - w[j]
            if delta:
                for i in range(len(resid)):
                    resid[i] -= COLS[j][i] * delta
                w[j] = updated
        d = sum(resid) / len(resid)
        intercept += d
        for i in range(len(resid)):
            resid[i] -= d
    return intercept, w


def predictor(m):
    return lambda x: m[0] + sum(a * b for a, b in zip(m[1], apply(x)))


PENALTIES = (3, 10, 30, 100, 300, 1000)
FORMS = (("ridge", 0.0), ("lasso", 1.0), ("elastic net", 0.5))
RESULTS = []
print(f"columns {P}, baseline on test {rmse(lambda x: MEAN, TEST):.4f} m3")
print(f"\n{'penalty form':<12} {'penalty':>8} {'exact zero':>10} {'train':>9} "
      f"{'validation':>10} {'test':>8}")
for name, share in FORMS:
    for c in PENALTIES:
        m = descent(c, share)
        g = predictor(m)
        zeros = sum(1 for v in m[1] if abs(v) < 1e-9)
        RESULTS.append((name, c, zeros, rmse(g, VAL), rmse(g, TEST)))
        print(f"{name:<12} {c:>8} {zeros:>10} {rmse(g, TRAIN):>9.4f} "
              f"{rmse(g, VAL):>10.4f} {rmse(g, TEST):>8.4f}")
print(f"\ncandidates tried {len(RESULTS)}")

columns 87, baseline on test 6.5324 m3

penalty form  penalty exact zero     train validation     test
ridge               3          0    5.0348     5.2982   5.2954
ridge              10          0    5.0367     5.3009   5.2968
ridge              30          0    5.0408     5.3088   5.3053
ridge             100          0    5.0581     5.3313   5.3342
ridge             300          0    5.1289     5.4180   5.4387
ridge            1000          0    5.3907     5.6894   5.7325
lasso               3         42    5.0371     5.3028   5.2933
lasso              10         65    5.0429     5.3092   5.2909
lasso              30         74    5.0515     5.3179   5.3077
lasso             100         79    5.0810     5.3486   5.3443
lasso             300         83    5.1447     5.4056   5.4420
lasso            1000         84    5.6036     5.8955   5.9568
elastic net         3         24    5.0360     5.3015   5.2947
elastic net        10         47    5.0405     5.3058   5.2941
elastic net        30         68    5.0477     5.3158   5.3022
elastic net       100         72    5.0722     5.3458   5.3483
elastic net       300         75    5.1431     5.4192   5.4509
elastic net      1000         83    5.4757     5.7450   5.8260

candidates tried 18

The third column is the answer of the lesson, and the difference between the three blocks is beyond dispute. In the rows of the ridge penalty that column is zero in all six settings: even when the penalty is raised a thousandfold, not one of the 87 coefficients becomes exactly zero. The lasso penalty rises from 42 to 84 over the same range. The elastic net penalty stands right in the middle: from 24 to 83.

The distinction comes from the force a coefficient feels as it approaches zero. The derivative of the ridge penalty is proportional to the coefficient; as the coefficient shrinks, the force trying to shrink it further shrinks too, and it never reaches zero. The derivative of the lasso penalty is independent of magnitude — however small the coefficient, it is pushed with the same constant force — and the moment that force exceeds the coefficient’s contribution to lowering the error, the coefficient stays at zero.

The last two columns read the same table from the error side, and there the difference is small. The best validation numbers of the three forms are 5.2982, 5.3028 and 5.3015; the distance between them is under 0.005 m³. The choice of penalty form is not visible in the error number; it is visible in the number of columns.

What Sparsity Is Worth

One of the eighteen candidates has to be chosen, and the price of sparsity is read right here.

best_val = min(RESULTS, key=lambda t: t[3])
best_test = min(RESULTS, key=lambda t: t[4])
print(f"chosen on validation {best_val[0]} penalty {best_val[1]}: "
      f"{best_val[2]} exact zero, test {best_val[4]:.4f}")
print(f"if chosen on test {best_test[0]} penalty {best_test[1]}: "
      f"{best_test[2]} exact zero, test {best_test[4]:.4f}")
print(f"difference between the two choices {best_val[4] - best_test[4]:.4f} m3\n")
ridge3, lasso10 = descent(3, 0.0), descent(10, 1.0)
print(f"{'threshold':>9} {'ridge penalty 3':>16} {'lasso penalty 10':>16}  "
      f"(coefficients below this absolute value)")
for e in (1e-9, 0.001, 0.01, 0.1):
    print(f"{e:>9} {sum(1 for v in ridge3[1] if abs(v) < e):>16} "
          f"{sum(1 for v in lasso10[1] if abs(v) < e):>16}")
kept = [NAMES[j] for j in range(P) if abs(lasso10[1][j]) > 1e-9]
print(f"\nlasso penalty 10 leaves {len(kept)} columns: {', '.join(kept[:9])} ...")
for c in (3000, 100000):
    for name, share in (("lasso", 1.0), ("ridge", 0.0)):
        mm = descent(c, share)
        print(f"{name} penalty {c}: exact zero {sum(1 for v in mm[1] if abs(v) < 1e-9)}/{P}, "
              f"intercept {mm[0]:.4f}, test {rmse(predictor(mm), TEST):.4f}")
print(f"train mean {MEAN:.4f}, baseline on test {rmse(lambda x: MEAN, TEST):.4f}")

chosen on validation ridge penalty 3: 0 exact zero, test 5.2954
if chosen on test lasso penalty 10: 65 exact zero, test 5.2909
difference between the two choices 0.0046 m3

threshold  ridge penalty 3 lasso penalty 10  (coefficients below this absolute value)
    1e-09                0               65
    0.001                0               65
     0.01                5               66
      0.1               29               71

lasso penalty 10 leaves 22 columns: area_south, area_east, area_west, area_central, p^1, s^1p^4, s^1p^5, s^3p^1, h^1 ...
lasso penalty 3000: exact zero 87/87, intercept 19.7665, test 6.5324
ridge penalty 3000: exact zero 0/87, intercept 19.7665, test 6.0935
lasso penalty 100000: exact zero 87/87, intercept 19.7665, test 6.5324
ridge penalty 100000: exact zero 0/87, intercept 19.7665, test 6.5086
train mean 19.7665, baseline on test 6.5324

The validation set chooses the ridge penalty, and that choice keeps all 87 coefficients. Had it been chosen from the test set, lasso 10 would have come: 65 coefficients exactly zero, the remaining 22 columns and 5.2909. The difference between the two choices is 0.0046 m³, and it is the price of searching eighteen times for the smallest on a 252-row set.

But this is not the difference the table is really telling. The same error level is reached with two different column counts: 5.2954 with 87 columns for the ridge penalty, 5.2909 with 22 columns for the lasso penalty. What sparsity is worth on the held-out set is zero on this set — 65 columns leave the model and the error number does not worsen. Sparsity is not a performance gain but a cost gain: the number of columns to be computed, stored and explained for a prediction falls to a quarter.

The threshold table gives the difference between the two penalties more sharply. Not one of the ridge penalty’s 87 coefficients is below $10^{-9}$, but 29 of them are below 0.1. That is, the ridge penalty also disables columns — it only keeps them in the model. The lasso penalty really drops 65 of them and rises only to 71 at the 0.1 threshold; it does not shrink the ones outside those it dropped.

The last four lines show the limit of the two forms. At lasso 3000 all 87 coefficients are zeroed, the intercept stays at the train mean, and the model gives exactly 6.5324 on test, that is, the baseline itself. The ridge penalty is at 6.0935 at the same setting, and at 6.5086 at a coefficient of a hundred thousand. The ridge penalty approaches the baseline but at a finite setting does not reach it; the lasso penalty does. The difference between the two forms shows not at any point of the error curve but at its ends.

Summary

• The form of the penalty shows not in the error number but in the number of columns: the best validation numbers of the three forms, 5.2982, 5.3028 and 5.3015, are within 0.005 m³ of one another.
• The ridge penalty does not set a single coefficient exactly to zero in any of the six settings; the lasso penalty rises from 42 to 84, the elastic net penalty from 24 to 83. The reason is the shape of the derivative: the force of the ridge penalty shrinks as the coefficient shrinks, that of the lasso penalty stays constant.
• The same error level is reached with two different column counts: 5.2954 with 87 columns for the ridge penalty, 5.2909 with 22 columns for the lasso penalty. What sparsity is worth on this set is a cost, not performance.
• The ridge penalty disables columns but does not drop them: 29 of its 87 coefficients are below 0.1 and none is exactly zero.
• At lasso 3000 all the coefficients are zeroed and the model returns exactly to the baseline, 6.5324; the ridge penalty stays at 6.5086 even at a coefficient of a hundred thousand and does not reach the baseline at a finite setting.
• The choice read from validation is ridge 3 (5.2954 on test); had it been read from test, lasso 10 would have been chosen (5.2909) and the price of the eighteen-candidate budget would have been 0.0046 m³.

Next Step

In three lessons the target was a continuous number and the model said that number directly. The real through-line problem of the course, though, is a classification: whether a subscriber carries illicit use or fault suspicion, and the baseline there is 0.7579. The same linear model cannot be used directly on this problem, because the output of the plane stretches from minus infinity to plus infinity and the label is only zero and one. The bridge between them is squeezing the output into a probability. The next lesson writes this squeeze, and its measure is this: the same model gives one number at a 0.5 threshold and other numbers at other thresholds; choosing the threshold is a separate job and is left to the evaluation-and-tuning topic. What will be shown here is that the output is continuous before the threshold and that this continuity carries more than a single class label.