Numpy slow at vector cross product?
BartC
bc at freeuk.com
Mon Nov 21 09:53:35 EST 2016
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
More information about the Python-list
mailing list