Numerical Hessian

Estimate the Hessian matrix of a scalar function with central finite differences.

The Hessian H[i, j] = d^2 f / (dx_i dx_j) describes the curvature of f. Newton's method uses it, and its eigenvalues tell you whether a point where the gradient is zero is a minimum, a maximum or a saddle.

Perturb two coordinates at a time by h. With e_i the unit vector along coordinate i:

H[i, j] = ( f(x + h e_i + h e_j) - f(x + h e_i - h e_j)
          - f(x - h e_i + h e_j) + f(x - h e_i - h e_j) ) / (4 h^2)

The same formula works for i == j, where it becomes a second difference with step 2h.

Function signature

def numerical_hessian(f, x, h=1e-4):

Example

Input

f = lambda v: v[0] ** 2 * v[1]
numerical_hessian(f, np.array([1.0, 2.0]))

Output

array([[4., 2.],
       [2., 0.]])

Reasoning

By hand, d^2f/dx0^2 = 2 * x1 = 4, d^2f/dx0 dx1 = 2 * x0 = 2, and d^2f/dx1^2 = 0. The finite differences match up to a rounding error of about 1e-7.

Constraints and notes

• x is a 1D float array. f returns a float.
• The result must be symmetric.
• Do not modify x.