Central Finite Differences for computing Gradient and Hessian

cfd_grad_hess.py

#!/usr/bin/env python

import numpy as np

def gradient(x, func, eps=1e-6):
    """
    Central finite differences for computing gradient

    INPUT:
       x: vector of parameters with shape (n, 1)
       func: function to compute the gradient on. Take x as a parameter.
       eps: parameter for the central differences

    OUTPUT:
       grad: vector with the values of the gradient with shape (n, 1)
    """
    # Size the problem, i.e. nbr of parameters
    n = len(x)

    # Prepare the vector for the gradient
    grad = np.zeros(n)

    # Prepare the array to add epsilon to.
    dx = np.zeros(n)

    # Go through all parameters
    for i in range(len(x)):
        # Add epsilon to variate a parameter
        dx[i] += eps

        # Central finite differences
        grad[i] = -(func(x+dx) - func(x-dx))/(2*eps)

        # Set back to 0
        dx[i] = 0

    return grad

def hessian(x, func, eps=1e-6):
    """
    Central finite differences for computing the hessian

    INPUT:
       x: vector of parameters with shape (n, 1)
       func: function to compute the gradient on. Take x as a parameter.
       eps: parameter for the central differences

    OUTPUT:
       grad: vector with the values of the gradient with shape (n, 1)
    """

    # Size the problem, i.e. nbr of parameters
    n = len(x)

    # Prepare the vector for the gradient
    hess = np.zeros((n,n))

    # Prepare the array to add epsilon to.
    dx = np.zeros(n)

    # Go through all parameters
    for i in range(n):
        # Add epsilon to variate a parameter
        dx[i] += eps

        # Compute the gradient with forward and backward difference
        grad_plus = gradient(x+dx, func, eps)
        grad_minus = gradient(x-dx, func, eps)

        # Central finite difference
        hess[i,:] = -(grad_plus - grad_minus)/(2*eps)

        # Set back to 0
        dx[i] = 0

    return hess

test.py

import time
import numpy as np
import numdifftools as nd
from cfd_grad_hess import gradient, hessian

# Create a function to test our implementation
func = lambda x: 3*x[0]**2 + 2*x[1]**3 + x[0]*x[1] + x[2]

# Parameter
x = [1,1,1]

print('Function value at [{}]: {}'.format(', '.join(map(str, x)), func(x)))

# Compute the gradient
start = time.time()
grad = gradient(x, func, eps=1e-4)
stop = time.time()

print('Gradient is equal to ({}) and has been computed in {:.2f}µs.'.format(', '.join(map(str, grad)), 1e6*(stop-start)))

# Compute the hessian
start = time.time()
hess = hessian(x, func, eps=1e-4)
stop = time.time()

print('Hessian is equal to ({}) and has been computed in {:.2f}µs.'.format(', '.join(map(str, hess)), 1e6*(stop-start)))

# Compare our implementation with numdifftools
start = time.time()
nd_hess = nd.Hessian(func)
nd_hess_vals = nd_hess(x)
stop = time.time()

print('Hessian from numdifftools is equal to ({}) and has been computed in {:.2f}µs.'.format(', '.join(map(str, nd_hess_vals)), 1e6*(stop-start)))

print('Norm of the difference between the two hessians: {:.2E}'.format(np.linalg.norm(hess - nd_hess_vals)))

After struggling to find a good finite difference implementation of a Hessian (and numdifftools being way too slow), I decided to share my own implementation.

You can use the file test.py in order to test the implementation and compare it to the Hessian computed from numdifftools. It's possible that you may have to change the parameter eps to get a better result. However, the smaller does not always mean the better result.