python - n维数组的numpy二阶导数

Question

我有一组模拟数据，我想在其中找到 n 维中的最低斜率。数据的间距在每个维度上都是恒定的，但并不完全相同(为了简单起见，我可以更改它)。

我可以忍受一些数值上的不准确，尤其是在边缘处。我非常不希望生成样条曲线并使用该导数；仅基于原始值就足够了。

可以使用 numpy.gradient() 函数计算 numpy 的一阶导数。

import numpy as np

data = np.random.rand(30,50,40,20)
first_derivative = np.gradient(data)
# second_derivative = ??? <--- there be kudos (:

---

这是关于拉普拉斯与海森矩阵的评论；这不再是一个问题，而是为了帮助 future 的读者理解。

我使用二维函数作为测试用例来确定低于阈值的“最平坦”区域。以下图片显示了使用 second_derivative_abs = np.abs(laplace(data)) 的最小值与以下最小值的结果差异:

second_derivative_abs = np.zeros(data.shape)
hess = hessian(data)
# based on the function description; would [-1] be more appropriate? 
for i in hess[0]: # calculate a norm
    for j in i[0]:
        second_derivative_abs += j*j

色标表示函数值，箭头表示一阶导数(梯度)，红点表示最接近零的点，红线表示阈值。

数据的生成器函数是 ( 1-np.exp(-10*xi**2 - yi**2) )/100.0 其中 xi、yi 是用 生成的np.meshgrid.

拉普拉斯:

黑森州:

Answer 1

二阶导数由 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)

请注意，如果您只对二阶导数的大小感兴趣，则可以使用 Laplace operator由 scipy.ndimage.filters.laplace 实现, 是Hessian矩阵的迹(对角线元素之和)。

取 Hessian 矩阵的最小元素可用于估计任何空间方向上的最低斜率。