这是轨迹的连续段(xy坐标)之间的点积的函数.结果如预期,但"for循环"使它非常慢.
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]
你的想法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
这是一个基准:在一个xy形状数组(einsum
9000,2 )上,显示使用比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与仅手动写出总和相比,仅仅2次迭代设置循环的成本太高了.