二阶导数由Hessian matrix给出。下面是ND数组的Python实现，包括应用np.gradient两次并适当存储输出，

import numpy as np

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):
        # iterate over dimensions
        # apply gradient again to every component of the first derivative.
        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)

字符串
注意，如果你只对二阶导数的大小感兴趣，你可以使用scipy.ndimage.filters.laplace实现的Laplace operator，它是Hessian矩阵的迹（对角元素之和）。
取Hessian矩阵的最小元素可用于估计任何空间方向上的最低斜率。

斜坡、黑森曲线和拉普拉斯曲线是相关的，但它们是三个不同的东西。
从2d开始：2个变量的函数（x，y），例如，一系列山丘的高度图，

• 斜率又称为梯度是方向向量，即在每个点x y处的方向和长度。

在笛卡尔坐标系中，这可以用2个数dx dy表示，在极坐标系中，这可以用Angular θ和长度sqrt( dx^2 + dy^2 )表示，在整个山丘范围内，我们得到vector field。

• 黑森描述了x y附近的曲率，例如抛物面或鞍形，有4个数字：dxx dxy dyx dyy。
• 拉普拉斯算子是一个数，dxx + dyy，在每个点x y。在一个山丘范围内，我们得到一个scalar field。（带有Laplacian = 0的函数或山丘是特别平滑的。）

斜率为线性拟合和Hessians二次拟合，用于点xy附近的微小步长h：

f(xy + h)  ~  f(xy)
        +  slope . h    -- dot product, linear in both slope and h
        +  h' H h / 2   -- quadratic in h

字符串
这里，xy、slope和h是2个数的向量，并且H是4个数的矩阵dxx dxy dyx dyy。
N-d是类似的：斜率是N个数的方向向量，海森是N^2个数的矩阵，拉普拉斯是1个数，在每个点。
(You可能会在math.stackexchange上找到更好的答案。）

你可以看到Hessian矩阵是一个梯度的梯度，你可以对第一个梯度的每个分量应用第二次梯度计算这里是一个wikipedia link definig Hessian矩阵，你可以清楚地看到这是一个梯度的梯度，这里是一个Python实现定义梯度然后hessian：

import numpy as np
#Gradient Function
def gradient_f(x, f):
  assert (x.shape[0] >= x.shape[1]), "the vector should be a column vector"
  x = x.astype(float)
  N = x.shape[0]
  gradient = []
  for i in range(N):
    eps = abs(x[i]) *  np.finfo(np.float32).eps 
    xx0 = 1. * x[i]
    f0 = f(x)
    x[i] = x[i] + eps
    f1 = f(x)
    gradient.append(np.asscalar(np.array([f1 - f0]))/eps)
    x[i] = xx0
  return np.array(gradient).reshape(x.shape)

#Hessian Matrix
def hessian (x, the_func):
  N = x.shape[0]
  hessian = np.zeros((N,N)) 
  gd_0 = gradient_f( x, the_func)
  eps = np.linalg.norm(gd_0) * np.finfo(np.float32).eps 
  for i in range(N):
    xx0 = 1.*x[i]
    x[i] = xx0 + eps
    gd_1 =  gradient_f(x, the_func)
    hessian[:,i] = ((gd_1 - gd_0)/eps).reshape(x.shape[0])
    x[i] =xx0
  return hessian

字符串
作为测试，（x^2 + y^2）的Hessian矩阵是2 * I_2，其中I_2是2维单位矩阵

hessians = np.asarray(np.gradient(np.gradient(f(X, Y))))
hessians[1:]

字符串
对三维函数f有用。

你有没有试过sympy hessian（二阶导数矩阵）或jacobian（一阶导数矩阵）？试试这个，对于你的函数：

import sympy as sp
sp.init_printing()

x_i,y_i = sp.symbols("x_i,y_i")
f = sp.Function("f")(x_i,y_i)
f = (1 - sp.exp(-10*x_i**2-y_i**2))/100
# Jacobian matrix - first derivative
J = sp.Matrix([f]).jacobian([x_i,y_i])
J
# Hessian matrix - second derivative
H = sp.hessian(f, (x_i, y_i))
H

字符串
它产生一个Hessian：