Level 6 · Machine Learning Fundamentals: Beyond the Basics

The first five levels built the core of machine learning: data as vectors, models with parameters, a loss, an optimiser, and honest measurement on unseen data. This level uses that core to step into six larger areas, each small enough to write with plain lists and check by eye: PCA, neural networks trained by backpropagation, embeddings and similarity search, anomaly detection, and two kinds of learning by trial and error, multi-armed bandits and Q-learning.

The thread of the path holds throughout. PCA maximises the variance along a direction. A neural network is trained by the gradient descent of Level 4, with backpropagation as an efficient way to get the gradient. Similarity search is Level 1's closest-vector search at scale. Bandits and Q-learners collect their own data by acting, and must balance using what they know against finding out more.

Links to the other paths: the covariance matrix comes from Statistics level 3, z-scores from Statistics level 2, and the confidence bounds of UCB echo the intervals of Statistics level 5; top-k search uses heaps from the Algorithms path, and a bandit is the budgeted black-box optimisation of Algorithms level 6 with noisy answers.

PCA by power iteration

Real datasets often have many features that move together. Twenty students' marks in maths, physics and art make a dataset with three features, but maths and physics are strongly related, so the data is closer to two-dimensional than three. Principal component analysis (PCA) finds the directions along which the data varies most. The first principal component is the unit vector \mathbf{v} for which the projections \mathbf{v}\cdot\mathbf{x} of the centred data have the largest possible variance: again an optimisation problem.

import math, random

rng = random.Random(7)
marks = []
for _ in range(20):
    s = rng.gauss(0, 1)                                   # shared "science" ability
    marks.append([round(60 + 10 * s + rng.gauss(0, 3)),  # maths
                  round(55 + 9 * s + rng.gauss(0, 3)),   # physics
                  round(65 + rng.gauss(0, 8))])           # art, unrelated
print(marks[:4])

n, d = len(marks), 3
means = [sum(col) / n for col in zip(*marks)]
cov = [[sum((r[i] - means[i]) * (r[j] - means[j]) for r in marks) / n for j in range(d)]
       for i in range(d)]
for row in cov:
    print([round(v, 1) for v in row])

[[59, 52, 62], [50, 50, 68], [71, 66, 66], [46, 42, 69]]
[111.7, 77.7, -7.6]
[77.7, 63.1, -9.9]
[-7.6, -9.9, 55.5]

The covariance matrix (Statistics level 3) holds each feature's variance on the diagonal and the covariance of each pair off it. The large positive entry for maths and physics shows them moving together; art barely relates to either. Write C for the covariance matrix. The variance of the data projected onto a unit vector \mathbf{v} is \mathbf{v}^\top C\,\mathbf{v}, and the \mathbf{v} that maximises it is the covariance matrix's top eigenvector: a direction that the matrix only stretches, C\mathbf{v} = \lambda\mathbf{v}, with \lambda (the eigenvalue) equal to the variance along it.

Power iteration finds it with nothing more than repeated matrix-vector products. Multiplying any vector by the matrix stretches its component along the top eigenvector more than the others, so multiplying again and again, and rescaling to length 1 each time, turns the vector towards that direction:

\mathbf{v} \leftarrow \frac{C\mathbf{v}}{\lVert C\mathbf{v} \rVert}

Each entry of C\mathbf{v} is a dot product of one row of C with \mathbf{v}, (C\mathbf{v})_i = \sum_j C_{ij} v_j, which is the comprehension in mat_vec:

def mat_vec(A, v):
    return [sum(a * b for a, b in zip(row, v)) for row in A]

def power_iteration(A, steps=100):
    v = [1.0] * len(A)
    for _ in range(steps):
        w = mat_vec(A, v)
        length = math.sqrt(sum(x * x for x in w))
        v = [x / length for x in w]
    return v, sum(a * b for a, b in zip(v, mat_vec(A, v)))    # direction and variance

v1, var1 = power_iteration(cov)
total = sum(cov[i][i] for i in range(d))
print([round(x, 3) for x in v1], round(var1, 1), f"{var1 / total:.0%} of the variance")

# remove the first component and repeat to find the second (deflation)
cov2 = [[cov[i][j] - var1 * v1[i] * v1[j] for j in range(d)] for i in range(d)]
v2, var2 = power_iteration(cov2)
print([round(x, 3) for x in v2], round(var2, 1), f"{var2 / total:.0%} of the variance")

[0.8, 0.591, -0.104] 170.1 74% of the variance
[0.126, 0.004, 0.992] 54.5 24% of the variance

The first component weights maths and physics together and ignores art: it is the "science ability" hidden in the data, and it carries almost three quarters of the total variance. The second component is essentially art. Projecting each student onto these two directions describes them with two numbers instead of three while keeping 98% of the variance:

def project(row, directions):
    centred = [x - m for x, m in zip(row, means)]
    return [round(sum(a * b for a, b in zip(centred, v)), 1) for v in directions]

for row in marks[:4]:
    print(row, "->", project(row, [v1, v2]))

[59, 52, 62] -> [3.3, -1.2]
[50, 50, 68] -> [-5.7, 3.6]
[71, 66, 66] -> [20.7, 4.4]
[46, 42, 69] -> [-13.8, 4.1]

PCA is used to compress data, to plot high-dimensional data in two dimensions, and to remove redundant features before another model. It only finds linear directions, and because it maximises variance it is sensitive to scale: standardise the features first (Level 3) unless they share units, as marks do here.

Common mistakes: forgetting to centre the data, so the first "component" just points at the mean. Starting power iteration from a vector perpendicular to the top eigenvector (it can never turn towards it; a random start avoids this). Running too few steps when the top two eigenvalues are close, which makes convergence slow.

Practice: The Direction Where the Data Spreads Most, How Many Directions Does the Panel Need?, Putting the Towns on a Two-Number Map

A small neural network with backpropagation

Every classifier so far drew a straight boundary (logistic regression, the perceptron) or cut the space into rectangles (trees). The simplest problem a straight line cannot solve is XOR: the output is 1 when exactly one of two inputs is 1. The points (0, 1) and (1, 0) must be on one side and (0, 0) and (1, 1) on the other, and no single line does that.

A neural network solves it by stacking simple units. Each hidden unit computes a weighted sum of the inputs plus a bias, as logistic regression does, and passes it through a nonlinear activation function (here tanh, which squashes to between −1 and 1). The output unit then does logistic regression on the hidden units' values. For an input \mathbf{x}, hidden unit j with weight vector \mathbf{w}_j and bias b_j, and output weights v_j with output bias c:

h_j = \tanh(\mathbf{w}_j\cdot\mathbf{x} + b_j), \qquad \hat{y} = \sigma\Bigl(\sum_j v_j h_j + c\Bigr)

A layer is all the hidden units together: the same input goes into each, and the list of their outputs h_j is the layer's output. In the code below \mathbf{w}_j is W1[j], b_j is b1[j], v_j is W2[j] and c is b2; each formula is one comprehension in forward. The hidden layer learns new features in which the problem is linearly separable. With 2 inputs, 3 hidden units and 1 output there are 3·2 + 3 + 3 + 1 = 13 parameters, and the loss is the log loss of Level 4.

X = [[0, 0], [0, 1], [1, 0], [1, 1]]
Y = [0, 1, 1, 0]

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

def init(seed, hidden=3):
    r = random.Random(seed)
    W1 = [[r.uniform(-1, 1) for _ in range(2)] for _ in range(hidden)]   # hidden weights
    b1 = [0.0] * hidden
    W2 = [r.uniform(-1, 1) for _ in range(hidden)]                        # output weights
    return {"W1": W1, "b1": b1, "W2": W2, "b2": 0.0}

def forward(p, x):
    h = [math.tanh(sum(w * xi for w, xi in zip(row, x)) + b) for row, b in zip(p["W1"], p["b1"])]
    out = sigmoid(sum(w * hi for w, hi in zip(p["W2"], h)) + p["b2"])
    return h, out

def loss(p):
    total = 0.0
    for x, y in zip(X, Y):
        o = forward(p, x)[1]
        total -= math.log(o) if y == 1 else math.log(1 - o)
    return total / len(X)

p = init(seed=0)
print(round(loss(p), 4), [round(forward(p, x)[1], 3) for x in X])

0.7009 [0.5, 0.612, 0.599, 0.669]

Training is gradient descent, so you need the gradient of the loss with respect to all 13 parameters. Backpropagation computes it with the chain rule, working backwards from the output. For one example:

At the output, the log loss with a sigmoid gives the same simple error as in Level 4: \delta = \hat{y} - y (delta = out - y). The gradient for each output weight v_j is \delta\,h_j, and for the output bias \delta.
Each hidden unit j influenced the output through the weight v_j (W2[j]), so its share of the error is \delta\,v_j, multiplied by the slope of \tanh at that point, which is 1 - h_j^2. Call that \delta_j = \delta\,v_j\,(1 - h_j^2). The gradient for hidden weight W1[j][i] is \delta_j\,x_i, and for the hidden bias \delta_j.

Sum over the examples and divide by their number:

def gradients(p):
    g = {"W1": [[0.0, 0.0] for _ in p["W1"]], "b1": [0.0] * len(p["b1"]),
         "W2": [0.0] * len(p["W2"]), "b2": 0.0}
    for x, y in zip(X, Y):
        h, out = forward(p, x)
        delta = out - y                                   # output error
        g["b2"] += delta / len(X)
        for j, hj in enumerate(h):
            g["W2"][j] += delta * hj / len(X)
            dj = delta * p["W2"][j] * (1 - hj ** 2)       # error passed back to hidden unit j
            g["b1"][j] += dj / len(X)
            for i, xi in enumerate(x):
                g["W1"][j][i] += dj * xi / len(X)
    return g

# check one entry against a numerical gradient (Level 4)
g = gradients(p)
eps = 1e-6
p["W1"][0][1] += eps; up = loss(p)
p["W1"][0][1] -= 2 * eps; down = loss(p)
p["W1"][0][1] += eps
print(round(g["W1"][0][1], 6), round((up - down) / (2 * eps), 6))

-0.013962 -0.013962

The two agree, so the backpropagation formulas are right. (This check is not optional: a wrong gradient still makes a network learn something, so bugs in backpropagation are easy to miss.) The cost matters: backpropagation gives all the partial derivatives for about the price of two forward passes, while the numerical method needs two forward passes per parameter. For a network with millions of parameters, that difference is what makes training possible.

Now train with the plain gradient descent loop of Level 4:

def step(p, g, lr):
    p["W1"] = [[w - lr * gw for w, gw in zip(row, grow)] for row, grow in zip(p["W1"], g["W1"])]
    p["b1"] = [b - lr * gb for b, gb in zip(p["b1"], g["b1"])]
    p["W2"] = [w - lr * gw for w, gw in zip(p["W2"], g["W2"])]
    p["b2"] -= lr * g["b2"]

p = init(seed=0)
for epoch in range(1, 2001):
    step(p, gradients(p), lr=0.5)
    if epoch in (1, 100, 500, 2000):
        print(epoch, round(loss(p), 4), [round(forward(p, x)[1], 3) for x in X])

1 0.695 [0.486, 0.594, 0.582, 0.652]
100 0.4905 [0.164, 0.611, 0.671, 0.59]
500 0.0207 [0.026, 0.986, 0.986, 0.028]
2000 0.0033 [0.004, 0.998, 0.998, 0.004]

The network has learned XOR: near 0 for (0, 0) and (1, 1), near 1 for the other two. Deep learning is this recipe at scale (more layers and units, better activations and optimisers, mini-batches), but still a loss, a gradient by backpropagation and a descent step. The loss is not convex, which is why the starting weights are random (equal weights would give every hidden unit the same gradient forever) and why a different seed can land in a different minimum.

The two activation functions of this network look alike but have different ranges: tanh runs from −1 to 1, the sigmoid from 0 to 1. Both are steepest at 0 and flatten out on either side, which is why the slope 1 - h_j^2 is small for a unit that is already close to −1 or 1. Press Run to compare them:

import math
zs = [z / 10 for z in range(-50, 51)]
plot(zs, [math.tanh(z) for z in zs], label="tanh", xlabel="z", title="Activation functions")
plot(zs, [1 / (1 + math.exp(-z)) for z in zs], label="sigmoid")

Common mistakes: initialising all weights to zero, so the hidden units stay identical. Updating some parameters before computing the gradients of others. Using the output of tanh where its derivative is needed, or the other way round. Skipping the gradient check.

Practice: What Does the Little Network Say?, One Step of Backpropagation, Training the Network From a Fixed Start

Embeddings and similarity search

An embedding represents each item (a word, a film, a product, a user) as a vector, chosen so that similar items get nearby vectors. Once items are embedded, "find things like this one" is the closest-vector search of Level 1, and cosine similarity is the usual measure. The interesting question is where the vectors come from. Hand-picked features work for produce; for films or words nobody can list the right features, so the vectors are learned from data.

A simple source is behaviour. Here twelve viewers have rated eight films from 1 to 5; each row is a film and each column a viewer:

import heapq

ratings = {
    "Rocket Run":     [2, 3, 4, 1, 2, 4, 4, 2, 3, 4, 3, 4],
    "Steel Storm":    [2, 2, 4, 1, 3, 5, 4, 3, 3, 3, 4, 4],
    "Night Chase":    [1, 2, 3, 1, 3, 3, 4, 3, 4, 2, 4, 3],
    "Summer Letters": [3, 3, 2, 4, 2, 3, 4, 1, 5, 4, 1, 3],
    "Paris Again":    [3, 4, 2, 4, 2, 3, 3, 1, 5, 3, 1, 3],
    "Two Hearts":     [3, 4, 1, 4, 3, 3, 2, 2, 5, 4, 1, 4],
    "Space Wedding":  [3, 3, 2, 3, 1, 4, 3, 2, 5, 3, 2, 4],
    "Quiet Garden":   [4, 4, 4, 4, 4, 1, 1, 3, 3, 3, 3, 2],
}
films = list(ratings)
viewers = len(ratings[films[0]])
# some viewers rate everything high: subtract each viewer's average first
avg = [sum(ratings[f][u] for f in films) / len(films) for u in range(viewers)]
emb = {f: [ratings[f][u] - avg[u] for u in range(viewers)] for f in films}

def cosine(u, v):
    return sum(a * b for a, b in zip(u, v)) / (math.hypot(*u) * math.hypot(*v))

def most_similar(query, k=2):
    others = (f for f in films if f != query)
    return [(f, round(cosine(emb[query], emb[f]), 2))
            for f in heapq.nlargest(k, others, key=lambda f: cosine(emb[query], emb[f]))]

for film in ("Rocket Run", "Paris Again", "Quiet Garden"):
    print(film, "->", most_similar(film))

Rocket Run -> [('Steel Storm', 0.78), ('Night Chase', 0.43)]
Paris Again -> [('Summer Letters', 0.81), ('Two Hearts', 0.68)]
Quiet Garden -> [('Two Hearts', 0.09), ('Paris Again', -0.02)]

Nobody told the program which films are action films and which are romances, yet the ratings reveal it: people who liked Rocket Run also liked Steel Storm. This is the idea behind collaborative filtering in recommender systems. Quiet Garden has no close match (its best similarity is near 0), which is itself useful information: it appeals to a different audience from everything else.

These vectors have one entry per viewer, so with a million viewers they would have a million entries. PCA from the first section compresses them: treat each film's vector as a data point, and its projections onto the top two principal components give every film a two-number embedding. For these ratings the first coordinate runs from action to romance and the second singles out Quiet Garden, and the two numbers keep over 80% of the variance. Word embeddings are built in the same spirit: represent each word by the contexts it appears in and compress, or train a small network (as in the previous section) whose hidden-layer values become the embedding.

Similarity search at scale is the Algorithms side of embeddings. Comparing a query with every stored vector costs one dot product per item: fine for thousands, too slow for billions. Normalise every vector once (Level 1) so that cosine similarity is a single dot product, keep only the top k with heapq.nlargest, and beyond that use approximate nearest-neighbour indexes, which group the vectors into clusters (as k-means does in Level 5) and search only the nearest clusters, trading a little accuracy for enormous speed.

Common mistakes: comparing raw ratings without removing each viewer's average, so generous viewers make everything look similar. Treating a missing rating as 0, which reads as "hated it". Forgetting that an embedding learned from few examples is noisy; similarities based on a handful of ratings deserve little trust.

Practice: Films Like This One, Paris Is to France as Rome Is to ..., Search Only the Promising Shelves

Anomaly detection

An anomaly is an example unlike the data a system normally sees: a failing machine, a fraudulent payment, a broken sensor. Anomalies are rare and each one is different, so there are usually too few labelled examples to train a classifier. Instead, learn what normal looks like from ordinary data, and flag whatever does not fit.

Here a machine reports its temperature (°C) and vibration (mm/s) every minute. Both rise with the load it is running at, so in normal operation they rise and fall together:

def normal_readings(n, seed):
    rng = random.Random(seed)
    readings = []
    for _ in range(n):
        load = rng.random()
        readings.append([round(40 + 20 * load + rng.gauss(0, 1.5), 1),
                         round(2 + 6 * load + rng.gauss(0, 0.4), 2)])
    return readings

normal = normal_readings(200, seed=1)
means = [sum(col) / len(col) for col in zip(*normal)]
sds = [math.sqrt(sum((v - m) ** 2 for v in col) / len(col)) for col, m in zip(zip(*normal), means)]
print([round(m, 2) for m in means], [round(s, 2) for s in sds])

[50.29, 5.0] [5.76, 1.69]

The simplest detector uses the z-score of Statistics level 2: how many standard deviations each feature is from its normal mean, z = (x - \mu)/\sigma. A reading with any |z| > 3 is flagged, since that is rare for roughly bell-shaped data.

def z_scores(r):
    return [(v - m) / s for v, m, s in zip(r, means, sds)]

new = [[50.2, 4.9], [58.1, 7.6], [75.0, 8.3], [59.0, 2.6], [43.0, 7.5]]
for r in new:
    z = z_scores(r)
    print(r, [round(v, 2) for v in z], "ANOMALY" if max(abs(v) for v in z) > 3 else "ok")

[50.2, 4.9] [-0.02, -0.06] ok
[58.1, 7.6] [1.36, 1.54] ok
[75.0, 8.3] [4.29, 1.95] ANOMALY
[59.0, 2.6] [1.51, -1.42] ok
[43.0, 7.5] [-1.27, 1.48] ok

It catches the overheating reading [75.0, 8.3], but it passes [59.0, 2.6] (hot but barely vibrating) and [43.0, 7.5] (cool but shaking hard). Each feature on its own is within the normal range; it is the combination that never happens in normal operation. Feature-by-feature checks cannot see that.

A distance-based detector can. Standardise the features (Level 3), then score each reading by its distance to its k-th nearest normal reading: in a crowded region of normal data that distance is small, and in empty space it is large. To pick the threshold, score every normal reading the same way (leaving itself out) and take a high percentile of those scores, so that only about 1% of normal readings would be flagged:

Z_normal = [z_scores(r) for r in normal]

def knn_score(z, k=5, skip=None):
    dists = (math.dist(z, p) for i, p in enumerate(Z_normal) if i != skip)
    return heapq.nsmallest(k, dists)[-1]                  # distance to the k-th nearest

train_scores = sorted(knn_score(z, skip=i) for i, z in enumerate(Z_normal))
threshold = train_scores[int(0.99 * len(train_scores))]
print("threshold:", round(threshold, 3))
for r in new:
    s = knn_score(z_scores(r))
    print(r, round(s, 2), "ANOMALY" if s > threshold else "ok")

threshold: 0.545
[50.2, 4.9] 0.1 ok
[58.1, 7.6] 0.14 ok
[75.0, 8.3] 2.55 ANOMALY
[59.0, 2.6] 1.7 ANOMALY
[43.0, 7.5] 1.66 ANOMALY

Now all three unusual readings stand out, with scores three to five times the threshold, while the two normal ones score far below it. As in Level 3, every detector trades false alarms against missed anomalies (precision against recall), and the percentile sets where you sit. Other detectors follow the same pattern of "model normal, measure misfit": the distance from the main PCA directions, or the distance to the nearest k-means centre.

Common mistakes: fitting the detector on data that already contains anomalies, so they count as normal. Forgetting to standardise, so one feature dominates the distance. Choosing the threshold by eye after seeing the anomalies you want to catch.

Practice: Which Sensor Raised the Alarm?, Readings Far From Any Normal Neighbour

Multi-armed bandits

In every model so far the data arrived first and learning came after. Now the learner collects its own data by acting. A website has three versions of its sign-up page and wants the one that converts best. Each visitor is shown one version and either signs up (reward 1) or not (reward 0). The true sign-up rates are unknown to the learner; the simulation below uses 30%, 50% and 40%. This is a multi-armed bandit, named after a row of slot machines ("one-armed bandits") with different, hidden payout rates.

The tension is between exploitation (show the version that has done best so far) and exploration (show others, in case the current favourite only looks good by luck). The score is the total reward, or equivalently the regret: how much reward was lost compared with always choosing the best arm, which here is 0.5 - \text{rate} per pull of a worse arm.

RATES = [0.30, 0.50, 0.40]

def run_bandit(choose, rounds, seed):
    rng = random.Random(seed)
    pulls, wins = [0] * len(RATES), [0] * len(RATES)
    for t in range(1, rounds + 1):
        arm = choose(pulls, wins, t, rng)
        pulls[arm] += 1
        wins[arm] += rng.random() < RATES[arm]          # True counts as 1
    regret = sum(n * (max(RATES) - r) for n, r in zip(pulls, RATES))
    return sum(wins), pulls, round(regret, 1)

def greedy(pulls, wins, t, rng):
    for arm, n in enumerate(pulls):
        if n == 0:
            return arm                                     # try every arm once
    return max(range(len(pulls)), key=lambda a: wins[a] / pulls[a])

def epsilon_greedy(pulls, wins, t, rng, epsilon=0.1):
    if rng.random() < epsilon:
        return rng.randrange(len(pulls))                  # explore
    return greedy(pulls, wins, t, rng)                     # exploit

for name, strategy in (("greedy", greedy), ("epsilon-greedy", epsilon_greedy)):
    for seed in (0, 2):
        print(name, seed, run_bandit(strategy, 2000, seed))

greedy 0 (575, [1998, 1, 1], 399.7)
greedy 2 (768, [1, 1, 1998], 200.0)
epsilon-greedy 0 (992, [65, 1830, 105], 23.5)
epsilon-greedy 2 (933, [162, 1710, 128], 45.2)

Pure greedy can lock onto the wrong arm for good: after one try of each page it settled on the 30% page with seed 0 and on the 40% page with seed 2, because the best page happened to fail its single early try, and greedy never looked at it again. Epsilon-greedy spends a fixed 10% of rounds on random choices, which guarantees every arm keeps being sampled, and it finds the best arm with either seed.

A fixed exploration rate wastes effort once the answer is clear. The upper confidence bound strategy (UCB1) explores where it is uncertain instead. It picks the arm with the highest optimistic estimate: the observed rate plus a bonus that is large for arms tried few times and shrinks as evidence accumulates. In round t, for an arm a pulled n_a times with w_a wins so far, the score is

\frac{w_a}{n_a} + \sqrt{\frac{2\ln t}{n_a}}

and the code below takes the max of it over the arms. It is the idea of a confidence interval from Statistics level 5, used as "optimism in the face of uncertainty": an arm is dropped only once even its optimistic estimate is poor.

def ucb1(pulls, wins, t, rng):
    for arm, n in enumerate(pulls):
        if n == 0:
            return arm
    return max(range(len(pulls)), key=lambda a: wins[a] / pulls[a] + math.sqrt(2 * math.log(t) / pulls[a]))

for rounds in (1000, 10000, 50000):
    for name, strategy in (("epsilon-greedy", epsilon_greedy), ("UCB1", ucb1)):
        print(rounds, name, run_bandit(strategy, rounds, 0)[1:])

1000 epsilon-greedy ([32, 897, 71], 13.5)
1000 UCB1 ([90, 686, 224], 40.4)
10000 epsilon-greedy ([339, 9287, 374], 105.2)
10000 UCB1 ([231, 9078, 691], 115.3)
50000 epsilon-greedy ([1651, 46616, 1733], 503.5)
50000 UCB1 ([392, 48417, 1191], 197.5)

Early on UCB1 explores more heavily and has the higher regret. But epsilon-greedy's regret grows in proportion to the number of rounds (10% of visitors always see a random page), while UCB1's grows only like the logarithm, because its bonus for the weaker arms keeps shrinking. By 50,000 rounds UCB1 has lost less than half as much. A bandit is an optimisation problem where each evaluation is noisy and costs something, the budgeted black-box setting of Algorithms level 6 with random answers.

Common mistakes: dividing by zero for an arm that has never been pulled (try each arm once first). Running greedy "because the data will sort it out"; it does not. Judging strategies on one run; regret is itself random, so compare over several seeds.

Practice: What Did the Exploring Cost?, Mostly Greedy, Sometimes Curious, Optimism About the Untested Page

Q-learning on a tiny grid

A bandit's choices have no lasting effect: the next visitor is the same whatever you showed the last one. In reinforcement learning actions change the situation, so a good action now may pay off only several steps later. The standard setting has states, actions, and after each action a reward and a new state. The learner's goal is a policy, a rule for choosing an action in each state, that maximises the total future reward.

Here is a robot on a 3 by 4 grid. It starts at S, the goal G is worth +10, the pit P costs 10, # is a wall, and every move costs 1, so shorter routes are better. Moves are the four directions; walking into the wall or the edge leaves the robot where it is.

GRID = ["...G",
        ".#.P",
        "S..."]
MOVES = {"^": (-1, 0), "v": (1, 0), "<": (0, -1), ">": (0, 1)}

def act(state, move):
    r, c = state
    dr, dc = MOVES[move]
    nr, nc = r + dr, c + dc
    if not (0 <= nr < 3 and 0 <= nc < 4) or GRID[nr][nc] == "#":
        nr, nc = r, c                                  # blocked: stay put
    cell = GRID[nr][nc]
    if cell == "G":
        return (nr, nc), 10, True                      # new state, reward, episode over
    if cell == "P":
        return (nr, nc), -10, True
    return (nr, nc), -1, False

print(act((2, 0), ">"), act((1, 2), ">"), act((2, 1), "^"))

((2, 1), -1, False) ((1, 3), -10, True) ((2, 1), -1, False)

The third call walks into the wall and stays put.

Q-learning learns a table Q[state, move]: an estimate of the total future reward from making that move in that state and playing well afterwards. The true values satisfy a consistency condition (the Bellman equation): the value of a move is its immediate reward plus the discounted value of the best move from wherever it leads, r + \gamma \max_{a'} Q(s', a'), where s' is the next state and the maximum runs over the moves a' available there. The discount \gamma (gamma, here 0.9) makes rewards sooner worth more than rewards later. Q-learning nudges each estimate towards that target every time it tries the move a in state s:

Q(s, a) \leftarrow Q(s, a) + \alpha\bigl(r + \gamma \max_{a'} Q(s', a') - Q(s, a)\bigr)

with learning rate \alpha. The update is gradient descent again, on the squared gap between the estimate and its target. The robot chooses moves epsilon-greedily, exactly as in the bandit section, so it keeps exploring while it learns:

Q = {}
def q(s, m):
    return Q.get((s, m), 0.0)

rng = random.Random(0)
alpha, gamma, epsilon = 0.5, 0.9, 0.2
steps_taken = []
for episode in range(300):
    s = (2, 0)
    for step in range(1, 51):
        if rng.random() < epsilon:
            m = rng.choice(list(MOVES))
        else:
            m = max(MOVES, key=lambda mv: q(s, mv))
        s2, reward, done = act(s, m)
        target = reward if done else reward + gamma * max(q(s2, mv) for mv in MOVES)
        Q[(s, m)] = q(s, m) + alpha * (target - q(s, m))
        s = s2
        if done:
            break
    steps_taken.append(step)
print("first episodes:", steps_taken[:8])
print("last episodes: ", steps_taken[-8:])

for r in range(3):
    print(" ".join(GRID[r][c] if GRID[r][c] in "G#P" else max(MOVES, key=lambda mv: q((r, c), mv))
                   for c in range(4)),
          "  ", " ".join(f"{max(q((r, c), mv) for mv in MOVES):5.2f}" for c in range(4)))

first episodes: [34, 12, 17, 19, 10, 15, 12, 9]
last episodes:  [4, 4, 5, 5, 5, 7, 5, 7]
> > > G     2.95  8.00 10.00  0.00
v # ^ P     1.81  0.00  8.00  0.00
> > ^ <     3.12  4.58  6.20  4.58

The first episodes wander for dozens of steps; the last ones mostly take five, the shortest route. The others are the price of exploration: 20% of moves are still random, which sometimes adds detours and sometimes (the 4-step episodes) tips the robot into the pit. The printed policy shows the learned route: right along the bottom, up the third column, and right into the goal. It has learned to avoid the pit: from the cell below the pit it goes left rather than up, and from the cell beside the pit it goes up rather than right. The values along the route are exactly what the Bellman equation predicts: 10 next to the goal, then 8, 6.2, 4.58 and 3.12 at the start, each one −1 + 0.9 × the next. Cells the robot rarely visits, such as the top-left corner, have values that have not yet converged, because Q-learning only learns about what it tries.

This is the path's optimisation idea again, with the data gathered by acting: the table is the model, the Bellman gap is the loss, and exploration decides which data exists. Replacing the table by a neural network that maps a state to its Q-values gives deep Q-learning, which learned to play video games from screen pixels.

Common mistakes: using the value of the move the robot takes next instead of the best move in the target (that is a different algorithm, SARSA). Forgetting that a finished episode has no future (the target is just the reward). Setting epsilon to 0 from the start, so the robot never discovers a better route than its first one. Not capping the episode length, so an early random walk runs forever.

Practice: Learning From the Robot's Diary, The Robot Learns Its Way Around the Pit

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.

Checklist

Build a covariance matrix, find its top principal components by power iteration and deflation, and project data onto them.
Explain why a single linear boundary cannot solve XOR, and train a 2-3-1 network with backpropagation and gradient descent.
Derive the backpropagation formulas for one hidden layer and verify them with a numerical gradient.
Build item embeddings from behaviour, find similar items by cosine similarity with a heap, and compress embeddings with PCA.
Detect anomalies with z-scores and with k-th nearest-neighbour distances, choose a threshold from normal data, and explain why combinations need the second method.
Simulate a multi-armed bandit and compare greedy, epsilon-greedy and UCB1 by regret.
Implement tabular Q-learning with epsilon-greedy exploration, read off the learned policy, and check its values against the Bellman equation.
Describe PCA, network training, bandits and Q-learning each as an optimisation problem over data.