[Numpy-discussion] Proposed new feature for numpy.einsum: repeated output subscripts as diagonal

Eelco Hoogendoorn hoogendoorn.eelco at gmail.com
Fri Aug 15 15:09:02 EDT 2014


here is a snippet I extracted from a project with similar aims (integrating
the functionality of einsum and numexpr, actually)

Not much to it, but in case someone needs a reminder on how to use striding
tricks:
http://pastebin.com/kQNySjcj


On Fri, Aug 15, 2014 at 5:20 PM, Eelco Hoogendoorn <
hoogendoorn.eelco at gmail.com> wrote:

> 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/1eb533a8/attachment.html>


More information about the NumPy-Discussion mailing list