numpy：按列点积答案

【问题标题】：numpy: column-wise dot productnumpy：按列点积
【发布时间】：2011-09-07 22:35:55
【问题描述】：

给定一个二维numpy 数组，我需要计算每一列与其自身的点积，并将结果存储在一维数组中。以下作品：

In [45]: A = np.array([[1,2,3,4],[5,6,7,8]])

In [46]: np.array([np.dot(A[:,i], A[:,i]) for i in xrange(A.shape[1])])
Out[46]: array([26, 40, 58, 80])

有没有一种简单的方法可以避免 Python 循环？以上几乎不是世界末日，但如果有一个numpy 原语，我想使用它。

edit 实际上，矩阵有很多行和相对较少的列。因此，我并不太热衷于创建大于 O(A.shape[1]) 的临时数组。我也无法原地修改A。

【问题讨论】：

标签： python numpy dot-product

【解决方案1】：

怎么样：

>>> A = np.array([[1,2,3,4],[5,6,7,8]])
>>> (A*A).sum(axis=0)
array([26, 40, 58, 80])

编辑：嗯，好吧，你不想要中间大对象。也许：

>>> from numpy.core.umath_tests import inner1d
>>> A = np.array([[1,2,3,4],[5,6,7,8]])
>>> inner1d(A.T, A.T)
array([26, 40, 58, 80])

这似乎有点快。这应该在幕后做你想做的事，因为 A.T 是一个视图（它不会制作自己的副本，IIUC），而 inner1d 似乎会按照它需要的方式循环。

很晚的更新：另一种选择是使用np.einsum：

>>> A = np.array([[1,2,3,4],[5,6,7,8]])
>>> np.einsum('ij,ij->j', A, A)
array([26, 40, 58, 80])
>>> timeit np.einsum('ij,ij->j', A, A)
100000 loops, best of 3: 3.65 us per loop
>>> timeit inner1d(A.T, A.T)
100000 loops, best of 3: 5.02 us per loop
>>> A = np.random.randint(0, 100, (2, 100000))
>>> timeit np.einsum('ij,ij->j', A, A)
1000 loops, best of 3: 363 us per loop
>>> timeit inner1d(A.T, A.T)
1000 loops, best of 3: 848 us per loop
>>> (np.einsum('ij,ij->j', A, A) == inner1d(A.T, A.T)).all()
True

【讨论】：

inner1d 似乎只是门票。谢谢。

【解决方案2】：

您可以计算所有元素的平方并使用逐列求和

np.sum(np.square(A),0);

（我不完全确定 sum 函数的第二个参数，它标识了取总和的轴，我目前没有安装 numpy。也许你必须尝试一下 :) 。 ..)

编辑

看DSM的帖子，好像应该用axis=0。使用square 函数可能比使用A*A 性能更高。

【讨论】：

【解决方案3】：

从线性代数来看，第 i 行与第 j 行的点积是 AA^T 的第 i,j 个条目。同样，第 i 列与第 j 列的点积是 (A^T)A 的第 i,j 个条目。

因此，如果您想要 A 的每个列向量与其自身的点积，您可以使用ColDot = np.dot(np.transpose(A), A).diagonal()。另一方面，如果你想要每一行的点积与它本身，你可以使用RowDot = np.dot(A, np.transpose(A)).diagonal()。

这两行都返回一个数组。

【讨论】：