How is the Hessian computed using ``Experimental`NumericalFunction``?

Question

Here and here it was explained how to use Experimental`NumericalFunction to compute the Hessian of a numerical function. I would like to know how this undocumented function works; in particular if it is as robust as Numdifftools, in which [...] Finite differences are used in an adaptive manner, coupled with a Richardson extrapolation methodology to provide a maximally accurate result. [...]

One could consider the following test code:

d = 3;
vecF = {1, 2, 3};
g[vec_?(VectorQ[#, NumericQ] &)] := vec.vec;
h[vec_] := vec.vec;

where g is a numerical function while h is suitable for symbolic evaluation. The Hessian is then obtained in the two cases using the following code:

Block[{e = 0},
  {f = Experimental`CreateNumericalFunction[
  Table[Subscript[x, i], {i, d}], g[Table[Subscript[x, i], {i, d}]], {}, 
  Hessian -> {Automatic, "DifferenceOrder" -> 2},EvaluationMonitor :> e++];
  f["Hessian"[vecF]], e}
]
(*{{{2., 0., 0.}, {0., 2.00001, 0.}, {0., 0., 2.00001}}, 24}*)

and

Block[{e = 0},
  {f = Experimental`CreateNumericalFunction[
  Table[Subscript[x, i], {i, d}], h[Table[Subscript[x, i], {i, d}]], {}, 
  Hessian -> {Automatic, "DifferenceOrder" -> 2},EvaluationMonitor :> e++];
  f["Hessian"[vecF]], e}
]
(*{{{2., 0., 0.}, {0., 2., 0.}, {0., 0., 2.}}, 0}*)

When evaluating the Hessian of g Mathematica uses FiniteDifference as method (with 24 evaluations of the function), while when evaluating the Hessian of h it uses Symbolic and DifferenceOrder is ignored. So, is it possible to know what is happening under the hood?

Answer 1

Converting former comments to an (incomplete) answer.

Well, in the first case, you supress symbolic evaluation by using the pattern _?(VectorQ[#, NumericQ] &). So Mathematica has to switch to numeric differentiation. In the second case, Mathematica is clever enough to differentiate with respect to the symbolic variables that you supplied, e.g. this way: D[h[Table[x[i], {i, d}]], {Table[x[i], {i, d}], 2}]. I am not completely sure, but I think that Experimental`CreateNumericalFunction compiles the Jacobian and the Hessian; since it probably uses reals (doubles) to do so, the result of f["Hessian"[vecF]] is automatically cast to reals.

The code in EvaluationMonitor gets evaluated only if the 0th derivative (e.g., f[vecF]) is evaluated. But when f["Hessian"[vecF]] is called, then f[vecF] will only be called if numerical differentiation is activated (which is only the case for your first example).