[Numpy-discussion] Allowing broadcasting of code dimensions in generalized ufuncs

Marten van Kerkwijk m.h.vankerkwijk at gmail.com
Mon Jun 11 10:59:49 EDT 2018

> Nathaniel:
> Output shape feels very similar to
> output dtype to me, so maybe the general way to handle this would be
> to make the first callback take the input shapes+dtypes and return the
> desired output shapes+dtypes?
> This hits on an interesting alternative to frozen dimensions - np.cross
> could just become a regular ufunc with signature np.dtype((float64, 3)),
> np.dtype((float64, 3)) → np.dtype((float64, 3))
As you note further down, the present proposal of just using numbers has
the advantage of being clear and easy. Another (small?) advantage is that I
can use `axis` to tell where my three coordinates are, rather than be stuck
with having them as the last dimension.

Indeed, in my trials for wrapping the Standards Of Fundamental Astronomy
routines, I started with just making every 3-vector and 3x3-matrix
structured arrays with the relevant single sub-array entry. That worked,
but I ended up disliking the casting to and fro.

> Furthermore, the expansion quickly becomes cumbersome. For instance, for
> the all_equal signature of (n|1),(n|1)->() …
> I think this is only a good argument when used in conjunction with the
> broadcasting syntax. I don’t think it’s a reason for matmul not to have
> multiple signatures. Having multiple signatures is an disincentive to
> introduced too many overloads of the same function, which seems like a good
> thing to me
But implementation for matmul is actually considerably trickier, since the
internal loop now has to check the number of distinct dimensions.

> Summarizing my overall opinions:
>    - I’m +0.5 on frozen dimensions. The use-cases seem reasonable, and it
>    seems like an easy-ish way to get them. Allowing ufuncs to natively support
>    subarray types might be a tidier solution, but that could come down the road
> Indeed, they are not mutually exclusive. My guess would be that the use
cases would be somewhat different.

>    - I’m -1 on optional dimensions: they seem to legitimize creating many
>    overloads of gufuncs. I’m already not a fan of how matmul has special cases
>    for lower dimensions that don’t generalize well. To me, the best way to
>    handle matmul would be to use the proposed __array_function__ to
>    handle the shape-based special-case dispatching, either by:
>       - Inserting dimensions, and calling the true gufunc
>       np.linalg.matmul_2d (which is a function I’d like direct access to
>       anyway).
>       - Dispatching to one of four ufuncs
> I must admit I wish that `@` was just pure matrix multiplication...  But
otherwise agree with Stephan as optional dimensions being the least-bad

Aside: do agree we should think about how to expose the `linalg` gufuncs.

>    - Broadcasting dimensions:
>       - I know you’re not suggesting this but: enabling broadcasting
>       unconditionally for all gufuncs would be a bad idea, masking linalg bugs.
>       (although einsum does support broadcasting…)
> Indeed, definitely *not* suggesting that!

>    -
>       - Does it really need a per-dimension flag, rather than a global
>       one? Can you give a case where that’s useful?
> Mostly simply that the implementation is easier given the optional
dimensions... Also, it has the benefit of being clear what the function can
handle by inspection of the signature, i.e., it self-documents better (one
of my main arguments in favour of frozen dimensions...).

>    -
>       - If we’d already made all_equal a gufunc, I’d be +1 on adding
>       broadcasting support to it
>       - I’m -0.5 on the all_equal path in the first place. I think we
>       either should have a more generic approach to combined ufuncs, or just
>       declare them numbas job.
> I am working on and off on a way to generically chain ufuncs (goal would
be to auto-create an inner loop that calls all the chained ufuncs loops in
turn). Not sure that short-circuiting will be all that easy.

I actually quite like the all_equal ufunc, but it is in part because I
remember discovering how painfully slow (a==b).all() was (and still have a
place where I would use it if it existed). And it does fit in the
(admittedly vague) plans to try to make `.reduce` a gufunc.

>    -
>       - Can you come up with a broadcasting use-case that isn’t just
>       chaining a reduction with a broadcasting ufunc?
> Perhaps the use is that it allows people to write gufuncs that are like
such functions... Absent a mechanism to chain ufuncs, more complicated
gufuncs are currently the easiest way to get fast more complicated algebra.

But perhaps a putative

weighted_mean(y, sigma) -> mean, sigma_mean

is a decent example? Its signature would be


but then you're forced to give individual sigmas for each point. With


you are no longer forced to do that (though the case of all y being the
same is less than useful here... I did at some point have an implementation
that worked by core dimension of each argument, but ended up feeling it was
not worth the extra complication)

-- Marten
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20180611/fcc7b4b2/attachment.html>

More information about the NumPy-Discussion mailing list