[Numpy-discussion] Consider improving numpy.outer's behavior with zero-dimensional vectors

Sebastian Berg sebastian at sipsolutions.net
Fri Apr 17 15:54:31 EDT 2015 

On Fr, 2015-04-17 at 20:56 +0200, Sebastian Berg wrote:
> On Fr, 2015-04-17 at 12:40 -0400, josef.pktd at gmail.com wrote:
> > On Fri, Apr 17, 2015 at 12:16 PM, Neil Girdhar <mistersheik at gmail.com> wrote:
> > >
> > >
> > > On Fri, Apr 17, 2015 at 12:09 PM, <josef.pktd at gmail.com> wrote:
> > >>
> > >> On Fri, Apr 17, 2015 at 11:22 AM, Neil Girdhar <mistersheik at gmail.com>
> > >> wrote:
> > >> >
> > >> >
> > >> > On Fri, Apr 17, 2015 at 10:47 AM, <josef.pktd at gmail.com> wrote:
> > >> >>
> > >> >> On Fri, Apr 17, 2015 at 10:07 AM, Sebastian Berg
> > >> >> <sebastian at sipsolutions.net> wrote:
> > >> >> > On Do, 2015-04-16 at 15:28 -0700, Matthew Brett wrote:
> > >> >> >> Hi,
> > >> >> >>
> > >> >> > <snip>
> > >> >> >>
> > >> >> >> So, how about a slight modification of your proposal?
> > >> >> >>
> > >> >> >> 1) Raise deprecation warning for np.outer for non 1D arrays for a
> > >> >> >> few
> > >> >> >> versions, with depraction in favor of np.multiply.outer, then
> > >> >> >> 2) Raise error for np.outer on non 1D arrays
> > >> >> >>
> > >> >> >
> > >> >> > I think that was Neil's proposal a bit earlier, too. +1 for it in any
> > >> >> > case, since at least for the moment I doubt outer is used a lot for
> > >> >> > non
> > >> >> > 1-d arrays. Possible step 3) make it work on higher dims after a long
> > >> >> > period.
> > >> >>
> > >> >> sounds ok to me
> > >> >>
> > >> >> Some random comments of what I remember or guess in terms of usage
> > >> >>
> > >> >> I think there are at most very few np.outer usages with 2d or higher
> > >> >> dimension.
> > >> >> (statsmodels has two models that switch between 2d and 1d
> > >> >> parameterization where we don't use outer but it has similar
> > >> >> characteristics. However, we need to control the ravel order, which
> > >> >> IIRC is Fortran)
> > >> >>
> > >> >> The current behavior of 0-D scalars in the initial post might be
> > >> >> useful if a numpy function returns a scalar instead of a 1-D array in
> > >> >> size=1. np.diag which is a common case, doesn't return a scalar (in my
> > >> >> version of numpy).
> > >> >>
> > >> >> I don't know any use case where I would ever want to have the 2d
> > >> >> behavior of np.multiply.outer.
> > >> >
> > >>
> > >> I only understand part of your example, but it looks similar to what
> > >> we are doing in statsmodels.
> > >>
> > >> >
> > >> > My use case is pretty simple.  Given an input vector x, and a weight
> > >> > matrix
> > >> > W, and a model y=Wx, I calculate the gradient of the loss L with respect
> > >> > W.
> > >> > It is the outer product of x with the vector of gradients dL/dy.  So the
> > >> > code is simply:
> > >> >
> > >> > W -= outer(x, dL_by_dy)
> > >>
> > >> if you sum/subtract over all the values, isn't this the same as
> > >> np.dot(x, dL_by_dy)
> > >>
> > >
> > > What?  Matrix subtraction is element-wise:
> > >
> > > In [1]: x = np.array([2,3,4])
> > >
> > > In [2]: dL_by_dy = np.array([7,9])
> > >
> > > In [5]: W = np.zeros((3, 2))
> > >
> > > In [6]: W -= np.outer(x, dL_by_dy)
> > >
> > > In [7]: W
> > > Out[7]:
> > > array([[-14., -18.],
> > >        [-21., -27.],
> > >        [-28., -36.]])
> > 
> > 
> > Ok, different use case
> > 
> > mine are more like variations on the following
> > 
> > >>> a1 = np.arange(18).reshape(6,3)
> > >>> a2 = np.arange(12).reshape(6, 2)
> > >>> index = [1, 2, 5]
> > 
> > 
> > text book version
> > >>> np.sum([np.outer(a1[i], a2[i]) for i in index], 0)
> > array([[180, 204],
> >        [196, 223],
> >        [212, 242]])
> > 
> > simpler
> > >>> np.dot(a1[index].T, a2[index])
> > array([[180, 204],
> >        [196, 223],
> >        [212, 242]])
> > 
> > 
> > >
> > >> >
> > >> > Sometimes, I have some x_indices and y_indices.  Now I want to do:
> > >> >
> > >> > W[x_indices, y_indices] -= outer(x[x_indices], dL_by_dy[y_indices])
> > >> >
> > >> > Unfortunately, if x_indices or y_indices are "int" or slice in some way
> > >> > that
> > >> > removes a dimension, the left side will have fewer dimensions than the
> > >> > right.  np.multipy.outer does the right thing without the ugly cases:
> > >> >
> > >> > if isinstance(x_indices, int): … # ugly hacks follow.
> > >>
> > >> My usual hacks are either to use np.atleast_1d or np.atleast_1d or
> > >> np.squeeze if there is shape mismatch in some cases.
> > >
> > >
> > > Yes, but in this case, the left side is the problem, which has too few
> > > dimensions.  So atleast_1d doesn't work.  I was conditionally squeezing, but
> > > that is extremely ugly.  Especially if you're conditionally squeezing based
> > > on both x_indices and y_indices.
> > 
> > I don't remember if I ever used something like this
> > 
> > >>> a1[0, 1]
> > 1
> > >>> a1[np.atleast_1d(0), np.atleast_1d(1)]
> > array([1])
> > 
> > >>> a1[np.atleast_1d(0), np.atleast_1d(1)] = [[100]]
> > 
> > >>> a1[0, 1]  = [[100]]
> > Traceback (most recent call last):
> >   File "<pyshell#314>", line 1, in <module>
> >     a1[0, 1]  = [[100]]
> > ValueError: setting an array element with a sequence.
> > 
> 
> Hehe, yeah, that difference. But if you really want that, you can
> usually do a1[0, 1, ...] if you don't mind the ugliness.
> 

Though actually I think I would usually prefer the other way around:
a1[None, None, 0, 1] = [[100]]
or instead
a1[0, 1] = np.array([[100]])[0, 0]


> > Josef
> > 
> > 
> > >
> > >>
> > >>
> > >> >
> > >> >> I guess we will or would have applications for outer along an axis,
> > >> >> for example if x.shape = (100, 10), then we have
> > >> >> x[:,None, :] * x[:, :, None]     (I guess)
> > >> >> Something like this shows up reasonably often in econometrics as
> > >> >> "Outer Product". However in most cases we can avoid constructing this
> > >> >> matrix and get the final results in a more memory efficient or faster
> > >> >> way.
> > >> >> (example an array of covariance matrices)
> > >> >
> > >> >
> > >> > Not sure I see this.  outer(a, b) should return something that has
> > >> > shape:
> > >> > (a.shape + b.shape).  If you're doing it "along an axis", you mean
> > >> > you're
> > >> > reshuffling the resulting shape vector?
> > >>
> > >> No I'm not reshaping the full tensor product.
> > >>
> > >> It's a vectorized version of looping over independent outer products
> > >>
> > >> np.array([outer(xi, yi) for xi,yi in zip(x, y)])
> > >> (which I would never use with outer)
> > >>
> > >> but I have code that works similar for a reduce (or reduce_at) loop over
> > >> this.
> > >>
> > >> Josef
> > >>
> > >>
> > >> >>
> > >> >>
> > >> >> Josef
> > >> >>
> > >> >>
> > >> >>
> > >> >>
> > >> >> >
> > >> >> > - Sebastian
> > >> >> >
> > >> >> >
> > >> >> >> Best,
> > >> >> >>
> > >> >> >> Matthew
> > >> >> >> _______________________________________________
> > >> >> >> NumPy-Discussion mailing list
> > >> >> >> NumPy-Discussion at scipy.org
> > >> >> >> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> > >> >> >>
> > >> >> >
> > >> >> >
> > >> >> > _______________________________________________
> > >> >> > NumPy-Discussion mailing list
> > >> >> > NumPy-Discussion at scipy.org
> > >> >> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> > >> >> >
> > >> >> _______________________________________________
> > >> >> NumPy-Discussion mailing list
> > >> >> NumPy-Discussion at scipy.org
> > >> >> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> > >> >
> > >> >
> > >> >
> > >> > _______________________________________________
> > >> > NumPy-Discussion mailing list
> > >> > NumPy-Discussion at scipy.org
> > >> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> > >> >
> > >> _______________________________________________
> > >> NumPy-Discussion mailing list
> > >> NumPy-Discussion at scipy.org
> > >> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> > >
> > >
> > >
> > > _______________________________________________
> > > NumPy-Discussion mailing list
> > > NumPy-Discussion at scipy.org
> > > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> > >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > NumPy-Discussion at scipy.org
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> 
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 819 bytes
Desc: This is a digitally signed message part
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20150417/b2879d6b/attachment.sig>

More information about the NumPy-Discussion mailing list