n维数组的numpy二阶导数

Question

我有一组模拟数据,我想在n维中找到最低的斜率.数据的间距沿着每个维度是恒定的,但不是全部相同(为了简单起见,我可以改变它).

我可以忍受一些数字不准确,尤其是边缘.我非常希望不生成样条并使用该衍生物; 只要原始价值就足够了.

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

import numpy as np

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

---

这是关于拉普拉斯与粗麻布矩阵的评论; 这不再是一个问题,而是为了帮助理解未来的读者.

我使用2D函数作为测试用例来确定阈值以下的"最平坦"区域.以下图片显示了使用以下最小值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矩阵给出.这是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,即Hessian矩阵的轨迹(对角元素之和).

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

Answer 2

斜率、Hessians 和 Laplacian 是相关的，但是是 3 个不同的东西。
\n从 2d 开始：2 个变量的函数( x, y )，例如一系列山丘的高度图，

\n\n

\n

• 斜率又名梯度是方向向量，每个点的方向和长度x y。
  \n这可以由笛卡尔坐标中的 2 个数字给出dx dy，\也可以由极坐标中的角度 \xce\xb8 和长度给出sqrt( dx^2 + dy^2 )。\n在整个山丘范围内，我们得到一个\n矢量场。

\n

• Hessians 描述 附近的曲率x y，例如抛物面或鞍形，\n用 4 个数字：dxx dxy dyx dyy。

\n

• 拉普拉斯算子dxx + dyy在每个点 处都是 1 个数字x y。\n在一系列山丘上，我们得到\n标量场。\n（拉普拉斯算子 = 0的函数或山丘\n 特别平滑。）

\n

\n\n

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

\n\n

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

\n\n

这里xy、slope和h是 2 个数字的向量，\nandH是 4 个数字的矩阵dxx dxy dyx dyy。

\n\n

Nd 类似：斜率是 N 个数字的方向向量，\nHessians 是 N^2 个数字的矩阵，而拉普拉斯矩阵是每个点上的 1 个数字。

\n\n

（您可能会在math.stackexchange上找到更好的答案。）

\n