<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Sep 17, 2014 at 2:30 AM, Sebastian Berg <span dir="ltr"><<a href="mailto:sebastian@sipsolutions.net" target="_blank">sebastian@sipsolutions.net</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class=""><div class="h5">On Di, 2014-09-16 at 16:51 -0400, Nathaniel Smith wrote:<br>
> On Tue, Sep 16, 2014 at 4:31 PM, Jaime Fernández del Río<br>
> <<a href="mailto:jaime.frio@gmail.com">jaime.frio@gmail.com</a>> wrote:<br>
> > If it is a bug, it is an extended one, because it is the same behavior of<br>
> > einsum:<br>
> ><br>
> >>>> np.einsum('i,i', [1,1,1], [1])<br>
> > 3<br>
> >>>> np.einsum('i,i', [1,1,1], [1,1])<br>
> > Traceback (most recent call last):<br>
> >   File "<stdin>", line 1, in <module><br>
> > ValueError: operands could not be broadcast together with remapped shapes<br>
> > [origi<br>
> > nal->remapped]: (3,)->(3,) (2,)->(2,)<br>
> ><br>
> > And I think it is a conscious design decision, there is a comment about<br>
> > broadcasting missing core dimensions here:<br>
> ><br>
> > <a href="https://github.com/numpy/numpy/blob/master/numpy/core/src/umath/ufunc_object.c#L1940" target="_blank">https://github.com/numpy/numpy/blob/master/numpy/core/src/umath/ufunc_object.c#L1940</a><br>
><br>
> "intentional" and "sensible" are not always the same thing :-). That<br>
> said, it isn't totally obvious to me what the correct behaviour for<br>
> einsum is in this case.<br>
><br>
> > and the code makes it very explicit that input argument dimensions with the<br>
> > same label are broadcast to a common shape, see here:<br>
> ><br>
> > <a href="https://github.com/numpy/numpy/blob/master/numpy/core/src/umath/ufunc_object.c#L1956" target="_blank">https://github.com/numpy/numpy/blob/master/numpy/core/src/umath/ufunc_object.c#L1956</a><br>
> ><br>
> > I kind of expect numpy to broadcast whenever possible, so this doesn't feel<br>
> > wrong to me.<br>
><br>
> The case Chuck is talking about is like if we allowed matrix<br>
> multiplication between an array with shape (n, 1) with an array with<br>
> (k, m), because (n, 1) can be broadcast to (n, k). This feels VERY<br>
> wrong to me, will certainly hide many bugs, and is definitely not how<br>
> it works right now (for np.dot, anyway; apparently it does work that<br>
> way for the brand-new gufunc np.linalg.matrix_multiply, but this must<br>
> be an accident).<br>
<br>
</div></div>Agreed, the only argument to not change it right away would be being<br>
afraid of breaking user code abusing the kind of thing Josef mentioned.<br>
<span class=""><font color="#888888"><br></font></span></blockquote><div><br></div><div>It *is* a big change. I think of the gufuncs as working with matrices and vectors, the array version of the matrix class. In that case the signature shapes must be preserved and there should be no broadcasting within the signature. The matrix class itself fails in that regard:<br><br>In [1]: matrix(eye(3)) + matrix([[1,1,1]])<br>Out[1]: <br>matrix([[ 2.,  1.,  1.],<br>        [ 1.,  2.,  1.],<br>        [ 1.,  1.,  2.]])<br><br></div><div>Which is all wrong for a matrix type.<br><br>It would also be nice if the order could be made part of the signature as DGEMM and friends like one of the argument axis to be contiguous, but I don't see a clean way to do that. The gufuncs do have an order parameter which should probably default to 'C'  if the arrays/vectors are stacked. I think the default is currently 'K'. Hmm, we could make 'K' refer to the last one or two dimensions in the inputs. OTOH, that isn't needed for types not handled by BLAS. Or it could be handled in the inner loops.<br><br></div><div><snip><br><br></div><div>Chuck<br></div></div></div></div>