Numpy slow at vector cross product?

[BartC]

On 21/11/2016 12:44, Peter Otten wrote:

> After a look into the source this is no longer a big surprise (numpy 1.8.2):
>
>     if axis is not None:
>         axisa, axisb, axisc=(axis,)*3
>     a = asarray(a).swapaxes(axisa, 0)
>     b = asarray(b).swapaxes(axisb, 0)

>
> The situation may be different when you process vectors in bulk, i. e.
> instead of [cross(a, bb) for bb in b] just say cross(a, b).

I thought that numpy was supposed to be fast? Also that the critical 
bits were not implemented in Python?

Anyway I tried your code (put into a function as shown below) in the 
same test program, and it was *still* 3 times as fast as numpy! 
(mycross() was 3 times as fast as np.cross().)

Explain that...


---------------------------------------

def mycross(a,b,axisa=-1, axisb=-1, axisc=-1, axis=None):
     if axis is not None:
         axisa, axisb, axisc=(axis,)*3
     a = np.asarray(a).swapaxes(axisa, 0)
     b = np.asarray(b).swapaxes(axisb, 0)
     msg = "incompatible dimensions for cross product\n"\
           "(dimension must be 2 or 3)"
     if (a.shape[0] not in [2, 3]) or (b.shape[0] not in [2, 3]):
         raise ValueError(msg)
     if a.shape[0] == 2:
         if (b.shape[0] == 2):
             cp = a[0]*b[1] - a[1]*b[0]
             if cp.ndim == 0:
                 return cp
             else:
                 return cp.swapaxes(0, axisc)
         else:
             x = a[1]*b[2]
             y = -a[0]*b[2]
             z = a[0]*b[1] - a[1]*b[0]
     elif a.shape[0] == 3:
         if (b.shape[0] == 3):
             x = a[1]*b[2] - a[2]*b[1]
             y = a[2]*b[0] - a[0]*b[2]
             z = a[0]*b[1] - a[1]*b[0]
         else:
             x = -a[2]*b[1]
             y = a[2]*b[0]
             z = a[0]*b[1] - a[1]*b[0]
     cp = np.array([x, y, z])
     if cp.ndim == 1:
         return cp
     else:
         return cp.swapaxes(0, axisc)

---------------------------------------
Tested as:

	x=np.array([1,2,3])
	y=np.array([4,5,6])

	start=time.clock()
	for i in range(loops):
	    z=mycross(x,y)
	print "Calc, %s loops: %.2g seconds" %(loops,time.clock()-start)


-- 
Bartc