Level 4 · Machine Learning Fundamentals: Learning as Optimisation

The models of Level 3 were each trained in their own way: kNN stores the data, the least-squares line has a formula, Naive Bayes counts. This level replaces those special cases with one idea that covers almost all of modern machine learning. A model has parameters (numbers you can adjust), a loss function measures how badly the model with those parameters fits the training data, and learning is the search for the parameters with the smallest loss. Machine learning is an optimisation problem over data, and this level, together with the two that follow, is built around that sentence.

The search method is gradient descent: work out which direction makes the loss fall fastest, take a small step that way, and repeat. You will compute gradients numerically and exactly, see how the step size decides between crawling, converging and exploding, train logistic regression (a classifier that outputs probabilities) with the same loop, meet the older mistake-driven perceptron, and ask when a local search like this is guaranteed to find the best answer.

If you have done Algorithms level 6, this will feel familiar: there you optimised routes and black-box functions by hill climbing and annealing, and fought local optima. Gradient descent is downhill hill climbing that uses the slope to choose each step, which is possible because here you can see inside the loss. The least-squares line of Statistics level 3 returns as a test case: gradient descent must find the same line the formula gives.

Parameters and the loss

Take the six deliveries from Level 2 and the model minutes ≈ w · km + b, or \hat{y} = wx + b. The model is the formula; the parameters are w (minutes per kilometre) and b (fixed minutes). Every choice of (w, b) is a different line, and the mean squared error turns each choice into one number:

The sum over the n deliveries is the generator expression inside loss:

km      = [1, 2, 3, 4, 5, 6]
minutes = [9, 11, 15, 16, 21, 24]

def loss(w, b):
    return sum((w * x + b - y) ** 2 for x, y in zip(km, minutes)) / len(km)

print(round(loss(0, 0), 2))       # predicts 0 for everything
print(round(loss(3, 6), 2))       # the guessed rule from Level 2
print(round(loss(3.03, 5.4), 2))  # close to the least-squares line

283.33
0.83
0.58

The loss is a function of the parameters, with the data held fixed. That is the change of viewpoint this level is about. Earlier you asked "given this model, how wrong is it?"; now you ask "which parameters make it least wrong?". Picture the loss as a landscape over the (w, b) plane, with height equal to the loss; training means finding the lowest point. Here is a slice through that landscape with b fixed at 5:

for w in [2.0, 2.5, 3.0, 3.5, 4.0]:
    print(w, round(loss(w, 5), 2), "#" * round(loss(w, 5)))

2.0 19.67 ####################
2.5 6.46 ######
3.0 0.83 #
3.5 2.79 ###
4.0 12.33 ############

A bowl, with its bottom near w = 3. With two parameters you could scan a fine grid and take the best point. With the thousands or millions of parameters of real models a grid is hopeless (a grid of just 10 values per parameter has 10¹⁰⁰⁰ points for 1,000 parameters), so you need a method that looks only at the neighbourhood of where it is and still heads for the bottom.

Common mistakes: confusing the model's input (km) with its parameters (w, b); in training the data is fixed and the parameters move. Scoring parameters on the test set; the loss you minimise is always the training loss.

Numerical gradients

The derivative of a function of one variable is its slope: how much the output changes per unit of input, for a tiny change. For a function of several parameters, the partial derivative with respect to one parameter is the slope in that parameter's direction with the others held still, and the gradient is the vector of all the partial derivatives. The gradient points uphill, in the direction of fastest increase, so its negative points downhill.

You can estimate a slope without any calculus: nudge the parameter a tiny amount h in each direction and see how much the loss changes. This central difference is accurate enough for most purposes:

The symbol \partial f / \partial p is the partial derivative of f with respect to the parameter p; the code applies the formula to each parameter in turn:

def numerical_gradient(f, params, h=1e-5):
    grad = []
    for i in range(len(params)):
        up = params[:i] + [params[i] + h] + params[i + 1:]
        down = params[:i] + [params[i] - h] + params[i + 1:]
        grad.append((f(up) - f(down)) / (2 * h))
    return grad

f = lambda p: loss(p[0], p[1])
print([round(g, 4) for g in numerical_gradient(f, [3, 5])])
print([round(g, 4) for g in numerical_gradient(f, [2, 5])])

[-3.6667, -1.0]
[-34.0, -8.0]

At (3, 5) both slopes are negative: increasing w or b a little would lower the loss, so both should go up slightly. At (2, 5) the w-slope is strongly negative: the line is too flat, and increasing w lowers the loss fast.

For the mean squared error you can also differentiate by hand. Each example contributes (w x_i + b - y_i)^2; by the chain rule its slope in w is 2(w x_i + b - y_i)\,x_i and in b is 2(w x_i + b - y_i). Average over the examples:

The code names the bracket errors once and reuses it in both sums:

def gradient(w, b):
    n = len(km)
    errors = [w * x + b - y for x, y in zip(km, minutes)]
    dw = 2 * sum(e * x for e, x in zip(errors, km)) / n
    db = 2 * sum(errors) / n
    return [dw, db]

print([round(g, 4) for g in gradient(3, 5)])
print([round(g, 4) for g in gradient(2, 5)])

[-3.6667, -1.0]
[-34.0, -8.0]

The two methods agree to four decimals. In practice you use the exact (analytic) gradient, because the numerical one needs two loss evaluations per parameter, which is far too slow for many parameters. But the numerical gradient is the perfect check: whenever you derive a gradient by hand, compare it with the numerical one at a few points before trusting it. You will do exactly that for a neural network in Level 6.

Common mistakes: a step h that is too large (the estimate is inaccurate) or too small (like 1e-12, where floating-point rounding swamps the difference). Using the one-sided difference (f(x + h) - f(x)) / h, which is noticeably less accurate than the central one. Modifying the parameter list in place and forgetting to restore it.

Practice: Slopes of a Black Box, Which Way Is Downhill for the Rent Model?, Catch the Wrong Derivative

Gradient descent

Gradient descent turns the gradient into an algorithm. Start with any parameters, then repeat: compute the gradient, and move each parameter a small step against it,

where the learning rate \eta (eta) sets the step size. In the code \eta is lr, and the arrow is the tuple assignment w, b = w - lr * dw, b - lr * db. Each step goes downhill, and near the bottom the gradient shrinks, so the steps shrink too and the parameters settle.

def gradient_descent(lr, steps, report_every):
    w, b = 0.0, 0.0
    for step in range(1, steps + 1):
        dw, db = gradient(w, b)
        w, b = w - lr * dw, b - lr * db
        if step % report_every == 0:
            print(step, round(w, 4), round(b, 4), round(loss(w, b), 4))
    return w, b

w, b = gradient_descent(lr=0.05, steps=1000, report_every=200)

200 3.0547 5.2881 0.5834
400 3.0292 5.3972 0.581
600 3.0286 5.3999 0.581
800 3.0286 5.4 0.581
1000 3.0286 5.4 0.581

Compare with the least-squares formula from Level 3 (and Statistics level 3), which gives the exact minimum directly:

mx, my = sum(km) / len(km), sum(minutes) / len(minutes)
slope = sum((x - mx) * (y - my) for x, y in zip(km, minutes)) / sum((x - mx) ** 2 for x in km)
print(round(slope, 4), round(my - slope * mx, 4))

3.0286 5.4

Gradient descent found the same line, to four decimals, without knowing the formula. That is the point. The formula works only for this particular model and loss; gradient descent works for any model whose loss you can differentiate, including the logistic regression below and the neural network of Level 6, which have no formula at all.

Each step above uses the whole training set (batch gradient descent). With millions of examples that is slow, so large models use stochastic gradient descent: each step uses one example or a small random mini-batch, a noisy but cheap estimate of the full gradient. The perceptron below learns one example at a time in this way.

Common mistakes: updating w and then using the new w to compute the gradient for b in the same step; compute the whole gradient first, then update every parameter (the tuple assignment above does that). Adding the gradient instead of subtracting it, which climbs the loss. Stopping after a fixed number of steps without checking that the loss has stopped falling.

Practice: Descend Until the Loss Stops Falling, How Long Until Descent Finds the Formula?, Mini-Batches in a Reproducible Order

The learning rate

The learning rate is the most important setting in gradient descent, and the easiest to get wrong. Run the same loop with three values and look at the loss after 100 steps:

def loss_after(lr, steps):
    w, b = 0.0, 0.0
    for _ in range(steps):
        dw, db = gradient(w, b)
        w, b = w - lr * dw, b - lr * db
    return loss(w, b)

for lr in (0.005, 0.05, 0.06, 0.065):
    print(lr, f"{loss_after(lr, 100):.4g}")

0.005 3.225
0.05 0.6768
0.06 0.6265
0.065 9.223e+08

The loss after every step, drawn as a loss curve, shows the difference at a glance. Press Run to see the loss fall slowly with the small rate and fast with the larger one (the vertical axis is logarithmic):

km, minutes = [1, 2, 3, 4, 5, 6], [9, 11, 15, 16, 21, 24]
for lr in (0.005, 0.05):
    w, b, losses = 0.0, 0.0, []
    for step in range(200):
        errors = [w * x + b - y for x, y in zip(km, minutes)]
        losses.append(sum(e * e for e in errors) / 6)
        w, b = w - lr * 2 * sum(e * x for e, x in zip(errors, km)) / 6, b - lr * 2 * sum(errors) / 6
    plot(losses, label=f"learning rate {lr}", xlabel="step", ylabel="training loss", logy=True)

Too small (0.005) and progress is slow: after 100 steps the loss is still far above its minimum of 0.581. In the useful range (0.05, 0.06) it is close. At 0.065 the loss is nearly a billion: each step overshoots the bottom of the bowl and lands higher up the other side than it started, and the overshoots grow until, given enough steps, the numbers overflow. The boundary between converging and exploding is sharp, and it depends on how steeply the bowl curves: the steeper the curve, the smaller the step must be.

Why does even a good rate need hundreds of steps? This bowl is a long, narrow valley: the loss changes fast in some directions (changing w moves every prediction, the long deliveries most) and slowly in others (pivoting the line about the middle of the data). The steep direction limits the step size, and the gentle direction then needs many small steps.

Feature scaling fixes this, for the same reason it fixed kNN in Level 3. Standardise the kilometres to mean 0 and standard deviation 1 and the valley becomes a round bowl:

mean_km = sum(km) / len(km)
sd_km = (sum((x - mean_km) ** 2 for x in km) / len(km)) ** 0.5
z = [(x - mean_km) / sd_km for x in km]

w, b = 0.0, 0.0
for step in range(1, 51):
    errors = [w * x + b - y for x, y in zip(z, minutes)]
    dw = 2 * sum(e * x for e, x in zip(errors, z)) / len(z)
    db = 2 * sum(errors) / len(z)
    w, b = w - 0.1 * dw, b - 0.1 * db
    if step in (10, 20, 50):
        print(step, round(w, 4), round(b, 4), round(sum(e * e for e in errors) / len(z), 4))
print("back in km:", round(w / sd_km, 4), round(b - w * mean_km / sd_km, 4))

10 4.6169 14.282 5.6746
20 5.1126 15.8155 0.6397
50 5.1722 15.9998 0.581
back in km: 3.0285 5.3999

In 50 steps with standardised input the loss reaches the minimum, which took about 400 steps before, and converting the parameters back to kilometres gives the same line again. Scaling features is standard practice before gradient descent for exactly this reason.

There is no universal best learning rate. A practical recipe: try values spaced by factors of about 3 (0.001, 0.003, 0.01, 0.03, 0.1, …), watch the loss, and pick the largest value for which the loss falls steadily. If the loss ever goes up, the rate is too big.

Common mistakes: choosing a tiny learning rate "to be safe" and then stopping long before convergence. Blaming the model when the real problem is unscaled features. Not printing the loss during training; a loss curve is the first thing to look at when training misbehaves.

Practice: The Largest Learning Rate That Behaves, Standardise, Descend, Translate Back, Where Gradient Descent Tips Over

Logistic regression

Gradient descent can train classifiers too, once you choose a model that outputs something smooth. A line w·x + b outputs any number; for classification you want a probability between 0 and 1. The sigmoid (logistic) function squashes any number into that range:

It is 0.5 at z = 0, close to 1 for large positive z, close to 0 for large negative z. Logistic regression predicts p = \sigma(wx + b): the probability that the label is 1. Here is a small dataset of hours studied and whether each student passed (1) or failed (0):

import math

def sigmoid(z):
    return 1 / (1 + math.exp(-z))

hours  = [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5]
passed = [0, 0, 0, 1, 0, 1, 0, 1, 1, 1]
print([round(sigmoid(z), 3) for z in (-4, -1, 0, 1, 4)])

[0.018, 0.269, 0.5, 0.731, 0.982]

The loss must reward confident correct probabilities and punish confident wrong ones. The standard choice is the log loss (also called cross-entropy): for an example with label 1 it is -\log p, for label 0 it is -\log(1 - p). Averaged over n examples with labels y_i and predicted probabilities p_i, both cases fit in one formula, because one of the two terms is always zero:

Predicting 0.99 for a student who passed costs 0.01; predicting 0.01 costs 4.6, and predicting exactly 0 would cost infinity. Minimising the average log loss is the same as maximising the probability the model assigns to the observed labels, the principle of maximum likelihood from Statistics level 6.

Its gradient turns out to be remarkably simple: for each example the error is p_i - y_i, and the partial derivatives are

the same shape as for the squared error of a line. The training loop is the familiar one:

def log_loss(w, b):
    total = 0.0
    for x, y in zip(hours, passed):
        p = sigmoid(w * x + b)
        total -= math.log(p) if y == 1 else math.log(1 - p)
    return total / len(hours)

w, b = 0.0, 0.0
for step in range(1, 1001):
    errors = [sigmoid(w * x + b) - y for x, y in zip(hours, passed)]
    dw = sum(e * x for e, x in zip(errors, hours)) / len(hours)
    db = sum(errors) / len(hours)
    w, b = w - 0.5 * dw, b - 0.5 * db
    if step in (1, 10, 100, 1000):
        print(step, round(w, 3), round(b, 3), round(log_loss(w, b), 4))

1 0.237 0.0 0.6458
10 0.349 -0.526 0.5802
100 1.021 -2.68 0.445
1000 1.353 -3.721 0.4335

Once trained, the model gives a probability for any number of hours, and a decision if you threshold it at 0.5. The threshold is crossed where w·x + b = 0, the decision boundary:

for h in (1, 2, 2.75, 3.5, 5):
    print(h, round(sigmoid(w * h + b), 3))
print("boundary at", round(-b / w, 2), "hours")

1 0.086
2 0.266
2.75 0.5
3.5 0.734
5 0.955
boundary at 2.75 hours

Press Run to see the students as dots at 0 and 1 and the fitted S-shaped curve passing 0.5 at the boundary (the weights are the ones just learned):

import math
hours  = [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5]
passed = [0, 0, 0, 1, 0, 1, 0, 1, 1, 1]
grid = [h / 10 for h in range(61)]
plot(hours, passed, kind="scatter", label="students", xlabel="hours studied", ylabel="P(pass)")
plot(grid, [1 / (1 + math.exp(-(1.353 * h - 3.721))) for h in grid], label="sigmoid(w·x + b)")

The data is noisy (a student with 2 hours passed and one with 3.5 failed), and the model reflects it: rather than drawing a hard line, it says a student with 3.5 hours passes about 73% of the time. With several features the model is p = \sigma(\mathbf{w}\cdot\mathbf{x} + b), with a weight vector and the dot product of Level 1, and its decision boundary is a straight line (or flat plane) through feature space. Each weight says how much one unit of its feature shifts the log-odds of the label.

Common mistakes: training logistic regression with squared error instead of log loss; it works poorly because the loss becomes flat for confident wrong predictions. Computing math.log(0) when a probability rounds to exactly 0 or 1 (clip p to a range like [1e-12, 1 − 1e-12] in real code). math.exp overflowing for large negative z; for big inputs, use 1 − sigmoid(−z) or clip z.

Practice: How Surprised Is the Spam Filter?, One Step Towards Better Click Predictions, When Does the Heating Come On?

The perceptron

The perceptron (1958) is the oldest learning algorithm for a linear boundary, and it learns in a different style: one example at a time, and only from mistakes. It predicts +1 if w·x + b > 0 and −1 otherwise. Whenever it gets an example wrong, it nudges the boundary towards getting that example right: \mathbf{w} \leftarrow \mathbf{w} + y\,\mathbf{x} and b \leftarrow b + y, where y is the true label +1 or -1.

X = [[2, 8], [5, 1], [1, 3], [8, 6], [6, 1], [5, 8], [4, 7], [9, 9]]   # sweetness, crunchiness
y = [-1, 1, -1, 1, 1, -1, -1, 1]                                          # fruit = +1, veg = -1

def dot(u, v):
    return sum(a * b for a, b in zip(u, v))

w, b = [0, 0], 0
for epoch in range(1, 21):
    mistakes = 0
    for x, label in zip(X, y):
        if label * (dot(w, x) + b) <= 0:          # wrong, or exactly on the boundary
            w = [wi + label * xi for wi, xi in zip(w, x)]
            b += label
            mistakes += 1
    print("epoch", epoch, "mistakes", mistakes, "w", w, "b", b)
    if mistakes == 0:
        break

epoch 1 mistakes 5 w [15, 0] b 1
epoch 2 mistakes 3 w [17, -7] b 0
epoch 3 mistakes 2 w [21, -6] b 0
epoch 4 mistakes 3 w [24, -8] b -1
epoch 5 mistakes 1 w [19, -16] b -2
epoch 6 mistakes 0 w [19, -16] b -2

After six epochs the line 19·sweetness − 16·crunchiness − 2 = 0 separates all eight items. The number of mistakes does not fall steadily (epoch 4 makes more than epoch 3), because fixing one example can break another.

A pass over the data is an epoch. Once an epoch makes no mistakes the perceptron has found a line that separates the two classes, and it stops. The perceptron convergence theorem guarantees this happens after finitely many updates whenever some line separates the classes perfectly (the data is linearly separable). If no such line exists, as with the overlapping produce of Levels 2 and 3, the perceptron never settles: it keeps fixing one mistake by making another.

The perceptron fits the optimisation view too: its update is stochastic gradient descent, with learning rate 1, on the loss max(0, −y·(w·x + b)), which is zero for correct examples and grows with how far a mistake is on the wrong side. The choice of loss shows. The perceptron stops at the first separating line, which may pass very close to some examples; logistic regression keeps pushing every example away from the boundary, gives probabilities, and copes with overlapping classes.

Common mistakes: using labels 0 and 1 instead of −1 and +1, so the update w + y·x does nothing for the 0 class. Forgetting the bias b, which forces the boundary through the origin. Looping forever on data that is not separable; always cap the number of epochs.

Practice: The Sorting Robot Learns from Its Mistakes, Keep the Best Line in Your Pocket

Convexity and local minima

Gradient descent only ever looks at the slope where it stands. It stops when the slope is zero, at a local minimum: a point lower than everything around it. Whether that is also the global minimum, the lowest point anywhere, depends on the shape of the loss.

A function is convex if the straight line between any two points on its graph never passes below the graph: a bowl with no dents. A convex function has no local minima except the global one, so gradient descent (with a sensible learning rate) finds the best answer from any starting point. The squared error of a linear model is convex, and so is the log loss of logistic regression. Start the line-fitting loop from three very different places:

for start in ([0.0, 0.0], [10.0, -20.0], [-5.0, 30.0]):
    w, b = start
    for _ in range(2000):
        dw, db = gradient(w, b)
        w, b = w - 0.05 * dw, b - 0.05 * db
    print(start, "->", round(w, 4), round(b, 4))

[0.0, 0.0] -> 3.0286 5.4
[10.0, -20.0] -> 3.0286 5.4
[-5.0, 30.0] -> 3.0286 5.4

All three roads lead to the least-squares line. For convex problems, the only questions are how fast you get there and how precisely.

Most interesting models are not convex. Here is a tiny one: the data comes from a wave, y = sin(2x), and the model sin(w·x) has a single parameter, the frequency w. Its loss landscape has several valleys:

xs = [0.5 * i for i in range(1, 11)]
ys = [round(math.sin(2 * x), 3) for x in xs]

def wave_loss(w):
    return sum((math.sin(w * x) - y) ** 2 for x, y in zip(xs, ys)) / len(xs)

def wave_gradient(w):
    return 2 * sum((math.sin(w * x) - y) * math.cos(w * x) * x for x, y in zip(xs, ys)) / len(xs)

for w in [0.5 * i for i in range(9)]:
    print(f"{w:4}", f"{wave_loss(w):.3f}", "|" + "#" * round(20 * wave_loss(w)))

 0.0 0.500 |##########
 0.5 1.021 |####################
 1.0 1.241 |#########################
 1.5 0.799 |################
 2.0 0.000 |
 2.5 0.748 |###############
 3.0 1.244 |#########################
 3.5 0.881 |##################
 4.0 1.146 |#######################

The true answer w = 2 sits at the bottom of a narrow valley, with hills on both sides and other, shallower valleys beyond them. Gradient descent rolls into whichever valley it starts in:

def descend(w, lr=0.01, steps=500):
    for _ in range(steps):
        w -= lr * wave_gradient(w)
    return w

for start in (1.0, 1.5, 2.5, 3.5):
    w = descend(start)
    print(start, "->", round(w, 3), "loss", round(wave_loss(w), 4))

1.0 -> 0.663 loss 0.9545
1.5 -> 2.0 loss 0.0
2.5 -> 2.0 loss 0.0
3.5 -> 3.591 loss 0.8669

Starting at 1.5 or 2.5 finds the true frequency with zero loss. Starting at 1.0 or 3.5 gets stuck in a local minimum with a loss near 1, and nothing in the gradient says that a better answer exists elsewhere.

This is the problem you met in Algorithms level 6, and the same remedies apply. The simplest is random restarts: run gradient descent from several random starting points and keep the result with the lowest loss.

import random

rng = random.Random(0)
results = [descend(rng.uniform(0, 4)) for _ in range(8)]
best = min(results, key=wave_loss)
print(sorted(round(w, 2) for w in results))
print("best:", round(best, 3), round(wave_loss(best), 4))

[0.66, 2.0, 2.0, 2.0, 2.0, 3.59, 3.59, 3.59]
best: 2.0 0.0

Other remedies carry over too: noise in the updates (stochastic gradient descent jitters out of small dents, much as annealing does), momentum (steps keep part of their previous direction and roll through shallow bumps), and good starting points. The neural networks of Level 6 have highly non-convex losses, yet their local minima tend to be nearly as good as each other, which is one reason gradient descent works for them at all.

Common mistakes: concluding that a model is wrong because one training run got stuck; try other starting points first. Assuming that a zero gradient means the best answer; for a non-convex loss it only means a flat spot. Starting a symmetric model (such as a network whose weights are all equal) at a point where the gradient cannot tell its parameters apart.

Practice: Which Peak Is the Right One?

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

• Slopes of a Black Box
• Which Way Is Downhill for the Rent Model?
• Catch the Wrong Derivative
• How Surprised Is the Spam Filter?
• One Step Towards Better Click Predictions
• How Long Until Descent Finds the Formula?

Checklist

• Describe a model as a formula with parameters, and training as minimising a loss over the training data with the data held fixed.
• Compute a numerical gradient with central differences, derive the analytic gradient of the squared error, and check one against the other.
• Implement batch gradient descent and recover the least-squares line without its formula.
• Diagnose a learning rate as too small, good or too large from the loss curve, and explain why standardised features speed up gradient descent.
• Train logistic regression with the sigmoid, the log loss and its (p − y)·x gradient, and find its decision boundary.
• Implement the perceptron's mistake-driven update, and say when it is guaranteed to stop.
• Explain convexity, show that a convex loss has one minimum, and use random restarts against local minima in a non-convex one.