我有一个3D数据数组A，形状分别为x、y和z。
我想求A在y维度上的梯度。这可以通过NumPy轻松实现：
dAdy = np.gradient(A, Y, axis=1)
其中Y是y维度中坐标的1D向量。
然而，如果Y是非结构化的，这就变得不平凡了。也就是说，固定位置(x, z) = (Xi, Zi)上的每一列数据都有一组唯一的y坐标。举例来说：

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)

上面的结果是一个3D数据集A，定义在X和Z坐标的结构化集合上，而Y坐标的值对于每个数据点都是唯一的（当然在y维度上是单调的）。我想通过有限差分来估计dA/dy。
本质上，我是在尝试取许多独立列的梯度。有没有一种方法可以用NumPy将其向量化？我尝试了下面的迭代方法，但它非常慢：

# 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?

我只得到了几个简单线性函数的grad1=grad2，而不是上面表示的函数。我现在想知道是否有一个理论上的原因，为什么链式法则不应该普遍适用于由有限差分估计的导数。