# How is the Hessian computed using ExperimentalNumericalFunction?

Here and here it was explained how to use ExperimentalNumericalFunction 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 = ExperimentalCreateNumericalFunction[
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 = ExperimentalCreateNumericalFunction[
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?

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