Caclulate second partial derivatives of a numerically defined function

Question

I am performing a chi square minimisation with respect to some parameters x,y,z that are themselves part of the kernel involved in a numeric integral. I therefore have a numerical function of three variables, say $\chi( x,y,z)$, which is only defined when $x,y,z$ are all numerical. As a test function, say we have

testfn[x_?NumericQ, y_?NumericQ, z_?NumericQ] := Sin[x y] + Cos[x z]

(obviously my actual function is much more complicated).

I've optimised my function numerically and found a solution set $(x_0,y_0,z_0)$, and want now I want to compute a correlation matrix in the space of parameters $\left\{x,y,z\right\}$ using the canonical statistical quantity $$a_{i,j} = \frac{1}{2}\left.\frac{\partial^2 \chi}{\partial x_i x_j}\right|_{(x=x_0,y=y_0,z=z_0)}$$

i.e differentiate with respect to the parameters that I've found the best fit values for. The problem is, these quantities were defined as numeric quantities in the above and I wish to produce a symbolic derivative whence evaluated at the optimal values. This leads to an error.

So, my question is how to circumvent this to allow me to obtain a correlation matrix through computing its $ij$th element through the above equation?

Note, I'm aware of alternate procedures such as NonlinearModelFit that computes a correlation matrix automatically however I want to proceed through this route, at least for now. Moreover, I've used NonlinearModelFit to obtain the same fit parameters but different correlation matrix in another exercise so I will post my questions about that at a later date.

Answer 1

There is a package called NumericalCalculus which has a function ND, which takes a numerical derivative. This is a bit odd (to me), in that it takes one-sided forward derivatives, but I'll use it here. This has two options, Terms and Scale, which control the order and the distance between sampling points (play around with the parameters in ND[f[x], x, 1, Terms -> 3,Scale->1] to figure out how this works).

However, for some reason this doesn't seem to have a nice generalisation to partial derivatives. This is a start of such an implementation.

First we load the package:

Needs["NumericalCalculus`"]

Then we define some functions for second differentiation with respect to the same variable and two different variables:

ndd[f_, {a_, a0_}, {a_, a0_}, {x0_, y0_, z0_}, terms_: 5, scale_: 0.1] :=
  ND[f[x, y, z], {a, 2}, a0, Terms -> terms, Scale -> scale] /. {x -> x0, y -> y0, z -> z0}

ndd[f_, {a_, a0_}, {b_, b0_}, {x0_, y0_, z0_}, terms_: 5, scale_: 0.1] /; (!SameQ[a, b]) := 
  ND[ND[f[x, y, z], a, a0, Terms -> terms, Scale -> scale], b, b0, 
       Terms -> terms, Scale -> scale] /. {x -> x0, y -> y0, z -> z0}

These are used in the following function that then takes the second derivative with respect to each variable at a point {x0,y0,z0}.

calculatePartialDerivatives[ff_, {x0_, y0_, z0_}, terms_: 5, scale_: 0.1] := 
 Table[ndd[ff, d1, d2, {x0, y0, z0}, terms, scale], 
    {d1, {{x, x0}, {y, y0}, {z, z0}}}, {d2, {{x, x0}, {y, y0}, {z, z0}}}]

We can test this on a function defined only on numerical input,

testfn[x_?NumericQ, y_?NumericQ, z_?NumericQ] := Sin[x y] + Cos[x z]
calculatePartialDerivatives[testfn, {1, 0.5, 0.2}] // Chop

{{-0.159059, 0.63787, -0.394683}, {0.63787, -0.479426, 0}, {-0.394683, 0, -0.980067}}

which is the same as the explicit result (to 10dp with these settings of terms and scale):

Table[D[Sin[x y] + Cos[x z], d1, d2], {d1, {x,y,z}}, {d2, {x,y,z}}]/. {x->1, y->0.5, z->0.2}

{{-0.159059, 0.63787, -0.394683}, {0.63787, -0.479426, 0}, {-0.394683, 0, -0.980067}}

It is important to either cache the function evaluations, or Simplify the output from calculatePartialDerivatives before plugging the function in. Doing so massively reduces the number of evaluations required:

Clear[testfn, testfn2, testfn3]
nEvals = 0; nEvals2 = 0; nEvals3 = 0; nTerms = 6;
testfn[x_?NumericQ, y_?NumericQ, z_?NumericQ] := (nEvals++; Sin[x y] + Cos[x z])
testfn2[x_?NumericQ, y_?NumericQ, z_?NumericQ] := (nEvals2++; Sin[x y] + Cos[x z])
testfn3[x_?NumericQ, y_?NumericQ, z_?NumericQ] := testfn3[x, y, z] = (nEvals3++; Sin[x y] + Cos[x z])
calculatePartialDerivatives[testfn, {1, 0.5, 0.2}, nTerms, 0.1] // Chop;
Simplify[calculatePartialDerivatives[q, {1, 0.5, 0.2}, nTerms, 0.1]] /. q -> testfn2 // Chop;
calculatePartialDerivatives[testfn3, {1, 0.5, 0.2}, nTerms, 0.1];

These give the same answer (precision dependent), but evaluate their numerical function very different numbers of times:

{nEvals, nEvals2, nEvals3}

(*   {24864, 318, 143}   *)

Adding the Simplify cancels out a lot of the repeated terms in the expansions, while caching the result (testfn3[x_?NumericQ, y_?NumericQ, z_?NumericQ] := testfn3[x, y, z] = ...) gets the same the same effect and is even better as some terms are repeated between the different partial derivatives. It is also marginally faster (as the Simplify takes some time to try and simplify the expression).