Optimization Algorithms

The previous lesson showed that the gradient is computed correctly: the number given by backpropagation was compared with the numerical derivative and the difference was printed. But the gradient gives a direction, not a step. In five lessons it was never asked by which rule the weights are moved in that direction; because it was never asked, there seems to be a single rule: subtract a constant multiple of the gradient from the weight.

This lesson fills that gap. With the same network, the same data, the same gradient and the same epoch budget, three different step rules are run. The measure is twofold: the job of the optimizer is to lower the training loss, the job of the model is to be right on the held-out set. The finding of this lesson is that the two can go in opposite directions.

• NN46. The set, the split and the seed are as in the previous lessons; the table is fictional. So that the run time stays within the lesson, the training set has been reduced from 756 rows to 300 rows; the validation and test sets did not change (252 and 252) and the baseline on the test set is 0.7579.
• NN47. The network is fixed: 12 inputs, 8 hidden units, 1 output. The only thing that changes is the step rule. There is the rectified linear unit in the hidden layer, a sigmoid at the output and cross-entropy as the loss; the initialization is gain-scaled and the biases start from zero.
• NN48. The epoch budget is 40 and the batch size is 20; all three spend the same number of updates.
• NN49. The three rules are: plain descent (a constant multiple of the gradient is subtracted from the weight), momentum descent (gradients accumulate with a decay factor of 0.9) and the adaptive step (a squared-gradient accumulation is kept for each weight and the step is divided by the root of that accumulation).
• NN50. The batch split order is the same in every run. The only thing that changes from run to run is the initialization.
• NN51. The number of candidates tried is 9: three optimizers times three learning rates, with the choice from the validation set.
• NN52. The comparison lines come from the Supervised Learning course and are not retrained: the linear model 0.8016, the decision tree 0.7857. Both had been measured on the same test set.
• NN53. The measure of speed is not time but the number of epochs: in which epoch the loss first falls below the 0.40 threshold. The threshold is arbitrary and the same for the three optimizers.

From Gradient to Step

The setup produces the fictional table, separates the sets and defines the network. The last coefficient of every weight row is the bias.

# optimization.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 forward(layers, x):                     # rectified linear unit in the hidden layers,
    a, trace = x, []                        # sigmoid at the output
    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]
        a = [squash(z[0])] if i == len(layers) - 1 else [t if t > 0 else 0.0 for t in z]
    trace.append(a)
    return trace


def gradient(layers, X, Y, idx):            # BACKPROPAGATION, cross-entropy loss
    G, loss = zeros(layers), 0.0
    for i in idx:
        trace = forward(layers, X[i])
        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))) * (a[j] > 0)
                     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 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), 8), (8, 1)])
print(f"train {len(TRAIN)} (reduced from 756), validation {len(VAL)}, "
      f"test {len(TEST)}, columns {len(FIELDS)}")
print(f"network {len(FIELDS)}-8-1, parameters {sum(len(w) for W in LAYERS0 for w in W)}, "
      f"epochs 40, batch 20, updates {40 * len(X_TR) // 20}")
print(f"baseline {BASELINE:.4f}; comparison lines linear model {LINEAR:.4f}, "
      f"decision tree {TREE:.4f}; initialized network {accuracy(LAYERS0, X_TE, Y_TE):.4f}")

train 300 (reduced from 756), validation 252, test 252, columns 12
network 12-8-1, parameters 113, epochs 40, batch 20, updates 600
baseline 0.7579; comparison lines linear model 0.8016, decision tree 0.7857; initialized network 0.3333

The last number gives the starting point: the untrained network scores 0.3333 on the test set, far below the baseline. All 113 parameters are still noise, and the job of the 600 updates is to carry this number upward.

Three Optimizers in the Same Epoch Budget

The three rules sit in a single function. Plain descent uses the gradient directly. Momentum descent accumulates past gradients; if the direction is consistent the step grows, if the direction oscillates the steps cancel each other. The adaptive step keeps a separate squared-gradient accumulation for each weight and divides the step of that weight by the root of the accumulation; a weight that keeps seeing large gradients moves with a small step, a weight that sees rare gradients with a large one.

def apply_step(layers, G, M, V, method, rate):  # a single update step
    for W, gW, mW, vW in zip(layers, G, M, V):
        for w, g, m, v in zip(W, gW, mW, vW):
            for j in range(len(w)):
                if method == "momentum":            # accumulated direction
                    m[j] = 0.9 * m[j] + g[j]
                    w[j] -= rate * m[j]
                elif method == "adaptive":          # squared-gradient accumulation
                    v[j] = 0.9 * v[j] + 0.1 * g[j] * g[j]
                    w[j] -= rate * g[j] / (math.sqrt(v[j]) + 1e-8)
                else:                               # plain descent
                    w[j] -= rate * g[j]


def train_model(method, rate, seed, epochs=40, batch=20):
    layers = init_weights(seed, [(len(FIELDS), 8), (8, 1)])
    M, V, curve = zeros(layers), zeros(layers), []
    for t in range(epochs):
        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])
            apply_step(layers, G, M, V, method, rate)
            total += k * batch
        curve.append(total / len(s))
    return layers, curve


def epochs_to_target(curve, threshold=0.40):    # first epoch the loss falls below the threshold
    for t, k in enumerate(curve):
        if k < threshold:
            return str(t + 1)
    return "never"


METHODS = ("plain", "momentum", "adaptive")
RESULTS = {}
print(f"{'optimizer':>12} {'rate':>6} {'train_loss':>13} {'epochs_to_target':>16} "
      f"{'validation':>10} {'test':>8}")
for m in METHODS:
    for rate in (0.02, 0.05, 0.15):
        layers, curve = train_model(m, rate, SEED + 1)
        RESULTS[(m, rate)] = (curve[-1], accuracy(layers, X_VA, Y_VA), accuracy(layers, X_TE, Y_TE))
        print(f"{m:>12} {rate:>6} {curve[-1]:>13.4f} {epochs_to_target(curve):>16} "
              f"{RESULTS[(m, rate)][1]:>10.4f} {RESULTS[(m, rate)][2]:>8.4f}")
CHOSEN = {m: max((0.02, 0.05, 0.15), key=lambda a: RESULTS[(m, a)][1]) for m in METHODS}
lowest = min(RESULTS.values(), key=lambda t: t[0])
highest = max(RESULTS.values(), key=lambda t: t[2])
print(f"\ncandidates tried {len(RESULTS)}; rates read from validation "
      + ", ".join(f"{m} {CHOSEN[m]}" for m in METHODS)
      + f"\nlowest loss {lowest[0]:.4f} (test {lowest[2]:.4f}), "
      + f"highest test {highest[2]:.4f} (loss {highest[0]:.4f})")

   optimizer   rate    train_loss epochs_to_target validation     test
       plain   0.02        0.4042            never     0.8135   0.7857
       plain   0.05        0.3664               19     0.8056   0.7897
       plain   0.15        0.3296                8     0.7817   0.7698
    momentum   0.02        0.3186                6     0.7976   0.7579
    momentum   0.05        0.2811                4     0.7817   0.7540
    momentum   0.15        0.2777                3     0.7817   0.7500
    adaptive   0.02        0.2948                4     0.7659   0.7341
    adaptive   0.05        0.2923                3     0.7500   0.7222
    adaptive   0.15        0.3038                8     0.7341   0.7460

candidates tried 9; rates read from validation plain 0.02, momentum 0.02, adaptive 0.02
lowest loss 0.2777 (test 0.7500), highest test 0.7897 (loss 0.3664)

The loss column shows the optimizers doing their job. In the same 600 updates momentum descent goes down to 0.2777 and the adaptive step to 0.2923; the best of plain descent stays at 0.3296. The difference is sharper in epochs: momentum descent crosses the 0.40 threshold in the third epoch, plain descent at the smallest rate does not cross it at all in forty epochs.

The test column does not give the counterpart. The setting that produces the lowest loss gives 0.7500 on the test set; the highest test number, 0.7897, is produced by plain descent, which brings the loss down only to 0.3664. The two columns go in opposite directions, because with 113 parameters on 300 rows there is a way to pull the loss lower: memorizing the training rows. Momentum descent and the adaptive step get there faster. The optimizer lowers the loss, it does not lower the generalization error. For this reason the choice is read from the validation set, and in all three the smallest rate, 0.02, wins.

Run Variability Swallows the Method Difference

Each row of the table above is a single run. The block below runs the three optimizers four times each with the rates read from validation; the data, the order, the budget and the rate are fixed, only the initialization changes.

print(f"{'optimizer':>12} {'loss 4 runs':>32} {'test 4 runs':>32}")
SPREAD = {}
for m in METHODS:
    runs = [train_model(m, CHOSEN[m], SEED + t) for t in (1, 2, 3, 4)]
    SPREAD[m] = ([curve[-1] for _, curve in runs], [accuracy(l, X_TE, Y_TE) for l, _ in runs])
    print(f"{m:>12} " + " ".join(f"{v:>7.4f}" for v in SPREAD[m][0])
          + "  " + " ".join(f"{v:>7.4f}" for v in SPREAD[m][1]))
for m in METHODS:
    s = SPREAD[m][1]
    print(f"{m:>12} test spread {max(s) - min(s):.4f} "
          f"({round((max(s) - min(s)) * len(TEST))} subscribers)")
best = max(METHODS, key=lambda m: SPREAD[m][1][0])
print(f"\nbaseline {BASELINE:.4f}, network number {SPREAD[best][1][0]:.4f} ({best}), "
      f"second run of the same network {SPREAD[best][1][1]:.4f}")
print("difference between the best and the worst optimizer per run: "
      + ", ".join(f"{max(SPREAD[m][1][t] for m in METHODS) - min(SPREAD[m][1][t] for m in METHODS):.4f}"
                  for t in range(4)))
print(f"runs beating the decision tree "
      f"{sum(v > TREE + 1e-9 for _, s in SPREAD.values() for v in s)}/12, the linear model "
      f"{sum(v > LINEAR + 1e-9 for _, s in SPREAD.values() for v in s)}/12")

   optimizer                      loss 4 runs                      test 4 runs
       plain  0.4042  0.3874  0.3980  0.3972   0.7857  0.7698  0.7897  0.7817
    momentum  0.3186  0.3373  0.3244  0.3177   0.7579  0.7659  0.7817  0.7738
    adaptive  0.2948  0.3064  0.2922  0.2976   0.7341  0.7619  0.7619  0.7659
       plain test spread 0.0198 (5 subscribers)
    momentum test spread 0.0238 (6 subscribers)
    adaptive test spread 0.0317 (8 subscribers)

baseline 0.7579, network number 0.7857 (plain), second run of the same network 0.7698
difference between the best and the worst optimizer per run: 0.0516, 0.0079, 0.0278, 0.0159
runs beating the decision tree 1/12, the linear model 0/12

The loss columns are narrow: each optimizer stays in its own band across four runs and the three bands do not touch each other. The effect of the optimizer on the loss is stable at this scale.

The test columns are not narrow. The adaptive step wanders between 0.7341 and 0.7659; this is 0.0317, that is, eight subscribers. The line before the last is the real finding of the lesson: the difference between the optimizers is 0.0516 in the first run and 0.0079 in the second. The same three methods, the same settings, the same data — only the initialization changed and the distance between them narrowed by more than a factor of six. To look at the first run and say “plain descent beats the adaptive step by 0.0516” is not a measured result.

The triple is as follows: baseline 0.7579, network number 0.7857, the second run of the same network 0.7698; the difference from the baseline is 0.0278, that is, seven subscribers. The last line gives the second claim of the course: only one of the twelve runs passes the 0.7857 of the decision tree and none passes the 0.8016 of the linear model. On this tabular data the network does not pass two simple families, and this is not hidden.

Summary

• Backpropagation gives a direction; the step rule is a separate decision, and three rules are compared in the same budget of 600 updates.
• Momentum descent brings the training loss down to 0.2777 and the adaptive step to 0.2923; plain descent stays at 0.3296 and cannot reach the 0.40 threshold at all at the smallest rate.
• The test column turns around: the lowest loss gives 0.7500, the highest test number gives 0.7897. The optimizer lowers the loss, it does not lower the generalization error.
• In four runs that change only the initialization, the test spread of a single optimizer is 0.0317, that is, eight subscribers; the difference between optimizers, on the other hand, moves between 0.0516 and 0.0079 depending on the run.
• The triple: baseline 0.7579, network number 0.7857, second run 0.7698. One of the twelve runs passes the decision tree, none passes the linear model.

Next Step

The network of this lesson had two layers and the gradient was chained only once as it came back. When the number of layers grows, the question does not stay clean: the gradient is filtered by a factor at every layer on its way from the output to the input, and the product of those factors either goes down to zero or overflows. The next lesson ties this to a number — as the network deepens, how many times smaller the gradient magnitude seen by the first layer becomes than that of the last layer, at which initialization gain it turns around and explodes, and which of these two extremes gradient clipping and the learning rate schedule correct, and by how much.