# Caclulate second partial derivatives of a numerically defined function

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.

• Your final code section doesn't work, it has an extra ]]. And rr and M are undefined. Can you ask your question as a minimal working example that gets to the question, which is how to take a numerical derivative. – SPPearce Jan 2 '18 at 16:52
• To be precise, your issue is that you have a function $f(x,y,z)$ that you are trying to find the second order derivatives of. – SPPearce Jan 2 '18 at 17:00
• @KraZug Yes, sorry as noted in the question, M can be set to one for convenience and rr[] represents the best fit parameters. But since I am just interested in how to compute the derivative of the chi square involving numerically defined variables, it doesn't matter what these best fit values are so can also just be set to arbitrary values. – CAF Jan 2 '18 at 22:19
• @KraZug Thanks, yes I have a function (namely the chi square) which I wish to take a derivative of wrt numerically defined variables a,b,NN . Is it clear? Or should I rephrase? – CAF Jan 2 '18 at 22:22
• I've heavily edited the question to make it general and helpful to other people. And I've also answered it. You should be able to figure out how to use my answer for your specific function. – SPPearce Jan 3 '18 at 10:13

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.

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).