[Numpy-discussion] @ operator

[Nathaniel Smith]

On Thu, Sep 11, 2014 at 10:12 PM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
>
> On Thu, Sep 11, 2014 at 8:01 PM, Nathaniel Smith <njs at pobox.com> wrote:
>>
>> On Thu, Sep 11, 2014 at 12:10 PM, Charles R Harris
>> <charlesr.harris at gmail.com> wrote:
>> >
>> > On Wed, Sep 10, 2014 at 10:08 PM, Nathaniel Smith <njs at pobox.com> wrote:
>> >>
>> >> My vote is:
>> >>
>> >> __matmul__/__rmatmul__ do the standard dispatch stuff that all __op__
>> >> methods do (so I guess check __array_priority__ or whatever it is we
>> >> always do). I'd also be okay with ignoring __array_priority__ on the
>> >> grounds that __numpy_ufunc__ is better, and there's no existing code
>> >> relying on __array_priority__ support in __matmul__.
>> >>
>> >> Having decided that we are actually going to run, they dispatch
>> >> unconditionally to np.newdot(a, b) (or newdot(a, b, out=a) for the
>> >> in-place version), similarly to how e.g. __add__ dispatches to np.add.
>> >>
>> >> newdot acts like a standard gufunc with all the standard niceties,
>> >> including __numpy_ufunc__ dispatch.
>> >>
>> >> ("newdot" here is intended as a placeholder name, maybe it should be
>> >> np.linalg.matmul or something else to be bikeshed later. I also vote
>> >> that eventually 'dot' become an alias for this function, but whether
>> >> to do that is an orthogonal discussion for later.)
>> >>
>> > 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,
>>
>> Couldn't we write a single inner loop like:
>>
>> void ufunc_loop(blah blah) {
>>     if (arg1_shape[0] == 1 && arg2_shape[1] == 1) {
>>         call DOT
>>     } else if (arg2_shape[0] == 1) {
>>         call GEMV
>>     } else if (...) {
>>         ...
>>     } else {
>>         call GEMM
>>     }
>> }
>> ?
>
> Not for generalized ufuncs, different signatures, or if linearized, more
> info on dimensions.

This sentence no verb, but I think the point you might be raising is:
we don't actually have the technical capability to define a single
gufunc for @, because the matmat, matvec, vecmat, and vecvec forms
have different gufunc signatures ("mn,nk->mk", "mn,n->m", "n,nk->k",
and "n,n->" respectively, I think)?

This is true, but Jaime has said he's willing to look at fixing this :-):
    http://thread.gmane.org/gmane.comp.python.numeric.general/58669/focus=58670

...and fundamentally, it's very difficult to solve this anywhere else
except the ufunc internals. When we first enter a function like
newdot, we need to check for special overloads like __numpy_ufunc__
*before* we do any array_like->ndarray coercion. But we have to do
array_like->ndarray coercion before we know what the shape/ndim of the
our inputs is. And in the wrapper-around-multiple-ufuncs approach, we
have to check the shape/ndim before we choose which ufunc to dispatch
to. So __numpy_ufunc__ won't get checked until after we've already
coerced to ndarray, oops... OTOH the ufunc internals already know how
to do this dance, so it's at least straightforward (if not necessarily
trivial) to insert the correct logic once in the correct place.

> What you show is essentially what dot does now for cblas
> enabled functions. But note, we need more than the simple '@', we also need
> stacks of vectors, and turning vectors into matrices, and then back into
> vectors seems unnecessarily complicated.

By the time we get into the data/shape/strides world that ufuncs work
with, converting vectors->matrices is literally just adding a single
entry to the shape+strides arrays. This feels trivial and cheap to me?

Or do you just mean that we do actually want broadcasting matvec,
vecmat, vecvec gufuncs? I agree with this too but it seems orthogonal
to me -- we can have those in any case, even if newdot doesn't
literally call them.

-n

-- 
Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org