矩阵乘法与 ndarray

Question

如何执行向量和矩阵乘法 其中 v 是一个向量 (v1, v2, v3) 而 A 是一个 3x3 矩阵？ Python 抱怨形状不对齐，可能是因为 v 是一个 ndarray。关于如何进行此操作的任何想法？最终结果应该是每个坐标点 (v1, v2, v3) 的标量。下面的基本代码在尝试进行乘法运算时会中断。

import numpy as np
a = np.linspace(0, 10, 21)
b = np.linspace(0, 20, 41)
a, b = np.meshgrid(a,b)
v = np.array([a*b, a+b, a])
A = np.ones([3,3])
s = v.T @ A @ v     # doesn't work

错误

----> 1 s = v.T @ A @ v    
ValueError: shapes (21,41,3) and (3,41,21) not aligned: 3 (dim 2) != 41 (dim 1)

编辑：矩阵运算应该在每个点v完成，其中v通常是一个大数组（向量）。例如取一个 1m 的立方体，其中心位于原点，并在每个网格点计算矩阵运算，例如在每个坐标轴上每隔 10cm。

编辑 2 在 (x,y,z) 处的单个点的示例

A = np.zeros([3,3])
A[0][0] = 1
A[1][1] = 2
A[2][2] = 3
x,y,z = 1, 1, 0
v = np.array([x, y, z])
s = v.T @ A @ v   # should give s=3

下一步是使代码适用于大量向量v。除了它有点复杂，因为矢量坐标 (x,y,z) 需要根据坐标 (a,b) 进行参数化。上面的原始代码试图做到这一点，但不起作用，可能不是最好的方法。还有其他想法吗？

Answer 1

似乎提到三个元素的向量v，你的意思是一个ndarray，它的第一个轴上有三个元素，每个元素都保存一个n-dim数组数据。对于列出的示例，您拥有与 3D 数组相同的内容。 对于三个元素的向量中的每一个，输出似乎也将减少为一个标量，即输出将是 2D。 所以，为了解决你的情况，我们需要对第一个使用张量乘法：V.T @ A sum-reducing 第一个轴，给我们一个 3D 数组。然后，使用einsum 保持前两个轴对齐并减少最后一个轴，就像这样 -

p1 = np.tensordot(v,A,axes=((0,0)))
out = np.einsum('jkl,ljk->jk',p1,v)

或者，使用einsum，我们可以一步完成所有操作，就像这样 -

out = np.einsum('ijk,il,ljk->jk',v,A,v)

我们可以通过 einsum's 可选参数使其更快：optimize 设置为 True：

np.einsum(..., optimize=True)

Answer 2

当您将两个 N 维矩阵与 numpy 相乘时，我想它会自动将最后两个维度相乘并保留第一个维度。 N 维矩阵乘法或多或少类似于 2D 乘法。除了二维之外，您的矩阵必须具有相同的形状。在这 2 个维度中，您应该遵守 2D 乘法的规则。例如，如果您有一个形状为 (a,b,c, ...,d,e,f) 的矩阵 A，并且您想将其与矩阵 B 相乘，则 B 的形状应为 (a,b,c,...,d,f,g)，结果的形状为 (a,b,c, ...,d,e,g)。

让我们忘记我们身处 4D 空间。如果你只有一个点，v^T*A*v 的形状应该是(1,3)x(3,3)x(3,1)。我们只是想将其应用于(41,21) 网格中的每个点。这为我们提供了需要相乘的每个组件的最后尺寸。为了保持一致，v^T*A*v 的形状应为(41,21,1,3)x(3,3)x(41,21,3,1)。

 import numpy as np
 a = np.linspace(0, 10, 21)
 b = np.linspace(0, 20, 41)
 a, b = np.meshgrid(a,b)
 a = np.expand_dims(a, axis=0)
 b = np.expand_dims(b, axis=0)
 print("Shape a = {}, Shape b = {}".format(a.shape, b.shape))
 v = np.array([a*b, a+b, a])
 print("Shape v = {}".format(v.shape))
 u1 = v.transpose((2,3,1,0))
 print("Shape u1 = {}".format(u1.shape))
 s = u1 @ A
 u2 = v.transpose((2,3,0,1))
 print("Shape u2 = {}".format(u2.shape))
 s = s @ u2
 print("{} x {} x {} = {} x {} = {}".format(u1.shape, A.shape, u2.shape, (u1 @ A).shape, u2.shape, s.shape))

返回：

Shape a = (1, 41, 21), Shape b = (1, 41, 21)
Shape v = (3, 1, 41, 21)
Shape u1 = (41, 21, 1, 3)
Shape u2 = (41, 21, 3, 1)
(41, 21, 1, 3) x (3, 3) x (41, 21, 3, 1) = (41, 21, 1, 3) x (41, 21, 3, 1) = (41, 21, 1, 1)

我向您推荐这个解决方案。首先将大小为 1 的维度添加到向量 a 和 b。它们的形状不是(41,21)，而是(1,41,21)。现在，当您构造v 时，您将获得(3,1,41,21) 的形状。现在，如果您使用通常的转置，您只需反转所有维度，这不是您想要的。您希望 v^T 可乘以 A，形状为 (3,3)。因此，您可以手动定义如何将向量的维度从(3,1,41,21) 转换为(41,21,1,3) 和(41,21,3,1)。最后，终于可以相乘了，而且是一致的。

注意 1 在理论上，您可以将除最后一个维度之外的其他维度相乘，只要您遵守这些维度的 2D 乘法规则，其他维度相同即可。但这是在 Python 中实现的方式。

注意 2 您可能想知道为什么我们可以将形状为 (41,21,1,3) 的矩阵与形状为 (3,3) 的矩阵相乘。这与将 2D 矩阵乘以标量时的机制完全相同。当你这样做时，你将标量的维度增加到二维（基本上是一个到处都有标量的矩阵），然后你执行逐元素乘法。同样，您创建一个形状为(41,21,3,3) 的矩阵，然后逐元素相乘，或“逐块”相乘（二维矩阵相乘）。元素给出乘法(1,3)x(3,3)。