Numeric calculation of Hessian

Question

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?

Thank you for any answers/comments.

Answer 1

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

Answer 2

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