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

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.])

标签： pythonpandasnumpypairwise

Answer 1

这是一个带有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模块并获得内存效率和性能。

改进#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