6 - A convolutional neural network from scratch

In chapter 6 Emma explains to the Minermont team the network that could help the dermatologists prioritize skin lesions: a convolutional neural network (CNN) that receives a 28×28 grayscale image and gives a probability for each class, benign, suspicious or malignant. The notebook Neural networks with TensorFlow and Keras builds it with a few lines of Keras. Here we write it from scratch, with NumPy only, to see what those lines hide:

1. Import the libraries.
2. The convolution by hand, with the book's example.
3. Padding, and the convolution as a matrix product.
4. Max-pooling.
5. The layers and their backward pass.
6. The network of the book: 105,411 parameters.
7. Checking the gradients.
8. The dermatology images.
9. Training with mini-batch gradient descent.
10. Evaluation on the test images.
11. What the filters have learned.

Step 1: Import the libraries

---

The network uses NumPy only. scikit-learn splits the data and draws the confusion matrix, as in the other notebooks.

import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, ConfusionMatrixDisplay

Step 2: The convolution by hand

---

Emma's first figure: a 3×3 kernel slides over a 5×5 input, one pixel at a time (stride 1). At each position it multiplies the window by the kernel, element by element, and adds up the nine products. That sum is one value of the feature map. With no padding, a 3×3 window fits in 3×3 positions of a 5×5 input.

book_input = np.array([[1, 3, 0, 2, 1],
                       [0, 2, 1, 3, 0],
                       [2, 1, 3, 1, 2],
                       [1, 0, 2, 0, 1],
                       [3, 1, 0, 2, 3]])
book_kernel = np.array([[0, 1, 0],
                        [1, 1, 1],
                        [0, 1, 0]])

def convolve(image, kernel):
    """Slide the kernel over the image (no padding, stride 1): one output per position."""
    k = kernel.shape[0]
    rows, cols = image.shape[0] - k + 1, image.shape[1] - k + 1
    feature_map = np.zeros((rows, cols))
    for i in range(rows):
        for j in range(cols):
            window = image[i:i + k, j:j + k]
            feature_map[i, j] = np.sum(window * kernel)
    return feature_map

print("First window:\n", book_input[:3, :3])
print("Window x kernel, element by element:\n", book_input[:3, :3] * book_kernel)
print("Sum:", np.sum(book_input[:3, :3] * book_kernel))
print("\nFeature map (3 x 3):\n", convolve(book_input, book_kernel))

First window:
 [[1 3 0]
 [0 2 1]
 [2 1 3]]
Window x kernel, element by element:
 [[0 3 0]
 [0 2 1]
 [0 1 0]]
Sum: 7

Feature map (3 x 3):
 [[7. 9. 7.]
 [8. 8. 9.]
 [5. 5. 6.]]

The first window gives 3 + 2 + 1 + 1 = 7, and the whole feature map is the one in Emma's figure. This kernel is a cross: it adds up each pixel and its four neighbours, so each value of the map says how much brightness there is in that cross-shaped neighbourhood of the input. A network uses the same operation, but it does not choose the kernels: it learns them.

Step 3: Padding, and the convolution as a matrix product

---

The book's network uses padding='same': a border of zeros around the image, so that the feature map keeps the size of the input (28×28 stays 28×28). And a real layer does not apply one kernel with a double loop, but many kernels to many images. The usual trick, called im2col, copies every 3×3 window into a row of a matrix. Then the whole layer, all positions, all images and all filters, is a single matrix product: (windows × 9·channels) times (9·channels × filters). We write it for images with any number of channels, stored as (images, height, width, channels), like Keras.

def im2col(x, k=3):
    """x: (n, h, w, c) -> one row per output position with its k*k*c input values ('same' padding)."""
    n, h, w, c = x.shape
    p = k // 2
    padded = np.pad(x, ((0, 0), (p, p), (p, p), (0, 0)))      # a border of zeros
    cols = np.empty((n, h, w, k, k, c), dtype=x.dtype)
    for i in range(k):
        for j in range(k):
            cols[:, :, :, i, j, :] = padded[:, i:i + h, j:j + w, :]
    return cols.reshape(n * h * w, k * k * c)

x = book_input[None, :, :, None].astype(float)                 # 1 image, 5 x 5, 1 channel
rows = im2col(x)
print(f"im2col matrix: {rows.shape[0]} windows x {rows.shape[1]} values")
same = (rows @ book_kernel.reshape(9, 1)).reshape(5, 5)
print("Feature map with 'same' padding (5 x 5):\n", same)
print("Its centre equals the feature map without padding:", np.array_equal(same[1:4, 1:4], convolve(book_input, book_kernel)))

im2col matrix: 25 windows x 9 values
Feature map with 'same' padding (5 x 5):
 [[4. 6. 6. 6. 3.]
 [5. 7. 9. 7. 6.]
 [4. 8. 8. 9. 4.]
 [6. 5. 5. 6. 6.]
 [5. 4. 5. 5. 6.]]
Its centre equals the feature map without padding: True

With a border of zeros the map has 25 values, one per pixel, and its 3×3 centre is exactly the map of step 2. The values on the edge are smaller, because part of the cross falls on the zeros. im2col turned the 25 windows into a 25×9 matrix, and a single matrix product computed all of them. In the first layer of the network the same product handles 32 images × 784 positions and 16 filters at once.

Step 4: Max-pooling

---

After each convolution the book's network applies a 2×2 max-pooling with stride 2: it scans each map with a 2×2 window that jumps 2 pixels at a time, never overlapping, and keeps only the largest value. The map is halved in each direction. We apply it to the top-left 4×4 corner of the feature map above, and we also remember where each maximum was, because the backward pass will need it.

def pool_forward(x):
    """2 x 2 max-pooling with stride 2. Returns the pooled maps and the position of each maximum (0 to 3)."""
    n, h, w, c = x.shape
    windows = x.reshape(n, h // 2, 2, w // 2, 2, c).transpose(0, 1, 3, 5, 2, 4).reshape(n, h // 2, w // 2, c, 4)
    where = windows.argmax(axis=-1)
    return np.take_along_axis(windows, where[..., None], axis=-1)[..., 0], where

corner = same[:4, :4]
pooled, where = pool_forward(corner[None, :, :, None])
print("4 x 4 corner of the feature map:\n", corner)
print("After 2 x 2 max-pooling:\n", pooled[0, :, :, 0])

4 x 4 corner of the feature map:
 [[4. 6. 6. 6.]
 [5. 7. 9. 7.]
 [4. 8. 8. 9.]
 [6. 5. 5. 6.]]
After 2 x 2 max-pooling:
 [[7. 9.]
 [8. 9.]]

Each 2×2 block keeps its maximum, 7, 9, 8 and 9, and the 4×4 map becomes 2×2: a quarter of the values, keeping the strongest response of each block. If the pattern moves by one pixel inside a block, the result does not change.

Step 5: The layers and their backward pass

---

To train, each layer needs two functions. The forward pass computes its output. The backward pass receives the gradient of the loss with respect to that output and returns the gradient with respect to its input and to its weights: it is backpropagation, the chain rule applied layer by layer from the loss to the image.

• Convolution. The forward pass is the matrix product of step 3 plus a bias per filter. Backwards, the gradient of the weights is the transposed im2col matrix times the incoming gradient, and the gradient of the input goes back through the same windows (col2im), adding up where windows overlap.
• ReLU, $\max(0, z)$, lets the gradient through where $z > 0$ and blocks it elsewhere.
• Max-pooling has no weights. Its gradient goes entirely to the position that held the maximum, and the other three get zero.
• Dense layer, $\mathbf{a}\,W + \mathbf{b}$, as in the perceptron notebook.
• Softmax with categorical cross-entropy. As in the multiclass notebook of chapter 4, the gradient with respect to the scores is the probability minus the one-hot label, $\hat{\mathbf{y}} - \mathbf{y}$.

def col2im(d_rows, shape, k=3):
    """Backward of im2col: send each window's gradient back to its pixels, adding up the overlaps."""
    n, h, w, c = shape
    p = k // 2
    d_rows = d_rows.reshape(n, h, w, k, k, c)
    d_padded = np.zeros((n, h + 2 * p, w + 2 * p, c), dtype=d_rows.dtype)
    for i in range(k):
        for j in range(k):
            d_padded[:, i:i + h, j:j + w, :] += d_rows[:, :, :, i, j, :]
    return d_padded[:, p:-p, p:-p, :]

def conv_forward(x, W, b):
    n, h, w, _ = x.shape
    rows = im2col(x)
    return (rows @ W.reshape(-1, W.shape[-1]) + b).reshape(n, h, w, -1), rows

def conv_backward(d_out, rows, x_shape, W):
    d = d_out.reshape(-1, W.shape[-1])
    dW = (rows.T @ d).reshape(W.shape)
    db = d.sum(axis=0)
    dx = col2im(d @ W.reshape(-1, W.shape[-1]).T, x_shape)
    return dx, dW, db

def pool_backward(d_out, where, shape):
    n, h, w, c = shape
    d_windows = np.zeros((n, h // 2, w // 2, c, 4), dtype=d_out.dtype)
    np.put_along_axis(d_windows, where[..., None], d_out[..., None], axis=-1)   # only the maximum gets the gradient
    return d_windows.reshape(n, h // 2, w // 2, c, 2, 2).transpose(0, 1, 4, 2, 5, 3).reshape(n, h, w, c)

def softmax(scores):
    scores = scores - scores.max(axis=1, keepdims=True)
    e = np.exp(scores)
    return e / e.sum(axis=1, keepdims=True)

def cross_entropy(probs, y):
    return -np.mean(np.log(probs[np.arange(len(y)), y] + 1e-12))

Step 6: The network of the book: 105,411 parameters

---

Emma's eight steps: input 28×28×1, convolution with 16 filters of 3×3 and ReLU, max-pooling to 14×14, convolution with 32 filters of 3×3 and ReLU, max-pooling to 7×7, flatten into 1,568 numbers, dense layer of 64 neurons with ReLU, and 3 outputs with softmax. The weights start at random values with the scale proposed by He for ReLU layers (standard deviation $\sqrt{2 / \text{inputs}}$), and the biases at zero. forward keeps everything the backward pass will need.

def init_network(seed=0, dtype=np.float32):
    rng = np.random.default_rng(seed)
    def he(shape, inputs):
        return rng.normal(0, np.sqrt(2 / inputs), shape).astype(dtype)
    return {'W1': he((3, 3, 1, 16), 9),       'b1': np.zeros(16, dtype),   # conv 1: 16 filters 3 x 3 x 1
            'W2': he((3, 3, 16, 32), 144),    'b2': np.zeros(32, dtype),   # conv 2: 32 filters 3 x 3 x 16
            'W3': he((1568, 64), 1568),       'b3': np.zeros(64, dtype),   # dense: 7 * 7 * 32 = 1568 inputs
            'W4': he((64, 3), 64),            'b4': np.zeros(3, dtype)}    # output: 3 classes

def forward(net, x):
    z1, rows1 = conv_forward(x, net['W1'], net['b1']); a1 = np.maximum(z1, 0); m1, where1 = pool_forward(a1)
    z2, rows2 = conv_forward(m1, net['W2'], net['b2']); a2 = np.maximum(z2, 0); m2, where2 = pool_forward(a2)
    flat = m2.reshape(len(x), -1)                                             # flatten
    z3 = flat @ net['W3'] + net['b3']; a3 = np.maximum(z3, 0)
    probs = softmax(a3 @ net['W4'] + net['b4'])
    return probs, (x, z1, rows1, a1, m1, where1, z2, rows2, a2, where2, m2, flat, z3, a3)

def backward(net, cache, probs, y):
    x, z1, rows1, a1, m1, where1, z2, rows2, a2, where2, m2, flat, z3, a3 = cache
    grads = {}
    d4 = probs.copy()
    d4[np.arange(len(y)), y] -= 1                                            # probability minus one-hot label
    d4 /= len(y)
    grads['W4'], grads['b4'] = a3.T @ d4, d4.sum(axis=0)
    d3 = (d4 @ net['W4'].T) * (z3 > 0)                                       # back through ReLU
    grads['W3'], grads['b3'] = flat.T @ d3, d3.sum(axis=0)
    d_m2 = (d3 @ net['W3'].T).reshape(m2.shape)                              # un-flatten
    d_z2 = pool_backward(d_m2, where2, a2.shape) * (z2 > 0)
    d_m1, grads['W2'], grads['b2'] = conv_backward(d_z2, rows2, m1.shape, net['W2'])
    d_z1 = pool_backward(d_m1, where1, a1.shape) * (z1 > 0)
    _, grads['W1'], grads['b1'] = conv_backward(d_z1, rows1, x.shape, net['W1'])
    return grads

net = init_network()
layers = [('Conv 1 (16 filters 3x3)', ['W1', 'b1']), ('MaxPooling 1', []), ('Conv 2 (32 filters 3x3)', ['W2', 'b2']),
          ('MaxPooling 2', []), ('Flatten', []), ('Dense (64)', ['W3', 'b3']), ('Output (3, softmax)', ['W4', 'b4'])]
for name, keys in layers:
    print(f"{name:<24}{sum(net[k].size for k in keys):>8,}")
print(f"{'Total':<24}{sum(v.size for v in net.values()):>8,}")

Conv 1 (16 filters 3x3)      160
MaxPooling 1                   0
Conv 2 (32 filters 3x3)    4,640
MaxPooling 2                   0
Flatten                        0
Dense (64)               100,416
Output (3, softmax)          195
Total                    105,411

The same 160, 4,640, 100,416 and 195 as in the book's table, 105,411 in total. Weight sharing is visible in these numbers: the first convolution looks at the whole image with only 160 parameters, because the same 16 small filters are reused at every position, while the dense layer needs 100,416, because each of its 1,568 inputs is connected to each of its 64 neurons.

Step 7: Checking the gradients

---

A mistake in a backward pass does not produce an error message: the network simply learns badly. So before training we compare each analytic gradient with a numerical one, obtained by nudging a single weight up and down by $\varepsilon$ and measuring how the loss changes: $\frac{L(w + \varepsilon) - L(w - \varepsilon)}{2\varepsilon}$. We use double precision, four images and five random weights of each array. The biases are set to small random values for the check: at exactly zero, some units sit exactly on the kink of the ReLU, where the derivative is not defined and the two methods can legitimately disagree.

check_rng = np.random.default_rng(1)
check_net = init_network(seed=1, dtype=np.float64)
for key in ['b1', 'b2', 'b3', 'b4']:
    check_net[key] = check_rng.normal(0, 0.1, check_net[key].shape)
check_x = check_rng.random((4, 28, 28, 1))
check_y = np.array([0, 1, 2, 1])

probs, cache = forward(check_net, check_x)
analytic = backward(check_net, cache, probs, check_y)
eps, worst = 1e-6, 0.0
for key, W in check_net.items():
    for _ in range(5):
        idx = tuple(check_rng.integers(0, s) for s in W.shape)
        original = W[idx]
        W[idx] = original + eps; loss_plus = cross_entropy(forward(check_net, check_x)[0], check_y)
        W[idx] = original - eps; loss_minus = cross_entropy(forward(check_net, check_x)[0], check_y)
        W[idx] = original
        numerical = (loss_plus - loss_minus) / (2 * eps)
        relative = abs(numerical - analytic[key][idx]) / max(abs(numerical) + abs(analytic[key][idx]), 1e-12)
        worst = max(worst, relative)
print(f"Largest relative difference between analytic and numerical gradients: {worst:.1e}")

Largest relative difference between analytic and numerical gradients: 9.5e-08

The largest relative difference is 9.5e-08, the size of rounding errors: the backward passes are right. An error in the backward pass of a layer would show up in the gradients of its weights and of all the layers before it, which are computed from it.

Step 8: The dermatology images

---

The same synthetic lesions as in the Keras notebook, with the same generator, seed and split, so that the results can be compared: 3,600 images of 28×28 pixels, 75% for training and 25% for the test. As validation_split=0.2 does in Keras, the last 20% of the training images are kept apart for validation.

def make_lesions(n, seed, shares=(0.5, 0.3, 0.2)):
    """Synthetic 28x28 grayscale images of skin lesions: 0 = benign, 1 = suspicious, 2 = malignant."""
    rng = np.random.default_rng(seed)
    labels = rng.choice(3, n, p=shares)
    yy, xx = np.mgrid[0:28, 0:28].astype(float)
    images = np.empty((n, 28, 28), dtype=np.float32)
    # the hidden rules: size (r), border irregularity (irr), heterogeneous pigment (het),
    # asymmetric extra lobe (lobe) and how blurred the border is (soft); the ranges overlap on purpose
    params = {0: dict(r=(4.0, 7.0), irr=(0.0, 0.10), het=(0.0, 0.06), lobe=0.05, soft=(1.2, 2.0)),
              1: dict(r=(5.0, 8.5), irr=(0.07, 0.20), het=(0.04, 0.12), lobe=0.35, soft=(0.8, 1.6)),
              2: dict(r=(6.0, 10.0), irr=(0.15, 0.32), het=(0.10, 0.22), lobe=0.75, soft=(0.5, 1.2))}
    for i, c in enumerate(labels):
        p = params[c]
        cx, cy = rng.uniform(11, 17, 2)
        r0, irr = rng.uniform(*p['r']), rng.uniform(*p['irr'])
        het, soft = rng.uniform(*p['het']), rng.uniform(*p['soft'])
        theta, dist = np.arctan2(yy - cy, xx - cx), np.hypot(xx - cx, yy - cy)
        radius = r0 * (1 + sum(irr / k ** 0.5 * rng.uniform(0.5, 1.0) * np.cos(k * theta + rng.uniform(0, 2 * np.pi))
                               for k in range(2, 7)))
        mask = 1 / (1 + np.exp(-(radius - dist) / soft))
        if rng.random() < p['lobe']:
            a = rng.uniform(0, 2 * np.pi)
            lx, ly = cx + 0.8 * r0 * np.cos(a), cy + 0.8 * r0 * np.sin(a)
            mask = np.maximum(mask, 1 / (1 + np.exp(-(0.55 * r0 - np.hypot(xx - lx, yy - ly)) / soft)))
        pigment = rng.uniform(0.45, 0.6) + het * rng.normal(0, 1, (28, 28))
        for _ in range(3):                                   # darker and lighter blotches
            bx, by = cx + rng.normal(0, r0 / 2), cy + rng.normal(0, r0 / 2)
            pigment += het * 2.5 * rng.choice([-1, 1]) * np.exp(-((xx - bx) ** 2 + (yy - by) ** 2) / (2 * (r0 / 3) ** 2))
        skin = 0.82 + 0.04 * (xx - 14) / 14 * rng.normal() + 0.03 * rng.normal(0, 1, (28, 28))
        images[i] = np.clip(skin - mask * pigment, 0, 1)
    return images, labels

class_names = ['benign', 'suspicious', 'malignant']
images, labels = make_lesions(3600, seed=0)
X_train, X_test, y_train, y_test = train_test_split(images[..., None], labels, test_size=0.25, random_state=42, stratify=labels)
split = int(len(y_train) * 0.8)
X_fit, y_fit, X_val, y_val = X_train[:split], y_train[:split], X_train[split:], y_train[split:]
print(f"Training: {len(y_fit)} images, validation: {len(y_val)}, test: {len(y_test)}")

fig, axes = plt.subplots(3, 6, figsize=(12, 6.5))
for c in range(3):
    for ax, idx in zip(axes[c], np.where(y_fit == c)[0][:6]):
        ax.imshow(X_fit[idx, :, :, 0], cmap='gray', vmin=0, vmax=1)
        ax.axis('off')
    axes[c, 0].set_title(class_names[c], loc='left')
plt.suptitle('Synthetic skin lesions (28 x 28 pixels)')
plt.show()

Training: 2160 images, validation: 540, test: 900

Step 9: Training with mini-batch gradient descent

---

The Keras notebook uses the Adam optimizer. Here we use the algorithm of chapter 2, mini-batch gradient descent: batches of 32 images, shuffled at the start of every epoch, and the update $w := w - \alpha \nabla L$ for each of the 105,411 parameters. After each epoch we measure the loss on the validation images. If it has not improved, the learning rate is halved, as in the gradient-descent notebook of chapter 2; if it has not improved for 4 epochs in a row, training stops. At most 20 epochs, starting with $\alpha = 0.05$, and we keep the weights of the epoch with the lowest validation loss (early stopping).

def predict_probs(net, x, batch=500):
    return np.concatenate([forward(net, x[i:i + batch])[0] for i in range(0, len(x), batch)])

net = init_network(seed=0)
learning_rate, batch_size, patience = 0.05, 32, 4
order_rng = np.random.default_rng(0)
history = {'loss': [], 'val_loss': [], 'accuracy': [], 'val_accuracy': [], 'learning_rate': []}
best_loss, best_net, best_epoch, waited = np.inf, None, 0, 0
start = time.time()
for epoch in range(1, 21):
    order = order_rng.permutation(len(y_fit))
    for first in range(0, len(order), batch_size):
        batch = order[first:first + batch_size]
        probs, cache = forward(net, X_fit[batch])
        grads = backward(net, cache, probs, y_fit[batch])
        for key in net:
            net[key] -= learning_rate * grads[key]
    p_fit, p_val = predict_probs(net, X_fit), predict_probs(net, X_val)
    history['loss'].append(cross_entropy(p_fit, y_fit)); history['val_loss'].append(cross_entropy(p_val, y_val))
    history['accuracy'].append(np.mean(p_fit.argmax(1) == y_fit)); history['val_accuracy'].append(np.mean(p_val.argmax(1) == y_val))
    history['learning_rate'].append(learning_rate)
    if history['val_loss'][-1] < best_loss:
        best_loss, best_net, best_epoch, waited = history['val_loss'][-1], {k: v.copy() for k, v in net.items()}, epoch, 0
    else:
        waited += 1
        learning_rate /= 2
        if waited == patience:
            break
print(f"Epochs run: {epoch}, best epoch: {best_epoch}, training time: {time.time() - start:.0f} s")
print(f"Best epoch: training accuracy {history['accuracy'][best_epoch - 1]:.3f}, "
      f"validation accuracy {history['val_accuracy'][best_epoch - 1]:.3f}, validation loss {best_loss:.4f}")
print("Learning rate in each epoch:", ", ".join(f"{rate:g}" for rate in history['learning_rate']))

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for ax, key, title in [(axes[0], 'loss', 'Loss (categorical cross-entropy)'), (axes[1], 'accuracy', 'Accuracy')]:
    epochs = range(1, len(history[key]) + 1)
    ax.plot(epochs, history[key], 'o-', label='training')
    ax.plot(epochs, history['val_' + key], 's-', label='validation')
    ax.axvline(best_epoch, color='gray', ls='--', label='best epoch')
    ax.set_xticks([1, 5, 10, 15, 20]); ax.set_xlabel('epoch'); ax.set_title(title); ax.legend(); ax.grid(True)
plt.show()
net = best_net

Epochs run: 20, best epoch: 20, training time: 20 s
Best epoch: training accuracy 0.948, validation accuracy 0.931, validation loss 0.1876
Learning rate in each epoch: 0.05, 0.05, 0.05, 0.05, 0.05, 0.025, 0.025, 0.0125, 0.0125, 0.00625, 0.00625, 0.003125, 0.0015625, 0.0015625, 0.0015625, 0.0015625, 0.0015625, 0.00078125, 0.000390625, 0.000390625

The loss falls fast in the first epochs, with a couple of bumps in the validation loss, at epochs 5 and 7, while the learning rate is still large. Each time the validation loss does not improve, the rate is halved, and from epoch 10 on the curves settle. Early stopping never triggers: the validation loss keeps improving, very slightly, until the 20th epoch, which is the best one, with a training accuracy of 0.948, a validation accuracy of 0.931 and a validation loss of 0.1876. The two curves stay close, so the network has not memorized the training images. On a laptop the whole training takes seconds; with larger images and networks it would take much longer, which is why Emma recommends GPUs.

Step 10: Evaluation on the test images

---

We evaluate once on the 900 test images. As in the Keras notebook, the number that matters most is the recall of the malignant class.

test_probs = predict_probs(net, X_test)
y_pred = test_probs.argmax(axis=1)
print(f"Test accuracy: {np.mean(y_pred == y_test):.3f}\n")
print(classification_report(y_test, y_pred, target_names=class_names, digits=3))
ConfusionMatrixDisplay.from_predictions(y_test, y_pred, display_labels=class_names, cmap='Blues')
plt.xlabel('Predicted class'); plt.ylabel('Real class'); plt.title('CNN from scratch on the test images')
plt.show()

Test accuracy: 0.926

              precision    recall  f1-score   support

      benign      0.954     0.961     0.957       457
  suspicious      0.869     0.873     0.871       259
   malignant      0.933     0.913     0.923       184

    accuracy                          0.926       900
   macro avg      0.919     0.915     0.917       900
weighted avg      0.926     0.926     0.926       900

The test accuracy is 0.926. The malignant lesions have a recall of 0.913: 168 of the 184 are flagged as malignant and the other 16 as suspicious, so none is taken for benign, and no benign lesion is taken for malignant either. The suspicious class, in the middle, is the hardest, with a recall of 0.873. For this problem, 105,411 parameters trained with the plain gradient descent of chapter 2 are enough.

Step 11: What the filters have learned

---

Nobody programmed the filters: they started as random numbers. Below are the 16 kernels of the first layer after training, and the 16 feature maps they produce for a malignant test lesion (after the ReLU). Light pixels in a map mean a strong response.

malignant = np.where((y_test == 2) & (y_pred == 2))[0][0]
maps = np.maximum(conv_forward(X_test[malignant:malignant + 1], net['W1'], net['b1'])[0][0], 0)

fig, axes = plt.subplots(2, 17, figsize=(17, 2.6))
axes[0, 0].imshow(X_test[malignant, :, :, 0], cmap='gray', vmin=0, vmax=1)
axes[0, 0].set_title('image', fontsize=9)
axes[1, 0].axis('off')
for f in range(16):
    axes[0, f + 1].imshow(net['W1'][:, :, 0, f], cmap='gray')
    axes[1, f + 1].imshow(maps[:, :, f], cmap='gray')
for ax in axes.ravel():
    ax.set_xticks([]); ax.set_yticks([])
axes[0, 1].set_ylabel('kernel'); axes[1, 1].set_ylabel('map')
plt.suptitle('First convolutional layer: the 16 learned kernels and their feature maps')
plt.show()

The kernels do not look like anything a person would design, and the maps show why they are still useful: some highlight the lesion, others the skin around it or parts of its border, and a few are almost black for this image. The second layer combines these 16 maps into 32 new ones, and the dense layer turns them into the three probabilities. As Emma warns, this is the network's own way of seeing: what each filter responds to does not necessarily match our human idea of edges and textures.

Conclusion

---

• A convolution is a small kernel that slides over the image and adds up products. With im2col, a whole layer becomes one matrix product.
• Max-pooling halves the maps by keeping the largest value of each 2×2 window, and backwards it sends the gradient only to that value.
• Backpropagation is the chain rule, layer by layer. A numerical gradient check is the way to know that it is right before training.
• The book's network has 105,411 parameters, almost all of them, 100,416, in the dense layer. Trained with plain mini-batch gradient descent and a learning rate that is halved when the validation loss stalls, it classifies the synthetic lesions about as well as the Keras version.
• Nobody designs the filters: they come out of gradient descent, and, as Emma warns, what they detect does not necessarily match our human idea of edges and textures.