计算多个向量的成对平方差

Question

我正在寻找计算两个向量 ((x1-x2)**2) 之间的平方差的最快方法，但是是成对的（所有组合或仅上三角）。

x1 = [1,3,5,6,8]
x2 = [3,6,7,9,12]

预期输出：

array([[   4.,   25.,   36.,   64.,  121.],
       [   0.,    9.,   16.,   36.,   81.],
       [   4.,    1.,    4.,   16.,   49.],
       [   9.,    0.,    1.,    9.,   36.],
       [  25.,    4.,    1.,    1.,   16.]])

或

array([[   4.,   25.,   36.,   64.,  121.],
       [   0.,    9.,   16.,   36.,   81.],
       [   0.,    0.,    4.,   16.,   49.],
       [   0.,    0.,    0.,    9.,   36.],
       [   0.,    0.,    0.,    0.,    16.]])

甚至（如果更快）：

array([   4.,   25.,   36.,   64.,  121.,    9.,   16.,   36.,   81.,
      4.,    1.,    4.,   16.,   49.,    9.,    1.,    9.,   36.,
     25.,    4.,    1.,    1.,   16.])

Answer 1

您可以使用广播：

x1 = np.asarray([1,3,5,6,8]).reshape(-1, 1)
x2 = np.asarray([3,6,7,9,12]).reshape(1, -1)
(x1 - x2)**2

输出：

array([[  4,  25,  36,  64, 121],
       [  0,   9,  16,  36,  81],
       [  4,   1,   4,  16,  49],
       [  9,   0,   1,   9,  36],
       [ 25,   4,   1,   1,  16]])

编码简单，但计算所有值，因此可以优化为仅计算上三角。

Answer 2

这是一个使用 broadcasting 和 masking 来获得上三角的，然后仅对那些进行平方以获得更好的性能效率的方法 -

def pairwise_squared_diff(x1, x2):
    x1 = np.asarray(x1)
    x2 = np.asarray(x2)
    diffs = x1[:,None] - x2
    mask = np.arange(len(x1))[:,None] <= np.arange(len(x2))
    return (diffs[mask])**2

样本运行-

In [85]: x1
Out[85]: array([1, 3, 5, 6, 8])

In [86]: x2
Out[86]: array([ 3,  6,  7,  9, 12])

In [87]: pairwise_squared_diff(x1, x2)
Out[87]: 
array([  4,  25,  36,  64, 121,   9,  16,  36,  81,   4,  16,  49,   9,
        36,  16])

可能的改进

改进#1：

我们也可以使用np.tri生成mask-

mask = ~np.tri(len(x1),len(x2),dtype=bool,k=-1)

改进#2：

如果我们接受 2D 输出并将下三角输出设置为 0s，那么与 mask 的简单元素乘法也可以解决它以获得最终输出 -

(diffs*mask)**2

这将与 numexpr module 一起很好地处理大数据并提高内存效率和性能。

改进#3：

我们还可以计算与 numexpr 的差异，因此也可以使用相同的 evaulate 方法计算屏蔽输出，从而为我们自己提供一个新的解决方案 -

def pairwise_squared_diff_numexpr(x1, x2):
    x1 = np.asarray(x1)
    x2 = np.asarray(x2)
    mask = ~np.tri(len(x1),len(x2),dtype=bool,k=-1)
    return ne.evaluate('mask*((x1D-x2)**2)',{'x1D':x1[:,None]})

改进时间

让我们研究一下这些关于大型阵列性能的建议 -

设置：

In [136]: x1 = np.random.randint(0,9,(1000))

In [137]: x2 = np.random.randint(0,9,(1000))

改进#1 :

In [138]: %timeit np.arange(len(x1))[:,None] <= np.arange(len(x2))
1000 loops, best of 3: 772 µs per loop

In [139]: %timeit ~np.tri(len(x1),len(x2),dtype=bool,k=-1)
1000 loops, best of 3: 243 µs per loop

改进 #2 :

In [140]: import numexpr as ne

In [141]: diffs = x1[:,None] - x2
     ...: mask = np.arange(len(x1))[:,None] <= np.arange(len(x2))

In [142]: %timeit (diffs[mask])**2
1000 loops, best of 3: 1.46 ms per loop

In [143]: %timeit ne.evaluate('(diffs*mask)**2')
1000 loops, best of 3: 1.05 ms per loop

改进 #3 完整解决方案：

In [170]: %timeit pairwise_squared_diff(x1, x2)
100 loops, best of 3: 3.66 ms per loop

In [171]: %timeit pairwise_squared_diff_numexpr(x1, x2)
1000 loops, best of 3: 1.54 ms per loop

循环一个

为了完整起见，这是一个利用 slicing 比纯 broadcasting 表现更好的循环，这归功于内存效率 -

def pairwise_squared_diff_loopy(x1,x2):
    n = len(x2)
    idx = np.concatenate(( [0], np.arange(n,0,-1).cumsum() ))
    start, stop = idx[:-1], idx[1:]
    L = n*(n+1)//2
    out = np.empty(L,dtype=np.result_type(x1,x2))
    for i,(s0,s1) in enumerate(zip(start,stop)):
        out[s0:s1] = x1[i] - x2[i:]
    return out**2

计时 -

In [300]: x1 = np.random.randint(0,9,(1000))
     ...: x2 = np.random.randint(0,9,(1000))

In [301]: %timeit pairwise_squared_diff(x1, x2)
100 loops, best of 3: 3.44 ms per loop

In [302]: %timeit pairwise_squared_diff_loopy(x1, x2)
100 loops, best of 3: 2.73 ms per loop

计算多个向量的成对平方差

Compute the pairwise squared difference of multiple vectors

python

numpy

pandas

pairwise

可能的改进

循环一个