CODEWITHSUNDEEP
pythonnumpy
subnivean
subnivean

Reputation: 1152

Numpy cross-product on rectangular grid

I have two numpy arrays holding 2d vectors:

import numpy as np
a = np.array([[ 0.999875,  0.015836],
              [ 0.997443,  0.071463],
              [ 0.686554,  0.727078],
              [ 0.93322 ,  0.359305]])

b = np.array([[ 0.7219  ,  0.691997],
              [ 0.313656,  0.949537],
              [ 0.507926,  0.861401],
              [ 0.818131,  0.575031],
              [ 0.117956,  0.993019]])

As you can see, a.shape is (4,2) and b.shape is (5,2).

Now, I can get the results I want thusly:

In [441]: np.array([np.cross(av, bv) for bv in b for av in a]).reshape(5, 4)
Out[441]: 
array([[ 0.680478,  0.638638, -0.049784,  0.386403],
       [ 0.944451,  0.924694,  0.423856,  0.773429],
       [ 0.85325 ,  0.8229  ,  0.222097,  0.621377],
       [ 0.562003,  0.515094, -0.200055,  0.242672],
       [ 0.991027,  0.982051,  0.595998,  0.884323]])

My question is: What's a more 'numpythonic' way of getting the above (i.e without the nested list comprehensions)? I've tried every combination of np.cross() I can think of, and I usually get results like this:

In [438]: np.cross(a, b.T, axisa=1, axisb=0)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-438-363c0765a7f9> in <module>()
----> 1 np.cross(a, b.T, axisa=1, axisb=0)

D:\users\ae4652t\Python27\lib\site-packages\numpy\core\numeric.p<snipped>
   1242     if a.shape[0] == 2:
   1243         if (b.shape[0] == 2):
-> 1244             cp = a[0]*b[1] - a[1]*b[0]
   1245             if cp.ndim == 0:
   1246                 return cp

ValueError: operands could not be broadcast together with shapes (4) (5)

Upvotes: 9

Views: 1691

Answers (2)

Daniel
Daniel

Reputation: 19547

I thought a bit more on this.

>>> a
array([[ 0.999875,  0.015836],
       [ 0.997443,  0.071463],
       [ 0.686554,  0.727078],
       [ 0.93322 ,  0.359305]])
>>> b
array([[ 0.7219  ,  0.691997],
       [ 0.313656,  0.949537],
       [ 0.507926,  0.861401],
       [ 0.818131,  0.575031],
       [ 0.117956,  0.993019]])
>>> c = np.tile(a, (b.shape[0], 1))
>>> d = np.repeat(b, a.shape[0], axis=0)
>>> np.cross(c, d).reshape(5,4)
array([[ 0.68047849,  0.63863842, -0.0497843 ,  0.38640316],
       [ 0.94445125,  0.92469424,  0.42385605,  0.77342875],
       [ 0.85324981,  0.82290048,  0.22209648,  0.62137629],
       [ 0.5620032 ,  0.51509455, -0.20005522,  0.24267187],
       [ 0.99102692,  0.98205036,  0.59599795,  0.88432301]])

Some timings:

import timeit

s="""
import numpy as np
a=np.random.random(100).reshape(-1, 2)
b=np.random.random(1000).reshape(-1, 2)
"""

ophion="""
np.cross(np.tile(a,(b.shape[0],1)),np.repeat(b,a.shape[0],axis=0))"""

subnivean="""
np.array([np.cross(av, bv) for bv in b for av in a]).reshape(b.shape[0], a.shape[0])"""

DSM="""
np.outer(b[:,1], a[:,0]) - np.outer(b[:,0], a[:,1])"""

Jamie="""
np.cross(a[None], b[:, None, :])"""

h=timeit.timeit(subnivean,setup=s,number=10)
m=timeit.timeit(ophion,setup=s,number=10)
d=timeit.timeit(DSM,setup=s,number=10)
j=timeit.timeit(Jamie,setup=s,number=10)

print "subnivean's method took",h,'seconds.'
print "Ophion's method took",m,'seconds.'
print "DSM's method took",d,'seconds.'

"
subnivean's method took 1.99507117271 seconds.
Ophion's method took 0.0149450302124 seconds.
DSM's method took 0.0040500164032 seconds.
Jamie's method took 0.00390195846558 seconds."

For when the length of a=10 and b=100:

"
subnivean's method took 0.0217308998108 seconds.
Ophion's method took 0.00046181678772 seconds.
DSM's method took 0.000531911849976 seconds.
Jamie's method took 0.000334024429321 seconds."

Hmm you switched the order of the cross product again, both answers are shown if you want (5,4) or (4,5).

Upvotes: 8

Jaime
Jaime

Reputation: 67427

I haven't timed my code, but I am almost certain there is no more numpythonic way of doing this than the nice and simple:

>>> np.cross(a[None], b[:, None])
array([[ 0.68047849,  0.63863842, -0.0497843 ,  0.38640316],
       [ 0.94445125,  0.92469424,  0.42385605,  0.77342875],
       [ 0.85324981,  0.82290048,  0.22209648,  0.62137629],
       [ 0.5620032 ,  0.51509455, -0.20005522,  0.24267187],
       [ 0.99102692,  0.98205036,  0.59599795,  0.88432301]])

Broadcasting is always the answer...

Upvotes: 7

Related Questions