CODEWITHSUNDEEP
pythonnumpy
Ferdinando Randisi
Ferdinando Randisi

Reputation: 4388

Efficient way of computing the cross products between two sets of vectors numpy

I have two sets of 2000 3D vectors each, and I need to compute the cross product between each possible pair. I currently do it like this

for tx in tangents_x:
    for ty in tangents_y:
         cross = np.cross(tx, ty)
         (... do something with the cross variable...)

This works, but it's pretty slow. Is there a way to make it faster?

If I was interested in the element-wise product, I could just do the following

# Define initial vectors
tx = np.array([np.random.randn(3) for i in range(2000)])
ty = np.array([np.random.randn(3) for i in range(2000)])
# Store them into matrices
X = np.array([tx for i in range(2000)])
Y = np.array([ty for i in range(2000)]).T
# Compute the element-wise product
ew = X * Y
# Use the element_wise product as usual
for i,tx in enumerate(tangents_x):
    for j,ty in enumerate(tangents_y):
        (... use the element wise product of tx and ty as ew[i,j])

How can I apply this to the cross product instead of the element-wise one? Or, do you see another alternative?

Thanks much :)

Upvotes: 4

Views: 4861

Answers (5)

Brenlla
Brenlla

Reputation: 1481

np.dot is almost always going to be faster. So you could convert one of the vectors into a matrix.

def skew(x):
    return np.array([[0, -x[2], x[1]],
                     [x[2], 0, -x[0]],
                     [-x[1], x[0], 0]])

On my machine this runs faster:

tx = np.array([np.random.randn(3) for i in range(100)])
ty = np.array([np.random.randn(3) for i in range(100)])

tt=time.clock()
for x in tx:
    for y in ty:
         cross = np.cross(x, y)
print(time.clock()-tt)

0.207 sec

tt=time.clock()
for x in tx:
    m=skew(x)
    for y in ty:
         cross = np.dot(m, y)
print(time.clock()-tt)

0.015 sec

This result may vary depending on the computer.

Upvotes: 2

Paul Panzer
Paul Panzer

Reputation: 53029

Like many numpy functions cross supports broadcasting, therefore you can simply do:

np.cross(tangents_x[:, None, :], tangents_y)

or - more verbose but maybe easier to read

np.cross(tangents_x[:, None, :], tangents_y[None, :, :])

This reshapes tangents_x and tangents_y to shapes 2000, 1, 3 and 1, 2000, 3. By the rules of broadcasting this will be interpreted like two arrays of shape 2000, 2000, 3 where tangents_x is repeated along axis 1 and tangents_y is repeated along axis 0.

Upvotes: 6

max9111
max9111

Reputation: 6482

Just write it out and compile it

import numpy as np
import numba as nb

@nb.njit(fastmath=True,parallel=True)
def calc_cros(vec_1,vec_2):
    res=np.empty((vec_1.shape[0],vec_2.shape[0],3),dtype=vec_1.dtype)
    for i in nb.prange(vec_1.shape[0]):
        for j in range(vec_2.shape[0]):
            res[i,j,0]=vec_1[i,1] * vec_2[j,2] - vec_1[i,2] * vec_2[j,1]
            res[i,j,1]=vec_1[i,2] * vec_2[j,0] - vec_1[i,0] * vec_2[j,2]
            res[i,j,2]=vec_1[i,0] * vec_2[j,1] - vec_1[i,1] * vec_2[j,0]
    
    return res

Performance

#create data
tx = np.random.rand(3000,3)
ty = np.random.rand(3000,3)
#don't measure compilation overhead
comb=calc_cros(tx,ty)

t1=time.time()
comb=calc_cros(tx,ty)
print(time.time()-t1)

This gives 0.08s for the two (3000,3) matrices.

Upvotes: 4

Alfe
Alfe

Reputation: 59426

You could use np.meshgrid() to build the combination matrix and then decompose the cross product. The rest is fiddling around with the axes etc:

# build two lists of 5 3D vecotrs as example values:
a_list = np.random.randint(0, 10, (5, 3))
b_list = np.random.randint(0, 10, (5, 3))

# here the original approach using slow list comprehensions:
slow = np.array([[ np.cross(a, b) for a in a_list ] for b in b_list ])

# now the faster proposed version:
g = np.array([ np.meshgrid(a_list[:,i], b_list[:,i]) for i in range(3) ])
fast = np.array([ g[1,0] * g[2,1] - g[2,0] * g[1,1],
                  g[2,0] * g[0,1] - g[0,0] * g[2,1],
                  g[0,0] * g[1,1] - g[1,0] * g[0,1] ]).transpose(1, 2, 0)

I tested this with 10000×10000 elements (instead of the 5×5 in the example above) and it took 6.4 seconds with the fast version. The slow version already took 27 seconds for 500 elements.

For your 2000×2000 elements the fast version takes 0.23s on my computer. Fast enough for you?

Upvotes: 1

lr1985
lr1985

Reputation: 1232

Use a cartesian product to get all possible pairs 

import itertools as it
all_pairs = it.product(tx, ty)

And then use map to loop over all pairs and compute the cross product:

map(lambda x: np.cross(x[0], x[1]), all_pairs)

Upvotes: 0

Related Questions