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没有实现这个功能。

补充说明1：我知道有np.diff(a, 2)。但这只是单边估计，所以我很好奇为什么np.gradient没有类似的选项。

补充说明2：np.sin()只是个玩具数据，真实的数据没有解析形式。

谢谢！

Answer 1

以下内容摘自原始文档（在写作时可以在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中新增。

返回值：
gradient : ndarray
返回N个数组，形状和f一样，表示f在每个维度上的导数。

Answer 2

因为我总是遇到这个问题，所以我决定写一个叫做 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()

Answer 3

双重梯度方法在处理一阶导数的不连续性时会出现问题。因为这个梯度函数会考虑左边和右边的一个数据点，所以当多次应用时，这种影响会继续扩散。

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

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

你可以在这里查看。这种方法的好处是只考虑两个相邻的像素。

在图片中，左边是双重np.gradient方法，右边是前面提到的公式，使用np.diff实现。因为f(x)在零点只有一个拐点，所以二阶导数（绿色）应该只在那个地方有一个峰值。由于双重梯度方法考虑了每个方向上的两个相邻点，这导致在+/- 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 4

我同意@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年的评论中指出的，这个答案（至少适用于Numpy 1.08）在Numpy 1.14及之后的版本中就不再适用了。请检查你的版本号和可用的选项。

Answer 5

对于数值梯度的计算，并没有一个放之四海而皆准的正确答案。在你计算样本数据的梯度之前，得先对生成这些数据的底层函数做一些假设。其实你可以用 np.diff 来计算梯度，使用 np.gradient 也是一个不错的方法。我觉得你现在的做法没有什么根本上的问题——这只是对一维函数的二阶导数的一种特定近似。