<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Sep 11, 2014 at 7:47 PM, Charles R Harris <span dir="ltr"><<a href="mailto:charlesr.harris@gmail.com" target="_blank">charlesr.harris@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote"><div><div class="h5">On Thu, Sep 11, 2014 at 10:10 AM, Charles R Harris <span dir="ltr"><<a href="mailto:charlesr.harris@gmail.com" target="_blank">charlesr.harris@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote"><div><div>On Wed, Sep 10, 2014 at 10:08 PM, Nathaniel Smith <span dir="ltr"><<a href="mailto:njs@pobox.com" target="_blank">njs@pobox.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div>On Wed, Sep 10, 2014 at 4:53 PM, Charles R Harris<br>
<<a href="mailto:charlesr.harris@gmail.com" target="_blank">charlesr.harris@gmail.com</a>> wrote:<br>
><br>
> On Wed, Sep 10, 2014 at 10:52 AM, Pauli Virtanen <<a href="mailto:pav@iki.fi" target="_blank">pav@iki.fi</a>> wrote:<br>
>><br>
>> 09.09.2014, 22:52, Charles R Harris kirjoitti:<br>
>> >   1. Should the operator accept array_like for one of the arguments?<br>
>> >   2. Does it need to handle __numpy_ufunc__, or will<br>
>> >   __array_priority__ serve?<br>
>><br>
>> I think the __matmul__ operator implementation should follow that of<br>
>> __mul__.<br>
>><br>
>> [clip]<br>
>> >    3. Do we want PyArray_Matmul in the numpy API?<br>
>> >    4. Should a matmul function be supplied by the multiarray module?<br>
>> ><br>
>> > If 3 and 4 are wanted, should they use the __numpy_ufunc__ machinery, or<br>
>> > will __array_priority__ serve?<br>
>><br>
>> dot() function deals with __numpy_ufunc__, and the matmul() function<br>
>> should behave similarly.<br>
>><br>
>> It seems dot() uses __array_priority__ for selection of output return<br>
>> subclass, so matmul() probably needs do the same thing.<br>
>><br>
>> > Note that the type number operators, __add__ and such, currently use<br>
>> > __numpy_ufunc__ in combination with __array_priority__, this in addition<br>
>> > to<br>
>> > the fact that they are by default using ufuncs that do the same. I'd<br>
>> > rather<br>
>> > that the __*__ operators simply rely on __array_priority__.<br>
>><br>
>> The whole business of __array_priority__ and __numpy_ufunc__ in the<br>
>> binary ops is solely about when __<op>__ should yield the execution to<br>
>> __r<op>__ of the other object.<br>
>><br>
>> The rule of operation currently is: "__rmul__ before __numpy_ufunc__"<br>
>><br>
>> If you remove the __numpy_ufunc__ handling, it becomes: "__rmul__ before<br>
>> __numpy_ufunc__, except if array_priority happens to be smaller than<br>
>> that of the other class and your class is not an ndarray subclass".<br>
>><br>
>> The following binops also do not IIRC respect __array_priority__ in<br>
>> preferring right-hand operand:<br>
>><br>
>> - in-place operations<br>
>> - comparisons<br>
>><br>
>> One question here is whether it's possible to change the behavior of<br>
>> __array_priority__ here at all, or whether changes are possible only in<br>
>> the context of adding new attributes telling Numpy what to do.<br>
><br>
> I was tempted to make it a generalized ufunc, which would take care of a lot<br>
> of things, but there is a lot of overhead in those functions. Sounds like<br>
> the easiest thing is to make it similar to dot, although having an inplace<br>
> versions complicates the type selection a bit.<br>
<br>
</div></div>Can we please fix the overhead instead of adding more half-complete<br>
implementations of the same concepts? I feel like this usually ends up<br>
slowing things down in the end, as optimization efforts get divided...<br>
<br>
My vote is:<br>
<br>
__matmul__/__rmatmul__ do the standard dispatch stuff that all __op__<br>
methods do (so I guess check __array_priority__ or whatever it is we<br>
always do). I'd also be okay with ignoring __array_priority__ on the<br>
grounds that __numpy_ufunc__ is better, and there's no existing code<br>
relying on __array_priority__ support in __matmul__.<br>
<br>
Having decided that we are actually going to run, they dispatch<br>
unconditionally to np.newdot(a, b) (or newdot(a, b, out=a) for the<br>
in-place version), similarly to how e.g. __add__ dispatches to np.add.<br>
<br>
newdot acts like a standard gufunc with all the standard niceties,<br>
including __numpy_ufunc__ dispatch.<br>
<br>
("newdot" here is intended as a placeholder name, maybe it should be<br>
np.linalg.matmul or something else to be bikeshed later. I also vote<br>
that eventually 'dot' become an alias for this function, but whether<br>
to do that is an orthogonal discussion for later.)<br>
<span><font color="#888888"><br></font></span></blockquote></div></div><div>If we went the ufunc route, I think we would want three of them, matxvec, vecxmat, and matxmat, because the best inner loops would be different in the three cases, but they couldn't be straight ufuncs themselves, as we don't need the other options, `reduce`, etc., but they can't be exactly like the linalg machinery, because we do want subclasses to be able to override. Hmm...<br><br></div><div>The ufunc machinery has some funky aspects. For instance, there are hardwired checks for `__radd__` and other such operators in PyUFunc_GenericFunction  that allows subclasses to overide the ufunc. Those options should really be part of the PyUFuncObject.</div></div></div></div></blockquote><div><br></div></div></div><div>Here are some timing results for the inner product of an integer array comparing inner, ufunc (inner1d), and einsum implementations. The np.long type was chosen because it is implemented for inner1d and to avoid the effect of using cblas.<br><br><span style="font-family:courier new,monospace">In [43]: a = ones(100000, dtype=np.long)<br><br>In [44]: timeit inner(a,a)<br>10000 loops, best of 3: 55.5 µs per loop<br><br>In [45]: timeit inner1d(a,a)<br>10000 loops, best of 3: 56.2 µs per loop<br><br>In [46]: timeit einsum('...i,...i',a,a)<br>10000 loops, best of 3: 43.8 µs per loop<br><br>In [47]: a = ones(1, dtype=np.long)<br><br>In [48]: timeit inner(a,a)<br>1000000 loops, best of 3: 484 ns per loop<br><br>In [49]: timeit inner1d(a,a)<br>1000000 loops, best of 3: 741 ns per loop<br><br>In [50]: timeit einsum('...i,...i',a,a)<br>1000000 loops, best of 3: 811 ns per loop<br></span><br></div><div>For big arrays, einsum has a better inner loop. For small arrays, einsum and inner1d suffer from call overhead. Note that the einsum overhead could be improved by special casing, and that the various loops could be consolidated into best of breed. <br><br></div><div>The ufunc implementation is easy to do for all cases.<br></div><div><br></div></div></div></div></blockquote><div><br></div><div>Same thing for doubles. The speedup due to cblas can be seen, and the iterator overhead becomes more visible.<br><br><span style="font-family:courier new,monospace">In [51]: a = ones(100000, dtype=np.double)<br><br>In [52]: timeit inner(a,a)<br>10000 loops, best of 3: 32.1 µs per loop<br><br>In [53]: timeit inner1d(a,a)<br>10000 loops, best of 3: 134 µs per loop<br><br>In [54]: timeit einsum('...i,...i',a,a)<br>10000 loops, best of 3: 42 µs per loop<br><br>In [55]: a = ones(1, dtype=np.double)<br><br>In [56]: timeit inner(a,a)<br>1000000 loops, best of 3: 336 ns per loop<br><br>In [57]: timeit inner1d(a,a)<br>1000000 loops, best of 3: 742 ns per loop<br><br>In [58]: timeit einsum('...i,...i',a,a)<br>1000000 loops, best of 3: 809 ns per loop</span><br><br></div><div>But stacked vectors make a difference<br><br><span style="font-family:courier new,monospace">In [59]: a = ones((2, 1), dtype=np.double)<br><br>In [60]: timeit for i in range(2): inner(a[i],a[i])<br>1000000 loops, best of 3: 1.39 µs per loop<br><br>In [61]: timeit inner1d(a,a)<br>1000000 loops, best of 3: 749 ns per loop<br><br>In [62]: timeit einsum('...i,...i',a,a)<br>1000000 loops, best of 3: 888 ns per loop<br></span><br></div><div>Chuck<br></div></div></div></div>