在 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 只能处理二维数组。在 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 方法是一致的。

您可能需要一段时间才能理解爱因斯坦求和，但它非常棒且功能强大。我强烈建议您进一步了解它，特别是如果这种线性代数对您来说是一个常见问题。