Numpy slow at vector cross product?
Peter Otten
__peter__ at web.de
Mon Nov 21 07:44:25 EST 2016
Steve D'Aprano wrote:
> On Mon, 21 Nov 2016 07:46 am, DFS wrote:
>
>> import sys, time, numpy as np
>> loops=int(sys.argv[1])
>>
>> x=np.array([1,2,3])
>> y=np.array([4,5,6])
>> start=time.clock()
>> for i in range(loops):
>> np.cross(x,y)
>> print "Numpy, %s loops: %.2g seconds" %(loops,time.clock()-start)
>
> [...]
>> $ python vector_cross.py
>> Numpy, 100000 loops: 2.5 seconds
>> Calc, 100000 loops: 0.13 seconds
>>
>>
>> Did I do something wrong, or is numpy slow at this?
>
> I can confirm similar results.
>
> However, your code is not a great way of timing code. Timing code is
> *very* difficult, and can be effected by many things, such as external
> processes, CPU caches, even the function you use for getting the time.
> Much of the time you are timing here will be in creating the range(loops)
> list, especially if loops is big.
>
> The best way to time small snippets of code is to use the timeit module.
> Open a terminal or shell (*not* the Python interactive interpreter, the
> operating system's shell: you should expect a $ or % prompt) and run
> timeit from that. Copy and paste the following two commands into your
> shell prompt:
>
>
> python2.7 -m timeit --repeat 5 -s "import numpy as np" \
> -s "x = np.array([1, 2, 3])" -s "y = np.array([4, 5, 6])" \
> -- "np.cross(x, y)"
>
>
> python2.7 -m timeit --repeat 5 -s "x = [1, 2, 3]" \
> -s "y = [4, 5, 6]" -s "z = [0, 0, 0]" \
> -- "z[0] = x[1]*y[2] - x[2]*y[1]; z[1] = x[2]*y[0] - \
> x[0]*y[2]; z[2] = x[0]*y[1] - x[1]*y[0]"
>
>
> The results I get are:
>
> 10000 loops, best of 5: 30 usec per loop
>
> 1000000 loops, best of 5: 1.23 usec per loop
>
>
> So on my machine, np.cross() is about 25 times slower than multiplying by
> hand.
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)
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 = array([x, y, z])
if cp.ndim == 1:
return cp
else:
return cp.swapaxes(0, axisc)
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).
More information about the Python-list
mailing list