Building MINST Neural Network only use Numpy

Building a neural network for MNIST digits from scratch using only NumPy requires pulling back the abstractions of modern frameworks. Below is the end-to-end mathematical breakdown—from raw grayscale pixels to calculus-driven gradient descent—before writing a single line of code.

Core Conceptual Pillars

Classifying digits ($0$ through $9$) with a feedforward architecture relies on three primary pillars:

• Architecture & Dimensions: Flattening a $28 \times 28$ grayscale image into vectors and defining explicit matrix dimensions for weights ($W$) and biases ($b$).

• Forward Propagation & Activations: Executing linear transformations ($Z = WX + b$) and non-linearities ($\text{ReLU}$ for hidden representations, $\text{Softmax}$ for multi-class probabilities).

• Loss & Backpropagation: Quantifying prediction divergence using Categorical Cross-Entropy, applying the multivariable chain rule, and updating parameters via Gradient Descent.

1. Network Architecture & Dimensions

A standard two-layer architecture (one hidden layer, one output layer) handles the MNIST classification task:

Layer	Input Dimension	Output Dimension	Weight Matrix (W)	Bias Vector (b)	Activation
Input ($X$)	—	$784 \times 1$	—	—	Flatten
Hidden Layer ($[1]$)	$784 \times 1$	$128 \times 1$	$W^{[1]} \in \mathbb{R}^{128 \times 784}$	$b^{[1]} \in \mathbb{R}^{128 \times 1}$	$\text{ReLU}$
Output Layer ($[2]$)	$128 \times 1$	$10 \times 1$	$W^{[2]} \in \mathbb{R}^{10 \times 128}$	$b^{[2]} \in \mathbb{R}^{10 \times 1}$	$\text{Softmax}$

The Input Layer

Each MNIST sample is a single-channel grayscale image of $28 \times 28$ pixels. Fully connected networks require a 1D column vector:

$28 \times 28 = 784 \implies X \in \mathbb{R}^{784 \times 1}$

The Hidden Layer

With $128$ hidden units, the pre-activation equation is:

$Z^{[1]} = W^{[1]} X + b^{[1]}$

To satisfy matrix dimension constraints $(m \times k) \times (k \times n) = (m \times n)$:

$(128 \times 784) \times (784 \times 1) + (128 \times 1) = (128 \times 1)$

Rather than viewing $W^{[1]}X$ as abstract arithmetic, interpret it as 128 parallel template matchers:

1. Row-Wise Decomposition: Each row $i$ of $W^{[1]}$ contains $784$ weights corresponding to neuron $i$:

   $z_i = \sum_{k=1}^{784} W_{i, k} x_k + b_i$

2. Neuron as a Spatial Filter: Reshaping that $1 \times 784$ row back into a $28 \times 28$ grid reveals a spatial blueprint. Computing the dot product acts as a cross-correlation:

   • Bright pixel $\times$ Positive weight: High positive activation.

   • Bright pixel $\times$ Negative weight: Penalty/inhibition.

   • Zero/Background pixel: No contribution.

3. Feature Projection ($784 \to 128$): The network shifts from raw sensory questions (“Is pixel (14, 12) lit?”) to semantic building blocks (“Does the image contain an upper loop or a vertical stem?“).

2. Forward Propagation & Activations

Hidden Activation: ReLU

The hidden layer maps linear combinations through a Rectified Linear Unit ($\text{ReLU}$) element-wise:

$A^{[1]} = \max(0, Z^{[1]}), \quad A^{[1]} \in \mathbb{R}^{128 \times 1}$

Output Transformation & Softmax

The output layer maps the $128$ features to $10$ class logits:

$Z^{[2]} = W^{[2]} A^{[1]} + b^{[2]}, \quad Z^{[2]} \in \mathbb{R}^{10 \times 1}$

To turn raw real-valued scores into a normalized probability distribution where $\sum_{j=0}^{9} \hat{y}_j = 1$, apply the Softmax function:

$\hat{y}_i = A^{[2]}_i = \frac{e^{Z^{[2]}_i}}{\sum_{j=0}^{9} e^{Z^{[2]}_j}}$

3. Target Encoding & Loss Function

One-Hot Encoding

Targets $y \in \{0, \dots, 9\}$ are converted to 10-dimensional unit vectors. For target digit $3$:

$y = \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}^T$

Categorical Cross-Entropy

Loss for a single sample measures the divergence between target $y$ and prediction $\hat{y}$:

$L(y, \hat{y}) = -\sum_{i=0}^{9} y_i \ln(\hat{y}_i)$

Because $y_i = 0$ for all classes except target $k$, the expression collapses to the negative log-likelihood of the correct class:

$L = -\ln(\hat{y}_k)$

4. Backpropagation & The Output Layer Gradients

Gradients trace how changes to intermediate variables affect total loss $L$ using the multivariable chain rule.

Output Error ($dZ^{[2]}$)

The derivative of the loss with respect to raw output logits $Z^{[2]}$ is defined as:

$dZ^{[2]} = \frac{\partial L}{\partial Z^{[2]}}$

Because changing a single logit $Z_i$ alters the normalization denominator $\sum_k e^{Z_k}$, every predicted probability $\hat{y}_j$ shifts. By the multivariable chain rule:

$\frac{\partial L}{\partial Z_i^{[2]}} = \sum_{j=0}^{9} \frac{\partial L}{\partial \hat{y}_j} \frac{\partial \hat{y}_j}{\partial Z_i^{[2]}}$

Evaluating this summation resolves with clean algebraic cancellation:

$dZ^{[2]} = \hat{y} - y \in \mathbb{R}^{10 \times 1}$

The gradient at the output is simply the error vector—predicted probabilities minus ground truth.

Output Parameters ($dW^{[2]}, db^{[2]}$)

Applying the chain rule through the linear transformation $Z^{[2]} = W^{[2]} A^{[1]} + b^{[2]}$:

• Weight Gradient ($dW^{[2]}$): $dW^{[2]} = dZ^{[2]} (A^{[1]})^T \implies (10 \times 1) \times (1 \times 128) = 10 \times 128$

• Bias Gradient ($db^{[2]}$): $db^{[2]} = dZ^{[2]} \implies 10 \times 1$

Examining a single weight scalar $W_{jk}$ connecting hidden unit $k$ to output unit $j$:

$\frac{\partial L}{\partial W_{jk}} = \frac{\partial L}{\partial Z_j} \cdot \frac{\partial Z_j}{\partial W_{jk}} = dz_j \cdot a_k$

$\text{Weight Gradient} = (\text{error at destination neuron}) \times (\text{incoming activation})$

To compute all pairwise products between error vector $dZ \in \mathbb{R}^{m \times 1}$ and activations $A \in \mathbb{R}^{n \times 1}$ simultaneously, matrix multiplication requires transposing $A$:

$dZ \cdot A^T = \begin{bmatrix} dz_1 \\ \vdots \\ dz_m \end{bmatrix} \begin{bmatrix} a_1 & \dots & a_n \end{bmatrix} = \begin{bmatrix} dz_1 a_1 & \dots & dz_1 a_n \\ \vdots & \ddots & \vdots \\ dz_m a_1 & \dots & dz_m a_n \end{bmatrix}$

Summary of Parameter Derivatives

Parameter	Scalar Form	Matrix Form	Physical Intuition
Weight ($W_{jk}$)	$dz_j \cdot a_k$	$dZ \cdot A^T$	Scales the downstream error by how actively the input unit was firing.
Bias ($b_j$)	$dz_j \cdot 1$	$dZ$	Adds directly to pre-activation ($z_j$), absorbing downstream error directly.

5. Backpropagation: The Hidden Layer Gradients

With the output layer gradients ($dW^{[2]}$, $db^{[2]}$) resolved, we must propagate the error backward across the hidden layer to update $W^{[1]}$ and $b^{[1]}$.

---

Step 1: Pulling Error Through the Weights ($dA^{[1]}$)

We need to determine how the total loss changes with respect to the post-activation hidden states:

$dA^{[1]} = \frac{\partial L}{\partial A^{[1]}}$

Recall the forward relationship: $Z^{[2]} = W^{[2]} A^{[1]} + b^{[2]}$. Because hidden unit $k$ connects to every output neuron $j \in \{0, \dots, 9\}$, any change in $a_k^{[1]}$ propagates downstream along 10 distinct pathways.

Applying the multivariable chain rule:

$\frac{\partial L}{\partial a_k^{[1]}} = \sum_{j=0}^{9} \frac{\partial L}{\partial z_j^{[2]}} \cdot \frac{\partial z_j^{[2]}}{\partial a_k^{[1]}} = \sum_{j=0}^{9} dz_j^{[2]} \cdot W_{jk}^{[2]}$

In matrix notation, this summation corresponds to multiplying the incoming error vector $dZ^{[2]}$ by the transpose of the weight matrix $(W^{[2]})^T$:

$dA^{[1]} = (W^{[2]})^T dZ^{[2]}$

• Matrix Dimensions: $(128 \times 10) \times (10 \times 1) = 128 \times 1$.

Intuition: The Transpose as a Reverse Router

During forward propagation, $W^{[2]}$ routes signals forward from $128$ hidden units to $10$ output units:

$Z^{[2]} = W^{[2]} A^{[1]} \quad [(10 \times 128) \times (128 \times 1)]$

During backward propagation, we run the network in reverse. Transposing the matrix to $(W^{[2]})^T$ inverts the mapping:

$(W^{[2]})^T dZ^{[2]} \quad [(128 \times 10) \times (10 \times 1)]$

Each hidden unit $k$ receives a fraction of the error from each output neuron $j$, directly proportional to the strength of the connection $W_{jk}^{[2]}$ that linked them during the forward pass.

---

Step 2: Passing Through the ReLU Gate ($dZ^{[1]}$)

Before we can calculate parameter updates, the error must pass backward through the activation function:

$A^{[1]} = \text{ReLU}(Z^{[1]}) = \max(0, Z^{[1]})$

Because ReLU is an element-wise operation, each neuron’s pre-activation $z_k^{[1]}$ affects only its own post-activation $a_k^{[1]}$:

$dz_k^{[1]} = \frac{\partial L}{\partial z_k^{[1]}} = \frac{\partial L}{\partial a_k^{[1]}} \cdot \frac{\partial a_k^{[1]}}{\partial z_k^{[1]}} = da_k^{[1]} \cdot \text{ReLU}'(z_k^{[1]})$

The derivative of ReLU is a simple binary indicator:

$\text{ReLU}'(z) = \begin{cases} 1 & \text{if } z > 0 \\ 0 & \text{if } z \le 0 \end{cases}$

Using the Hadamard (element-wise) product $\odot$:

$dZ^{[1]} = dA^{[1]} \odot \mathbf{1}_{Z^{[1]} > 0}$

• Matrix Dimensions: $(128 \times 1) \odot (128 \times 1) = 128 \times 1$.
• The “Dead Gate” Mechanism: If a neuron had a negative pre-activation during forward propagation ($z \le 0$), its derivative is $0$. It completely zeros out the gradient, stopping backpropagation dead in its tracks for that neuron.

---

Step 3: Gradients for Layer 1 Parameters ($dW^{[1]}, db^{[1]}$)

Finally, we project the filtered error $dZ^{[1]}$ back to the network’s input features $X$:

$Z^{[1]} = W^{[1]} X + b^{[1]}$

Using the identical outer-product relationship established for Layer 2:

• Weight Gradient ($dW^{[1]}$): $dW^{[1]} = dZ^{[1]} X^T$ $\text{Dimensions: } (128 \times 1) \times (1 \times 784) = 128 \times 784$

• Bias Gradient ($db^{[1]}$): $db^{[1]} = dZ^{[1]}$ $\text{Dimensions: } 128 \times 1$

---

Complete Gradient Computational Graph

Tracing backwards from loss to raw pixels:

Loss L
  │
  ▼
dZ[2] = ŷ - y                    (10 × 1)
  ├──► dW[2] = dZ[2] · (A[1])ᵀ    (10 × 128)
  ├──► db[2] = dZ[2]             (10 × 1)
  │
  ▼
dA[1] = (W[2])ᵀ · dZ[2]          (128 × 1)
  │
  ▼
dZ[1] = dA[1] ⊙ (Z[1] > 0)       (128 × 1)
  ├──► dW[1] = dZ[1] · Xᵀ        (128 × 784)
  └──► db[1] = dZ[1]             (128 × 1)

6. Bridging Math to Code: Vectorization, Initialization, and Numerical Stability

Before translating our derivations into NumPy, we must address three practical engineering realities: processing multiple samples in parallel, initializing weights to preserve gradient scale, and preventing floating-point overflow.

---

1. Batch Vectorization ($m$ Samples)

Training one image at a time via pure Python loops is computationally prohibitive. Instead, we stack $m$ training examples side-by-side into a single batch matrix:

$X = \begin{bmatrix} \vert & \vert & & \vert \\ x^{(1)} & x^{(2)} & \cdots & x^{(m)} \\ \vert & \vert & & \vert \end{bmatrix} \in \mathbb{R}^{784 \times m}$

Because linear algebra distributes over columns, the forward propagation equations remain identical in matrix notation:

$Z^{[1]} = W^{[1]} X + b^{[1]}, \quad Z^{[1]} \in \mathbb{R}^{128 \times m}$ $Z^{[2]} = W^{[2]} A^{[1]} + b^{[2]}, \quad Z^{[2]} \in \mathbb{R}^{10 \times m}$

Broadcasting Note: The bias vector $b^{[l]} \in \mathbb{R}^{n^{[l]} \times 1}$ is automatically copied across all $m$ columns by NumPy during addition.

When computing gradients over a batch, the loss is defined as the mean error across all $m$ examples. Thus, parameter gradients must be averaged by dividing by $m$:

Gradient	Vectorized Formula	NumPy Equivalent
Output Error	$dZ^{[2]} = \hat{Y} - Y$	dZ2 = A2 - Y
Layer 2 Weights	$dW^{[2]} = \frac{1}{m} dZ^{[2]} (A^{[1]})^T$	dW2 = (1 / m) * np.dot(dZ2, A1.T)
Layer 2 Bias	$db^{[2]} = \frac{1}{m} \sum_{i=1}^m dZ^{[2](i)}$	db2 = (1 / m) * np.sum(dZ2, axis=1, keepdims=True)
Hidden Error	$dZ^{[1]} = \left((W^{[2]})^T dZ^{[2]}\right) \odot \mathbf{1}_{Z^{[1]} > 0}$	dZ1 = np.dot(W2.T, dZ2) * (Z1 > 0)
Layer 1 Weights	$dW^{[1]} = \frac{1}{m} dZ^{[1]} X^T$	dW1 = (1 / m) * np.dot(dZ1, X.T)
Layer 1 Bias	$db^{[1]} = \frac{1}{m} \sum_{i=1}^m dZ^{[1](i)}$	db1 = (1 / m) * np.sum(dZ1, axis=1, keepdims=True)

keepdims=True ensures the bias shape stays $(n \times 1)$ rather than collapsing into a 1D rank-1 array $(n,)$, which prevents subtle broadcasting bugs.

---

2. Parameter Initialization (He / Kaiming Normal)

How weights are initialized dictates whether signals survive propagation across deep layers.

• The Zero Trap ($W = 0$): If weights are initialized to zeros, every neuron in hidden layer 1 computes $z = 0 + b$, outputs identical activations, and receives identical gradients. The neurons remain symmetric and cannot learn distinct features.
• Large Random Values: Initializing with unscaled random normals (e.g., standard normal $\mathcal{N}(0, 1)$) causes the variance of layer outputs to grow with the number of input connections ($n_{\text{in}}$), triggering exploding activations.

Because Layer 1 uses ReLU, which zeroes out negative inputs and halves output variance, we use He (Kaiming) Initialization. Weights are sampled from a zero-mean normal distribution scaled by $\sqrt{2 / n_{\text{in}}}$:

$W^{[l]} \sim \mathcal{N}\left(0, \sqrt{\frac{2}{n_{\text{in}}}}\right)$

# Layer 1: n_in = 784, n_out = 128
W1 = np.random.randn(128, 784) * np.sqrt(2.0 / 784)
b1 = np.zeros((128, 1))

# Layer 2: n_in = 128, n_out = 10
W2 = np.random.randn(10, 128) * np.sqrt(2.0 / 128)
b2 = np.zeros((10, 1))

Biases contribute no symmetry issues and are safely initialized to zeros.

3. Numerically Stable Softmax

The standard Softmax formula involves computing exponentials directly:

$\text{Softmax}(Z)_{ij} = \frac{e^{Z_{ij}}}{\sum_{k} e^{Z_{kj}}}$

In 64-bit floating point arithmetic (float64), $e^{x}$ overflows to inf when $x \gtrapprox 709.78$. In 32-bit (float32), this occurs around $x \gtrapprox 88.72$. Once an exponential overflows, division produces NaN.

To resolve this, apply the identity:

$\frac{e^{Z_{ij}}}{\sum_k e^{Z_{kj}}} = \frac{e^{-C} \cdot e^{Z_{ij}}}{e^{-C} \cdot \sum_k e^{Z_{kj}}} = \frac{e^{Z_{ij} - C}}{\sum_k e^{Z_{kj} - C}}$

By choosing $C = \max_k(Z_{kj})$ (the maximum logit for each column sample):

• The largest exponent becomes $e^{\max(Z) - \max(Z)} = e^0 = 1$.

• All other terms are strictly $\le 0$, bounding exponents to the interval $(0, 1]$ and eliminating overflow risks entirely.

Python

def softmax(Z):
    # Shift Z along class axis (axis=0) by the maximum logit of each sample
    shifted_Z = Z - np.max(Z, axis=0, keepdims=True)
    exp_Z = np.exp(shifted_Z)
    return exp_Z / np.sum(exp_Z, axis=0, keepdims=True)

4. Parameter Updates (SGD)

With the gradients averaged across the batch, update parameters using learning rate $\alpha$:

$W^{[1]} := W^{[1]} - \alpha \, dW^{[1]}$

$b^{[1]} := b^{[1]} - \alpha \, db^{[1]}$

$W^{[2]} := W^{[2]} - \alpha \, dW^{[2]}$

$b^{[2]} := b^{[2]} - \alpha \, db^{[2]}$

<FollowUp label="Want to assemble the complete, runnable NumPy training script?" query="Write the complete, runnable Python/NumPy training script for MNIST based on these notes."/>

import pandas as pd
import kagglehub

# Download the Kaggle dataset
path = kagglehub.dataset_download("oddrationale/mnist-in-csv")

train_path = f"{path}/mnist_train.csv"
test_path = f"{path}/mnist_test.csv"
train_df = pd.read_csv(train_path)
test_df = pd.read_csv(test_path)

x_train = train_df.drop('label', axis=1).values / 255.0
y_train = train_df['label'].values
x_test = test_df.drop('label', axis=1).values / 255.0
y_test = test_df['label'].values

np.random.seed(7882)

W1 = np.random.rand(128, 784) * 0.01  # Scale from [0,1) to [0, 0.01)
b1 = np.random.rand(128, 1) * 0.01
W2 = np.random.rand(10, 128) * 0.01
b2 = np.random.rand(10, 1) * 0.01

# W1 = np.random.randn(128, 784) * np.sqrt(2.0 / 784)
# b1 = np.zeros((128, 1))
# W2 = np.random.randn(10, 128) * np.sqrt(2.0 / 128)
# b2 = np.zeros((10, 1))


def relu(x):
    return np.maximum(0, x)


def softmax(logits):
    # res = np.zeros(10)
    # for i, logit in enumerate(logits.ravel()):
    #     res[i] = np.exp(logit) / np.sum(np.exp(logits))
    # return res

    # Shift logits by max to avoid overflow
    shifted = logits - np.max(logits)
    exp_scores = np.exp(shifted)
    return exp_scores / np.sum(exp_scores)


def forward(x):
    # Ensure x is a column vector
    # in our dataloader, x.shape is (784, ) which's 1d numpy array
    if x.ndim == 1:
        x = x.reshape(-1, 1)

    Z1 = W1 @ x + b1
    A1 = relu(Z1)
    Z2 = W2 @ A1 + b2
    A2 = softmax(Z2)
    return A2, A1, Z1


def backward(x, y):
    if x.ndim == 1:
        x = x.reshape(-1, 1)

    pred, A1, Z1 = forward(x)
    y_one_hot = np.zeros_like(pred)
    y_one_hot[y] = 1.0

    loss = -np.log(pred[y])
    dZ2 = pred - y_one_hot  # dim(dZ2) = (10,)
    dZ2 = dZ2.reshape(-1, 1)

    dW2 = dZ2 @ A1.T  # (10, 128)
    db2 = dZ2

    # Hidden Layer
    dA1 = W2.T @ dZ2  # (128,1)
    dZ1 = dA1 * (Z1 > 0)  # (128,1) element-wise
    dW1 = dZ1 @ x.T  # (128,784)
    db1 = dZ1  # (128,1)
    return dW1, db1, dW2, db2, loss


def gdescent(lr=0.01):
    global W1, b1, W2, b2
    for i in range(5):
        total_loss = 0
        for sample_i in range(len(x_train)):
            x = x_train[sample_i]
            y = y_train[sample_i]
            dW1, db1, dW2, db2, loss = backward(x, y)
            W1 = W1 - dW1 * lr
            b1 = b1 - db1 * lr
            W2 = W2 - dW2 * lr
            b2 = b2 - db2 * lr
            total_loss += loss

        print(f"Epoch {i} loss: {total_loss / len(x_train)}")


import time

# Start timer
start_time = time.time()
gdescent()
# End timer
end_time = time.time()
training_time = end_time - start_time
print(f"Total training time: {training_time:.2f} seconds")

correct = 0

for i in range(x_test.shape[0]):
    probs, _, _ = forward(x_test[i])
    pred_label = np.argmax(probs)
    ground_truth = y_test[i]
    if pred_label == ground_truth:
        correct += 1

print(f"Acc: {correct}/{x_test.shape[0]} -> {correct / x_test.shape[0]}")