Handout #5 — Analytical Backpropagation

Handout #4 established losses, gradient descent, and the chain rule. Now we derive an efficient backward sweep for a layered network, then implement and check its parameter gradients.

1Deriving analytical backpropagation

Slides 90–109 · Direct gradients, cached sensitivities, and parameter updates

1.1  Name the forward computations

Use the multilayer network from the slides, with biases made explicit. For one example, let $h_0=x\in\mathbb R^{d_0}$ and, for $i=0,\ldots,L-1$, define

Here $W_i\in\mathbb R^{d_{i+1}\times d_i}$, $b_i\in\mathbb R^{d_{i+1}}$, and $u\in\mathbb R^{d_L}$. The hidden activations $f_i$ act elementwise. Setting biases to zero recovers the slides' forward equations. Cache the activations and the values needed to evaluate activation derivatives during the forward pass.

We use column gradients: $g_i=\nabla_{h_i}\ell$. The slides' cached upstream sensitivities are written as rows; $g_i$ is their transpose. Below, $\delta_{i+1}$ means the column gradient with respect to the preactivation $z_{i+1}$, not $h_{i+1}$.

1.2  Start at the loss and output layer

First compute the scalar $q=\partial\ell/\partial\hat y$. For the convenient squared-error loss $\ell=\tfrac12(\hat y-t)^2$, we have $q=\hat y-t$. Since $\hat y=\sum_j u_j h_{L,j}+c$, the chain rule gives

The first two are parameter gradients to save for the optimizer. The last is the upstream sensitivity to send into the last hidden layer. A different loss changes the starting derivative $q$; the hidden-layer rules below remain the same.

1.3  Derive one hidden layer, then reuse it

Suppose $g_{i+1}$ has already arrived from the output side. The elementwise activation gives

This is a Jacobian-transpose product: the activation Jacobian is diagonal, so multiplication reduces to the elementwise product $\odot$. For sigmoid, $f'(z)=f(z)(1-f(z))$; for tanh, $f'(z)=1-f(z)^2$; for ReLU, it is 1 at positive inputs and 0 at negative inputs (we use 0 at zero).

Now expand the affine operation coordinate by coordinate:

Changing $W_{i,jk}$ directly changes only $z_{i+1,j}$, with local derivative $h_{i,k}$. Changing $h_{i,k}$ affects every output coordinate, so those paths must be summed:

Collecting these entries gives the three reusable rules:

1.4  Why the backward sweep saves work

Let $A_i=J_{h_{i+1}}(h_i)$. Expanding the chain rule independently for every hidden layer would repeatedly compute the same suffix of the network:

Instead, cache $g_{i+1}$ and compute $g_i=A_i^\top g_{i+1}$ once. In practice we evaluate this as $W_i^\top(g_{i+1}\odot f_i'(z_{i+1}))$, without constructing the full Jacobian. The saved upstream gradient also supplies that layer's weight and bias gradients.

With comparable layer widths, independently recomputing all suffixes takes $1+2+\cdots+L=O(L^2)$ layer operations. A backward sweep takes $O(L)$, comparable to the forward pass. For dense layers, a more informative cost is $O(\sum_i d_i d_{i+1})$ per example for either pass, plus elementwise work and the output layer. Cached activations cost memory.

1. Run the forward pass and cache intermediate values.
2. Differentiate the loss and scalar output to initialize $g_L$ and the output parameter gradients.
3. For $i=L-1,\ldots,0$, compute $\delta_{i+1}$, save weight and bias gradients, and pass $g_i$ backward.
4. After all gradients are computed, update each parameter using the chosen optimizer.

1.5  A complete numerical backward pass

Take $x=(1,2)^\top$, target $t=1$, ReLU hidden units, and

The forward pass gives $z=(1,-1)^\top$, $h=(1,0)^\top$, $\hat y=2$, and $\ell=\tfrac12(2-1)^2=0.5$. Before opening the answer, compute every parameter gradient and the gradient with respect to $x$.

2Code lab: implement and check the derivation

2.1  Manual gradients for a minibatch

For row-batched inputs, $Z=XW^\top+b$, $H=\tanh Z$, and $\hat y=Hu+c$. Use $J=\tfrac1{2B}\sum_n(\hat y_n-t_n)^2$. Then $q=(\hat y-t)/B$ already includes averaging. With $D=(q u^\top)\odot(1-H^2)$, the weight gradient is $D^\top X$ and the bias gradient sums the rows of $D$. Do not divide by $B$ again.

This runnable example implements the backward formulas explicitly. Autograd is used only as an independent comparison, and finite differences check every parameter of this small network.

import torch

torch.manual_seed(7)
B, d, h = 3, 2, 4
X = torch.randn(B, d, dtype=torch.float64)
t = torch.randn(B, dtype=torch.float64)
W = (0.2 * torch.randn(h, d, dtype=torch.float64)).requires_grad_()
b = torch.zeros(h, dtype=torch.float64, requires_grad=True)
u = (0.2 * torch.randn(h, dtype=torch.float64)).requires_grad_()
c = torch.zeros((), dtype=torch.float64, requires_grad=True)
params = [W, b, u, c]

def objective():
    H = torch.tanh(X @ W.T + b)
    pred = H @ u + c
    return 0.5 * ((pred - t) ** 2).mean()

# Our backward pass: no autograd operations are used to find derivatives.
with torch.no_grad():
    H = torch.tanh(X @ W.T + b)
    pred = H @ u + c
    q = (pred - t) / B
    D = (q[:, None] * u[None, :]) * (1 - H ** 2)
    manual = [D.T @ X, D.sum(dim=0), H.T @ q, q.sum()]

reference = torch.autograd.grad(objective(), params)
for p, actual, expected in zip(params, manual, reference):
    assert actual.shape == p.shape
    torch.testing.assert_close(actual, expected, rtol=1e-10, atol=1e-12)

# Central differences: deliberately expensive, useful for a tiny check.
eps = 1e-6
with torch.no_grad():
    for p, expected in zip(params, manual):
        numeric = torch.empty_like(p)
        flat = p.view(-1)
        for k in range(p.numel()):
            old = flat[k].item()
            flat[k] = old + eps
            plus = objective().item()
            flat[k] = old - eps
            minus = objective().item()
            flat[k] = old
            numeric.view(-1)[k] = (plus - minus) / (2 * eps)
        torch.testing.assert_close(numeric, expected, rtol=1e-5, atol=1e-8)
    before = objective().item()
    for p, grad in zip(params, manual):
        p -= 0.01 * grad
    after = objective().item()
assert after < before  # checked for this example, not a general guarantee
print("All gradients agree; loss:", before, "->", after)

2.2  Connecting to a classification output

The scalar regression head makes the derivation easy to follow. For classification, replace it with logits $s=Vh+a$, probabilities $p=\operatorname{softmax}(s)$, and a one-hot target $y$. For one example, cross-entropy can be rewritten as

Thus the output sensitivity is $r=p-y$, and the same affine-layer derivation gives $\nabla_V\ell=rh^\top$, $\nabla_a\ell=r$, and $g_h=V^\top r$. For a mean loss over a batch, divide each example's output sensitivity by the batch size before summing contributions. The hidden-layer recurrence is unchanged.

2.3  Code lab: manual gradients for a classifier

Keep the tanh hidden layer, but replace the scalar output with $K$ class logits. For a batch of $B$ examples, let $X\in\mathbb R^{B\times d}$, $W\in\mathbb R^{h\times d}$, $V\in\mathbb R^{K\times h}$, and let $Y\in\mathbb R^{B\times K}$ contain one-hot targets:

Apply §2.2 to every row. The mean reduction enters once, at the logits:

Bias gradients sum over the rows that shared the bias. The input gradient is useful as a check; it does not mean we update the training inputs. Predict each shape before running this independent example. Manual derivatives use tensor arithmetic under no_grad(); autograd is only the reference calculation.

import torch
import torch.nn.functional as F

torch.manual_seed(11)
B, d, h, K = 4, 3, 5, 3
X = torch.randn(B, d, dtype=torch.float64, requires_grad=True)
target = torch.tensor([0, 2, 1, 2])
W = (0.2 * torch.randn(h, d, dtype=torch.float64)).requires_grad_()
b = torch.zeros(h, dtype=torch.float64, requires_grad=True)
V = (0.2 * torch.randn(K, h, dtype=torch.float64)).requires_grad_()
a = torch.zeros(K, dtype=torch.float64, requires_grad=True)
params = [W, b, V, a]

def objective():
    H = torch.tanh(X @ W.T + b)
    return F.cross_entropy(H @ V.T + a, target)

with torch.no_grad():
    H = torch.tanh(X @ W.T + b)
    logits = H @ V.T + a
    P = logits.softmax(dim=-1)
    Y = F.one_hot(target, num_classes=K).to(dtype=X.dtype)
    R = (P - Y) / B
    D = (R @ V) * (1 - H.square())
    manual = [D.T @ X, D.sum(0), R.T @ H, R.sum(0)]
    input_grad = D @ W
    torch.testing.assert_close(R.sum(-1), torch.zeros(B, dtype=X.dtype))

reference = torch.autograd.grad(objective(), params + [X])
for tensor, actual, expected in zip(params + [X], manual + [input_grad], reference):
    assert actual.shape == tensor.shape
    torch.testing.assert_close(actual, expected, rtol=1e-10, atol=1e-12)

# Also compare a directional derivative with two loss evaluations.
directions = [torch.ones_like(p) for p in params]
directions[0] = torch.arange(W.numel(), dtype=W.dtype).reshape_as(W) / W.numel()
analytic = sum((g * v).sum().item() for g, v in zip(manual, directions))
eps = 1e-6
with torch.no_grad():
    originals = [p.clone() for p in params]
    for p, original, v in zip(params, originals, directions):
        p.copy_(original + eps * v)
    plus = objective().item()
    for p, original, v in zip(params, originals, directions):
        p.copy_(original - eps * v)
    minus = objective().item()
    for p, original in zip(params, originals):
        p.copy_(original)
assert abs((plus - minus) / (2 * eps) - analytic) < 1e-7
print("Classifier: every parameter and input gradient agrees with autograd.")

The directional check compares $\sum_p\langle\nabla_pJ,v_p\rangle$ with $[J(\theta+\epsilon v)-J(\theta-\epsilon v)]/(2\epsilon)$. It is a cheap additional check of one direction, not a replacement for checking every component as in §2.1.