Lab: Neural Networks from Scratch (MNIST)

This is the "Moment of Truth." After exploring forward propagation, backpropagation, and optimization, it's time to assemble everything into a fully functional Deep Neural Network—built entirely from scratch using only NumPy. We will apply it to the MNIST dataset, the "Hello World" of computer vision.

---

The Objective

Build a multi-layer perceptron (MLP) to classify 70,000 images of handwritten digits (0-9). Each image is $28 \times 28$ pixels, which we flatten into a 784-dimensional vector.

Step 1: Data Preparation

We use fetch_openml to load the dataset and perform standard preprocessing: normalization (0-1 range) and One-Hot Encoding for the labels.

train_mnist.py

import numpy as np
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import OneHotEncoder

def load_and_prepare_data():
    print("Loading MNIST dataset...")
    X, y = fetch_openml('mnist_784', version=1, return_X_y=True, as_frame=False, parser='auto')
    
    # 1. Normalize pixel values (0-255 -> 0.0-1.0)
    X = X / 255.0
    
    # 2. One-hot encode the labels
    encoder = OneHotEncoder(sparse_output=False)
    Y_onehot = encoder.fit_transform(y.reshape(-1, 1))

    # 3. Split the data (90% Train, 10% Test)
    X_train, X_test, Y_train, Y_test = train_test_split(X, Y_onehot, test_size=0.1, random_state=42)
    
    # Transpose to match our math: (features, samples)
    return X_train.T, X_test.T, Y_train.T, Y_test.T, y.astype(int)

---

Step 2: The Core Neural Network Class

Here is our engine. It implements Vectorized Forward Pass, Backpropagation, and Mini-batch Gradient Descent.

nn.py

import numpy as np

class NeuralNetwork:
    def __init__(self, layer_sizes):
        self.L = len(layer_sizes) - 1
        # Xavier/Glorot Initialization
        self.weights = [np.random.randn(layer_sizes[l], layer_sizes[l-1]) * np.sqrt(1/layer_sizes[l-1]) 
                        for l in range(1, self.L + 1)]
        self.biases = [np.zeros((layer_sizes[l], 1)) for l in range(1, self.L + 1)]

    def sigmoid(self, z): return 1 / (1 + np.exp(-z))
    
    def softmax(self, z):
        exp_z = np.exp(z - np.max(z, axis=0, keepdims=True))
        return exp_z / np.sum(exp_z, axis=0, keepdims=True)

    def forward(self, X):
        cache = {"A0": X}
        for l in range(1, self.L):
            Z = np.dot(self.weights[l-1], cache[f"A{l-1}"]) + self.biases[l-1]
            cache[f"A{l}"] = self.sigmoid(Z)
        
        # Output layer with Softmax
        ZL = np.dot(self.weights[self.L-1], cache[f"A{self.L-1}"]) + self.biases[self.L-1]
        cache[f"A{self.L}"] = self.softmax(ZL)
        return cache[f"A{self.L}"], cache

    def backward(self, AL, Y, cache):
        m = Y.shape[1]
        grads = {}
        
        # Initial dZ for Softmax + Cross-Entropy
        dZ = AL - Y 
        
        for l in range(self.L, 0, -1):
            grads[f"dW{l}"] = (1/m) * np.dot(dZ, cache[f"A{l-1}"].T)
            grads[f"db{l}"] = (1/m) * np.sum(dZ, axis=1, keepdims=True)
            if l > 1:
                dA_prev = np.dot(self.weights[l-1].T, dZ)
                # Sigmoid derivative: A * (1 - A)
                dZ = dA_prev * cache[f"A{l-1}"] * (1 - cache[f"A{l-1}"])
        return grads

    def update(self, grads, lr):
        for l in range(1, self.L + 1):
            self.weights[l-1] -= lr * grads[f"dW{l}"]
            self.biases[l-1] -= lr * grads[f"db{l}"]

    def train(self, X, Y, epochs=50, lr=0.1, batch_size=64):
        m = X.shape[1]
        for i in range(epochs):
            permutation = np.random.permutation(m)
            X_shuffled, Y_shuffled = X[:, permutation], Y[:, permutation]
            
            for j in range(0, m, batch_size):
                X_batch = X_shuffled[:, j:j+batch_size]
                Y_batch = Y_shuffled[:, j:j+batch_size]
                
                AL, cache = self.forward(X_batch)
                grads = self.backward(AL, Y_batch, cache)
                self.update(grads, lr)
            
            if i % 5 == 0:
                AL_full, _ = self.forward(X)
                cost = -np.mean(np.sum(Y * np.log(AL_full + 1e-15), axis=0))
                print(f"Epoch {i}: Cost {cost:.4f}")

---

Step 3: Execution and Evaluation

main.py

X_train, X_test, Y_train, Y_test, y_orig = load_and_prepare_data()

# Architecture: 784 -> 128 -> 64 -> 10
model = NeuralNetwork([784, 128, 64, 10])
model.train(X_train, Y_train, epochs=30, lr=0.1)

# Prediction
AL_test, _ = model.forward(X_test)
predictions = np.argmax(AL_test, axis=0)
accuracy = np.mean(predictions == np.argmax(Y_test, axis=0))
print(f"Test Accuracy: {accuracy * 100:.2f}%")

---

Visualizing Results

A professional lab doesn't just print numbers; it visualizes the network's behavior.

1. Training Curve

If you plot the costs stored during training, you should see a smooth exponential decay, indicating successful convergence.

2. Sample Predictions

Let's see the model in action:

import matplotlib.pyplot as plt

def visualize_prediction(index):
    img = X_test[:, index].reshape(28, 28)
    plt.imshow(img, cmap='gray')
    pred = model.predict(X_test[:, index:index+1])
    plt.title(f"Label: {np.argmax(Y_test[:, index])}, Prediction: {pred[0]}")
    plt.show()

Pro Tip: Why Xavier Initialization?

Notice the np.sqrt(1/layer_sizes[l-1]) in the weight initialization. This is Xavier Initialization. It prevents gradients from exploding or vanishing in the first few iterations by keeping the variance of activations constant across layers.

---

Notebook Implementation

For an interactive, step-by-step experience with plots, check the MNIST from Scratch Google Colab.