numpy的einsum()函数能否在轨迹的不同段之间执行叉积运算？

Question

我使用以下脚本执行轨迹（xy坐标）的连续段的叉积：

In [129]:
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.cross(p1-p2, p4-p3)
    return out

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]

    tmp1 = p1-p2
    tmp2 = p4-p3
    return tmp1[:, 0] * tmp2[:, 1] - tmp2[:, 0] * tmp1[:, 1]


In [136]:
xy = np.array([[1,2],[2,3],[3,4],[5,6],[7,8],[2,4],[5,2],[9,9],[1,1]])
func2(xy, 2)

Out[136]:
array([ 0, -3, 16,  1, 22])

func1由于内部循环而特别缓慢，因此我自己重新编写了叉积函数（func2），这个函数的速度比原来快了好几个数量级。

是否可能使用numpy的einsum函数进行相同的计算？

- user3329302

---

在我的测试中，您的 func2 比其他替代方案更快，甚至比新的 cross 更快。 - hpaulj

Answer 1

einsum函数只能计算积的和，但是可以通过将tmp2的列翻转并更改第一列的符号，将叉积强行转化为积的和:

def func3(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]
    tmp2 = tmp2[:, ::-1]
    tmp2[:, 0] *= -1
    return np.einsum('ij,ij->i', tmp1, tmp2)

---

但是func3比func2慢。

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

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

In [105]: %timeit func2(xy, 2)
10000 loops, best of 3: 73.2 µs per loop

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

合理性检查：

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

我认为func2之所以能在这里胜过einsum的原因是，仅仅为了进行两次迭代而在einsum中设置循环的成本太高了，相比手动写出求和，反转和乘法也需要一些时间。

- unutbu

---

谢谢unutbu，我知道使用基于common view的方法翻转矩阵列也相当慢。 - user3329302

Answer 2

np.cross 是一个非常智能的小家伙，可以处理广播而不会出现任何问题。因此，您可以将您的 func2 重写为：

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.cross(p1-p2, p4-p3)

并且它将会产生正确的结果：

>>> func2(xy, 2)
array([ 0, -3, 16,  1, 22])

在最新的numpy中，它可能比你的代码运行速度稍快，因为它被重写以最小化中间数组的创建。您可以查看源代码（纯Python）此处。

- Jaime

---

你所提到的 cross 版本在2014年2月和6月有重大更改，版本号为v1.9.1。将感兴趣的轴滚动到末尾，并在可能的情况下使用 multiply(...,out=c)。 - hpaulj