Mailman 3 March 2014 - NumPy-Discussion

[help needed] associativity and precedence of '@'
by Nathaniel Smith 20 Mar '14

20 Mar '14

Hi all, Here's the main blocker for adding a matrix multiply operator '@' to Python: we need to decide what we think its precedence and associativity should be. I'll explain what that means so we're on the same page, and what the choices are, and then we can all argue about it. But even better would be if we could get some data to guide our decision, and this would be a lot easier if some of you all can help; I'll suggest some ways you might be able to do that. So! Precedence and left- versus right-associativity. If you already know what these are you can skim down until you see CAPITAL LETTERS. We all know what precedence is. Code like this: a + b * c gets evaluated as: a + (b * c) because * has higher precedence than +. It "binds more tightly", as they say. Python's complete precedence able is here: http://docs.python.org/3/reference/expressions.html#operator-precedence Associativity, in the parsing sense, is less well known, though it's just as important. It's about deciding how to evaluate code like this: a * b * c Do we use a * (b * c) # * is "right associative" or (a * b) * c # * is "left associative" ? Here all the operators have the same precedence (because, uh... they're the same operator), so precedence doesn't help. And mostly we can ignore this in day-to-day life, because both versions give the same answer, so who cares. But a programming language has to pick one (consider what happens if one of those objects has a non-default __mul__ implementation). And of course it matters a lot for non-associative operations like a - b - c or a / b / c So when figuring out order of evaluations, what you do first is check the precedence, and then if you have multiple operators next to each other with the same precedence, you check their associativity. Notice that this means that if you have different operators that share the same precedence level (like + and -, or * and /), then they have to all have the same associativity. All else being equal, it's generally considered nice to have fewer precedence levels, because these have to be memorized by users. Right now in Python, every precedence level is left-associative, except for '**'. If you write these formulas without any parentheses, then what the interpreter will actually execute is: (a * b) * c (a - b) - c (a / b) / c but a ** (b ** c) Okay, that's the background. Here's the question. We need to decide on precedence and associativity for '@'. In particular, there are three different options that are interesting: OPTION 1 FOR @: Precedence: same as * Associativity: left My shorthand name for it: "same-left" (yes, very creative) This means that if you don't use parentheses, you get: a @ b @ c -> (a @ b) @ c a * b @ c -> (a * b) @ c a @ b * c -> (a @ b) * c OPTION 2 FOR @: Precedence: more-weakly-binding than * Associativity: right My shorthand name for it: "weak-right" This means that if you don't use parentheses, you get: a @ b @ c -> a @ (b @ c) a * b @ c -> (a * b) @ c a @ b * c -> a @ (b * c) OPTION 3 FOR @: Precedence: more-tightly-binding than * Associativity: right My shorthand name for it: "tight-right" This means that if you don't use parentheses, you get: a @ b @ c -> a @ (b @ c) a * b @ c -> a * (b @ c) a @ b * c -> (a @ b) * c We need to pick which of which options we think is best, based on whatever reasons we can think of, ideally more than "hmm, weak-right gives me warm fuzzy feelings" ;-). (In principle the other 2 possible options are tight-left and weak-left, but there doesn't seem to be any argument in favor of either, so we'll leave them out of the discussion.) Some things to consider: * and @ are actually not associative (in the math sense) with respect to each other, i.e., (a * b) @ c and a * (b @ c) in general give different results when 'a' is not a scalar. So considering the two expressions 'a * b @ c' and 'a @ b * c', we can see that each of these three options gives produces different results in some cases. "Same-left" is the easiest to explain and remember, because it's just, "@ acts like * and /". So we already have to know the rule in order to understand other non-associative expressions like a / b / c or a - b - c, and it'd be nice if the same rule applied to things like a * b @ c so we only had to memorize *one* rule. (Of course there's ** which uses the opposite rule, but I guess everyone internalized that one in secondary school; that's not true for * versus @.) This is definitely the default we should choose unless we have a good reason to do otherwise. BUT: there might indeed be a good reason to do otherwise, which is the whole reason this has come up. Consider: Mat1 @ Mat2 @ vec Obviously this will execute much more quickly if we do Mat1 @ (Mat2 @ vec) because that results in two cheap matrix-vector multiplies, while (Mat1 @ Mat2) @ vec starts out by doing an expensive matrix-matrix multiply. So: maybe @ should be right associative, so that we get the fast behaviour without having to use explicit parentheses! /If/ these kinds of expressions are common enough that having to remember to put explicit parentheses in all the time is more of a programmer burden than having to memorize a special associativity rule for @. Obviously Mat @ Mat @ vec is more common than vec @ Mat @ Mat, but maybe they're both so rare that it doesn't matter in practice -- I don't know. Also, if we do want @ to be right associative, then I can't think of any clever reasons to prefer weak-right over tight-right, or vice-versa. For the scalar multiplication case, I believe both options produce the same result in the same amount of time. For the non-scalar case, they give different answers. Do people have strong intuitions about what expressions like a * b @ c a @ b * c should do actually? (I'm guessing not, but hey, you never know.) And, while intuition is useful, it would be really *really* nice to be basing these decisions on more than *just* intuition, since whatever we decide will be subtly influencing the experience of writing linear algebra code in Python for the rest of time. So here's where I could use some help. First, of course, if you have any other reasons why one or the other of these options is better, then please share! But second, I think we need to know something about how often the Mat @ Mat @ vec type cases arise in practice. How often do non-scalar * and np.dot show up in the same expression? How often does it look like a * np.dot(b, c), and how often does it look like np.dot(a * b, c)? How often do we see expressions like np.dot(np.dot(a, b), c), and how often do we see expressions like np.dot(a, np.dot(b, c))? This would really help guide the debate. I don't have this data, and I'm not sure the best way to get it. A super-fancy approach would be to write a little script that uses the 'ast' module to count things automatically. A less fancy approach would be to just pick some code you've written, or a well-known package, grep through for calls to 'dot', and make notes on what you see. (An advantage of the less-fancy approach is that as a human you might be able to tell the difference between scalar and non-scalar *, or check whether it actually matters what order the 'dot' calls are done in.) -n -- Nathaniel J. Smith Postdoctoral researcher - Informatics - University of Edinburgh http://vorpus.org

16 87

OT: job opening at STScI
by Perry Greenfield 20 Mar '14

20 Mar '14

The Science Software Branch at the Space Telescope Science Institute is seeking a developer that enjoys and is good at digging into and understanding the details of Python internals and its libraries--particularly those related to scientific computing--in order to support of our development of software tools for Python and the astronomical community, particularly for the Hubble Space Telescope and the next telescope under construction (JWST). Details on the position and how to apply can be found at this link: https://rn11.ultipro.com/SPA1004/JobBoard/JobDetails.aspx?__ID=*5ECC2DFF670… (applications are still being taken despite the date given…) Perry Greenfield

1 0

[RFC] should we argue for a matrix power operator, @@?
by Nathaniel Smith 19 Mar '14

19 Mar '14

Hi all, Here's the second thread for discussion about Guido's concerns about PEP 465. The issue here is that PEP 465 as currently written proposes two new operators, @ for matrix multiplication and @@ for matrix power (analogous to * and **): http://legacy.python.org/dev/peps/pep-0465/ The main thing we care about of course is @; I pushed for including @@ because I thought it was nicer to have than not, and I thought the analogy between * and ** might make the overall package more appealing to Guido's aesthetic sense. It turns out I was wrong :-). Guido is -0 on @@, but willing to be swayed if we think it's worth the trouble to make a solid case. Note that question now is *not*, how will @@ affect the reception of @. @ itself is AFAICT a done deal, regardless of what happens with @@. For this discussion let's assume @ can be taken for granted, and that we can freely choose to either add @@ or not add @@ to the language. The question is: which do we think makes Python a better language (for us and in general)? Some thoughts to start us off: Here are the interesting use cases for @@ that I can think of: - 'vector @@ 2' gives the squared Euclidean length (because it's the same as vector @ vector). Kind of handy. - 'matrix @@ n' of course gives the matrix power, which is of marginal use but does come in handy sometimes, e.g., when looking at graph connectivity. - 'matrix @@ -1' provides a very transparent notation for translating textbook formulas (with all their inverses) into code. It's a bit unhelpful in practice, because (a) usually you should use solve(), and (b) 'matrix @@ -1' is actually more characters than 'inv(matrix)'. But sometimes transparent notation may be important. (And in some cases, like using numba or theano or whatever, 'matrix @@ -1 @ foo' could be compiled into a call to solve() anyway.) (Did I miss any?) In practice it seems to me that the last use case is the one that's might matter a lot practice, but then again, it might not -- I'm not sure. For example, does anyone who teaches programming with numpy have a feeling about whether the existence of '@@ -1' would make a big difference to you and your students? (Alan? I know you were worried about losing the .I attribute on matrices if switching to ndarrays for teaching -- given that ndarray will probably not get a .I attribute, how much would the existence of @@ -1 affect you?) On a more technical level, Guido is worried about how @@'s precedence should work (and this is somewhat related to the other thread about @'s precedence and associativity, because he feels that if we end up giving @ and * different precedence, then that makes it much less clear what to do with @@, and reduces the strength of the */**/@/@@ analogy). In particular, if we want to argue for @@ then we'll need to figure out what expressions like a @@ b @@ c and a ** b @@ c and a @@ b ** c should do. A related question is what @@ should do if given an array as its right argument. In the current PEP, only integers are accepted, which rules out a bunch of the more complicated cases like a @@ b @@ c (at least assuming @@ is right-associative, like **, and I can't see why you'd want anything else). OTOH, in the brave new gufunc world, it technically would make sense to define @@ as being a gufunc with signature (m,m),()->(m,m), and the way gufuncs work this *would* allow the "power" to be an array -- for example, we'd have: mat = randn(m, m) pow = range(n) result = gufunc_matrix_power(mat, pow) assert result.shape == (n, m, m) for i in xrange(n): assert np.all(result[i, :, :] == mat ** i) In this case, a @@ b @@ c would at least be a meaningful expression to write. OTOH it would be incredibly bizarre and useless, so probably no-one would ever write it. As far as these technical issues go, my guess is that the correct rule is that @@ should just have the same precedence and the same (right) associativity as **, and in practice no-one will ever write stuff like a @@ b @@ c. But if we want to argue for @@ we need to come to some consensus or another here. It's also possible the answer is "ugh, these issues are too complicated, we should defer this until later when we have more experience with @ and gufuncs and stuff". After all, I doubt anyone else will swoop in and steal @@ to mean something else! OTOH, if e.g. there's a strong feeling that '@@ -1' will make a big difference in pedagogical contexts, then putting that off for years might be a mistake. -n -- Nathaniel J. Smith Postdoctoral researcher - Informatics - University of Edinburgh http://vorpus.org

16 19

PR with changes to triangular array functions
by Jaime Fernández del Río 18 Mar '14

18 Mar '14

I submitted a PR that makes some improvements to the numpy functions dealing with triangular arrays. Aside from a general speed-up of about 2x for most functions, there are some minor changes to the public API. In case anyone is concerned about them, here's a list: * 'np.tri' now accepts a boolean 'invert' kwarg that is equivalent to '1 - np.tri' only faster. * 'np.mask_indices' is no longer used by any of the triangular array functions. While it is part of the public API, it is not even mentioned in the documentation AFAICT. It may be a candidate for deprecation IMO. * 'np.tril_indices' and 'np.triu_indices' now accept an 'm' kwarg to indicate the number of columns of the array, so they are no longer restricted to square arrays. The weird thing is that, to preserve the order of the existing arguments, the signature is '(n, k=0, m=None)', while other similar functions, such as 'np.tri', have signature '(n, m=None, k=0)'. * 'np.triu_indices_from' and 'np.tril_indices_from' now also accept rectangular arrays. The PR can be found here: https://github.com/numpy/numpy/pull/4509 Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

2 1

_gufuncs_linalg module
by Jay Bourque 18 Mar '14

18 Mar '14

I was just about to submit some pull requests for fixes to the _gufuncs_linalg module and discovered that it no longer exists. It looks like it was removed in this commit<https://github.com/numpy/numpy/commit/0f516827dd081625b8b2262297be57ac855a9…>. Is there any reason why it was removed without any apparent discussion? It looks like it was originally added in this PR<https://github.com/numpy/numpy/pull/3220>after a looooong discussion. Thanks, -Jay

2 2

Re: [Numpy-discussion] NumPy-Discussion Digest, Vol 90, Issue 56
by Colin J. Williams 17 Mar '14

17 Mar '14

Julian, I can see the need to recognize both column and row vectors, but why not with np.matrix? I can see no need for a new operator and hope to be able to comment more fully on PEP 465 in a few days. Colin W. On 17-Mar-2014 7:19 PM, numpy-discussion-request(a)scipy.org wrote: > Send NumPy-Discussion mailing list submissions to > numpy-discussion(a)scipy.org > > To subscribe or unsubscribe via the World Wide Web, visit > http://mail.scipy.org/mailman/listinfo/numpy-discussion > or, via email, send a message with subject or body 'help' to > numpy-discussion-request(a)scipy.org > > You can reach the person managing the list at > numpy-discussion-owner(a)scipy.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of NumPy-Discussion digest..." > > > Today's Topics: > > 1. Re: [help needed] associativity and precedence of '@' > (Nathaniel Smith) > 2. Re: GSoC project: draft of proposal (Julian Taylor) > 3. Re: [help needed] associativity and precedence of '@' > (Christophe Bal) > 4. Re: [help needed] associativity and precedence of '@' > (Alexander Belopolsky) > 5. Re: [help needed] associativity and precedence of '@' (Bago) > 6. Re: [help needed] associativity and precedence of '@' > (Christophe Bal) > 7. Re: [help needed] associativity and precedence of '@' > (Christophe Bal) > 8. Re: [help needed] associativity and precedence of '@' > (Nathaniel Smith) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Mon, 17 Mar 2014 22:02:33 +0000 > From: Nathaniel Smith <njs(a)pobox.com> > Subject: Re: [Numpy-discussion] [help needed] associativity and > precedence of '@' > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > Message-ID: > <CAPJVwB=zBazN+fiYWJeiWOL=4a9Bf2XGxJGoTT8GFTt-kDUDZw(a)mail.gmail.com> > Content-Type: text/plain; charset=UTF-8 > > On Mon, Mar 17, 2014 at 9:38 PM, Christophe Bal <projetmbc(a)gmail.com> wrote: >> Here is the translation. ;-) >> >> Hello, >> and what about something like that ? >> >> a @ b @ c -> (a @ b) @ c >> a * b @ c -> (a * b) @ c >> a @ b * c -> a @ (b * c) >> >> Easy to remember: the *-product has priority regarding to the @-product, and >> we just do @-product from left to right. > In the terminology we've been using in this thread, this is "weak-left". > >> An advantage of this is that most parsers do analyze from left to right. >> >> So I really think that it is a better choice than the weak-right one. > We've mostly ignored this option because of assuming that if we want > left-associativity, we should go with "same-left" instead of > "weak-left". Same-left is: > > a @ b @ c -> (a @ b) @ c > a * b @ c -> (a * b) @ c > a @ b * c -> (a @ b) * c > > i.e., even more left-to-right than weak-left :-) > > Do you think weak-left is better than same-left? >

1 0

GSoC project: draft of proposal
by Leo Mao 17 Mar '14

17 Mar '14

Hi, The attachment is my draft of proposal. The project is "vector math library integration". I think I need some feedback to make it solider. Any comment will be appreciated. Thanks in advance. Regards, Leo Mao

9 16

ANN: Scipy 0.14.0 beta 1 release
by Ralf Gommers 16 Mar '14

16 Mar '14

Hi, I'm pleased to announce the availability of the first beta release of Scipy0.14.0. Please try this beta and report any issues on the scipy-dev mailing list. Source tarballs, binaries and the full release notes can be found at http://sourceforge.net/projects/scipy/files/scipy/0.14.0b1/. Part of the release notes copied below. A big thank you to everyone who contributed to this release! Ralf SciPy 0.14.0 is the culmination of 8 months of hard work. It contains many new features, numerous bug-fixes, improved test coverage and better documentation. There have been a number of deprecations and API changes in this release, which are documented below. All users are encouraged to upgrade to this release, as there are a large number of bug-fixes and optimizations. Moreover, our development attention will now shift to bug-fix releases on the 0.14.x branch, and on adding new features on the master branch. This release requires Python 2.6, 2.7 or 3.2-3.4 and NumPy 1.5.1 or greater. New features ============ ``scipy.interpolate`` improvements ---------------------------------- A new wrapper function `scipy.interpolate.interpn` for interpolation on regular grids has been added. `interpn` supports linear and nearest-neighbor interpolation in arbitrary dimensions and spline interpolation in two dimensions. Faster implementations of piecewise polynomials in power and Bernstein polynomial bases have been added as `scipy.interpolate.PPoly` and `scipy.interpolate.BPoly`. New users should use these in favor of `scipy.interpolate.PiecewisePolynomial`. `scipy.interpolate.interp1d` now accepts non-monotonic inputs and sorts them. If performance is critical, sorting can be turned off by using the new ``assume_sorted`` keyword. Functionality for evaluation of bivariate spline derivatives in ``scipy.interpolate`` has been added. The new class `scipy.interpolate.Akima1DInterpolator` implements the piecewise cubic polynomial interpolation scheme devised by H. Akima. Functionality for fast interpolation on regular, unevenly spaced grids in arbitrary dimensions has been added as `scipy.interpolate.RegularGridInterpolator` . ``scipy.linalg`` improvements ----------------------------- The new function `scipy.linalg.dft` computes the matrix of the discrete Fourier transform. A condition number estimation function for matrix exponential, `scipy.linalg.expm_cond`, has been added. ``scipy.optimize`` improvements ------------------------------- A set of benchmarks for optimize, which can be run with ``optimize.bench()``, has been added. `scipy.optimize.curve_fit` now has more controllable error estimation via the ``absolute_sigma`` keyword. Support for passing custom minimization methods to ``optimize.minimize()`` and ``optimize.minimize_scalar()`` has been added, currently useful especially for combining ``optimize.basinhopping()`` with custom local optimizer routines. ``scipy.stats`` improvements ---------------------------- A new class `scipy.stats.multivariate_normal` with functionality for multivariate normal random variables has been added. A lot of work on the ``scipy.stats`` distribution framework has been done. Moment calculations (skew and kurtosis mainly) are fixed and verified, all examples are now runnable, and many small accuracy and performance improvements for individual distributions were merged. The new function `scipy.stats.anderson_ksamp` computes the k-sample Anderson-Darling test for the null hypothesis that k samples come from the same parent population. ``scipy.signal`` improvements ----------------------------- ``scipy.signal.iirfilter`` and related functions to design Butterworth, Chebyshev, elliptical and Bessel IIR filters now all use pole-zero ("zpk") format internally instead of using transformations to numerator/denominator format. The accuracy of the produced filters, especially high-order ones, is improved significantly as a result. The new function `scipy.signal.vectorstrength` computes the vector strength, a measure of phase synchrony, of a set of events. ``scipy.special`` improvements ------------------------------ The functions `scipy.special.boxcox` and `scipy.special.boxcox1p`, which compute the Box-Cox transformation, have been added. ``scipy.sparse`` improvements ----------------------------- - Significant performance improvement in CSR, CSC, and DOK indexing speed. - When using Numpy >= 1.9 (to be released in MM 2014), sparse matrices function correctly when given to arguments of ``np.dot``, ``np.multiply`` and other ufuncs. With earlier Numpy and Scipy versions, the results of such operations are undefined and usually unexpected. - Sparse matrices are no longer limited to ``2^31`` nonzero elements. They automatically switch to using 64-bit index data type for matrices containing more elements. User code written assuming the sparse matrices use int32 as the index data type will continue to work, except for such large matrices. Code dealing with larger matrices needs to accept either int32 or int64 indices. Deprecated features =================== ``anneal`` ---------- The global minimization function `scipy.optimize.anneal` is deprecated. All users should use the `scipy.optimize.basinhopping` function instead. ``scipy.stats`` --------------- ``randwcdf`` and ``randwppf`` functions are deprecated. All users should use distribution-specific ``rvs`` methods instead. Probability calculation aliases ``zprob``, ``fprob`` and ``ksprob`` are deprecated. Use instead the ``sf`` methods of the corresponding distributions or the ``special`` functions directly. ``scipy.interpolate`` --------------------- ``PiecewisePolynomial`` class is deprecated. Backwards incompatible changes ============================== scipy.special.lpmn ------------------ ``lpmn`` no longer accepts complex-valued arguments. A new function ``clpmn`` with uniform complex analytic behavior has been added, and it should be used instead. scipy.sparse.linalg ------------------- Eigenvectors in the case of generalized eigenvalue problem are normalized to unit vectors in 2-norm, rather than following the LAPACK normalization convention. The deprecated UMFPACK wrapper in ``scipy.sparse.linalg`` has been removed due to license and install issues. If available, ``scikits.umfpack`` is still used transparently in the ``spsolve`` and ``factorized`` functions. Otherwise, SuperLU is used instead in these functions. scipy.stats ----------- The deprecated functions ``glm``, ``oneway`` and ``cmedian`` have been removed from ``scipy.stats``. ``stats.scoreatpercentile`` now returns an array instead of a list of percentiles. scipy.interpolate ----------------- The API for computing derivatives of a monotone piecewise interpolation has changed: if `p` is a ``PchipInterpolator`` object, `p.derivative(der)` returns a callable object representing the derivative of `p`. For in-place derivatives use the second argument of the `__call__` method: `p(0.1, der=2)` evaluates the second derivative of `p` at `x=0.1`. The method `p.derivatives` has been removed. Authors ======= * Marc Abramowitz + * andbo + * Vincent Arel-Bundock + * Petr Baudis + * Max Bolingbroke * François Boulogne * Matthew Brett * Lars Buitinck * Evgeni Burovski * CJ Carey + * Thomas A Caswell + * Pawel Chojnacki + * Phillip Cloud + * Stefano Costa + * David Cournapeau * Dapid + * Matthieu Dartiailh + * Christoph Deil + * Jörg Dietrich + * endolith * Francisco de la Peña + * Ben FrantzDale + * Jim Garrison + * André Gaul * Christoph Gohlke * Ralf Gommers * Robert David Grant * Alex Griffing * Blake Griffith * Yaroslav Halchenko * Andreas Hilboll * Kat Huang * Gert-Ludwig Ingold * jamestwebber + * Dorota Jarecka + * Todd Jennings + * Thouis (Ray) Jones * Juan Luis Cano Rodríguez * ktritz + * Jacques Kvam + * Eric Larson + * Justin Lavoie + * Denis Laxalde * Jussi Leinonen + * lemonlaug + * Tim Leslie * Alain Leufroy + * George Lewis + * Max Linke + * Brandon Liu + * Benny Malengier + * Matthias Kümmerer + * Cimarron Mittelsteadt + * Eric Moore * Andrew Nelson + * Niklas Hambüchen + * Joel Nothman + * Clemens Novak * Emanuele Olivetti + * Stefan Otte + * peb + * Josef Perktold * pjwerneck * Andrew Sczesnak + * poolio * Jérôme Roy + * Carl Sandrock + * Shauna + * Fabrice Silva * Daniel B. Smith * Patrick Snape + * Thomas Spura + * Jacob Stevenson * Julian Taylor * Tomas Tomecek * Richard Tsai * Joris Vankerschaver + * Pauli Virtanen * Warren Weckesser A total of 78 people contributed to this release. People with a "+" by their names contributed a patch for the first time. This list of names is automatically generated, and may not be fully complete.

1 0

It looks like Py 3.5 will include a dedicated infix matrix multiply operator
by Nathaniel Smith 16 Mar '14

16 Mar '14

Well, that was fast. Guido says he'll accept the addition of '@' as an infix operator for matrix multiplication, once some details are ironed out: https://mail.python.org/pipermail/python-ideas/2014-March/027109.html http://legacy.python.org/dev/peps/pep-0465/ Specifically, we need to figure out whether we want to make an argument for a matrix power operator ("@@"), and what precedence/associativity we want '@' to have. I'll post two separate threads to get feedback on those in an organized way -- this is just a heads-up. -n -- Nathaniel J. Smith Postdoctoral researcher - Informatics - University of Edinburgh http://vorpus.org

12 16

Re: [Numpy-discussion] NumPy-Discussion Digest, Vol 90, Issue 45
by Colin J. Williams 16 Mar '14

16 Mar '14

I would like to see the case made for @. Yes, I know that Guido has accepted the idea, but he has changed his mind before. The PEP seems neutral to retaining both np.matrix and @. Nearly ten years ago, Tim Peters <http://legacy.python.org/dev/peps/pep-0020/> gave us: /There should be one-- and preferably only one --obvious way to do it. / W/e now have: / /C= A * B C becomes an instance of the Matrix class (m, p) When A and B are matrices a matrix of (m, n) and (n, p) respectively. Actually, the rules are a little more general than the above. / The PEP proposes that /C= /A @ B where the types or classes of A, B and C are not clear. We also have A.I for the inverse, for the square matrix) or A.T for the transpose of a matrix. One way is recommended in the Zen of Python, of the two, which is the obvious way? Colin W. / / On 15-Mar-2014 9:25 PM, numpy-discussion-request(a)scipy.org wrote: > Send NumPy-Discussion mailing list submissions to > numpy-discussion(a)scipy.org > > To subscribe or unsubscribe via the World Wide Web, visit > http://mail.scipy.org/mailman/listinfo/numpy-discussion > or, via email, send a message with subject or body 'help' to > numpy-discussion-request(a)scipy.org > > You can reach the person managing the list at > numpy-discussion-owner(a)scipy.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of NumPy-Discussion digest..." > > > Today's Topics: > > 1. Re: [help needed] associativity and precedence of '@' > (josef.pktd(a)gmail.com) > 2. Re: [RFC] should we argue for a matrix power operator, @@? > (josef.pktd(a)gmail.com) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Sat, 15 Mar 2014 21:20:40 -0400 > From: josef.pktd(a)gmail.com > Subject: Re: [Numpy-discussion] [help needed] associativity and > precedence of '@' > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > Message-ID: > <CAMMTP+Ahag9fN3XPtS4uDRThBknVXzudc0G8TtJ7G3w3dWbBWw(a)mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > On Fri, Mar 14, 2014 at 11:41 PM, Nathaniel Smith <njs(a)pobox.com> wrote: > >> Hi all, >> >> Here's the main blocker for adding a matrix multiply operator '@' to >> Python: we need to decide what we think its precedence and associativity >> should be. I'll explain what that means so we're on the same page, and what >> the choices are, and then we can all argue about it. But even better would >> be if we could get some data to guide our decision, and this would be a lot >> easier if some of you all can help; I'll suggest some ways you might be >> able to do that. >> >> So! Precedence and left- versus right-associativity. If you already know >> what these are you can skim down until you see CAPITAL LETTERS. >> >> We all know what precedence is. Code like this: >> a + b * c >> gets evaluated as: >> a + (b * c) >> because * has higher precedence than +. It "binds more tightly", as they >> say. Python's complete precedence able is here: >> http://docs.python.org/3/reference/expressions.html#operator-precedence >> >> Associativity, in the parsing sense, is less well known, though it's just >> as important. It's about deciding how to evaluate code like this: >> a * b * c >> Do we use >> a * (b * c) # * is "right associative" >> or >> (a * b) * c # * is "left associative" >> ? Here all the operators have the same precedence (because, uh... they're >> the same operator), so precedence doesn't help. And mostly we can ignore >> this in day-to-day life, because both versions give the same answer, so who >> cares. But a programming language has to pick one (consider what happens if >> one of those objects has a non-default __mul__ implementation). And of >> course it matters a lot for non-associative operations like >> a - b - c >> or >> a / b / c >> So when figuring out order of evaluations, what you do first is check the >> precedence, and then if you have multiple operators next to each other with >> the same precedence, you check their associativity. Notice that this means >> that if you have different operators that share the same precedence level >> (like + and -, or * and /), then they have to all have the same >> associativity. All else being equal, it's generally considered nice to have >> fewer precedence levels, because these have to be memorized by users. >> >> Right now in Python, every precedence level is left-associative, except >> for '**'. If you write these formulas without any parentheses, then what >> the interpreter will actually execute is: >> (a * b) * c >> (a - b) - c >> (a / b) / c >> but >> a ** (b ** c) >> >> Okay, that's the background. Here's the question. We need to decide on >> precedence and associativity for '@'. In particular, there are three >> different options that are interesting: >> >> OPTION 1 FOR @: >> Precedence: same as * >> Associativity: left >> My shorthand name for it: "same-left" (yes, very creative) >> >> This means that if you don't use parentheses, you get: >> a @ b @ c -> (a @ b) @ c >> a * b @ c -> (a * b) @ c >> a @ b * c -> (a @ b) * c >> >> OPTION 2 FOR @: >> Precedence: more-weakly-binding than * >> Associativity: right >> My shorthand name for it: "weak-right" >> >> This means that if you don't use parentheses, you get: >> a @ b @ c -> a @ (b @ c) >> a * b @ c -> (a * b) @ c >> a @ b * c -> a @ (b * c) >> >> OPTION 3 FOR @: >> Precedence: more-tightly-binding than * >> Associativity: right >> My shorthand name for it: "tight-right" >> >> This means that if you don't use parentheses, you get: >> a @ b @ c -> a @ (b @ c) >> a * b @ c -> a * (b @ c) >> a @ b * c -> (a @ b) * c >> >> We need to pick which of which options we think is best, based on whatever >> reasons we can think of, ideally more than "hmm, weak-right gives me warm >> fuzzy feelings" ;-). (In principle the other 2 possible options are >> tight-left and weak-left, but there doesn't seem to be any argument in >> favor of either, so we'll leave them out of the discussion.) >> >> Some things to consider: >> >> * and @ are actually not associative (in the math sense) with respect to >> each other, i.e., (a * b) @ c and a * (b @ c) in general give different >> results when 'a' is not a scalar. So considering the two expressions 'a * b >> @ c' and 'a @ b * c', we can see that each of these three options gives >> produces different results in some cases. >> >> "Same-left" is the easiest to explain and remember, because it's just, "@ >> acts like * and /". So we already have to know the rule in order to >> understand other non-associative expressions like a / b / c or a - b - c, >> and it'd be nice if the same rule applied to things like a * b @ c so we >> only had to memorize *one* rule. (Of course there's ** which uses the >> opposite rule, but I guess everyone internalized that one in secondary >> school; that's not true for * versus @.) This is definitely the default we >> should choose unless we have a good reason to do otherwise. >> >> BUT: there might indeed be a good reason to do otherwise, which is the >> whole reason this has come up. Consider: >> Mat1 @ Mat2 @ vec >> Obviously this will execute much more quickly if we do >> Mat1 @ (Mat2 @ vec) >> because that results in two cheap matrix-vector multiplies, while >> (Mat1 @ Mat2) @ vec >> starts out by doing an expensive matrix-matrix multiply. So: maybe @ >> should be right associative, so that we get the fast behaviour without >> having to use explicit parentheses! /If/ these kinds of expressions are >> common enough that having to remember to put explicit parentheses in all >> the time is more of a programmer burden than having to memorize a special >> associativity rule for @. Obviously Mat @ Mat @ vec is more common than vec >> @ Mat @ Mat, but maybe they're both so rare that it doesn't matter in >> practice -- I don't know. >> >> Also, if we do want @ to be right associative, then I can't think of any >> clever reasons to prefer weak-right over tight-right, or vice-versa. For >> the scalar multiplication case, I believe both options produce the same >> result in the same amount of time. For the non-scalar case, they give >> different answers. Do people have strong intuitions about what expressions >> like >> a * b @ c >> a @ b * c >> should do actually? (I'm guessing not, but hey, you never know.) >> >> And, while intuition is useful, it would be really *really* nice to be >> basing these decisions on more than *just* intuition, since whatever we >> decide will be subtly influencing the experience of writing linear algebra >> code in Python for the rest of time. So here's where I could use some help. >> First, of course, if you have any other reasons why one or the other of >> these options is better, then please share! But second, I think we need to >> know something about how often the Mat @ Mat @ vec type cases arise in >> practice. How often do non-scalar * and np.dot show up in the same >> expression? How often does it look like a * np.dot(b, c), and how often >> does it look like np.dot(a * b, c)? How often do we see expressions like >> np.dot(np.dot(a, b), c), and how often do we see expressions like np.dot(a, >> np.dot(b, c))? This would really help guide the debate. I don't have this >> data, and I'm not sure the best way to get it. A super-fancy approach would >> be to write a little script that uses the 'ast' module to count things >> automatically. A less fancy approach would be to just pick some code you've >> written, or a well-known package, grep through for calls to 'dot', and make >> notes on what you see. (An advantage of the less-fancy approach is that as >> a human you might be able to tell the difference between scalar and >> non-scalar *, or check whether it actually matters what order the 'dot' >> calls are done in.) >> >> -n >> >> -- >> Nathaniel J. Smith >> Postdoctoral researcher - Informatics - University of Edinburgh >> http://vorpus.org >> >> >> _______________________________________________ >> NumPy-Discussion mailing list >> NumPy-Discussion(a)scipy.org >> http://mail.scipy.org/mailman/listinfo/numpy-discussion >> >> > I'm in favor of same-left because it's the easiest to remember. > with scalar factors it is how I read formulas. > > Both calculating dot @ first or calculating elementwise * first sound > logical, but I wouldn't know which should go first. (My "feeling" would be > @ first.) > > > two cases I remembered in statsmodels > H = np.dot(results.model.pinv_wexog, scale[:,None] * > results.model.pinv_wexog.T) > se = (exog * np.dot(covb, exog.T).T).sum(1) > > we are mixing * and dot pretty freely in all combinations AFAIR > > my guess is that I wouldn't trust any sequence without parenthesis for a > long time. > (and I don't trust a sequence of dots @ without parenthesis either, in our > applications.) > > x @ (W.T @ W) @ x ( W.shape = (10000, 5) ) > or > x * (W.T @ W) * x > > (w * x) @ x weighted sum of squares > > Josef >

2 1