一组点之间的成对位移向量

Question

我有一个包含 N 个点的数组，这些点在 d 维空间中，形状是 (N, d)。我想为每一对点生成一个新的数组，这些点之间的位移向量数量是 (N 选择 2, d)。如果我只想要这些向量的大小，我可以使用 pdist，这个函数来自 scipy.spatial.distance。

如果我能直接这样做就好了：

pdist(points, lambda u, v: u - v)

但是 metric 函数必须返回一个标量值，所以会出现错误（ValueError: setting an array element with a sequence.）

---

我的解决方案是使用 np.triu_indices：

i, j = np.triu_indices(len(points), 1)
displacements = points[i] - points[j]

这个方法比使用 pdist 慢大约 20 到 30 倍（我通过计算 displacements 的大小来比较，虽然这不是耗时的部分，我猜真正耗时的是生成上三角矩阵和进行复杂的索引操作）。

Answer 1

如果你计算所有差值的完整笛卡尔积，把得到的二维数组压扁，然后自己创建索引来提取上三角部分，你可以让这个过程“仅仅”比pdist慢6倍：

In [39]: points = np.random.rand(1000, 2)

In [40]: %timeit pdist(points)
100 loops, best of 3: 5.81 ms per loop

In [41]: %%timeit
    ...: n = len(points)
    ...: rng = np.arange(1, n)
    ...: idx = np.arange(n *(n-1) // 2) + np.repeat(np.cumsum(rng), rng[::-1])
    ...: np.take((points[:, None] - points).reshape(-1, 2), idx, axis=0)
    ...: 
10 loops, best of 3: 33.9 ms per loop

你还可以通过自己创建索引来加速你的解决方案，并使用取值方法，而不是花哨的索引：

In [75]: %%timeit
    ...: n = len(points)
    ...: rng = np.arange(1, n)
    ...: idx1 = np.repeat(rng - 1, rng[::-1])
    ...: idx2 = np.arange(n*(n-1)//2) + np.repeat(n - np.cumsum(rng[::-1]), rng[::-1])
    ...: np.take(points, idx1, axis=0) - np.take(points, idx2, axis=0)
    ...: 
10 loops, best of 3: 38.8 ms per loop

Answer 2

直接的方法是

dis_vectors = [l - r for l, r in itertools.combinations(points, 2)]

不过我怀疑这个方法快不快。实际上，%timeit显示：

对于3个点：

list : 13 us
pdist: 24 us

但是对于27个点：

list : 798 us
pdist: 35.2 us

我们这里讨论的是多少个点呢？

还有一种可能性是类似于

import numpy
from operator import mul
from fractions import Fraction

def binomial_coefficient(n,k):
    # credit to http://stackoverflow.com/users/226086/nas-banov
    return int( reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1) )

def pairwise_displacements(a):
    n = a.shape[0]
    d = a.shape[1]
    c = binomial_coefficient(n, 2)

    out = numpy.zeros( (c, d) )

    l = 0
    r = l + n - 1
    for sl in range(1, n): # no point1 - point1!
        out[l:r] = a[:n-sl] - a[sl:]
        l = r
        r += n - (sl + 1)
    return out

这个方法简单来说就是让数组在各个维度上“滑动”自己，并在每一步进行可广播的减法。注意，这里没有考虑重复的情况，也就是说不会计算相同的点（比如点1减去点1）。

这个函数在处理1000个点时表现还不错，耗时31.3毫秒，而pdist则更快，耗时20.7毫秒，列表推导法则排在第三，耗时1.23秒。