Backpropagation Is Not Magic. It Is a Chain of Jacobian Matrices. Here Are All of Them.

You have trained a thousand models. You have never seen the actual matrix that flows backward through a layer. Until now.

— -

You call loss.backward().

Gradients appear. Weights update. The model gets better.

But what just happened? What actually flowed backward through the network?

Not numbers. Not scalars. Matrices.

At every layer, backpropagation computes a Jacobian — a matrix of partial derivatives that tells you how every output of that layer changes with respect to every input. Then it chains these matrices together using the chain rule.

The chain rule you learned in calculus class was for scalars: dy/dx = dy/du × du/dx.

The chain rule in a neural network is for vectors: dL/dx = J₃ × J₂ × J₁ — a product of Jacobian matrices, one per layer.

Nobody shows you these matrices. Let me show you every single one.

— -

## The Network

The tiniest possible network that still shows you everything.

Input: x = [x₁, x₂] (2 neurons)

Hidden: h = [h₁, h₂, h₃] (3 neurons, ReLU activation)

Output: y = [y₁] (1 neuron, no activation)

Loss: L = (y₁ — target)² (MSE)

Two layers of weights:

W₁ (3×2): transforms input → hidden

b₁ (3×1): hidden bias

W₂ (1×3): transforms hidden → output

b₂ (1×1): output bias

Let’s assign concrete values.

import numpy as np

# Weights

W1 = np.array([[ 0.3, -0.5],

[ 0.8, 0.2],

[-0.1, 0.6]]) # shape (3, 2)

b1 = np.array([0.1, -0.2, 0.05]) # shape (3,)

W2 = np.array([[0.4, -0.3, 0.7]]) # shape (1, 3)

b2 = np.array([0.1]) # shape (1,)

# Input and target

x = np.array([1.0, 2.0]) # shape (2,)

target = 1.0

— -

## Forward Pass (By Hand)

Layer 1: Linear

z₁ = W₁ @ x + b₁

z₁[0] = 0.3×1.0 + (-0.5)×2.0 + 0.1 = 0.3–1.0 + 0.1 = -0.6

z₁[1] = 0.8×1.0 + 0.2×2.0 + (-0.2) = 0.8 + 0.4–0.2 = 1.0

z₁[2] = (-0.1)×1.0 + 0.6×2.0 + 0.05 = -0.1 + 1.2 + 0.05 = 1.15

z₁ = [-0.6, 1.0, 1.15]

Layer 2: ReLU

h = ReLU(z₁)

h[0] = max(0, -0.6) = 0.0 ← killed

h[1] = max(0, 1.0) = 1.0 ← alive

h[2] = max(0, 1.15) = 1.15 ← alive

h = [0.0, 1.0, 1.15]

Neuron 0 is dead. Its gradient will be zero. This is ReLU doing its job — sparsifying the representation.

Layer 3: Linear

z₂ = W₂ @ h + b₂

z₂ = 0.4×0.0 + (-0.3)×1.0 + 0.7×1.15 + 0.1

= 0.0–0.3 + 0.805 + 0.1

= 0.605

y = 0.605

Loss:

L = (y — target)² = (0.605–1.0)² = (-0.395)² = 0.1560

Now let’s go backward. This is where the Jacobians live.

— -

## Backward Pass: The Jacobian at Each Layer

The chain rule says:

dL/dW₁ = dL/dy × dy/dh × dh/dz₁ × dz₁/dW₁

Each of those “d/d” terms is a Jacobian matrix. Let’s compute every one.

— -

### Jacobian 1: Loss → Output

dL/dy = 2(y — target) = 2(0.605–1.0) = -0.790

This one is a scalar because both L and y are scalars. Simple.

— -

### Jacobian 2: Output → Hidden (The Linear Layer Jacobian)

The output is y = W₂ @ h + b₂. How does y change when h changes?

dy/dh = W₂ = [0.4, -0.3, 0.7]

This is a (1×3) Jacobian. Each entry tells you: if I nudge h[j] by a tiny amount, how much does y change?

• Nudge h[0] by +0.01 → y changes by +0.004 (weight 0.4)

• Nudge h[1] by +0.01 → y changes by −0.003 (weight −0.3)

• Nudge h[2] by +0.01 → y changes by +0.007 (weight 0.7)

The Jacobian of a linear layer is just the weight matrix. This is why linear layers are simple to backpropagate through.

— -

### Jacobian 3: ReLU (The Diagonal Jacobian)

h = ReLU(z₁). How does h change when z₁ changes?

ReLU is element-wise. Each h[i] depends only on z₁[i]. So the Jacobian is diagonal:

dh/dz₁ = diag(ReLU’(z₁))

ReLU’(z₁[0]) = ReLU’(-0.6) = 0 (input was negative → gradient is 0)

ReLU’(z₁[1]) = ReLU’(1.0) = 1 (input was positive → gradient is 1)

ReLU’(z₁[2]) = ReLU’(1.15) = 1 (input was positive → gradient is 1)

dh/dz₁ = [[0, 0, 0],

[0, 1, 0],

[0, 0, 1]]

This is a (3×3) diagonal matrix. The zero in position (0,0) means: no gradient flows through the dead neuron. This is the gradient-killing property of ReLU, visible as a zero on the diagonal of the Jacobian.

— -

### Jacobian 4: Hidden Linear → Input

z₁ = W₁ @ x + b₁. How does z₁ change when x changes?

dz₁/dx = W₁ = [[ 0.3, -0.5],

[ 0.8, 0.2],

[-0.1, 0.6]]

A (3×2) Jacobian. Each row tells you how one hidden neuron responds to changes in the input.

— -

## Chaining the Jacobians

Now multiply them. The chain rule for vectors:

dL/dx = dL/dy × dy/dh × dh/dz₁ × dz₁/dx

Let’s go step by step.

Step 1: dL/dh = dL/dy × dy/dh

dL/dh = -0.790 × [0.4, -0.3, 0.7]

= [-0.316, 0.237, -0.553]

Shape: (1×3). This is the gradient of the loss with respect to the hidden activations.

Step 2: dL/dz₁ = dL/dh × dh/dz₁

dL/dz₁ = [-0.316, 0.237, -0.553] × [[0, 0, 0],

[0, 1, 0],

[0, 0, 1]]

= [0.0, 0.237, -0.553]

The dead neuron (z₁[0] = −0.6) zeroed out. Its gradient is exactly 0.0. The signal that was flowing backward through the network hit the dead ReLU and stopped.

This is the Jacobian chain in action. One zero on the diagonal of the ReLU Jacobian killed one dimension of the gradient. Information that could have flowed to the input along that path is gone.

Step 3: dL/dx = dL/dz₁ × dz₁/dx

dL/dx = [0.0, 0.237, -0.553] × [[ 0.3, -0.5],

[ 0.8, 0.2],

[-0.1, 0.6]]

dL/dx[0] = 0.0×0.3 + 0.237×0.8 + (-0.553)×(-0.1)

= 0.0 + 0.1896 + 0.0553

= 0.2449

dL/dx[1] = 0.0×(-0.5) + 0.237×0.2 + (-0.553)×0.6

= 0.0 + 0.0474–0.3318

= -0.2844

dL/dx = [0.2449, -0.2844]

That is the gradient of the loss with respect to the input. In a deeper network, this would continue flowing backward to the previous layer.

— -

## The Weight Gradients (What Actually Updates)

You do not update inputs. You update weights. The weight gradients use the same Jacobian chain, but stop at the weight matrix.

Gradient for W₂:

dL/dW₂ = dL/dy × dy/dW₂

dy/dW₂ = hᵀ = [0.0, 1.0, 1.15]

dL/dW₂ = -0.790 × [0.0, 1.0, 1.15]

= [0.0, -0.790, -0.909]

Gradient for W₁:

dL/dW₁ = (dL/dz₁)ᵀ × xᵀ

dL/dz₁ = [0.0, 0.237, -0.553] (column vector)

x = [1.0, 2.0] (row vector)

dL/dW₁ = [[0.0], × [1.0, 2.0]

[0.237],

[-0.553]]

= [[0.0, 0.0 ],

[0.237, 0.474 ],

[-0.553, -1.106]]

Look at row 0 of dL/dW₁: all zeros. The dead ReLU neuron receives zero gradient. Its weights will not update. It is frozen until a different input wakes it up.

— -

## The Complete Code

import numpy as np

# === Network setup ===

W1 = np.array([[ 0.3, -0.5],

[ 0.8, 0.2],

[-0.1, 0.6]])

b1 = np.array([0.1, -0.2, 0.05])

W2 = np.array([[0.4, -0.3, 0.7]])

b2 = np.array([0.1])

x = np.array([1.0, 2.0])

target = 1.0

# === Forward pass ===

z1 = W1 @ x + b1

h = np.maximum(0, z1) # ReLU

z2 = W2 @ h + b2

y = z2[0]

loss = (y — target) ** 2

print(“=== FORWARD PASS ===”)

print(f”z1 (pre-ReLU): {z1}”)

print(f”h (post-ReLU): {h}”)

print(f”y (output): {y:.4f}”)

print(f”Loss: {loss:.4f}”)

print()

# === Backward pass (Jacobian chain) ===

# Jacobian 1: dL/dy

dL_dy = 2 * (y — target)

print(“=== BACKWARD PASS ===”)

print(f”dL/dy = {dL_dy:.4f}”)

# Jacobian 2: dy/dh = W2

J2 = W2 # shape (1, 3)

dL_dh = dL_dy * J2 # shape (1, 3)

print(f”dy/dh (= W2): {J2}”)

print(f”dL/dh: {np.round(dL_dh, 4)}”)

# Jacobian 3: dh/dz1 = diag(ReLU’)

relu_grad = (z1 > 0).astype(float)

J3 = np.diag(relu_grad) # shape (3, 3)

dL_dz1 = dL_dh @ J3 # shape (1, 3)

print(f”ReLU gradients: {relu_grad} ← neuron 0 is DEAD”)

print(f”dL/dz1: {np.round(dL_dz1, 4)}”)

# Jacobian 4: dz1/dx = W1

J4 = W1 # shape (3, 2)

dL_dx = dL_dz1 @ J4 # shape (1, 2)

print(f”dL/dx: {np.round(dL_dx, 4)}”)

print()

# === Weight gradients ===

dL_dW2 = dL_dy * h.reshape(1, -1) # outer product

dL_db2 = np.array([dL_dy])

dL_dW1 = dL_dz1.T @ x.reshape(1, -1) # outer product

dL_db1 = dL_dz1.flatten()

print(“=== WEIGHT GRADIENTS ===”)

print(f”dL/dW2:\n{np.round(dL_dW2, 4)}”)

print(f”dL/db2: {np.round(dL_db2, 4)}”)

print(f”dL/dW1:\n{np.round(dL_dW1, 4)}”)

print(f”dL/db1: {np.round(dL_db1, 4)}”)

print()

print(“Note: Row 0 of dL/dW1 is all zeros — the dead ReLU neuron gets no gradient update.”)

# === Verify with PyTorch ===

print(“\n=== PYTORCH VERIFICATION ===”)

try:

import torch

import torch.nn as nn

W1_t = torch.tensor(W1, dtype=torch.float64, requires_grad=True)

b1_t = torch.tensor(b1, dtype=torch.float64, requires_grad=True)

W2_t = torch.tensor(W2, dtype=torch.float64, requires_grad=True)

b2_t = torch.tensor(b2, dtype=torch.float64, requires_grad=True)

x_t = torch.tensor(x, dtype=torch.float64)

z1_t = W1_t @ x_t + b1_t

h_t = torch.relu(z1_t)

y_t = (W2_t @ h_t + b2_t)[0]

loss_t = (y_t — target) ** 2

loss_t.backward()

print(f”PyTorch dL/dW1:\n{np.round(W1_t.grad.numpy(), 4)}”)

print(f”Our dL/dW1:\n{np.round(dL_dW1, 4)}”)

print(f”Match: {np.allclose(W1_t.grad.numpy(), dL_dW1)}”)

except ImportError:

print(“(PyTorch not available — manual computation stands on its own)”)

— -

## What You Just Saw

Four Jacobian matrices. One per operation.

dL/dy → scalar (loss derivative)

dy/dh → W₂ (1×3) (linear layer = its own Jacobian)

dh/dz₁ → diag (3×3) (ReLU = diagonal mask)

dz₁/dx → W₁ (3×2) (linear layer = its own Jacobian)

The backward pass multiplied them together, right to left:

dL/dx = dL/dy × W₂ × diag(ReLU’) × W₁

= scalar × (1×3) × (3×3) × (3×2)

= (1×2)

Every loss.backward() you have ever called did exactly this. Jacobian after Jacobian, multiplied in sequence, from the loss back to the first layer.

The dead neuron showed up as a zero on the diagonal. The weight gradient for that neuron was zero. No magic. Just a matrix with a zero in the right place.

That is backpropagation. Not an algorithm. A matrix product.

— -

References

Rumelhart, D.E., Hinton, G.E., and Williams, R.J. (1986). Learning representations by back-propagating errors. Nature, 323, 533–536. https://doi.org/10.1038/323533a0

Goodfellow, I., Bengio, Y., and Courville, A. (2016). Deep Learning, Chapter 6.5: Back-Propagation. MIT Press. https://www.deeplearningbook.org/