numpy中的二阶梯度

Question

我试图在numpy数值计算数组的二阶梯度.

a = np.sin(np.arange(0, 10, .01))
da = np.gradient(a)
dda = np.gradient(da)

这就是我的想法.应该这样做吗？

我问这个,因为在numpy中没有选项说np.gradient(a,order = 2).我担心这种用法是否错误,这就是为什么numpy没有实现这一点.

PS1:我确实知道有np.diff(a,2).但这只是单边估计,所以我很好奇为什么np.gradient没有类似的关键字.

PS2:np.sin()是玩具数据 - 真实数据没有分析形式.

谢谢!

Answer 1

我将第二句@jrennie的第一句话 - 它可以全部依赖.numpy.gradient函数要求数据均匀分布(尽管如果是多维的,每个方向允许不同的距离).如果您的数据不符合这一点,那么numpy.gradient将不会有太大用处.实验数据可能具有(OK,将具有)噪声,除了不必均匀间隔.在这种情况下,最好使用scipy.interpolate样条函数(或对象)之一.这些可以采用不均匀间隔的数据,允许平滑,并且可以返回到k-1的导数,其中k是所请求的样条拟合的顺序.k的默认值是3,所以二阶导数就好了.例:

spl = scipy.interpolate.splrep(x,y,k=3) # no smoothing, 3rd order spline
ddy = scipy.interpolate.splev(x,spl,der=2) # use those knots to get second derivative 

像scipy.interpolate.UnivariateSpline这样的面向对象的样条线具有导数的方法.注意,导数方法在Scipy 0.13中实现,并且不存在于0.12中.

请注意,正如@JosephCottham在2018年的评论中指出的那样,这个答案(至少对Numpy 1.08有利)已不再适用,因为(至少)Numpy 1.14.检查您的版本号和呼叫的可用选项.

Answer 2

对于一阶导数的不连续性，双梯度方法会失败。由于梯度函数考虑向左和向右的一个数据点，因此在多次应用它时，这种情况会继续/传播。

另一方面，二阶导数可以通过以下公式计算

d^2 f(x[i]) / dx^2 = (f(x[i-1]) - 2*f(x[i]) + f(x[i+1])) / h^2

在这里比较一下。这样做的优点是只考虑两个相邻像素。

在图中，比较了由np.diff实现的双np.gradient方法（左）和上述公式（右） 。由于 f(x) 在零处只有一个扭结，因此二阶导数（绿色）应该只在那里有一个峰值。由于双梯度解考虑了每个方向上的 2 个相邻点，这导致二阶导数值有限，为 +/- 1。

然而，在某些情况下，您可能更喜欢双梯度解决方案，因为这对噪声更稳健。

我不确定为什么存在np.gradient和np.diff，但原因可能是， np.gradient 的第二个参数定义了像素距离（对于每个维度），并且对于图像，它可以同时应用于两个维度gy, gx = np.gradient(a)。

代码

import numpy as np
import matplotlib.pyplot as plt

xs = np.arange(-5,6,1)
f = np.abs(xs)
f_x = np.gradient(f)
f_xx_bad = np.gradient(f_x)
f_xx_good = np.diff(f, 2)
test = f[:-2] - 2* f[1:-1] + f[2:]





# lets plot all this

fig, axs = plt.subplots(1, 2, figsize=(9, 3), sharey=True)

ax = axs[0]
ax.set_title('bad: double gradient')
ax.plot(xs, f, marker='o', label='f(x)')
ax.plot(xs, f_x, marker='o', label='d f(x) / dx')
ax.plot(xs, f_xx_bad, marker='o', label='d^2 f(x) / dx^2')
ax.legend()

ax = axs[1]
ax.set_title('good: diff with n=2')
ax.plot(xs, f, marker='o', label='f(x)')
ax.plot(xs, f_x, marker='o', label='d f(x) / dx')
ax.plot(xs[1:-1], f_xx_good, marker='o', label='d^2 f(x) / dx^2')
ax.plot(xs[1:-1], test, marker='o', label='test', markersize=1)
ax.legend()

Answer 3

数值梯度计算没有通用的正确答案。在计算样本数据的梯度之前，您必须对生成该数据的基础函数做出一些假设。您可以在技术上np.diff用于梯度计算。使用np.gradient是一种合理的方法。我看不出你所做的事情有什么根本性的错误——它是一维函数的二阶导数的一种特殊近似。