将 (a x a) 方阵乘以 (a x c x d) 矩阵

Question

假设我有一个(c x d) 矩阵列表。比如说我有a。每个矩阵都有a 系数。

在 NumPy 中是否有一种快速的方法可以一次将每个矩阵乘以它的系数，同时仍然保持张量数据结构，或者我是否需要在 for 循环中手动完成，即X = np.array([np.multiply(coefs[i], X[i]) for i in range(len(coefs))])

即X.shape = (3, 4, 5), coefs.shape = (3).

Answer 1

X = np.array([[[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]],
              [[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]],
              [[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]]])

coeffs = np.array([2,4,6])

您需要将坐标轴添加到coeffs，以便它将在您想要的维度中广播。

>>> X * coeffs[:, np.newaxis, np.newaxis]
array([[[2, 2, 2, 2, 2],
        [2, 2, 2, 2, 2],
        [2, 2, 2, 2, 2],
        [2, 2, 2, 2, 2]],

       [[4, 4, 4, 4, 4],
        [4, 4, 4, 4, 4],
        [4, 4, 4, 4, 4],
        [4, 4, 4, 4, 4]],

       [[6, 6, 6, 6, 6],
        [6, 6, 6, 6, 6],
        [6, 6, 6, 6, 6],
        [6, 6, 6, 6, 6]]])
>>> 

np.newaxis 允许coeffs 的值与X 的第一个维度对齐，然后它们在其余维度中为broadcast。

Answer 2

解决方案：

这是采用pythonic方式的经典问题。下面的 1-liner 使用强大的 einstein sum、np.einsum() 函数和丰富的下标标签框架来解决问题，这些框架允许 broadcasting 的维度和控制输出：

np.einsum('ijk,i...->ijk', X, coeffs)

下标字符串：

• 逗号, 将左侧第一个操作数的下标标签与右侧第二个操作数的下标标签分开。
• -> 符号后的下标标签给出了输出的尺寸。
• 省略号... 将允许将coeffs 向量广播到X 矩阵的额外2 维中（... 实际上对于更明确的广播规范是可选的）

这需要一点时间来适应，但有许多矩阵向量运算，如trace、diag、inner/outer product、element-wise multiplication 等，都集中在这个强大的公式中。

言下之意：

所以字符串本质上是说取X矩阵的所有维度并将其与coeffs向量逐元素相乘（在X的额外维度上广播）并生成一个相同维度的输出矩阵和X一样。

输出：

>>> np.einsum( 'ijk,i->ijk', np.ones((2,2,2)), np.array([2,4]))

array([[[ 2.,  2.],
        [ 2.,  2.]],

       [[ 4.,  4.],
        [ 4.,  4.]]])

在二维情况下，沿任一轴逐元素：

>>> np.einsum( 'ij,i->ij', np.ones((2,2)), np.array([2,4]))
array([[ 2.,  2.],
       [ 4.,  4.]]) 

>>> np.einsum( 'ij,j->ij', np.ones((2,2)), np.array([2,4]))
array([[ 2.,  4.],
       [ 2.,  4.]])

替代方案：

当然，您可以使用np.newaxis 参数和简单的乘法在额外维度上手动广播向量：

X * coeff[:, np.newaxis, np.newaxis]

时间：

使用 ipython line magic 函数%timeit 我们实际上可以看到einsum 对于大型数据集可能会贵一点：

>>> %timeit np.einsum('ijk,i...->ijk',np.ones((10,100,100)),np.ones(10))
1000 loops, best of 3: 209 µs per loop

>>> %timeit np.ones((10,100,100))*np.ones(10)[:,np.newaxis, np.newaxis]
1000 loops, best of 3: 129 µs per loop

然后权衡einstein sum 提供的多功能性。

P.S.如果你是一个痴迷于它提供的丰富索引框架的 Matlab 用户，你可能有兴趣查看这个 Numpy-Matlab comparison 页面