# 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?

## 1 Answer

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

• thanks but what i would like to know is how the numerical differentiation works, especially regarding the choice of the steps for the finite differences evaluation Commented May 30, 2018 at 12:00
• That's good question. Since it is about internals, I cannot say more. But maybe somebody else can provide more detail... Commented May 30, 2018 at 12:28