Regularization

The previous two lessons pointed in the same direction: as the training loss fell, the held-out set number tended to fall. This is a condition with a known name — when a few hundred parameters meet a training set of 300 rows, the network memorizes the rows. Overfitting itself was measured in the Generalization topic of Introduction to Machine Learning and is not repeated here. The question here is narrower: which countermeasures are specific to the network, and how many subscribers does each one earn?

Weight decay is not the first answer to this question and it is not swept again in this course. Regularization Intuition (04) in the Generalization topic of Introduction to Machine Learning had moved the penalty coefficient from 0.001 to 1000 and measured the weight magnitude going from 10.040 to 0.601, and Regularization Penalties (03) in Supervised Learning had measured the squared, absolute and mixed penalties separating on sparsity. Here the penalty is used with a single coefficient and only as a comparison line. The two methods this lesson adds are dropout and early stopping.

A distinction of terms: dropout in this course is temporarily silencing a part of the hidden units during training. The downsampling in the Chart Design topic of Visualization, Reporting and Analysis at Scale was a drawing technique and reduced the number of points; the two share only a name.

• NN61. The set, the split, the seed and the scale are as in the previous lessons; the table is fictional, the training set has been reduced to 300 rows and the baseline on the test set is 0.7579.
• NN62. In this lesson the network has 16 hidden units and carries 225 parameters. It was widened because dropout needs units to close.
• NN63. Training is 45 epochs, batch size 20, plain descent, learning rate 0.05. These settings were not searched for; they come from the previous lesson.
• NN64. Dropout is applied with a mask redrawn for every mini-batch: the output of a closed unit is zero and the output of the remaining ones is scaled up by $1/(1-p)$. Because the scale is applied in training, no correction is needed at inference; outside training the network runs without a mask.
• NN65. The mask is per mini-batch, not per row. This is a decision and it keeps the number of masks at 45 times 15, that is, 675.
• NN66. The weight decay coefficient is 0.01 and it is added to the gradient as the coefficient times the weight. Biases are not penalized: the bias carries the position of the decision boundary, not its curvature.
• NN67. The criterion of early stopping is the validation loss, not accuracy; accuracy measures in steps of 0.0040 on a set of 252 rows and leaves the bottom unclear.
• NN68. The number of settings tried is 5; this is not a search but the four countermeasures placed side by side with the bare network. The comparison lines come from Supervised Learning: the linear model 0.8016, the decision tree 0.7857.

Two Methods Specific to the Network

The setup is the same as in the previous lessons; what is added is that the forward pass takes a mask and backpropagation can put a penalty term on the gradient. The output of a masked unit is zero, so its gradient is zero as well; the gradient of the remaining units is scaled by the scaling factor in the mask.

# regularization.py — MODEL: the same FICTIONAL subscriber table, the same seed, the same split;
# the training set is reduced to 300 rows in this topic.
import math
import statistics

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


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:
        continue
    sub_id = 10001 + i
    r, v = rng(SEED + 7000 + sub_id), []
    for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            b[2] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - d * 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 + sub_id)
    z = (0.052 * (mean - 20) + 0.85 * (min(v) == 0.0) + 0.026 * vol + 0.24 * (household >= 5)
         - 0.20 * (rating >= 4) + b[3] + 0.9 * (sub_id - 10001) / 1399
         + (r() + r() + r() - 1.5) * 1.30)
    x = {"household": household, "satisfaction": rating, "period": len(v),
         "mean_consumption": round(mean, 2), "volatility": round(vol, 2),
         "zero_readings": int(min(v) == 0.0),
         "suspicious": int(z > 1.35)}
    x["per_person"] = round(x["mean_consumption"] / x["household"], 3)
    for q in REGION:
        x["area_" + q[0]] = int(b[0] == q[0])
    DATA.append(x)

r, S = rng(SEED + 90000), 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]
K = [DATA[i] for i in S]
TRAIN, VAL, TEST = K[:300], K[756:1008], K[1008:]      # training REDUCED from 756 to 300
FIELDS = ["mean_consumption", "volatility", "household", "satisfaction", "period",
          "zero_readings", "per_person", "area_north", "area_south", "area_east",
          "area_west", "area_central"]
MEAN = {a: statistics.fmean(x[a] for x in TRAIN) for a in FIELDS}
STD = {a: max(1e-9, statistics.pstdev([x[a] for x in TRAIN])) for a in FIELDS}


def scale(s):                               # LEARNING from training, APPLICATION to every set
    return ([[(x[a] - MEAN[a]) / STD[a] for a in FIELDS] for x in s],
            [float(x["suspicious"]) for x in s])


X_TR, Y_TR = scale(TRAIN)
X_VA, Y_VA = scale(VAL)
X_TE, Y_TE = scale(TEST)
BASELINE = sum(y == 0.0 for y in Y_TE) / len(Y_TE)
TREE, LINEAR = 198 / len(TEST), 202 / len(TEST)     # the two numbers of the Supervised Learning course


def squash(z):
    if z >= 0:
        return 1.0 / (1.0 + math.exp(-z)) if z < 700 else 1.0
    e = math.exp(z) if z > -700 else 0.0
    return e / (1.0 + e)


def init_weights(seed, shape, gain=1.0):    # gain-scaled; the last column of every row is the BIAS
    r, layers = rng(seed), []
    for fan_in, fan_out in shape:
        limit = math.sqrt(6.0 * gain / fan_in)
        layers.append([[(r() * 2 - 1) * limit for _ in range(fan_in)] + [0.0]
                       for _ in range(fan_out)])
    return layers


def zeros(layers):
    return [[[0.0] * len(w) for w in W] for W in layers]


def make_mask(layers, rate, seed):          # DROPOUT: a closed unit is 0, the rest 1/(1-rate)
    r = rng(seed)
    return [[0.0 if r() < rate else 1.0 / (1.0 - rate) for _ in W] for W in layers[:-1]]


def forward(layers, x, mask=None):
    a, trace = x, []
    for i, W in enumerate(layers):
        trace.append(a)
        z = [sum(p * q for p, q in zip(w, a)) + w[-1] for w in W]
        if i == len(layers) - 1:
            a = [squash(z[0])]
        else:
            a = [t if t > 0 else 0.0 for t in z]
            if mask:
                a = [v * m for v, m in zip(a, mask[i])]
    trace.append(a)
    return trace


def gradient(layers, X, Y, idx, mask=None, penalty=0.0):      # BACKPROPAGATION
    G, loss = zeros(layers), 0.0
    for i in idx:
        trace = forward(layers, X[i], mask)
        p = min(max(trace[-1][0], 1e-12), 1 - 1e-12)
        loss -= Y[i] * math.log(p) + (1 - Y[i]) * math.log(1 - p)
        d = [p - Y[i]]
        for L in range(len(layers) - 1, -1, -1):
            W, a = layers[L], trace[L]
            for k, dk in enumerate(d):
                gk = G[L][k]
                for j, aj in enumerate(a):
                    gk[j] += dk * aj
                gk[-1] += dk
            if L:
                mk = mask[L - 1] if mask else None
                d = [sum(d[k] * W[k][j] for k in range(len(W)))
                     * (0.0 if a[j] <= 0 else (mk[j] if mk else 1.0))
                     for j in range(len(a))]
    for W in G:
        for w in W:
            w[:] = [v / len(idx) for v in w]
    if penalty:                             # WEIGHT DECAY: the bias is not penalized
        for W, gW in zip(layers, G):
            for w, g in zip(W, gW):
                for j in range(len(w) - 1):
                    g[j] += penalty * w[j]
    return G, loss / len(idx)


def loss_of(layers, X, Y):
    t = 0.0
    for x, y in zip(X, Y):
        p = min(max(forward(layers, x)[-1][0], 1e-12), 1 - 1e-12)
        t -= y * math.log(p) + (1 - y) * math.log(1 - p)
    return t / len(X)


def accuracy(layers, X, Y):
    return sum((forward(layers, x)[-1][0] > 0.5) == (y > 0.5) for x, y in zip(X, Y)) / len(X)


def order(epoch):                           # batch split order is the SAME in every run
    r, s = rng(SEED + 500 + epoch), list(range(len(X_TR)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return s


LAYERS0 = init_weights(SEED + 1, [(len(FIELDS), 16), (16, 1)])
print(f"network {len(FIELDS)}-16-1, parameters {sum(len(w) for W in LAYERS0 for w in W)}, "
      f"train {len(TRAIN)}, validation {len(VAL)}, test {len(TEST)}, epochs 45, batch 20")
print(f"baseline {BASELINE:.4f}; comparison lines linear model {LINEAR:.4f}, "
      f"decision tree {TREE:.4f}")

network 12-16-1, parameters 225, train 300, validation 252, test 252, epochs 45, batch 20
baseline 0.7579; comparison lines linear model 0.8016, decision tree 0.7857

Dropout and Weight Decay

The training loop is in a single function and takes both countermeasures as switches. The measure is three columns: training accuracy, test accuracy and the gap between the two. The gap is the direct measure of overfitting; if a regularization method is working it narrows the gap first.

def train_model(seed, rate=0.0, penalty=0.0, watch=False, step=0.05, epochs=45, batch=20,
                hidden=16):
    layers, curve = init_weights(seed, [(len(FIELDS), hidden), (hidden, 1)]), []
    for t in range(epochs):
        s = order(t)
        for bi, b in enumerate(range(0, len(s), batch)):
            mk = make_mask(layers, rate, SEED + 900 + t * 100 + bi) if rate else None
            G, _ = gradient(layers, X_TR, Y_TR, s[b:b + batch], mk, penalty)
            for W, gW in zip(layers, G):
                for w, g in zip(W, gW):
                    for j in range(len(w)):
                        w[j] -= step * g[j]
        if watch:                           # the held-out set is read at the end of every epoch
            curve.append((loss_of(layers, X_TR, Y_TR), loss_of(layers, X_VA, Y_VA),
                          accuracy(layers, X_TE, Y_TE)))
    return layers, curve


SETTINGS = (("bare", 0.0, 0.0), ("dropout 0.2", 0.2, 0.0), ("dropout 0.4", 0.4, 0.0),
            ("weight decay", 0.0, 0.01), ("dropout 0.4 + decay", 0.4, 0.01))
print(f"{'setting':>21} {'train':>15} {'test':>15} {'gap':>7} {'run':>7}")
OUTCOME = {}
for name, rate, penalty in SETTINGS:
    nets = [train_model(SEED + s, rate, penalty)[0] for s in (1, 2)]
    tr = [accuracy(k, X_TR, Y_TR) for k in nets]
    OUTCOME[name] = [accuracy(k, X_TE, Y_TE) for k in nets]
    print(f"{name:>21} {tr[0]:>7.4f} {tr[1]:>7.4f} {OUTCOME[name][0]:>7.4f} {OUTCOME[name][1]:>7.4f} "
          f"{tr[0] - OUTCOME[name][0]:>7.4f} {abs(OUTCOME[name][0] - OUTCOME[name][1]):>7.4f}")
print(f"\nsettings tried {len(SETTINGS)}; baseline {BASELINE:.4f}, network number "
      f"{OUTCOME['dropout 0.4'][0]:.4f}, second run of the same network {OUTCOME['dropout 0.4'][1]:.4f}")
print(f"runs beating the decision tree {sum(v > TREE + 1e-9 for n in OUTCOME.values() for v in n)}/10, "
      f"the linear model {sum(v > LINEAR + 1e-9 for n in OUTCOME.values() for v in n)}/10")

              setting           train            test     gap     run
                 bare  0.8467  0.8467  0.7579  0.7817  0.0887  0.0238
          dropout 0.2  0.8267  0.8300  0.7778  0.7937  0.0489  0.0159
          dropout 0.4  0.8267  0.8133  0.7937  0.7738  0.0330  0.0198
         weight decay  0.8400  0.8433  0.7698  0.7937  0.0702  0.0238
  dropout 0.4 + decay  0.8233  0.8133  0.7897  0.7738  0.0337  0.0159

settings tried 5; baseline 0.7579, network number 0.7937, second run of the same network 0.7738
runs beating the decision tree 4/10, the linear model 0/10

The gap column behaves regularly. In the bare network there are 0.0887 between training and test — twenty-two subscribers. With dropout rate 0.2 this goes down to 0.0489, with 0.4 to 0.0330. A narrowing of almost three times shows that dropout is doing its job: because every mini-batch works with different units, no unit can memorize a single training row on its own. Weight decay brings the same gap down only to 0.0702; in a setting where the penalty is not swept this is an expected result, and what the sweep earns was measured in Introduction to Machine Learning.

The test column confirms the gap in the first run: from 0.7579 to 0.7937, nine subscribers. But the second run does not confirm it. There the bare network gives 0.7817 and the dropout network 0.7738 — the ranking flips. Someone looking at the first run says “dropout 0.4 beats the bare network by 0.0358”, someone looking at the second run says “it trails by 0.0079”. Both ran the same code. The last line sustains the second claim of the course: four of the ten runs beat the decision tree, none beats the linear model.

Which Set Early Stopping Is Read From

Early stopping is not a setting but a reading point: after how many epochs the training will be cut. The real question is which set is looked at to find the epoch. The block below puts three sources side by side.

for s in (1, 2):
    _, curve = train_model(SEED + s, watch=True)
    val = min(range(len(curve)), key=lambda t: curve[t][1])     # the bottom of the validation loss
    tst = max(range(len(curve)), key=lambda t: curve[t][2])     # the peak of the test
    print(f"run {s}: training loss {curve[0][0]:.4f} at epoch 1, {curve[-1][0]:.4f} at the last "
          f"epoch; validation loss {curve[0][1]:.4f} -> {curve[-1][1]:.4f}")
    print(f"  stopping epoch read from training 45: test {curve[-1][2]:.4f}")
    print(f"  stopping epoch read from validation {val + 1}: validation loss "
          f"{curve[val][1]:.4f}, test {curve[val][2]:.4f}")
    print(f"  stopping epoch read from test {tst + 1}: test {curve[tst][2]:.4f}, "
          f"inflation {curve[tst][2] - curve[val][2]:.4f} "
          f"(subscribers {round((curve[tst][2] - curve[val][2]) * len(TEST))})")

run 1: training loss 0.5183 at epoch 1, 0.3549 at the last epoch; validation loss 0.5244 -> 0.4505
  stopping epoch read from training 45: test 0.7579
  stopping epoch read from validation 10: validation loss 0.4475, test 0.7698
  stopping epoch read from test 11: test 0.7738, inflation 0.0040 (subscribers 1)
run 2: training loss 0.6826 at epoch 1, 0.3478 at the last epoch; validation loss 0.6677 -> 0.4244
  stopping epoch read from training 45: test 0.7817
  stopping epoch read from validation 14: validation loss 0.4129, test 0.7976
  stopping epoch read from test 11: test 0.8016, inflation 0.0040 (subscribers 1)

The two loss curves part ways. The training loss never rises across forty-five epochs: in the first run it goes from 0.5183 to 0.3549. The validation loss, on the other hand, sees its bottom at the tenth epoch with 0.4475, turns upward and rises to 0.4505 at the last epoch. A stopping epoch cannot be chosen by looking at the training set, because there is no place to stop there — the bottom of the curve is always at the last epoch.

The epoch read from the validation set is 10 in the first run and 14 in the second. Their counterparts on the test are 0.7698 and 0.7976; had forty-five epochs been run they would have been 0.7579 and 0.7817. The contribution of early stopping is 0.0119 in the first run and 0.0159 in the second — that is, three and four subscribers. In the first run the bare network had fallen exactly to the baseline; early stopping lifts it above the baseline.

The last lines are the real warning of the lesson. Had the epoch been chosen by looking at the test set, it would have been 11 in both runs and the reported numbers would have been 0.7738 and 0.8016. The inflation is 0.0040, that is, one subscriber. It is small, but it is an unmeasured gain: the moment the test set is read once, it stops being held out. The same logic was built in the pipelines lesson of Data Preparation and Feature Engineering to keep a transformation from being learned from the whole set; the early stopping epoch carries as much information as a transformation.

Summary

• Weight decay had been swept in Introduction to Machine Learning and Supervised Learning; here it is applied with a single coefficient and brings the gap between training and test from 0.0887 down to 0.0702.
• Dropout brings the same gap down to 0.0330 at rate 0.4; because every mini-batch works with different units, a single unit cannot memorize a single row.
• Because the scaling is done in training, no correction is needed at inference; outside training the network runs without a mask.
• The test column shows dropout ahead by 0.0358 in the first run and behind by 0.0079 in the second; four of the ten runs beat the decision tree, none beats the linear model.
• The validation loss sees its bottom at the 10th and the 14th epoch; when the stopping epoch is read from validation the test comes out as 0.7698 and 0.7976, from test as 0.7738 and 0.8016 — the inflation is 0.0040 in both runs, one subscriber.

Next Step

Both methods of this lesson touched the output of the network: one silenced units, the other shrank weights. The input of the layers was never touched. Yet input scaling was the only preparation step used throughout the course and it was applied only to the first layer; the numbers the second layer sees drift as training goes on and nobody holds them in place. The next lesson ties this to a measure: when layer inputs are rescaled at every batch, what does the network do at a learning rate that is unstable without normalization. The same lesson also opens a source of leakage, because normalization does not use the same number in training and in inference.