【问题标题】：Python - Matrix outer productPython - 矩阵外积
【发布时间】：2014-09-10 10:42:41
【问题描述】：

给定两个矩阵

A: m * r
B: n * r

我想生成另一个矩阵C: m * n，每个条目C_ij 是由A_i 和B_j 的外积计算的矩阵。

例如，

A: [[1, 2],
    [3, 4]]

B: [[3, 1],
    [1, 2]]

给了

C: [[[3, 1],  [[1 ,2],
     [6, 2]],  [2 ,4]],
     [9, 3],  [[3, 6],
     [12,4]],  [4, 8]]]

我可以使用 for 循环来做到这一点，比如

    for i in range (A.shape(0)):
      for j in range (B.shape(0)):
         C_ij = np.outer(A_i, B_j)

我想知道是否有一种矢量化的方式来做这个计算来加速它？

【问题讨论】：

您想要一个 4D、(m, n, r, r)-shape 数组，还是想要一个 2D、(m, n)-shape object dtype 数组，其中每个元素都是另一个数组？我强烈推荐第一个选项，但您的描述听起来更接近第二个。
对不起，我更喜欢第一个，4D (m, n, r, r)-shape 数组。

标签： python numpy matrix

【解决方案1】：

使用 Numpy 数组广播的简单解决方案

既然你想要C_ij = A_i * B_j，这可以简单地通过对column-vector-A和row-vector-B的element-wise-product进行numpy广播来实现，如下所示：

# import numpy as np
# A = [[1, 2], [3, 4]]
# B = [[3, 1], [1, 2]]
A, B = np.array(A), np.array(B)
C = A.reshape(-1,1) * B.reshape(1,-1)
# same as: 
# C = np.einsum('i,j->ij', A.flatten(), B.flatten())
print(C)

输出：

array([[ 3,  1,  1,  2],
       [ 6,  2,  2,  4],
       [ 9,  3,  3,  6],
       [12,  4,  4,  8]])

然后您可以使用numpy.dsplit() 或numpy.array_split() 获得所需的四个子矩阵，如下所示：

np.dsplit(C.reshape(2, 2, 4), 2)
# same as:
# np.array_split(C.reshape(2,2,4), 2, axis=2)

输出：

[array([[[ 3,  1],
         [ 6,  2]],

        [[ 9,  3],
         [12,  4]]]), 
array([[[1, 2],
         [2, 4]],

        [[3, 6],
         [4, 8]]])]

【讨论】：

【解决方案2】：

爱因斯坦符号很好地表达了这个问题

In [85]: np.einsum('ac,bd->abcd',A,B)
Out[85]: 
array([[[[ 3,  1],
         [ 6,  2]],

        [[ 1,  2],
         [ 2,  4]]],


       [[[ 9,  3],
         [12,  4]],

        [[ 3,  6],
         [ 4,  8]]]])

【讨论】：

伙计，我应该学习那个符号。每次有人使用它发布答案时，它总是比我想出的要短得多。如果您知道爱因斯坦求和约定，可能也更容易理解。

【解决方案3】：

temp = numpy.multiply.outer(A, B)
C = numpy.swapaxes(temp, 1, 2)

NumPy ufunc，例如 multiply，有一个 outer 方法，几乎可以满足您的需求。以下：

temp = numpy.multiply.outer(A, B)

产生一个结果，如temp[a, b, c, d] == A[a, b] * B[c, d]。你想要C[a, b, c, d] == A[a, c] * B[b, d]。 swapaxes 调用重新排列 temp 以按您想要的顺序放置。

【讨论】：

【解决方案4】：

使用numpy；

In [1]: import numpy as np

In [2]: A = np.array([[1, 2], [3, 4]])

In [3]: B = np.array([[3, 1], [1, 2]])

In [4]: C = np.outer(A, B)

In [5]: C
Out[5]: 
array([[ 3,  1,  1,  2],
       [ 6,  2,  2,  4],
       [ 9,  3,  3,  6],
       [12,  4,  4,  8]])

一旦你得到想要的结果，你就可以使用numpy.reshape()将它塑造成你想要的几乎任何形状；

In [6]: C.reshape([4,2,2])
Out[6]: 
array([[[ 3,  1],
        [ 1,  2]],

       [[ 6,  2],
        [ 2,  4]],

       [[ 9,  3],
        [ 3,  6]],

       [[12,  4],
        [ 4,  8]]])

【讨论】：

我有点惊讶这在视觉上如此接近 OP 想要的，但尽管相似，它仍然不太正确。如果您将所需的结果可视化为一个大的 2D 网格，单元格中有一些小的 2D 网格，这就是您将所有的小网格在边缘粘合在一起形成一个网格的结果。
是的，同意你的看法。也许先做np.outer(A, B)，然后把它分成更小的网格？
但这会产生不同的结果。我希望内部矩阵为[[3, 1], [6, 2]]、[[1, 2], [2, 4]]、[[9, 3], [12, 4]] 和[[3, 6], [4, 8]]。