矩阵乘法与 ndarray

Question

如何执行向量和矩阵乘法 $\bf v^T \bf A v$ 其中 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 阵列相同的内容。对于三个元素的向量中的每一个，输出似乎也将减少为标量，即输出将是二维的。因此，为了解决您的问题，我们需要对第一个使用张量乘法：V.T @ A 对第一个轴求和，从而为我们提供一个 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 维矩阵乘法或多或少类似于二维乘法。除了 2 个维度之外，您的矩阵必须具有相同的形状。在这两个维度中，您应该遵守二维乘法的规则。例如，如果您有一个形状为 (a,b,c, ...,d,e,f) 的矩阵 A，并且您想将其与矩阵 B 相乘，则 B 的形状应为 (a,b,c,...,d,f,g)，结果的形状将为 (a,b,c, ...,d,e,g) .

让我们忘记我们处于 4D space。如果只有一个点，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))

returns :

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)

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

注意 1 理论上，除了最后一个维度，你可以乘以其他维度，只要你尊重这些维度的二维乘法规则，其他维度相同方面。但这是 Python.

中的方法

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

矩阵乘法与 ndarray

matrix multiplication with ndarray

python

arrays

numpy

matrix

linear-algebra