Backpropagation — Why Your 10-Layer Network Stops Learning

Every time you unlock your phone with your face, or get a eerily accurate Netflix recommendation, or watch GPT-4 complete your sentence — backpropagation is the algorithm that made those models smart. It's the engine inside every gradient-based deep learning model ever trained. Without it, neural networks would be untrained noise generators, not intelligent systems. It's not an exaggeration to say backpropagation is the most important algorithm in modern AI.

The core problem backpropagation solves is credit assignment: when a network of thousands or millions of parameters makes a wrong prediction, which parameters are responsible, and by how much? Tweaking weights randomly is computationally hopeless. You need an efficient, mathematically principled way to propagate blame backwards through a computational graph — assigning each weight a gradient that tells you exactly which direction to push it. Backpropagation does this in a single backwards pass using the chain rule of calculus, turning what would be an O(n²) problem into O(n).

By the end of this article you'll understand the chain rule derivation from scratch, implement forward and backward passes without any autograd framework, recognize the vanishing and exploding gradient problems and know exactly what causes them at the weight-initialization level, and understand what PyTorch's autograd is actually doing under the hood when you call loss.backward(). You'll also walk away with the kind of nuanced understanding that separates candidates who get ML engineering offers from those who don't.

What is Backpropagation Explained?

Backpropagation is essentially the application of the Chain Rule from calculus to a directed acyclic graph (DAG) of computations. In a neural network, we calculate the partial derivative of the loss function $L$ with respect to every weight $w$ in the network. Mathematically, for a single weight at layer $l$, the gradient is:

$$\frac{\partial L}{\partial w^{(l)}} = \frac{\partial L}{\partial a^{(l)}} \cdot \frac{\partial a^{(l)}}{\partial z^{(l)}} \cdot \frac{\partial z^{(l)}}{\partial w^{(l)}}$$

By caching these intermediate derivatives during the backward pass, we avoid redundant calculations, making deep learning computationally feasible.

import numpy as np

# io.thecodeforge implementation of a basic neuron backprop
def forge_backprop_demo():
    # 1. Forward Pass
    x = 1.0       # Input
    w = -2.0      # Weight
    b = 3.0       # Bias
    
    # Linear transformation: z = w*x + b
    z = w * x + b 
    # Activation: a = tanh(z)
    a = np.tanh(z)
    
    # 2. Backward Pass (Manual Gradient Calculation)
    # dL/da = 1.0 (assuming this is the end of the graph)
    da = 1.0
    
    # da/dz = 1 - tanh^2(z)
    dz = (1 - a**2) * da
    
    # dz/dw = x
    dw = x * dz
    # dz/db = 1
    db = 1.0 * dz
    
    print(f"Activation: {a:.4f}")
    print(f"Weight Gradient: {dw:.4f}")
    print(f"Bias Gradient: {db:.4f}")

if __name__ == "__main__":
    forge_backprop_demo()

The Chain Rule in Depth: From Single Neuron to Multilayer Networks

When you stack layers, the chain rule chains together. For a two-layer network with hidden layer $h = f(W_1 x + b_1)$ and output $\hat{y} = f(W_2 h + b_2)$, the gradient of loss $L$ with respect to $W_1$ is:

$$\frac{\partial L}{\partial W_1} = \frac{\partial L}{\partial \hat{y}} \cdot \frac{\partial \hat{y}}{\partial h} \cdot \frac{\partial h}{\partial z_1} \cdot \frac{\partial z_1}{\partial W_1}$$

Each new layer multiplies in an extra derivative. The key insight: we compute these from output back to input, reusing intermediate values. In code, this means storing activations ($h$, $\hat{y}$) during the forward pass. Let's extend our earlier example to two layers.

import numpy as np

def forge_two_layer_backprop():
    # Two-layer network: input -> hidden (tanh) -> output (linear)
    x = np.array([1.0, 2.0])  # Input (2-dim)
    W1 = np.random.randn(3, 2) * 0.5  # 3 hidden neurons
    b1 = np.zeros(3)
    W2 = np.random.randn(1, 3) * 0.5
    b2 = np.zeros(1)
    
    # Forward pass
    z1 = W1.dot(x) + b1          # (3,)
    h = np.tanh(z1)              # activation
    z2 = W2.dot(h) + b2          # (1,)
    y_pred = z2                  # linear output
    loss = 0.5 * (y_pred - 1.0)**2  # MSE with target=1

    # Backward pass (dL/dW1, dL/dW2, etc.)
    dloss = y_pred - 1.0         # dL/dy_pred
    dz2 = dloss                  # linear: dL/dz2 = dL/dy_pred
    dW2 = dz2 * h                # (1,3)
    db2 = dz2                    # (1,)
    dh = W2.T.dot(dz2)           # (3,) backprop through linear
    dz1 = (1 - h**2) * dh        # tanh derivative
    dW1 = np.outer(dz1, x)       # (3,2)
    db1 = dz1
    
    print(f"Loss: {loss:.6f}")
    print(f"grad W1 shape: {dW1.shape}, norm: {np.linalg.norm(dW1):.4f}")
    print(f"grad W2 shape: {dW2.shape}, norm: {np.linalg.norm(dW2):.4f}")

forge_two_layer_backprop()

Vanishing and Exploding Gradients: The Deep Network Killer

As depth increases, gradients can either shrink to zero (vanishing) or grow exponentially (exploding). Vanishing happens when activation derivatives are small (sigmoid max 0.25, tanh max 1.0 but saturates). With 10 layers using sigmoid, the gradient can shrink to $0.25^{10} \approx 9.5 \times 10^{-7}$. Exploding occurs with poor weight initialization – large weights compound multiplicatively.

This isn't just a theory – it's the reason deep learning didn't work before ReLU, batch norm, and He initialization. The 2015 paper 'Delving Deep into Rectifiers' showed that proper initialization alone can make 30-layer networks trainable.

Let's simulate both conditions to see the effect.

import numpy as np

def simulate_gradient_flow(num_layers=10, activation='sigmoid', init_std=1.0):
    x = np.random.randn(100)
    W = [np.random.randn(100, 100) * init_std for _ in range(num_layers)]
    b = [np.zeros(100) for _ in range(num_layers)]
    
    # Forward with stored activations
    a = x
    activations = [a]
    for i in range(num_layers):
        z = W[i].dot(a) + b[i]
        if activation == 'sigmoid':
            a = 1/(1+np.exp(-z))
        elif activation == 'tanh':
            a = np.tanh(z)
        else:  # relu
            a = np.maximum(0, z)
        activations.append(a)
    
    # Backward
    grad = np.ones(100)
    norm_per_layer = []
    for i in reversed(range(num_layers)):
        z = W[i].dot(activations[i]) + b[i]
        if activation == 'sigmoid':
            da = activations[i+1] * (1 - activations[i+1])
        elif activation == 'tanh':
            da = 1 - activations[i+1]**2
        else:
            da = (activations[i+1] > 0).astype(float)
        grad = (W[i].T.dot(grad * da))
        norm_per_layer.insert(0, np.linalg.norm(grad))
    return norm_per_layer

norms_sigmoid = simulate_gradient_flow(num_layers=10, activation='sigmoid')
norms_relu = simulate_gradient_flow(num_layers=10, activation='relu')
print("Gradient norms per layer (sigmoid):", [f"{n:.3e}" for n in norms_sigmoid])
print("Gradient norms per layer (ReLU):", [f"{n:.3e}" for n in norms_relu])

Numerical Stability: When Your Gradients Become NaN or Inf

Floating-point arithmetic is finite precision. During backprop, you'll often compute log, exp, or division operations that can produce NaNs or infinities. Common culprits:

• Cross-entropy loss: $L = -\log(\hat{y})$ when $\hat{y} = 0$ produces $\infty$
• Softmax: $e^{z_i}$ where $z_i$ is large (e.g., 1000) causes overflow
• Sigmoid: $1 / (1 + e^{-z})$ for $z \approx -1000$ underflows to 0
• Division by a small gradient can produce Inf

Production engineers learn to use numerically stable alternatives: the log-softmax function, adding epsilons to denominators, and gradient clipping. Let's see how PyTorch handles this and where it can still fail.

import torch
import torch.nn.functional as F

def stable_cross_entropy(logits, target):
    # Logits: (batch, classes), target: (batch,)
    # Using built-in which is numerically stable
    loss = F.cross_entropy(logits, target)
    return loss

# Demonstrate unstable vs stable
unstable_logits = torch.tensor([[1000.0, -1000.0, 0.0]])
target = torch.tensor([0])

try:
    # Manual unstable softmax and log
    softmax = torch.exp(unstable_logits) / torch.exp(unstable_logits).sum(dim=1, keepdim=True)
    loss_manual = -torch.log(softmax[0, 0])
    print(f"Manual loss: {loss_manual.item()}")  # Likely inf or nan
except Exception as e:
    print(f"Manual calculation failed: {e}")

# Stable version
loss_stable = stable_cross_entropy(unstable_logits, target)
print(f"Stable loss: {loss_stable.item():.4f}")

# Backward works
loss_stable.backward()
print(f"Gradient exists: {unstable_logits.grad is not None}")

Autograd: What Happens When You Call loss.backward()

PyTorch's autograd builds a computational graph on the fly. When you call loss.backward(), it traverses the graph in reverse, computing gradients using the chain rule. The graph is dynamic – it's rebuilt every forward pass. Under the hood:

1. GradFn nodes store pointers to input tensors and the operation
2. Each tensor has a .grad attribute that accumulates gradients
3. After backward, the graph is freed by default (to save memory)
4. Setting retain_graph=True keeps it for multiple backward calls

TensorFlow uses a similar model but with a static graph (eager mode makes it dynamic). The key difference: PyTorch's tape is imperative – you can put Python control flow inside the model. TensorFlow's graph mode compiles the whole graph first, which enables optimisations but makes debugging harder.

Let's inspect a simple autograd graph to see what's happening.

import torch

x = torch.tensor([3.0], requires_grad=True)
w = torch.tensor([2.0], requires_grad=True)
b = torch.tensor([1.0], requires_grad=True)

z = w * x + b
a = torch.tanh(z)
loss = (a - 0.5)**2

print("Forward values:")
print(f"z = {z.item():.4f}")
print(f"a = {a.item():.4f}")
print(f"loss = {loss.item():.4f}")

# Backward
loss.backward()

print("\nGradients:")
print(f"dl/dx: {x.grad.item():.4f}")
print(f"dl/dw: {w.grad.item():.4f}")
print(f"dl/db: {b.grad.item():.4f}")

# Inspect the graph nodes
print("\nComputation graph nodes:")
def print_grad_fn(tensor, depth=0):
    if tensor.grad_fn:
        print(' ' * depth + f'{tensor.grad_fn.__class__.__name__}: {tensor.grad_fn}')
        for next_fn, _ in tensor.grad_fn.next_functions:
            if next_fn:
                print(next_fn.grad_fn.__class__.__name__ if hasattr(next_fn, 'grad_fn') else str(next_fn))
    else:
        print(' ' * depth + 'leaf tensor')

print_grad_fn(loss)