[Numpy-discussion] Preliminary thoughts on implementing __matmul__

[Jaime Fernández del Río]

On Wed, Aug 6, 2014 at 5:31 AM, Nathaniel Smith <njs at pobox.com> wrote:

> I think the other obvious strategy to consider, is defining a 'dot'
> gufunc, with semantics identical to @. (This would be useful for
> backcompat as well: adding/dropping compatibility with older python
> versions would be as simple as mechanically replacing a @ b with
> newdot(a, b) or vice-versa.) This would require one new feature in the
> gufunc machinery: support for "optional core axes", to get the right
> semantics for 1d arrays.


Can you elaborate on what those optional core axes would look like? If I am
understanding you correctly, this is what now is solved by having more than
one gufunc defined, and choosing which one to use based on the input's
shapes in a thin Python wrapper. There are several examples in the linalg
module you are certainly well aware of.

Say you could define the matmul signature as "(i)j,j(k)->(ik)", with
dimensions in parenthesis being "optional." Say we modified the gufunc
machinery to detect which optional core axes are present and which not. It
seems to me that you would then still need to write 4 traditional gufuncs
(ij,jk->ik, j,jk->k, ij,j->i, j,j->) and dispatch to one of them. I haven't
thought it through, but are there really a set of universal dispatch rules
that will apply to any optional core axes problem? Would we not be losing
flexibility in doing so?

When I looked into gufuncs several months ago, what I missed was a way of
defining signatures like n,m->n*(n-1), which would come in very handy if
computing all pairwise distances. You can work around this by making the
signature n,m->p and always calling the gufunc from a Python wrapper that
passes in an out parameter of the right shape. But if someone gets a hold
of the gufunc handle and calls it directly without an out parameter, the p
defaults to 1 and you are probably in for a big crash. So it would be nice
if you could provide a pointer to a function to produce the output shape
based on the inputs'.

On my wish list for gufunc signatures there is also frozen dimensions, e.g.
a gufunc to compute greater circle distances on a sphere can be defined as
m,m->, but m has to be 2, and since you don't typically want to be raising
errors in the kernel, a Python wrapper is once more necessary. And again an
unwrapped call to the gufunc is potentially catastrophic.

Sorry for hijacking the thread, but I wouldn't mind spending some time
working on expanding this functionality to include the optional axes and my
wish-list, if the whole thing makes sense.

Jaime
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20140806/4dd952d3/attachment.html>