6 - Neural networks with TensorFlow and Keras

In the previous notebook we wrote a neural network by hand, derivatives included. That is what Emma calls level 0 in her notes (section 6.8): perfect for understanding, far too laborious for real networks. Here we climb the other levels of her "building":

• Level 1, automatic differentiation (TensorFlow): we write the operations and the library computes every gradient.
• Level 2, high-level API (Keras): we combine standard layers instead of individual operations.
• Level 3, predefined models (scikit-learn): a neural network as one more fit/predict model.

With Keras we build the convolutional network that Emma designs for the Minermont dermatology project (section 6.4), layer by layer and parameter by parameter.

1. Import the libraries.
2. Level 1: automatic differentiation, and Martha's example again.
3. Level 2: XOR with Keras.
4. The dermatology images.
5. The book's convolutional network.
6. Evaluation.
7. What the network sees.
8. Level 3: a scikit-learn network on the same images.

Step 1: Import the libraries

---

We fix the random seeds and ask TensorFlow for deterministic operations, so that the results can be reproduced. On a different machine or with a GPU the numbers may still vary slightly.

import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import keras
from keras import layers
from sklearn.model_selection import train_test_split
from sklearn.neural_network import MLPClassifier
from sklearn.metrics import accuracy_score, recall_score, classification_report, ConfusionMatrixDisplay

keras.utils.set_random_seed(42)
tf.config.experimental.enable_op_determinism()
np.set_printoptions(precision=3, suppress=True)
print(f"TensorFlow {tf.__version__}, Keras {keras.__version__}")

TensorFlow 2.20.0, Keras 3.13.2

Step 2: Level 1, automatic differentiation

---

TensorFlow represents every calculation as a computational graph: each node is a variable, an input or an operation. While a GradientTape is "recording", TensorFlow remembers the graph and can then apply the chain rule to it automatically. Emma's example is $y = 3x + 2$, whose derivative is 3:

x = tf.Variable(2.0)
with tf.GradientTape() as tape:
    y = 3 * x + 2
print(f"y = {y.numpy()}, dy/dx = {tape.gradient(y, x).numpy()}")

y = 8.0, dy/dx = 3.0

Now Martha's diagnostic network (section 6.7), the one we differentiated by hand in the previous notebook. We only write the forward pass; the gradients come from the tape:

x_in = tf.constant([0.5, 0.7])
W1 = tf.Variable([[0.8, 0.5], [0.3, 0.9], [1.2, 0.4]]);                b1 = tf.Variable([0.1, -0.2, 0.0])
W2 = tf.Variable([[1.1, 0.3, 0.5], [0.2, 1.0, 0.7], [0.6, 0.4, 1.3]]); b2 = tf.Variable([0.0, 0.1, -0.1])
y_true = tf.constant([0.0, 1.0, 0.0])                                   # flu

with tf.GradientTape() as tape:
    hidden = tf.nn.relu(tf.linalg.matvec(W1, x_in) + b1)
    probs = tf.nn.softmax(tf.linalg.matvec(W2, hidden) + b2)
    loss = -tf.reduce_sum(y_true * tf.math.log(probs))

grads = dict(zip(['W1', 'b1', 'W2', 'b2'], tape.gradient(loss, [W1, b1, W2, b2])))
print(f"Probabilities (cold, flu, pneumonia): {probs.numpy()}, loss {loss.numpy():.4f}")
for name, g in grads.items():
    print(f"dL/d{name} =\n{g.numpy()}")

Probabilities (cold, flu, pneumonia): [0.314 0.289 0.398], loss 1.2423
dL/dW1 =
[[ 0.221  0.309]
 [-0.229 -0.321]
 [ 0.088  0.123]]
dL/db1 =
[ 0.441 -0.458  0.176]
dL/dW2 =
[[ 0.267  0.182  0.276]
 [-0.605 -0.413 -0.626]
 [ 0.338  0.231  0.35 ]]
dL/db2 =
[ 0.314 -0.711  0.398]

The same probabilities (0.314, 0.289, 0.398), the same loss (1.2423) and the same gradients that we derived by hand in the previous notebook, for example −0.605 and −0.321, the two elements Martha checks. But here we did not write a single derivative: TensorFlow recorded the operations and applied the chain rule to them. That is what makes it possible to train networks of any architecture with millions of parameters.

Step 3: Level 2, XOR with Keras

---

With Keras we no longer write operations: we stack layers, choose a loss and an optimizer, and call fit. The XOR network of the previous notebook takes four lines. We use the Adam optimizer, a refined version of gradient descent that adapts the step of each weight.

X_xor = np.array([[0, 0], [0, 1], [1, 0], [1, 1]], dtype='float32')
y_xor = np.array([0, 1, 1, 0], dtype='float32')

xor_model = keras.Sequential([
    keras.Input(shape=(2,)),
    layers.Dense(4, activation='tanh'),
    layers.Dense(1, activation='sigmoid'),
])
xor_model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.05), loss='binary_crossentropy')
history = xor_model.fit(X_xor, y_xor, epochs=300, verbose=0)
print(f"Final loss: {history.history['loss'][-1]:.4f}")
print(f"Probabilities for (0,0), (0,1), (1,0), (1,1): {xor_model.predict(X_xor, verbose=0).ravel()}")

Final loss: 0.0029
Probabilities for (0,0), (0,1), (1,0), (1,1): [0.    0.995 0.998 0.005]

The network solves XOR: probabilities close to 0 for (0, 0) and (1, 1), and close to 1 for (0, 1) and (1, 0). Everything we wrote by hand in the previous notebook, the forward pass, backpropagation and the updates, is now inside fit.

Step 4: The dermatology images

---

In the book, Minermont's dermatologists need help to prioritize skin lesions, and defining by hand features like "symmetry" or "border irregularity" turns out to be too subjective. A real dataset of clinical images is confidential, so we generate synthetic images of 28×28 pixels in grayscale, the simplified format Emma uses to explain her network. Each image is a dark lesion on lighter skin. Following the usual warning signs, benign lesions tend to be small, round and uniform, and malignant ones larger, asymmetric, with an irregular border and a heterogeneous pigment; suspicious ones are in between. The rules overlap on purpose, as they do in reality. As in real clinics, benign lesions are the most frequent.

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)
print("Images per class:", dict(zip(class_names, np.bincount(labels).tolist())))

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

X_train, X_test, y_train, y_test = train_test_split(images[..., None], labels, test_size=0.25, random_state=42, stratify=labels)
Y_train = keras.utils.to_categorical(y_train, 3)                 # one-hot labels for the softmax output
print(f"Training images: {len(y_train)}, test images: {len(y_test)}, shape of each image: {X_train.shape[1:]}")

Images per class: {'benign': 1828, 'suspicious': 1036, 'malignant': 736}

Training images: 2700, test images: 900, shape of each image: (28, 28, 1)

Step 5: The book's convolutional network

---

Emma's network, exactly as she describes it in section 6.4:

1. Input of $28 \times 28 \times 1$ (one channel, grayscale).
2. Convolution with 16 filters of $3 \times 3$ and ReLU. Each filter slides over the image looking for a local pattern.
3. Max-pooling $2 \times 2$: keeps the largest value of each window and halves the size, to $14 \times 14$.
4. Convolution with 32 filters of $3 \times 3$, which combine the patterns of the first layer.
5. Max-pooling $2 \times 2$, down to $7 \times 7 \times 32 = 1568$ values.
6. Flatten into a list of 1568 numbers.
7. Dense layer of 64 neurons with ReLU.
8. Output of 3 neurons with softmax: benign, suspicious and malignant.

padding='same' keeps the size of the image in each convolution, as in the book. summary() lists the parameters of each layer, which we can compare with the book's table:

cnn = keras.Sequential([
    keras.Input(shape=(28, 28, 1)),
    layers.Conv2D(16, 3, padding='same', activation='relu'),
    layers.MaxPooling2D(2),
    layers.Conv2D(32, 3, padding='same', activation='relu'),
    layers.MaxPooling2D(2),
    layers.Flatten(),
    layers.Dense(64, activation='relu'),
    layers.Dense(3, activation='softmax'),
])
cnn.summary()
print(f"Total parameters: {cnn.count_params():,}")

Model: "sequential_1"

┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┓
┃ Layer (type)                    ┃ Output Shape           ┃       Param # ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━┩
│ conv2d (Conv2D)                 │ (None, 28, 28, 16)     │           160 │
├─────────────────────────────────┼────────────────────────┼───────────────┤
│ max_pooling2d (MaxPooling2D)    │ (None, 14, 14, 16)     │             0 │
├─────────────────────────────────┼────────────────────────┼───────────────┤
│ conv2d_1 (Conv2D)               │ (None, 14, 14, 32)     │         4,640 │
├─────────────────────────────────┼────────────────────────┼───────────────┤
│ max_pooling2d_1 (MaxPooling2D)  │ (None, 7, 7, 32)       │             0 │
├─────────────────────────────────┼────────────────────────┼───────────────┤
│ flatten (Flatten)               │ (None, 1568)           │             0 │
├─────────────────────────────────┼────────────────────────┼───────────────┤
│ dense_2 (Dense)                 │ (None, 64)             │       100,416 │
├─────────────────────────────────┼────────────────────────┼───────────────┤
│ dense_3 (Dense)                 │ (None, 3)              │           195 │
└─────────────────────────────────┴────────────────────────┴───────────────┘

 Total params: 105,411 (411.76 KB)

 Trainable params: 105,411 (411.76 KB)

 Non-trainable params: 0 (0.00 B)

Total parameters: 105,411

The parameters match the book's table layer by layer: 160 in the first convolution, 4,640 in the second, 100,416 in the dense layer and 195 in the output, 105,411 in total. Almost all of them are in the dense layer. The convolutions need very few, because each small filter is shared across the whole image.

We train with the categorical cross-entropy, as Emma proposes, and the Adam optimizer. 20% of the training images are set aside for validation: after each epoch (a full pass over the training images) we measure the loss on them, and early stopping halts the training when it stops improving for 4 epochs, keeping the best weights. It is the same idea as choosing the depth of a tree with validation data: stop before memorizing.

cnn.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy'])
early_stopping = keras.callbacks.EarlyStopping(monitor='val_loss', patience=4, restore_best_weights=True)
history = cnn.fit(X_train, Y_train, validation_split=0.2, epochs=40, batch_size=32, callbacks=[early_stopping], verbose=0)

h = history.history
best = int(np.argmin(h['val_loss']))
print(f"Epochs run: {len(h['loss'])}, best epoch: {best + 1}")
print(f"Best epoch: training accuracy {h['accuracy'][best]:.3f}, validation accuracy {h['val_accuracy'][best]:.3f}")

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')]:
    ax.plot(np.arange(1, len(h[key]) + 1), h[key], 'o-', label='training')
    ax.plot(np.arange(1, len(h[key]) + 1), h['val_' + key], 's-', label='validation')
    ax.axvline(best + 1, color='gray', linestyle='--')
    ax.set_xlabel('Epoch')
    ax.set_title(title)
    ax.grid(True)
    ax.legend()
plt.show()

Epochs run: 12, best epoch: 8
Best epoch: training accuracy 0.919, validation accuracy 0.900

Early stopping ended the training after 12 epochs and kept the weights of epoch 8, the one with the lowest validation loss: an accuracy of 0.919 on the training images and 0.900 on the validation ones. The two curves stay close, so the network has not memorized the training images.

Step 6: Evaluation

---

We evaluate once on the test images. In dermatology the critical number is the recall of the malignant class: how many malignant lesions the network flags.

probabilities = cnn.predict(X_test, verbose=0)
y_pred = probabilities.argmax(axis=1)
print(f"Test accuracy: {accuracy_score(y_test, y_pred):.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('Convolutional network on the test images')
plt.show()

Test accuracy: 0.928

              precision    recall  f1-score   support

      benign      0.917     0.996     0.955       457
  suspicious      0.918     0.822     0.868       259
   malignant      0.971     0.908     0.938       184

    accuracy                          0.928       900
   macro avg      0.935     0.909     0.920       900
weighted avg      0.929     0.928     0.926       900

The network classifies 92.8% of the test images correctly. It recognizes almost every benign lesion (recall 0.996) and 90.8% of the malignant ones, and when it says "malignant" it is right 97.1% of the time. The suspicious class, which lies between the other two, is the hardest (recall 0.822). In a real clinic, the roughly 9% of malignant lesions that are missed is the number that matters most, and it is why, as in the book, the system only helps to prioritize and the dermatologist always makes the decision.

The network does not only give a class: softmax gives a probability for each one, like the "benign 0.12, suspicious 0.21, malignant 0.67" of Emma's example. These are some test images with the three probabilities:

examples = [np.where(y_test == c)[0][k] for c in range(3) for k in range(2)]
fig, axes = plt.subplots(1, 6, figsize=(16, 3.6))
for ax, i in zip(axes, examples):
    ax.imshow(X_test[i, :, :, 0], cmap='gray', vmin=0, vmax=1)
    ax.axis('off')
    ax.set_title(f"real: {class_names[y_test[i]]}\n" + "\n".join(f"{n} {p:.2f}" for n, p in zip(class_names, probabilities[i])), fontsize=9)
plt.show()

Step 7: What the network sees

---

Ethan asks in the book how we know what the filters are really detecting. We can look at them directly: these are the 16 filters of $3 \times 3$ that the first convolutional layer has learned, and the 16 feature maps they produce for one malignant lesion (how strongly each filter responds at each position).

filters = cnn.layers[0].get_weights()[0][:, :, 0, :]              # 3 x 3 x 16
first_layer = keras.Model(cnn.inputs, cnn.layers[0].output)
image = X_test[np.where(y_test == 2)[0][0]]
maps = first_layer.predict(image[None], verbose=0)[0]             # 28 x 28 x 16

fig, axes = plt.subplots(2, 17, figsize=(19, 2.8))
for row in range(2):
    axes[row, 0].imshow(image[:, :, 0], cmap='gray', vmin=0, vmax=1)
for k in range(16):
    axes[0, k + 1].imshow(filters[:, :, k], cmap='bwr')
    axes[1, k + 1].imshow(maps[:, :, k], cmap='viridis')
for ax in axes.ravel():
    ax.axis('off')
axes[0, 0].set_title('image')
axes[0, 1].set_title('the 16 filters', loc='left')
axes[1, 1].set_title('their 16 feature maps', loc='left')
plt.show()

Nobody told the network what to look for. Some feature maps light up at the border of the lesion, others on the dark pigment, on the lighter skin around it or on the texture inside. As Emma warns, this is a simplified picture: the filters do not correspond exactly to human concepts such as "edge" or "texture", and the deeper layers combine them in ways that are much harder to interpret.

Step 8: Level 3, a scikit-learn network

---

At the top of Emma's building, scikit-learn's MLPClassifier is a neural network with the same interface as any other model. It only accepts dense layers, so each image has to be flattened into 784 numbers from the start, losing its spatial structure. We give it a hidden layer of 64 neurons, like the dense part of the CNN:

mlp = MLPClassifier(hidden_layer_sizes=(64,), early_stopping=True, max_iter=300, random_state=0)
mlp.fit(X_train.reshape(len(X_train), -1), y_train)
mlp_pred = mlp.predict(X_test.reshape(len(X_test), -1))
mlp_params = sum(w.size for w in mlp.coefs_) + sum(b.size for b in mlp.intercepts_)

for name, pred, n_params in [('CNN (Keras)', y_pred, cnn.count_params()), ('MLP (scikit-learn)', mlp_pred, mlp_params)]:
    rec = recall_score(y_test, pred, average=None)
    print(f"{name:<20} parameters {n_params:>7,} | accuracy {accuracy_score(y_test, pred):.3f} | "
          f"recall benign {rec[0]:.3f}, suspicious {rec[1]:.3f}, malignant {rec[2]:.3f}")

CNN (Keras)          parameters 105,411 | accuracy 0.928 | recall benign 0.996, suspicious 0.822, malignant 0.908
MLP (scikit-learn)   parameters  50,435 | accuracy 0.670 | recall benign 0.917, suspicious 0.456, malignant 0.359

The scikit-learn network, with 50,435 parameters, only reaches 0.670 and misses almost two out of three malignant lesions (recall 0.359). The CNN reaches 0.928. The difference is not the number of parameters, since most of the CNN's are in its dense layer, but the architecture. By keeping the image as a grid and reusing the same small filters everywhere, the CNN detects borders and textures wherever they are, which is exactly Emma's argument for convolutional networks.

Conclusion

---

• Automatic differentiation (TensorFlow) computes the same gradients we derived by hand, for any network we write: that is what made deep learning practical.
• Keras builds a network by stacking layers. Emma's CNN takes a dozen lines, and its 105,411 parameters are exactly those of the book's table.
• A convolutional network keeps the spatial structure of the image: small filters, shared across the whole image, detect local patterns such as edges and textures, and pooling makes them robust to small shifts. On our images it clearly beats a dense network that sees the pixels as a flat list.
• Validation data and early stopping control overfitting, and the evaluation must look at each class, above all at the malignant lesions that the network misses.
• These images are synthetic and very simple. As the book insists, a real system needs large, representative and well-labeled proprietary data, careful validation, and dermatologists who always make the final decision.