def hessian(x):
    """
    Calculate the hessian matrix with finite differences
    Parameters:
       - x : ndarray
    Returns:
       an array of shape (x.dim, x.ndim) + x.shape
       where the array[i, j, ...] corresponds to the second derivative x_ij
    """
x_grad = np.gradient(x)
hessian = np.empty((x.ndim, x.ndim) + x.shape, dtype=x.dtype)
for k, grad_k in enumerate(x_grad):
    tmp_grad = np.gradient(grad_k)
for l, grad_kl in enumerate(tmp_grad):
    hessian[(k), (l), :, :] = grad_kl
return hessian
x = np.random.randn(100, 100, 100)
hessian(x)
numpy second derivative of a ndimensional array