成对平方差异的有效Numpy计算

Question

下面的代码完全符合我的要求,即计算向量元素之间差异的平方和(本例中长度为3),其中我有一个长序列(此处限制为5).所需的结果显示在底部.但由于两个原因,实施感觉很糟糕:

1)需要添加一个幻像维度,将形状从(5,3)改为(5,1,3)以避免广播问题,以及

2)显然需要一个明确的'for'循环,我敢肯定,为什么需要花费数小时来执行我更大的数据集(一百万个长度为2904的向量).

是否有更高效和/或pythonic的方法来实现相同的结果？

a = np.array([[ 4,  2,  3], [-1, -5,  4], [ 2,  1,  4], [-5, -1,  4], [6, -3,  3]])
a = a.reshape((5,1,3))

m = a.shape[0]
n = a.shape[2]
d = np.zeros((n,n))
for i in range(m):
    c = a[i,:] - np.transpose(a[i,:])
    c = c**2
    d += c

print d

[[   0.  118.  120.]
 [ 118.    0.  152.]
 [ 120.  152.    0.]]

Answer 1

您可以使用以下方法消除 for 循环：

In [48]: ((a - a.swapaxes(1,2))**2).sum(axis=0)
Out[48]: 
array([[  0, 118, 120],
       [118,   0, 152],
       [120, 152,   0]])

请注意，如果a具有 shape(N, 1, M)则(a - a.swapaxes(1,2))具有 shape (N, M, M)。确保您有足够的 RAM 来容纳这种大小的数组。页面交换也可以将计算减慢到爬行。

如果你的内存太少，你将不得不将计算分成块：

m, _, n = a.shape
chunksize = 10**4
d = np.zeros((n,n))
for i in range(0, m, chunksize):
    b = a[i:i+chunksize]
    d += ((b - b.swapaxes(1,2))**2).sum(axis=0)

这是对整个数组执行计算和逐行计算之间的折衷。如果有一百万行，并且块大小为 10**4，那么循环将只有 100 次迭代而不是一百万次。因此，它应该比逐行计算快得多。选择最大的 chunksize 值，它允许在 RAM 中执行计算。

Answer 2

如果你不介意依赖scipy,你可以使用scipy.spatial.distance库中的函数:

In [17]: from scipy.spatial.distance import pdist, squareform

In [18]: a = np.array([[ 4,  2,  3], [-1, -5,  4], [ 2,  1,  4], [-5, -1,  4], [6, -3,  3]])

In [19]: d = pdist(a.T, metric='sqeuclidean')

In [20]: d
Out[20]: array([ 118.,  120.,  152.])

In [21]: squareform(d)
Out[21]: 
array([[   0.,  118.,  120.],
       [ 118.,    0.,  152.],
       [ 120.,  152.,    0.]])