Computation Graph and Automatic Differentiation

In the previous ten lessons the architecture was chosen: convolution, pooling, skip connection, recurrent state, gated unit, attention, embedding, transfer, autoencoder. In each one the structure of the network changed and the gradient computation was re-derived and rewritten; all of them ran in a small hand-written training loop. The machinery of real training was never addressed.

This topic addresses that machinery, and its first piece is the abstraction that separates gradient computation from the architecture. In the catalog this lesson mentions deep learning libraries; what is described here is not libraries but the abstraction libraries sell — the computation graph and automatic differentiation — and the abstraction is written within the lesson. The measure is three numbers: the difference between the gradient given by the hand-written backpropagation and by automatic differentiation, how many lines the two routes take, and how many nodes the graph builds.

• TP1. The dataset is the fictional subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering, Supervised Learning and Unsupervised Learning courses, seed 20260218, 1,260 subscribers; the label is suspected illicit use or a fault and is fictional. The baseline is saying the most frequent class: 0.7579 on the test set, inherited from the Introduction to Machine Learning, Data Preparation and Feature Engineering and Supervised Learning courses.
• TP2. The training set in this topic is 300 rows instead of 756; it is a scale decision, and every block of the topic runs in a few seconds. The test set is 252 rows and one subscriber is worth 0.0040.
• TP3. The network is fixed throughout this topic: eight features, a single hidden layer of eight units (rectified linear unit), a single output; the loss is cross-entropy. The constant term of the layer is called the bias and is distinct from the bias in statistics. Scaling is learned from the training set only (the Data Preparation and Feature Engineering course, pipeline).
• TP4. The third number of the course is the second run: the same network, the same data, the same example order; only the initial values change. The measure is not raw time: the cost of the graph is measured in nodes, the cost of the code in lines.
• TP5. The loss node is fused: sigmoid and cross-entropy are held in a single node, because held separately they overflow on large outputs. Its local derivative is the same on both routes. The machinery is a model, and no framework, library or vendor name is mentioned.

The Dataset, Scaling and the Baseline

# COMMON — FICTIONAL. The subscriber table of the Introduction to Machine Learning through
# Unsupervised Learning courses is generated with the same seed (20260218);
# the label is suspected illicit use or a fault. The training set is REDUCED in this topic
# from 756 to 300 rows: every block of the topic runs in a few seconds.
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)]
FIELDS = ["area_south", "area_east", "area_west", "area_central",
          "household", "satisfaction", "period", "mean"]


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


def mix(items, seed):
    r, s = rng(seed), list(items)
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return s


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
    subscriber_id = 10001 + i
    r, v = rng(SEED + 7000 + subscriber_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(
            b[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 + subscriber_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 * (subscriber_id - 10001) / 1399 + (r() + r() + r() - 1.5) * 1.30)
    x = {"household": household, "satisfaction": rating, "period": len(v),
         "mean": round(mean, 2), "suspicious": int(z > 1.35)}
    for q in REGION[1:]:
        x["area_" + q[0]] = int(b[0] == q[0])
    DATA.append(x)

K = mix(DATA, SEED + 90000)
TRAIN, TEST = K[:300], K[1008:]


def table(rows, stats=None):
    X = [[float(s[a]) for a in FIELDS] for s in rows]
    if stats is None:                     # scaling is learned from the training set only
        stats = []
        for j in range(4, 8):
            m = sum(v[j] for v in X) / len(X)
            stats.append((j, m, math.sqrt(sum((v[j] - m) ** 2 for v in X) / len(X))))
    for v in X:
        for j, m, sd in stats:
            v[j] = (v[j] - m) / sd
    return X, [s["suspicious"] for s in rows], stats


X_train, y_train, STATS = table(TRAIN)
X_test, y_test, _ = table(TEST, STATS)
BASELINE = sum(v == 0 for v in y_test) / len(y_test)
print(f"fictional subscribers {len(DATA)}; train {len(X_train)} (reduced from 756), test "
      f"{len(X_test)}, features {len(FIELDS)}, suspicious share in train "
      f"{sum(y_train) / len(y_train):.4f}")
print(f"baseline (predicting the most frequent class, on the test set) {BASELINE:.4f}")

fictional subscribers 1260; train 300 (reduced from 756), test 252, features 8, suspicious share in train 0.2200
baseline (predicting the most frequent class, on the test set) 0.7579

A Graph Is a Data Structure

The trouble with deriving by hand is not that the computation is hard but that it is tied to the structure: when a layer is added, the derivative is rewritten. The computation graph breaks that tie. Every intermediate value is kept as a node, and a node stores its own value, its parents and its local derivative with respect to its parents; the backward pass walks the graph in reverse and distributes the gradient of each node to its parents. The chain rule is this distribution, and it is written only once.

# Computation graph machinery — MODEL. Four operations: addition, multiplication, rectification
# and fused loss. Each node holds its value, its gradient, its parents and its LOCAL DERIVATIVE
# WITH RESPECT TO ITS PARENTS; the backward step uses only these.
import inspect


class Value:
    __slots__ = ("v", "grad", "parents", "local")

    def __init__(self, v, parents=(), local=()):
        self.v, self.grad, self.parents, self.local = float(v), 0.0, parents, local

    def __add__(self, o):
        o = o if isinstance(o, Value) else Value(o)
        return Value(self.v + o.v, (self, o), (1.0, 1.0))

    def __mul__(self, o):
        o = o if isinstance(o, Value) else Value(o)
        return Value(self.v * o.v, (self, o), (o.v, self.v))

    def rectify(self):
        return Value(self.v if self.v > 0 else 0.0, (self,),
                     (1.0 if self.v > 0 else 0.0,))

    def loss(self, y):                   # cross-entropy, directly on the output
        p = 1.0 / (1.0 + math.exp(-self.v)) if self.v > -30 else 0.0
        k = max(0.0, self.v) + math.log(1.0 + math.exp(-abs(self.v))) - self.v * y
        return Value(k, (self,), (p - y,))


def backward(root):
    order, seen = [], set()

    def visit(d):
        if id(d) in seen:
            return
        seen.add(id(d))
        for a in d.parents:
            visit(a)
        order.append(d)

    visit(root)
    root.grad = 1.0
    for d in reversed(order):
        for a, y in zip(d.parents, d.local):
            a.grad += d.grad * y
    return len(order)


LINES = sum(len([s for s in inspect.getsource(o).splitlines() if s.strip()])
            for o in (Value, backward))
print(f"machinery {LINES} lines, 4 operations;  L = 5 * rectify(a*b + c)")
for b_value in (3.0, -3.0):
    a, b, c = Value(2.0), Value(b_value), Value(4.0)
    L = (a * b + c).rectify() * Value(5.0)
    n = backward(L)
    print(f"  a=2 b={b_value:+.0f} c=4 -> L={L.v:5.1f}, {n} nodes;  dL/da {a.grad:>6.1f}"
          f"   dL/db {b.grad:>6.1f}   dL/dc {c.grad:>6.1f}")

machinery 32 lines, 4 operations;  L = 5 * rectify(a*b + c)
  a=2 b=+3 c=4 -> L= 50.0, 8 nodes;  dL/da   15.0   dL/db   10.0   dL/dc    5.0
  a=2 b=-3 c=4 -> L=  0.0, 8 nodes;  dL/da    0.0   dL/db    0.0   dL/dc    0.0

Thirty-two lines and four operations. The derivatives can be checked by hand: $L = 5\,(ab+c)$ gives $\partial L/\partial a = 5b = 15$, $\partial L/\partial b = 5a = 10$ and $\partial L/\partial c = 5$ in the positive region; the output says so. In the second call only $b$ became negative three: the sum went negative, the rectifier gate closed, and all three derivatives became zero. The graph is the same, the code is the same; what changed is only the path the data opens in the graph.

The Same Gradient, Two Routes

# The gradient of the same network by two routes: walking the graph and writing the chain
# rule BY HAND. The parameters are one flat list: p[j*G+i] weight, p[H*G+j] bias, p[H*G+H+j]
# and p[-1] output layer.
G, H = len(FIELDS), 8


def init_weights(seed):                      # gain-scaled initialization
    r = rng(seed)
    return [(r() * 2 - 1) * math.sqrt(2.0 / G) for _ in range(H * G)] + [0.0] * H \
        + [(r() * 2 - 1) * math.sqrt(2.0 / H) for _ in range(H)] + [0.0]


def by_graph(p, X, y, idx):
    nodes, total = [Value(v) for v in p], None
    for k in idx:
        x, o = X[k], nodes[-1]
        for j in range(H):
            z = nodes[H * G + j]
            for i in range(G):
                z = z + nodes[j * G + i] * x[i]
            o = o + nodes[H * G + H + j] * z.rectify()
        L = o.loss(y[k])
        total = L if total is None else total + L
    root = total * (1.0 / len(idx))
    return root, backward(root), [u.grad for u in nodes]


def by_hand(p, X, y, idx):
    g, loss_sum = [0.0] * len(p), 0.0
    for k in idx:
        x, hh = X[k], [0.0] * H
        for j in range(H):
            s = p[H * G + j] + sum(p[j * G + i] * x[i] for i in range(G))
            hh[j] = s if s > 0 else 0.0
        o = p[-1] + sum(p[H * G + H + j] * hh[j] for j in range(H))
        loss_sum += max(0.0, o) + math.log(1.0 + math.exp(-abs(o))) - o * y[k]
        d = 1.0 / (1.0 + math.exp(-o)) - y[k]    # local derivative of the loss node
        g[-1] += d
        for j in range(H):
            g[H * G + H + j] += d * hh[j]
            if hh[j] > 0:                        # a closed rectifier gate does not pass it
                g[H * G + j] += d * p[H * G + H + j]
                for i in range(G):
                    g[j * G + i] += d * p[H * G + H + j] * x[i]
    return loss_sum / len(idx), [v / len(idx) for v in g]


P0, BATCH = init_weights(SEED + 1), list(range(30))
root, NODES, auto = by_graph(P0, X_train, y_train, BATCH)
loss_value, manual = by_hand(P0, X_train, y_train, BATCH)
print(f"parameters {len(P0)}; for a 30-example batch the graph has {NODES} nodes, "
      f"{NODES / 30:.1f} per example")
print(f"loss  graph {root.v:.10f}   by hand {loss_value:.10f}")
print(f"largest absolute gradient difference {max(abs(u - v) for u, v in zip(auto, manual)):.2e}"
      f"  (largest gradient {max(abs(v) for v in auto):.4f})")
print(f"hand-written route {len(inspect.getsource(by_hand).splitlines())} lines "
      f"and belongs to this network only; machinery {LINES} lines and knows nothing of the network")

parameters 81; for a 30-example batch the graph has 6622 nodes, 220.7 per example
loss  graph 0.7366651135   by hand 0.7366651135
largest absolute gradient difference 2.78e-17  (largest gradient 0.3097)
hand-written route 18 lines and belongs to this network only; machinery 32 lines and knows nothing of the network

This is the lesson’s first number. The gradient of eighty-one parameters was computed by two routes and the largest absolute difference came out as 2.78e-17: a difference produced by the order of floating-point addition. Since the largest gradient is 0.3097, the two routes give the same number to fifteen digits; equality is checked not with == but with the ratio of the difference to the largest gradient.

The second number says the opposite. The hand-written route took eighteen lines, and these lines belong to this network only; if a layer is added they are rewritten. The thirty-two-line machinery knows nothing of the network’s structure; whoever writes the forward pass, the same backward call gives the gradient. Its cost is in the third number: for a batch of thirty examples the graph built 6,622 nodes, 220.7 per example. The hand-written route holds no nodes; it accumulates the intermediate values in place.

The Course’s Three Numbers

# Training runs with the HAND-WRITTEN backpropagation; the graph was built for verification. The two
# runs differ ONLY IN THE INITIAL VALUES: the example order is the same in both.
STEP, EPOCHS, BATCH_SIZE, ORDER_SEED = 0.30, 40, 30, SEED + 700


def train_model(seed):
    p = init_weights(seed)
    for t in range(EPOCHS):
        indices = mix(range(len(X_train)), ORDER_SEED + t)
        for start in range(0, len(indices), BATCH_SIZE):
            _, g = by_hand(p, X_train, y_train, indices[start:start + BATCH_SIZE])
            p = [a - STEP * b for a, b in zip(p, g)]
    return p


def accuracy(p, X, y):
    n = 0
    for k in range(len(X)):
        x = X[k]
        o = p[-1] + sum(p[H * G + H + j] * max(0.0, p[H * G + j]
                        + sum(p[j * G + i] * x[i] for i in range(G)))
                        for j in range(H))
        n += (o >= 0.0) == (y[k] == 1)
    return n / len(X)


S = [accuracy(train_model(SEED + t), X_test, y_test) for t in range(1, 11)]
print(f"network {G}-{H}-1, {len(P0)} parameters, {EPOCHS} epochs, batch size {BATCH_SIZE}, "
      f"learning rate {STEP}, {EPOCHS * (len(X_train) // BATCH_SIZE)} updates")
print(f"  baseline          {BASELINE:.4f}")
print(f"  network's number  {S[0]:.4f}   (difference {S[0] - BASELINE:+.4f})")
print(f"  second run        {S[1]:.4f}   (difference {S[1] - BASELINE:+.4f})")
print(f"range over ten initializations {min(S):.4f} - {max(S):.4f}, width "
      f"{max(S) - min(S):.4f} ({round((max(S) - min(S)) * len(X_test))} subscribers)")
TREE, LINEAR = 198, 202         # the Supervised Learning course's numbers, on the same 252-row test set
low, high = round(min(S) * len(X_test)), round(max(S) * len(X_test))
print(f"comparison lines (with the 756-row training set): decision tree {TREE / len(X_test):.4f}, "
      f"linear model {LINEAR / len(X_test):.4f}; the network's range {low}-{high} subscribers"
      f" -> {sum(low <= v <= high for v in (TREE, LINEAR))}/2 inside")
print(f"had the same training run through the graph, "
      f"{EPOCHS * (len(X_train) // BATCH_SIZE) * NODES:,} nodes would have been built")

network 8-8-1, 81 parameters, 40 epochs, batch size 30, learning rate 0.3, 400 updates
  baseline          0.7579
  network's number  0.8095   (difference +0.0516)
  second run        0.7937   (difference +0.0357)
range over ten initializations 0.7857 - 0.8095, width 0.0238 (6 subscribers)
comparison lines (with the 756-row training set): decision tree 0.7857, linear model 0.8016; the network's range 198-204 subscribers -> 2/2 inside
had the same training run through the graph, 2,648,800 nodes would have been built

Three numbers side by side: the baseline 0.7579, the network’s number 0.8095, the second run 0.7937. The difference between them is 0.0158, that is, four subscribers; this difference between the two runs is close to a third of the gain over the baseline (+0.0516). Over ten initializations the range spreads between 0.7857 and 0.8095, six subscribers wide. With nothing changing in capacity, data, learning rate or number of epochs, the network’s number can land anywhere in this range.

The Supervised Learning course’s decision tree got 198 and its linear model 202 subscribers right on the test set; the network’s range over ten runs is 198–204 and both are inside. The comparison is not exact, since those numbers came from the 756-row training set; but the direction is clear: saying that the network passed the two older families with the number it produced in a single run is reading run-to-run variability as a finding. This is why the course’s rule exists — a result that does not state that the same capacity gave two different numbers in two runs counts as unmeasured.

Summary

• The computation graph is a data structure: each node holds its value, its parents and its local derivative; the backward pass walks the graph in reverse and the chain rule is written only once.
• The thirty-two-line machinery with four operations gives the gradient of eighty-one parameters with a difference of 2.78e-17 from the eighteen-line hand-written backpropagation; the largest gradient is 0.3097.
• Thirty-two lines do not know the structure of the network, eighteen lines belong to this network only; the cost is 6,622 nodes for a batch of thirty examples, 220.7 per example.
• The course’s three numbers: baseline 0.7579, the network’s number 0.8095, second run 0.7937; the difference between them is 0.0158, close to a third of the gain over the baseline.
• The range of ten initializations is 0.7857–0.8095, and the numbers of the Supervised Learning course’s decision tree (0.7857) and linear model (0.8016) fall within this range; an advantage read from a single run cannot be separated from run-to-run variability.

Next Step

The gradient was computed correctly, but which examples were fed to the network in which order passed in a single line: the training set was shuffled every epoch and split into batches of thirty. That line is a decision. The next lesson addresses the data loading pipeline: how much the batch generation order moves the number, what data augmentation adds, and in which case augmentation produces leakage into the held-out set are measured.