Backpropagation becomes much easier to use—and debug—when you understand what the framework is calculating for you. This guide shows how to implement backpropagation without PyTorch or TensorFlow using only Python and NumPy. You will build a two-layer classifier, derive its gradients with the chain rule, train it in mini-batches, and verify the implementation numerically.
This is an educational implementation, not a replacement for production frameworks. For real systems, automatic differentiation, GPU kernels, mixed precision, and distributed training are valuable. But a small NumPy network is an excellent foundation for understanding model behaviour and for prototyping custom learning rules. The same discipline is useful when building scalable machine learning pipelines in Python or neural models for domain-specific datasets such as Indian agriculture data.
What backpropagation actually computes
A neural network is a composition of functions. For a two-layer network:
- Hidden pre-activation:
z1 = XW1 + b1 - Hidden activation:
a1 = ReLU(z1) - Output logits:
z2 = a1W2 + b2 - Predicted probabilities:
p = softmax(z2) - Loss: cross-entropy between
pand the target labels
Backpropagation applies the chain rule from the loss backwards. It produces the derivative of the loss with respect to every parameter, after which gradient descent updates each parameter:
parameter = parameter - learning_rate × gradient
The forward pass must cache intermediate values such as z1, a1, and p. Without those values, the backward pass would need to recompute the graph and would be harder to inspect.
Set up a small NumPy network
Install NumPy in a virtual environment:
python -m venv .venv
source .venv/bin/activate # Windows: .venv\\Scripts\\activate
pip install numpyThe example below classifies two-dimensional points into two classes. It uses a hidden layer with ReLU and a two-unit softmax output. Unlike the original sigmoid-and-mean-squared-error pattern, this combination is a standard choice for classification and gives a particularly clean output gradient.
import numpy as np
rng = np.random.default_rng(42)
# XOR-like training data
X = np.array([
[0., 0.], [0., 1.], [1., 0.], [1., 1.],
[0., 0.2], [0.2, 0.], [1., 0.8], [0.8, 1.]
])
y = np.array([0, 1, 1, 0, 0, 0, 1, 1])
n_inputs, n_hidden, n_classes = 2, 8, 2
W1 = rng.normal(0, np.sqrt(2 / n_inputs), (n_inputs, n_hidden))
b1 = np.zeros((1, n_hidden))
W2 = rng.normal(0, np.sqrt(2 / n_hidden), (n_hidden, n_classes))
b2 = np.zeros((1, n_classes))He-style initialisation keeps ReLU activations at a useful scale. Biases start at zero, while weights should not: identical weights would make hidden neurons learn the same function.
Write stable activation and loss functions
Numerical stability matters even in a teaching implementation. Subtracting the largest logit before exponentiation prevents overflow in exp.
def relu(z):
return np.maximum(0, z)
def softmax(logits):
shifted = logits - np.max(logits, axis=1, keepdims=True)
exp_values = np.exp(shifted)
return exp_values / np.sum(exp_values, axis=1, keepdims=True)
def cross_entropy(probs, targets):
n = targets.shape[0]
# Clip only for the logarithm; probabilities still come from softmax.
return -np.mean(np.log(probs[np.arange(n), targets] + 1e-12))ReLU has derivative 1 when z > 0 and 0 otherwise. At exactly zero, either convention is acceptable for this example; the code uses zero.
Implement the forward pass
The function returns both the prediction and the cache needed by backpropagation:
def forward(X, W1, b1, W2, b2):
z1 = X @ W1 + b1
a1 = relu(z1)
logits = a1 @ W2 + b2
probs = softmax(logits)
cache = (X, z1, a1, probs)
return probs, cacheFor an input matrix with shape (batch_size, features), W1 has shape (features, hidden_units) and W2 has shape (hidden_units, classes). Checking these shapes early prevents many silent broadcasting errors.
Derive and code the backward pass
With softmax followed by cross-entropy, the derivative of the loss with respect to the logits is:
d_logits = (probs - one_hot_targets) / batch_size
The remaining gradients follow the matrix chain rule:
dW2 = a1.T @ d_logitsdb2 = sum(d_logits, axis=0)da1 = d_logits @ W2.Tdz1 = da1 * (z1 > 0)dW1 = X.T @ dz1db1 = sum(dz1, axis=0)
def backward(y, cache, W2):
X, z1, a1, probs = cache
batch_size = X.shape[0]
one_hot = np.zeros_like(probs)
one_hot[np.arange(batch_size), y] = 1
d_logits = (probs - one_hot) / batch_size
dW2 = a1.T @ d_logits
db2 = np.sum(d_logits, axis=0, keepdims=True)
da1 = d_logits @ W2.T
dz1 = da1 * (z1 > 0)
dW1 = X.T @ dz1
db1 = np.sum(dz1, axis=0, keepdims=True)
return dW1, db1, dW2, db2Notice that the gradients are calculated separately from the parameter update. This makes the code easier to test and avoids accidental use of already-updated weights during the same backward pass.
Train with mini-batch gradient descent
def accuracy(probs, y):
return np.mean(np.argmax(probs, axis=1) == y)
learning_rate = 0.08
batch_size = 8
for epoch in range(3000):
indices = rng.permutation(len(X))
for start in range(0, len(X), batch_size):
batch_idx = indices[start:start + batch_size]
X_batch, y_batch = X[batch_idx], y[batch_idx]
probs, cache = forward(X_batch, W1, b1, W2, b2)
dW1, db1, dW2, db2 = backward(y_batch, cache, W2)
W1 -= learning_rate * dW1
b1 -= learning_rate * db1
W2 -= learning_rate * dW2
b2 -= learning_rate * db2
if epoch % 500 == 0:
probs, _ = forward(X, W1, b1, W2, b2)
print(f"epoch={epoch:4d} loss={cross_entropy(probs, y):.4f} "
f"accuracy={accuracy(probs, y):.2%}")A healthy run should show loss declining and accuracy improving. Do not expect every run to produce identical numbers: random initialisation and batch order affect optimisation. Keep the random seed fixed while debugging, then vary it to assess robustness.
Verify gradients before trusting training
A decreasing loss is not proof that every gradient is correct. Use finite differences to compare an analytical gradient with a numerical approximation:
def numerical_gradient(param, loss_fn, epsilon=1e-5):
gradient = np.zeros_like(param)
iterator = np.nditer(param, flags=["multi_index"], op_flags=["readwrite"])
while not iterator.finished:
index = iterator.multi_index
original = param[index]
param[index] = original + epsilon
plus = loss_fn()
param[index] = original - epsilon
minus = loss_fn()
param[index] = original
gradient[index] = (plus - minus) / (2 * epsilon)
iterator.iternext()
return gradientCompare only a few entries on a small batch. The relative error should usually be close to machine precision for a correct implementation:
abs(analytical - numerical) / max(1, abs(analytical), abs(numerical))
Common causes of failure include forgetting the division by batch size, using the wrong transpose, applying the ReLU derivative to a1 instead of z1, and modifying a parameter before all gradients are computed.
Practical improvements and limitations
This implementation is intentionally small. Before using it for a real project, add:
- Feature scaling: standardise numeric inputs and handle missing values explicitly.
- Validation data: monitor validation loss, not just training accuracy, and stop when it deteriorates.
- Regularisation: add L2 penalties or dropout when overfitting appears.
- Better data loading: shuffle training data and use reproducible train-validation-test splits.
- Monitoring: log loss, accuracy, gradient norms, and parameter norms.
- Unit tests: test each activation, loss, shape, and gradient independently.
For larger models, manual NumPy backpropagation becomes slow and error-prone. Frameworks are appropriate when you need convolutional layers, GPUs, automatic differentiation, or deployment tooling. For example, a production object-detection project may eventually move to custom object detection models with PyTorch. The point of this exercise is not to avoid those tools permanently; it is to know what their training engine is doing.
The same forward-cache-backward-update pattern extends to deeper networks, embeddings, regression losses, and custom layers. It also gives AI teams a stronger basis for reviewing model behaviour before integrating it into enterprise workflows such as private LLMs for faculty research data.
FAQ
Why avoid PyTorch or TensorFlow for this exercise?
Because writing the derivatives yourself exposes matrix shapes, cached activations, numerical stability, and the chain rule. Use a framework once the learning goal becomes a production requirement.
Can this code train more than one hidden layer?
Yes. Store each layer's pre-activation and activation in a list, then traverse the list in reverse order. Each layer passes its gradient to the preceding layer.
Why use cross-entropy instead of mean squared error?
Cross-entropy with softmax is generally better suited to classification and produces a simple, well-scaled gradient. Mean squared error remains useful for regression and for studying basic calculus.
What should I try after this example?
Add a validation split, implement L2 regularisation, compare ReLU with tanh, and write gradient checks for every parameter matrix. Then compare your results with a framework implementation on the same data.
Apply for AI Grants India
If you are building an AI product, research prototype, or India-focused deployment, apply to AI Grants India for potential funding, mentorship, and ecosystem support. A clear implementation plan— including data, evaluation, compute requirements, and responsible deployment—will make your proposal more credible.