python中不同形状数组的元素乘法答案

【问题标题】：Elementwise multiplication of arrays of different shapes in pythonpython中不同形状数组的元素乘法
【发布时间】：2012-05-20 21:53:55
【问题描述】：

假设我有两个数组 a 和 b，

  a.shape = (5,2,3)
  b.shape = (2,3)

然后c = a * b 会给我一个数组c，形状为(5,2,3) 和c[i,j,k] = a[i,j,k]*b[j,k]。

现在的情况是，

  a.shape = (5,2,3)
  b.shape = (2,3,8)

我希望c 具有(5,2,3,8) 和c[i,j,k,l] = a[i,j,k]*b[j,k,l] 的形状。

如何有效地做到这一点？我的a 和b 实际上很大。

【问题讨论】：

标签： python arrays numpy multiplication

【解决方案1】：

这应该可行：

a[..., numpy.newaxis] * b[numpy.newaxis, ...]

用法：

In : a = numpy.random.randn(5,2,3)

In : b = numpy.random.randn(2,3,8)

In : c = a[..., numpy.newaxis]*b[numpy.newaxis, ...]

In : c.shape
Out: (5, 2, 3, 8)

参考：Array Broadcasting in numpy

编辑：更新参考网址

【讨论】：

我发现None也可以，除了numpy.newaxis；事实上，numpy.newaxis is None 对我来说是正确的。有理由不在这里使用None 吗？
@senderle。显然，它们是可以互换的：numpy.newaxis。但是，将来它可能会改变。我想这就是为什么他们把newaxis 作为None 的同义词放在那里。
np.newaxis is None 返回True。我倾向于使用newaxis 只是为了在我的代码中明确说明。另见stackoverflow.com/questions/944863/…
@Avaris 如果你对我的评论表示赞赏，我将不胜感激。

【解决方案2】：

我认为以下应该可行：

import numpy as np
a = np.random.normal(size=(5,2,3))
b = np.random.normal(size=(2,3,8))
c = np.einsum('ijk,jkl->ijkl',a,b)

和：

In [5]: c.shape
Out[5]: (5, 2, 3, 8)

In [6]: a[0,0,1]*b[0,1,2]
Out[6]: -0.041308376453821738

In [7]: c[0,0,1,2]
Out[7]: -0.041308376453821738

np.einsum 使用起来可能有点棘手，但对于这类索引问题来说非常强大：

http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html

另请注意，这需要 numpy >= v1.6.0

我不确定您的特定问题的效率，但如果它的性能不如所需，请务必考虑使用带有显式 for 循环的 Cython，并可能使用 prange 并行化它

更新

In [18]: %timeit np.einsum('ijk,jkl->ijkl',a,b)
100000 loops, best of 3: 4.78 us per loop

In [19]: %timeit a[..., np.newaxis]*b[np.newaxis, ...]
100000 loops, best of 3: 12.2 us per loop


In [20]: a = np.random.normal(size=(50,20,30))
In [21]: b = np.random.normal(size=(20,30,80))

In [22]: %timeit np.einsum('ijk,jkl->ijkl',a,b)
100 loops, best of 3: 16.6 ms per loop

In [23]: %timeit a[..., np.newaxis]*b[np.newaxis, ...]
100 loops, best of 3: 16.6 ms per loop

In [2]: a = np.random.normal(size=(500,20,30))
In [3]: b = np.random.normal(size=(20,30,800))

In [4]: %timeit np.einsum('ijk,jkl->ijkl',a,b)
1 loops, best of 3: 3.31 s per loop

In [5]: %timeit a[..., np.newaxis]*b[np.newaxis, ...]
1 loops, best of 3: 2.6 s per loop

【讨论】：

嗯，好吧...对于更大的数组numpy.einsum 似乎有点慢。
@Avaris 我一般同意您对时间的评估，尽管在相当大的数组大小范围内，np.einsum 与广播解决方案一样快或更快。在某些方面，我更喜欢 einsum 解决方案的可读性和明确性，尽管其他人可能不同意。但可以肯定的是，如果对于用例，广播速度更快，请使用它。