Numpy slow at vector cross product?

[Peter Otten]

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