Automatic differentiation

In mathematics and computer algebra, automatic differentiation (auto-differentiation, autodiff, or AD), also called algorithmic differentiation, computational differentiation, and differentiation arithmetic is a set of techniques to evaluate the partial derivative of a function specified by a computer program. Automatic differentiation enables the simultaneous computation of the numerical values of arbitrarily complex functions and their derivatives using only the function itself; no symbolic derivative is required. Auto-differentiation is thus neither numeric nor symbolic, nor is it a combination of both. In contrast to the more traditional numerical methods based on finite differences, auto-differentiation is 'in theory' exact, and in comparison to symbolic algorithms, it is computationally inexpensive. Automatic differentiation exploits the fact that every computer calculation, no matter how complicated, executes a sequence of elementary arithmetic operations (addition, subtraction, multiplication, division, etc.) and elementary functions (exp, log, sin, cos, etc.). By applying the chain rule repeatedly to these operations, partial derivatives of arbitrary order can be computed automatically, accurately to working precision, and using at most a small constant factor of more arithmetic operations than the original program.

In mathematics and computer algebra, automatic differentiation (auto-differentiation, autodiff, or AD), also called algorithmic differentiation, computational differentiation, and differentiation arithmetic is a set of techniques to evaluate the partial derivative of a function specified by a computer program. Automatic differentiation enables the simultaneous computation of the numerical values of arbitrarily complex functions and their derivatives using only the function itself; no symbolic derivative is required. Auto-differentiation is thus neither numeric nor symbolic, nor is it a combination of both. In contrast to the more traditional numerical methods based on finite differences, auto-differentiation is 'in theory' exact, and in comparison to symbolic algorithms, it is computationally inexpensive. Automatic differentiation exploits the fact that every computer calculation, no matter how complicated, executes a sequence of elementary arithmetic operations (addition, subtraction, multiplication, division, etc.) and elementary functions (exp, log, sin, cos, etc.). By applying the chain rule repeatedly to these operations, partial derivatives of arbitrary order can be computed automatically, accurately to working precision, and using at most a small constant factor of more arithmetic operations than the original program. \begin{align} y &= f(g(h(x))) = f(g(h(w_0))) = f(g(w_1)) = f(w_2) = w_3 \\ w_0 &= x \\ w_1 &= h(w_0) \\ w_2 &= g(w_1) \\ w_3 &= f(w_2) = y \end{align} the chain rule gives \frac{\partial y}{\partial x} = \frac{\partial y}{\partial w_2} \frac{\partial w_2}{\partial w_1} \frac{\partial w_1}{\partial x} = \frac{\partial f(w_2)}{\partial w_2} \frac{\partial g(w_1)}{\partial w_1} \frac{\partial h(w_0)}{\partial x} Forward accumulation specifies that one traverses the chain rule from inside to outside (that is, first compute \frac{\partial w_1}{\partial x} and then \frac{\partial w_2}{\partial w_1} and lastly \frac{\partial y}{\partial w_2} ), while reverse accumulation traverses from outside to inside (first compute \frac{\partial y}{\partial w_2} and then \frac{\partial w_2}{\partial w_1} and lastly \frac{\partial w_1}{\partial x} ). More succinctly, Forward accumulation computes the recursive relation: \frac{\partial w_i}{\partial x} = \frac{\partial w_i}{\partial w_{i-1}} \frac{\partial w_{i-1}}{\partial x} \quad\text{with } w_3 = y, Reverse accumulation computes the recursive relation: \frac{\partial y}{\partial w_i} = \frac{\partial y}{\partial w_{i+1}} \frac{\partial w_{i+1}}{\partial w_{i}} \quad\text{with } w_0 = x. The value of the partial derivative, called the seed, is propagated forward or backward and is initially \frac{\partial x}{\partial x}=1 or \frac{\partial y}{\partial y}=1 . Forward accumulation evaluates the function and calculates the derivative with respect to one independent variable in one pass. For each independent variable x_1,x_2,\dots,x_n a separate pass is therefore necessary in which the derivative with respect to that independent variable is set to one ( \frac{\partial x_1}{\partial x_1}=1 ) and of all others to zero ( \frac{\partial x_2}{\partial x_1}= \dots = \frac{\partial x_n}{\partial x_1} = 0 ). In contrast, reverse accumulation requires the evaluated partial functions for the partial derivatives. Reverse accumulation therefore evaluates the function first and calculates the derivatives with respect to all independent variables in an additional pass. Which of these two types should be used depends on the sweep count. The computational complexity of one sweep is proportional to the complexity of the original code. Forward accumulation is more efficient than reverse accumulation for functions f : \R^{n} \to \R^{m} with n \ll m as only n sweeps are necessary, compared to m sweeps for reverse accumulation. Reverse accumulation is more efficient than forward accumulation for functions f : \R^{n} \to \R^{m} with n \gg m as only m sweeps are necessary, compared to n sweeps for forward accumulation. Backpropagation of errors in multilayer perceptrons, a technique used in machine learning, is a special case of reverse accumulation. Forward accumulation was introduced by R. E. Wengert in 1964. According to Andreas Griewank, reverse accumulation has been suggested since the late 1960s, but the inventor is unknown. Seppo Linnainmaa published reverse accumulation in 1976¹. \begin{align} \frac{\partial y}{\partial x} &= \frac{\partial y}{\partial w_{n-1}} \frac{\partial w_{n-1}}{\partial x} \\[6pt] &= \frac{\partial y}{\partial w_{n-1}} \left(\frac{\partial w_{n-1}}{\partial w_{n-2}} \frac{\partial w_{n-2}}{\partial x}\right) \\[6pt] &= \frac{\partial y}{\partial w_{n-1}} \left(\frac{\partial w_{n-1}}{\partial w_{n-2}} \left(\frac{\partial w_{n-2}}{\partial w_{n-3}} \frac{\partial w_{n-3}}{\partial x}\right)\right) \\[6pt] &= \cdots \end{align} This can be generalized to multiple variables as a matrix product of Jacobians. Compared to reverse accumulation, forward accumulation is natural and easy to implement as the flow of derivative information coincides with the order of evaluation. Each variable w_i is augmented with its derivative \dot w_i (stored as a numerical value, not a symbolic expression), \dot w_i = \frac{\partial w_i}{\partial x} as denoted by the dot. The derivatives are then computed in sync with the evaluation steps and combined with other derivatives via the chain rule. Using the chain rule, if w_i has predecessors in the computational graph: \dot w_i = \sum_{j \in \{\text{predecessors of }i\}} \frac{\partial w_i}{\partial w_j} \dot w_j As an example, consider the function: \begin{align} y &= f(x_1, x_2) \\ &= x_1 x_2 + \sin x_1 \\ &= w_1 w_2 + \sin w_1 \\ &= w_3 + w_4 \\ &= w_5 \end{align} For clarity, the individual sub-expressions have been labeled with the variables w_i . The choice of the independent variable to which differentiation is performed affects the seed values and . Given interest in the derivative of this function with respect to , the seed values should be set to: \begin{align} \dot w_1 = \frac{\partial w_1}{\partial x_1} = \frac{\partial x_1}{\partial x_1} = 1 \\ \dot w_2 = \frac{\partial w_2}{\partial x_1} = \frac{\partial x_2}{\partial x_1} = 0 \end{align} With the seed values set, the values propagate using the chain rule as shown. Figure 2 shows a pictorial depiction of this process as a computational graph. To compute the gradient of this example function, which requires not only \frac{\partial y}{\partial x_1} but also \frac{\partial y}{\partial x_2} , an additional sweep is performed over the computational graph using the seed values \dot w_1 = 0 ; \dot w_2 = 1 . tuple evaluateAndDerive(Expression Z, Variable V) { if isVariable(Z) if (Z = V) return {valueOf(Z), 1}; else return {valueOf(Z), 0}; else if (Z = A + B) {a, a'} = evaluateAndDerive(A, V); {b, b'} = evaluateAndDerive(B, V); return {a + b, a' + b'}; else if (Z = A - B) {a, a'} = evaluateAndDerive(A, V); {b, b'} = evaluateAndDerive(B, V); return {a - b, a' - b'}; else if (Z = A * B) {a, a'} = evaluateAndDerive(A, V); {b, b'} = evaluateAndDerive(B, V); return {a * b, b * a' + a * b'}; } include struct ValueAndPartial { float value, partial; }; struct Variable; struct Expression { virtual ValueAndPartial evaluateAndDerive(Variable &variable) = 0; }; struct Variable: public Expression { float value; Variable(float value): value(value) {} ValueAndPartial evaluateAndDerive(Variable &variable) { float partial = (this == &variable) ? 1.0f : 0.0f; return {value, partial}; } }; struct Plus: public Expression { Expression &a, &b; Plus(Expression &a, Expression &b): a(a), b(b) {} ValueAndPartial evaluateAndDerive(Variable &variable) { auto [valueA, partialA] = a.evaluateAndDerive(variable); auto [valueB, partialB] = b.evaluateAndDerive(variable); return {valueA + valueB, partialA + partialB}; } }; struct Multiply: public Expression { Expression &a, &b; Multiply(Expression &a, Expression &b): a(a), b(b) {} ValueAndPartial evaluateAndDerive(Variable &variable) { auto [valueA, partialA] = a.evaluateAndDerive(variable); auto [valueB, partialB] = b.evaluateAndDerive(variable); return {valueA * valueB, valueB * partialA + valueA * partialB}; } }; int main () { // Example: Finding the partials of z = x * (x + y) + y * y at (x, y) = (2, 3) Variable x(2), y(3); Plus p1(x, y); Multiply m1(x, p1); Multiply m2(y, y); Plus z(m1, m2); float xPartial = z.evaluateAndDerive(x).partial; float yPartial = z.evaluateAndDerive(y).partial; std::cout \begin{align} \frac{\partial y}{\partial x} &= \frac{\partial y}{\partial w_1} \frac{\partial w_1}{\partial x}\\[6px] &= \left(\frac{\partial y}{\partial w_2} \frac{\partial w_2}{\partial w_1}\right) \frac{\partial w_1}{\partial x}\\[6px] &= \left(\left(\frac{\partial y}{\partial w_3} \frac{\partial w_3}{\partial w_2}\right) \frac{\partial w_2}{\partial w_1}\right) \frac{\partial w_1}{\partial x}\\[6px] &= \cdots \end{align} In reverse accumulation, the quantity of interest is the adjoint, denoted with a bar \bar w_i ; it is a derivative of a chosen dependent variable with respect to a subexpression w_i : \bar w_i = \frac{\partial y}{\partial w_i} Using the chain rule, if w_i has successors in the computational graph: \bar w_i = \sum_{j \in \{\text{successors of }i\}} \bar w_j \frac{\partial w_j}{\partial w_i} Reverse accumulation traverses the chain rule from outside to inside, or in the case of the computational graph in Figure 3, from top to bottom. The example function is scalar-valued, and thus there is only one seed for the derivative computation, and only one sweep of the computational graph is needed to calculate the (two-component) gradient. This is only half the work when compared to forward accumulation, but reverse accumulation requires the storage of the intermediate variables as well as the instructions that produced them in a data structure known as a "tape" or a Wengert list (however, Wengert published forward accumulation, not reverse accumulation ), which may consume significant memory if the computational graph is large. This can be mitigated to some extent by storing only a subset of the intermediate variables and then reconstructing the necessary work variables by repeating the evaluations, a technique known as rematerialization. Checkpointing is also used to save intermediary states. We consider again the example function y=f(x_1,x_2)=x_1 x_2 + \sin x_1 with subexpressions w_1=x_1, w_2=x_2, w_3=w_1 w_2, w_4=\sin(w_1), w_5=w_3+w_4 . The operations to compute the adjoint values using reverse accumulation are shown in the table below (note the reversed order): \begin{align} \bar w_5 &= 1 \quad\text{(seed)} \\ \bar w_4 &= \bar w_5 \cdot 1 \\ \bar w_3 &= \bar w_5 \cdot 1 \\ \bar w_2 &= \bar w_3 \cdot w_1 \\ \bar w_1 &= \bar w_3 \cdot w_2 + \bar w_4 \cdot \cos w_1 \end{align} The sought-after partial derivatives are then \partial y/\partial x_1 = \bar{w}_1 and \partial y/\partial x_2 = \bar{w}_2 . The data flow graph of a computation can be manipulated to calculate the gradient of its original calculation. This is done by adding an adjoint node for each primal node, connected by adjoint edges which parallel the primal edges but flow in the opposite direction. The nodes in the adjoint graph represent multiplication by the derivatives of the functions calculated by the nodes in the primal. For instance, addition in the primal causes fanout in the adjoint; fanout in the primal causes addition in the adjoint; a unary function in the primal causes in the adjoint; etc. void derive(Expression Z, float seed) { if isVariable(Z) partialDerivativeOf(Z) += seed; else if (Z = A + B) derive(A, seed); derive(B, seed); else if (Z = A - B) derive(A, seed); derive(B, -seed); else if (Z = A * B) derive(A, valueOf(B) * seed); derive(B, valueOf(A) * seed); } include struct Expression { float value; virtual void evaluate() = 0; virtual void derive(float seed) = 0; }; struct Variable: public Expression { float partial; Variable(float value) { this->value = value; partial = 0.0f; } void evaluate() {} void derive(float seed) { partial += seed; } }; struct Plus: public Expression { Expression &a, &b; Plus(Expression &a, Expression &b): a(a), b(b) {} void evaluate() { a.evaluate(); b.evaluate(); value = a.value + b.value; } void derive(float seed) { a.derive(seed); b.derive(seed); } }; struct Multiply: public Expression { Expression &a, &b; Multiply(Expression &a, Expression &b): a(a), b(b) {} void evaluate() { a.evaluate(); b.evaluate(); value = a.value * b.value; } void derive(float seed) { a.derive(b.value * seed); b.derive(a.value * seed); } }; int main () { // Example: Finding the partials of z = x * (x + y) + y * y at (x, y) = (2, 3) Variable x(2), y(3); Plus p1(x, y); Multiply m1(x, p1); Multiply m2(y, y); Plus z(m1, m2); z.evaluate(); std::cout Replace every number \,x with the number x + x'\varepsilon , where x' is a real number, but \varepsilon is an abstract number with the property \varepsilon^2=0 (an infinitesimal; see Smooth infinitesimal analysis). Using only this, regular arithmetic gives \begin{align} (x + x'\varepsilon) + (y + y'\varepsilon) &= x + y + (x' + y')\varepsilon \\[4px] (x + x'\varepsilon) - (y + y'\varepsilon) &= x - y + (x' - y')\varepsilon \\[4px] (x + x'\varepsilon) \cdot (y + y'\varepsilon) &= xy + xy'\varepsilon + yx'\varepsilon + x'y'\varepsilon^2 = xy + (x y' + yx')\varepsilon \\[8px] \frac{x + x'\varepsilon}{y + y'\varepsilon} &= \frac{\frac{x}{y} + \frac{x'\varepsilon}{y}} {1 + \frac{y'\varepsilon}{y}} = \left(\frac{x}{y} + \frac{x'\varepsilon}{y}\right) \cdot \left(1 - \frac{y'\varepsilon}{y}\right) = \frac{x}{y} + \left(\frac{x'y-xy'}{y^2} \right)\varepsilon \end{align} using the fact that \left(1 + \frac{y'\varepsilon}{y}\right) \cdot \left(1 - \frac{y'\varepsilon}{y}\right) = 1. Now, polynomials can be calculated in this augmented arithmetic. If P(x) = p_0 + p_1 x + p_2x^2 + \cdots + p_n x^n, then \begin{align} P(x + x'\varepsilon) &= p_0 + p_1(x + x'\varepsilon) + \cdots + p_n (x + x'\varepsilon)^n \\ &= p_0 + p_1 x + \cdots + p_n x^n + p_1x'\varepsilon + 2p_2xx'\varepsilon + \cdots + np_n x^{n-1} x'\varepsilon \\ &= P(x) + P^{(1)}(x)x'\varepsilon \end{align} where P^{(1)} denotes the derivative of P with respect to its first argument, and x' , called a seed, can be chosen arbitrarily. The new arithmetic consists of ordered pairs, elements written \langle x, x' \rangle , with ordinary arithmetics on the first component, and first order differentiation arithmetic on the second component, as described above. Extending the above results on polynomials to analytic functions gives a list of the basic arithmetic and some standard functions for the new arithmetic: \begin{align} \left\langle u,u'\right\rangle + \left\langle v,v'\right\rangle &= \left\langle u + v, u' + v' \right\rangle \\[4px] \left\langle u,u'\right\rangle - \left\langle v,v'\right\rangle &= \left\langle u - v, u' - v' \right\rangle \\[4px] \left\langle u,u'\right\rangle \cdot \left\langle v,v'\right\rangle &= \left\langle u v, u'v + uv' \right\rangle \\[8px] \frac{\left\langle u,u'\right\rangle}{\left\langle v,v'\right\rangle} &= \left\langle \frac{u}{v}, \frac{u'v - uv'}{v^2} \right\rangle && ( v\ne 0) \\[8px] \sin\left\langle u,u'\right\rangle &= \left\langle \sin(u) , u' \cos(u) \right\rangle \\[4px] \cos\left\langle u,u'\right\rangle &= \left\langle \cos(u) , -u' \sin(u) \right\rangle \\[4px] e^{\left\langle u,u'\right\rangle} &= \left\langle e^u , u' e^u \right\rangle \\[8px] \log\left\langle u,u'\right\rangle &= \left\langle \log(u) , \frac{u'}{u} \right\rangle && (u>0) \\[8px] \left\langle u,u'\right\rangle^k &= \left\langle u^k , u' k u^{k - 1} \right\rangle && (u \ne 0) \\[8px] \left| \left\langle u,u'\right\rangle \right| &= \left\langle \left| u \right| , u' \operatorname{sgn} u \right\rangle && (u \ne 0) \end{align} and in general for the primitive function g , g(\langle u,u' \rangle , \langle v,v' \rangle ) = \langle g(u,v) , g_u(u,v) u' + g_v(u,v) v' \rangle where g_u and g_v are the derivatives of g with respect to its first and second arguments, respectively. When a binary basic arithmetic operation is applied to mixed arguments—the pair \langle u, u' \rangle and the real number c —the real number is first lifted to \langle c, 0 \rangle . The derivative of a function f : \R\to\R at the point x_0 is now found by calculating f(\langle x_0, 1 \rangle) using the above arithmetic, which gives \langle f ( x_0 ) , f' ( x_0 ) \rangle as the result. {{pre|1= Dual plus(Dual A, Dual B) { return { realPartOf(A) + realPartOf(B), infinitesimalPartOf(A) + infinitesimalPartOf(B) }; } Dual minus(Dual A, Dual B) { return { realPartOf(A) - realPartOf(B), infinitesimalPartOf(A) - infinitesimalPartOf(B) }; } Dual multiply(Dual A, Dual B) { return { realPartOf(A) * realPartOf(B), realPartOf(B) * infinitesimalPartOf(A) + realPartOf(A) * infinitesimalPartOf(B) }; } X = {x, 0}; Y = {y, 0}; Epsilon = {0, 1}; xPartial = infinitesimalPartOf(f(X + Epsilon, Y)); yPartial = infinitesimalPartOf(f(X, Y + Epsilon)); }} include struct Dual { float realPart, infinitesimalPart; Dual(float realPart, float infinitesimalPart=0): realPart(realPart), infinitesimalPart(infinitesimalPart) {} Dual operator+(Dual other) { return Dual( realPart + other.realPart, infinitesimalPart + other.infinitesimalPart ); } Dual operator*(Dual other) { return Dual( realPart * other.realPart, other.realPart * infinitesimalPart + realPart * other.infinitesimalPart ); } }; // Example: Finding the partials of z = x * (x + y) + y * y at (x, y) = (2, 3) Dual f(Dual x, Dual y) { return x * (x + y) + y * y; } int main () { Dual x = Dual(2); Dual y = Dual(3); Dual epsilon = Dual(0, 1); Dual a = f(x + epsilon, y); Dual b = f(x, y + epsilon); std::cout ===Vector arguments and functions=== Multivariate functions can be handled with the same efficiency and mechanisms as univariate functions by adopting a directional derivative operator. That is, if it is sufficient to compute y' = \nabla f(x)\cdot x' , the directional derivative y' \in \R^m of f:\R^n\to\R^m at x \in \R^n in the direction x' \in \R^n may be calculated as (\langle y_1,y'_1\rangle, \ldots, \langle y_m,y'_m\rangle) = f(\langle x_1,x'_1\rangle, \ldots, \langle x_n,x'_n\rangle) using the same arithmetic as above. If all the elements of \nabla f are desired, then n function evaluations are required. Note that in many optimization applications, the directional derivative is indeed sufficient. ===High order and many variables=== The above arithmetic can be generalized to calculate second order and higher derivatives of multivariate functions. However, the arithmetic rules quickly grow complicated: complexity is quadratic in the highest derivative degree. Instead, truncated Taylor polynomial algebra can be used. The resulting arithmetic, defined on generalized dual numbers, allows efficient computation using functions as if they were a data type. Once the Taylor polynomial of a function is known, the derivatives are easily extracted. ==See also== Differentiable programming ==Notes== ==References== ==External links== www.autodiff.org, An "entry site to everything you want to know about automatic differentiation" Automatic Differentiation of Parallel OpenMP Programs Automatic Differentiation, C++ Templates and Photogrammetry Automatic Differentiation, Operator Overloading Approach Compute analytic derivatives of any Fortran77, Fortran95, or C program through a web-based interface Automatic Differentiation of Fortran programs Description and example code for forward Automatic Differentiation in Scala finmath-lib stochastic automatic differentiation, Automatic differentiation for random variables (Java implementation of the stochastic automatic differentiation). Adjoint Algorithmic Differentiation: Calibration and Implicit Function Theorem C++ Template-based automatic differentiation article and implementation Tangent Source-to-Source Debuggable Derivatives Exact First- and Second-Order Greeks by Algorithmic Differentiation Adjoint Algorithmic Differentiation of a GPU Accelerated Application Adjoint Methods in Computational Finance Software Tool Support for Algorithmic Differentiationop More than a Thousand Fold Speed Up for xVA Pricing Calculations with Intel Xeon Scalable Processors Sparse truncated Taylor series implementation with VBIC95 example for higher order derivatives

Text derived from https://en.wikipedia.org/w/index.php?oldid=1373407057 under CC BY-SA 4.0; modifications have been made.