对称函数的 numpy.gradient 的分量不同答案

【问题标题】：The components of numpy.gradient of a symmetric function are different对称函数的 numpy.gradient 的分量不同
【发布时间】：2019-02-26 18:13:17
【问题描述】：

对称函数的梯度应该在所有维度上都具有相同的导数。 numpy.gradient 提供不同的组件。

这是一个 MWE。

import numpy as np
x = (-1,0,1)
y = (-1,0,1)
X,Y = np.meshgrid(x,y)
f = 1/(X*X + Y*Y +1.0)

print(f)
>> [[0.33333333 0.5        0.33333333]
   [0.5        1.         0.5       ]
   [0.33333333 0.5        0.33333333]]

这在两个维度上具有相同的值。

但是np.gradient(f)给了

[array([[ 0.16666667,  0.5       ,  0.16666667],
    [ 0.        ,  0.        ,  0.        ],
    [-0.16666667, -0.5       , -0.16666667]]),

 array([[ 0.16666667,  0.        , -0.16666667],
    [ 0.5       ,  0.        , -0.5       ],
    [ 0.16666667,  0.        , -0.16666667]])]

渐变的两个分量是不同的。

为什么会这样？我在解释输出时遗漏了什么？

【问题讨论】：

因为导数是沿方向计算的，所以在比较它们时也需要考虑这个方向。沿 Y 轴查看的 Y 分量与沿 X 轴查看的 X 分量相同。即你可以通过将Y导函数顺时针旋转90度来看到原始函数的旋转对称性。
@meowgoesthedog 那么这样的3维数组的结果应该怎么看呢？
在 3D 中，您还必须在更改坐标系时保持轴的循环顺序，即在 XYZ 框架中查看 X 分量，在 YZX 中查看 Y，在 ZXY 中查看 Z .
@meowgoesthedog 是的，似乎是这样。我认为这应该在 numpy.gradient 文档中提及。
@meowgoesthedog 是否有任何方法可以撤消这种“循环”并获得 n 维数组的导数？

标签： python numpy gradient

【解决方案1】：

让我们一步一步来。首先，正如 meowgoesthedog 正确提到的那样 numpy 计算一个方向的导数。

Numpy 计算梯度的方式

需要注意的是，np.gradient 使用中心差异的含义（为简单起见，我们只关注一个方向）：

grad_f[i] = (f[i+1] - f[i])/2 + (f[i] - f[i-1])/2 =  (f[i+1] - f[i-1])/2

在边界处numpy计算（以min为例）

grad_f[min] = f[min+1] - f[min]
grad_f[max] = f[max] - f[max-1]

在您的情况下，边界是 0 和 2。

二维案例

如果您使用多个维度，我们需要将导数的方向考虑在内。 np.gradient 计算所有可能方向的导数。让我们重现您的结果：

让我们沿着列移动，因此我们使用行向量

进行计算

f[1,:] - f[0,:]

输出

array([0.16666667, 0.5       , 0.16666667])

这正是渐变的第一个元素的第一行。

该行是使用中心导数计算的，因此：

(f[2,:]-f[1,:])/2 + (f[1,:]-f[0,:])/2

输出

array([0., 0., 0.])

第三行：

f[2,:] - f[1,:]

输出

array([-0.16666667, -0.5       , -0.16666667])

对于另一个方向，只需交换: 和数字并记住您现在正在计算列向量。在对称函数的情况下，这直接导致转置导数，就像在你的情况下一样。

3D案例

x_ = (-1,0,4)
y_ = (-3,0,1)
z_ = (-1,0,12)

x, y, z = np.meshgrid(x_, y_, z_, indexing='ij')
f = 1/(x**2 + y**2 + z**2 + 1)
np.gradient(f)[1]

输出

array([[[ *2.50000000e-01,  4.09090909e-01,  3.97702165e-04*],
        [ 8.33333333e-02,  1.21212121e-01,  1.75554093e-04],
        [-8.33333333e-02, -1.66666667e-01, -4.65939801e-05]],

       [[ **4.09090909e-01,  9.00000000e-01,  4.03045231e-04**],
        [ 1.21212121e-01,  2.00000000e-01,  1.77904287e-04],
        [-1.66666667e-01, -5.00000000e-01, -4.72366556e-05]],

       [[ ***1.85185185e-02,  2.03619910e-02,  3.28827183e-04***],
        [ 7.79727096e-03,  8.54700855e-03,  1.45243282e-04],
        [-2.92397661e-03, -3.26797386e-03, -3.83406181e-05]]])

此处给出的梯度是沿行计算的（0 沿矩阵计算，1 沿行计算，2 沿列计算）。

这个可以计算出来

(f[:,1,:] - f[:,0,:])

输出

array([[*2.50000000e-01, 4.09090909e-01, 3.97702165e-04*],
       [**4.09090909e-01, 9.00000000e-01, 4.03045231e-04**],
       [***1.85185185e-02, 2.03619910e-02, 3.28827183e-04***]])

我添加了星号，以便清楚在哪里可以找到相应的行向量。由于我们计算了方向1 的梯度，我们必须寻找行向量。

如果想要重现整个渐变，这是由

完成的

np.stack(((f[:,1,:] - f[:,0,:]), (f[:,2,:] - f[:,0,:])/2, (f[:,2,:] - f[:,1,:])), axis=1)

n-dim 大小写

我们可以概括我们在这里学到的东西来计算任意函数沿方向的梯度。

def grad_along_axis(f, ax):
    f_grad_ind = []
    for i in range(f.shape[ax]):
        if i == 0:
            f_grad_ind.append(np.take(f, i+1, ax) - np.take(f, i, ax))
        elif i == f.shape[ax] -1:
            f_grad_ind.append(np.take(f, i, ax) - np.take(f, i-1, ax))
        else:
            f_grad_ind.append((np.take(f, i+1, ax) - np.take(f, i-1, ax))/2)
    f_grad = np.stack(f_grad_ind, axis=ax)
    return f_grad

在哪里

np.take(f, i, ax) = f[:,...,i,...,:]

而i 位于索引ax。

【讨论】：

【解决方案2】：

通常梯度和雅可比是函数上的运算符

如果你需要f = 1/(X*X + Y*Y +1.0) 的梯度，那么你必须象征性地计算它。或者使用使用该函数的数值方法来估计它。

我不知道常量 3d 数组的梯度是什么。 numpy.gradient 是一维概念。

Python 有 sympy 包，可以自动计算雅可比符号。

如果second order derivative of a scalar 3d field 指的是拉普拉斯算子，那么您可以使用标准的 4 点模板进行估算。

【讨论】：