Can numpy einsum () perform cross-product between path segments

Question

I am doing a cross-product of contiguous path segments (xy coordinates) using the following script:

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 is especially slow due to the inner loop, so I rewrote the cross product myself (func2), which is orders of magnitude faster.

Is it possible to use the numpy einsum function to perform the same calculation?

Answer 1

einsum only calculates the sum of the pieces, but you can retell the cross-product into the sum of the pieces by changing the tmp2 columns and changing the sign of the first column:

 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) 

---

But func3 slower than 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 

Health Check:

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

I think the reason func2 beat einsum here is because the cost of setting up a loop in einsum just 2 iterations is too expensive compared to just manually writing out the amount, and reversing and multiplying eats up some time.

Answer 2

np.cross is a smart little beast that can handle broadcasts without any problems. Therefore, you can rewrite your func2 as:

 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) 

and it will give the correct result:

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

In the latest numpy, it will most likely run a little faster than your code, as it has been overwritten to minimize the creation of an intermediate array. You can see the source code (pure Python) here .