特定维度矩阵的乘法

Question

假设我有两个矩阵，A 和 B。

A 是一个三维数组，或者说是一个张量。

[1,2,3,4]
[5,6,7,8]
[1,2,3,4]
[5,6,7,8]

在 A 的第三个维度上，有 4 个不同的二维矩阵，像上面那样。

B 是一个矩阵。

[1,2,3,4]

B 中也有 4 个这样的矩阵，每一个都是相同的。

我想知道如何把 B 中的每一个向量和 A 中的每个二维矩阵相乘。

[1,2,3,4]*[1,2,3,4]*[1;2;3;4]
          [5,6,7,8]
          [1,2,3,4]
          [5,6,7,8]

这样会有四次这样的乘法，得到四个 4x1 的向量。我尝试用 numpy 来实现：

y = numpy.arange(4).reshape(1,4)
z = numpy.arange(64).reshape(4,4,4)
y.dot(z).dot(numpy.transpose(y))
------
Output:
array([[[ 420],
        [ 996],
        [1572],
        [2148]]])

这样做的效果正是我想要的。但我对 numpy 的广播机制并不太了解，我想学习一下这个知识点。此外，不同的库在处理矩阵时，广播的方式也可能不同。我尝试用不同的方法来调整 B，以达到相同的结果，但都没有成功。如果我解释得不够清楚，请告诉我。

另外，我更希望从 numpy 得到 4x1 的结果，而不是三维的返回值。

Answer 1

为了保持维度的清晰，我会使用变量。在测试时，我喜欢在每个维度使用不同的大小（比如 z=...reshape(3,4,5)），这样错误就能很明显地显示出来。

k,j,l = 4
y = numpy.arange(4).reshape(1,j)
z = numpy.arange(64).reshape(k,j,l)

y.dot(z) 会把 y 的最后一个维度和 z 的倒数第二个维度（j）结合起来（可以查看它的文档）。结果会是 (1,k,l)。

np.einsum 可以帮助确认这一点：

print y.dot(z)
print np.einsum('ij,kjl',y,z) # note the repeated j

这两者的结果都是

array([[[ 56,  62,  68,  74],
        [152, 158, 164, 170],
        [248, 254, 260, 266],
        [344, 350, 356, 362]]])

y.T（转置）是 (j,1)。所以第二次点乘会把 (1,k,l) 和 (j,1) 结合，最后一个维度（l）和第一个维度（j）结合，结果是 (1,k,1)，这就是你的三维数组：

array([[[ 420],
        [ 996],
        [1572],
        [2148]]])

np.einsum('ijl,lm', np.einsum('ij,kjl',y,z), y.T) # note the repeated l

这个 einsum 可以合并成一个调用（使用相同的索引）：

np.einsum('ij,kjl,ml',y,z,y)

这里没有涉及广播。也就是说，没有维度被添加或扩展。你只需要系统地跟踪维度在 dot 里是如何变化的。

如果想得到一个维度更小的结果，你需要把 1 的维度挤掉或者重新调整结果的形状。如果你一开始用的是一维的 y，那么结果也是一维的。

y1=y.squeeze()
np.einsum('j,kjl,l',y1,z,y1)
y1.dot(z).dot(y1.T)
# array([ 420,  996, 1572, 2148])

Answer 2

注意，np.dot在处理不同维度的数组时表现是不同的，具体可以参考这里：

对于二维数组，它相当于矩阵乘法。

对于一维数组，它是向量的内积（不涉及复共轭）。

对于N维数组，它是在数组a的最后一个维度和数组b的倒数第二个维度上进行求和乘积。

在你的例子中，z有三维。

另外，还有一个重要的区别在于numpy数组和numpy矩阵：

你在例子中使用的是数组。Numpy矩阵只有二维，矩阵乘法是通过*运算符来实现的。（而在数组中，这个运算符是逐元素相乘）

Answer 3

你可以看看 np.einsum 这个功能。举个例子：

>>> mat = np.arange(80).reshape(4, 4, 5)
>>> vec = np.arange(12).reshape(3, 4)
>>> np.einsum('ij,jkl,ik->il', vec, mat, vec)
array([[ 2100,  2136,  2172,  2208,  2244],
       [20900, 21384, 21868, 22352, 22836],
       [58900, 60344, 61788, 63232, 64676]])

如果我没有搞错索引的话，我有5个形状为4x4的矩阵，还有3个长度为4的向量，我正在计算每个矩阵和每个向量之间的3x5的二次形式。