Numeric mixed derivatives

Question

In short

I need NumericCalculus`ND extended to mixed partial derivatives. This can be done by nesting them, but care has to be taken to evaluate underling ND only when the parameters become numeric. There are examples in the answers on "Numerical partial derivative" and "Numeric calculation of Hessian" (make sure to notice the discussion on the efficiency). The challenge is to build such patterns with arbitrary number of parameters.

Long story

I am forced to take the derivative of my 200k leafcount expressions numerically. I am aware of NumericCalculus`ND which works great, but does not allow mixed derivatives. There are some examples on how to do it efficiently in this and this answer. I would like to play around with a greater number of parameters and use a Scale option for each of them. So I would like to have a nice top level function.

An example of what I would like it to be capable of:

nd[Sin[x Conjugate[y]] Cos[u v] /.{u->u0, v->v0}, {x, x0}, {y, y0}, Scale -> {xScale, yScale}]

To define such a function with variable number of parameters and understanding the syntax in some of the examples is beyond my capabilities. Could someone provide the function? It is a pity that Mathematica does not ship with one.

Update to the long story: Both of the related answers I linked discuss how the underling ND should be called only, when all parameters are finally numeric. This avoids big number of ND calls with partially numeric expressions. I would like the method to have the same feature and not just nesting the ND calls.

Update 2 As I did not consider it relevant, I did not mention that I have complex variables. In this case finite difference method is ill defined and some simple reimplementation of ND would not do.

Answer 1

Initial answer

I have similar situation (but I numerically compute only first derivatives) and found that the most reliable and efficient solution is do not rely on ND but implement numerical derivatives by yourself using the definition

$$f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}$$

Such approach requires at the first step finding optimal $h$ for each of the parameters which gives maximum precision. It cannot be too small because in this case the precision of the difference $f(a+h)-f(a)$ will vanish. In my case defining $h$ as

$$h=\frac{\log_2(|a|)}{2^n}$$

reduces the problem to finding optimal value of $n$ which is equal for all the parameters. For this purpose I compute numerical approximations for all partial derivatives with $n$ running from 15 to 60 (I compute with high WorkingPrecision) and see where the precision of these approximations becomes sufficient or maximal for all the partial derivatives.

UPDATE: efficient use of ND

Thanks to acl's great answer now it is clear that ND with Terms->1 uses exactly the same finite difference formula:

Needs["NumericalCalculus`"]
ND[f[u], u, x, Terms -> 1, Scale -> h, WorkingPrecision -> Infinity]

(-f[x] + f[h + x])/h

But one should be careful here: increasing the number of terms leads to division of the step $h$ by powers of 2 in successive higher order terms:

ND[f[u], u, x, Terms -> 2, Scale -> h, WorkingPrecision -> Infinity] // Simplify

-((3 f[x] - 4 f[h/2 + x] + f[h + x])/h)

ND[f[u], u, x, Terms -> 3, Scale -> h, WorkingPrecision -> Infinity] // Simplify

(-21 f[x] + 32 f[h/4 + x] - 12 f[h/2 + x] + f[h + x])/(3 h)

It means that for getting the best possible approximation with ND one should estimate the step $h$ by the method described in the initial section of this answer, then multiply it by $2^{t-1}$ (where $t$ is the value of the Terms option of ND) and use as value for the Scale option.

These considerations can be compiled into a function with more intuitive and consistent behavior:

Options[myND] = {Terms -> 1, WorkingPrecision -> MachinePrecision};
myND[expr_, x_, x0_, h_, opts : OptionsPattern[]] := 
 ND[expr, x, x0, Scale -> h*2^(OptionValue[Terms]-1), opts]

Answer 2

If I understand correctly, you're trying to figure out how to generate expressions like

ND[ND[ND[g[x, y, z], x, 1], y, 1], z, 1]

(taken from the first link you provided). If so, I think what you want is something along these lines:

ClearAll[nd];
SetAttributes[nd, HoldFirst];
Options[nd] = {Scale -> Automatic};
nd[func_, params__List, OptionsPattern[]] :=
  With[
   {vars = {params}[[All, 1]], inits = {params}[[All, 2]]},
   Fold[
    NDFunc[#1, Sequence @@ MapAt[Scale -> # &, #2, -1]] &,
    Apply[func, vars],
    Transpose@{
      vars,
      inits,
      OptionValue[Scale] /. 
       Automatic :> ConstantArray[1, Length@vars]}]];

So:

nd[f, {x, x0}, {y, y0}, {z, z0}]

NDFunc[NDFunc[NDFunc[f[x, y, z], x, x0, Scale -> 1], y, y0, Scale -> 1], z, z0, Scale -> 1]

nd[f, {x, x0}, {y, y0}, {z, z0}, Scale -> {xscale, yscale, zscale}]

NDFunc[NDFunc[NDFunc[f[x, y, z], x, x0, Scale -> xscale], y, y0, Scale -> yscale], z, z0, Scale -> zscale]