How to do automatic differentiation?

Question

Does Mathematica have any AD functionality or does it only support symbolic differentiation? If not, are there any packages or other third party implementations available?

---

Edit (J. M.)

Here is an example that might clarify how AD works. Consider the Horner method for evaluating a polynomial from its coefficients:

h[coeffs_?VectorQ, x_] := Module[{np1 = Length[coeffs], p, k},
   p = 0;
   Do[p = x p + coeffs[[k]], {k, np1, 1, -1}];
   p
   ]

An AD method should be able to "differentiate" this procedure and automatically produce this:

hp[coeffs_?VectorQ, x_] := Module[{np1 = Length[coeffs], p, p$, k},
       p$ = 0; p = 0;
   Do[p$ = x p$ + p;
      p = x p + coeffs[[k]], {k, np1, 1, -1}];
   p$
   ]

Note the use of the product rule to differentiate the line p = x p + coeffs[[k]], as well as the ordering of the variable assignments in the "differentiated procedure".

There are more complicated examples, but traditionally, AD methods have dealt with routines that have been written procedurally. A challenge, then, is if AD methods can cope with more "idiomatic" Mathematica programs; that is, functional/rule-based algorithms. I believe that will be more difficult.

Answer 1

This is not a package, just an idea how such a package might work. You could define a custom value type autoDiffValue that stores a value and it's first derivative (with respect to some variable), then define arithmetic operations for that type:

Clear[autoDiffValue]
autoDiffValue /: autoDiffValue[a_, adx_] + autoDiffValue[b_, bdx_] := 
 autoDiffValue[a + b, adx + bdx]
autoDiffValue /: autoDiffValue[a_, adx_] * autoDiffValue[b_, bdx_] := 
 autoDiffValue[a*b, adx*b + a*bdx]
autoDiffValue /: autoDiffValue[a_, adx_] + b_ := 
 autoDiffValue[a + b, adx]
autoDiffValue /: autoDiffValue[a_, adx_]*b_ := 
 autoDiffValue[a*b, adx*b]
autoDiffValue /: autoDiffValue[a_, adx_]^b_ := 
 autoDiffValue[a^b, a^(b - 1)*adx*b]

Using that and your function h, you can then call:

h[{1, 2, 3}, autoDiffValue[5, 1]]

and get

autoDiffValue[86, 32]

i.e. the value and the first derivative of h[{1, 2, 3}, x] at x=5

You could automate generating derivative definitions:

Clear[defineFunc]
defineFunc[fn_] := 
 autoDiffValue /: fn[autoDiffValue[a_, adx_]] = 
  autoDiffValue[fn[a], D[fn[a + x*adx], x] /. x -> 0]
defineFunc /@ {Sin, Cos, Exp}; 

(it's probably possible to extend defineFunc so it can create definitions for Plus, Times, ... too)

This is more or less the same technique you would use to define AD in e.g. C++ using operator overloading. It should work in principle, although a lot of work would still be needed to implement multivariable and higher derivatives and to support vector and matrix arguments (so calculations can be done with packed arrays).

Answer 2

I am not sure I understood the advantages of AD from your example. It would be nice to have a function that you can't differentiate properly or fast with Derivative.

This is not an answer, but too long for a comment. I suppose the automatic accuracy tracking is some sort of very limited built-in automatic differentiation, right? At least for combinations of real monotonically increasing functions.

makeDual[n_] := SetAccuracy[n, 1]
getDual[n_] := {N@n, 10^(1 - Accuracy[n])}

f[x_] := If[x > 0, Log[Sin[x]] + Tan[x]^3 - Sin[Cos[Sin[x]]]^17, x - 3]
r = 0.234;
{f[r], f'[r]} == getDual@f[makeDual[r]]
(* True *)