如何矢量化轨迹段之间的点积

Question

这是轨迹（xy 坐标）的连续段之间的点积函数。结果符合预期，但 "for loop" 使其非常慢。

In [94]:
def func1(xy, s):
    size = xy.shape[0]-2*s
    out = np.zeros(size)
    for i in range(size):
        p1, p2 = xy[i], xy[i+s]     #segment 1
        p3, p4 = xy[i+s], xy[i+2*s] #segment 2
        out[i] = np.dot(p1-p2, p4-p3)
    return out

xy = np.array([[1,2],[2,3],[3,4],[5,6],[7,8],[2,4],[5,2],[9,9],[1,1]])
func1(xy, 2)

Out[94]:
array([-16.,  15.,  32.,  31., -14.])

我正在寻找一种对上述内容进行矢量化处理的方法，希望能使其更快。这是我想出的：

In [95]:
def func2(xy, s):
    size = xy.shape[0]-2*s
    p1 = xy[0:size]
    p2 = xy[s:size+s]
    p3 = p2
    p4 = xy[2*s:size+2*s]
    return np.diagonal(np.dot((p1-p2), (p4-p3).T))

func2(xy, 2)

Out[95]:
array([-16,  15,  32,  31, -14])

不幸的是，点积产生了一个方阵，我必须从中取对角线：

In [96]:
print np.dot((p1-p2), (p4-p3).T)
np.diagonal(np.dot((p1-p2), (p4-p3).T))

[[-16  10  16 -24  10]
 [-24  15  24 -36  15]
 [-32  20  32 -48  20]
 [ 20 -13 -18  31 -14]
 [ 32 -18 -40  44 -14]]

Out[96]:
array([-16,  15,  32,  31, -14])

我的解决方案真的很糟糕。它仅将速度提高了 2 倍，更重要的是，它现在不可扩展。我的平均轨迹有几万个点，这意味着我将不得不处理巨大的矩阵。

你们知道更好的方法吗？谢谢

编辑：惊人的！ einsum 绝对是解决方案。在我沮丧的时候，我自己写了点积。我知道，可读性不是很好，它违背了使用优化库的目的，但无论如何 (func4)。速度与 einsum 相当。

def func4(xy, s):
    size = xy.shape[0]-2*s
    tmp1 = xy[0:size] - xy[s:size+s]
    tmp2 = xy[2*s:size+2*s] - xy[s:size+s]
    return tmp1[:, 0] * tmp2[:, 0] + tmp1[:, 1] * tmp2[:, 1]

Answer 1

您在 func2 中的想法自然会导致使用 np.einsum。

func2 的优点是它只计算 p1、p2、p3、p4 一次更大的阵列而不是 func1.

中的小块

func2 的缺点是它做了很多你不会做的点积关心.

这就是 einsum 的用武之地。它是 np.dot 的更灵活的版本。每当计算乘积总和时，请考虑使用 np.einsum。它可能是最快的（如果不是最快的）计算使用 NumPy 的数量。

def func3(xy, s):
    size = xy.shape[0]-2*s
    p1 = xy[0:size]
    p2 = xy[s:size+s]
    p3 = p2
    p4 = xy[2*s:size+2*s]
    return np.einsum('ij,ij->i', p1-p2, p4-p3)

下标串'ij,ij->i'含义如下：

下标字符串有两部分'ij,ij->i'：在箭头之前 (->)，在左侧，ij,ij 和箭头后，i。

左边逗号前的ij指的是p1-p2的下标，逗号后的ij指的是p4-p3.

的下标

爱因斯坦求和符号对未出现在箭头后的重复下标求和。在这种情况下，j 重复出现并且不会出现在箭头之后。

因此，对于每个 i，计算总和 (p1-p2)[i,j]*(p4-p3)[i,j]，总和遍及所有 j。结果是一个由 i.

索引的数组

完整性检查：

In [90]: np.allclose(func1(xy, 2), func3(xy, 2))
Out[90]: True

这是一个基准：在形状为 (9000, 2) 的数组 xy 上，显示使用 einsum 比 func1 快 450 倍，比 func2:

快 7470 倍

In [13]: xy = np.tile(xy, (1000,1))

In [14]: %timeit func1(xy, 2)
10 loops, best of 3: 42.1 ms per loop

In [15]: %timeit func2(xy, 2)
1 loops, best of 3: 686 ms per loop

In [16]: %timeit func3(xy, 2)
10000 loops, best of 3: 91.8 µs per loop

OP 的 func4 比 func3 表现得更好！

In [92]: %timeit func4(xy, 2)
10000 loops, best of 3: 74.1 µs per loop

我认为 func4 在这里击败 einsum 的原因是因为与手动编写相比，在 einsum 中设置循环仅两次迭代的成本太高求和。

Answer 2

einsum 是推广 dot 产品的好工具。玩弄它，我可以重现你的数字：

np.einsum('ij,ij->i',p1-p2,p4-p3)

'ij,kj' 产生 dot(p1-p2, (p4-p3).T)； 'i...,i...->i' 做对角线 - 一步完成。

作为我试过的交叉产品问题的副产品

tmp11,tmp21),tmp11[:,0]*tmp21[:,0]+tmp11[:,1]*tmp21[:,1])

对于 5000 行数组，它几乎是 einsum 计算速度的 2 倍。

如何矢量化轨迹段之间的点积

How to vectorize the dot product between segments of a trajectory

python

numpy