<div dir="ltr"><div>The following write up was triggered by issue <a href="https://github.com/numpy/numpy/issues/7265">#7265</a>, you may want to review that discussion for context:</div><div><br></div>In the docs, ufunc reduction operations are explained as "cumulatively applying the ufunc along an axis." This basically limits just about any reduction operation to ufuncs with two inpus and one output. Another take on the same reduction idea is that, if the inputs are <font face="monospace, monospace">a</font> and <font face="monospace, monospace">b</font> and the output is <font face="monospace, monospace">c</font>, <font face="monospace, monospace">c</font> is the result of updating <font face="monospace, monospace">a</font> with the value in <font face="monospace, monospace">b</font>. With this idea in mind, in e.g.<div><br></div><div><font face="monospace, monospace">    add(a, b) -> c = a + b</font>,</div><div><br></div><div><font face="monospace, monospace">c</font> would be thought of as the updated value of <font face="monospace, monospace">a</font>, after applying <font face="monospace, monospace">b</font> to it. One would expect that, for any registered loop suitable for this interpretation, the types of <font face="monospace, monospace">a</font> and <font face="monospace, monospace">c</font> would be identical, but b could be of some other type. As an example consider</div><div><br></div><div><font face="monospace, monospace">    count(a, b) -> c = a + 1</font>,</div><div><br></div><div>where <font face="monospace, monospace">a</font> and <font face="monospace, monospace">c</font> would typically be <font face="monospace, monospace">intp</font>, but <font face="monospace, monospace">b</font> could have just about any type.</div><div><br></div><div>The power of this description of ufuncs suited for reduction is that it is very easy to generalize beyond the current "two inputs, one output" restriction. E.g. a "sum of squared differences" ufunc defined as:</div><div><br></div><div><font face="monospace, monospace">    ssqd(sd2, x, y) -> sd2 += (x - y)**2</font>,</div><div><br></div><div>could be used to compute squared euclidean distances between vectors doing:</div><div><br></div><div><font face="monospace, monospace">    ssqd.reduce((x, y))</font>,<br></div><div><br></div><div>without the memory overhead of current available approaches.</div><div><br></div><div>In general, a reduction of a ufunc with <font face="monospace, monospace">m</font> inputs and <font face="monospace, monospace">n</font> outputs, <font face="monospace, monospace">m > n</font>, would produce <font face="monospace, monospace">n</font> results out of <font face="monospace, monospace">m - n</font> inputs. Such a ufunc would have to define a suitable identity to initialize each of those <font face="monospace, monospace">m - n</font> outputs at the beginning of any reduction, rather than the single identity ufuncs now hold.</div><div><br></div><div>It can be argued that this generalization of reduction is just a redefinition of a subset of what gufuncs can do, e.g. a "sum of squared differences" gufunc with signature <font face="monospace, monospace">(n),(n)->()</font> would do the same as the above reduction, probably faster. And it is not clear that getting <font face="monospace, monospace">accumulate</font> and <font face="monospace, monospace">reduceat</font> for free, or reduction over multiple axes, justifies the extra complexity. Especially since it is likely that something similarly powerful could be built on top of gufuncs.</div><div><br></div><div>Where an approach such as this would shine is if we had a reduction method that did not require a single strided memory segment to act on. Setting aside the generalization of reduction described above, a <font face="monospace, monospace">groupby</font> or <font face="monospace, monospace">reduceby</font> method would take an array of <font face="monospace, monospace">values</font> and an array of <font face="monospace, monospace">groups</font>, and compute a reduction value for each of the unique <font face="monospace, monospace">groups</font>. With the above generalizations, one could compute, e.g. the variance over each group by doing something like:</div><div><br></div><div><font face="monospace, monospace">    # how exactly does keepdims work here is not clear at all!</font></div><div><font face="monospace, monospace">    counts = count.reduceby(values, groups, keepdims=True)</font></div><div><font face="monospace, monospace">    mean = np.add.reduceby(values, groups, keepdims=True) / counts</font></div><div><font face="monospace, monospace">    var = np.ssqd.reduceby((values, mean), groups) / counts</font></div><div><br></div><div>I am not fully sure of whether this is just useful for this particular little example of computing variance using a ufunc method that doesn't exist, or if it is valid in other settings too. Implementing this would add complexity to an already complex system, so it better be good for something!</div><div><br></div><div>One of the weakest points I see is the need to rely on a set of predefined identities to start the reduction. Think of trying to generalize reduction to gufuncs, another worthy effort, and take matrix multiplication with signature <font face="monospace, monospace">(i,j),(j,k)->(i,k)</font>. In the current paradigm you would impose that both inputs and the output parameters be interchangeable, which leads rather naturally to the condition <font face="monospace, monospace">i = j = k</font> for this reduction to be doable. But with the new paradigm you only need that the first input and output be interchangeable, which only provides you with <font face="monospace, monospace">j = k</font> as a condition. And while this is a more general behavior, the questions of what value to set <font face="monospace, monospace">i</font> to in the output, and what to init the output to, when trying to do a reduction over a stack of square arrays, make it fairly unusable. So a backup to the "two inputs, one output, initialize to some item in the reduction array" would have to be kept in place</div><div><br></div><div>It is also important to note that expanding reduction to gufuncs would probably require that we impose some iteration order guarantees on ourselves, as the poster child for this (matrix multiplication) is not in general a commutative operation. Would we want to do this to ourselves? Would the features today justify the burden for ever after?</div><div><br></div><div>Any thoughts on where, if anywhere, to go with this, are very welcome!</div><div><div><br></div><div>Jaime</div><div><br></div>-- <br><div>(\__/)<br>( O.o)<br>( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.</div>
</div></div>