5

我正在尝试从此矩阵乘法中删除循环(并了解有关优化代码的更多信息),并且我认为我需要某种形式的np.broadcastingor np.einsum,但是在阅读了它们之后,我仍然不确定如何将它们用于我的问题。

A = np.array([[1, 2, 3, 4, 5],
         [6, 7, 8, 9, 10],
         [11,12,13,14,15]])
#A is a 3x5 matrix, such that the shape of A is (3, 5) (and A[0] is (5,))

B = np.array([[1,0,0],
         [0,2,0],
         [0,0,3]])
#B is a 3x3 (diagonal) matrix, with a shape of (3, 3)

C = np.zeros(5)
for i in range(5):
    C[i] = np.linalg.multi_dot([A[:,i].T, B, A[:,i]])

#Each row of matrix math is [1x3]*[3x3]*[3x1] to become a scaler value in each row
#C becomes a [5x1] matrix with a shape of (5,)

我知道我不能自己做np.multidot,因为这会产生一个 (5,5) 数组。

我还发现了这一点:将矩阵乘以 Numpy 中另一个矩阵的每一行,但我不知道它是否与我的问题相同。

4

3 回答 3

4
In [601]: C
Out[601]: array([436., 534., 644., 766., 900.])

这是很自然的einsum。我i和你一样使用,来表示结果的索引。 j并且k是在产品总和中使用的索引。

In [602]: np.einsum('ji,jk,ki->i',A,B,A)
Out[602]: array([436, 534, 644, 766, 900])

它可能也可以用 来完成mutmul,尽管它可能需要添加维度和后期挤压。

dot使用的方法diag做的工作比必要的要多得多。diag抛出了很多值。

要使用matmul,我们必须使i维度成为 3d 数组的第一个。这就是“被动”的结果:

In [603]: A.T[:,None,:]@B@A.T[:,:,None]
Out[603]: 
array([[[436]],     # (5,1,1) result

       [[534]],

       [[644]],

       [[766]],

       [[900]]])
In [604]: (A.T[:,None,:]@B@A.T[:,:,None]).squeeze()
Out[604]: array([436, 534, 644, 766, 900])

或者索引额外的维度:(A.T[:,None,:]@B@A.T[:,:,None])[:,0,0]

于 2018-12-07T20:53:50.203 回答
1

您可以将调用链接dot在一起,然后得到对角线:

# your original output:
# >>> C
# array([436., 534., 644., 766., 900.])

>>> np.diag(np.dot(np.dot(A.T,B), A))
array([436, 534, 644, 766, 900])

或者等效地,使用您最初multi_dot的思路,但采用生成的 5x5 数组的对角线。这可能会提高一些性能(根据文档

>>> np.diag(np.linalg.multi_dot([A.T, B, A]))
array([436, 534, 644, 766, 900])
于 2018-12-07T20:43:45.830 回答
0

在添加到答案。如果您想将矩阵相乘,您可以使用广播。编辑:注意这是元素乘法,而不是点积。为此,您可以使用点方法。

 B [...,None] * A

给出:

array([[[ 1,  2,  3,  4,  5],
        [ 0,  0,  0,  0,  0],
        [ 0,  0,  0,  0,  0]],

       [[ 0,  0,  0,  0,  0],
        [12, 14, 16, 18, 20],
        [ 0,  0,  0,  0,  0]],

       [[ 0,  0,  0,  0,  0],
        [ 0,  0,  0,  0,  0],
        [33, 36, 39, 42, 45]]])
于 2018-12-07T21:05:17.417 回答