Automatic Differentiation - Kenjudd

AUTOMATIC DIFFERENTIATIONPHILIPP MÜLLERUNIVERSITY OF ZURICH

MOTIVATIONDerivatives are omnipresent in numerical algorithms1st order derivatives Solving non-linear equations E.g., by Newton’s method(Un-)constrained optimization Gradient-Based optimization algorithms Especially difficult for high dimensional variables, i.e., objective function 𝑓: 𝑅 𝑛 𝑅 Structural sparsity can be key2nd order derivatives (Un-)constrained optimizationHigher order derivatives Higher-order differential equations

MOTIVATIONSuppose we want to solve the unconstrained optimization problemmin 𝑓(𝑥)𝑥with 𝑓: 𝑅 𝑅 and 𝑥 𝑅.Gradient-based optimization requires the gradient 𝑓 𝑥

FINITE DIFFERENCESRecall: The Taylor series expansion of a real-valued function 𝑓 𝒞 𝑛 , 𝑛 2, around 𝑥 and evaluated at 𝑎 reads𝑓 𝑎 𝛼 𝑛1 𝛼 𝑓 𝑥 𝑎 𝑥𝛼! 𝑓 𝑥 𝑓 𝑥𝑥 𝑎 𝑥 𝛼 𝑅1 2 𝑓2! 𝑥 2 O( a x 3 )𝑥 𝑎 𝑥2

TAYLOR APPROXIMATION OF 𝑒 𝑥 AROUND 0𝑓 𝑥 𝑒𝑥

TAYLOR APPROXIMATION OF 𝑒 𝑥 AROUND 0𝑥2 𝑥3 𝑥4𝒯𝑓 𝑥 1 𝑥 26 24

TAYLOR APPROXIMATION OF 𝑒 𝑥 AROUND 0𝑥2 𝑥3 𝑥4𝒯𝑓 𝑥 1 𝑥 26 24

FINITE DIFFERENCESRecall (Taylor Series):𝑓 𝑎 𝑓 𝑥 𝑓 𝑥𝑥 𝑎 𝑥 1 2 𝑓2! 𝑥 2𝑥 𝑎 𝑥2Truncate the Taylor series and set a 𝑥 ℎ with a small ℎ yields: 𝑓𝑓 𝑥 ℎ 𝑓 𝑥 𝑥 𝑥 ℎ 𝑥 𝑂 ℎ2 𝑥 𝑓𝑓 𝑥 ℎ 𝑓 𝑥֞𝑥 𝑂(ℎ) 𝑥ℎThis results in the well-known forward difference equation: 𝑓𝑓 𝑥 ℎ 𝑓 𝑥𝑥 𝑥ℎ

WHY WOULD WE NEED ANYTHING ELSE? 𝑓We derivedby truncating the Taylor series resulting in the error 𝑂 ℎ . This truncation error decreases in the 𝑥step size. 𝑓𝑓 𝑥 ℎ 𝑓 𝑥ℎℎ 0Accurate and efficient approximation of 𝑥 by choosing a very small h, i.e., lim?

PROBLEM SOLVED?Apply forward differences to𝑓 𝑥 𝑥3and increase the step size from 10 16 to 0.1. 𝑓 3𝑥 2 𝑥 𝐹𝐷ℎ

PROBLEM SOLVED?Apply forward differences to𝑓 𝑥 𝑥3and increase the step size from 10 16 to 0.1.CD: 𝑓 𝑥1 1𝑓 𝑥 2ℎ 𝑓 𝑥 2ℎℎ 𝑓 3𝑥 2 𝑥 𝐹𝐷ℎ

NUMERICAL ERRORS IN FINITE PRECISION ARITHMETICRounding error Intermediate results are rounded Any rounding error propagates and amplifiesTruncation error Even if an algorithm is converging to the true solution, we are stopping it after some finite time. Mitigated by the appropriate convergence criteria as introduced by Ken.Cancellation error

CANCELLATION ERROR𝑎 ℎ 𝑛𝑜𝑖𝑠𝑒𝑛𝑜𝑖𝑠𝑒1 ℎℎ 𝑛𝑜𝑖𝑠𝑒𝑎h 1E-12a 1-True valueNoise /h b a h# add h to ac b–a# c should be equal to hd c/h# c h, thus d should 1d 0.999200722162641

FINITE DIFFERENCES FOR 𝑓: 𝑅𝑛 𝑅Let’s consider a n-dimensional unconstrained optimization problemmin 𝑓(𝑥)𝑥with 𝑓: 𝑅𝑛 𝑅 and 𝑥 𝑅𝑛 .The finite difference quotients resemble directional derivatives for 𝑓: 𝑅𝑛 𝑅𝑚 and n 1 : 𝑓𝑓 𝑥 ℎ𝒆𝒊 𝑓 𝑥𝑥 . 𝑥𝑖ℎThe cost of FD scales with n - O(n) * cost(f).

SCALING𝑛 Let’s consider the Rosenbrock function 𝑓: 𝑅 𝑛 𝑅 asbenchmark𝑛 1𝑓 𝑥 10 𝑥𝑖 1 2𝑥𝑖2 1 𝑥𝑖𝑖 1 The runtimes are averaged across 1000 runs.2.𝑱𝑭𝑫 [s]𝒇 [s]106.7931e-05 (54x)1.2472e-061006.0959e-04 (332x)1.8355e-0610001.4839e-02 (1500x)9.9629e-06100001.3282 (20000x)6.6663e-05100000?9.0872e-04

AUTOMATIC DIFFERENTIATIONBasic Idea: Every computer program is a composition of differentiable elementary operations as, basic arithmetic operations as, e.g., , -, and *, and basic functions as, e.g., sin, cos and tan.Automatic differentiation can transform the source code of your function into the source code of the gradient.

TOY EXAMPLEConsider the function 𝑓: 𝑅2 𝑅𝑓 𝑥1 , 𝑥2 𝑥1 𝑥2 sin 𝑥1This function can be discomposedin differentiable elementaryoperations:𝑤1 𝑥1𝑤2 𝑥2𝑤3 𝑤1 𝑤2𝑤4 sin(𝑤1 )𝑤5 𝑤3 𝑤4𝑓 𝑤5𝑥1𝑤1sin𝑤4 𝑥2𝑤2 𝑤3𝑤5𝑓

FORWARD MODEConsider the function 𝑓: 𝑅2 𝑅𝑓 𝑥1 , 𝑥2 𝑥1 𝑥2 sin 𝑥1To calculate the Gradient, calculate 𝑓 𝑥1 , 𝑥2 𝑥1Choose input variable 𝑥1 andcalculate the sensitivity of eachintermediate value as 𝑤𝑖𝑤ሶ 𝑖 𝑥𝑗 𝑤4 𝑤1𝑤ሶ 4 cos(𝑤1 )𝑤ሶ 1 𝑤1 𝑥𝑗 𝑤1𝑤ሶ 1 𝑥𝑗𝑥1𝑤1sin𝑤4𝑤ሶ 5 𝑤ሶ 3 𝑤ሶ 4 𝑥2𝑤2 𝑤2𝑤ሶ 2 𝑥𝑗 𝑤3𝑤ሶ 3 𝑤ሶ 1 𝑤2 𝑤1 𝑤ሶ 2𝑤5𝑓𝑓𝑥𝑗 𝑤ሶ 3 𝑤ሶ 4

FORWARD MODE: EVALUATIONSuppose: 𝑗 1, 𝑥1 1, 𝑥2 2Consider the function 𝑓: 𝑅2 𝑅𝑓 𝑥1 , 𝑥2 𝑥1 𝑥2 sin 𝑥1To calculate the Gradient, calculate 𝑓 𝑥1 , 𝑥2 𝑥1Choose input variable 𝑥1 andcalculate the sensitivity of eachintermediate value as 𝑤𝑖𝑤ሶ 𝑖 𝑥1𝑤ሶ 1 𝑥1 𝑤1 1 𝑥1𝑤1𝑤ሶ 4 cos 𝑤1 𝑤ሶ 1 cos(1)sin𝑤4𝑤ሶ 5 𝑤ሶ 3 𝑤ሶ 4 2 cos(1) 𝑥2𝑤2𝑤ሶ 2 𝑤2 0 𝑥1 𝑤3𝑤ሶ 3 𝑤ሶ 1 𝑤2 𝑤1 𝑤ሶ 2 2Accurate up to working precision, but still scales linearly in n. 𝐶𝑜𝑠𝑡 𝐽𝑓 𝑛 𝑐 𝐶𝑜𝑠𝑡(𝑓)𝑤5𝑓𝑓𝑥1 2 cos(1)

REVERSE MODE (ADJOINT MODE) – PRIMAL TRACESuppose: 𝑥1 1, 𝑥2 2Consider the function 𝑓: 𝑅2 𝑅𝑓 𝑥1 , 𝑥2 𝑥1 𝑥2 sin 𝑥1𝑤4 sin 1 0.175𝑤1 1𝑥1𝑤1sin𝑤4 Calculate the sensitivity of the outputw.r.t. each intermediate value 𝑓𝑤ഥ𝑖 𝑤𝑗All intermediate values are stored. Thisleads to a high memory consumption,mitigated by good AD software𝑥2𝑤2𝑤2 2 𝑤3𝑤3 2𝑤5𝑤5 2.175𝑓

REVERSE MODE (ADJOINT MODE) 𝑓 𝑤4 𝑓 𝑤3 𝑤4 𝑤1 𝑤3 𝑤1 𝑤ഥ4 cos 𝑤1 𝑤ഥ3 𝑤2𝑤ഥ1 Consider the function 𝑓: 𝑅2 𝑅𝑓 𝑥1 , 𝑥2 𝑥1 𝑥2 sin 𝑥1𝑥1𝑤1sin 𝑓 𝑓 𝑤5𝑤ഥ4 𝑤ഥ5 1 𝑤4 𝑤5 𝑤4𝑤4 Calculate the sensitivity of the outputw.r.t. each intermediate value 𝑓𝑤ഥ𝑖 𝑤𝑗𝑥2𝑤2𝑤ഥ2 𝑓 𝑤3 𝑤ഥ3 𝑤1 𝑤3 𝑤2𝑤3𝑤ഥ3 𝑓𝑤5𝑤ഥ5 𝑓 𝑤5 𝑤ഥ5 1 𝑤5 𝑤3 𝑓 𝑤5

REVERSE MODE (ADJOINT MODE) – DUAL TRACESuppose: 𝑥1 1, 𝑥2 2Consider the function 𝑓: 𝑅2 𝑅𝑓 𝑥1 , 𝑥2 𝑥1 𝑥2 sin 𝑥1𝑤1 1𝑤ഥ1 𝑤ഥ4 cos 𝑤1 𝑤ഥ3 𝑤2 cos 1 2𝑥1𝑤1sin𝑤4 sin 1 0.175𝑤ഥ4 𝑤ഥ5 1𝑤4 Calculate the sensitivity of the outputw.r.t. each intermediate value 𝑓𝑤ഥ𝑖 𝑤𝑗𝑥2𝑤2𝑤2 2𝑤ഥ2 1 𝑤3𝑤3 2𝑤ഥ3 𝑤ഥ5 1Accurate up to working precision, scales linearly in m. 𝐶𝑜𝑠𝑡 𝐽𝑓 𝑚 𝑐2 𝐶𝑜𝑠𝑡(𝑓)𝑓𝑤5𝑤5 2.175 𝑓𝑤ഥ5 𝑤 15

SUMMARYFinite differences The approximation error decreases as 𝑂 ℎ for forward finite differences. BUT, the error due to the finite precision arithmetic cannot be neglected. The time required to compute the Jacobian of 𝑓: 𝑅 𝑛 𝑅 𝑚 scales with 𝑂 𝑛 𝑐𝑜𝑠𝑡 𝑓 .AD - Forward mode The gradients are accurate up to machine precision. The time required to compute the Jacobian of 𝑓: 𝑅 𝑛 𝑅 𝑚 scales with 𝑂 𝑛 𝑐𝑜𝑠𝑡 𝑓AD - Reverse mode The gradients are accurate up to machine precision. The memory requirement may be huge depending on the underlyingimplementation. The time required to compute the Jacobian of 𝑓: 𝑅 𝑛 𝑅 𝑚 scales with 𝑂 𝑚 𝑐𝑜𝑠𝑡 𝑓

AD TOOLSCasADi Available for Python, Matlab, Octave and C Includes interfaces to a lot of free as well as commercial optimizers (as e.g., IPOPT (IP), KNITRO (IP & SQP), WORHP (SQP),SNOPT (SQP)) Structural sparsity detectionADiMat Available for MatlabPyTorch / Tensorflow

TUTORIAL SESSIONImplementation of the Rosenbrock function1.𝑛 1𝑓 𝑥 10 𝑥𝑖 1 𝑥𝑖22 1 𝑥𝑖2,𝑖 12.Implementation of the finite difference approximation and reverse mode AD of f(x). Comparison of theirruntimes for 𝑛 10𝑖 , i 1, 2, 3, 4, 3.Optimization of the Rosenbrock function using fminunc using1.the finite difference approximation of f(x), and2.the reverse mode approximation of f(x).

CASADI Include the casadi directory in the Matlab pathimport casadi.*x MX MX.sym(‘some name’, size rows, size columns)% create symbolic variabled rosenbrock jacobian(rosenbrock(x MX), x MX)% differentiate rosenbrockd rosenbrock Function(‘some name’, {x MX}, {d rosenbrock })d rosenbrock(x)% create callable function

AUTOMATIC DIFFERENTIATION Basic Idea: Every computer program is a composition of differentiable elementary operations as, basic arithmetic operations as, e.g., , -, and *, and basic functions as, e.g., sin, cos and tan. Automatic differentiation can transform the source code of your function into the source code of the gradient.