The Chain Rule, Computational Graphs, and Automatic Differentiation

Which mode is cheaper depends on how many inputs and how many outputs you need derivatives for. Forward mode carries input-direction derivatives alongside the values. It computes Jacobian–vector products, is natural for directional sensitivity, and costs one sweep per chosen input direction. That makes it efficient for few inputs and many outputs. Reverse mode pushes output sensitivities backward instead. It computes vector–Jacobian products and is efficient for many inputs and few outputs — the shape of a scalar training loss. It also needs access to the intermediate values. Finite differences do neither: they perturb the inputs and run the program again. They are approximate and step-size sensitive, and the cost grows with the input dimension. They also work through some black-box code, which makes them a useful independent check.

The two modes have a published price. For a map f:Rn→Rm, the m×n Jacobian costs n·c·ops(f) in forward mode and m·c·ops(f) in reverse mode. The overhead constant is small: “where c is a constant guaranteed to be c < 6 and typically c ∼ [2, 3] (Griewank and Walther, 2008)”. That is from the 2018 survey of automatic differentiation by Baydin and three colleagues. The two formulas differ in one symbol. A training loss has m = 1. That is why a single reverse sweep returns the whole gradient, however many parameters n contains.

The gradient case has an exact bound. Griewank's 2012 history of the reverse mode gives it as OPS{∇f(x)} ≤ ω·OPS{f(x)}, independent of n, with ω = 3 when only polynomial operations are counted by multiplications. His example is the determinant of a dense symmetric positive matrix computed by Cholesky, whose gradient is the adjugate. The guarantee stops at gradients. For F(x) = b·sin(aᵀx) the function costs about n+m multiplications, while F′(x) costs n·m, because no method can produce an n×m array of distinct numbers cheaply. Griewank draws the conclusion in one line: “Hence, the cheap gradient principle does not extend to a cheap Jacobian principle.”

Keeping the intermediate values is not a footnote either. It is the mode's bill, and it has been measured. The survey records that “The advantages of reverse mode AD, however, come with the cost of increased storage requirements growing (in the worst case) in proportion to the number of operations in the evaluated function”, and names checkpointing as the remedy. Griewank derives MEM{F′(x)} ∼ OPS{F(x)} and states that checkpointing buys “a logarithmic increase for both spatial and temporal complexity”. Chen and three colleagues put numbers on the trade in April 2016: O(√n) memory for an n-layer network at the cost of one extra forward pass per mini-batch, and O(log n) memory at O(n log n) extra forward computation in the extreme. Their abstract reports the measured case: “Our experiments show that we can reduce the memory cost of a 1,000-layer deep residual network from 48G to 7G with only 30 percent additional running time cost on ImageNet problems.” Reverse mode is not free. It is cheap in time and expensive in memory, and the exchange rate is known.

A neural network is a composition of simple operations: matrix products, additions, normalizations, activations, and reductions. The chain rule tells us how local derivatives combine. A computational graph records values and dependencies. Automatic differentiation traverses that graph while propagating derivative information according to the program's executed operations. Backpropagation is reverse-mode automatic differentiation specialized to scalar objectives and layered computations. It is not a separate law of calculus. The record of who invented what, when, and for which purpose is the shortest way to show it.

Neither mode was invented for neural networks. The 2018 survey dates them: forward mode “was essentially discovered by Wengert (1964)”, and “Prior to Werbos, the work by Linnainmaa (1970, 1976) is often cited as the first published description of the reverse mode.” The 1976 item is a real, findable paper about something else entirely. Seppo Linnainmaa's paper on the Taylor expansion of the accumulated rounding error appeared in BIT, received 22 January 1976, and it propagates coefficients backward through a sequence of elementary operations u_i = Q_i(u_j, u_k). Griewank's history opens with where the idea came from: “Seppo Linnainmaa (Lin76) of Helsinki says the idea came to him on a sunny afternoon in a Copenhagen park in 1970. He used it as a tool for estimating the effects of arithmetic rounding errors on the results of complex expressions.” The machinery that trains modern networks was built to track floating-point error. That was six years before the paper, and sixteen before backpropagation reached Nature.

The famous paper is four pages long. It appeared in Nature on 9 October 1986, by Rumelhart and two co-authors, describing “a new learning procedure, back-propagation, for networks of neurone-like units” that “repeatedly adjusts the weights of the connections in the network so as to minimize a measure of the difference between the actual output vector of the net and the desired output vector”. Schmidhuber's 2015 survey of deep learning — 88 pages, 888 references — places that letter as popularisation rather than invention. It “significantly contributed to the popularisation of BP for NNs”, while the first description belongs to Linnainmaa's 1970 master's thesis, “albeit without reference to NNs”. The same survey states the identification directly: “BP is also known as the reverse mode of automatic differentiation (Griewank, 2012), where the costs of forward activation spreading essentially equal the costs of backward derivative calculation.” Same algorithm. Different decade, different field, different purpose.

Let u=x² and v=3x, then L=u+v.

For L=(wx−y)², the forward pass records intermediate values and the reverse pass accumulates sensitivities. First, multiply: a = wx, recording x and w for the local derivatives. Second, subtract: e = a−y, whose local derivatives are 1 and −1. Third, square: L = e², whose local derivative is 2e. Fourth, reverse: propagate ∂L/∂e back to a, then to w and x. Fifth, accumulate: if a value influences several paths, add all incoming contributions.

Every stored intermediate in step one is the memory bill from the comparison above, on the smallest possible scale — a and e here, 48G of activations in a 1,000-layer residual network. It is the same trade. The forward values are kept because the backward sweep needs them, and checkpointing recomputes them instead when keeping them is too expensive.

Follow the gradient path rather than inspecting only the final norm. First, confirm participation: check that every intended parameter influences the loss and requires gradients. Second, inspect local scales: log activation, Jacobian, and gradient magnitudes at key boundaries. The λ₁ < 1/γ condition tells you what you are looking for, and for sigmoid units the number to compare against is 4. Third, check branches: verify accumulation, masking, detachment, and conditional execution. That includes any straight-through construction of the y_hard - y_soft.detach() + y_soft form, which is designed to make the forward and backward graphs disagree. Fourth, compare a small case. Fifth, validate the surrogate: ask whether the differentiable objective represents the actual discrete or operational goal.

Step four has tools with published defaults, so it need not be improvised. PyTorch's torch.autograd.gradcheck runs with eps=1e-06, atol=1e-05, rtol=0.001 and nondet_tol=0.0. Its documentation attaches a condition to those numbers: “The default values are designed for input of double precision. This check will likely fail if input is of less precision, e.g., FloatTensor.” The same page keeps a failed check from being read as a proven bug: “Gradcheck may fail when evaluated on non-differentiable points because the numerically computed gradients via finite differencing may differ those computed analytically (not necessarily because either is incorrect)”. JAX offers jax.test_util.check_grads(f, args, order, modes=('fwd','rev'), atol=None, rtol=None, eps=None) to “Check gradients from automatic differentiation against finite differences.” It solves the dimension problem by sampling: “Gradients are only checked in a single randomly chosen direction, which ensures that the finite difference calculation does not become prohibitively expensive even for large input/output spaces.” A cross-check in double precision, in one random direction, against a documented tolerance, is a stronger report than a gradient norm that looks plausible.

If z=g(x) and y=f(z), then the derivative of y with respect to x composes the derivative maps. In Jacobian notation, Jᵧₓ = Jᵧ_z J_zₓ under a common convention. Reverse-mode systems often work with row or column cotangents and compute vector–Jacobian products without forming the full Jacobian. Frameworks differ in how they do it, but the invariant idea is composition of local linear maps.

Not forming the full Jacobian is the entire economy of the method, not an implementation detail. Griewank's counterexample makes the point: F(x) = b·sin(aᵀx) costs about n+m multiplications, while the derivative array costs n·m. A system that insisted on materialising Jᵧₓ would forfeit the cheap-gradient bound OPS{∇f(x)} ≤ ω·OPS{f(x)} with ω = 3. Forward mode carries a Jacobian–vector product. Reverse mode carries a vector–Jacobian product. In both cases what moves through the graph is a vector the size of the data, never the matrix the notation names.