Build a Neural Network from Scratch in Python

Every major deep learning framework (PyTorch 2.10, JAX 0.9, TensorFlow) hides thousands of lines of optimized C++ and CUDA behind a single model.fit() call. That abstraction is wonderful for shipping products, but terrible for understanding what actually happens when a neural network learns. If you can't build one from scratch in NumPy, you're flying blind when debugging vanishing gradients, choosing learning rates, or diagnosing why your model won't converge.

In this guide, we'll build a fully functional multi-layer neural network, piece by piece, and train it to classify handwritten digits from the sklearn digits dataset (1,797 images of 8x8 pixels). Every line of code runs directly in your browser. By the end, you'll have a working classifier that hits 97.78% accuracy, and you'll understand every matrix multiplication and gradient computation that makes it possible.

The Single Neuron as a Building Block

Before we build a full network, we need to understand the smallest unit: a single neuron.

Think of a neuron like a voting machine. It takes in a bunch of numbers (inputs), assigns each one an importance score (weight), adds them all up, and then decides whether to "fire" or stay quiet. That's it. Every neural network, from a simple perceptron to GPT-class large language models, is built from millions of these tiny voters working together.

Here's the math for what a single neuron computes:

$z = \sum_{i=1}^{n} w_i x_i + b = \mathbf{w}^T \mathbf{x} + b$

Where:

• $\mathbf{w}$ is the weight vector (one weight per input feature)
• $\mathbf{x}$ is the input vector (e.g., 64 pixel values from an 8x8 digit image)
• $b$ is the bias term, which shifts the decision boundary
• $z$ is the pre-activation value (also called the logit)

The next step is the activation function. Without it, stacking layers would be pointless because multiple linear transformations collapse into a single linear transformation (matrix multiplication is associative). The activation function introduces nonlinearity, which is what gives neural networks the ability to learn curved decision boundaries instead of just straight lines.

For hidden layers, ReLU remains the default choice in 2026. It simply outputs the input if it's positive, and zero if it's negative. For the output layer of a multi-class classifier, we use softmax, which squashes a vector of raw scores into probabilities that sum to 1. Deep dives into activation choices are covered in Activation Functions: ReLU, Sigmoid, and Beyond.

Click to expandSingle neuron computing a weighted sum of pixel inputs through activation to produce output

Weight Initialization Determines Whether Training Succeeds

Now that we know what a single neuron does, let's talk about how we set up the weights before training begins. This might sound like a minor detail, but it's actually the difference between a network that learns in 50 epochs and one that never converges at all.

Why does it matter so much? Consider three scenarios:

All weights set to zero. Every neuron in a layer computes the same output. During backpropagation (which we'll cover soon), every neuron receives the same gradient. They all update identically, forever. Your 128-neuron layer effectively behaves like a single neuron. This is called the symmetry problem.

Weights too large. The outputs grow exponentially as they pass through layers. By layer 5, your numbers overflow to inf. Gradients explode, and training diverges.

Weights too small. The outputs shrink toward zero at each layer. By layer 10, your signal is essentially zero. Gradients vanish, and the early layers stop learning entirely.

The solution is to initialize weights from a carefully chosen random distribution. Two strategies dominate in practice:

Since our network uses ReLU in hidden layers, we'll use He initialization (He et al., 2015):

$W \sim \mathcal{N}\left(0,\ \sqrt{\frac{2}{n_{in}}}\right)$

Where:

• $W$ is the weight matrix for a given layer
• $n_{in}$ is the number of neurons in the previous layer (fan-in)
• $\mathcal{N}$ denotes a normal (Gaussian) distribution

Here's our initialization code. Notice that we set np.random.seed(42) for reproducibility, so you'll get the exact same weights every time you run this:

W1 shape: (64, 128), std: 0.1782
W2 shape: (128, 64), std: 0.1248
W3 shape: (64, 10), std: 0.1728

Let's unpack what just happened. The architecture = [64, 128, 64, 10] defines our network shape: 64 input features (one per pixel), then a hidden layer with 128 neurons, then another hidden layer with 64 neurons, and finally 10 output neurons (one per digit class, 0 through 9).

Notice how the standard deviation adapts to each layer's fan-in. The first layer ($n_{in}=64$) gets a std of 0.1782, while the second layer ($n_{in}=128$) gets a smaller std of 0.1248. More inputs means each individual weight should be smaller, so their combined effect doesn't blow up the output.

Forward Pass Through Multiple Layers

Now that our weights are initialized, we need a way to push data through the network and get predictions out the other end. This is called the forward pass.

Think of it like an assembly line. Raw materials (pixel values) enter at one end. Each station (layer) transforms them: multiply by weights, add bias, apply activation. The finished product (a probability distribution over 10 digit classes) comes out the other end.

For our digit classifier, an 8x8 image (flattened to 64 values) enters the network and produces 10 output probabilities:

Click to expandForward pass data flow from input pixels through hidden layers to digit class probabilities

Let's implement this. We need three functions: relu for the hidden layer activations, softmax for the output layer, and forward to wire everything together.

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

def softmax(z):
    # Subtract max for numerical stability
    exp_z = np.exp(z - np.max(z, axis=1, keepdims=True))
    return exp_z / np.sum(exp_z, axis=1, keepdims=True)

def forward(X, params):
    """Forward pass through our 3-layer network."""
    cache = {'A0': X}

    # Hidden layer 1: Linear -> ReLU
    cache['Z1'] = X @ params['W1'] + params['b1']
    cache['A1'] = relu(cache['Z1'])

    # Hidden layer 2: Linear -> ReLU
    cache['Z2'] = cache['A1'] @ params['W2'] + params['b2']
    cache['A2'] = relu(cache['Z2'])

    # Output layer: Linear -> Softmax
    cache['Z3'] = cache['A2'] @ params['W3'] + params['b3']
    cache['A3'] = softmax(cache['Z3'])

    return cache['A3'], cache

Let's walk through what each line does:

• X @ params['W1'] + params['b1'] — This is the linear transformation. The @ operator is matrix multiplication. For each sample, it computes the weighted sum of all inputs for every neuron in layer 1. The result Z1 contains the raw (pre-activation) values.
• relu(cache['Z1']) — This applies the ReLU activation: keep positive values as-is, replace negatives with zero. The result A1 is what gets passed to the next layer.
• The pattern repeats for layer 2, except now the input is A1 (the output of layer 1) instead of X.
• The output layer uses softmax instead of relu. Softmax converts the raw logits into probabilities that sum to 1.0, so we can interpret them as "the network thinks this is a 7 with 85% confidence."

The softmax function deserves a closer look. Subtracting np.max(z) before exponentiating doesn't change the output probabilities (the subtraction cancels in the numerator and denominator), but it prevents overflow when logits are large. Without this trick, np.exp(1000) returns inf, which would propagate NaN through the entire computation.

Cross-Entropy Loss Measures Prediction Quality

We can now push data through our network and get predictions. But how do we know if those predictions are any good? We need a way to measure the gap between what the network predicts and what the true labels are. That measurement is the loss function.

Cross-entropy loss is the standard choice for classification tasks. It quantifies how far the predicted probability distribution sits from the true labels. The reason it's preferred over MSE (mean squared error) for classification: cross-entropy produces much stronger gradients when predictions are confidently wrong, which means faster learning.

$\mathcal{L} = -\frac{1}{N} \sum_{i=1}^{N} \sum_{c=1}^{C} y_{i,c} \log(\hat{y}_{i,c})$

Where:

• $N$ is the number of samples in the batch
• $C$ is the number of classes (10 digits in our example)
• $y_{i,c}$ is 1 if sample $i$ belongs to class $c$, otherwise 0 (one-hot encoded)
• $\hat{y}_{i,c}$ is the predicted probability for sample $i$ belonging to class $c$

def cross_entropy_loss(y_pred, y_true):
    """Compute cross-entropy loss with numerical stability."""
    N = y_true.shape[0]
    # Clip predictions to avoid log(0)
    y_pred_clipped = np.clip(y_pred, 1e-12, 1 - 1e-12)
    loss = -np.sum(y_true * np.log(y_pred_clipped)) / N
    return loss

def one_hot_encode(y, num_classes):
    """Convert integer labels to one-hot vectors."""
    one_hot = np.zeros((len(y), num_classes))
    one_hot[np.arange(len(y)), y] = 1
    return one_hot

The one_hot_encode function converts integer labels (like 7) into vectors (like [0,0,0,0,0,0,0,1,0,0]). This format matches what our softmax layer outputs: a probability for each class. We need both in the same format to compute the loss.

Backpropagation Computes Gradients Layer by Layer

So far we can push data forward through the network and measure how bad our predictions are. The next question is: how do we improve? Which weights should increase, which should decrease, and by how much?

This is the job of backpropagation. It uses the chain rule from calculus to compute, for every single weight in the network, the answer to: "If I nudge this weight up by a tiny amount, how much does the loss change?" That answer is the gradient, and it tells us the direction to adjust each weight to reduce the loss.

For a thorough mathematical treatment, see Backpropagation: Engine of Deep Learning. Here, we'll focus on the practical implementation.

The gradient of cross-entropy loss combined with softmax simplifies to something beautifully simple:

$\frac{\partial \mathcal{L}}{\partial \mathbf{Z}^{[L]}} = \hat{\mathbf{Y}} - \mathbf{Y}$

Where:

• $\mathbf{Z}^{[L]}$ is the pre-activation output of the last layer
• $\hat{\mathbf{Y}}$ is the predicted probability matrix (from softmax)
• $\mathbf{Y}$ is the one-hot encoded true label matrix

For hidden layers, we propagate the gradient backward through the chain rule:

$\frac{\partial \mathcal{L}}{\partial \mathbf{W}^{[l]}} = \frac{1}{N} \mathbf{A}^{[l-1]T} \cdot \delta^{[l]}$

Where:

• $\mathbf{A}^{[l-1]}$ is the activation from the previous layer (the input to the current layer)
• $\delta^{[l]}$ is the error signal at layer $l$
• $N$ is the batch size

Now let's implement it:

def backward(y_pred, y_true, cache, params):
    """Backpropagation through our 3-layer network."""
    N = y_true.shape[0]
    grads = {}

    # Output layer gradient (softmax + cross-entropy shortcut)
    dZ3 = (y_pred - y_true) / N
    grads['dW3'] = cache['A2'].T @ dZ3
    grads['db3'] = np.sum(dZ3, axis=0, keepdims=True)

    # Hidden layer 2
    dA2 = dZ3 @ params['W3'].T
    dZ2 = dA2 * (cache['Z2'] > 0)  # ReLU derivative
    grads['dW2'] = cache['A1'].T @ dZ2
    grads['db2'] = np.sum(dZ2, axis=0, keepdims=True)

    # Hidden layer 1
    dA1 = dZ2 @ params['W2'].T
    dZ1 = dA1 * (cache['Z1'] > 0)  # ReLU derivative
    grads['dW1'] = cache['A0'].T @ dZ1
    grads['db1'] = np.sum(dZ1, axis=0, keepdims=True)

    return grads

Let's trace through the key lines:

1. dZ3 = (y_pred - y_true) / N — The output gradient. Predicted minus actual, divided by batch size to get the average.
2. grads['dW3'] = cache['A2'].T @ dZ3 — How much each weight in layer 3 contributed to the error. We multiply the transposed activations from layer 2 (what went into layer 3) by the error signal.
3. dA2 = dZ3 @ params['W3'].T — We propagate the error backward through layer 3's weights to find how much layer 2's output contributed to the error.
4. dZ2 = dA2 * (cache['Z2'] > 0) — The ReLU derivative. This is elegantly simple: 1 where the pre-activation input was positive, 0 where it was negative. If a neuron was "off" (ReLU zeroed it out), it gets zero gradient, meaning it contributed nothing to the error and won't be updated.
5. The pattern repeats identically for layer 1.

This is why we saved those intermediate values in cache during the forward pass. Without A2, A1, Z2, and Z1, we couldn't compute any of these gradients.

The Training Loop Ties Everything Together

We now have all four building blocks: weight initialization, forward pass, loss computation, and backpropagation. The training loop is the engine that ties them all together.

Here's the rhythm: grab a batch of training data, push it through the network (forward pass), measure how wrong the predictions are (loss), compute which direction to adjust each weight (backward pass), then actually adjust the weights (parameter update). Repeat this thousands of times, and the network gradually gets better.

For a deeper look at optimizer choices beyond vanilla SGD, see Deep Learning Optimizers: SGD to AdamW.

Click to expandTraining loop pipeline from data batching through forward pass, loss computation, backpropagation, and weight update

def train(X_train, y_train, X_val, y_val, architecture, epochs=200,
          lr=0.1, batch_size=64):
    """Train neural network with mini-batch gradient descent."""
    params = initialize_weights(architecture)
    y_train_oh = one_hot_encode(y_train, 10)
    N = X_train.shape[0]
    history = {'train_loss': [], 'val_acc': []}

    for epoch in range(epochs):
        # Shuffle training data each epoch
        indices = np.random.permutation(N)
        X_shuffled = X_train[indices]
        y_shuffled = y_train_oh[indices]

        epoch_loss = 0
        num_batches = 0

        for start in range(0, N, batch_size):
            end = min(start + batch_size, N)
            X_batch = X_shuffled[start:end]
            y_batch = y_shuffled[start:end]

            # Forward pass
            y_pred, cache = forward(X_batch, params)

            # Compute loss
            loss = cross_entropy_loss(y_pred, y_batch)
            epoch_loss += loss
            num_batches += 1

            # Backward pass
            grads = backward(y_pred, y_batch, cache, params)

            # Update weights (vanilla SGD)
            for key in params:
                params[key] -= lr * grads[f'd{key}']

        # Track metrics
        avg_loss = epoch_loss / num_batches
        val_pred, _ = forward(X_val, params)
        val_acc = np.mean(np.argmax(val_pred, axis=1) == y_val)
        history['train_loss'].append(avg_loss)
        history['val_acc'].append(val_acc)

        if (epoch + 1) % 50 == 0:
            print(f"Epoch {epoch+1:3d} | Loss: {avg_loss:.4f} | Val Acc: {val_acc:.4f}")

    return params, history

A few things worth highlighting:

Why shuffle every epoch? Without shuffling, the network sees the same batch ordering each time, which can create oscillating gradients that prevent convergence. Think of it like studying for an exam: reviewing topics in a different order each day produces better retention than the same sequence every time.

Why mini-batches instead of the full dataset? Using all 1,437 training samples at once gives you a precise gradient, but it's slow (one update per epoch). Using one sample at a time gives you a noisy gradient, but you get 1,437 updates per epoch. Mini-batches (64 samples here) are the sweet spot: reasonably accurate gradients with frequent updates.

The weight update rule is vanilla SGD (stochastic gradient descent): params[key] -= lr * grads[f'd{key}']. This says "move each weight in the opposite direction of its gradient, scaled by the learning rate." If the gradient says "increasing this weight increases the loss," we decrease the weight. The learning rate (0.1 here) controls the step size.

Putting It All Together on Real Data

Let's train our network on the sklearn digits dataset. This is a clean, manageable dataset (1,797 samples, 10 classes, 64 features) that's perfect for validating a from-scratch implementation before moving on to larger challenges.

The code below is fully self-contained. It includes every function we've built so far, loads the data, preprocesses it, trains the network, and prints the results. You can run it directly in your browser:

Training samples: 1437
Validation samples: 360
Features: 64, Classes: 10
Epoch  50 | Loss: 0.0038 | Val Acc: 0.9722
Epoch 100 | Loss: 0.0015 | Val Acc: 0.9722
Epoch 150 | Loss: 0.0009 | Val Acc: 0.9778
Epoch 200 | Loss: 0.0006 | Val Acc: 0.9778

Final validation accuracy: 0.9778

97.78% validation accuracy with a from-scratch implementation. Let's break down what the output tells us:

The loss drops dramatically. From the first epoch to epoch 50, the loss already falls to 0.0038. By epoch 200, it's at 0.0006. This means the network is getting extremely confident about its correct predictions (remember, cross-entropy loss approaches zero when the predicted probability for the correct class approaches 1.0).

The accuracy plateaus around epoch 150. We reach 97.78% at epoch 150 and stay there through epoch 200. This means additional training isn't helping (or hurting). In a larger project, you'd use early stopping to save the model at this point and avoid wasting compute.

8 out of 360 validation samples are misclassified. For 8x8 pixel images, some digits genuinely look ambiguous (a sloppy 4 can look like a 9), so ~98% is excellent for this dataset.

How Our NumPy Network Compares to PyTorch

The question everyone asks: "How does a from-scratch implementation compare to a real framework?" Here's the same architecture expressed in PyTorch for comparison:

import torch
import torch.nn as nn
import torch.optim as optim

class DigitClassifier(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(64, 128),
            nn.ReLU(),
            nn.Linear(128, 64),
            nn.ReLU(),
            nn.Linear(64, 10)
        )

    def forward(self, x):
        return self.net(x)

model = DigitClassifier()
criterion = nn.CrossEntropyLoss()
optimizer = optim.SGD(model.parameters(), lr=0.1)

The PyTorch version is ~20 lines vs. ~80 for our NumPy implementation. But both express identical math: the same matrix multiplications, the same ReLU gating, the same cross-entropy gradients. On a small dataset like digits (1,797 samples), accuracy will be comparable regardless of framework because the underlying algorithm is the same.

PyTorch's real advantage isn't accuracy on toy problems. It's automatic differentiation (torch.autograd), GPU acceleration, and torch.compile (which in PyTorch 2.10 fuses operations for up to 2x training speedups). For production models with millions of parameters trained on GPU clusters (or when applying transfer learning from pretrained checkpoints), these matter enormously. For understanding what those frameworks actually do, building from scratch matters more.

When to Build from Scratch vs. Use a Framework

This decision comes up in interviews, coursework, and early-stage prototyping. Here's a practical framework:

Build from scratch when:

• Learning fundamentals (university courses, self-study, interview prep)
• Debugging framework behavior (you suspect a bug in your model, not the framework)
• Implementing novel architectures not supported by existing layers
• Working in constrained environments (embedded, no Python dependencies)
• Teaching others (nothing beats walking through the math and code together)

Use PyTorch/JAX when:

• Training on GPU or TPU (no way around it for large models)
• Any model with > 100K parameters (manual backprop becomes error-prone)
• Production deployment (ONNX export, TorchScript, torch.compile)
• Research requiring automatic differentiation through complex control flow
• Anything involving convolutional networks, recurrent architectures, or attention mechanisms

Production Considerations

Conclusion

Building a neural network from scratch strips away every abstraction and forces you to confront the math directly. You've implemented He initialization, ReLU activations, softmax output, cross-entropy loss, and full backpropagation through three layers. The result, a 97.78% digit classifier built entirely in NumPy, proves that the core algorithm is straightforward once you see each piece clearly.

The same bias-variance tradeoff that governs simpler models applies here: more hidden units increase capacity (lower bias) but risk overfitting (higher variance). Techniques like dropout, batch normalization, and learning rate scheduling, all built on top of the foundation we coded today, address this balance in modern architectures. To validate your model properly, pair this knowledge with cross-validation techniques rather than relying on a single train/test split.

The next step is clear. Pick a problem you care about, swap in your own data, and train. Break things on purpose: try zero initialization, remove the softmax, use a learning rate of 10. Watching a network fail teaches more than watching it succeed. And once you're comfortable with dense layers, move on to the transformer architecture to see how attention mechanisms replaced recurrence entirely.

Interview Questions

Why can't you initialize all weights to zero in a neural network?

All neurons in a layer would compute identical outputs and receive identical gradients during backpropagation. They'd update identically every step, meaning the network effectively has one neuron per layer regardless of width. This is called the symmetry problem. Random initialization breaks this symmetry so each neuron can specialize in detecting a different pattern.

Explain the difference between MSE and cross-entropy loss for classification.

MSE treats outputs as continuous values and penalizes errors quadratically, producing small gradients when predictions are near 0 or 1 (where sigmoid saturates). Cross-entropy directly measures the divergence between predicted and true probability distributions, generating strong gradients even for confidently wrong predictions. For classification, cross-entropy converges faster and reaches better solutions.

What happens if you use sigmoid activations throughout a 20-layer network?

Sigmoid outputs are bounded between 0 and 1, so each layer squashes its inputs. During backpropagation, the sigmoid derivative (maximum value 0.25) multiplies at every layer. After 20 layers, gradients shrink by roughly $0.25^{20} \approx 10^{-12}$, making early layers virtually untrainable. This is the vanishing gradient problem. ReLU avoids it because its derivative is 1 for positive inputs, allowing gradients to flow unchanged through active neurons.

Why does mini-batch gradient descent often outperform full-batch gradient descent?

Mini-batches introduce noise into gradient estimates, which acts as implicit regularization and helps escape shallow local minima. Full-batch gradients are more accurate but point to the nearest local minimum, which may not generalize well. Mini-batches also allow more frequent weight updates per epoch, so the network sees more learning steps in the same wall-clock time. Typical batch sizes range from 32 to 256.

Your neural network's training loss decreases but validation accuracy stops improving after 50 epochs. What's happening?

The network is overfitting to the training data. After 50 epochs, it starts memorizing training examples rather than learning general patterns. Practical remedies include early stopping (stop training at the lowest validation loss), dropout (randomly zeroing neurons during training), L2 regularization on weights, or simply collecting more training data. You should also verify your network capacity isn't too large for the dataset size.

How does the learning rate interact with batch size during training?

Larger batch sizes produce lower-variance gradient estimates, which can tolerate higher learning rates. A common heuristic (the linear scaling rule from Goyal et al., 2017) is to multiply the learning rate by $k$ when you multiply the batch size by $k$. However, this breaks down for very large batch sizes (above ~8K), where specialized techniques like learning rate warmup become necessary.

What's the computational complexity of backpropagation through a fully connected layer?

For a layer with $n_{in}$ inputs and $n_{out}$ outputs and batch size $B$, the forward pass costs $O(B \cdot n_{in} \cdot n_{out})$ for the matrix multiplication. Backpropagation requires two matrix multiplications of similar size (one for the weight gradient, one for the input gradient), so it costs roughly 2x the forward pass. Total training cost per layer per batch is $O(3 \cdot B \cdot n_{in} \cdot n_{out})$.