Python - 矩阵外积

Question

给定两个矩阵

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 i in range (A.shape(0)):
      for j in range (B.shape(0)):
         C_ij = np.outer(A_i, B_j)

我想知道有没有更快的方法来做这个计算？

Answer 1

使用 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]]])

Answer 2

你可以使用

C = numpy.tensordot(A, B, axes=0)

tensordot 正好能满足你的需求。axes 参数的作用是可以在某些维度上进行求和（对于超过2维的张量，默认值是 2，这会把每个数组的两个维度进行求和），但如果把它设置为 0，那么就不会减少维度，而是保留整个外积。

Answer 3

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

因为你想要的结果是 C_ij = A_i * B_j，这可以通过对列向量A和行向量B进行逐元素相乘，利用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]]])]

Answer 4

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

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]]]])

Answer 5

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

NumPy的ufuncs，比如说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，让它变成你想要的顺序。