Numpy slow at vector cross product?

[BartC]

On 22/11/2016 16:48, Steve D'Aprano wrote:
> On Tue, 22 Nov 2016 11:45 pm, BartC wrote:
>
>
>> I will have a look. Don't forget however that all someone is trying to
>> do is to multiply two vectors. They're not interested in axes
>> transformation or making them broadcastable, whatever that means.
>
> You don't know that.

When I did vector arithmetic at school, the formula for the 
cross-product of (x1,y1,z1) with (x2,y2,z2) was quite straightforward. 
This is what you expect given a library function to do the job. 
Especially in an implementation that could be crippled by running as 
unoptimised byte-code.

If someone was given the task of cross-multiplying two three-element 
vectors, the formula in the text-book (or Wikipedia) is what they would 
write.

> Could numpy optimize the single pair of vectors case a bit better? Perhaps
> they could. It looks like there's a bunch of minor improvements which could
> be made to the code, a few micro-optimizations that shave a microsecond or
> two off the execution time. And maybe they could even detect the case where
> the arguments are a single pair of vectors, and optimize that. Even
> replacing it with a naive pure-Python cross product would be a big win.

The code it ends up executing for individual pairs of vectors is a joke. 
97% of what it does is nothing to do with the cross-product!

> But for the big array of vectors case, you absolutely have to support doing
> fast vectorized cross-products over a huge number of vectors.

That's how I expected numpy to work. I said something along those lines 
in my first post:

BC:
 > Maybe numpy has extra overheads, and the arrays being operated on are
 > very small, but even so, 30 times slower than CPython? (2.5 to 0.083
 > seconds.)


> py> a = np.array([
> ...     [1, 2, 3],
> ...     [4, 5, 6],
> ...     [7, 8, 9],
> ...     ]*1000
> ...     )
> py> b = np.array([
> ...     [9, 8, 7],
> ...     [6, 5, 4],
> ...     [3, 2, 1],
> ...     ]*1000
> ...     )
> py> a.shape
> (3000, 3)
> py> result = np.cross(a, b)
> py> result.shape
> (3000, 3)
>
> On my computer, numpy took only 10 times longer to cross-multiply 3000 pairs
> of vectors than it took to cross-multiply a single pair of vectors. If I
> did that in pure Python, it would take 3000 times longer, or more, so numpy
> wins here by a factor of 300.

I tested this with 3 million pairs on my machine. It managed 8 million 
products per second (about the same as the simplest Python code, but run 
with pypy), taking nearly 0.4 seconds. Except it took over 4 seconds to 
create the special numpy arrays!

Setting up ordinary arrays was much faster, but then the cross-product 
calculation was slower.

-- 
bartc