用numpy

如何解决用numpy

我正在尝试计算坐标转换（从球坐标到笛卡尔坐标）的协方差矩阵。该数据集包含504个不同的球坐标。

我通过使用矩阵乘积实现的误差传播来计算变换的协方差矩阵。

此矩阵乘积为A dot CR dot B。这些矩阵的一般形状可以在这里看到：

A = [[dx/dr,dx/dt,dx/dp]
     [dy/dr,dy/dt,dy/dp]
     [dz/dr,dz/dt,dz/dp]]

B = [[dx/dr,dy/dr,dz/dr]
     [dx/dt,dz/dt]
     [dx/dp,dy/dp,dz/dp]]

CR = [[0.001 ** 2,0]
      [0,0.003 ** 2,0.005 ** 2]]

现在，出于性能原因，我想以不用于循环的方式来实现此目的。 实际代码中的A和B（在文章的底部）是包含504个3x3矩阵的数组。 因此，阵列的形状为（508、3、3）。我想做的是对所有对应的矩阵A[0] dot CR dot B[0],A[0] dot CR dot B[0],...采用上述的点积，并使所得数组的形状为（508，3，3）。我目前已使用列表推导（请参阅CX的计算）来实现此功能，但我想完全删除for循环。

我尝试将（508，3，3）数组传递给np.dot()-> CX = np.dot(np.dot(A,CR),B)，但这导致（508，3，508，3）数组不是我想要的。似乎没有关键字参数可以指定点积应使用哪个轴，例如np.transpose()。有人知道我该怎么实现吗？

SIGMA_R,SIGMA_THETA,SIGMA_PHI = 0.001,0.003,0.005
CR = np.array([[SIGMA_R ** 2,0],[0,SIGMA_THETA ** 2,SIGMA_PHI ** 2]])

A = np.transpose(np.array(([derivative_r(r,theta,phi),derivative_theta(r,derivative_phi(r,phi)])))

B = np.transpose(A,axes=(0,2,1))
CX = np.array([np.dot(np.dot(A[i],B[i]) for i in range(504)])

解决方法

您似乎需要np.einsum

CX = np.einsum('ijk,kl,ilm -> ijm',A,CR,B) 

甚至：

CX = np.einsum('ijk,iml -> ijm',A)