[Numpy-discussion] Proposed new feature for numpy.einsum: repeated output subscripts as diagonal
Eelco Hoogendoorn
hoogendoorn.eelco at gmail.com
Fri Aug 15 11:20:16 EDT 2014
Well, there is the numpy-API C level, and then there is the arcane macro C
level. The two might as well be a completely different language.
Indeed, it should be doing something similar for the inputs. Actually, I
think I wrote a wrapper around einsum/numexpr once that performed this
generalized indexing once... ill see if I can dig that up.
On Fri, Aug 15, 2014 at 5:01 PM, Sebastian Berg <sebastian at sipsolutions.net>
wrote:
> On Fr, 2014-08-15 at 16:42 +0200, Eelco Hoogendoorn wrote:
> > Agreed; this addition occurred to me as well. Note that the
> > implemenatation should be straightforward: just allocate an enlarged
> > array, use some striding logic to construct the relevant view, and let
> > einsums internals act on the view. hopefully, you wont even have to
> > touch the guts of einsum at the C level, because id say that isn't for
> > the faint of heart...
> >
>
> I am not sure if einsum isn't pure C :). But even if, it should be doing
> something identical already for duplicate indices on the inputs...
>
> - Sebastian
>
> >
> > On Fri, Aug 15, 2014 at 3:53 PM, Sebastian Berg
> > <sebastian at sipsolutions.net> wrote:
> > On Do, 2014-08-14 at 12:42 -0700, Stephan Hoyer wrote:
> > > I think this would be very nice addition.
> > >
> > >
> > > On Thu, Aug 14, 2014 at 12:21 PM, Benjamin Root
> > <ben.root at ou.edu>
> > > wrote:
> > > You had me at Kronecker delta... :-) +1
> > >
> >
> >
> > Sounds good to me. I don't see a reason for not relaxing the
> > restriction, unless there is some technical issue, but I doubt
> > that.
> >
> > - Sebastian
> >
> > >
> > >
> > > On Thu, Aug 14, 2014 at 3:07 PM, Pierre-Andre Noel
> > > <noel.pierre.andre at gmail.com> wrote:
> > > (I created issue 4965 earlier today on this
> > topic, and
> > > I have been
> > > advised to email to this mailing list to
> > discuss
> > > whether it is a good
> > > idea or not. I include my original post
> > as-is,
> > > followed by additional
> > > comments.)
> > >
> > > I think that the following new feature would
> > make
> > > `numpy.einsum` even
> > > more powerful/useful/awesome than it already
> > is.
> > > Moreover, the change
> > > should not interfere with existing code, it
> > would
> > > preserve the
> > > "minimalistic" spirit of `numpy.einsum`, and
> > the new
> > > functionality would
> > > integrate in a seamless/intuitive manner for
> > the
> > > users.
> > >
> > > In short, the new feature would allow for
> > repeated
> > > subscripts to appear
> > > in the "output" part of the `subscripts`
> > parameter
> > > (i.e., on the
> > > right-hand side of `->`). The corresponding
> > dimensions
> > > in the resulting
> > > `ndarray` would only be filled along their
> > diagonal,
> > > leaving the off
> > > diagonal entries to the default value for
> > this `dtype`
> > > (typically zero).
> > > Note that the current behavior is to raise
> > an
> > > exception when repeated
> > > output subscripts are being used.
> > >
> > > This is simplest to describe with an example
> > involving
> > > the dual behavior
> > > of `numpy.diag`.
> > >
> > > ```python
> > > # Extracting the diagonal of a 2-D array.
> > > A = arange(16).reshape(4,4)
> > > print(diag(A)) # Output: [ 0 5 10 15 ]
> > > print(einsum('ii->i', A)) # Same as previous
> > line
> > > (current behavior).
> > >
> > > # Constructing a diagonal 2-D array.
> > > v = arange(4)
> > > print(diag(v)) # Output: [[0 0 0 0] [0 1 0
> > 0] [0 0 2
> > > 0] [0 0 0 3]]
> > > print(einsum('i->ii', v)) # New behavior
> > would be same
> > > as previous line.
> > > # The current behavior of the previous line
> > is to
> > > raise an exception.
> > > ```
> > >
> > > By opposition to `numpy.diag`, the approach
> > > generalizes to higher
> > > dimensions: `einsum('iii->i', A)` extracts
> > the
> > > diagonal of a 3-D array,
> > > and `einsum('i->iii', v)` would build a
> > diagonal 3-D
> > > array.
> > >
> > > The proposed behavior really starts to shine
> > in more
> > > intricate cases.
> > >
> > > ```python
> > > # Dummy values, these should be
> > probabilities to make
> > > sense below.
> > > P_w_ab = arange(24).reshape(3,2,4)
> > > P_y_wxab = arange(144).reshape(3,3,2,2,4)
> > >
> > > # With the proposed behavior, the following
> > two lines
> > > should be equivalent.
> > > P_xyz_ab = einsum('wab,xa,ywxab,zy->xyzab',
> > P_w_ab,
> > > eye(2), P_y_wxab,
> > > eye(3))
> > > also_P_xyz_ab = einsum('wab,ywaab->ayyab',
> > P_w_ab,
> > > P_y_wxab)
> > > ```
> > >
> > > If this is not convincing enough, replace
> > `eye(2)` by
> > > `eye(P_w_ab.shape[1])` and replace `eye(3)`
> > by
> > > `eye(P_y_wxab.shape[0])`,
> > > then imagine more dimensions and repeated
> > indices...
> > > The new notation
> > > would allow for crisper codes and reduce the
> > > opportunities for dumb
> > > mistakes.
> > >
> > > For those who wonder, the above computation
> > amounts to
> > > $P(X=x,Y=y,Z=z|A=a,B=b) = \sum_w P(W=w|
> > A=a,B=b) P(X=x|
> > > A=a)
> > > P(Y=y|W=w,X=x,A=a,B=b) P(Z=z|Y=y)$ with
> > $P(X=x|A=a)=
> > > \delta_{xa}$ and
> > > $P(Z=z|Y=y)=\delta_{zy}$ (using LaTeX
> > notation, and
> > > $\delta_{ij}$ is
> > > [Kronecker's
> > >
> > delta](http://en.wikipedia.org/wiki/Kronecker_delta)).
> > >
> > > (End of original post.)
> > >
> > > I have been told by @jaimefrio that "The
> > best way of
> > > getting a new
> > > feature into numpy is putting it in
> > yourself." Hence,
> > > if discussions
> > > here do reveal that this is a good idea,
> > then I may
> > > give a try at coding
> > > it myself. However, I currently know nothing
> > of the
> > > inner workings of
> > > numpy/ndarray/einsum, and I have higher
> > priorities
> > > right now. This means
> > > that it could take a long while before I
> > contribute
> > > any code, if I ever
> > > do. Hence, if anyone feels like doing it,
> > feel free to
> > > do so!
> > >
> > > Also, I am aware that storing a lot of zeros
> > in an
> > > `ndarray` may not, a
> > > priori, be a desirable avenue. However,
> > there are
> > > times where you have
> > > to do it: think of `numpy.eye` as an
> > example. In my
> > > case of application,
> > > I use such diagonal structures in the
> > initialization
> > > of an `ndarray`
> > > which is later updated through an iterative
> > process.
> > > After these
> > > iterations, most of the zeros will be gone.
> > Do other
> > > people see a use
> > > for such capabilities?
> > >
> > > Thank you for your time and have a nice day.
> > >
> > > Sincerely,
> > >
> > > Pierre-André Noël
> > >
> > _______________________________________________
> > > 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 --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20140815/f315ae48/attachment.html>
More information about the NumPy-Discussion
mailing list