Training Stability

The previous lesson compared three step rules in a network with two layers. There the gradient was chained once on its way back, so every layer saw a gradient of a similar magnitude. When the number of layers grows this stops being true: while the error at the output moves toward the input it is filtered by a factor at every layer, and the product of those factors decides what is left for the first layer. If the factors are smaller than one the product goes down to zero, if they are larger than one it overflows.

This lesson measures that product. The two extremes are counted separately — the vanishing gradient and the exploding gradient — and then two corrections are tried: gradient clipping and the learning rate schedule. The measure is again twofold: rescuing the training and improving the held-out set number are not the same thing.

• NN54. The set, the split, the seed and the scale are the same as in the previous lesson; the table is fictional, the training set has been reduced to 300 rows and the baseline on the test set is 0.7579.
• NN55. Every hidden layer has 8 units: one layer comes to 113 parameters, two layers to 185, four layers to 329, six layers to 473.
• NN56. The gradient magnitude is the root of the mean of the squares of the gradient entries in a layer; because layers carry different numbers of weights, the mean is taken, not the sum.
• NN57. The measurement is made right after initialization, on a single batch of 60 rows; the question is not the result of the training but how the gradient is distributed across the layers. Two activations stand side by side: the sigmoid (saturating) and the rectified linear unit (non-saturating). Gain is the factor that sets the width of the initialization range.
• NN58. In the stability experiment the network has three hidden layers, makes 450 updates with 30 epochs and a batch size of 20, and uses plain descent. The number of settings tried is 5; this is not a search.
• NN59. The gradient clipping threshold is 0.10 and is applied once to the gradient magnitude of the whole network: if the magnitude exceeds the threshold the gradient is reduced by the same ratio, and the direction is kept. The schedule halves the learning rate every six epochs.
• NN60. A dead unit is a unit that gives an output greater than zero on no training row; across the three hidden layers there are 24 units in total.

The Gradient on Its Way Through the Layers

The setup is the same as in the previous lesson; the only difference is that the forward pass and backpropagation take the activation as a variable. The derivative of the sigmoid is computed from the output itself, the derivative of the rectified linear unit is zero or one.

# stability.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_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 forward(layers, x, activation="rectified"):
    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])]
        elif activation == "sigmoid":
            a = [squash(t) for t in z]
        else:
            a = [t if t > 0 else 0.0 for t in z]
    trace.append(a)
    return trace


def gradient(layers, X, Y, idx, activation="rectified"):      # BACKPROPAGATION
    G, loss = zeros(layers), 0.0
    slope = (lambda v: v * (1 - v)) if activation == "sigmoid" else (lambda v: v > 0)
    for i in idx:
        trace = forward(layers, X[i], activation)
        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:
                d = [sum(d[k] * W[k][j] for k in range(len(W))) * slope(a[j])
                     for j in range(len(a))]
    for W in G:
        for w in W:
            w[:] = [v / len(idx) for v in w]
    return G, loss / len(idx)


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 layer_shape(depth, hidden=8):           # depth hidden layers, then a single output
    return [(len(FIELDS), hidden)] + [(hidden, hidden)] * (depth - 1) + [(hidden, 1)]


print(f"train {len(TRAIN)} (reduced from 756), test {len(TEST)}, columns {len(FIELDS)}")
print("hidden layers and parameters: " + ", ".join(
    f"{d} -> {sum(len(w) for W in init_weights(SEED, layer_shape(d)) for w in W)}"
    for d in (1, 2, 3, 4, 6)))
print(f"baseline {BASELINE:.4f}; comparison lines linear model {LINEAR:.4f}, "
      f"decision tree {TREE:.4f}")

train 300 (reduced from 756), test 252, columns 12
hidden layers and parameters: 1 -> 113, 2 -> 185, 3 -> 257, 4 -> 329, 6 -> 473
baseline 0.7579; comparison lines linear model 0.8016, decision tree 0.7857

The Numbers of Vanishing and Exploding Gradients

At every depth the network is initialized, the gradient is computed on a single batch and the gradient magnitude of the first layer is printed next to that of the last layer. The last column gives how many times the first layer falls relative to the single-layer network of the same family.

def magnitude(W):                           # root of the mean of squares
    return math.sqrt(sum(v * v for w in W for v in w) / sum(len(w) for w in W))


print(f"{'activation':>14} {'gain':>7} {'depth':>9} {'first_layer':>12} "
      f"{'last_layer':>12} {'last/first':>12} {'vs_first':>10}")
FIRST = {}
for activation, gain in (("sigmoid", 1.0), ("rectified", 2.0), ("rectified", 16.0)):
    for depth in (1, 2, 4, 6):
        layers = init_weights(SEED + 1, layer_shape(depth), gain)
        G, _ = gradient(layers, X_TR, Y_TR, range(60), activation)
        b = [magnitude(W) for W in G]
        FIRST.setdefault((activation, gain), b[0])
        print(f"{activation:>14} {gain:>7} {depth:>9} {b[0]:>12.3e} {b[-1]:>12.3e} "
              f"{b[-1] / b[0]:>12.1f} {b[0] / FIRST[(activation, gain)]:>10.4f}")

    activation    gain     depth  first_layer   last_layer   last/first   vs_first
       sigmoid     1.0         1    1.274e-02    2.393e-01         18.8     1.0000
       sigmoid     1.0         2    2.820e-03    4.822e-02         17.1     0.2213
       sigmoid     1.0         4    1.930e-04    1.318e-01        682.6     0.0151
       sigmoid     1.0         6    5.539e-05    2.203e-01       3977.2     0.0043
     rectified     2.0         1    9.817e-02    5.694e-01          5.8     1.0000
     rectified     2.0         2    5.689e-02    2.176e-01          3.8     0.5796
     rectified     2.0         4    2.513e-02    5.219e-01         20.8     0.2560
     rectified     2.0         6    6.893e-02    3.273e-01          4.7     0.7022
     rectified    16.0         1    3.065e-01    1.707e+00          5.6     1.0000
     rectified    16.0         2    5.330e-01    2.059e+00          3.9     1.7387
     rectified    16.0         4    1.769e+00    3.761e+01         21.3     5.7719
     rectified    16.0         6    4.138e+01    1.813e+02          4.4   134.9930

The first four rows are the number of the vanishing gradient. With saturating activation the gradient magnitude seen by the first layer is 1.274e-02 with one layer and 5.539e-05 with six layers — by the last column 0.0043, that is, 230 times lower; the ratio to the last layer rises from 18.8 to 3,977.2. The last layer of the six-layer network sees a gradient of the same magnitude as that of the single-layer network; what disappears is only the share of the layers near the input. If this network is trained, the greater part of the 473 parameters stays where it was initialized.

In the middle four rows the same depths are measured with non-saturating activation and at six layers the ratio is 0.7022; there is almost no vanishing. What makes it vanish is not depth but the magnitude of the factor. The last four rows are the opposite extreme: when the gain is raised to 16 the factors become larger than one and the product runs upward. At six layers the gradient magnitude of the first layer is 134.99 times that of the single-layer one. A vanishing gradient freezes the weights, an exploding gradient throws them away in a single update.

What Clipping and the Schedule Correct

The explosion does not come only from initialization; a large learning rate produces the same result. The block below trains a network with three hidden layers first with a safe step and then with a step two hundred times larger, and tries the two corrections one by one and together.

def norm(G):                                # gradient magnitude of the whole network
    return math.sqrt(sum(v * v for W in G for w in W for v in w))


def dead_units(layers):                     # a unit that never activates on any training row
    alive = [[False] * len(W) for W in layers[:-1]]
    for x in X_TR:
        trace = forward(layers, x)
        for L in range(len(layers) - 1):
            for j, v in enumerate(trace[L + 1]):
                if v > 0:
                    alive[L][j] = True
    return sum(not c for W in alive for c in W)


def train_model(rate, seed, clip=0.0, schedule=False, depth=3, epochs=30, batch=20):
    layers, clipped, peak, curve = init_weights(seed, layer_shape(depth), 2.0), 0, 0.0, []
    for t in range(epochs):
        a = rate * 0.5 ** (t // 6) if schedule else rate
        s, total = order(t), 0.0
        for b in range(0, len(s), batch):
            G, k = gradient(layers, X_TR, Y_TR, s[b:b + batch])
            n = norm(G)
            peak = max(peak, n)
            if clip and n > clip:           # direction is kept, magnitude is pulled to the threshold
                clipped += 1
                for W in G:
                    for w in W:
                        w[:] = [v * clip / n for v in w]
            for W, gW in zip(layers, G):
                for w, g in zip(W, gW):
                    for j in range(len(w)):
                        w[j] -= a * g[j]
            total += k * batch
        curve.append(total / len(s))
    return layers, curve[-1], clipped, peak


SETTINGS = (("rate 0.05", 0.05, 0.0, False), ("rate 10 bare", 10.0, 0.0, False),
            ("rate 10 clipping", 10.0, 0.10, False), ("rate 10 schedule", 10.0, 0.0, True),
            ("rate 10 both", 10.0, 0.10, True))
print(f"{'setting':>16} {'loss 2 runs':>16} {'test 2 runs':>16} {'dead':>7} "
      f"{'clipped':>9} {'max_gradient':>14}")
OUTCOME = {}
for name, rate, clip, sched in SETTINGS:
    runs = [train_model(rate, SEED + s, clip, sched) for s in (1, 2)]
    OUTCOME[name] = [sum((forward(k, x)[-1][0] > 0.5) == (y > 0.5)
                         for x, y in zip(X_TE, Y_TE)) / len(X_TE) for k, _, _, _ in runs]
    print(f"{name:>16} {runs[0][1]:>7.4f} {runs[1][1]:>8.4f} {OUTCOME[name][0]:>8.4f} "
          f"{OUTCOME[name][1]:>7.4f} {dead_units(runs[0][0]):>6}/24 {runs[0][2]:>5}/450 "
          f"{runs[0][3]:>14.3e}")
print(f"\nbaseline {BASELINE:.4f}, network number {OUTCOME['rate 0.05'][0]:.4f}, "
      f"second run of the same network {OUTCOME['rate 0.05'][1]:.4f}")
for name in ("rate 10 bare", "rate 10 clipping", "rate 10 both"):
    print(f"{name:>16} run difference {abs(OUTCOME[name][0] - OUTCOME[name][1]):.4f} "
          f"(subscribers {round(abs(OUTCOME[name][0] - OUTCOME[name][1]) * len(TEST))})")

         setting      loss 2 runs      test 2 runs    dead   clipped   max_gradient
       rate 0.05  0.3666   0.3776   0.7817  0.7778      0/24     0/450      4.599e+00
    rate 10 bare  0.6210   0.6210   0.7579  0.7579     15/24     0/450      2.800e+07
rate 10 clipping  0.4333   0.3578   0.6984  0.6944      1/24   400/450      7.982e+00
rate 10 schedule  0.5287   0.5287   0.7579  0.7579     15/24     0/450      2.800e+07
    rate 10 both  0.2722   0.2473   0.7698  0.7183      0/24   422/450      6.445e+00

baseline 0.7579, network number 0.7817, second run of the same network 0.7778
    rate 10 bare run difference 0.0000 (subscribers 0)
rate 10 clipping run difference 0.0040 (subscribers 1)
    rate 10 both run difference 0.0516 (subscribers 13)

The second row is the explosion within the training. The gradient magnitude shoots up to 2.800e+07 at some point, that update throws the weights away and 15 of the 24 units never activate again on any training row. The result is the same in both runs: 0.7579, exactly the baseline. The network now says “not suspicious” to everyone. That the two initializations give the same number is not good news here — both have fallen to the same dead point.

The fourth row says what the schedule does by itself: nothing. The step is halved every six epochs but the explosion happens in the first epoch. A schedule does not undo damage that took place before it.

In the third row gradient clipping does its job: 400 of the 450 updates are clipped, the largest gradient it sees falls to 7.982 and the dead units go from 15 to 1. But the test number is 0.6984, clearly below the baseline. Clipping rescued the training, it did not rescue the model. In the fifth row the two together bring the loss to the lowest of the table, 0.2722; the test rises to 0.7698 but in the second run it falls to 0.7183 — 0.0516, thirteen subscribers. The widest run gap is in the setting that looks best.

The triple is as follows: baseline 0.7579, network number 0.7817, second run 0.7778. The shallow setting with the safe step is higher than all of the rescued deep settings, and none of them comes close to the 0.8016 of the linear model.

Summary

• With saturating activation the gradient magnitude of the first layer at six layers falls to 0.0043 times that of the single-layer one, that is, 230 times lower; the ratio to the last layer rises from 18.8 to 3,977.2.
• With non-saturating activation the ratio at the same depth stays at 0.7022; at gain 16 the same measure rises to 134.99 times and explodes. What makes it vanish is not depth but the magnitude of the factor.
• At a learning rate of 10 the gradient magnitude shoots up to 2.800e+07, 15 of the 24 units die and both runs give exactly the baseline; the schedule by itself does not correct this.
• Gradient clipping clips 400 of the 450 updates and brings the dead units down to 1 but the test falls to 0.6984; clipping and the schedule together produce the lowest loss and the widest run gap.

Next Step

In these two lessons the network went wrong in the same direction every time: as the training loss fell, the held-out set number tended to fall. With 300 rows and a few hundred parameters the name of this is known. The next lesson ties that name to a number. Weight decay was measured in Introduction to Machine Learning and Supervised Learning and will not be repeated; what will be added are two methods specific to the network — dropout, which randomly silences a part of the units during training, and early stopping, which cuts the training short. In the second the real question is not the method but the source: which set the stopping epoch is read from.