Numpy cross-product on rectangular grid

Question

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) 

Answer 1

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).