Backprop by Hand: Nine Gradients and a Verified XOR Network

0 · Backprop is a local sensitivity calculation

Backpropagation applies the chain rule to the computation that actually ran, reusing intermediate derivatives rather than perturbing every weight separately. Rumelhart, Hinton and Williams's 1986 Nature article concerns learning internal representations; it occupies pages 533–536. The original claims about mailed floppy disks, a particular workstation and an eight-page letter were not needed or supported. Karpathy's 2019 training recipe is a different article from his 2016 essay about understanding backprop.

Here we recover the nine-parameter ReLU derivation from DL009 and preserve the independently executable 17-parameter tanh XOR experiment. They illustrate the same chain-rule contract with different activation derivatives; they are not interchangeable numerical examples.

1 · What you'll build today

Derive every weight and bias derivative of the 2-2-1 ReLU network; test every coordinate with centered differences; implement a stable logit loss; train a small nonlinear network; compare NumPy to CPU PyTorch at exactly copied parameters; and explain sum/mean reduction, shared edges, checkpointing and accumulation. No fixed seed-independent training time or accuracy is promised.

2 · Checkpoint · you should already know

Bring the chain rule, broadcasting, sigmoid/BCE cancellation and the forward cache (X,Z1,H,z,p). Review DL004, DL007, DL008 and DL009. Optional pre-reading: backprop intuition, calculus, micrograd lecture and the training recipe above.

3 · Recap — the forward pass from S009

For N examples and one binary output, use the stable mean loss

This equals binary cross entropy in exact arithmetic without log(0). Adding epsilon inside log(p+epsilon) changes the objective and its derivative; it is not an exact match to the fused p-y formula.

4 · The four cheat-sheet formulas — derived from coordinates

For Z=X@W+b and upstream G=dL/dZ:

The first contracts the batch axis; the second contracts the output-feature axis. Hence dW=X.T@G, dX=G@W.T, db=G.sum(0). ReLU multiplies by (Z>0), selecting zero as the derivative convention at zero. For sigmoid multiply by p*(1-p); for tanh by 1-H**2. The fused mean BCE endpoint gives (p-y)/N; multiclass normalized-target CE gives (P-target)/N under the same simple mean convention. Weighted/ignored targets can change the denominator.

The factory-blame analogy means signed local sensitivity, not moral fault or a guaranteed finite-step improvement. Every gradient has its variable's shape, but shape checks do not catch missing factors, wrong data or a transpose hidden by square dimensions.

5 · Backward pass — the full derivation, walked

5.1 Output logit and reduction

Differentiate -y log p-(1-y)log(1-p) with respect to p, then multiply by p(1-p); the terms simplify to p-y. Divide by N once for mean loss. For a sum, remove that division; for a scaled loss, scale every derivative identically.

5.2 Output affine layer

Set dz=(p-y)/N. Then dW2=H.T@dz, db2=dz.sum(0), and dH=dz@W2.T. Use the forward W2; do not update it before computing earlier gradients.

5.3 ReLU node

Set dZ1=dH*(Z1>0). For finite real values from the same forward pass, (H>0) and (Z1>0) are equivalent. The original supposed boundary disagreement was false. A mask from another batch, or a different activation, is not equivalent. Zero activation on one example does not mean the unit is dead on the entire dataset.

5.4 Input affine layer

Compute dW1=X.T@dZ1, db1=dZ1.sum(0) and, if needed, dX=dZ1@W1.T. Input gradients matter for saliency, adversarial optimization and input reconstruction even though the input is not an optimizer parameter. This network uses three parameter/backprop matmuls plus one extra for dX; other architectures add different local derivatives.

6 · By-hand computation on one example

Use DL009's x=[1,2], hidden values [.8,1.1], output logit1.26 and target1. With d=0.7790261077798122−1=−0.2209738922201878:

These are the nine parameter derivatives, plus two optional input derivatives. At rate .1, W2[0,0] becomes .5−.1*(−.1767791138)=.5176779114. A derivative predicts the first-order change for a small perturbation; it is not itself the new weight. In our (in,hidden) layout, a hidden unit corresponds to a column of W1, not a row.

7 · Executable ReLU backward and all-coordinate checks

Run this CPU-only program in order. It includes a frozen-snapshot finite difference for all nine parameters, reduction controls and a deliberately wrong gradient. Centered differences have truncation and roundoff error; choose float64 and an h range away from ReLU kinks rather than demanding a universal 1e-8 result.

The helper accepts one binary output per example: logits (N,1) and targets (N,) or (N,1), with nonempty batches and targets in [0,1]. Convert only the explicitly supported vector layout; blindly reshaping a (1,N) target can hide an axis mistake. Reject an unknown reduction instead of silently treating a misspelled mean as sum. A multi-label extension needs its own reduction contract (DL007), not just a wider final matrix.

import numpy as np
 
def binary_targets(y, logits):
    if logits.ndim != 2 or logits.shape[1] != 1 or len(logits) == 0:
        raise ValueError('nonempty (N,1) logits required')
    y = np.asarray(y)
    if y.shape == (len(logits),):
        y = y[:, None]
    if y.shape != logits.shape or not np.isfinite(y).all() or not ((y >= 0) & (y <= 1)).all():
        raise ValueError('targets must have shape (N,) or (N,1) and lie in [0,1]')
    return y
 
def relu_eval(X, y, params, reduction='mean'):
    if reduction not in ('mean', 'sum'):
        raise ValueError('reduction must be mean or sum')
    W1,b1,W2,b2 = params
    Z1 = X@W1+b1
    H = np.maximum(0,Z1)
    z = H@W2+b2
    y = binary_targets(y, z)
    e = np.exp(-np.abs(z))
    prob = np.where(z>=0, 1/(1+e), e/(1+e))
    terms = np.maximum(z,0)-y*z+np.log1p(np.exp(-np.abs(z)))
    denom = len(X) if reduction == 'mean' else 1
    dz = (prob-y)/denom
    dZ1 = (dz@W2.T)*(Z1>0)
    grads = (X.T@dZ1,dZ1.sum(0),H.T@dz,dz.sum(0))
    assert all(g.shape==p.shape for g,p in zip(grads,params))
    return float(terms.sum()/denom),grads,prob,dZ1@W1.T
 
R = (np.array([[.1,.2],[.3,.4]]),np.array([.1,.1]),
     np.array([[.5],[.6]]),np.array([.2]))
RX, Ry = np.array([[1.,2.]]), np.array([[1.]])
loss, grad, prob, dx = relu_eval(RX,Ry,R)
d = 1/(1+np.exp(-1.26))-1
expected = (np.array([[.5*d,.6*d],[d,1.2*d]]),
            np.array([.5*d,.6*d]), np.array([[.8*d],[1.1*d]]),np.array([d]))
assert all(np.allclose(g,e,atol=1e-14) for g,e in zip(grad,expected))
errors=[]
for h in (1e-4,1e-5,1e-6):
    for k,p0 in enumerate(R):
        for ix in np.ndindex(p0.shape):
            plus,minus=[p.copy() for p in R],[p.copy() for p in R]
            plus[k][ix]+=h;minus[k][ix]-=h
            n=(relu_eval(RX,Ry,plus)[0]-relu_eval(RX,Ry,minus)[0])/(2*h)
            errors.append(abs(n-grad[k][ix]))
assert max(errors)<1e-7
assert not np.allclose(grad[0]*2,expected[0])
BX,BY=np.repeat(RX,5,axis=0),np.repeat(Ry,5,axis=0)
gmean=relu_eval(BX,BY,R)[1];gsum=relu_eval(BX,BY,R,'sum')[1]
assert all(np.allclose(a,b) for a,b in zip(gmean,grad))
assert all(np.allclose(a,5*b) for a,b in zip(gsum,grad))
print('nine-gradient check',max(errors),'loss',loss,'dX',dx)
# Rejection controls: these must fail before a loss or gradient is returned.
for bad_y, reduction in [(BY.T, 'mean'), (BY, 'meen'), (BY * 2, 'mean')]:
    try:
        relu_eval(BX, bad_y, R, reduction)
    except ValueError:
        pass
    else:
        raise AssertionError('invalid target layout/value or reduction accepted')
assert np.allclose(relu_eval(BX, BY.ravel(), R)[2], relu_eval(BX, BY, R)[2])

8 · Numeric gradient check — diagnostic limits

For each coordinate use (L(theta+h)-L(theta-h))/(2h) with fresh copies, identical inputs and deterministic operations. Report max absolute error and a relative error with a nonzero floor. A ReLU kink can make a centered slope disagree with the selected subgradient even when backward is correct; try nearby smooth points. Too-small h creates subtraction cancellation, especially in float32/bfloat16. Verify the mathematical implementation in float64, then separately test the training dtype's stability; the original advice to perform tight finite differences directly in bf16 was misleading.

A zero gradient can be correct at a saturated point or a disconnected path. Test several batch sizes and inputs, not only one lucky example. Use a wrong denominator, missing shared-edge contribution and early W2 update as negative controls: the checker must reject each relevant defect.

9 · Full training loop — the protected tanh XOR experiment

The current verified core is retained exactly below. Its widths are 2→4→1, giving 17 parameters; tanh changes the local derivative from a ReLU mask to 1-H*H. The seed-specific fit is an implementation/expressivity test on all four truth-table rows, not held-out generalization.

import numpy as np
 
rng = np.random.default_rng(3)
X = np.array([[0., 0.], [0., 1.], [1., 0.], [1., 1.]])
y = np.array([[0.], [1.], [1.], [0.]])
p = {'W1': rng.normal(0, .7, (2, 4)), 'b1': np.zeros(4),
     'W2': rng.normal(0, .7, (4, 1)), 'b2': np.zeros(1)}
 
def evaluate(p):
    H = np.tanh(X @ p['W1'] + p['b1'])
    z = H @ p['W2'] + p['b2']
    e = np.exp(-np.abs(z))
    prob = np.where(z >= 0, 1/(1+e), e/(1+e))
    loss = np.mean(np.maximum(z, 0) - y*z + np.log1p(np.exp(-np.abs(z))))
    dz = (prob-y) / len(X)
    da = (dz @ p['W2'].T) * (1-H*H)
    grads = {'W2': H.T @ dz, 'b2': dz.sum(axis=0),
             'W1': X.T @ da, 'b1': da.sum(axis=0)}
    return float(loss), grads, prob
 
initial, grads, _ = evaluate(p)
errors = []
for name, value in p.items():
    for index in np.ndindex(value.shape):
        old = value[index]
        value[index] = old + 1e-5
        plus = evaluate(p)[0]
        value[index] = old - 1e-5
        minus = evaluate(p)[0]
        value[index] = old
        errors.append(abs((plus-minus)/2e-5 - grads[name][index]))
assert max(errors) < 1e-7
 
for _ in range(5000):
    loss, grads, prob = evaluate(p)
    for name in p:
        p[name] -= .1 * grads[name]
final, _, prob = evaluate(p)
assert final < .03 and ((prob >= .5) == y).all()
print('parameters', sum(v.size for v in p.values()), 'gradient error', max(errors))
print('loss', round(initial, 6), '->', round(final, 6), 'probabilities', prob.ravel())

Its original smaller hand trace remains useful: at x=(1,0), W1=identity, zero biases and W2=(1,−1)ᵀ, H=(tanh(1),0), z≈.761594 and p≈.681700. Target1 gives dz≈−.318300, dW2[0]≈−.24242 and dZ1[0]≈−.13368. Do not compare those values to §6: both parameters and activation differ.

10 · Forward and backward together

Two sweeps, one cache: each local rule consumes its upstream derivative and saved forward values. Fan-out adds contributions. A residual u=f(x)+x sends both du directly to x and du through f; replacing addition by assignment drops a path.

11 · NumPy to PyTorch — copied parameters, not matching seeds

Matching NumPy and PyTorch seeds does not generate identical weights. Copy each tensor explicitly, transpose for nn.Linear, use float64 and compare the same loss reduction. PyTorch rebuilds an eager graph during forward and accumulates leaf gradients during backward. A scalar output supplies a unit cotangent implicitly; it does not normally populate nonleaf loss.grad. Non-scalar outputs can accept an explicit cotangent.

import torch
import torch.nn.functional as F
torch.set_num_threads(1)
tp=[torch.tensor(a,dtype=torch.float64,requires_grad=True) for a in R]
tx=torch.tensor(RX,dtype=torch.float64);ty=torch.tensor(Ry,dtype=torch.float64)
tz=torch.relu(tx@tp[0]+tp[1])@tp[2]+tp[3]
tl=F.binary_cross_entropy_with_logits(tz,ty);tl.backward()
assert abs(tl.item()-relu_eval(RX,Ry,R)[0])<1e-14
assert all(np.allclose(t.grad.numpy(),g,atol=1e-14) for t,g in zip(tp,grad))
a=torch.tensor(3.,requires_grad=True)
(a*a+a).backward()
assert a.grad.item()==7
(a*a+a).backward() # fresh graph, same leaf, accumulates
assert a.grad.item()==14

12 · Custom kernels, checkpointing and accumulation

For scaled attention, S=QKᵀ/sqrt(d), P=softmax(S) and O=PV. Given dO, derive dV=PᵀdO, dP=dOVᵀ, dS=P*(dP−sum(P*dP,last_axis)), dQ=dSK/sqrt(d) and dK=dSᵀQ/sqrt(d). These are the same affine contractions plus the coupled softmax pullback, not just an elementwise activation. FlashAttention-2 reorganizes work and tiling to reduce non-matmul work and memory movement; its backward recomputes blocks rather than storing the full probability matrix. It is not a proof that all derivatives fit into one universally faster kernel.

Checkpointing discards selected saved activations and recomputes during backward. For an ideal uniform L-layer chain with segment length k, storing boundaries plus one segment costs roughly L/k+k activation units, minimized near sqrt(L). Parameters, optimizer states, nonuniform layers and temporary workspaces remain. Preserve RNG behavior when recomputing dropout; CPU checkpoint tests do not establish GPU memory savings.

Microbatch accumulation retains one optimizer step after all pieces. With mean losses of sizes n_j, accumulate (n_j/N)*L_j, not L_j/number_of_pieces when sizes differ. BatchNorm and stochastic state can break full-batch equivalence. Recompute may cost extra arithmetic; smaller microbatches change utilization and launches, so neither method has a universal wall-clock winner.

13 · Failure cases and repair rules

A missing /N multiplies the gradient by N at the same parameter snapshot, not necessarily every future step's loss improvement. A dead unit has zero positive-activation fraction across the relevant data, not merely a fraction below .5; five dead units out of 64 is 7.8125%, not 5%. A mask from H or Z is equivalent for the same ReLU; a stale mask is wrong. A zero row of W1's gradient may reflect an unused input rather than a dead neuron.

Without clearing leaf gradients between optimizer updates, new gradients add to old ones evaluated at different parameters; this is not a legitimate giant batch at one fixed point. Clearing must respect intentional accumulation. Store-all, checkpoint and microbatch strategies should be chosen from measured peak memory and time, not solely accelerator utilization.

14 · Practice — solutions and controlled variants

1. All nine derivatives: §7 checks all W1,b1,W2,b2 coordinates across three h values. Extend to widths 3 and 5 with nonsquare batches and avoid zero preactivations.
2. Two moons: generate local make_moons, split before fitting, compare a logistic baseline and hidden16 MLP on the same data/seeds, and plot their held-out boundaries. A curved boundary is possible; 85% versus99% is not a guaranteed result. Inspect class overlap caused by noise.
3. Break it: a very large rate can overshoot; removing ReLU yields an affine classifier incapable of perfectly fitting XOR but it can classify three of four points, so “always50%” is false. All-zero ReLU parameters with zero derivative at zero give no hidden gradient; tanh's issue is symmetry, not that same mask.
4. Three-layer backward: compute dz, (dW3,db3), dH2, dZ2, (dW2,db2), dH1, dZ1, (dW1,db1) using the same contraction rules. “Eight formulas” is a grouping choice, not a new differentiation algorithm.
5. Framework comparison: §11 copies the exact NumPy weights and validates every gradient. Reinitializing both libraries with the same integer is not this test.

Three affine layers — a complete worked extension

Each hidden cache holds its layer input, preactivation and activation. Walk those caches backward with the forward weights frozen; only the last layer is linear. This2→3→2→1 example checks every parameter numerically, so the extension is not merely an instruction to add another layer.

def deep_relu_eval(X,y,layers):
    caches=[];a=X
    for i,(W,b) in enumerate(layers):
        z=a@W+b
        caches.append((a,z))
        a=np.maximum(0,z) if i<len(layers)-1 else z
    logits=a;y=binary_targets(y,logits)
    e=np.exp(-np.abs(logits))
    prob=np.where(logits>=0,1/(1+e),e/(1+e))
    value=float(np.mean(np.maximum(logits,0)-y*logits+np.log1p(e)))
    dz=(prob-y)/len(X);grads=[]
    for i in range(len(layers)-1,-1,-1):
        inp,z=caches[i];W,b=layers[i]
        grads.append((inp.T@dz,dz.sum(0)))
        if i:
            dz=(dz@W.T)*(caches[i-1][1]>0)
    return value,list(reversed(grads))
rng_deep=np.random.default_rng(31)
deep_layers=[(rng_deep.normal(size=(a,b))*.2,np.full(b,.4)) for a,b in [(2,3),(3,2),(2,1)]]
deep_X=np.array([[.3,-.2],[.7,.1],[-.1,.6]])
deep_y=np.array([1.,0.,1.]);_,deep_g=deep_relu_eval(deep_X,deep_y,deep_layers)
deep_errors=[]
for li,layer in enumerate(deep_layers):
    for pi,arr in enumerate(layer):
        for ix in np.ndindex(arr.shape):
            old=arr[ix];arr[ix]=old+1e-5;plus=deep_relu_eval(deep_X,deep_y,deep_layers)[0]
            arr[ix]=old-1e-5;minus=deep_relu_eval(deep_X,deep_y,deep_layers)[0];arr[ix]=old
            deep_errors.append(abs((plus-minus)/2e-5-deep_g[li][pi][ix]))
assert max(deep_errors)<1e-7
print('three-layer ReLU all-coordinate check',len(deep_errors),max(deep_errors))

Complete ReLU training and a local two-moons comparison

The protected tanh core is not a substitute for practicing the ReLU derivative. This additional program trains a2→8→1 ReLU network on XOR and a2→16→1 network on locally generated moons. It logs actual held-out results against a linear-logit baseline rather than assuming99% versus85%. The boundary SVG uses blue/red cells for predictions and includes no external plotting dependency.

def relu_train(X,y,hidden=8,lr=.5,steps=2000,seed=0):
    rng=np.random.default_rng(seed)
    theta=(rng.normal(size=(X.shape[1],hidden))*np.sqrt(2/X.shape[1]),
           np.zeros(hidden),rng.normal(size=(hidden,1))/np.sqrt(hidden),np.zeros(1))
    trace=[]
    for _ in range(steps):
        value,g,_,_=relu_eval(X,y,theta)
        trace.append(value)
        theta=tuple(a-lr*b for a,b in zip(theta,g))
    return theta,np.asarray(trace)
 
xor_X=np.array([[0.,0.],[0.,1.],[1.,0.],[1.,1.]])
xor_y=np.array([[0.],[1.],[1.],[0.]])
xor_p,xor_trace=relu_train(xor_X,xor_y)
xor_loss,_,xor_prob,_=relu_eval(xor_X,xor_y,xor_p)
assert xor_loss<xor_trace[0] and np.array_equal(xor_prob>.5,xor_y>0)
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
mx,my=make_moons(n_samples=300,noise=.2,random_state=0)
xa,xt,ya,yt=train_test_split(mx,my,test_size=.3,random_state=0,stratify=my)
moon_p,moon_trace=relu_train(xa,ya,hidden=16,lr=.3,steps=2000)
moon_prob=relu_eval(xt,yt,moon_p)[2].ravel()
linear=LogisticRegression(random_state=0).fit(xa,ya)
print('ReLU XOR fit',xor_loss,'moons held-out',((moon_prob>.5)==yt).mean(),
      'logistic held-out',linear.score(xt,yt))
assert np.isfinite(moon_trace).all()
 
def boundary_svg(X,theta,n=32):
    lo=X.min(0)-.2;hi=X.max(0)+.2
    pts=np.array([[a,b] for b in np.linspace(lo[1],hi[1],n) for a in np.linspace(lo[0],hi[0],n)])
    probs=relu_eval(pts,np.zeros(len(pts)),theta)[2].ravel()
    out=['<svg xmlns="http://www.w3.org/2000/svg" width="320" height="320" viewBox="0 0 32 32">']
    for i,pred in enumerate(probs):
        color='#93c5fd' if pred>.5 else '#fca5a5'
        out.append(f'<rect x="{i%n}" y="{n-1-i//n}" width="1" height="1" fill="{color}"/>')
    return ''.join(out)+'</svg>'
assert '<rect' in boundary_svg(mx,moon_p)

The logistic baseline uses scikit-learn 1.7's default L2 regularization (C=1, solver lbfgs), whereas this teaching MLP has none; the output is a small comparison, not an isolated proof that depth caused a gain. Both use the same fixed split. The review saves the synthetic boundary and loss values as evidence; no MNIST data is involved.

15 · Recap and retention answers

1. Linear rules: dX=G@W.T, dW=X.T@G, db=G.sum(0); coordinate sums explain each axis.
2. Activation rules: ReLU G*(Z>0) with a zero convention; sigmoid G*p*(1-p); fused mean BCE (p-y)/N.
3. Nine-parameter trace: dW1 is [[-.1104869461,-.1325843353],[-.2209738922,-.2651686707]].
4. Cache: dW2 needs H, dW1 needs X, and the activation derivative needs its forward state; recomputation is an alternative.
5. Shape rule: gradient and differentiated tensor shapes agree; still test semantics and factors independently.
6. Sum loss: remove the mean denominator; repeated identical examples multiply the sum gradient by N.
7. Three layers: repeat the §14 chain at frozen parameters and numerically check all coordinates, not only a selected matrix.
8. Attention backward: §12 contains the matmul pullbacks and softmax row-centering term. FA2 changes the execution schedule and memory traffic, not those derivatives.

Tomorrow reconstruct the nine numbers and the flow; next week reproduce the protected XOR run and explain why fit is not generalization. MNIST comes next as a local-data recipe, not a promised accuracy deadline.

Sources and execution boundary

PyTorch 2.9 autograd mechanics, BCEWithLogitsLoss, checkpoint API, Chen et al. 2016, FlashAttention-2, and Baydin et al.'s AD survey. Exact passages and raw hashes are in the review artifacts. CPU NumPy/PyTorch experiments and independent derivative checks are separate from source support; no FlashAttention kernel, GPU or large-model run is claimed.