Backpropagation

In four lessons the gradient always arrived ready-made. How the number at the output layer is found was written in 04, but how that number is distributed over the weights of the hidden layer was never opened. This lesson writes that distribution and tests what it wrote with a separate method.

The rule of the distribution is the chain rule. The effect of a weight on the loss is the product of two parts: the value coming back to the unit it feeds and the input it multiplies. The value coming back is also a product — the sum of what comes from the next layer weighted by the weights, times the derivative of the activation. When a single rule is applied layer by layer, the gradients of all parameters come out.

• NN29. The set, the split, the seed, the twelve columns and the training set reduced to 300 rows are as in 01–04; the network is 12–8–1, 113 parameters. The baseline on the test set is 0.7579.
• NN30. The loss is cross-entropy. The value arriving at the output layer is $a-y$, because the derivative of the sigmoid cancels against the derivative of the loss; this was shown in 04.
• NN31. The activation of the hidden layer is the rectified linear unit. Its derivative is 1 when the input is positive and 0 otherwise; it is undefined at zero and 0 is chosen in the code.
• NN32. The gradients are summed over the examples. Because the numerical derivative is computed over the mean loss, its result is multiplied by the number of examples and the two methods come to the same scale.
• NN33. The numerical derivative is the central difference: the parameter is nudged once up and once down and the difference is divided by two steps. It is more accurate than the one-sided difference and asks for twice the passes.
• NN34. The verification is made on a subset of 20 rows; the reason is the number of passes.
• NN35. The number of passes is in terms of forward passes. The backward pass of backpropagation is of the same order as a forward pass and is not counted separately; the ratio is for this reason a lower bound.

The Chain Rule over the Network

Backpropagation consists of two runs. In the forward pass the intermediate values are stored, because the weight gradient needs them. In the backward pass one starts from the output and does two jobs at every layer: the gradients of that layer are written, then the value coming back is carried to the previous layer.

# backpropagation.py — MODEL: the same FICTIONAL subscriber table, the same seed and the same split;
# the training set is reduced to 300 rows as in 01.
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


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
    mean, vol, r = sum(v) / len(v), max(v) - min(v), rng(SEED + 51000 + k["id"])
    z = (0.052 * (mean - 20) + 0.85 * (min(v) == 0.0) + 0.026 * vol
         + 0.24 * (k["household"] >= 5) - 0.20 * (k["satisfaction"] >= 4) + k["region"][3]
         + 0.9 * (k["id"] - 10001) / 1399 + (r() + r() + r() - 1.5) * 1.30)
    x = {"household": k["household"], "satisfaction": k["satisfaction"], "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 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:]
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}
BASELINE = sum(x["suspicious"] == 0 for x in TEST) / len(TEST)


def make_set(data):
    return ([[(x[a] - MEAN[a]) / STD[a] for a in FIELDS] for x in data],
            [float(x["suspicious"]) for x in data])


X_TR, Y_TR = make_set(TRAIN[:300])          # REDUCED training set
X_VA, Y_VA = make_set(VAL)
X_TE, Y_TE = make_set(TEST)


def squash(z):                              # S-shaped function, overflow-protected
    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 rectified(z):
    return z if z > 0.0 else 0.0


def init_weights(seed, sizes):              # gain-scaled initialization, bias zero
    r, net = rng(seed), []
    for a, b in zip(sizes, sizes[1:]):
        scale = math.sqrt(2.0 / a)
        net.append(([[(r() * 2 - 1) * scale for _ in range(a)] for _ in range(b)], [0.0] * b))
    return net


def matmul(W, bias, a):                     # MATRIX MULTIPLICATION by hand: each row is a unit
    return [sum(w * v for w, v in zip(row, a)) + b
            for row, b in zip(W, bias)]


def forward(net, x):                        # FORWARD PASS: layer by layer
    a = x
    for k, (W, bias) in enumerate(net):
        z = matmul(W, bias, a)
        a = [squash(z[0])] if k == len(net) - 1 else [rectified(t) for t in z]
    return a


def loss(a, y):                             # cross-entropy
    return -(math.log(max(a, 1e-12)) if y else math.log(max(1 - a, 1e-12)))


def backpropagation(net, X, Y, trace=False):    # MODEL: chain rule layer by layer
    gW = [[[0.0] * len(W[0]) for _ in W] for W, _ in net]
    gb = [[0.0] * len(bias) for _, bias in net]
    total = 0.0
    for i in range(len(X)):
        a, acts = X[i], [X[i]]              # FORWARD PASS: intermediate values are stored
        for k, (W, bias) in enumerate(net):
            z = matmul(W, bias, a)
            a = [squash(z[0])] if k == len(net) - 1 else [rectified(q) for q in z]
            acts.append(a)
        total += loss(acts[-1][0], Y[i])
        d = [acts[-1][0] - Y[i]]            # output layer: loss and sigmoid cancel
        for k in range(len(net) - 1, -1, -1):
            W, a = net[k][0], acts[k]       # a: the INPUT of this layer
            if trace:
                print(f"  layer {k + 1}: incoming width {len(d)}, first two "
                      f"{[round(v, 5) for v in d[:2]]}")
            for j, dj in enumerate(d):      # weight gradient = incoming x input
                s = gW[k][j]
                for p in range(len(a)):
                    s[p] += dj * a[p]
                gb[k][j] += dj              # bias gradient = incoming
            if k:                           # carry to the previous layer, with the activation derivative
                prev = [0.0] * len(a)
                for j, dj in enumerate(d):
                    s = W[j]
                    for p in range(len(a)):
                        prev[p] += dj * s[p]
                d = [prev[p] if a[p] > 0 else 0.0 for p in range(len(a))]
    return gW, gb, total / len(X)


net = init_weights(SEED + 1, [12, 8, 1])
print(f"single subscriber, true label {int(Y_TR[0])}, network output {forward(net, X_TR[0])[0]:.5f}")
gW, gb, k0 = backpropagation(net, X_TR[:1], Y_TR[:1], trace=True)
print(f"  loss {k0:.5f}")
print(f"  output layer weight gradient, first three "
      f"{[round(v, 5) for v in gW[1][0][:3]]}")
print(f"  hidden layer unit 1 weight gradient, first three "
      f"{[round(v, 5) for v in gW[0][0][:3]]}")
print(f"  hidden layer bias gradients, first three {[round(v, 5) for v in gb[0][:3]]}")

single subscriber, true label 0, network output 0.49369
  layer 2: incoming width 1, first two [0.49369]
  layer 1: incoming width 8, first two [0.0, 0.00155]
  loss 0.68061
  output layer weight gradient, first three [0.0, 0.15433, 0.18243]
  hidden layer unit 1 weight gradient, first three [0.0, 0.0, 0.0]
  hidden layer bias gradients, first three [0.0, 0.00155, -0.03503]

The value arriving at the output layer is 0.49369, that is, the output of the network minus the true label. Eight values go down to the hidden layer and the first is exactly 0.0. The reason had been measured in 02: that unit is not active for this subscriber, the rectified function has cut its input, its derivative is 0 and the chain breaks there. The result shows in the last two lines — the twelve weight gradients of the first unit are zero and so is its bias gradient. This subscriber cannot teach that unit anything.

Verification with the Numerical Derivative

A derivative written by hand can be silently wrong: the code runs, the loss falls, numbers are produced, and the error shows itself only as performance staying a little low. The independent test is to measure the loss twice with a small nudge of the parameter and take the difference. It is slow, but because it knows nothing about the structure of the network, it does not make the same error as the chain rule.

def mean_loss(net, X, Y):
    return sum(loss(forward(net, x)[0], y) for x, y in zip(X, Y)) / len(X)


def numerical_derivative(net, X, Y, eps):   # central difference, two passes per parameter
    gW = [[[0.0] * len(W[0]) for _ in W] for W, _ in net]
    gb = [[0.0] * len(bias) for _, bias in net]
    passes = 0
    for k, (W, bias) in enumerate(net):
        for j in range(len(W)):
            for p in range(len(W[j]) + 1):
                target, index = (W[j], p) if p < len(W[j]) else (bias, j)
                old = target[index]
                target[index] = old + eps
                plus = mean_loss(net, X, Y)
                target[index] = old - eps
                minus = mean_loss(net, X, Y)
                target[index] = old
                value = (plus - minus) / (2 * eps) * len(X)
                passes += 2 * len(X)
                if p < len(W[j]):
                    gW[k][j][p] = value
                else:
                    gb[k][j] = value
    return gW, gb, passes


def check(net, X, Y, eps):                  # hand-written gradient against the numerical derivative
    aW, ab, _ = backpropagation(net, X, Y)
    nW, nb, passes = numerical_derivative(net, X, Y, eps)
    diffs = []
    for k in range(len(net)):
        for j in range(len(aW[k])):
            for p in range(len(aW[k][j])):
                diffs.append((abs(aW[k][j][p] - nW[k][j][p]), abs(aW[k][j][p])))
            diffs.append((abs(ab[k][j] - nb[k][j]), abs(ab[k][j])))
    return len(diffs), max(a for a, _ in diffs), max(a / max(b, 1e-12) for a, b in diffs), passes


def train_model(net, X, Y, rate, epochs):   # MODEL: full-batch gradient descent
    n = len(X)
    for _ in range(epochs):
        gW, gb, _ = backpropagation(net, X, Y)
        for k, (W, bias) in enumerate(net):
            for j in range(len(W)):
                s, g = W[j], gW[k][j]
                for p in range(len(s)):
                    s[p] -= rate * g[p] / n
                bias[j] -= rate * gb[k][j] / n
    return net


SUB_X, SUB_Y = X_TR[:20], Y_TR[:20]         # verification subset, 20 rows
TRAINED = train_model(init_weights(SEED + 1, [12, 8, 1]), X_TR, Y_TR, 0.3, 100)
print(f"{'network':<24} {'eps':>7} {'params':>10} {'max absolute':>16} "
      f"{'max relative':>15}")
for label, net in (("at initialization", init_weights(SEED + 1, [12, 8, 1])),
                   ("after 100 epochs", TRAINED)):
    for eps in (1e-3, 1e-5, 1e-7):
        count, max_abs, max_rel, passes = check(net, SUB_X, SUB_Y, eps)
        print(f"{label:<24} {eps:>7g} {count:>10} {max_abs:>16.3e} {max_rel:>15.3e}")

network                      eps     params     max absolute    max relative
at initialization          0.001        113        2.199e-07       1.583e-07
at initialization          1e-05        113        2.193e-10       4.567e-08
at initialization          1e-07        113        2.591e-08       4.932e-06
after 100 epochs           0.001        113        2.518e-01       2.881e+00
after 100 epochs           1e-05        113        6.441e-11       1.227e-09
after 100 epochs           1e-07        113        5.694e-09       2.697e-07

The first three rows show the error of the numerical derivative itself. When the step is large the truncation error dominates, when it is small the rounding error; the best agreement is in the middle, at 1e-05, with 2.193e-10. That is, across all 113 parameters the two methods give the same number up to about ten significant digits. This is the proof that the hand-written chain rule is correct.

The fourth row is a trap and must not be passed over without measuring. In the trained network, when the step is chosen as 0.001, the largest absolute difference jumps to 2.518e-01 and the relative difference to 2.881. This is not a calculation error; the same network at the same step gives 6.441e-11 at 1e-05. The cause is the kink of the rectified unit: a nudge of 0.001 is enough to carry the weighted sum from one side of zero to the other, and the two ends of the central difference fall on two different linear pieces of the function. When a verification fails, the first place to look is not the code but the step size.

Cost: How Many Passes

If the two methods produce the same number, which one to use depends on the cost. Backpropagation produces the gradients of all parameters in one forward and one backward pass. The numerical derivative asks for two separate passes for each parameter.

print(f"{'units':>6} {'params':>10} {'backpropagation':>15} {'numerical derivative':>20} "
      f"{'ratio':>7}")
for h in (2, 8, 32):
    sizes = [12, h, 1]
    P = sum(a * b + b for a, b in zip(sizes, sizes[1:]))
    print(f"{h:>6} {P:>10} {len(SUB_X):>15} {2 * P * len(SUB_X):>20} {2 * P:>7}")

print(f"\n{'run':>6} {'train_loss':>10} {'validation':>10} {'test':>8}")
TEST_ACCS = []
for run in (1, 2):
    net = train_model(init_weights(SEED + run, [12, 8, 1]), X_TR, Y_TR, 0.3, 100)
    _, _, loss_value = backpropagation(net, X_TR, Y_TR)
    val_acc = sum((forward(net, x)[0] > 0.5) == (y > 0.5) for x, y in zip(X_VA, Y_VA)) / len(X_VA)
    test_acc = sum((forward(net, x)[0] > 0.5) == (y > 0.5) for x, y in zip(X_TE, Y_TE)) / len(X_TE)
    TEST_ACCS.append(test_acc)
    print(f"{run:>6} {loss_value:>10.4f} {val_acc:>10.4f} {test_acc:>8.4f}")
print(f"baseline {BASELINE:.4f}, network number {TEST_ACCS[0]:.4f}, "
      f"second run {TEST_ACCS[1]:.4f}, difference "
      f"{round(abs(TEST_ACCS[0] - TEST_ACCS[1]) * len(X_TE))} subscribers")

 units     params backpropagation numerical derivative   ratio
     2         29              20                 1160      58
     8        113              20                 4520     226
    32        449              20                17960     898

   run train_loss validation     test
     1     0.3597     0.8214   0.7857
     2     0.3646     0.8056   0.7698
baseline 0.7579, network number 0.7857, second run 0.7698, difference 4 subscribers

The backpropagation column is 20 in all three rows — as many as the number of examples. The numerical derivative column grows together with the parameter count: 1,160, 4,520, 17,960. The ratio is twice the parameter count and rises to 898 in the network with 32 units. The cost of backpropagation is independent of the parameter count; the work is done by the intermediate values stored in the forward pass. For this reason the numerical derivative is not a training method but a tool of inspection.

The last row confirms the number of the four lessons. In this section the training was carried out with backpropagation and the triple is the same output as in 01: 0.7579 — 0.7857 — 0.7698. The reason is plain; the code running since 01 was already backpropagation, only it had not been named.

Summary

• The gradient of a weight is the product of the value coming back to the unit it feeds and the input that weight multiplies; the value coming back is carried from the next layer and multiplied by the derivative of the activation.
• All gradients of an inactive unit are zero: in the example subscriber, the twelve weight gradients and the bias gradient of the first unit of the hidden layer all come out as 0.0.
• Verification with the central difference agrees in all 113 parameters: the largest absolute difference is 2.193e-10 at initialization and 6.441e-11 after a hundred epochs of training. When the step is chosen too large or too small, the error of the numerical derivative itself grows.
• In the trained network the difference jumps to 2.518e-01 at step 0.001; this is the kink of the rectified unit, not a calculation error. The same network gives 6.441e-11 at 1e-05.
• Backpropagation produces all 113 gradients with one forward and one backward pass per example; the numerical derivative asks for 4,520 passes for the same job and the ratio grows with the parameter count: 58, 226, 898.

Next Step

The gradient is now computed correctly and cheaply. But what to do with it was never asked in five lessons: in every epoch the gradient was multiplied by the same constant number and subtracted from the weights, always of the same size, without ever looking at what the previous epochs had done. Using the same gradient in another way — remembering the direction of the previous steps or giving each parameter its own step size — what loss it reaches in the same epoch budget and how much it moves between two runs was not measured. The next lesson treats the gradient not as a result but as an input.