我正在尝试以数字方式计算 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() 是一个玩具数据 - 真实数据没有解析形式。
谢谢！
【问题讨论】：
我将第二个 @jrennie 的第一句话 - 这一切都取决于。 numpy.gradient 函数要求数据均匀分布（尽管如果是多维的，则允许每个方向上的距离不同）。如果您的数据不遵守这一点，那么 numpy.gradient 就没有多大用处。实验数据可能有（好的，会有）噪音，除了不一定都是均匀分布的。在这种情况下，最好使用 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 年在 cmets 中指出的那样，这个答案（至少对 Numpy 1.08 有好处）自（至少）Numpy 1.14 起不再适用。检查您的版本号和调用的可用选项。
【讨论】：
数值梯度计算没有通用的正确答案。在计算样本数据的梯度之前，您必须对生成该数据的底层函数做出一些假设。从技术上讲，您可以使用np.diff 进行梯度计算。使用np.gradient 是一种合理的方法。我认为您所做的事情没有任何根本性的错误——它是一维函数二阶导数的一种特殊近似。
【讨论】：
对于一阶导数的不连续性，双梯度方法失败。 由于梯度函数将一个数据点向左和向右考虑在内，因此在多次应用时会继续/传播。
另一方面，二阶导数可以通过公式计算
d^2 f(x[i]) / dx^2 = (f(x[i-1]) - 2*f(x[i]) + f(x[i+1])) / h^2
比较here。这样做的好处是只需考虑两个相邻像素。
在图片中，比较了双np.gradient 方法（左）和上面提到的公式（右），由np.diff 实现。由于 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()
【讨论】：
当我一次又一次地以一种或另一种形式解决这个问题时，我决定编写一个函数gradient_n，它为np.gradient 添加了一个微分或其他功能。并非 np.gradient 的所有功能都受支持，例如多轴微分。
与np.gradient 一样，gradient_n 以与输入相同的形状返回微分结果。还支持像素距离参数 (d)。
import numpy as np
def gradient_n(arr, n, d=1, axis=0):
"""Differentiate np.ndarray n times.
Similar to np.diff, but additional support of pixel distance d
and padding of the result to the same shape as arr.
If n is even: np.diff is applied and the result is zero-padded
If n is odd:
np.diff is applied n-1 times and zero-padded.
Then gradient is applied. This ensures the right output shape.
"""
n2 = int((n // 2) * 2)
diff = arr
if n2 > 0:
a0 = max(0, axis)
a1 = max(0, arr.ndim-axis-1)
diff = np.diff(arr, n2, axis=axis) / d**n2
diff = np.pad(diff, tuple([(0,0)]*a0 + [(1,1)] +[(0,0)]*a1),
'constant', constant_values=0)
if n > n2:
assert n-n2 == 1, 'n={:f}, n2={:f}'.format(n, n2)
diff = np.gradient(diff, d, axis=axis)
return diff
def test_gradient_n():
import matplotlib.pyplot as plt
x = np.linspace(-4, 4, 17)
y = np.linspace(-2, 2, 9)
X, Y = np.meshgrid(x, y)
arr = np.abs(X)
arr_x = np.gradient(arr, .5, axis=1)
arr_x2 = gradient_n(arr, 1, .5, axis=1)
arr_xx = np.diff(arr, 2, axis=1) / .5**2
arr_xx = np.pad(arr_xx, ((0, 0), (1, 1)), 'constant', constant_values=0)
arr_xx2 = gradient_n(arr, 2, .5, axis=1)
assert np.sum(arr_x - arr_x2) == 0
assert np.sum(arr_xx - arr_xx2) == 0
fig, axs = plt.subplots(2, 2, figsize=(29, 21))
axs = np.array(axs).flatten()
ax = axs[0]
ax.set_title('x-cut')
ax.plot(x, arr[0, :], marker='o', label='arr')
ax.plot(x, arr_x[0, :], marker='o', label='arr_x')
ax.plot(x, arr_x2[0, :], marker='x', label='arr_x2', ls='--')
ax.plot(x, arr_xx[0, :], marker='o', label='arr_xx')
ax.plot(x, arr_xx2[0, :], marker='x', label='arr_xx2', ls='--')
ax.legend()
ax = axs[1]
ax.set_title('arr')
im = ax.imshow(arr, cmap='bwr')
cbar = ax.figure.colorbar(im, ax=ax, pad=.05)
ax = axs[2]
ax.set_title('arr_x')
im = ax.imshow(arr_x, cmap='bwr')
cbar = ax.figure.colorbar(im, ax=ax, pad=.05)
ax = axs[3]
ax.set_title('arr_xx')
im = ax.imshow(arr_xx, cmap='bwr')
cbar = ax.figure.colorbar(im, ax=ax, pad=.05)
test_gradient_n()
【讨论】：
这是原始文档的摘录（在撰写本文时可在http://docs.scipy.org/doc/numpy/reference/generated/numpy.gradient.html 找到）。它指出，除非采样距离为 1，否则您需要包含一个包含距离作为参数的列表。
numpy.gradient(f, *varargs, **kwargs)返回一个 N 维数组的梯度。
梯度是使用内部的二阶精确中心差和边界处的一阶精确单边（前向或后向）差异计算的。因此，返回的梯度与输入数组的形状相同。
参数：
f : array_like 包含标量函数样本的 N 维数组。varargs : 标量列表，可选 N 个标量，指定每个维度的样本距离，即 dx、dy、dz、... 默认距离：1。
edge_order : {1, 2}，可选 梯度是使用边界处的 N 阶精确差异计算的。默认值：1。 1.9.1 版中的新功能。
返回：
渐变：ndarray 与 f 形状相同的 N 个数组，给出 f 对每个维度的导数。
【讨论】：