将矩阵与NumPy中的向量数组相乘

Question

我正在尝试以几何方式旋转NumPy中的向量数组。首先，我生成网格的坐标矢量。

width = 128
height = 128

x_axis = np.linspace(-1, 1, width)
y_axis = np.linspace(-1, 1, height)

x, y = np.meshgrid(x_axis, y_axis)
z = np.full((width, height), 0)
vectors = np.stack((x, y, z), axis=2)

所以'向量'的形状为（128，128，3）

我已经准备了旋转矩阵，其中a，b和c为沿轴的旋转角度。

rotation_matrix = np.array([
    [np.cos(b) * np.cos(c),
     - np.cos(b) * np.sin(c),
     np.sin(b)],

    [np.sin(a) * np.sin(b) * np.cos(c) + np.cos(a) * np.sin(c),
     - np.sin(a) * np.sin(b) * np.sin(c) + np.cos(a) * np.cos(c),
     - np.sin(a) * np.cos(b)],

    [- np.cos(a) * np.sin(b) * np.cos(c) + np.sin(a) * np.sin(c),
     np.cos(a) * np.sin(b) * np.sin(c) + np.sin(a) * np.cos(c),
     np.cos(a) * np.cos(b)]
])

现在，我希望将数组的每个向量都矩阵乘以'rotation_matrix'，

vector_rotated = rotation_matrix @ vector

因此，所得数组也应具有形状（128、128、3）。我在处理此3维数组时遇到一些问题。 Matmul仅能处理2d数组。 NumPy中是否有用于此用例的优雅方法，还是我必须使用for循环来解决此问题？

非常感谢，祝您有美好的一天！

Answer 1

有几种不同的方法可以解决此问题。

选项1：

最直接的方法是对数组vectors进行整形，使其具有形状(3, 128 * 128)，然后调用内置的np.dot函数，并将结果整形为所需的形状。

（（请注意，阵列形状的(128, 128)部分与旋转实际上无关；这是一种解释，您可能想使问题更清楚，但对要应用的线性变换没有影响。方式是旋转3个向量，其中有128 * 128 == 16384个，它们恰好像上面一样被组织成3D数组。）

这种方法看起来像：

>>> v = vectors.reshape(-1, 3).T
>>> np.dot(rotation_matrix, v).shape
(3, 16384)
>>> rotated = np.dot(rotation_matrix, v).T.reshape(vectors.shape)
>>> rotated.shape == vectors.shape
True

选项2：

不涉及任何重塑的另一种方法是使用NumPy的Einstein summation。爱因斯坦求和是非常灵活的，需要一段时间才能理解，但是其力量证明了其复杂性。以其最简单的形式，“标记”您要相乘在一起的轴。省略的轴是“收缩的”，这意味着将计算该轴上的和。对于您的情况，它将是：

>>> np.einsum('ij,klj->kli', rotation_matrix, vectors).shape
(128, 128, 3)
>>> np.allclose(rotated, np.einsum('ij,klj->kli', rotation_matrix_vectors))
True

这里是对索引的快速说明。我们标记旋转矩阵i和j的轴，以及矢量k，l和j的轴。重复的j表示那些是相乘的轴。这等效于将上面的变形数组与旋转矩阵右乘（即，是旋转）。

输出轴标记为kli。这意味着我们将保留矢量的k和l轴。由于j不在输出标签中，因此该轴上存在一个求和。相反，我们有轴i，因此最终形状为(128, 128, 3)。您可以在上方看到点积方法和einsum方法一致。

将您的头放在爱因斯坦的总和上可能要花一些时间，但是它很棒而且功能强大。我强烈建议您了解更多有关它的信息，特别是如果这种线性代数是您的常见问题。