Calculating "element-wise" the angles between two lists of vectors

Question

Let's say you have two lists of vectors:

v1s = [a, b, c]

v2s = [d, e, f]

I am interested in generating the following result:

angles = [angleBetween(a, d), angleBetween(b, e), angleBetween(c, f)]

Here is my function to do this (it uses vector dot products in order to calculate the angles):

import numpy as np

def findAnglesBetweenTwoVectors(v1s, v2s):
    dot_v1_v2 = np.einsum('ij,ij->i', v1s, v2s)
    dot_v1_v1 = np.einsum('ij,ij->i', v1s, v1s)
    dot_v2_v2 = np.einsum('ij,ij->i', v2s, v2s)

    return np.arccos(dot_v1_v2/(np.sqrt(dot_v1_v1)*np.sqrt(dot_v2_v2)))

I call this function a lot, so optimizing it would be very helpful for me. Is this possible?

Answer 1

Here's a way of packing the calculation into on call to einsum.

def findAnglesBetweenTwoVectors1(v1s, v2s):
    dot = np.einsum('ijk,ijk->ij',[v1s,v1s,v2s],[v2s,v1s,v2s])
    return np.arccos(dot[0,:]/(np.sqrt(dot[1,:])*np.sqrt(dot[2,:])))

for random vectors of length 10, the speed for this function and yours is basically the same. But with 10000, (corrected times)

In [62]: timeit findAnglesBetweenTwoVectors0(v1s,v2s)
100 loops, best of 3: 2.26 ms per loop

In [63]: timeit findAnglesBetweenTwoVectors1(v1s,v2s)
100 loops, best of 3: 4.24 ms per loop

Yours is clearly faster.

We could try to figure out why the more complex einsum is slower (assuming the issue is there rather than the indexing in the last line). But this result is consistent with other einsum testing that I've seen and done. einsum on simple 2D cases is as fast as np.dot, and in some cases faster. But in more complex cases (3 or more indexes) it is slower than the equivalent constructed from multiple calls to np.dot or einsum.

---

I tried tried factoring out the array construction, and got a speed improvement.

In [12]: V1=np.array(v1s)  # shape (3,10000,3)
In [13]: V2=np.array(v2s)
In [22]: timeit findAnglesBetweenTwoVectors2(V1,V2)
100 loops, best of 3: 2.63 ms per loop

---

I've corrected the times, using v1s.shape==(10000,3)