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 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矩阵的最小元素可用于估计任何空间方向上的最低斜率。

Answer 2

您可以看到Hessian矩阵为渐变渐变，其中您为此处计算的第一个渐变的每个分量第二次应用渐变是维基百科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的单位矩阵

Answer 3

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

• 斜率aka梯度是方向向量，每个点x y的方向和长度。 这可以通过笛卡尔坐标中的2个数字dx dy或极坐标中的角度θ和长度sqrt( dx^2 + dy^2 )给出。在一系列的山丘上，我们得到了一个vector field。
• Hessians描述了x y附近的曲率，例如：一个抛物面或一个马鞍，有4个数字：dxx dxy dyx dyy。
• 拉普拉斯算子是1个数字，dxx + dyy，在每个点x y。在一系列的山丘上，我们得到了一个scalar field。 （Laplacian = 0的功能或山丘特别光滑。）

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

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个数的方向向量，Hessians是N ^ 2个数的矩阵，拉普拉斯数是1个数。

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