A slope without a subtraction

The slope of a single point defined the derivative as the limit of secant slopes, and then noticed that a computer cannot follow the definition. The quotient $(f(x + h) - f(x))/h$ subtracts two numbers that agree in their leading digits, and the digits they share are destroyed: shrink $h$ and the truncation error falls while the rounding error rises, and the best compromise, near $h = 10^{-8}$ in double precision, still loses half the digits the machine carries. That essay mentioned, in one sentence, that real software does not compute slopes this way, and named three methods that never form the fatal difference.

Two of those methods are worth more than a sentence, because they are not tricks for computing derivatives so much as other definitions of the derivative — ones in which no limit is taken and nothing small is divided by. One uses the complex numbers; the other uses a number system of its own, in which a nonzero number squares to nought. The second, run backwards, is how every large optimisation in present use gets its slopes.

The test function, and the subtraction’s price #

The function in the figures is $f(x) = e^x/\sqrt{\sin^3 x + \cos^3 x}$, which William Squire and George Trapp used in 1998 to show off the first of the two methods, evaluated at $x = 1.5$. Its slope there is $4.053427893898620$, computed from the formula for $f'$ that ordinary calculus gives; every error is measured against that.

The forward difference behaves as the first essay predicted. As $h$ shrinks from $10^{-1}$ its error falls like $h$, reaches a minimum of $2.6 \times 10^{-9}$ near $h = 10^{-8}$, and then rises in a jagged staircase as more and more of $f(x + h)$ and $f(x)$ coincide; below $h = 10^{-16}$ the two values are the same floating-point number and the quotient is nought, an error of 100%. The central difference $(f(x + h) - f(x - h))/2h$ does better, bottoming out at $6.4 \times 10^{-12}$ near $h = 3 \times 10^{-6}$, because its truncation error falls like $h^2$ and the balance with rounding is struck at a larger step. Then it climbs the same staircase. The same fight between a vanishing correction and a growing rounding error spoiled computing $e$ from compound interest, and it is the general hazard of any limit taken in floating point, as the orbit a computer draws found for chaos. No choice of step rescues either one, because the loss is in the subtraction, and the subtraction is the definition.

Stepping sideways #

The complex step evaluates $f$ at $x + ih$ — a step of size $h$ in the imaginary direction — and divides the imaginary part of the answer by $h$:

$$f'(x) \approx \frac{\operatorname{Im} f(x + ih)}{h}.$$

Taylor’s theorem says why. For a function given by a power series with real coefficients,

$$f(x + ih) = f(x) + ih\,f'(x) - \frac{h^2}{2} f''(x) - \frac{ih^3}{6} f'''(x) + \cdots,$$

so the real parts collect the even terms and the imaginary parts the odd ones, and $\operatorname{Im} f(x + ih)/h = f'(x) - h^2 f'''(x)/6 + \cdots$. The error is of order $h^2$, like the central difference, and the crucial difference is what is not there: $f(x)$ sits in the real part and never meets $f'(x)$, so nothing is subtracted from anything of similar size. Rounding error stays at the level of a single evaluation however small $h$ is.

The figure shows exactly that. The complex step’s error falls like $h^2$, reaches the rounding level of about $10^{-16}$ at $h = 10^{-8}$, and stays there for every smaller step plotted — and at $h = 10^{-200}$, far past the edge of the figure, it is still correct to the last digit. The only requirement is that the program computing $f$ accept complex arguments and do nothing that is not analytic: the sine of a complex number, its exponential, its square root, its product and quotient are all fine; an absolute value or a comparison of sizes is not, because those are not given by power series. Squire and Trapp’s paper made the method known in 1998, though James Lyness and Cleve Moler had found the underlying idea in 1967, and it is now the standard way to check derivatives that other methods have produced.

A number whose square is nought #

The second method changes the arithmetic itself. Take numbers of the form $a + b\varepsilon$, where $\varepsilon$ is a new symbol subject to the single rule $\varepsilon^2 = 0$. They add in the obvious way, and they multiply by

$$(a + b\varepsilon)(c + d\varepsilon) = ac + (ad + bc)\varepsilon,$$

which is the product rule — the corner of the rectangle grown on two sides is the $\varepsilon^2$ term, and here it is nought by decree rather than by a limit. These are dual numbers. For any polynomial, $p(a + \varepsilon) = p(a) + p'(a)\varepsilon$ exactly, because expanding the powers leaves only the terms with no $\varepsilon$ or one. Define every other function by the same rule — $\sin(a + b\varepsilon) = \sin a + b\cos a\,\varepsilon$, and so on — and every composition obeys the chain rule automatically, because the chain rule is what composing these definitions does.

The trace shows the whole computation. The input is $1.5 + 1\varepsilon$. The sine becomes $0.99749 + 0.07074\varepsilon$, its square and cube follow by the product rule, the cosine’s cube likewise, the sum, the square root by its rule, the exponential by its own, and the final quotient by the quotient rule. The last line’s $\varepsilon$ part is $4.05342789$, the slope, and in the hero figure the dual number’s error is at rounding level with no step to choose. Nothing in the calculation is a limit, and nothing is divided by anything small. The derivative has become arithmetic in a slightly larger number system, the way complex multiplication became turning once $i^2 = -1$ was adopted.

The two methods are closer relatives than they look. A complex number $x + ih$ with $h$ tiny is almost a dual number: its imaginary part squares to $-h^2$ rather than to nought, and $-h^2$ is too small to matter once $h$ is below about $10^{-8}$. So the complex step is dual-number arithmetic with a small leak, and the leak is exactly its $h^2$ error term — the second-order part of the Taylor series reaching the imaginary part through $i^2 = -1$. Shrinking $h$ to $10^{-200}$ closes the leak to far below rounding, which is why the method works with absurdly small steps: it is simulating a number system in which the square of the step is nought, using one in which it is merely negligible. The two differ in what they ask of the program. The complex step needs nothing but complex arithmetic, which most languages already provide, and that is its appeal for checking derivatives; dual numbers need every operation redefined, and in exchange they are exact and they never need a step.

There is a cost, and it is honest: every operation now does about three times the work, because the $\varepsilon$ part must be updated alongside the value. For a function of one variable that is a bargain. For a function of a thousand variables, one dual sweep gives one partial derivative, and a thousand sweeps are needed for the gradient.

Second derivatives cost twice #

The subtraction does more damage one derivative higher. The second difference $(f(x + h) - 2f(x) + f(x - h))/h^2$ subtracts three nearly equal numbers and divides by $h^2$, so rounding error grows like $\epsilon/h^2$ while truncation falls like $h^2$; the balance is struck near $h = 10^{-4}$, and the figure finds the best error at $2.5 \times 10^{-9}$ — about half the digits again, but from a step a thousand times larger, and with no margin around it. The complex step does not help here: the second derivative lives in its real part, $f(x) - \tfrac12 h^2 f''(x)$, and extracting it subtracts $f(x)$ after all, so it fares slightly worse, at $1.7 \times 10^{-8}$.

Dual numbers extend without difficulty. Nest them — let the coefficients of a dual number themselves be dual numbers, with two independent symbols $\varepsilon_1$ and $\varepsilon_2$ each squaring to nought — and the coefficient of $\varepsilon_1\varepsilon_2$ in $f(x + \varepsilon_1 + \varepsilon_2)$ is the second derivative. One evaluation in this arithmetic gives $f''(1.5) = 9.46307368$ to rounding, agreeing with a fourth-order difference of the exact first derivative to ten places. Any order of derivative can be had the same way, at a cost that grows with the order but never with any loss of digits.

Running the chain rule backwards #

A derivative in many variables is a linear map, and the chain rule says the derivative of a composition is the product of the derivatives of its steps. A computation of $f$ from $n$ inputs is a long chain of elementary steps, so its gradient is a long product of small matrices, one per step. Forward mode multiplies that product starting from the inputs: it pushes one direction at a time through the chain, which is why it needs $n$ sweeps for $n$ partial derivatives. But matrix products can be multiplied in either order. Starting from the output — one number — and multiplying backwards, each step receives a row vector of sensitivities, “how much the final answer changes per unit change in this intermediate quantity”, and passes it to its own inputs. When the sweep reaches the start, the row vector is the whole gradient.

That backward multiplication is the transpose of the forward one, which is to say it is the adjoint of the derivative: the map that moves across an inner product, applied step by step. It needs the evaluation to be recorded first, so the backward sweep knows each step’s local slopes, and it does about one multiplication and one addition for every link it passes back along. The figure measures it on a function of $n$ variables built from $n - 1$ terms of Rosenbrock’s valley, each with a sine twist. Forward mode’s cost, counted as three real operations for each dual one, grows exactly in proportion to $n$: three thousand evaluations’ worth for a thousand variables. The reverse sweep costs 4.8 evaluations’ worth at every size from 2 variables to 1,000, and one of its thousand partial derivatives agrees with a forward dual sweep to rounding at every size.

That constant is a theorem, Walter Baur and Volker Strassen’s of 1983, and it is the fact on which modern machine learning runs. Training a network with a billion parameters needs a gradient with a billion components at every step, and the reverse sweep — called backpropagation there — delivers it for a small fixed multiple of the cost of running the network once. Forward differences would need a billion and one evaluations, each losing half its digits.

Who ran it backwards first #

Forward mode was programmed in 1964 by Robert Wengert, whose short paper computed derivatives by carrying them alongside values through a list of elementary operations — the trace in the second figure is his idea. The reverse sweep was found by Seppo Linnainmaa in 1970, in a master’s thesis about the accumulation of rounding errors: he wanted to know how sensitive a long computation’s result was to the error committed at every step, which is exactly the vector of sensitivities the backward pass computes. Bert Speelpenning built a program that produced reverse-mode code automatically in 1980, and David Rumelhart, Geoffrey Hinton and Ronald Williams made the same calculation the training method for neural networks in 1986 under the name backpropagation, without, at first, connecting it to the earlier work.

The repeated rediscovery is a sign of how natural the idea is once the chain rule is read as a statement about a computation rather than a formula. The same reading explains why Newton’s method can be applied to functions defined only by programs, and why the slope of an inverse can be had by differentiating the iteration that computes the inverse rather than by any formula — which is the subject of the last figure.

Differentiating an algorithm #

The methods differentiate programs, not formulas, and a program can be an algorithm that never mentions a derivative. Newton’s iteration $x \leftarrow (x + a/x)/2$ computes $\sqrt a$; run it on the dual number $a + \varepsilon$, starting from $x = 1$, and each iterate carries an approximation to $\sqrt a$ and its slope with respect to $a$. The figure follows both for $a = 2$. The value converges in the usual way, its correct digits doubling each step, and the slope converges to $1/(2\sqrt2)$, the derivative of the square root, doubling its digits too, one step behind. After five steps both are exact to rounding. The derivative of the algorithm has converged to the derivative of the function the algorithm computes.

That is not guaranteed in general, and the exceptions are instructive. A program that branches — “if $x$ is zero, return nought” — has derivatives only on each side of the branch, and dual numbers report whichever side was taken; a slope that swings like that of $x^2\sin(1/x)$ at the origin cannot be computed this way at all, because the formula divides by $x$. An iteration stopped after a fixed number of steps has a derivative that is the derivative of that finite computation, which can differ from the derivative of the limit; for Newton’s method they agree because the iteration converges so fast that its slope inherits the convergence. Differentiation is exact for the program as written, and whether the program’s derivative is the one wanted is a question about the program.

What the picture cannot show #

Every figure here is one function at one point. The error curves of the differences depend on the function and the point — on the sizes of $f''$ and $f'''$ and on how $f$’s values round — so the minima move from case to case, though their orders, $10^{-8}$ and $10^{-11}$ in double precision, do not. The operation counts for the reverse sweep are a model, three real operations to a dual one and two to each backward link, measured on one family of functions; Baur and Strassen’s theorem bounds the ratio for every function, with a constant that depends on how the elementary operations are counted.

The exactness of dual numbers is exactness up to rounding: each operation rounds, as an evaluation of $f$ does, and the derivative inherits the same few units in the last place. What the figures show is that no additional loss occurs — the derivative is as accurate as the function value — and that is the whole point.

Still open: the cheapest way to multiply a chain #

Forward and reverse are the two extreme orders for multiplying the chain of small matrices; any bracketing is allowed, and the cost depends on the order. Finding the cheapest order for computing a full matrix of derivatives — a Jacobian — is the optimal Jacobian accumulation problem, and Uwe Naumann proved in 2008 that it is NP-complete in a natural formulation: no fast method is known for finding the cheapest elimination order, though exhaustive or heuristic orders can beat both pure sweeps by large factors on particular programs.

A second question is about memory rather than time. The reverse sweep needs the forward computation’s intermediate values, and storing them all can cost as much memory as the computation had steps. Checkpointing schemes store only some and recompute the rest, and Andreas Griewank found the optimal schedule for a chain of fixed length in 1992; for programs whose length is not known in advance, or that run on many processors at once, the best trade between recomputation and storage is not settled, and it is what decides how large a model can be trained on a given machine.

The slope as arithmetic #

A difference quotient computes a slope by subtracting nearly equal numbers, and loses half the machine’s digits however the step is chosen; the complex step and dual numbers compute the same slope with no subtraction, and lose nothing. The complex step hides the slope in the imaginary part, where the function’s value cannot cancel it; dual numbers make the product rule a law of arithmetic by declaring $\varepsilon^2 = 0$, so that every composition carries its own slope through the chain rule. Running the chain rule from the output backwards — the adjoint of the derivative, step by step — gives every partial derivative of a function of a thousand variables for under five times the cost of the function, and differentiating Newton’s square-root iteration converges to the derivative of the square root. The definition by limits says what a slope is; these say how to have one.