有没有办法用 Numpy 计算非结构化坐标上的梯度？

Question

我有一个 3D 数据数组

A

，形状分别为

(NX, NY, NZ)

、

x

和

y

维度。我想求

z

在

A

维度上的梯度。这可以使用 NumPy 轻松完成：

y

其中

dAdy = np.gradient(A, Y, axis=1)

是

Y

维度中的一维坐标向量。但是，如果

y

是非结构化的，这就变得很重要。也就是说，固定位置

Y

处的每个数据“列”都有一组唯一的

(x, z) = (Xi, Zi)

坐标。例如：y

上面的结果是一个 3D 数据集 

A = np.random.random((10, 10, 10))

X = np.arange(10)
Y = np.sort(np.random.random((10, 10, 10)), axis=1)
Z = np.arange(10)

，定义在一组结构化的

A

和

X

坐标上，而

Z

坐标的值对于每个数据点都是唯一的（但在

Y

尺寸）。我想通过有限差分来估计

y

。本质上，我正在尝试获取许多独立列的梯度。有没有办法用 NumPy 对其进行矢量化？我尝试了以下迭代方法，但速度非常慢：

dA/dy

我还认为通过实施链式法则我可以变得聪明：

# A is the 3D dataset
# Y is the 3D dataset with shape matching that of A; gives the y-position of each datapoint in A
NX, NY, NZ = A.shape[0], A.shape[1], A.shape[2]
dA_dy = np.zeros((NX, NY, NZ))
for i in range(NX):
    for k in range(NZ):
        dA_dy[i, :, k] = np.gradient(A[i,:,k], Y[i,:,k])

但是对于下面的简单测试，这种方法不起作用：

dA_dy = np.gradient(A, axis=1) / np.gradient(Y, axis=1)

我只得到了一些简单的线性函数的

g = np.array([1, 5, 6, 10])  # an unstructured coordinate
f = g**2                     # function value on the points x
grad1 = np.gradient(f, g)                # df/dg
grad2 = np.gradient(f) / np.gradient(g)  # df/dg?

，而不是上面表示的函数。我现在想知道是否存在理论上的原因为什么链式法则对于由有限差分估计的导数通常不成立。

Answer 1

（不是解决问题的答案）

对于一些简单的线性函数，我只得到 grad1=grad2

当然：

grad1=grad2

如果 

# np.gradient(f) is equivalent to:
>>> np.gradient(f, np.arange(f.size))
array([24. , 17.5, 37.5, 64. ])

# np.gradient(x) is equivalent to:
>>> np.gradient(x, np.arange(x.size))
array([4. , 2.5, 2.5, 4. ])

# so np.gradient(f) / np.gradient(x) is equivalent to:
>>> np.gradient(f, np.arange(f.size)) / np.gradient(x, np.arange(f.size))
array([ 6.,  7., 15., 16.])

均匀分布，则

x

等于

grad1

，即使函数

grad2

不是线性的：f

输出：

x = np.array([1, 3, 5, 7])
f = x**2
grad1 = np.gradient(f, x)
grad2 = np.gradient(f) / np.gradient(x)