0

我正在寻找最快的方法来计算两个向量((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.])
4

2 回答 2

4

这是一个带有broadcastingmasking获得上三角形的,然后只对那些进行平方以获得更好的性能效率 -

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
于 2018-10-24T09:37:47.677 回答
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]])

这很容易编码,但会计算所有值,因此可以对其进行优化以仅计算上三角形。

于 2018-10-24T09:33:31.733 回答