[Numpy-discussion] Bringing order to higher dimensional operations

[Nathaniel Smith]

On Thu, Jul 18, 2013 at 2:23 PM, Sebastian Berg
<sebastian at sipsolutions.net> wrote:
> On Thu, 2013-07-18 at 13:52 +0100, Nathaniel Smith wrote:
>> Hi all,
>>
> <snip>
>>
>> So:
>>
>> QUESTION 1: does that sound right: that in a perfect world, the
>> current gufunc convention would be the only one, and that's what we
>> should work towards, at least in the cases where that's possible?
>>
>
> Sounds right to me, ufunc/gufunc broadcasting assumes the "inner"
> dimensions are the right-most. Since we are normally in C-order arrays,
> this also seems the sensible way if you consider the memory layout.
>
>> QUESTION 2: Assuming that's right, it would be *really nice* if we
>> could at least get np.dot onto our new convention, for consistency
>> with the rest of np.linalg, and to allow it to be overridden. I'm sort
>> of terrified to touch np.dot's API, but the only cases where it would
>> act differently is when *both* arguments have *3 or more dimensions*,
>> and I guess there are very very few users who fall into that category.
>> So maybe we could start raising some big FutureWarnings for this case
>> in the next release, and eventually switch?
>>
>
> It is noble to try to get do to use the gufunc convention, but if you
> look at the new gufunc linalg functions, they already have to have some
> weird tricks in the case of np.linalg.solve.
> It is so difficult because of the fact that dot is basically a
> combination of many functions:
>   o vector * vector -> vector
>   o vector * matrix -> matrix (add dimensions to vector on right)
>   o matrix * vector -> matrix (add dimensions to vector on left)
>   o matrix * matrix -> matrix
> plus scalar cases.

Oh ugh, I forgot about all those special cases.

> I somewhat believe we should not touch dot, or deprecate anything but
> the most basic dot functionality. Then we can point to matrix_multiply,
> inner1d, etc. which are gufuncs (even if they are not exposed at this
> time). The whole dance that is already done for np.linalg.solve right
> now is not pretty there, and it will be worse for dot. Because dot is
> basically overloaded, marrying it with the broadcasting machinery in a
> general way is impossible.

While it would be kind of nice if we could eventually make
isinstance(np.dot, np.ufunc) be True, I'm not so bothered if it
remains a wrapper around gufuncs (like the np.linalg wrappers
currently are). Most of these special cases, while ugly, can be
handled perfectly well by this sort of mechanism. What I'm most
bothered about is pseudo-outer case:
  np.dot(array with ndim >2, array with ndim >2)
This simply *can't* be emulated with a gufunc. And as long as that's
true it's going to be very hard to get 'dot' to play along with the
general ufunc machinery.

So that's specifically the case I was talking about in question 2.

>> (I'm even more terrified of trying to mess with np.sum or
>> np.add.reduce, so I'll leave that alone for now -- maybe we're just
>> stuck with them. And at least they do already go through the ufunc
>> machinery.)
>
> I did not understand where the inconsistency/problem for the reductions
> is.

What I mean is: Suppose we wrote a gufunc for 'sum', where the
intrinsic operation took a vector and returned a scalar. (E.g. we want
to implement one of the specialized algorithms for vector summation,
like Kahan summation, which can be more accurate than applying scalar
addition repeatedly.)

Then we'd have:

np.sum(ones((2, 3))).shape == ()
np.add.reduce(ones((2, 3))).shape == (3,)
gufunc_sum(ones((2, 3))).shape == (2,)

These are three names for exactly the same underlying function... but
they all have different defaults for how they vectorize.

-n