Batch Normalization

Both methods of the previous lesson touched the output of the network; the input of the layers was never touched. Yet scaling was the only preparation step of this course and it was applied only to the first layer: twelve columns were scaled with the mean and variability learned from the training set. The numbers the hidden layer sees have no such guarantee; because the weights change at every update, the input scale of the second layer drifts as well.

Batch normalization fills this gap: the inputs of a layer are rescaled at every step with the mean and variability of the mini-batch they are in. The measure of this lesson is the learning rate that normalization allows. The second measure is its cost: normalization does not use the same number in training and in inference, and this difference is a source of leakage.

• NN69. 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.
• NN70. The network is 12-16-1 with 225 parameters, as in the previous lesson. Training is 30 epochs, batch size 20, plain descent. The only things that change are the learning rate and whether there is normalization.
• NN71. Normalization is applied to the pre-activation values of the hidden layer; the learned scale and shift coefficients are not used. To avoid dividing by zero, $10^{-8}$ is added to the variability and this threshold moves nothing at the visible digit.
• NN72. The derivative of the normalization itself is written into backpropagation: because the mean and the deviation come from all the rows of the batch, the gradient of one row depends on the other rows as well.
• NN73. At inference there are no batch statistics; a running average with coefficient 0.9 is accumulated during training and that is used at inference. The dead unit is the definition of the previous lesson.
• NN74. The number of settings tried is 8: two normalization states times four learning rates. This is not a search; the rates were chosen to show the limit of stability. The comparison lines come from Supervised Learning: the linear model 0.8016, the decision tree 0.7857.

Scaling Within the Mini-Batch

The code of this lesson writes the forward pass and backpropagation batch by batch, not row by row, because the statistic to be normalized comes from the whole of the batch. The normalization derivative in backpropagation is a consequence of this: the gradient of one row contains the values of the other rows of the same batch.

# normalization.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, TEST = K[:300], 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_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 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


def forward_batch(W1, W2, Xb, stats=None, norm=False):   # THE WHOLE OF ONE MINI-BATCH
    Z = [[sum(p * q for p, q in zip(w, x)) + w[-1] for w in W1] for x in Xb]
    mu = sd = None
    if norm:
        n, m = len(Z), len(W1)
        if stats is None:                   # TRAINING: the batch's own statistics
            mu = [sum(z[j] for z in Z) / n for j in range(m)]
            sd = [math.sqrt(sum((z[j] - mu[j]) ** 2 for z in Z) / n + 1e-8)
                  for j in range(m)]
        else:                               # INFERENCE: running averages
            mu, sd = stats
        Z = [[(z[j] - mu[j]) / sd[j] for j in range(m)] for z in Z]
    A = [[t if t > 0 else 0.0 for t in z] for z in Z]
    P = [squash(sum(p * q for p, q in zip(W2[0], a)) + W2[0][-1]) for a in A]
    return Z, mu, sd, A, P


def gradient_batch(W1, W2, Xb, Yb, norm=False):      # BACKPROPAGATION
    Z, mu, sd, A, P = forward_batch(W1, W2, Xb, norm=norm)
    n, m, loss = len(Xb), len(W1), 0.0
    g1 = [[0.0] * len(W1[0]) for _ in W1]
    g2 = [[0.0] * len(W2[0])]
    dZ = [[0.0] * m for _ in range(n)]
    for i in range(n):
        p = min(max(P[i], 1e-12), 1 - 1e-12)
        loss -= Yb[i] * math.log(p) + (1 - Yb[i]) * math.log(1 - p)
        d = p - Yb[i]
        for j in range(m):
            g2[0][j] += d * A[i][j] / n
            dZ[i][j] = d * W2[0][j] * (Z[i][j] > 0)
        g2[0][-1] += d / n
    if norm:                                # the derivative of the normalization itself
        for j in range(m):
            s1 = sum(dZ[i][j] for i in range(n))
            s2 = sum(dZ[i][j] * Z[i][j] for i in range(n))
            for i in range(n):
                dZ[i][j] = (n * dZ[i][j] - s1 - Z[i][j] * s2) / (n * sd[j])
    for i in range(n):
        for j in range(m):
            if dZ[i][j]:
                for k, xk in enumerate(Xb[i]):
                    g1[j][k] += dZ[i][j] * xk / n
                g1[j][-1] += dZ[i][j] / n
    return g1, g2, loss / n, mu, sd


def predict(W1, W2, X, stats=None, norm=False, batch=0):
    if not batch:
        return forward_batch(W1, W2, X, stats, norm)[4]
    return [p for b in range(0, len(X), batch)
            for p in forward_batch(W1, W2, X[b:b + batch], None, norm)[4]]


def accuracy(P, Y):
    return sum((p > 0.5) == (y > 0.5) for p, y in zip(P, Y)) / len(Y)


def dead_units(W1, W2, run_mean, run_std, norm):    # a unit that never activates on any training row
    A = forward_batch(W1, W2, X_TR, (run_mean, run_std) if norm else None, norm)[3]
    return sum(all(a[j] <= 0 for a in A) for j in range(len(W1)))


W1, W2 = init_weights(SEED + 1, [(len(FIELDS), 16), (16, 1)])
print(f"network {len(FIELDS)}-16-1, parameters {sum(len(w) for w in W1 + W2)}, train {len(TRAIN)}, "
      f"test {len(TEST)}, epochs 30, 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, test 252, epochs 30, batch 20
baseline 0.7579; comparison lines linear model 0.8016, decision tree 0.7857

The Learning Rate That Normalization Allows

The training loop also accumulates the running averages, because at inference there is no batch.

def train_model(seed, rate, norm=False, epochs=30, batch=20, hidden=16):
    W1, W2 = init_weights(seed, [(len(FIELDS), hidden), (hidden, 1)])
    run_mean, run_std, last = [0.0] * hidden, [1.0] * hidden, []
    for t in range(epochs):
        s, total = order(t), 0.0
        for b in range(0, len(s), batch):
            i = s[b:b + batch]
            g1, g2, k, mu, sd = gradient_batch(W1, W2, [X_TR[q] for q in i],
                                               [Y_TR[q] for q in i], norm)
            for w, g in zip(W1, g1):
                for j in range(len(w)):
                    w[j] -= rate * g[j]
            for j in range(len(W2[0])):
                W2[0][j] -= rate * g2[0][j]
            if norm:                        # running averages are accumulated for inference
                run_mean = [0.9 * a + 0.1 * b for a, b in zip(run_mean, mu)]
                run_std = [0.9 * a + 0.1 * b for a, b in zip(run_std, sd)]
            total += k * len(i)
        last.append(total / len(s))
    return W1, W2, run_mean, run_std, last[-1]


print(f"{'normalization':>15} {'rate':>6} {'loss 2 runs':>16} {'test 2 runs':>16} "
      f"{'dead':>11}")
OUTCOME = {}
for norm in (False, True):
    for rate in (0.05, 5.0, 20.0, 60.0):
        r = [train_model(SEED + s, rate, norm) for s in (1, 2)]
        n = [accuracy(predict(a, b, X_TE, (c, d), norm), Y_TE) for a, b, c, d, _ in r]
        OUTCOME[(norm, rate)] = n
        print(f"{'yes' if norm else 'no':>15} {rate:>6} {r[0][4]:>7.4f} {r[1][4]:>8.4f} "
              f"{n[0]:>8.4f} {n[1]:>7.4f} {dead_units(*r[0][:4], norm):>7}/16 "
              f"{dead_units(*r[1][:4], norm):>2}/16")
print(f"\nsettings tried {len(OUTCOME)}; baseline {BASELINE:.4f}, network number "
      f"{OUTCOME[(True, 0.05)][0]:.4f}, second run of the same network {OUTCOME[(True, 0.05)][1]:.4f}")

  normalization   rate      loss 2 runs      test 2 runs        dead
             no   0.05  0.3740   0.3671   0.7659  0.7817       0/16  0/16
             no    5.0  0.3841   0.3717   0.6905  0.5992       1/16  2/16
             no   20.0  8.6577   8.1051   0.6270  0.6190       0/16  0/16
             no   60.0  6.0788   6.0788   0.7579  0.7579      14/16 16/16
            yes   0.05  0.4020   0.3922   0.7698  0.7857       0/16  0/16
            yes    5.0  0.3809   0.4014   0.7381  0.7579       0/16  0/16
            yes   20.0  0.7206   0.6698   0.6984  0.6865       0/16  0/16
            yes   60.0  2.3263   1.9161   0.6429  0.6587       0/16  0/16

settings tried 8; baseline 0.7579, network number 0.7698, second run of the same network 0.7857

The loss column gives the number in the title of the lesson. The network without normalization is still standing at 5, at 20 the loss goes up to 8.6577; the network with normalization stays at 0.7206 at the same 20 and does not exceed 2.3263 even at 60. The largest rate that keeps the loss under 1.0 is 5 without normalization and 20 with it: fourfold. This is what normalization allows.

The dead unit column shows the difference more sharply. The network without normalization loses 14 of 16 units at 60 and all of them in the second run, and both runs give exactly the baseline. The network with normalization loses not a single unit at any rate: rescaling the input at every step keeps the pre-activation values from being pushed permanently below zero.

The test column carries a warning: normalization makes the large rate possible but does not make it useful, because the highest number in both families is at the smallest rate. The triple is this — baseline 0.7579, network number 0.7698, second run of the same network 0.7857; the second is 0.0159 higher than the first and lands exactly on the number of the decision tree.

Different in Training, Different at Inference

The cost of normalization is the splitting of behavior in two. In training the output of a row depends on the nineteen rows next to it; at inference it should not. The block below makes the inference in three ways.

W1, W2, RUN_MEAN, RUN_STD, _ = train_model(SEED + 1, 0.05, True)
P_A = predict(W1, W2, X_TE, (RUN_MEAN, RUN_STD), True)  # the right way: running averages
P_B = predict(W1, W2, X_TE, norm=True, batch=20)        # the batch's own statistics
r, S2 = rng(SEED + 4242), list(range(len(X_TE)))
for i in range(len(S2) - 1, 0, -1):
    j = int(r() * (i + 1))
    S2[i], S2[j] = S2[j], S2[i]
P_C = [0.0] * len(X_TE)
for k, i in enumerate(predict(W1, W2, [X_TE[q] for q in S2], norm=True, batch=20)):
    P_C[S2[k]] = i
print(f"inference with running averages   {accuracy(P_A, Y_TE):.4f}")
print(f"inference with batch statistics   {accuracy(P_B, Y_TE):.4f}  "
      f"inflation {accuracy(P_B, Y_TE) - accuracy(P_A, Y_TE):.4f} "
      f"(subscribers {round((accuracy(P_B, Y_TE) - accuracy(P_A, Y_TE)) * len(TEST))})")
print(f"same path, test order changed     {accuracy(P_C, Y_TE):.4f}")
print(f"subscribers who changed class: between running averages and batch statistics "
      f"{sum((a > 0.5) != (b > 0.5) for a, b in zip(P_A, P_B))}, "
      f"only when the order changes {sum((a > 0.5) != (b > 0.5) for a, b in zip(P_B, P_C))}")
print(f"largest movement in the probability of a single subscriber "
      f"{max(abs(a - b) for a, b in zip(P_B, P_C)):.4f}")

inference with running averages   0.7698
inference with batch statistics   0.7937  inflation 0.0238 (subscribers 6)
same path, test order changed     0.7817
subscribers who changed class: between running averages and batch statistics 6, only when the order changes 7
largest movement in the probability of a single subscriber 0.2273

The second line looks like good news at first glance: 0.7937 instead of 0.7698, six more subscribers correct. But those six subscribers do not come from learning. The mean of the test batch was computed from the test rows themselves; the prediction of every subscriber contains the numbers of the other nineteen subscribers in the same batch.

The third line proves it. The model is the same, the weights are the same, the test set is the same — only the order of the rows changed and the number went down from 0.7937 to 0.7817. Seven subscribers changed class, the probability of one subscriber moved by 0.2273. The prediction of a model should depend on no one but the person being predicted; here it does. The pipelines lesson of Data Preparation and Feature Engineering had measured that learning a transformation from the whole set produced an inflation between 0.0397 and 0.0476; batch normalization makes the same mistake possible silently. There is only one right way: running averages at inference.

Summary

• Batch normalization scales the pre-activation values of a layer with the mean and variability of the mini-batch it is in; its derivative also comes from the whole of the batch.
• The largest learning rate that keeps the loss under 1.0 is 5 without normalization and 20 with it — fourfold. The network without normalization carries the loss to 8.6577 at 20.
• At rate 60, 14 of the 16 units of the network without normalization die, and in the second run all of them, and the result is exactly the baseline; with normalization no unit dies at any rate.
• The large rate becomes possible but not useful: the highest test number is at the smallest rate in both families. The triple is baseline 0.7579, network number 0.7698, second run 0.7857.
• When inference is made with the statistics of the batch itself the number inflates to 0.7937; when only the test order changes it goes down to 0.7817 and one probability moves by 0.2273.

Next Step

Nine lessons built a network piece by piece: nonlinear separation, activation, forward propagation, loss, backpropagation, optimizer, stability, regularization and normalization. In every lesson three numbers stood side by side and the third — the second run of the same network — never came out the same as the first. The rule of the course was confirmed in the last four lessons: in each lesson the capacity was held constant and only one setting was changed, and the optimizer, the learning rate, dropout and normalization moved the number between 0.5992 and 0.7937. The second claim of the course also holds: no run beat 0.8016 of the linear model.

There is one side of this network that did not change in nine lessons. Every input column was connected to every hidden unit and every connection carried its own weight — the network stayed fully connected. This is behaving as if the input columns had no relation to one another: if the order of the columns were shuffled the network would give the same result. Yet most real inputs have structure; in a series of meter readings neighboring periods, in a coverage grid neighboring cells are neighbors. The next topic writes this structure into the network itself and asks the first question with a parameter count: how many weights does a fully connected layer need to do the same job, and how many does a layer that knows about neighborhood need.