Build a Two-Layer MNIST Network

Build a NumPy MNIST classifier with manual backpropagation and gradient checks.

Introduction

30 Second Summary

A handwritten 7 is obvious to you, but a computer sees only a grid of light and dark squares. Its first guesses are little better than chance.

In this project, you will build a two-hidden-layer MNIST digit classifier in Google Colab using NumPy. You will implement manual backpropagation so every gradient that teaches the network remains visible.

What You'll Build

You will run one notebook that turns unreliable guesses into a falling loss curve, a measured test score, and a grid showing each predicted digit beside its true label.

By the end of this project, you'll have:

  • An untrained prediction demo that passes eight handwritten digits through a 784 -> 128 -> 64 -> 10 network to expose the unreliable guesses produced by random weights.
  • A set of inspectable gradients for W1, b1, W2, b2, W3, and b3, plus a one-batch update that visibly lowers the loss.
  • A trained classifier report that prints eight epoch summaries, plots mean squared error, reports held-out accuracy, and displays ten labeled predictions.
  • Secret Mission: Use finite differences to compare one analytical W3 gradient with an independent numerical estimate.

Are there any prerequisites?

You should be comfortable writing Python and performing basic NumPy array operations.

The project runs on Google Colab's free standard hosted runtime, so you do not need a local Python setup or paid compute.

Before We Start

Before the matrix work begins, this is your moment to commit to a handwritten digit classifier and your reason for implementing its gradients manually. Your blank Google Colab notebook gives that goal a clear starting point.

Pin and Verify the Colab Environment

Every gradient in this project depends on consistent numerical behavior. A shifting Google Colab runtime can make a package mismatch look like a bug in your network.

Pinning runtime 2026.07 fixes Python at 3.12.13 with NumPy 2.0.2. The package cell fixes scikit-learn at 1.9.0 plus Matplotlib at 3.11.1.

In this step, get ready to:
  • Pin the hosted runtime to the project version.
  • Install the required library versions.
  • Create the initial network configuration in a verified environment.
Pin the hosted runtime

The runtime version controls the Python base plus the preinstalled NumPy release. Switching now keeps the later gradient results reproducible.

A runtime change starts a fresh hosted machine. Your notebook is blank, so there is no runtime state to lose.

  • Click Runtime in the Google Colab menu bar.
  • Click Change runtime type.
  • Select 2026.07 from the Runtime Version dropdown.
  • Keep the standard hosted CPU runtime selected.
  • Confirm the change using the button in the lower-right corner.

You'll return to the notebook while Google Colab connects to the selected runtime. The active virtual machine now uses the project's pinned foundation.

Install the pinned packages

Runtime 2026.07 already supplies NumPy 2.0.2. The installation cell aligns the loader plus plotting libraries with the code you build next.

  • Add the first code cell using the notebook toolbar.
  • Use the first cell to apply both package pins by running this command:
%pip install scikit-learn==1.9.0 matplotlib==3.11.1

What does this command do?

  • The %pip magic installs packages inside the active notebook kernel.
  • The first pin installs scikit-learn 1.9.0.
  • The second pin installs Matplotlib 3.11.1.
  • Wait for the installation cell to finish.

You'll see the cell finish without an installation error. The pinned packages are now available to the active kernel.

Installation cell not completing?

  • Confirm that Google Colab still shows an active hosted runtime connection.
  • Rerun the installation cell if the runtime connection was interrupted.
  • Use the installation output when you ask for help diagnosing the package installation.
Add imports and configuration

The first project cell defines the libraries plus dimensions used throughout the classifier. Keeping these values together gives every later matrix shape one source of truth.

  • Add a code cell below the installation output using the notebook toolbar.
  • Paste the library imports into the new cell:
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml

What do these imports provide?

  • The NumPy import provides arrays plus vectorized matrix operations.
  • The Matplotlib import provides the loss curve plus digit grids.
  • The fetch_openml import provides the fixed MNIST dataset loader used in the next step.
  • Run the import cell.

You'll see the cell complete without an exception. All three project imports now load in the pinned runtime.

Imports not completing?

  • Rerun the installation cell if scikit-learn or Matplotlib is unavailable.
  • Confirm that each import matches the displayed code exactly.
  • Use the import traceback when you ask for help diagnosing the failed import.
  • Place the cursor on a blank line beneath the imports.
  • Paste the network configuration below the imports:
SEED = 42
INPUT_SIZE = 784
HIDDEN_1 = 128
HIDDEN_2 = 64
OUTPUT_SIZE = 10
TRAIN_SIZE = 12_000
TEST_SIZE = 2_000
BATCH_SIZE = 128
EPOCHS = 8
LEARNING_RATE = 0.1

What does this configuration control?

  • The SEED value makes later random sampling plus parameter initialization reproducible.
  • The layer sizes define the 784 -> 128 -> 64 -> 10 network shape.
  • The training constants control subset sizes, batch size, epoch count, plus update strength.
  • Run the updated project cell.

You'll see the cell complete without an exception. The imports plus configuration values now exist in the active kernel.

Configuration cell showing an error?

  • Confirm that every constant uses uppercase letters.
  • Check that the import lines remain above the constants.
  • Use the failing line when you ask for help checking the configuration cell.

✔️ Awesome, I've got everything!

Great. Double-check that the project cell has run successfully.

ⓧ I'd like to double check the full code

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml

SEED = 42
INPUT_SIZE = 784
HIDDEN_1 = 128
HIDDEN_2 = 64
OUTPUT_SIZE = 10
TRAIN_SIZE = 12_000
TEST_SIZE = 2_000
BATCH_SIZE = 128
EPOCHS = 8
LEARNING_RATE = 0.1

A successful project cell proves the imports load. The final diagnostic checks the exact versions behind those imports.

Before you run the diagnostic, do you expect every version to match its pin?

  • Add a temporary code cell below the project cell using the notebook toolbar.
  • Print the active language plus library versions by running this code:
import sys
import sklearn
import matplotlib

# Read version details from the active notebook kernel
print("Python", f"{sys.version_info.major}.{sys.version_info.minor}.{sys.version_info.micro}")
print("NumPy", np.__version__)
print("scikit-learn", sklearn.__version__)
print("Matplotlib", matplotlib.__version__)

What does this diagnostic check?

  • The Python line reads the interpreter version from the active runtime.
  • Each library line reads the version loaded by the current kernel.
  • The four printed values expose runtime drift before it reaches the network code.

You'll see Python 3.12.13. You'll also see NumPy 2.0.2, scikit-learn 1.9.0, plus Matplotlib 3.11.1.

✔️ I see the four pinned versions

That's the environment locked down. Your runtime now matches every version used by the project.

ⓧ Python or NumPy differs

A different Python or NumPy version means the notebook is using another runtime release.

  • Return to Runtime in the Google Colab menu bar.
  • Select Change runtime type.
  • Set Runtime Version to 2026.07.
  • Confirm the runtime change using the button in the lower-right corner.
  • Rerun the installation cell.
  • Rerun the temporary version-check cell.

Still seeing a different runtime?

Use the selected runtime version plus diagnostic output when you ask for help correcting the Colab runtime.

ⓧ A package is missing or different

The runtime is active, but one pinned package did not load at the expected version.

  • Rerun the package installation cell.
  • Wait for the installation cell to finish.
  • Rerun the project import cell.
  • Rerun the temporary version-check cell.

Package version still differs?

Use the installation output plus printed versions when you ask for help resolving the package mismatch.

Your notebook now has a reproducible runtime plus the network's starting configuration. Next, you'll load MNIST and watch the untrained network make its first unreliable predictions.

Build the Untrained Forward Pass

Your pinned Google Colab runtime now gives every matrix operation a stable base.

Random parameters have no learned digit patterns. Forward propagation lets you see that unreliable baseline before training changes a weight.

In this step, get ready to:
  • Fetch reproducible training and test subsets from MNIST.
  • Initialize the parameters for three dense layers.
  • Run eight test images through the untrained network.
Fetch and prepare MNIST

The MNIST dataset contains handwritten digit images with 784 pixel values per image. A fixed OpenML version keeps the input data consistent across runs.

  • Select the first blank code cell beneath LEARNING_RATE.
  • Define load_data(rng) by pasting this code:
def load_data(rng):
    X, y = fetch_openml(
        "mnist_784",
        version=1,
        return_X_y=True,
        as_frame=False,
    )

    order = rng.permutation(X.shape[0])
    selected = order[: TRAIN_SIZE + TEST_SIZE]
    X_selected = X[selected].astype(np.float32) / 255.0
    y_selected = y[selected].astype(np.int64)

    X_train = X_selected[:TRAIN_SIZE]
    y_train = y_selected[:TRAIN_SIZE]
    X_test = X_selected[TRAIN_SIZE:]
    y_test = y_selected[TRAIN_SIZE:]
    return X_train, y_train, X_test, y_test

What does this function do?

  • The fixed mnist_784 dataset version makes each fetch target the same dataset release.
  • The seeded permutation selects 14,000 images in a reproducible order.
  • Dividing the pixel values by 255.0 scales them before they enter the network.
  • The final slices create 12,000 training images and 2,000 test images.
  • Run the cell that defines load_data(rng).
  • Create a temporary verification cell beneath it.
  • Test the dataset download by running this code:
rng = np.random.default_rng(SEED)
X_train, y_train, X_test, y_test = load_data(rng)
print("Train shape:", X_train.shape, "Test shape:", X_test.shape)

What does this check prove?

  • The seeded rng makes the selected subset reproducible.
  • The printed shapes confirm that the download produced the required training and test arrays.
  • The image arrays use float32 values. The label arrays use int64 values.

You should see Train shape: (12000, 784) Test shape: (2000, 784). That successful fetch gives your classifier 14,000 normalized images to work with.

Dataset download not completing?

Confirm that the hosted runtime still has an internet connection. Reconnect the runtime if the download repeatedly stops.

Rerun the function cell before rerunning the temporary verification cell. Help me troubleshoot this OpenML dataset fetch.

Initialize the network and activations

Each dense layer needs a weight matrix plus a bias array. A seeded NumPy generator gives the network the same random starting point each time.

  • Select the next blank code cell.
  • Define the six parameter arrays by pasting this code:
def initialize_parameters(rng):
    # Scale each weight matrix by its input width for sigmoid activations
    return {
        "W1": rng.standard_normal((INPUT_SIZE, HIDDEN_1))
        / np.sqrt(INPUT_SIZE),
        "b1": np.zeros((1, HIDDEN_1)),
        "W2": rng.standard_normal((HIDDEN_1, HIDDEN_2))
        / np.sqrt(HIDDEN_1),
        "b2": np.zeros((1, HIDDEN_2)),
        "W3": rng.standard_normal((HIDDEN_2, OUTPUT_SIZE))
        / np.sqrt(HIDDEN_2),
        "b3": np.zeros((1, OUTPUT_SIZE)),
    }

How are the parameters organized?

  • The W1 matrix connects 784 input pixels to 128 first-layer activations.
  • The W2 matrix connects 128 first-layer activations to 64 second-layer activations.
  • The W3 matrix connects 64 second-layer activations to 10 digit outputs.
  • Dividing each random matrix by the square root of its input width helps keep sigmoid activations away from their flat extremes.
  • Run the parameter initialization cell.
  • Create another temporary verification cell beneath it.
  • Inspect the parameter shapes by running this code:
params = initialize_parameters(rng)
print("Parameter shapes:", {name: value.shape for name, value in params.items()})

What does this check reveal?

The dictionary keeps every trainable array under its matching name. The printed dimensions trace the full 784 -> 128 -> 64 -> 10 architecture.

You should see six entries named W1, b1, W2, b2, W3, and b3. Their dimensions should connect each neighboring layer.

Parameter shapes look different?

Confirm that the constants cell ran before the initialization cell. Check that each dimension uses the matching size constant.

Rerun the seeded generator cell if rng is missing. Help me diagnose these parameter shapes.

Each layer needs a nonlinear activation between its matrix operations. Sigmoid maps every value into the range from zero to one, which lets the network compare its ten outputs with one-hot targets.

  • Select the next blank code cell.
  • Define sigmoid(values) by pasting this code:
def sigmoid(values):
    # Keep exponentiation stable for large pre-activation values
    clipped = np.clip(values, -500.0, 500.0)
    negative_exponentials = np.exp(-clipped)
    return 1.0 / (1.0 + negative_exponentials)

How does sigmoid shape the outputs?

Sigmoid converts each pre-activation into a value between zero and one. Each output can therefore move toward the matching zero or one in a one-hot target.

Clipping limits the input to np.exp. This prevents overflow without changing ordinary activation values.

  • Run the cell that defines sigmoid(values).

The cell should complete without an exception. The sigmoid function is now available to all three layers.

Sigmoid cell raising an exception?

Check that the NumPy import cell ran earlier. Confirm that the return line remains indented inside sigmoid().

Compare the parentheses in the np.clip and np.exp calls. Help me fix this sigmoid function.

Run the untrained forward pass

The forward pass transforms one batch through three sigmoid layers before producing ten output activations. Its cache preserves the intermediate arrays that backpropagation needs later.

  • Select the next blank code cell.
  • Define forward(X, params) by pasting this code:
def forward(X, params):
    # Transform pixels through the first sigmoid hidden layer
    Z1 = X @ params["W1"] + params["b1"]
    A1 = sigmoid(Z1)

    # Build the second hidden representation
    Z2 = A1 @ params["W2"] + params["b2"]
    A2 = sigmoid(Z2)

    # Produce one bounded output for each digit class
    Z3 = A2 @ params["W3"] + params["b3"]
    outputs = sigmoid(Z3)

    cache = {
        "X": X,
        "A1": A1,
        "A2": A2,
        "outputs": outputs,
    }
    return outputs, cache

How does data cross the network?

  • The first matrix operation produces Z1. Sigmoid converts it into A1.
  • The second matrix operation produces Z2. A second sigmoid converts it into A2.
  • The output matrix produces Z3. Sigmoid converts each class score into a bounded output.
  • The cache keeps the input plus every activation required by the backward pass.
  • Run the cell that defines forward(X, params).
  • Create a temporary cell beneath it.
  • Confirm that eight images cross the full network by running this expression:
forward(X_test[:8], params)

What should this return?

The first returned array contains eight rows of ten sigmoid outputs. The second returned object is the cache of intermediate activations.

You should see a tuple containing the output array plus the cache. This proves that all three matrix operations connect without a shape error.

Seeing a matrix shape error?

Rerun the data and parameter verification cells so X_test plus params exist in memory. Compare the parameter dimensions with the architecture constants.

Check that each bias is added to the matrix from its matching layer. Help me trace this forward-pass shape error.

The classifier needs a numeric signal that rewards outputs close to the correct digit target. Mean squared error compares the ten sigmoid outputs with a one-hot row containing one for the true class and zero elsewhere.

  • Select the next blank code cell.
  • Define one_hot(y) and mean_squared_error(outputs, y) by pasting this code:
def one_hot(y):
    # Create one target column for each digit class
    targets = np.zeros((y.shape[0], OUTPUT_SIZE))
    targets[np.arange(y.shape[0]), y] = 1.0
    return targets


def mean_squared_error(outputs, y):
    # Average the squared gap across every image and class
    targets = one_hot(y)
    squared_errors = (outputs - targets) ** 2
    return np.mean(squared_errors)

What does the loss measure?

  • one_hot() creates one row of ten targets for each label.
  • The true class receives 1.0. Every other class stays at 0.0.
  • mean_squared_error() averages every squared output-to-target gap into one batch loss.
  • Run the cell that defines one_hot(y) and mean_squared_error(outputs, y).

The cell should complete without an exception. The network can now measure how far its untrained outputs are from the one-hot targets.

A small temporary plotting helper turns the raw predictions into images you can inspect. Matplotlib places each predicted label beside its true label.

  • Create a temporary plotting cell beneath the loss function.
  • Add the plotting helper by pasting this code:
def show_predictions(X, y, predictions, count):
    columns = count // 2
    _, axes = plt.subplots(2, columns)

    for index in range(count):
        ax = axes[index // columns, index % columns]
        ax.imshow(X[index].reshape(28, 28), cmap="gray")
        ax.set_title(f"P:{predictions[index]} T:{y[index]}")
        ax.axis("off")

    plt.show()

How does the grid help?

Each panel reshapes one row of 784 pixels into a 28 by 28 image. Its title places the predicted digit after P: plus the true digit after T:.

  • Run the temporary plotting cell.

The helper is now available for the designed obstacle. The next cell exposes how random weights behave before training.

Plotting helper not defined?

Confirm that you ran the whole helper cell. Check that every plotting line remains indented inside the loop where shown.

Rerun the Matplotlib import cell if plt is missing. Help me fix this prediction grid.

Before you run the final check, do you think random parameters can match all eight true digit labels?

  • Create one final temporary cell beneath the plotting helper.
  • Reveal the untrained predictions by running this code:
initial_outputs, _ = forward(X_test[:8], params)
initial_predictions = np.argmax(initial_outputs, axis=1)
print("Untrained MSE:", mean_squared_error(initial_outputs, y_test[:8]))
show_predictions(X_test, y_test, initial_predictions, count=8)

What should you see?

The output prints one untrained MSE value. You should also see a grid of eight handwritten digits.

Compare each P: value with its T: value. Several guesses should be incorrect because no parameter has learned from a label yet.

You have now exposed the intended failure: the complete forward pass runs, but its random guesses are unreliable. That baseline gives the next update something concrete to improve.

Prediction grid not appearing?

Rerun the plotting helper before rerunning the final check. Confirm that the earlier cells created X_test, y_test, and params.

Check that the final call uses count=8. Help me diagnose the missing prediction grid.

  • Delete the five temporary verification cells after capturing the grid.

Your runtime keeps the arrays plus parameters in memory. Your saved project code now contains only the required reusable functions.

✔️ Awesome, I've got everything!

Great. Double-check that the notebook has saved your reusable data, parameter, activation, forward-pass, and loss functions.

ⓧ I'd like to double check the full code

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml

SEED = 42
INPUT_SIZE = 784
HIDDEN_1 = 128
HIDDEN_2 = 64
OUTPUT_SIZE = 10
TRAIN_SIZE = 12_000
TEST_SIZE = 2_000
BATCH_SIZE = 128
EPOCHS = 8
LEARNING_RATE = 0.1


def load_data(rng):
    X, y = fetch_openml(
        "mnist_784",
        version=1,
        return_X_y=True,
        as_frame=False,
    )

    # Select reproducible training and test subsets
    order = rng.permutation(X.shape[0])
    selected = order[: TRAIN_SIZE + TEST_SIZE]
    X_selected = X[selected].astype(np.float32) / 255.0
    y_selected = y[selected].astype(np.int64)

    X_train = X_selected[:TRAIN_SIZE]
    y_train = y_selected[:TRAIN_SIZE]
    X_test = X_selected[TRAIN_SIZE:]
    y_test = y_selected[TRAIN_SIZE:]
    return X_train, y_train, X_test, y_test


def initialize_parameters(rng):
    # Scale each weight matrix by its input width for sigmoid activations
    return {
        "W1": rng.standard_normal((INPUT_SIZE, HIDDEN_1))
        / np.sqrt(INPUT_SIZE),
        "b1": np.zeros((1, HIDDEN_1)),
        "W2": rng.standard_normal((HIDDEN_1, HIDDEN_2))
        / np.sqrt(HIDDEN_1),
        "b2": np.zeros((1, HIDDEN_2)),
        "W3": rng.standard_normal((HIDDEN_2, OUTPUT_SIZE))
        / np.sqrt(HIDDEN_2),
        "b3": np.zeros((1, OUTPUT_SIZE)),
    }


def sigmoid(values):
    # Keep exponentiation stable for large pre-activation values
    clipped = np.clip(values, -500.0, 500.0)
    negative_exponentials = np.exp(-clipped)
    return 1.0 / (1.0 + negative_exponentials)


def forward(X, params):
    # Transform pixels through both sigmoid hidden layers
    Z1 = X @ params["W1"] + params["b1"]
    A1 = sigmoid(Z1)
    Z2 = A1 @ params["W2"] + params["b2"]
    A2 = sigmoid(Z2)

    # Produce one bounded output for each digit class
    Z3 = A2 @ params["W3"] + params["b3"]
    outputs = sigmoid(Z3)

    cache = {
        "X": X,
        "A1": A1,
        "A2": A2,
        "outputs": outputs,
    }
    return outputs, cache


def one_hot(y):
    # Create one target column for each digit class
    targets = np.zeros((y.shape[0], OUTPUT_SIZE))
    targets[np.arange(y.shape[0]), y] = 1.0
    return targets


def mean_squared_error(outputs, y):
    # Average the squared gap across every image and class
    targets = one_hot(y)
    squared_errors = (outputs - targets) ** 2
    return np.mean(squared_errors)

Your random network now turns MNIST pixels into inspectable sigmoid outputs. Next, you will carry its MSE gradient backward through both hidden layers so one update can lower the loss.

Backpropagate Through Both Hidden Layers

Your forward pass turned random parameters into sigmoid outputs. The eight digit guesses showed how unreliable those parameters are.

Backpropagation carries the output error through both hidden layers. The chain rule determines how much each parameter contributed to that error.

In this step, get ready to:
  • Form the output gradient from MSE and the output sigmoid.
  • Carry gradients backward through both sigmoid hidden layers.
  • Apply one parameter update that lowers the batch loss.
Form the output gradient

Mean squared error measures the gap between each sigmoid output and its one-hot target. The chain rule then carries that error through the output sigmoid before it reaches the final weight matrix.

  • Return to the code cell that ends with the mean_squared_error() function.
  • Place your cursor at the end of its final return line.
  • Add the accuracy calculation plus the MSE output gradient by pasting this code:


def accuracy(outputs, y):
    # Treat the largest output as the predicted digit
    predictions = np.argmax(outputs, axis=1)
    return np.mean(predictions == y)


def backward(y, params, cache):
    # Recover the activations used by the forward pass
    X = cache["X"]
    A1 = cache["A1"]
    A2 = cache["A2"]
    outputs = cache["outputs"]

    # Differentiate the mean squared error through the output sigmoid
    targets = one_hot(y)
    scale = 2.0 / (y.shape[0] * OUTPUT_SIZE)
    dA3 = scale * (outputs - targets)
    dZ3 = dA3 * outputs * (1.0 - outputs)

    # Measure how the output parameters affect the loss
    dW3 = A2.T @ dZ3
    db3 = np.sum(dZ3, axis=0, keepdims=True)

What does this first section do?

  • accuracy() selects the largest sigmoid output and compares its class index with each true label.
  • The cached activations recover the values needed for the reverse pass.
  • dA3 measures the MSE gradient with respect to every output activation.
  • dZ3 applies the output sigmoid derivative before dW3 and db3 measure how the final parameters affect the loss.
Carry gradients through both hidden layers

Each hidden sigmoid contributes the derivative activation * (1.0 - activation). The resulting gradients travel through W2 and W1 in reverse order.

  • Place your cursor at the end of the db3 line in the same code cell.
  • Complete the reverse pass plus the parameter updater by pasting this code:


    # Carry the gradient through the second sigmoid hidden layer
    dA2 = dZ3 @ params["W3"].T
    dZ2 = dA2 * A2 * (1.0 - A2)
    dW2 = A1.T @ dZ2
    db2 = np.sum(dZ2, axis=0, keepdims=True)

    # Carry the gradient through the first sigmoid hidden layer
    dA1 = dZ2 @ params["W2"].T
    dZ1 = dA1 * A1 * (1.0 - A1)
    dW1 = X.T @ dZ1
    db1 = np.sum(dZ1, axis=0, keepdims=True)

    return {
        "W1": dW1,
        "b1": db1,
        "W2": dW2,
        "b2": db2,
        "W3": dW3,
        "b3": db3,
    }


def update_parameters(params, gradients, learning_rate):
    # Move every parameter opposite its loss gradient
    for name in params:
        params[name] -= learning_rate * gradients[name]

How does the error travel backward?

  • dA2 carries the output-layer gradient through W3.
  • A2 * (1.0 - A2) applies the second hidden sigmoid derivative.
  • A1 * (1.0 - A1) applies the first hidden sigmoid derivative.
  • update_parameters() moves every parameter opposite its gradient.

✔️ Awesome, I've got everything!

Great. Double-check that the completed cell defines accuracy(), backward(), and update_parameters().

ⓧ I'd like to double check the full code

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml

SEED = 42
INPUT_SIZE = 784
HIDDEN_1 = 128
HIDDEN_2 = 64
OUTPUT_SIZE = 10
TRAIN_SIZE = 12_000
TEST_SIZE = 2_000
BATCH_SIZE = 128
EPOCHS = 8
LEARNING_RATE = 0.1


def load_data(rng):
    # Fetch the fixed MNIST release as NumPy arrays
    X, y = fetch_openml(
        "mnist_784",
        version=1,
        return_X_y=True,
        as_frame=False,
    )

    # Select reproducible training and test subsets
    order = rng.permutation(X.shape[0])
    selected = order[: TRAIN_SIZE + TEST_SIZE]
    X_selected = X[selected].astype(np.float32) / 255.0
    y_selected = y[selected].astype(np.int64)

    X_train = X_selected[:TRAIN_SIZE]
    y_train = y_selected[:TRAIN_SIZE]
    X_test = X_selected[TRAIN_SIZE:]
    y_test = y_selected[TRAIN_SIZE:]
    return X_train, y_train, X_test, y_test


def initialize_parameters(rng):
    # Scale each weight matrix by its input width for sigmoid activations
    return {
        "W1": rng.standard_normal((INPUT_SIZE, HIDDEN_1))
        / np.sqrt(INPUT_SIZE),
        "b1": np.zeros((1, HIDDEN_1)),
        "W2": rng.standard_normal((HIDDEN_1, HIDDEN_2))
        / np.sqrt(HIDDEN_1),
        "b2": np.zeros((1, HIDDEN_2)),
        "W3": rng.standard_normal((HIDDEN_2, OUTPUT_SIZE))
        / np.sqrt(HIDDEN_2),
        "b3": np.zeros((1, OUTPUT_SIZE)),
    }


def sigmoid(values):
    # Keep exponentiation stable for large pre-activation values
    clipped = np.clip(values, -500.0, 500.0)
    return 1.0 / (1.0 + np.exp(-clipped))


def forward(X, params):
    # Transform pixels through both sigmoid hidden layers
    Z1 = X @ params["W1"] + params["b1"]
    A1 = sigmoid(Z1)
    Z2 = A1 @ params["W2"] + params["b2"]
    A2 = sigmoid(Z2)

    # Produce one bounded output for each digit class
    Z3 = A2 @ params["W3"] + params["b3"]
    outputs = sigmoid(Z3)
    cache = {"X": X, "A1": A1, "A2": A2, "outputs": outputs}
    return outputs, cache


def one_hot(y):
    # Create one target column for each digit class
    targets = np.zeros((y.shape[0], OUTPUT_SIZE))
    targets[np.arange(y.shape[0]), y] = 1.0
    return targets


def mean_squared_error(outputs, y):
    # Average the squared gap across every image and class
    targets = one_hot(y)
    return np.mean((outputs - targets) ** 2)


def accuracy(outputs, y):
    # Treat the largest output as the predicted digit
    predictions = np.argmax(outputs, axis=1)
    return np.mean(predictions == y)


def backward(y, params, cache):
    # Recover the activations used by the forward pass
    X = cache["X"]
    A1 = cache["A1"]
    A2 = cache["A2"]
    outputs = cache["outputs"]

    # Differentiate MSE through the output sigmoid
    targets = one_hot(y)
    scale = 2.0 / (y.shape[0] * OUTPUT_SIZE)
    dA3 = scale * (outputs - targets)
    dZ3 = dA3 * outputs * (1.0 - outputs)
    dW3 = A2.T @ dZ3
    db3 = np.sum(dZ3, axis=0, keepdims=True)

    # Carry the gradient through both hidden sigmoids
    dA2 = dZ3 @ params["W3"].T
    dZ2 = dA2 * A2 * (1.0 - A2)
    dW2 = A1.T @ dZ2
    db2 = np.sum(dZ2, axis=0, keepdims=True)
    dA1 = dZ2 @ params["W2"].T
    dZ1 = dA1 * A1 * (1.0 - A1)
    dW1 = X.T @ dZ1
    db1 = np.sum(dZ1, axis=0, keepdims=True)

    return {
        "W1": dW1,
        "b1": db1,
        "W2": dW2,
        "b2": db2,
        "W3": dW3,
        "b3": db3,
    }


def update_parameters(params, gradients, learning_rate):
    # Move every parameter opposite its loss gradient
    for name in params:
        params[name] -= learning_rate * gradients[name]

How to use this reference

Compare this reference with your cumulative code cell. Each gradient name must match its corresponding parameter name so the updater can pair them correctly.

Apply one update and compare the loss

A gradient descent update should reduce the loss on the same batch used to calculate its gradients. This small test isolates the learning step before the full training loop adds shuffling and repeated epochs.

  • Run the completed definitions cell using its play button.

Before you run the diagnostic, predict whether the loss will move up or down.

  • Apply one update to 128 training images by running this code in a new cell:
# Reuse one batch before and after a single gradient update
X_batch = X_train[:BATCH_SIZE]
y_batch = y_train[:BATCH_SIZE]
before_outputs, before_cache = forward(X_batch, params)
before_loss = mean_squared_error(before_outputs, y_batch)
gradients = backward(y_batch, params, before_cache)
update_parameters(params, gradients, LEARNING_RATE)
after_outputs, _ = forward(X_batch, params)
after_loss = mean_squared_error(after_outputs, y_batch)
print("One-batch MSE before update:", before_loss)
print("One-batch MSE after one update:", after_loss)

What does this diagnostic prove?

  • The first forward pass measures the loss before any parameter changes.
  • The backward pass calculates gradients for all six parameter arrays.
  • The update changes each parameter using LEARNING_RATE.
  • The second forward pass measures whether that update improved the same batch.

You will see One-batch MSE before update followed by One-batch MSE after one update. The second value is lower than the first.

That is the learning mechanism working. One manual update has reduced the error across the same 128 images.

Is the second MSE higher?

  • Check that dA3 divides by the batch size multiplied by OUTPUT_SIZE.
  • Confirm that dZ3, dZ2, and dZ1 each include their matching sigmoid derivative.
  • Verify that every parameter update uses the gradient with the same dictionary key.

Use the two MSE values to narrow down the failing stage. Help me debug why my one-batch MSE did not decrease.

Your gradients now reach every weight and bias in the network. Next, you will repeat these updates across shuffled mini-batches to train the classifier.

Train, Evaluate, and Inspect Predictions

Your NumPy network now lowers its loss after one manual gradient update. That result proves your backpropagation can improve one batch.

Shuffled mini-batch gradient descent repeats that update across the training images for eight epochs. Held-out test images then reveal how well the learned parameters handle digits they did not train on.

In this step, get ready to:
  • Train the parameters with shuffled mini-batches for eight epochs.
  • Measure progress with epoch MSE summaries plus a final test accuracy.
  • Inspect ten predicted labels beside their true labels.
Train with shuffled mini-batches

A single batch covers only a small slice of the handwriting styles in MNIST. Shuffling the indices before each epoch gives every batch a different mix of images.

  • In the notebook code cell from earlier, place your cursor below the existing update_parameters() function.
  • Define the training loop at that location by pasting this code:
def train(X_train, y_train, params, rng):
    # Record one comparable MSE value after each epoch
    loss_history = []

    for epoch in range(EPOCHS):
        order = rng.permutation(X_train.shape[0])

        # Update the parameters with shuffled mini-batches
        for start in range(0, X_train.shape[0], BATCH_SIZE):
            batch_indices = order[start : start + BATCH_SIZE]
            X_batch = X_train[batch_indices]
            y_batch = y_train[batch_indices]
            outputs, cache = forward(X_batch, params)
            gradients = backward(y_batch, params, cache)
            update_parameters(params, gradients, LEARNING_RATE)

        # Measure progress on the same monitoring slice
        monitor_outputs, _ = forward(X_train[:2_000], params)
        epoch_loss = mean_squared_error(monitor_outputs, y_train[:2_000])
        epoch_accuracy = accuracy(monitor_outputs, y_train[:2_000])
        loss_history.append(epoch_loss)
        print(
            f"Epoch {epoch + 1}/{EPOCHS} | "
            f"mse={epoch_loss:.4f} | accuracy={epoch_accuracy:.3f}"
        )

    return loss_history

What does the training loop do?

  • The outer loop repeats training for the eight epochs stored in EPOCHS.
  • The shuffled order changes how training images are grouped during each epoch.
  • Each mini-batch passes through forward() before backward() calculates all six gradients.
  • The fixed monitoring slice produces one comparable MSE plus one accuracy value after each epoch.
  • Run the notebook code cell to register train() in the current Google Colab runtime.
  • Confirm the cell finishes without an exception.

Training function not registering?

Check that return loss_history aligns with the outer for loop. It must remain inside train().

Confirm every mini-batch line stays inside the inner loop. Help me check the indentation in my train function.

Build the reusable run

The finished notebook needs one reusable path from fresh parameters to visible predictions. A plotting helper keeps the image layout separate from the training logic.

  • Place your cursor below train() in the same notebook code cell.
  • Define the prediction-grid helper at that location by pasting this code:
def show_predictions(X, y, predictions, count):
    columns = count // 2
    _, axes = plt.subplots(2, columns)

    for index in range(count):
        ax = axes[index // columns, index % columns]
        ax.imshow(X[index].reshape(28, 28), cmap="gray")
        ax.set_title(f"P:{predictions[index]} T:{y[index]}")
        ax.axis("off")

    plt.show()

How does the prediction grid work?

  • The requested image count becomes two rows of equal-width columns.
  • Each flattened image is reshaped into its original 28 by 28 pixel layout.
  • Every title pairs a P: prediction with its T: true label.
  • The axes are hidden so the digit plus its labels remain the focus.
  • Run the notebook code cell to register show_predictions().
  • Confirm the cell completes without an exception.

Plotting helper showing an exception?

Check that plt.subplots(2, columns) remains inside show_predictions(). Confirm the loop body is indented beneath for index in range(count):.

Make sure matplotlib.pyplot is still imported as plt near the top of the notebook. Help me diagnose my prediction plotting helper.

The main() function rebuilds the seeded experiment from the beginning. Its first section also preserves the untrained prediction plus the one-batch diagnostic from your earlier checks.

  • Place your cursor below show_predictions() in the same notebook code cell.
  • Start the reusable experiment function at that location by pasting this code:
def main():
    # Rebuild the same data split and initial parameters
    rng = np.random.default_rng(SEED)
    X_train, y_train, X_test, y_test = load_data(rng)
    print("Train shape:", X_train.shape, "Test shape:", X_test.shape)
    params = initialize_parameters(rng)
    print("Parameter shapes:", {name: value.shape for name, value in params.items()})

    # Capture the untrained baseline
    initial_outputs, _ = forward(X_test[:8], params)
    initial_predictions = np.argmax(initial_outputs, axis=1)
    print("Untrained MSE:", mean_squared_error(initial_outputs, y_test[:8]))
    show_predictions(X_test, y_test, initial_predictions, count=8)

    # Confirm that one manual update lowers MSE
    X_batch = X_train[:BATCH_SIZE]
    y_batch = y_train[:BATCH_SIZE]
    before_outputs, before_cache = forward(X_batch, params)
    before_loss = mean_squared_error(before_outputs, y_batch)
    gradients = backward(y_batch, params, before_cache)
    update_parameters(params, gradients, LEARNING_RATE)
    after_outputs, _ = forward(X_batch, params)
    after_loss = mean_squared_error(after_outputs, y_batch)
    print("One-batch MSE before update:", before_loss)
    print("One-batch MSE after one update:", after_loss)

    # Train across all shuffled mini-batches
    loss_history = train(X_train, y_train, params, rng)

What does this part of main do?

  • The seeded random generator makes the data selection plus parameter initialization reproducible.
  • The untrained prediction grid records the unreliable random-weight baseline.
  • The one-batch diagnostic confirms that one manual update still lowers loss.
  • The call to train() applies repeated updates across shuffled mini-batches.
  • Run the notebook code cell to register the current main() definition.
  • Confirm the cell finishes without a syntax or indentation exception.

Main function not registering?

Check that every line after def main(): has one level of indentation. Keep the blank lines between logical sections.

Confirm loss_history = train(X_train, y_train, params, rng) remains inside the function. Help me find the syntax problem in my main function.

Evaluate the trained classifier

Training loss measures progress on images used during learning. Test accuracy measures whether those updates transfer to the held-out subset.

A Matplotlib loss curve makes the eight recorded checkpoints visible. The final digit grid exposes individual successes plus mistakes behind the overall score.

  • In main(), locate this final line:
    loss_history = train(X_train, y_train, params, rng)

Where does the evaluation belong?

This line finishes the parameter updates before evaluation begins. The trained params dictionary is ready for the held-out test subset at this point.

  • Add the evaluation plus visualization logic immediately below that line by pasting this code:
    # Evaluate the trained parameters on held-out images
    test_outputs, _ = forward(X_test, params)
    test_predictions = np.argmax(test_outputs, axis=1)
    test_accuracy = accuracy(test_outputs, y_test)
    print(f"Final test accuracy: {test_accuracy:.3f}")

    # Plot the recorded mean squared error after each epoch
    _, loss_ax = plt.subplots()
    loss_ax.plot(np.arange(1, EPOCHS + 1), loss_history)
    loss_ax.set_xlabel("Epoch")
    loss_ax.set_ylabel("Mean squared error")
    loss_ax.set_title("Training MSE")
    plt.show()

    # Display individual predictions behind the aggregate score
    show_predictions(X_test, y_test, test_predictions, count=10)
    return params, X_train, y_train


params, X_train, y_train = main()

What does the final evaluation do?

  • A final forward pass produces ten sigmoid outputs for every held-out test image.
  • The highest output in each row becomes the predicted digit.
  • The loss plot maps each recorded epoch to its monitored mean squared error.
  • The last grid displays ten test images with predicted and true labels.

Training on the CPU runtime may leave an epoch line on screen briefly while the next mini-batches run. That pause means the matrix updates are still being calculated.

Before you run the finished notebook, do you expect the trained prediction grid to contain more correct matches than the untrained grid?

  • Run the complete notebook code cell.

You will see eight epoch summaries followed by a final measured test accuracy. You will also see the training loss curve plus a ten-image grid with P: and T: labels.

That is the full classifier working from random initialization to held-out predictions. Your manual gradients now drive a complete training run.

Training run not completing?

If an epoch fails with a missing name, rerun the earlier definition cells so forward(), backward(), update_parameters() plus train() are available in the current runtime.

If the plots do not appear, confirm both plt.show() calls remain inside their functions. Help me diagnose my final MNIST training run.

✔️ Awesome, I've got everything!

Great work. Confirm the notebook has saved your complete classifier plus its outputs.

ⓧ I'd like to double check the full code

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml

SEED = 42
INPUT_SIZE = 784
HIDDEN_1 = 128
HIDDEN_2 = 64
OUTPUT_SIZE = 10
TRAIN_SIZE = 12_000
TEST_SIZE = 2_000
BATCH_SIZE = 128
EPOCHS = 8
LEARNING_RATE = 0.1


def load_data(rng):
    # Fetch the fixed MNIST release as NumPy arrays
    X, y = fetch_openml(
        "mnist_784", version=1, return_X_y=True, as_frame=False
    )
    order = rng.permutation(X.shape[0])
    selected = order[: TRAIN_SIZE + TEST_SIZE]
    X_selected = X[selected].astype(np.float32) / 255.0
    y_selected = y[selected].astype(np.int64)
    X_train = X_selected[:TRAIN_SIZE]
    y_train = y_selected[:TRAIN_SIZE]
    X_test = X_selected[TRAIN_SIZE:]
    y_test = y_selected[TRAIN_SIZE:]
    return X_train, y_train, X_test, y_test


def initialize_parameters(rng):
    # Scale weights by their input widths for sigmoid activations
    return {
        "W1": rng.standard_normal((INPUT_SIZE, HIDDEN_1)) / np.sqrt(INPUT_SIZE),
        "b1": np.zeros((1, HIDDEN_1)),
        "W2": rng.standard_normal((HIDDEN_1, HIDDEN_2)) / np.sqrt(HIDDEN_1),
        "b2": np.zeros((1, HIDDEN_2)),
        "W3": rng.standard_normal((HIDDEN_2, OUTPUT_SIZE)) / np.sqrt(HIDDEN_2),
        "b3": np.zeros((1, OUTPUT_SIZE)),
    }


def sigmoid(values):
    # Keep exponentiation stable for large pre-activation values
    clipped = np.clip(values, -500.0, 500.0)
    return 1.0 / (1.0 + np.exp(-clipped))


def forward(X, params):
    # Transform pixels through three sigmoid layers
    A1 = sigmoid(X @ params["W1"] + params["b1"])
    A2 = sigmoid(A1 @ params["W2"] + params["b2"])
    outputs = sigmoid(A2 @ params["W3"] + params["b3"])
    cache = {"X": X, "A1": A1, "A2": A2, "outputs": outputs}
    return outputs, cache


def one_hot(y):
    # Create one target column for each digit class
    targets = np.zeros((y.shape[0], OUTPUT_SIZE))
    targets[np.arange(y.shape[0]), y] = 1.0
    return targets


def mean_squared_error(outputs, y):
    # Average the squared gap across every image and class
    return np.mean((outputs - one_hot(y)) ** 2)


def accuracy(outputs, y):
    # Treat the largest output as the predicted digit
    return np.mean(np.argmax(outputs, axis=1) == y)


def backward(y, params, cache):
    # Differentiate MSE through all three sigmoid layers
    X = cache["X"]
    A1 = cache["A1"]
    A2 = cache["A2"]
    outputs = cache["outputs"]
    dA3 = 2.0 * (outputs - one_hot(y)) / (y.shape[0] * OUTPUT_SIZE)
    dZ3 = dA3 * outputs * (1.0 - outputs)
    dW3 = A2.T @ dZ3
    db3 = np.sum(dZ3, axis=0, keepdims=True)
    dA2 = dZ3 @ params["W3"].T
    dZ2 = dA2 * A2 * (1.0 - A2)
    dW2 = A1.T @ dZ2
    db2 = np.sum(dZ2, axis=0, keepdims=True)
    dA1 = dZ2 @ params["W2"].T
    dZ1 = dA1 * A1 * (1.0 - A1)
    dW1 = X.T @ dZ1
    db1 = np.sum(dZ1, axis=0, keepdims=True)
    return {"W1": dW1, "b1": db1, "W2": dW2, "b2": db2, "W3": dW3, "b3": db3}


def update_parameters(params, gradients, learning_rate):
    # Move every parameter opposite its loss gradient
    for name in params:
        params[name] -= learning_rate * gradients[name]


def train(X_train, y_train, params, rng):
    # Record one comparable MSE value after each epoch
    loss_history = []
    for epoch in range(EPOCHS):
        order = rng.permutation(X_train.shape[0])
        for start in range(0, X_train.shape[0], BATCH_SIZE):
            batch_indices = order[start : start + BATCH_SIZE]
            X_batch = X_train[batch_indices]
            y_batch = y_train[batch_indices]
            outputs, cache = forward(X_batch, params)
            gradients = backward(y_batch, params, cache)
            update_parameters(params, gradients, LEARNING_RATE)

        # Monitor the same training slice after each epoch
        monitor_outputs, _ = forward(X_train[:2_000], params)
        epoch_loss = mean_squared_error(monitor_outputs, y_train[:2_000])
        epoch_accuracy = accuracy(monitor_outputs, y_train[:2_000])
        loss_history.append(epoch_loss)
        print(
            f"Epoch {epoch + 1}/{EPOCHS} | "
            f"mse={epoch_loss:.4f} | accuracy={epoch_accuracy:.3f}"
        )
    return loss_history


def show_predictions(X, y, predictions, count):
    # Display predictions beside their true digit labels
    columns = count // 2
    _, axes = plt.subplots(2, columns)
    for index in range(count):
        ax = axes[index // columns, index % columns]
        ax.imshow(X[index].reshape(28, 28), cmap="gray")
        ax.set_title(f"P:{predictions[index]} T:{y[index]}")
        ax.axis("off")
    plt.show()
def main():
    # Rebuild the seeded experiment
    rng = np.random.default_rng(SEED)
    X_train, y_train, X_test, y_test = load_data(rng)
    print("Train shape:", X_train.shape, "Test shape:", X_test.shape)
    params = initialize_parameters(rng)
    print("Parameter shapes:", {name: value.shape for name, value in params.items()})

    # Capture the untrained prediction baseline
    initial_outputs, _ = forward(X_test[:8], params)
    initial_predictions = np.argmax(initial_outputs, axis=1)
    print("Untrained MSE:", mean_squared_error(initial_outputs, y_test[:8]))
    show_predictions(X_test, y_test, initial_predictions, count=8)

    # Confirm that one manual update lowers MSE
    X_batch = X_train[:BATCH_SIZE]
    y_batch = y_train[:BATCH_SIZE]
    before_outputs, before_cache = forward(X_batch, params)
    before_loss = mean_squared_error(before_outputs, y_batch)
    gradients = backward(y_batch, params, before_cache)
    update_parameters(params, gradients, LEARNING_RATE)
    after_outputs, _ = forward(X_batch, params)
    after_loss = mean_squared_error(after_outputs, y_batch)
    print("One-batch MSE before update:", before_loss)
    print("One-batch MSE after one update:", after_loss)

    # Train and evaluate the classifier
    loss_history = train(X_train, y_train, params, rng)
    test_outputs, _ = forward(X_test, params)
    test_predictions = np.argmax(test_outputs, axis=1)
    test_accuracy = accuracy(test_outputs, y_test)
    print(f"Final test accuracy: {test_accuracy:.3f}")

    # Plot the recorded MSE and show individual predictions
    _, loss_ax = plt.subplots()
    loss_ax.plot(np.arange(1, EPOCHS + 1), loss_history)
    loss_ax.set_xlabel("Epoch")
    loss_ax.set_ylabel("Mean squared error")
    loss_ax.set_title("Training MSE")
    plt.show()
    show_predictions(X_test, y_test, test_predictions, count=10)
    return params, X_train, y_train


params, X_train, y_train = main()

Secret mission

Numerically Check a Backpropagation Gradient

Use finite differences to test one analytical W3 gradient independently. You will compare the numerical slope with the value produced by manual backpropagation.

Clean Up Your Resources

Clean Up Your Resources

Your Google Colab project uses no paid resources. Decide whether to keep the project, disconnect its runtime for later, or delete its saved notebook.

Resources you used:

  • Your Google Drive notebook containing the completed NumPy classifier.
  • The temporary Google Colab hosted runtime containing your trained parameters.
  • The temporary OpenML cache containing the MNIST dataset.

Keep everything running

No cleanup is needed while you are still working. Colab can replace the hosted runtime after its temporary session expires.

  • Leave the project notebook in Google Drive so its code and saved outputs remain available.
  • Keep the environment setup cells at the top of the notebook so they can reconstruct the runtime.
  • Rerun the notebook from the first cell when Colab provides a new runtime.

Pause - I'll come back to this later

Disconnecting the runtime pauses your active session. The notebook remains saved in Google Drive.

  • Switch back to the Google Colab notebook from earlier.
  • Disconnect the current runtime using Colab's runtime connection controls.

Your code and saved outputs remain in the notebook. A future session reconstructs the temporary data when you rerun the cells.

Delete - I don't want to use this again

Full removal clears the temporary runtime data. It also permanently removes the notebook containing your code and saved outputs.

  • Switch back to the Google Colab notebook from earlier.
  • Select Runtime in the top menu.
  • Select Disconnect and delete runtime.
  • Complete the confirmation dialog.

You have cleared the temporary environment. Its cached MNIST data and trained parameters are gone.

  • Open Google Drive in a new browser tab.
  • Locate the Colab notebook used for this project.
  • Delete the project notebook from Google Drive.
  • Remove the project notebook permanently from Google Drive's trash.

Cleanup complete. You have removed the saved notebook and its temporary runtime data.

Nice Work!

Nice Work!

You did it! Your Google Colab notebook now proves that a two-hidden-layer MNIST classifier can learn through gradients you implemented with NumPy.

You've learned how to:

  • Build a 784 -> 128 -> 64 -> 10 classifier that turns pixel arrays into ten bounded outputs through two sigmoid hidden layers.
  • Implement manual backpropagation through the output layer and both hidden layers. Inspect the gradients for W1, b1, W2, b2, W3, and b3.
  • Train with shuffled mini-batches and track a falling mean squared error curve. Measure accuracy on held-out images and inspect individual predictions.
  • Secret Mission: Use finite differences to compare an analytical W3 gradient with an independently estimated numerical gradient.

Ready to quiz yourself?