The Elements of Differentiable Programming
Part III Differentiating through programs

Chapter 7Finite differences

One of the simplest ways to numerically compute derivatives is to use finite differences, which approximate the infinitesimal definition of derivatives. Finite differences only require function evaluations, and can therefore work with blackbox functions: they ignore the compositional structure of functions. Without loss of generality, our exposition focuses on computing directional derivatives ∂f(w)[v]\partial f(\w)[\v], for a function f ⁣:E→Ff \colon \cE \to \cF, evaluated at w∈E\w \in \cE, in the direction v∈E\v \in \cE.

7.1 Forward differences

From Definition 2.4 and Definition 2.13, the directional derivative and more generally the JVP are defined as a limit, ∂f(w)[v]≔lim⁡δ→0f(w+δv)−f(w)δ. \partial f(\w)[\v] \coloneqq \lim_{\delta \rightarrow 0} \frac{f(\w + \delta\v) -f(\w)}{\delta}. This suggests that we can approximate the directional derivative and the JVP using ∂f(w)[v]≈f(w+δv)−f(w)δ, \partial f(\w)[\v] \approx \frac{f(\w + \delta\v) - f(\w)}{\delta}, for some 0<δ≪10 < \delta \ll 1. This formula is called a forward difference.

To measure the error incurred by this approximation, we use Landau’s OO notation.

Remark 7.1 (Big OO notation)

For two functions gg and ff, we write g(δ)=O(f(δ)) as δ→γ, g(\delta) = O(f(\delta)) \ \text{as} \ \delta \to \gamma, if the ratio ∣g(δ)∣/∣f(δ)∣|g(\delta)|/|f(\delta)| is bounded as δ→γ\delta \to \gamma. For example, we write δcos⁡(δ)=O(δ)\delta \cos(\delta) = O(\delta) as δ→0\delta \to 0. In the following, unless otherwise specified, we compare limits as δ→0\delta \to 0. We observe that g(δ)=O(δk)g(\delta) = O(\delta^k) implies g(δ)=o(δk−1)g(\delta) = o(\delta^{k-1}). However, Landau’s OO notation provides finer practical interpretations than Landau’s oo notation. For example, if g(δ)=O(δ2)g(\delta) = O(\delta^2), then we know its upper bound scales quadratically, meaning that reducing δ\delta by a factor of 2 reduces this bound by a factor of 4. Knowing only that g(δ)=o(δ)g(\delta) = o(\delta) does not allow us to make such a concrete statement, because it merely dictates that g(δ)g(\delta) decays strictly faster than δ\delta, without providing a specific bounding rate or scaling factor for that decay. For numerical analysis, Landau’s OO notation is therefore preferable.

From the Taylor expansion in Section 2.5.4, we indeed have f(w+δv)−f(w)=δ∂f(w)[v]+δ22∂2f(w)[v,v]+δ33!∂3f(w)[v,v,v]+… f(\w+\delta\v) - f(\w) = \delta \partial f(\w)[\v] + \frac{\delta^2}{2} \partial^2 f(\w)[\v, \v] + \frac{\delta^3}{3!} \partial^3 f(\w)[\v, \v, \v] + \dots so that f(w+δv)−f(w)δ=∂f(w)[v]+δ2∂2f(w)[v,v]+δ23!∂3f(w)[v,v,v]+…=∂f(w)[v]+O(δ).\begin{align*} \frac{f(\w+\delta\v) - f(\w)}{\delta} &= \partial f(\w)[\v] + \frac{\delta}{2} \partial^2 f(\w)[\v, \v] + \frac{\delta^2}{3!} \partial^3 f(\w)[\v, \v, \v] + \dots \\ &= \partial f(\w)[\v] + O(\delta). \end{align*} The error incurred by choosing a finite rather than infinitesimal δ\delta in the forward difference formula is called the truncation error. The Taylor approximation above shows that this error is of the order of O(δ)O(\delta).

However, we cannot choose a too small value of δ\delta, because the evaluation of the function ff on a computer rounds the value of ff to machine precision. Mathematically, a scalar-valued function ff evaluated on a computer becomes a function f~\tilde{f} such that f~(w)≈[f(w)/ε]ε\tilde{f}(\w) \approx [f(\w)/\varepsilon] \varepsilon, where [f(w)/ε][f(\w)/\varepsilon] denotes the closest integer of f(w)/ε∈Rf(\w)/\varepsilon \in \R and ε\varepsilon is the machine precision, i.e., the smallest non-zero real number encoded by the machine. This means that the difference f(w+δv)−f(w)f(\w+\delta\v) - f(\w) evaluated on a computer is prone to round-off error of the order of O(ε)O(\varepsilon). We illustrate the trade-off between truncation and round-off errors in Fig. 7.1.

Figure 7.1: Numerical differentiation of f(x)≔softplus(x)=log⁡(1+exp⁡(x))f(x) \coloneqq \mathrm{softplus}(x) = \log(1+\exp(x)), to approximate f′(x)=logistic(x)f'(x) = \mathrm{logistic}(x) at x=1x=1. The forward and central difference methods induce both truncation error (for large δ\delta) and round-off error (for small δ\delta). The complex step method enjoys smaller round-off error.

7.2 Backward differences

As an alternative, we can approximate the directional derivative and the JVP by ∂f(w)[v]≈f(w)−f(w−δv)δ, \partial f(\w)[\v] \approx \frac{f(\w) - f(\w - \delta\v)}{\delta}, for some 0<δ≪10 < \delta \ll 1. This formula is called a backward difference. From the Taylor expansion in Section 2.5.4, we easily verify that (f(w)−f(w−δv))/δ=∂f(w)[v]+O(δ)(f(\w) - f(\w-\delta\v))/\delta = \partial f(\w)[\v] + O(\delta), so that the truncation error is the same as for the forward difference.

7.3 Central differences

Rather than using an asymmetric formula to approximate the derivative, as in forward and backward differences, we can use a symmetric formula ∂f(w)[v]≈f(w+δv)−f(w−δv)2δ, \partial f(\w)[\v] \approx \frac{f(\w + \delta\v) - f(\w-\delta\v)}{2\delta}, for some 0<δ≪10 < \delta \ll 1. This formula is called a central difference. From the Taylor expansion in Section 2.5.4, we have f(w+δv)−f(w−δv)=2δ∂f(w)[v]+2δ33!∂3f(w)[v,v,v]+2δ55!∂5f(w)[v,…,v]+…\begin{align*} f(\w+\delta\v) - f(\w-\delta\v) =& 2\delta \partial f(\w)[\v] + \frac{2\delta^3}{3!} \partial^3 f(\w)[\v, \v, \v] \\ &+ \frac{2\delta^5}{5!} \partial^5 f(\w)[\v, \dots, \v] + \dots \end{align*} so that f(w+δv)−f(w−δv)2δ=∂f(w)[v]+δ23!∂3f(w)[v,v,v]+…=∂f(w)[v]+O(δ2).\begin{align*} \frac{f(\w+\delta\v) - f(\w-\delta\v)}{2\delta} =& \partial f(\w)[\v] + \frac{\delta^2}{3!} \partial^3 f(\w)[\v, \v, \v] + \ldots\\ =& \partial f(\w)[\v] + O(\delta^2). \end{align*} We see that the terms corresponding to derivatives of even order canceled out, allowing the formula to achieve O(δ2)O(\delta^2) truncation error. For any δ<1\delta<1, the truncation error of the central difference is much smaller than the one of the forward or backward differences as confirmed empirically in Fig. 7.1.

7.4 Higher-accuracy finite differences

The truncation error can be further reduced by making use of additional function evaluations. One can generalize the forward difference scheme by a formula of the form ∂f(w)[v]≈∑t=0patδf(w+tδv) \partial f(\w)[\v] \approx \sum_{t=0}^p \frac{a_t}{\delta} f(\w + t\delta \v) requiring p+1p+1 evaluations. To select the ata_t and reach a truncation error of order O(δp)O(\delta^p), we can use a Taylor expansion on each term of the sum to get ∑t=0patδf(w+tδv)=∑t=0pat∑j=0ptjδj−1j!∂jf(w)[v,…,v]+O(δp).\begin{align*} \sum_{t=0}^p \frac{a_t}{\delta} f(\w + t\delta \v) & = \sum_{t=0}^p a_t \sum_{j=0}^p \frac{t^j\delta^{j-1}}{j!} \partial^j f(\w)[\v, \ldots, \v] + O(\delta^p). \end{align*} By grouping the terms in the sum for each order of derivative, we obtain a set of p+1p+1 equations to be satisfied by the p+1p+1 coefficients a0,…,apa_0, \ldots, a_p, that is, a0+a1+…+ap=0a1+2a2+…+pap=1a1+2ja2+…+pjap=0∀j∈{2,…,p}.\begin{align*} a_0 + a_1 + \ldots + a_p & = 0 \\ a_1 + 2 a_2 + \ldots + p a_p & = 1 \\ a_1 + 2^j a_2 + \ldots + p^j a_p & = 0 \quad \forall j \in \{2, \ldots, p\}. \end{align*} This system of equations can be solved analytically to derive the coefficients. Backward differences can be generalized similarly by using ∂f(w)[v]≈∑t=0patδf(w−tδv)\partial f(\w)[\v] \approx \sum_{t=0}^p \frac{a_t}{\delta} f(\w - t\delta \v). Similarly, the central difference scheme can be generalized by using ∂f(w)[v]≈∑t=−ppatδf(w+tδv), \partial f(\w)[\v] \approx \sum_{t=-p}^p \frac{a_t}{\delta} f(\w + t\delta \v), to reach a truncation error of order O(δ2p)O(\delta^{2p}). Solving for the coefficients a−p,…,apa_{-p}, \ldots, a_p as above reveals that a0=0a_0 = 0. Therefore, only 2p2p evaluations are necessary.

7.5 Higher-order finite differences

To approximate higher order derivatives, we can follow a similar reasoning. Namely, we can generalize the forward difference scheme to approximate the derivative of order kk by ∂kf(w)[v,…,v]≈∑t=0patδkf(w+tδv). \partial^k f(\w)[\v, \ldots, \v] \approx \sum_{t=0}^p \frac{a_t}{\delta^k}f(\w + t \delta \v). As before, we can expand the terms in the sum. For the approximation to capture only the kkth derivative, we now require the coefficients ata_t to satisfy 0ja0+1ja1+2ja2+…+pjap=0∀j∈{0,…,k−1}.0ka0+1ka1+2ka2+…+pkap=k!0ja0+1ja1+2ja2+…+pjap=0∀j∈{k+1,…,p}.\begin{align*} 0^j a_0 + 1^j a_1 + 2^j a_2 + \ldots + p^j a_p & = 0 \quad \forall j \in \{0, \ldots, k-1\}. \\ 0^k a_0 + 1^k a_1 + 2^k a_2 + \ldots + p^k a_p & = k! \\ 0^j a_0 + 1^j a_1 + 2^j a_2 + \ldots + p^j a_p & = 0 \quad \forall j \in \{k+1, \ldots, p\}. \end{align*} With the resulting coefficients, we obtain a truncation error of order O(δp−k+1)O(\delta^{p - k +1}), while making p+1p+1 evaluations. For example, for p=k=2p=k=2, we can approximate the second-order derivative as ∂2f(w)[v,v]≈−(3/2)f(w)+2f(w+δv)−(1/2)f(w+2δv)δ2, \partial^2 f(\w)[\v, \v] \approx \frac{-(3/2)f(\w) +2 f(\w+ \delta \v) - (1/2) f(\w + 2\delta \v)}{\delta^2}, with a truncation error of order O(δ)O(\delta).

The central difference scheme can be generalized similarly by ∂kf(w)[v,…,v]≈∑t=−ppatδkf(w+tδv), \partial^k f(\w)[\v, \ldots, \v] \approx \sum_{t=-p}^p \frac{a_t}{\delta^k} f(\w + t\delta \v), to reach truncation errors of order O(δ2p+2−2⌈(k+1)/2⌉)O(\delta^{2p +2 - 2\lceil (k+1)/2\rceil}). For example, for k=2k=2, p=1p=1, we obtain the second-order central difference ∂2f(w)[v,v]≈f(w+δv)+f(w−δv)−2f(w)δ2. \partial^2 f(\w)[\v, \v] \approx \frac{f(\w+\delta\v) + f(\w-\delta\v) - 2f(\w)}{\delta^2}. By using a Taylor expansion we see that, this time, the terms corresponding to derivatives of odd order cancel out and the truncation error is O(δ2)O(\delta^2) while requiring 33 evaluations.

7.6 Complex-step derivatives

Suppose ff is well defined on CP\CC^P, the space of PP-dimensional complex numbers. Let us denote the imaginary unit by i=−1i = \sqrt{-1}. Then, the Taylor expansion of ff reads f(w+(iδ)v)=f(w)+(iδ)∂f(w)[v]+(iδ)22∂2f(w)[v,v]+(iδ)33!∂3f(w)[v,v,v]+…=f(w)+(iδ)∂f(w)[v]−δ22∂2f(w)[v,v]−iδ33!∂3f(w)[v,v,v]+….\begin{align*} f(\w+(i \delta)\v) =& f(\w) + (i\delta) \partial f(\w)[\v] + \frac{(i\delta)^2}{2} \partial^2 f(\w)[\v, \v] \\ &+ \frac{(i\delta)^3}{3!} \partial^3 f(\w)[\v, \v, \v] + \dots \\ =& f(\w) + (i\delta) \partial f(\w)[\v] - \frac{\delta^2}{2} \partial^2 f(\w)[\v, \v] \\ &- \frac{i\delta^3}{3!} \partial^3 f(\w)[\v, \v, \v] + \dots \quad . \end{align*} We see that the real part corresponds to even-degree terms and the imaginary part corresponds to odd-degree terms. We therefore obtain Re⁡(f(w+(iδ)v))=f(w)+O(δ2) \operatorname{Re}(f(\w+(i \delta)\v)) = f(\w) + O(\delta^2) and Im⁡(f(w+(iδ)v)δ)=∂f(w)[v]+O(δ2). \operatorname{Im}\left(\frac{f(\w+(i \delta)\v)}{\delta}\right) = \partial f(\w)[\v] + O(\delta^2). This suggests that we can compute directional derivatives using the approximation ∂f(w)[v]≈Im⁡(f(w+(iδ)v)δ), \partial f(\w)[\v] \approx \operatorname{Im}\left(\frac{f(\w+(i \delta)\v)}{\delta}\right), for 0<δ≪10 < \delta \ll 1. This is called the complex-step derivative approximation (Martins et al., 2003; Squire and Trapp, 1998).

Contrary to forward, backward and central differences, we see that only a single function call is necessary. A function call on complex numbers may take roughly twice the cost of a function call on real numbers. However, thanks to the fact that a difference of functions is no longer needed, the complex-step derivative approximation usually enjoys smaller round-off error as illustrated in Fig. 7.1. That said, one drawback of the method is that all elementary operations within the program implementing the function ff must be well-defined on complex numbers, e.g., using overloading.

7.7 Complexity

We now discuss the computational complexity in terms of function evaluations of finite differences and complex-step derivatives. For concreteness, as this is the most common use case in machine learning, we discuss the case of a single M=1M=1 output, i.e., we want to differentiate a function f ⁣:RP→Rf \colon \RR^P \to \RR. Whether we use forward, backward or central differences, the computational complexity of computing the directional derivative ∂f(w)[v]\partial f(\w)[\v] in any direction v\v amounts to two calls to ff. For computing the gradient ∇f(w)\nabla f(\w), we can use (see Definition 2.7) that [∇f(w)]j=⟨∇f(w),ej⟩=∂f(w)[ej], [\nabla f(\w)]_j = \langle \nabla f(\w), \e_j \rangle = \partial f(\w)[\e_j], for j∈[P]j \in [P]. For forward and backward differences, we therefore need P+1P+1 function calls to compute the gradient, while we need 2P2P function calls for central differences. For the complex step approximation, we need PP complex function calls. We summarize the complexities in Table 7.1.

Table 7.1: Computational complexity in number of function evaluations for computing the directional derivative and the gradient of a function f ⁣:RP→Rf \colon \RR^P \to \RR by finite differences and complex step derivatives.

Directional derivative Gradient
Forward difference 22 P+1P+1
Backward difference 22 P+1P+1
Central difference 22 2P2P
Complex step 11 PP

7.8 Summary