# Numeric calculation of Hessian

I want to calculate the Hessian matrix for a function that can only be evaluated numerically. So far, I have the following (where f is just for testing):

Needs["NumericalCalculus"]

np = 3;  (* number of parameters *)
f0[{x_, y_, z_}] := x^2*y*z^3 + z + 1
(* numeric version; c to count evaluations *)
f[{pars__?NumericQ}] := (c += 1; f0[{pars}])

(* construct symmetric hessian *)
hessian[{pars__?NumericQ}] := Module[{diag, urt, temp1, temp2},
(* diagonal *)
diag = Table[
ND[f[ReplacePart[{pars}, i -> temp1]], {temp1,
2}, {pars}[[i]]], {i, 1, np}
];
(* upper right triangle *)
urt = Table[
If[j > i,
ND[ND[f[ReplacePart[{pars}, {i -> temp1, j -> temp2}]],
temp1, {pars}[[i]]], temp2, {pars}[[j]]],
0
],
{i, 1, np}, {j, 1, np}
];
(* result *)
Table[
If[i == j, diag[[i]],
If[j > 1,
urt[[i, j]],
urt[[j, i]]
]
],
{i, 1, np}, {j, 1, np}
]
];


This seems to work:

 c = 0;
hessian[{3, 4, 5}] // MatrixForm


yields the expected result. But the number of function evaluations appears to be rather large (I got c = 3099). Is this normal for ND, or can the above calculation be improved?

• The culprit is ND[ND[f[..]..]..] -- the inside ND especially, for which the j-th parameter temp2 is not numeric. I don't have an solution, yet, though. – Michael E2 Jul 21 '13 at 16:05
• Closely related: Numerical partial derivative – Jens Jul 21 '13 at 23:28

As I mentioned in a comment the culprit is the nested ND in ND[ND[f[..],..],..]. The inside ND is of the form below. The function g below (f in your code) does not evaluate to a numerical result because at least one of the arguments remains a symbol. ND returns an expression that can be used, but it has 128 function calls to g in it.

g[x0_?NumericQ, y0_?NumericQ, z0_?NumericQ] := (foo++; x0 + y0 + z0);
Count[ND[g[x, y, 2.], x, 1.], _g, Infinity]
(* 128 *)


With a normal numeric function, ND does eight evaluations:

foo = 0;
ND[g[x, 1., 2.], x, 1.]
foo
(* 1. *)
(* 8  *)


So for a mixed partial derivative, one might hope for 64 evaluations from nested ND. One might hope for better since second order, non-mixed partials require 9 function evaluations. The numerical calculus package does not provide such methods, and I don't know if they exist. Mixed partials are higher dimensional and estimating them may be intrinsically more difficult.

The key is to prevent ND from being evaluated until the function has all numeric arguments. To do that we set up partial derivative functions that hold up evaluation of ND until they have all numeric arguments. This can be done using UpSetDelayed to associate such a partial derivative of f with an "upvalue" of f.

c = 0;
Needs["NumericalCalculus"];
Clear[f, fx, fy, fz, f0, x, y, z];
f0[{x_, y_, z_}] := x^2*y*z^3 + z + 1
f[{x0_?NumericQ, y0_?NumericQ, z0_?NumericQ}] := (c++; f0[{x0, y0, z0}]);

fx[{x0_?NumericQ, y0_?NumericQ, z0_?NumericQ}] := ND[f[{x, y0, z0}], x, x0];
ND[f[{x_Symbol, y_Symbol, z_} | {x_Symbol, y_, z_Symbol}], x_, x0_?NumericQ] ^:=
fx[{x0, y, z}];

fy[{x0_?NumericQ, y0_?NumericQ, z0_?NumericQ}] := ND[f[{x0, y, z0}], y, x0];
ND[f[{x_Symbol, y_Symbol, z_} | {x_, y_Symbol, z_Symbol}], y_, y0_?NumericQ] ^:=
fy[{x, y0, z}];

fz[{x0_?NumericQ, y0_?NumericQ, z0_?NumericQ}] :=  ND[f[{x0, y0, z}], z, z0];
ND[f[{x_Symbol, y_, z_Symbol} | {x_, y_Symbol, z_Symbol}], z_, z0_?NumericQ] ^:=
fz[{x, y, z0}];


Then the number of function calls is reduce to 219 = 3 x 64 + 3 x 9 (according to the two kinds of second-order partial derivatives, mixed and non-mixed):

c = 0;
hessian[{3, 4, 5}] // MatrixForm
c
(* 1000.  750.  1800.
750.     0.   675.
1800.    0   1080.  *)
(* 219 *)

• Thank you for the improvement and explanation. I wasn't aware of the UpSetDelayed command. This should be quite useful. – muton Jul 21 '13 at 20:16
• Did you forget to define the hessian in your answer? Or is it just me? – chris Jul 21 '13 at 20:54
• @chris Use the OP's hessian. I left it out because it was long. (Hope that's ok.) – Michael E2 Jul 21 '13 at 21:58
• @Michael E2 would it be possible to generalize your answer to np dimensions? Also, how can I use the OP's hessian if fx etc are not defined? – Valerio Sep 1 '17 at 19:42

For what it's worth, here's what I wrote for myself. (It doesn't use ND.)

NHessian::usage = "NHessian[f, x] computes a numerical approximation \
to the Hessian matrix evaluated at f[x]. NHessian take the option \
Scale, which can be a scalar or a vector (matching the length of the \
vector x). The default value is Scale -> \!$$\*SuperscriptBox[\(10$$, \
$$-3$$]\)."

Options[NHessian] = {Scale -> 10^-3}

NHessian[f_, x_?(VectorQ[#, NumericQ] &), opts___?OptionQ] :=
Module[{n, h, norm, z, mat, f0},
n = Length[x];
h = Scale /. {opts} /. Options[NHessian];
norm = If[VectorQ[h], Outer[Times, 2 h, 2 h], 4 h^2];
z = If[VectorQ[h], DiagonalMatrix[h], h*IdentityMatrix[n]];
mat = ConstantArray[0., {n, n}];
f0 = f[x];
Do[
mat[[i, j]] =
If[i == j,
(* then *)
.5 (f[x + 2 * z[[i]]] - 2 f0 + f[x - 2 * z[[i]]]),
(* else *)
f[x + z[[i]] + z[[j]]] - f[x + z[[i]] - z[[j]]] -
f[x - z[[i]] + z[[j]]] + f[x - z[[i]] - z[[j]]]
],
{i, n}, {j, i, n}
];
(mat + Transpose[mat])/norm
]


Let's apply this to the given function:

f0[{x_, y_, z_}] := (c++; x^2*y*z^3 + z + 1)

c = 0;
NHessian[f0, {3, 4, 5}]\\MatrixForm
c
(* 1000.  750.  1800.
750.    0.   675.
1800.  675.  1080.  *)
(* 19 *)

• Really nice. This is an even faster solution. Many thanks. – muton Jul 24 '13 at 21:40
• I liked the simplicity of this solution! Would it make sense to use as default for Scale something proportional to the vector vec at which the Hessian is being evaluated? Eg, Scale->Abs[vec] 10^-5. – Valerio Apr 6 '17 at 20:27