n维数组的numpy二阶导数 [英] numpy second derivative of a ndimensional array

问题描述

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

I have a set of simulation data where I would like to find the lowest slope in n dimensions. The spacing of the data is constant along each dimension, but not all the same (I could change that for the sake of simplicity).

我可以忍受一些数字不准确，尤其是边缘.我非常不想生成样条并使用该导数；仅在原始值上就足够了.

I can live with some numerical inaccuracy, especially towards the edges. I would heavily prefer not to generate a spline and use that derivative; just on the raw values would be sufficient.

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

It is possible to calculate the first derivative with numpy using the numpy.gradient() function.

import numpy as np

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

<小时>

这是关于拉普拉斯矩阵与黑森矩阵的注释；这不再是一个问题，而是为了帮助理解未来的读者.

---

This is a comment regarding laplace versus the hessian matrix; this is no more a question but is meant to help understanding of future readers.

我使用 2D 函数作为测试用例来确定阈值以下的最平坦"区域.下图显示了使用second_derivative_abs = np.abs(laplace(data))的最小值和以下的最小值的结果差异:

I use as a testcase a 2D function to determine the 'flattest' area below a threshold. The following pictures show the difference in results between using the minimum of second_derivative_abs = np.abs(laplace(data)) and the minimum of the following:

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

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

The color scale depicts the functions values, the arrows depict the first derivative (gradient), the red dot the point closest to zero and the red line the threshold.

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

The generator function for the data was ( 1-np.exp(-10*xi**2 - yi**2) )/100.0 with xi, yi being generated with np.meshgrid.

拉普拉斯:

黑森州:

推荐答案

二阶导数由 Hessian 给出矩阵.这是 ND 数组的 Python 实现，包括应用 np.gradient 两次并适当地存储输出，

The second derivatives are given by the Hessian matrix. Here is a Python implementation for ND arrays, that consists in applying the np.gradient twice and storing the output appropriately,

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 矩阵的迹(对角线元素的和).

Note that if you are only interested in the magnitude of the second derivatives, you could use the Laplace operator implemented by scipy.ndimage.filters.laplace, which is the trace (sum of diagonal elements) of the Hessian matrix.

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

Taking the smallest element of the the Hessian matrix could be used to estimate the lowest slope in any spatial direction.

这篇关于n维数组的numpy二阶导数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！