Numerical Jacobian matrix

Is there a simple way to evaluate the numerical Jacobian of a "black-box" function (defined only for lists of numerical values, for example)?

I the meanwhile I am using the following based on ND:

Needs["NumericalCalculus"]
jac[f_, x0_] := Table[xtmp = x0; xtmp[[i]] = x; ND[f[xtmp], x, 1], {i, Length@x0}]


Here is a very dirty trick that is somewhat related to automatic differentiation (note that automatic differentiation $\neq$ automated differentiation). It works for twice differentiable, real-valued functions $f$ that

• are implemented for complex arguments but
• do not involve Re and Im and such things.

Examples are polynomials and everything that is numerically approximated by Taylor expansions or Padé approximations. That's quite a lot but not everything, of course.

The strength of the method is that it needs only a single evaluation of the function per directional derivative. The idea is to emulate a monade (in the sense of Leibniz) by using $\varepsilon = \operatorname{i}\, 10^{-8}$ as this squares approximately to $0$ in machine precision. Here is the code:

monade = I 1. 10^-8
nD[f_, x : {__?NumericQ}] :=
f,
]]


A short test:

f[{x_?NumericQ, y_?NumericQ}] := Sin[x y];
nD[f, {1., 2.}] - (D[Sin[x y], {{x, y}, 1}] /. {x -> 1., y -> 2.})


{-1.11022*10^-16, -5.55112*10^-17}

The following Taylor expansion shows why it works

nDWith1Point[f_, x_] := 1/h (f[x + I h])
ComplexExpand[Im@Normal[Series[nDWith1Point[f, x] , {h, 0, 4}]]]


Experimenting a bit, one can find a 4-point method that has even 4th order accuracy:

nDWith4Point[f_, x_] := 1/(4 h) (f[x + h] - I f[x + I h] + I f[x - I h] - f[x - h])
Normal[Series[nDWith4Point[f, x] , {h, 0, 9}]]


The main advantages: 1) The monade can be chosen much larger for this, e.g. $10^{-4}$ such that this is much less vulnerable to cancellation errors. 2) It can also be applied to mixed derivatives, leading to a 16-point rule for second derivatives.

I am sure that @J.M. even knows how this differential quadrature is called in the literature...

Oh, and see what we can do, a $2^k$-point rule of order $2^k$!

nDWith2PowerKPoints[f_, x_, k_, h_] := With[{e = Table[Exp[2 Pi I/2^k j], {j, 0, 2^k - 1}]},
1/(2^k h) Total[Times[1/e, f /@ (x + h e)]]
];
k = 5;
Simplify[Series[nDWith2PowerKPoints[f, x, k, h], {h, 0, 2^k}]]


Actually, this works also with other $N$, not only for $N = 2^k$, but Mathematica does not simplify the above expression enough.

Look how precise it is for h = 0.1:

x = 1.100;
h = 0.1100;
Re[nDWith2PowerKPoints[Sin, x, 5, h] - Sin'[x]]


6.22230269105447896725778959568*10^-70

Edit

I found out that the reason is the following formula which holds for every analytic function (and in a truncated variant if the function has at least sufficiently many derivatives):

$$\frac{1}{N \, h} \sum_{k=1}^{N} f \Big( z + h \, \operatorname{e}^{2 \pi \operatorname{i} \tfrac{k}{N}} \Big) \, \operatorname{e}^{- 2 \pi\operatorname{i}\tfrac{k}{N}} = \sum_{j = 0}^\infty \frac{f^{(N j +1)}(z)}{(N j +1)!} \, h^{Nj}.$$

• What Henrik has done here is called "complex step differentiation", and is valuable if your black-box function can be evaluated at complex arguments. See this and this for more information. – J. M. will be back soon Nov 10 '17 at 12:05
• You flatter me, Henrik :P but the new formula you have does resemble the extensions studied by Lai and Crassidis; see e.g. this. – J. M. will be back soon Nov 10 '17 at 13:21
• @J.M. Hehe. Thanks for the pointer to that paper. The next step I took is precisely into its direction. Look, I just found a $2^k$-point rule of order $2^k$, completely by chance (added in answer). And of course that pretty much looks like a Cauchy integral... – Henrik Schumacher Nov 10 '17 at 13:42
• ...and as the paper itself notes, you can even do Richardson extrapolation if you're for some reason not satisfied with the complex step evaluations themselves. – J. M. will be back soon Nov 10 '17 at 13:56

Update

This update corrects my confusion in conflating the Jacobian with the Hessian matrix.

Derivative already has support for numerical differentiation of black box functions (although it doesn't work for mixed derivatives). First, here's a function to compute the numerical gradient:

grad[f_, x:{__?NumericQ}] := With[{d = Through @ (Derivative@@@IdentityMatrix[Length[x]]) @ f},
Through @ Apply[d] @ x
]


Here is an example using grad:

Clear[f]
f[x_?NumericQ, y_?NumericQ] := Sin[x y]



{-0.832294, -0.416147}

Let's check this with the usual definition:

D[Sin[x y], {{x, y}}] /. {x->1., y->2.}


{-0.832294, -0.416147}

The nice thing about relying on the numerical differentiation support of Derivative is that it works with extended precision numbers as well. Here is an example:

grad[f, {130, 230}]
D[Sin[x y], {{x, y}}] /. {x->130, y->230}


{-0.83229367309428477399513646, -0.41614683654714238699756823}

{-0.83229367309428477399513645900, -0.41614683654714238699756822950}

Very good agreement! (Note however, as pointed out by @HenrikSchumacher, Mathematica does use quite a lot of black box evaluations of the function to arrive at this answer!).

Hessian

(My previous revision unnecessarily provided an implementation of the Hessian, so I won't remove it).

To compute the Hessian we will need to deal with mixed derivatives, e.g., Derivative[1, 1][f][1., 2.]. Unfortunately, while:

Derivative[1, 0][f][1., 2.]


-0.832294

works fine, the mixed version does not:

Derivative[1, 1][f][1., 2.]


(f^(1,1))[1.,2.]

To deal with this, the following function creates a set of intermediate gradient functions, which it then differentiates:

hess[f_, x_:{__?NumericQ}] := Module[{g, len=Length[x], a, h},
(#[a__?NumericQ] := (Derivative@@#2)[f][a])&,
{h, IdentityMatrix@len}
];
]


Here is an example:

hess[f, {1., 2.}]
D[Sin[x y], {{x, y}, 2}] /. {x->1., y->2.}


{{-3.63719, -2.23474}, {-2.23474, -0.909297}}

{{-3.63719, -2.23474}, {-2.23474, -0.909297}}

The same example at higher precision:

hess[f, {130, 230}]
D[Sin[x y], {{x, y}, 2}] /. {x->130, y->230}


{{-3.637189707302726781584079, -2.234741690198505777789608}, \ {-2.234741690198505777789608, -0.909297426825681695396020}}

{{-3.6371897073027267815840794636, -2.2347416901985057777896079613}, \ {-2.2347416901985057777896079613, -0.90929742682568169539601986591}}

• Could it be that "Jacobian" and "Hessian" get mixed up here? Very confusing... Jacobians are first derivatives of vector valued functions while Hessians are (essentially) second derivatives of scalar valued functions. For each scalar valued function, the Jacobian of the gradient equals the Hessian. – Henrik Schumacher Nov 10 '17 at 8:52
• Moreover, the following code snippet tells me, that even the simple function f in two variables from your example gets evaluated 102 times(!) to obtain a numerical derivate. That's quite often, isn't it? Here is the code: Clear[g];c = 0;g[x_?NumericQ, y_?NumericQ] := (c++; Sin[x y]);grad[g, {1., 2.}];c`. – Henrik Schumacher Nov 10 '17 at 9:02
• @HenrikSchumacher Thanks for the correction! – Carl Woll Nov 15 '17 at 23:07
• You're welcome! – Henrik Schumacher Nov 15 '17 at 23:31