44.03.08 · optimization-control / unconstrained-optimization

Derivative Computation: Finite Differences and Automatic Differentiation

shipped3 tiersLean: none

Anchor (Master): Griewank & Walther 2008 Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation 2e (SIAM) Ch. 3-4 (forward and reverse accumulation), Ch. 12 (reversal schedules and checkpointing); Nocedal & Wright 2006 Numerical Optimization 2e (Springer) Ch. 8 (finite-difference and automatic derivatives for optimization); Baydin, Pearlmutter, Radul & Siskind 2018 Automatic Differentiation in Machine Learning: a Survey (JMLR 18) for the ML-training context

Intuition Beginner

Almost every fast optimizer needs to know which way is downhill, and how the slope is bending. That information is the derivative: the gradient (the vector of slopes in each direction) and sometimes the second derivatives (how the slopes change). The trouble is that you often have a function only as a piece of code — a long recipe of additions, multiplications, square roots, and so on — with no tidy formula for its slope. So the practical question is: given the code for $f$ , how do you get the code for its derivative?

The oldest trick is to nudge and compare. Move the input a tiny amount $h$ , see how much the output changes, and divide: the change in output over the change in input is an estimate of the slope. This is the finite-difference idea, and it needs nothing but the ability to run $f$ . Its weakness is the size of the nudge. Too big a nudge and the straight-line estimate is crude; too small a nudge and the two outputs are so close that rounding noise in the computer swamps the tiny difference. There is a sweet spot in the middle, and it is not as accurate as you might hope.

The sharper trick is to differentiate the recipe itself. Each basic step — add, multiply, take a sine — has a known derivative rule, and the chain rule stitches those rules together. A computer can apply the chain rule mechanically as it runs the code, producing the exact derivative (up to ordinary rounding), with no nudge and no guessing of $h$ . This is automatic differentiation, and it comes in two directions of travel: forward, alongside the computation, and backward, retracing it.

The backward direction holds the prize that powers modern machine learning: it returns the slope in every input direction at once, for a cost close to a single run of $f$ — no matter whether there are ten inputs or ten million.

Visual Beginner

Picture the finite-difference estimate first. On the curve $y = f (x)$ , mark your point and a second point a little to the right at distance $h$ . Draw the straight line through the two marks: its steepness is the finite-difference slope, a chord standing in for the true tangent. A symmetric version uses one mark on each side, left and right, and the chord through those two hugs the tangent more closely.

The step-size sweet spot is the heart of the method. Shrink $h$ and the chord lines up better with the tangent — until $h$ gets so small that the two output values agree in almost all their digits, and the subtraction throws away accuracy. The error first falls, bottoms out, then climbs again.

method	what it needs	cost of a full gradient ( $n$ inputs)	accuracy
forward difference	run $f$ only	about $n + 1$ runs of $f$	modest (about half the digits)
central difference	run $f$ only	about $2 n$ runs of $f$	better (about two-thirds of the digits)
forward-mode AD	the code of $f$	about $n$ runs' worth	full (exact rule, only rounding)
reverse-mode AD	the code of $f$	a few runs' worth, any $n$	full (exact rule, only rounding)

The last row is the surprise: the cost of the whole gradient does not grow with the number of inputs. That single fact is why training a network with millions of parameters is even possible.

Worked example Beginner

Take $f (x) = x^{3}$ at the point $x = 2$ , where the true slope is $f^{'} (2) = 3 \cdot 2^{2} = 12$ . We estimate it three ways.

Forward difference with $h = 0.1$ . Compute $f (2.1) = 9.261$ and $f (2) = 8$ . The chord slope is

\frac{f ( 2.1 ) - f ( 2 )}{0.1} = \frac{9.261 - 8}{0.1} = \frac{1.261}{0.1} = 12.61.

The error is $12.61 - 12 = 0.61$ — too high, because the chord leans up over the interval.

Central difference with $h = 0.1$ . Compute $f (2.1) = 9.261$ and $f (1.9) = 6.859$ . The symmetric chord slope is

\frac{f ( 2.1 ) - f ( 1.9 )}{2 \cdot 0.1} = \frac{9.261 - 6.859}{0.2} = \frac{2.402}{0.2} = 12.01.

The error is only $0.01$ — the two-sided chord cancels most of the lean.

Automatic differentiation. The code for $f$ is "multiply $x$ by $x$ to get $x^{2}$ , then multiply by $x$ again." Differentiating step by step with the product rule and feeding the slope through gives $3 x^{2} = 12$ at $x = 2$ , exactly, with no $h$ at all.

What this tells us: halving the lean by going symmetric turned an error of $0.61$ into $0.01$ , yet both finite-difference numbers are still approximations tied to the choice of $h$ . The automatic-differentiation answer is the exact $12$ , because it differentiates the recipe rather than sampling the output. The next sections make the step-size sweet spot precise and explain how the backward sweep gets every slope for nearly free.

Check your understanding Beginner

Exercise (easy, multiple-choice).

What makes reverse-mode automatic differentiation special for problems with very many inputs?

It needs a separate run of $f$ for each input.
It returns the slope in every input direction at once, for a cost close to a single run of $f$ .
It avoids the chain rule entirely.
It only works when there is exactly one input.

Hint

Think about the cost row in the table that does not grow with the number of inputs.

Answer

B. A single backward sweep produces the whole gradient — every input's slope — at a cost of a few function evaluations, independent of how many inputs there are. Feedback-correct: this is the cheap-gradient property that powers training with millions of parameters. Feedback-wrong: A describes finite differences and forward mode; C is false since reverse mode is the chain rule run backward; D misreads the method.

Formal definition Intermediate+

Let $f : R^{n} \to R$ be smooth, evaluated in floating-point arithmetic with unit roundoff $u$ (the floating-point model of 43.01.01). Write $e_{i}$ for the $i$ -th coordinate vector and $h > 0$ for a step.

Finite differences. The forward-difference approximation to a partial derivative is

D_{h}^{+} f (x)_{i} = \frac{f ( x + h e _{i} ) - f ( x )}{h},

and the central-difference approximation is

D_{h}^{0} f (x)_{i} = \frac{f ( x + h e _{i} ) - f ( x - h e _{i} )}{2 h} .

By Taylor's theorem 02.05.05, $D_{h}^{+} f (x)_{i} = \partial_{i} f (x) + \frac{1}{2} h \partial_{i}^{2} f (ξ) = \partial_{i} f (x) + O (h)$ , while the odd-order terms cancel in the symmetric formula, giving $D_{h}^{0} f (x)_{i} = \partial_{i} f (x) + \frac{1}{6} h^{2} \partial_{i}^{3} f (ξ) = \partial_{i} f (x) + O (h^{2})$ . A Hessian-vector product is obtained by differencing the gradient,

\nabla^{2} f (x) v \approx \frac{\nabla f ( x + h v ) - \nabla f ( x )}{h},

so a Newton-Krylov solver can apply $\nabla^{2} f$ to a vector using only gradient evaluations, never forming the Hessian. A full forward-difference gradient costs $n + 1$ evaluations of $f$ ; a central-difference gradient costs $2 n$ . The symbols $u$ (unit roundoff), $\partial_{i}$ (partial derivative), $\nabla$ (gradient), $\nabla^{2}$ (Hessian), $J$ (Jacobian), $\overset{v}{˙}$ (forward tangent), and $\overset{v}{ˉ}$ (reverse adjoint) are recorded in _meta/NOTATION.md.

The evaluation procedure. Automatic differentiation acts not on a formula but on an evaluation procedure (a Wengert list): the function is a finite straight-line program of intermediate variables $v_{1}, \dots, v_{L}$ , each $v_{k} = φ_{k} (v_{j} : j ≺ k)$ an elementary operation ( $+, \times, sin, exp, \dots$ ) of earlier variables, with the inputs $v_{1}, \dots, v_{n} = x_{1}, \dots, x_{n}$ and the output $y = v_{L}$ . The partial order $≺$ is the computational graph: a directed acyclic graph with the inputs as sources and $y$ as the sink.

Forward (tangent) mode. Fix a direction $\overset{x}{˙} \in R^{n}$ . Carry with each $v_{k}$ a tangent $\overset{v}{˙}_{k}$ and propagate it forward by the chain rule alongside the primal values:

\overset{v}{˙}_{k} = j ≺ k \sum \frac{\partial φ _{k}}{\partial v _{j}} \overset{v}{˙}_{j} .

Equivalently, replace each value $v$ by the dual number $v + \overset{v}{˙} ε$ with $ε^{2} = 0$ ; running the unchanged program on duals computes $y + (J \overset{x}{˙}) ε$ , so one forward sweep yields the Jacobian-vector product $J \overset{x}{˙}$ (a single directional derivative). For a scalar $f$ this is one component of the gradient per sweep, hence $n$ sweeps for the full gradient.

Reverse (adjoint) mode. Run the primal forward once, recording the graph (the tape). Then propagate adjoints $\overset{v}{ˉ}_{k} = \partial y / \partial v_{k}$ backward from $\overset{y}{ˉ} = 1$ :

\overset{v}{ˉ}_{j} = k : j ≺ k \sum \overset{v}{ˉ}_{k} \frac{\partial φ _{k}}{\partial v _{j}},

accumulating into each node from its successors. After the backward sweep, $\overset{v}{ˉ}_{1}, \dots, \overset{v}{ˉ}_{n}$ are exactly the components of $\nabla f (x)$ . One forward run plus one backward run yields the entire gradient — the vector-Jacobian product $\overset{y}{ˉ}^{⊤} J$ — independent of $n$ . Backpropagation is precisely reverse-mode AD applied to a scalar loss.

Counterexamples to common slips Intermediate+

"AD is just finite differences done carefully." No nudge is taken and no $h$ is chosen; AD applies exact differentiation rules to each elementary operation, so its only error is ordinary floating-point rounding, not truncation. The step-size dilemma below does not apply to AD at all.
"AD is symbolic differentiation." Symbolic differentiation manipulates expression trees and suffers expression swell: the derivative of a product of $m$ factors expands to $m$ terms, and nested compositions blow up combinatorially. AD instead evaluates derivatives numerically at the given point while reusing every shared subexpression through the graph, keeping the work proportional to the original program.
"Forward and reverse mode cost the same." For $f : R^{n} \to R^{m}$ , forward mode costs $O (n)$ sweeps (one per input direction) and reverse mode costs $O (m)$ sweeps (one per output). For a scalar loss ( $m = 1$ and $n$ large) reverse mode wins decisively; for $m ≫ n$ (a tall Jacobian) forward mode is cheaper.
"Reverse mode is free." It is cheap in time but pays in memory: the tape must store the intermediate values needed by the backward sweep, a cost proportional to the number of operations. Checkpointing trades some of this memory back for recomputation.

Key theorem with proof Intermediate+

The signature result is the error decomposition that fixes the optimal finite-difference step, because it is the single statement that quantifies the truncation-versus-roundoff dilemma and pins the best attainable accuracy to the machine precision $u$ of 43.01.01.

Theorem (optimal finite-difference step and accuracy floor). Let $f : R \to R$ be $C^{2}$ near $x$ with $∣ f^{''} ∣ \leq M_{2}$ on the relevant interval, and suppose each evaluation of $f$ is computed with a relative error at most $u$ , so the stored value satisfies $∣ fl (f (x)) - f (x) ∣ \leq u ∣ f (x) ∣ =: δ$ . Then the computed forward difference $D_{h}^{+} f (x)$ satisfies

D_{h}^{+} f (x) - f^{'} (x) \leq truncation \frac{1}{2} M_{2} h + roundoff \frac{2 δ}{h},

and the right-hand side is minimised at

h^{⋆} = 2 δ / M_{2} \sim u, with minimal bound 2 M_{2} δ \sim u .

The central difference, with $∣ f^{'''} ∣ \leq M_{3}$ , obeys $∣ D_{h}^{0} f (x) - f^{'} (x) ∣ \leq \frac{1}{6} M_{3} h^{2} + δ / h$ , minimised at $h^{⋆} \sim u^{1/3}$ with accuracy floor $\sim u^{2/3}$ .

Proof. Taylor's theorem 02.05.05 with Lagrange remainder gives $f (x + h) = f (x) + h f^{'} (x) + \frac{1}{2} h^{2} f^{''} (ξ)$ for some $ξ$ between $x$ and $x + h$ , so the exact forward difference satisfies

\frac{f ( x + h ) - f ( x )}{h} - f^{'} (x) = \frac{1}{2} h ∣ f^{''} (ξ) ∣ \leq \frac{1}{2} M_{2} h,

the truncation term. The computed numerator replaces $f (x + h), f (x)$ by stored values differing from the exact ones by at most $δ$ each, so the computed quotient differs from the exact quotient by at most $(δ + δ) / h = 2 δ / h$ , the roundoff term. The triangle inequality adds the two bounds. Let $g (h) = \frac{1}{2} M_{2} h + 2 δ / h$ ; then $g^{'} (h) = \frac{1}{2} M_{2} - 2 δ / h^{2}$ vanishes at $h^{⋆} = 2 δ / M_{2}$ , where $g^{''} (h^{⋆}) = 4 δ / (h^{⋆})^{3} > 0$ , so $h^{⋆}$ is the minimiser. Substituting, $g (h^{⋆}) = \frac{1}{2} M_{2} \cdot 2 δ / M_{2} + 2 δ / (2 δ / M_{2}) = 2 M_{2} δ$ . Since $δ = u ∣ f (x) ∣$ , both $h^{⋆}$ and the floor scale as $u$ up to the $O (1)$ factor $∣ f (x) ∣/ M_{2}$ . For the central difference the truncation term is $\frac{1}{6} M_{3} h^{2}$ (the odd-order Taylor terms cancel) and the roundoff term is $2 δ / (2 h) = δ / h$ ; minimising $\frac{1}{6} M_{3} h^{2} + δ / h$ gives $h^{⋆} = (3 δ / M_{3})^{1/3} \sim u^{1/3}$ and floor $\sim u^{2/3}$ . $□$

Bridge. This decomposition is the foundational reason finite differences can never reach full machine accuracy: the subtraction of two nearby values forfeits about half the available digits for a forward difference and a third for a central difference, so the best a gradient-by-differencing can do is $u$ or $u^{2/3}$ , and this is exactly the cancellation phenomenon the floating-point model of 43.01.01 predicts. The $u$ rule of thumb builds toward the contrast with automatic differentiation, where no subtraction of nearby quantities occurs and the only error is the ordinary rounding of 43.01.01, so AD returns derivatives to nearly full precision. The same Taylor cancellation that costs accuracy is dual to the cancellation that buys the higher order of the central formula: killing the $O (h)$ term is the central insight that improves the floor from $u$ to $u^{2/3}$ . Putting these together, the step-size analysis generalises directly to Hessian-vector products formed by differencing the gradient, where the accuracy floor degrades by another factor, and the dilemma appears again in 44.03.04 wherever a quasi-Newton or Newton-Krylov method must choose between differenced and automatically differentiated curvature information; the bridge is that AD removes the dilemma entirely by differentiating the program rather than sampling its output.

Exercises Intermediate+

Exercise 4 (medium, symbolic).

Show that running an evaluation procedure on dual numbers $v + \overset{v}{˙} ε$ with $ε^{2} = 0$ computes the directional derivative. Verify the product rule: if $a = uv$ with duals $u + \overset{u}{˙} ε$ , $v + \overset{v}{˙} ε$ , then the $ε$ -coefficient of the product is $\overset{u}{˙} v + u \overset{v}{˙}$ .

Hint

Multiply $(u + \overset{u}{˙} ε) (v + \overset{v}{˙} ε)$ and discard the $ε^{2}$ term.

Answer

Expanding, $(u + \overset{u}{˙} ε) (v + \overset{v}{˙} ε) = uv + (\overset{u}{˙} v + u \overset{v}{˙}) ε + \overset{u}{˙} \overset{v}{˙} ε^{2}$ . Since $ε^{2} = 0$ , this is $uv + (\overset{u}{˙} v + u \overset{v}{˙}) ε$ , whose real part is the primal product and whose $ε$ -coefficient is exactly the product-rule derivative $\overset{u}{˙} v + u \overset{v}{˙}$ . By the same mechanism each elementary operation $φ$ maps $v + \overset{v}{˙} ε \mapsto φ (v) + φ^{'} (v) \overset{v}{˙} ε$ , so composing the whole program sends $x + \overset{x}{˙} ε$ to $f (x) + (J \overset{x}{˙}) ε$ : the $ε$ -coefficient is the Jacobian-vector product, the forward-mode directional derivative. Full credit for the expansion and the closure of the rule under composition.

Exercise 5 (medium, symbolic).

For $f (x_{1}, x_{2}) = x_{1} x_{2} + sin x_{1}$ , trace reverse mode. Forward: $v_{1} = x_{1}$ , $v_{2} = x_{2}$ , $v_{3} = v_{1} v_{2}$ , $v_{4} = sin v_{1}$ , $y = v_{5} = v_{3} + v_{4}$ . Starting from $\overset{y}{ˉ} = 1$ , compute the adjoints and confirm $\overset{v}{ˉ}_{1} = \partial f / \partial x_{1}$ , $\overset{v}{ˉ}_{2} = \partial f / \partial x_{2}$ .

Hint

$\overset{v}{ˉ}_{3} = \overset{v}{ˉ}_{4} = \overset{v}{ˉ}_{5} = 1$ from the sum. Then $v_{1}$ feeds both $v_{3}$ and $v_{4}$ , so accumulate $\overset{v}{ˉ}_{1} = \overset{v}{ˉ}_{3} \cdot v_{2} + \overset{v}{ˉ}_{4} \cdot cos v_{1}$ .

Answer

From $v_{5} = v_{3} + v_{4}$ : $\overset{v}{ˉ}_{3} = \overset{v}{ˉ}_{4} = \overset{v}{ˉ}_{5} = 1$ . From $v_{4} = sin v_{1}$ : contribute $\overset{v}{ˉ}_{4} cos v_{1} = cos x_{1}$ to $\overset{v}{ˉ}_{1}$ . From $v_{3} = v_{1} v_{2}$ : contribute $\overset{v}{ˉ}_{3} v_{2} = x_{2}$ to $\overset{v}{ˉ}_{1}$ and $\overset{v}{ˉ}_{3} v_{1} = x_{1}$ to $\overset{v}{ˉ}_{2}$ . Accumulating, $\overset{v}{ˉ}_{1} = x_{2} + cos x_{1}$ and $\overset{v}{ˉ}_{2} = x_{1}$ . These are exactly $\partial f / \partial x_{1} = x_{2} + cos x_{1}$ and $\partial f / \partial x_{2} = x_{1}$ : a single backward sweep produced the whole gradient. The node $v_{1}$ correctly summed contributions from both of its successors $v_{3}$ and $v_{4}$ . Full credit for both adjoints and the accumulation at the fan-out node.

Exercise 6 (hard, symbolic).

Let $f : R^{n} \to R^{m}$ be evaluated by a procedure of $OPS$ elementary operations. Give the per-evaluation cost of computing the full Jacobian $J$ by (a) forward mode, (b) reverse mode, (c) forward differences, in terms of $n$ , $m$ , and $OPS$ , and state the regime in which each is preferred.

Hint

Forward mode gets one Jacobian column (one $J \overset{x}{˙}$ ) per sweep; reverse mode gets one Jacobian row (one $\overset{y}{ˉ}^{⊤} J$ ) per sweep; a forward difference gets one column per extra function evaluation.

Answer

(a) Forward mode propagates one tangent per input direction, so the full $J$ (its $n$ columns) costs $O (n \cdot OPS)$ , with a small constant ( $\leq 2.5$ ) per sweep. (b) Reverse mode propagates one adjoint per output, so $J$ (its $m$ rows) costs $O (m \cdot OPS)$ , with constant $\leq 4$ per sweep but additional $O (OPS)$ memory for the tape. (c) Forward differences need one extra evaluation of $f$ per column, so $O (n \cdot OPS)$ in time like forward mode, but with truncation/roundoff error capped at $u$ instead of full precision. Preference: reverse mode when $m ≪ n$ (a scalar loss, $m = 1$ , is the extreme case — the cheap-gradient regime); forward mode when $n ≪ m$ (a tall Jacobian); finite differences only when the source code is unavailable. Full credit for the three cost orders and the $m$ versus $n$ crossover.

Exercise 7 (hard, symbolic).

Derive the Hessian-vector product accuracy when $\nabla^{2} f (x) v$ is formed by differencing an automatically differentiated gradient: $\nabla^{2} f (x) v \approx (\nabla f (x + h v) - \nabla f (x)) / h$ , with $\nabla f$ computed to relative error $u$ by reverse mode. Find the optimal $h$ and the floor, and contrast with applying forward-over-reverse AD instead.

Hint

Treat $\nabla f$ as a vector-valued function known to accuracy $δ \sim u ∥\nabla f ∥$ ; reuse the forward-difference balance $\frac{1}{2} M h + 2 δ / h$ componentwise.

Answer

Componentwise the differenced gradient has truncation error $\frac{1}{2} h ∥ \nabla^{3} f ∥ ∥ v ∥^{2}$ and roundoff error $\sim 2 δ / h$ with $δ \sim u ∥\nabla f ∥$ , so by the Key theorem's balance the optimum is $h^{⋆} \sim u$ and the accuracy floor is $\sim u$ relative — only half the digits, even though $\nabla f$ itself was nearly exact, because the outer difference still subtracts nearby vectors. Applying forward-over-reverse AD instead — pushing a dual through the reverse sweep, i.e. evaluating the directional derivative of $\nabla f$ in direction $v$ — computes $\nabla^{2} f (x) v$ to full precision with no step size, at cost a small constant times one gradient evaluation. So differencing an exact gradient still loses half the digits; only nesting AD removes the loss. Full credit for $h^{⋆} \sim u$ , floor $\sim u$ , and the forward-over-reverse alternative.

Advanced results Master

The finite-difference floor is set by cancellation and the optimal step by a one-dimensional minimisation; the results below place the exact-derivative alternative on its complexity-theoretic footing — the cheap-gradient principle, the duality of the two AD modes, the memory cost of reversal and its checkpointing remedy, and the position of AD between the symbolic and numerical extremes.

Theorem 1 (the cheap-gradient principle). Let $f : R^{n} \to R$ be computed by an evaluation procedure whose runtime, measured in a fixed cost model counting additions, multiplications, and library-call evaluations, is $TIME (f)$ . Reverse-mode AD computes the complete gradient $\nabla f (x)$ in time

TIME (\nabla f) \leq ω_{grad} \cdot TIME (f), ω_{grad} \in [3, 4],

independent of $n$ . The bound holds because each elementary operation $v_{k} = φ_{k} (v_{j})$ contributes to the backward sweep only its local partials $\partial φ_{k} / \partial v_{j}$ times the incoming adjoint, a constant multiple of the work already done in the forward sweep; summing over all $L$ operations gives a constant overhead factor. The contrast with the $n + 1$ evaluations of a forward-difference gradient is the entire reason gradient-based training of models with $n \sim 1 0^{9}$ parameters is feasible ^{[Griewank, A. & Walther, A. — Evaluating Derivatives]}.

Theorem 2 (forward-reverse duality and the Jacobian). For $f : R^{n} \to R^{m}$ with Jacobian $J \in R^{m \times n}$ , forward mode evaluates the matrix-free product $v \mapsto J v$ (a tangent) and reverse mode evaluates $w \mapsto J^{⊤} w$ (an adjoint/cotangent), each at cost a small constant times $TIME (f)$ . The full Jacobian costs $min (n, m)$ such sweeps up to the optimal-Jacobian-accumulation problem — finding the cheapest order in which to multiply the elementary local Jacobians along the graph — which is NP-hard in general, so practical AD uses the forward or reverse extreme or hand-tuned cross-country schedules. Forward and reverse are formal transposes: reverse mode on $f$ is forward mode on the adjoint program, the algorithmic shadow of $(J^{⊤})^{⊤} = J$ ^{[Baydin, A. G., Pearlmutter, B. A., Radul, A. A. & Siskind, J. M. — Automatic Differentiation in Machine Learning: a Survey]}.

Theorem 3 (memory cost of reversal and logarithmic checkpointing). Reverse mode must make the intermediate values of the forward sweep available to the backward sweep. The naive tape stores all $L$ intermediates, $SPACE = O (L)$ . For an evaluation that decomposes into a chain of $L$ time steps (a deep network, a time-stepping PDE solver), the Revolve checkpointing schedule stores only $s$ checkpoints and recomputes intermediate segments on demand, achieving the optimal tradeoff: with $s$ checkpoints the chain of length $L$ is reversed using at most $t$ extra forward recomputations where $(s s + t) \geq L$ . Choosing $s = t = O (lo g L)$ gives memory $O (lo g L)$ at the price of an $O (lo g L)$ factor in time — the binomial reversal schedule, optimal among all schedules ^{[Griewank, A. & Walther, A. — Evaluating Derivatives]}.

Theorem 4 (AD versus symbolic and numerical differentiation). Symbolic differentiation returns a closed-form expression but suffers expression swell: the symbolic derivative of a length- $L$ program can have size exponential in $L$ when common subexpressions are not shared, and evaluating it then costs far more than the cheap-gradient bound. Numerical (finite) differentiation is cheap to implement but caps accuracy at $u$ (forward) or $u^{2/3}$ (central) by the Key theorem and costs $n + 1$ evaluations for a gradient. Automatic differentiation occupies the remaining corner: it evaluates derivatives to full machine precision (no truncation), reuses every shared subexpression (no swell), and computes a scalar gradient at $O (1)$ function evaluations (no $n$ -fold blowup). The three methods are distinguished not by the chain rule — all three obey it — but by what they carry: symbolic carries expressions, numerical carries perturbed values, AD carries derivative values attached to the live computation ^{[Baydin, A. G., Pearlmutter, B. A., Radul, A. A. & Siskind, J. M. — Automatic Differentiation in Machine Learning: a Survey]}.

Synthesis. Derivative computation splits into two regimes, and the foundational reason they differ is whether one samples the function or differentiates its program. Finite differences sample: the Key theorem shows this forfeits half the digits to the cancellation of 43.01.01, the $u$ floor an irreducible tax on subtracting nearby values, and the cost scales with the input dimension $n$ . Automatic differentiation differentiates the program: the cheap-gradient principle is exactly the statement that the reverse sweep reuses the forward computation operation for operation, so a scalar gradient costs a fixed multiple of the function regardless of $n$ — this is the central insight that makes large-scale optimization and the training of billion-parameter models possible. The two AD modes are dual, forward computing $J v$ and reverse computing $J^{⊤} w$ , and the choice between them is governed by the shape of the Jacobian just as the choice of finite-difference order is governed by the balance of Taylor terms; backpropagation is simply reverse mode at a scalar loss.

The memory cost of reversal generalises the time-space tradeoff that pervades scientific computing, and logarithmic checkpointing is its optimal resolution, putting these together so that even an arbitrarily deep computation can be differentiated in bounded memory. The bridge from this unit to the rest of optimization is that every Newton, quasi-Newton, and stochastic-gradient method downstream consumes the gradients and Hessian-vector products produced here, and their efficiency is inherited from the cheap-gradient bound — AD is dual to the optimizer in the sense that the optimizer chooses the search direction and AD supplies, almost for free, the curvature information that direction needs.

Full proof set Master

Proposition 1 (forward-difference error decomposition and optimal step). Under the Key theorem's hypotheses, the computed forward difference obeys $∣ D_{h}^{+} f (x) - f^{'} (x) ∣ \leq \frac{1}{2} M_{2} h + 2 δ / h$ with $δ = u ∣ f (x) ∣$ , minimised at $h^{⋆} = 2 δ / M_{2} \sim u$ with floor $2 M_{2} δ \sim u$ .

Proof. This is the Key theorem. Taylor with Lagrange remainder 02.05.05 bounds the exact-quotient truncation by $\frac{1}{2} M_{2} h$ ; the floating-point model of 43.01.01 bounds each stored value's error by $δ$ , contributing $2 δ / h$ after division; the triangle inequality sums them. Calculus on $g (h) = \frac{1}{2} M_{2} h + 2 δ / h$ gives the unique positive critical point $h^{⋆} = 2 δ / M_{2}$ with $g^{''} > 0$ , and $g (h^{⋆}) = 2 M_{2} δ$ . $□$

Proposition 2 (central-difference order and floor). For $f \in C^{3}$ with $∣ f^{'''} ∣ \leq M_{3}$ , $∣ D_{h}^{0} f (x) - f^{'} (x) ∣ \leq \frac{1}{6} M_{3} h^{2} + δ / h$ , minimised at $h^{⋆} = (3 δ / M_{3})^{1/3} \sim u^{1/3}$ with floor $\sim u^{2/3}$ .

Proof. Taylor to third order at $x \pm h$ gives $f (x \pm h) = f (x) \pm h f^{'} (x) + \frac{1}{2} h^{2} f^{''} (x) \pm \frac{1}{6} h^{3} f^{'''} (ξ_{\pm})$ ; subtracting cancels the even-order terms, so $\frac{f ( x + h ) - f ( x - h )}{2 h} - f^{'} (x) = \frac{1}{12} h^{2} (f^{'''} (ξ_{+}) + f^{'''} (ξ_{-}))$ , bounded by $\frac{1}{6} M_{3} h^{2}$ . The two stored values each carry error $δ$ , contributing $2 δ / (2 h) = δ / h$ . Minimising $\frac{1}{6} M_{3} h^{2} + δ / h$ : the derivative $\frac{1}{3} M_{3} h - δ / h^{2} = 0$ gives $h^{3} = 3 δ / M_{3}$ , so $h^{⋆} = (3 δ / M_{3})^{1/3} \sim u^{1/3}$ , and back-substitution gives a floor of order $δ^{2/3} \sim u^{2/3}$ . $□$

Proposition 3 (forward mode computes $J \overset{x}{˙}$ via dual numbers). Running an evaluation procedure for $f : R^{n} \to R^{m}$ on dual numbers over the ring $R [ε] / (ε^{2})$ with input $x + \overset{x}{˙} ε$ returns $f (x) + J (x) \overset{x}{˙} ε$ .

Proof. The dual numbers $R [ε] / (ε^{2})$ form a commutative ring in which $ε^{2} = 0$ . Each elementary smooth operation $φ$ extends to duals by its first-order Taylor jet, $φ (v + \overset{v}{˙} ε) = φ (v) + φ^{'} (v) \overset{v}{˙} ε$ , since higher terms carry $ε^{2} = 0$ ; for a binary $φ (u, v)$ this gives $φ (u, v) + (\partial_{u} φ \overset{u}{˙} + \partial_{v} φ \overset{v}{˙}) ε$ , the chain rule on the $ε$ -coefficient. By induction on the length of the Wengert list, after processing $v_{k}$ its dual carries $(v_{k}, \overset{v}{˙}_{k})$ with $\overset{v}{˙}_{k} = \sum_{j ≺ k} (\partial φ_{k} / \partial v_{j}) \overset{v}{˙}_{j}$ , which is the $k$ -th component of the propagated directional derivative. At the outputs the $ε$ -coefficients assemble into $J (x) \overset{x}{˙}$ . $□$

Proposition 4 (reverse mode computes $\nabla f$ in $O (1)$ function evaluations). For $f : R^{n} \to R$ given by a procedure of $L$ elementary operations each of bounded arity, one forward sweep followed by one reverse adjoint sweep yields the exact $\nabla f (x)$ in total work $\leq ω_{grad} TIME (f)$ with $ω_{grad}$ an absolute constant.

Proof. Initialise $\overset{y}{ˉ} = 1$ and $\overset{v}{ˉ}_{k} = 0$ for $k < L$ . Process the operations in reverse order; when handling $v_{k} = φ_{k} (v_{j} : j ≺ k)$ , add $\overset{v}{ˉ}_{k} (\partial φ_{k} / \partial v_{j})$ to each $\overset{v}{ˉ}_{j}$ . By reverse induction this maintains the invariant $\overset{v}{ˉ}_{k} = \partial y / \partial v_{k}$ (total derivative through all paths from $v_{k}$ to $y$ ), because every path from $v_{k}$ to the sink passes through exactly one of $v_{k}$ 's immediate successors and the partials multiply along edges and add over paths — the multivariate chain rule. At termination $\overset{v}{ˉ}_{1}, \dots, \overset{v}{ˉ}_{n} = \partial y / \partial x_{1}, \dots, \partial y / \partial x_{n} = \nabla f (x)$ . Each operation is touched once in the forward sweep and once in the reverse sweep, and its reverse contribution is a fixed number of multiply-adds determined by its bounded arity, so the reverse work is at most a constant multiple of the forward work; with the forward sweep itself this gives total work $\leq ω_{grad} TIME (f)$ , the constant absorbing the per-operation local-partial cost and independent of $n$ . $□$

Proposition 5 (binomial checkpointing bound). A computation structured as a chain of $L$ steps can be reversed using $s$ stored checkpoints and at most $t$ extra forward recomputations whenever $(s s + t) \geq L$ ; choosing $s = t = O (lo g L)$ reverses the chain in memory $O (lo g L)$ and time $O (L lo g L)$ .

Proof. Let $c (s, t)$ be the largest chain length reversible with $s$ checkpoints and $t$ recomputation passes. Place the first checkpoint after $c (s, t - 1)$ steps: the suffix is reversed with $s$ checkpoints and $t - 1$ passes, the prefix with the remaining $s - 1$ checkpoints and $t$ passes once its end state is recomputed, giving the recurrence $c (s, t) = c (s - 1, t) + c (s, t - 1)$ with $c (s, 0) = c (0, t) = 1$ , whose solution is the binomial coefficient $c (s, t) = (s s + t)$ . Hence any $L \leq (s s + t)$ is reversible within the stated resources. Taking $s = t = k$ gives $(k 2 k) \sim 4^{k} / π k \geq L$ for $k = O (lo g L)$ , so memory $s = O (lo g L)$ and total recomputation work $O (L lo g L)$ . The schedule is optimal: no reversal of a length- $L$ chain with $s$ checkpoints and $t$ passes can exceed $(s s + t)$ steps, by the same recurrence read as an upper bound. $□$

Connections Master

Multivariable Taylor expansion and extrema 02.05.05 supplies the remainder estimates that make the finite-difference analysis quantitative: the $O (h)$ forward and $O (h^{2})$ central truncation terms are the Lagrange remainders of the first- and second-order expansions, and the same second-order Taylor data — gradient and Hessian — is exactly what the optimizers downstream need and what AD delivers exactly. The vanishing-gradient condition for an extremum from that unit is the target every gradient-based method of this chapter drives toward, and the curvature (Hessian) it classifies critical points with is the object Hessian-vector products approximate or, through forward-over-reverse AD, compute exactly.
Floating-point arithmetic and the IEEE model 43.01.01 is the source of the accuracy floor: the unit roundoff $u$ that bounds each evaluation's relative error is precisely what forces the $2 δ / h$ roundoff branch and caps differenced derivatives at $u$ . The catastrophic cancellation that the IEEE model predicts when subtracting nearby quantities is the exact mechanism that limits finite differences and that automatic differentiation sidesteps by never forming such a difference, so this unit is the applied payoff of that unit's error model.
The conjugate gradient method 43.07.04 and Newton-Krylov optimization consume Hessian-vector products of the kind built here: a truncated-Newton method solves the Newton system $\nabla^{2} f (x) p = - \nabla f (x)$ by CG, and CG needs only the action $v \mapsto \nabla^{2} f (x) v$ , supplied either by differencing the gradient (with the $u$ floor of the Key theorem) or, to full precision, by forward-over-reverse AD. The cheap-gradient principle is what makes the per-CG-iteration cost a small multiple of one function evaluation, so the matrix-free curvature access of that unit and the matrix-free derivative access of this one compose into large-scale second-order optimization.
The artificial-intelligence and machine-learning survey 25.09.01 pending is the large-scale consumer of reverse-mode AD: backpropagation is reverse-mode AD applied to a scalar training loss, and the cheap-gradient bound proved here is the reason a network with billions of parameters can be trained at all. This unit is the theorem-level cost-and-accuracy treatment that the survey references at a higher level; the survey covers the modelling and architectural context while this unit owns the derivative-computation cost theory, so the two are complementary rather than overlapping.

Historical & philosophical context Master

The finite-difference quotient is as old as the calculus, but its numerical-analysis treatment — the explicit balance of truncation against roundoff and the resulting optimal step near $u$ — belongs to the mid-twentieth-century error-analysis tradition associated with James H. Wilkinson, whose backward-error framework of the 1960s made the $u / h$ cancellation term a quantitative object. The forward mode of automatic differentiation was given its modern form by Robert E. Wengert in 1964, whose "simple automatic derivative evaluation program" introduced the list of intermediate operations now called the Wengert list ^{[Wengert 1964]}.

The reverse mode has a more tangled lineage. Seppo Linnainmaa's 1970 master's thesis and 1976 paper described the reverse accumulation of rounding errors over a computational graph — the algorithm now called backpropagation — predating its 1986 popularization in neural networks by Rumelhart, Hinton, and Williams ^{[Linnainmaa 1976]}. Bert Speelpenning's 1980 dissertation proved that the reverse mode computes a gradient at a cost independent of the number of inputs, the cheap-gradient principle. The systematic theory — complexity bounds, optimal Jacobian accumulation, and the binomial checkpointing schedule — was consolidated by Andreas Griewank and Andrea Walther, the latter's Revolve algorithm of 1996 establishing the optimal logarithmic time-memory tradeoff for reversal ^{[Griewank, A. & Walther, A. — Evaluating Derivatives]}. The 2018 survey of Baydin, Pearlmutter, Radul, and Siskind documented how reverse-mode AD, rediscovered repeatedly across numerical analysis and machine learning, became the computational substrate of deep learning ^{[Baydin, A. G., Pearlmutter, B. A., Radul, A. A. & Siskind, J. M. — Automatic Differentiation in Machine Learning: a Survey]}.

Bibliography Master

@book{NocedalWright2006,
  author    = {Nocedal, Jorge and Wright, Stephen J.},
  title     = {Numerical Optimization},
  edition   = {2},
  series    = {Springer Series in Operations Research and Financial Engineering},
  publisher = {Springer},
  year      = {2006}
}

@book{GriewankWalther2008,
  author    = {Griewank, Andreas and Walther, Andrea},
  title     = {Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation},
  edition   = {2},
  publisher = {Society for Industrial and Applied Mathematics},
  year      = {2008}
}

@article{BaydinEtAl2018,
  author  = {Baydin, At{\i}l{\i}m G\"{u}ne\c{s} and Pearlmutter, Barak A. and Radul, Alexey Andreyevich and Siskind, Jeffrey Mark},
  title   = {Automatic Differentiation in Machine Learning: a Survey},
  journal = {Journal of Machine Learning Research},
  volume  = {18},
  number  = {153},
  year    = {2018},
  pages   = {1--43}
}

@article{Wengert1964,
  author  = {Wengert, Robert E.},
  title   = {A Simple Automatic Derivative Evaluation Program},
  journal = {Communications of the ACM},
  volume  = {7},
  number  = {8},
  year    = {1964},
  pages   = {463--464}
}

@article{Linnainmaa1976,
  author  = {Linnainmaa, Seppo},
  title   = {Taylor Expansion of the Accumulated Rounding Error},
  journal = {BIT Numerical Mathematics},
  volume  = {16},
  number  = {2},
  year    = {1976},
  pages   = {146--160}
}

@article{GriewankWalther2000revolve,
  author  = {Griewank, Andreas and Walther, Andrea},
  title   = {Algorithm 799: Revolve: An Implementation of Checkpointing for the Reverse or Adjoint Mode of Computational Differentiation},
  journal = {ACM Transactions on Mathematical Software},
  volume  = {26},
  number  = {1},
  year    = {2000},
  pages   = {19--45}
}

@book{Wilkinson1963,
  author    = {Wilkinson, James H.},
  title     = {Rounding Errors in Algebraic Processes},
  publisher = {Prentice-Hall},
  year      = {1963}
}

Prerequisites

02.05.05
43.01.01

Tier anchors

beginner: Nocedal & Wright 2006 Numerical Optimization 2e (Springer) §8.1 (the forward and central finite-difference formulas, the step-size dilemma between truncation and roundoff); 3Blue1Brown — 'Backpropagation calculus' (the chain rule run backwards to get every derivative in one sweep)
intermediate: Nocedal & Wright 2006 Numerical Optimization 2e (Springer) §8.1-8.2 (the $O (h)$ forward and $O (h^{2})$ central difference truncation error, the optimal step $h \sim u$ for gradients and $h \sim u^{1/3}$ for second differences, gradient and Hessian-vector products by differencing); Griewank & Walther 2008 Evaluating Derivatives 2e (SIAM) §3-4 (the forward and reverse modes of automatic differentiation on the computational graph and the cheap-gradient principle)
master: Griewank & Walther 2008 Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation 2e (SIAM) Ch. 3-4 (forward and reverse accumulation), Ch. 12 (reversal schedules and checkpointing); Nocedal & Wright 2006 Numerical Optimization 2e (Springer) Ch. 8 (finite-difference and automatic derivatives for optimization); Baydin, Pearlmutter, Radul & Siskind 2018 Automatic Differentiation in Machine Learning: a Survey (JMLR 18) for the ML-training context

References

Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.) · Springer Series in Operations Research and Financial Engineering (2006). Chapter 8 'Calculating Derivatives'. §8.1 develops finite-difference approximations: the forward-difference gradient $\partial f / \partial x_{i} \approx (f (x + h e_{i}) - f (x)) / h$ with truncation error $O (h)$ and the central-difference variant $(f (x + h e_{i}) - f (x - h e_{i})) / (2 h)$ with error $O (h^{2})$ ; the analysis of the total error as the sum of a truncation term growing like $h$ (or $h^{2}$ ) and a roundoff term shrinking like $u / h$ (with $u$ the unit roundoff), giving the optimal step $h^{⋆} \sim u$ for forward differences and $h^{⋆} \sim u^{1/3}$ for central differences and second differences, and the corresponding best attainable accuracy $u$ and $u^{2/3}$ ; Hessian-vector products $\nabla^{2} f (x) v \approx (\nabla f (x + h v) - \nabla f (x)) / h$ by differencing the gradient, and full finite-difference Hessians exploiting sparsity by graph colouring. §8.2 introduces automatic (algorithmic) differentiation: the computational graph of a function, the forward mode propagating directional derivatives (tangents) and the reverse mode (adjoint / backpropagation) propagating sensitivities backward, with the cost theorem that one gradient of a scalar function of $n$ variables costs a small constant multiple (typically $\leq 5$ ) of one function evaluation independent of $n$ , in contrast to the $n + 1$ evaluations of a finite-difference gradient; the contrast with symbolic differentiation and the memory cost of storing the forward trace for reversal.
Griewank, A. & Walther, A. — Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation (2nd ed.) · SIAM (2008). Chapters 3-4 give the formal account of algorithmic differentiation on the evaluation procedure (Wengert list / computational graph): the forward (tangent) mode as the propagation of dual-number / directional-derivative pairs $(v, \overset{v}{˙})$ obeying the elementary rules $\dot{(uv)} = \overset{u}{˙} v + u \overset{v}{˙}$ etc., computing a Jacobian-vector product $J \overset{x}{˙}$ at cost a small multiple of the function; the reverse (adjoint, cotangent) mode propagating adjoints $\overset{v}{ˉ}_{i} = \partial y / \partial v_{i}$ backward through the graph, computing a vector-Jacobian product $\overset{y}{ˉ}^{⊤} J$ and hence a full gradient of a scalar at cost $\leq ω_{grad}$ times the function cost with $ω_{grad} \in [3, 4]$ independent of the input dimension $n$ — the 'cheap gradient principle'; the temporal-complexity and spatial-complexity (memory) bounds; and Chapter 12 on reversal schedules and checkpointing, where storing a logarithmic number of intermediate states trades a factor $O (lo g L)$ in time for an $O (lo g L)$ reduction in the memory needed to reverse a chain of length $L$ (the binomial / Revolve schedule). The contrast with symbolic differentiation (expression swell) and numerical differentiation (truncation/roundoff) is drawn explicitly.
Baydin, A. G., Pearlmutter, B. A., Radul, A. A. & Siskind, J. M. — Automatic Differentiation in Machine Learning: a Survey · Journal of Machine Learning Research 18 (2018), 1-43. Surveys forward- and reverse-mode automatic differentiation, distinguishes AD from both symbolic and numerical differentiation, explains backpropagation as the special case of reverse-mode AD applied to a scalar loss, and documents the cheap-gradient property ( $O (1)$ function evaluations for the gradient of a scalar regardless of the number of parameters) that makes reverse-mode AD foundational to large-scale neural-network training; covers operator-overloading and source-transformation implementations, dual numbers for the forward mode, and the taping / Wengert-list mechanism for the reverse mode.

Estimated time

beginner: 18m
intermediate: 48m
master: 90m