Building Neural Networks from Scratch in NumPy

Every deep learning framework is hiding the same ten lines of math. Before you reach for PyTorch or TensorFlow, writing a network in raw NumPy forces you to understand what actually happens during training. This post builds a working XOR classifier step by step.

The Math You Actually Need

A neural network is a composition of functions. Each layer applies a linear transformation followed by a non-linearity:

Z = X @ W + b       # linear step
A = activation(Z)   # non-linear step

Training is gradient descent on a loss function. To compute gradients, you apply the chain rule backwards through the computation graph — that's all backpropagation is.

Forward Pass

For XOR, we need at least one hidden layer. A 2→8→1 architecture works:

import numpy as np
 
def sigmoid(z):
    return 1 / (1 + np.exp(-z))
 
def sigmoid_deriv(a):
    return a * (1 - a)
 
def relu(z):
    return np.maximum(0, z)
 
def relu_deriv(z):
    return (z > 0).astype(float)
 
# Xavier initialization
W1 = np.random.randn(2, 8) * np.sqrt(1 / 2)
b1 = np.zeros((1, 8))
W2 = np.random.randn(8, 1) * np.sqrt(1 / 8)
b2 = np.zeros((1, 1))
 
def forward(X):
    Z1 = X @ W1 + b1
    A1 = relu(Z1)
    Z2 = A1 @ W2 + b2
    A2 = sigmoid(Z2)
    cache = (X, Z1, A1, Z2, A2)
    return A2, cache

The cache stores intermediate values — you'll need them during backprop.

Loss Function

Binary cross-entropy for a binary classifier:

def bce_loss(y_pred, y_true):
    eps = 1e-9
    return -np.mean(
        y_true * np.log(y_pred + eps) + (1 - y_true) * np.log(1 - y_pred + eps)
    )

The gradient of BCE with respect to the output (before sigmoid) simplifies to A2 - y, which is why sigmoid + BCE is such a common pairing.

Backward Pass

This is where most tutorials get vague. Every gradient is computed by asking: "how does this variable affect the loss?"

def backward(cache, y_true, lr=0.05):
    global W1, b1, W2, b2
    X, Z1, A1, Z2, A2 = cache
    m = X.shape[0]
 
    # Output layer gradient
    dA2 = (A2 - y_true) / m          # dL/dA2 (BCE + sigmoid shortcut)
    dW2 = A1.T @ dA2
    db2 = np.sum(dA2, axis=0, keepdims=True)
 
    # Hidden layer gradient
    dA1 = dA2 @ W2.T
    dZ1 = dA1 * relu_deriv(Z1)       # chain rule through ReLU
    dW1 = X.T @ dZ1
    db1 = np.sum(dZ1, axis=0, keepdims=True)
 
    # Gradient descent step
    W2 -= lr * dW2
    b2 -= lr * db2
    W1 -= lr * dW1
    b1 -= lr * db1

The key line is dZ1 = dA1 * relu_deriv(Z1). The derivative of ReLU is 1 where the input was positive and 0 where it was negative — gradients don't flow through dead neurons.

Training Loop

X = np.array([[0,0],[0,1],[1,0],[1,1]], dtype=np.float32)
y = np.array([[0],[1],[1],[0]], dtype=np.float32)
 
for epoch in range(3000):
    pred, cache = forward(X)
    loss = bce_loss(pred, y)
    backward(cache, y)
    if epoch % 500 == 0:
        print(f"Epoch {epoch}: loss={loss:.4f}")
 
# Results after training
pred, _ = forward(X)
print(pred.round(2))  # → [[0.02], [0.97], [0.97], [0.03]]

What Neural Nexus Adds

Neural Nexus wraps this pattern in a block system so you don't repeat the boilerplate. Affine(2, 8) is the W1/b1 pair. ReLU() caches its input for backward. NeuralNet.fit() runs the loop. But the math is identical — the library just orchestrates it.

Understanding this from scratch means you can debug a training run when the loss explodes (usually a learning rate or initialization problem), when gradients vanish (usually depth + saturating activations), or when the model memorizes but doesn't generalize (usually lack of regularization).