一个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 (:

---

这是一条关于拉普拉斯和黑塞矩阵的评论；这不再是一个问题，而是旨在帮助未来读者理解。

我使用一个二维函数作为测试用例来确定低于阈值的“最平坦”区域。以下图片展示了使用 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生成。

拉普拉斯：

海森矩阵：

- Faultier

---

我想知道的是斜率的大小，而不是方向。如果我计算绝对一阶导数列表条目之和的梯度，是否足够？second_derivative = np.gradient(sum([df*df for d in first_derivative]))（为了论证而保留形状的sum）。 - Faultier

好的，我现在明白你想要什么了。你只是想得到最平坦的区域，或者说在N维空间中"最平坦"的意思是什么。我建议不要使用二阶导数，而是计算所有点的绝对梯度（对np.gradient结果的第一维的平方求和，就像你在评论中提到的那样），然后从中找到阈值区域，并在阈值区域内找到最小值（如果函数足够复杂，全局最小值可能很难找到）。我会尝试一下，如果我找到了什么东西，就会发布另一个答案。 - Carsten

@Carsten 所有梯度的总和不会产生最平坦的区域；在这个测试案例中，它将产生2D高斯的中心，这绝不是最平坦的区域。因此，我认为需要使用适当的二阶导数来完成，而不是一阶导数。 - Faultier

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)

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

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

- rth

---

我目前正在对2D测试案例进行两种提议解决方案的实验（我喜欢看东西）。当我说到斜率时，我的意思是在所有方向上最小的斜率；基本上是我能找到的局部最平坦的区域。目前看起来很有前途，但还需要进行一些测试 (: - Faultier

我没有足够的空间来发布2D测试用例，以免污染问题；拉普拉斯和黑塞矩阵之间的结果差异似乎在于它们产生不同的点。对于黑塞矩阵，我评估x.dim、x.ndim沿着平方和的最小值或总和。据我了解，后者考虑了混合导数，因此对于我的目的应该更准确？ - Faultier

我添加了Laplace和Hessian的2D测试用例的图片。虽然我认为这两种算法都做得很好，但我认为我更喜欢Hessian并将其扩展到可变步长；正如我们讨论的那样，它可以更正确地完成工作。编辑：我刚刚看到您在回答中的编辑；因此，您能否请查看我的图片/代码，并告诉我是否真的做对了(: - Faultier

对于已经实现的海森矩阵计算算法，请查看 Numdifftools 项目: https://pypi.org/project/numdifftools/ - Dr_Zaszuś

Answer 2

斜率、黑塞矩阵和拉普拉斯算子是相关的，但它们是三个不同的概念。
从二维开始：一个关于两个变量x和y的函数，例如一片山区的高度图。

• 斜率（又称梯度）是方向向量，表示每个点x y的方向和长度。
  这可以通过笛卡儿坐标系中的2个数字dx dy， 或者是极坐标系中的角度θ和长度sqrt(dx^2+dy^2)来表示。 在整个山区范围内，我们得到一个矢量场。

• 海森矩阵描述了接近x y处的曲率，例如抛物面或鞍点，具有4个数字：dxx dxy dyx dyy。

• 拉普拉斯算子是每个点x y的1个数字，即dxx+dyy。 在一定范围内的山丘上，我们得到一个标量场。 （具有Laplacian = 0的函数或山丘特别平滑。）

斜率是线性拟合，海森矩阵是二次拟合，对于在点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个数字组成的方向向量， Hessian是由N^2个数字组成的矩阵，每个点的Laplacian为1个数字。

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

- denis

Answer 3

您可以将Hessian矩阵视为梯度的梯度，其中对于此处计算的第一个梯度的每个分量，您会再次应用梯度。这是维基百科link定义Hessian矩阵的清晰说明，您可以清楚地看到它是梯度的梯度。以下是定义梯度然后Hessian的Python实现：

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) 的海森矩阵是 2 * I_2，其中 I_2 是 2 维单位矩阵。

- SENHAJI RHAZI Hamza

---

我认为Hessian矩阵是梯度的雅可比矩阵，而不是梯度的梯度。在这里查看维基百科链接：https://en.wikipedia.org/wiki/Hessian_matrix - Srikiran

@Srikiran 这有点语义上的问题。雅可比矩阵实际上只是为向量值函数定义的梯度。 - Tyberius

Answer 4

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

曾经为三维函数 f 工作。

- MaizeMaze

Answer 5

你试过用sympy的海森矩阵（二阶导数矩阵）或雅可比矩阵（一阶导数）了吗？试试这个，针对你的函数：

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

它生成了一个赫塞矩阵：

- Carlos Camargo