Chapter 32

Building and Running Our Reader

At the beginning of the laboratory we drew a digit and watched a machine answer. We now know what its stored numbers mean, how they divide image space, and how calculus can help us choose them. Let's actually choose them.

We will start with every parameter equal to zero and write the computations ourselves. Python will keep track of the lists and loops. NumPy will speed up the repeated arithmetic across thousands of images. Neither will decide what a score, a loss, or a derivative ought to be; those are the formulas we have already worked out.

32.1Follow one update all the way through

Before giving the computer thousands of images, give it one very small one. Take an invented two-pixel input $⃗ 𝑥 = (0.5, 1)$ , labeled three. We still use ten digits, so the miniature reader has ten parameter vectors. Each vector has three entries: a bias and two weights.

Append the bias coordinate to get $̃ ⃗ 𝑥 = (1, 0.5, 1)$ . At the zero reader all ten scores are zero and all ten probabilities are $0.1$ . The three's residual is $- 0.9$ ; the other nine residuals are $0.1$ . Therefore the image contributes

\nabla 𝜃 3 ℓ = (- 0.9, - 0.45, - 0.9), \nabla 𝜃 𝑘 ℓ = (0.1, 0.05, 0.1) (𝑘 \neq 3) .

The penalty contributes zero at the zero reader. With step size $𝜂 = 0.2$ , subtracting the gradient gives

𝜃 n e w 3 = (0.18, 0.09, 0.18), 𝜃 n e w 𝑘 = (- 0.02, - 0.01, - 0.02) (𝑘 \neq 3) .

Score the same image again. The three's score is $0.405$ and each competing score is $- 0.045$ . Its probability has risen from $0.1$ to

𝑒 0.405 𝑒 0.405 + 9 𝑒 - 0.045 \approx 0.1484 .

One small step has made the reader less wrong about this image. It has not taught it all handwriting. With the full training set, we will add all the image gradients before changing any parameters.

Figure 32.1 Walk through the zero reader, its probabilities, the residuals, the gradient, the parameter update, and the new probabilities. The ten classes stay present throughout. Only the image dimension has been reduced to make every number inspectable.

32.2Ten lists and a dot product

In the program, theta[k][j] is parameter $𝑗$ for digit $𝑘$ . Index zero holds the bias. The real reader uses a ten-by-65 array, while an augmented image is a list of 65 entries.

Here is the scoring calculation written with ordinary Python loops:

def scores_one(theta, image):
    return [sum(weight * pixel for weight, pixel in zip(row, image))
            for row in theta]

zip pairs each weight with its image coordinate. The inner expression multiplies these pairs and adds them. The outer expression repeats the dot product for the ten parameter rows. This is the same calculation we uncovered in the first chapter.

For the probabilities, subtract the largest score before exponentiating:

def probabilities_one(scores):
    largest = max(scores)
    exponentials = [np.exp(score - largest) for score in scores]
    total = sum(exponentials)
    return [value / total for value in exponentials]

Why change the formula? We have not changed its value. Subtracting the largest score multiplies the numerator and denominator of every probability by the same factor. But now every exponential is at most one, which prevents overflow when scores are large. At least one is exactly one, so their sum cannot be zero.

The derivative for one labeled image is just as short:

def one_image_gradient(theta, image, label):
    probabilities = probabilities_one(scores_one(theta, image))
    return [[(probabilities[k] - int(k == label)) * pixel
             for pixel in image]
            for k in range(10)]

The comparison k == label is true for the known digit; converting it to an integer gives the indicator $𝟏 𝑘 = 𝑦$ . Each residual then multiplies the whole augmented image. The output has the same shape as theta, as a gradient should: one derivative for each parameter.

32.3Do the same arithmetic for all the images

Repeating these Python loops for every image at every step would spend much of our time managing loops. NumPy can perform the same sums and products on whole arrays. We will use two inputs:

images has shape $(𝑁, 65)$ : one augmented image per row.
theta has shape $(10, 65)$ : one parameter vector per digit.

The expression

scores = images @ theta.T

produces an $(𝑁, 10)$ array. Its entry in row $𝑖$ , column $𝑘$ is the dot product of image $𝑖$ with parameter vector $𝑘$ . The symbol @ asks for matrix multiplication and .T transposes the array, arranging those dot products. If matrix notation is new, the entry-by-entry rule is enough: multiply corresponding coordinates and add.

We compute the probabilities separately in each image row:

largest = scores.max(axis=1, keepdims=True)
exponentials = np.exp(scores - largest)
totals = exponentials.sum(axis=1, keepdims=True)
probabilities = exponentials / totals

Here axis=1 means “across the ten digit columns.” Keeping the result as a column lets NumPy subtract or divide by that row's number in all ten places. Images share the parameter array; they do not share their softmax denominator.

To evaluate the objective, use the stable logarithmic expression directly:

losses = (largest[:, 0] + np.log(totals[:, 0])
          - scores[np.arange(len(labels)), labels])
loss = losses.sum() + lam * np.sum(theta * theta) / 2

The indexing in the last term selects the known digit's score from each image row. Computing -log(probability) would be algebraically equivalent, but a computer can round an extremely small probability to zero. The log-sum expression avoids taking the logarithm of that rounded zero.

Now make the residual array, subtracting one in each known-label position:

residuals = probabilities.copy()
residuals[np.arange(len(labels)), labels] -= 1
grad = residuals.T @ images + lam * theta

Check a single entry of that last line. Entry $(𝑘, 𝑗)$ of the product is

𝑁 \sum 𝑖 = 1 (𝑝 (𝑖) 𝑘 - 𝟏 𝑘 = 𝑦 𝑖) ̃ 𝑥 (𝑖) 𝑗 .

Adding lam * theta supplies the penalty derivative. The short array expression is exactly the long gradient formula from the previous chapter.

32.4Check a derivative before trusting it

A program can run without computing the function we intended. Before training, the supplied file checks its derivative on a small problem with ten classes and two image coordinates.

First it sums the explicit one-image gradients and compares them with the NumPy result. Then it varies each parameter by a small amount $ℎ$ and checks

𝐽 (Θ + ℎ 𝐸 𝑘, 𝑗) - 𝐽 (Θ - ℎ 𝐸 𝑘, 𝑗) 2 ℎ \approx 𝜕 𝐽 𝜕 𝜃 𝑘, 𝑗,

where $𝐸 𝑘, 𝑗$ changes just that parameter. This is a symmetric difference quotient. It is a useful independent check of the formula implemented by the program, although it would be an expensive way to compute the training gradient at every step.

On the supplied check, using $ℎ = 10 - 5$ , the largest absolute discrepancy is about $2.4 \times 10 - 10$ . A discrepancy that small is consistent with the differentiation and arithmetic being right. It does not tell us whether a linear reader is an adequate model of every style of handwriting; that is a different experiment.

32.5Turn the gradient into a new reader

The essential update is one line:

candidate = theta - trial * grad

Before accepting it, the program evaluates its objective and applies the decrease test from the previous chapter. In the following code, loss_and_gradient(..., need_gradient=False) computes only the objective:

trial = min(2 * eta, 0.05)
for attempt in range(40):
    candidate = theta - trial * grad
    new_loss = loss_and_gradient(
        candidate, images, labels, lam, need_gradient=False)
    if new_loss <= loss - 1e-4 * trial * np.sum(grad * grad):
        break
    trial /= 2
else:
    raise ArithmeticError("No acceptable step found")
theta, eta = candidate, trial

The previous accepted size is eta. Trying twice that size gives the method a chance to increase its steps again; the cap limits the initial proposal. A rejected trial costs another objective evaluation, but it does not require a new gradient. Every trial starts from the same current parameters and uses the same current gradient.

We repeat this procedure, stopping after 4,000 accepted updates or earlier if the gradient norm falls below $10 - 5$ . The iteration limit is a computational budget, not a proof that we have reached the exact minimum. The program reports the final gradient norm so we can see how close the run came to its stopping tolerance.

32.6Run the file

The complete Python program contains the calculations above, the data loader, derivative checks, training loop, evaluation, and export. Download it together with the all-digit CSV data, put them in one folder, and run:

python3 -m pip install numpy
python3 train_reader.py --data optdigits-all.csv --output my-reader.json

From this book's repository root, the equivalent command is:

python3 book/digit-reader-laboratory/building-and-running-our-reader/code/train_reader.py

The CSV contains the original training/test designation and 64 block counts for each image. The loader divides each count by 16, prepends a one for the bias, and keeps the two groups separate. Its labels are the digits zero through nine. No images are shuffled between the official groups.

The images come from E. Alpaydin and C. Kaynak's Optical Recognition of Handwritten Digits data in the UCI Machine Learning Repository. The data notes preserve the source citation, license, extraction procedure, and relationship between the block counts and original bitmaps.

A run of this file gives the following results:

Quantity	Result
Accepted updates	4,000
Final objective $𝐽$	$593.034582$
Final gradient norm	approximately $3.52 \times 10 - 4$
Correct training answers	$3748 / 3823$
Correct test answers	$1705 / 1797$
Test accuracy	approximately $94.88 %$

This run stopped at the iteration limit, before reaching the requested gradient tolerance. The objective agrees closely with that of the demonstration reader, and the two runs make the same number of correct training and test predictions. Small numerical differences across machines are possible.

32.7Look at the mistakes

An accuracy of $94.88 %$ leaves 92 test images answered incorrectly. That is enough success to see the model working, and enough failure to investigate its limits.

A confusion table records the mistakes by digit. Row $𝑦$ , column $𝑘$ counts images written as $𝑦$ that the reader calls $𝑘$ . The diagonal contains correct answers. The other cells tell us which kinds of confusion occurred; clicking a cell lets us inspect the images rather than speculate from a single percentage.

Figure 32.2 These counts are computed from the parameter file currently loaded in this chapter. Choose a cell to inspect its test images and all ten probabilities. You can send an example to the drawing-pad figure below. Loading your own trained file there also recomputes this table.

A confidently wrong answer is particularly revealing. Softmax compares the scores this reader knows how to produce; it does not check whether its geometric model is appropriate for the image. Giving nearly all the probability to a five does not make a handwritten three into a five. It tells us how decisively this parameter list prefers one region to the others.

The test set is useful because it did not enter the gradient calculations. If we repeatedly change our model after studying its test mistakes, those images begin to influence our choices. For a larger project we would use a separate validation set for that experimentation and reserve a final test set. Here we use the fixed recipe above and inspect what it produced.

32.8Put your parameters inside the machine

The program saves my-reader.json. Its central entry is the list of 650 parameters, in digit order, with the bias first in each group. The file also records the coordinate convention so the browser can interpret the numbers correctly. There is no hidden training operation when we load it. The drawing pad simply replaces its ten lists and runs the dot products again.

Figure 32.3 The default parameters were produced by the Python file in this chapter. Load your own my-reader.json, try the stored images, and draw a digit. The active filename is displayed above the pad. Both the pad and the test-results table use that file's parameters.

New drawings pass through the same kind of measurement process used in Chapter 1: normalize the stroke image, count ink in blocks, divide by 16, and flatten the grid. The preprocessing matters. A reader trained on centered, size-normalized digits cannot be expected to treat an arbitrarily positioned mark as the same input. A freehand drawing on this pad may still differ from the handwriting in the dataset. Try changing the size and shape of a digit, and watch which changes the reader tolerates.

We have arrived back at the machine we began with, but now its numbers have a history. Each training image contributed to a scalar function on the space of readers. The gradient told us how that function changed. Repeated vector updates moved us to a parameter list that assigns useful regions to the ten digits. To read a new image, all that training becomes ten dot products.

32.9Further experiments

Train for 250 or 1,000 updates by adding --iterations 250 or --iterations 1000 to the command. Save each run under a different filename. Compare their objectives and training behavior first, then load the files to see the resulting readers. Increasing the number of updates need not improve test accuracy at every stage.

Another extension is to approximate the full gradient with a small randomly chosen group of training images. If a batch contains $𝐵$ images sampled uniformly, then

̂ \nabla 𝐽 = 𝑁 𝐵 \sum 𝑖 i n b a t c h \nabla ℓ 𝑖 + 𝜆 Θ

has the full gradient as its expectation. The factor $𝑁 / 𝐵$ matters because our objective uses a sum. Each estimate is cheaper but fluctuates from batch to batch; the full-objective decrease test and its guarantees cannot simply be transferred to these noisy updates. This is a starting point for stochastic methods, beyond the complete implementation here.

We could also try reading several separated digits in order. If their image boxes have already been supplied, we can apply the same reader to each box and concatenate its answers. Finding those boxes in a connected handwritten string is an additional problem. Our ten-way classifier answers “which digit is in this image?”; it has not learned where one digit ends and the next begins.