1381 lines
67 KiB
Plaintext
1381 lines
67 KiB
Plaintext
PEP: 465
|
||
Title: A dedicated infix operator for matrix multiplication
|
||
Version: $Revision$
|
||
Last-Modified: $Date$
|
||
Author: Nathaniel J. Smith <njs@pobox.com>
|
||
Status: Final
|
||
Type: Standards Track
|
||
Content-Type: text/x-rst
|
||
Created: 20-Feb-2014
|
||
Python-Version: 3.5
|
||
Post-History: 13-Mar-2014
|
||
|
||
Abstract
|
||
========
|
||
|
||
This PEP proposes a new binary operator to be used for matrix
|
||
multiplication, called ``@``. (Mnemonic: ``@`` is ``*`` for
|
||
mATrices.)
|
||
|
||
|
||
Specification
|
||
=============
|
||
|
||
A new binary operator is added to the Python language, together
|
||
with the corresponding in-place version:
|
||
|
||
======= ========================= ===============================
|
||
Op Precedence/associativity Methods
|
||
======= ========================= ===============================
|
||
``@`` Same as ``*`` ``__matmul__``, ``__rmatmul__``
|
||
``@=`` n/a ``__imatmul__``
|
||
======= ========================= ===============================
|
||
|
||
No implementations of these methods are added to the builtin or
|
||
standard library types. However, a number of projects have reached
|
||
consensus on the recommended semantics for these operations; see
|
||
`Intended usage details`_ below for details.
|
||
|
||
For details on how this operator will be implemented in CPython, see
|
||
`Implementation details`_.
|
||
|
||
|
||
Motivation
|
||
==========
|
||
|
||
Executive summary
|
||
-----------------
|
||
|
||
In numerical code, there are two important operations which compete
|
||
for use of Python's ``*`` operator: elementwise multiplication, and
|
||
matrix multiplication. In the nearly twenty years since the Numeric
|
||
library was first proposed, there have been many attempts to resolve
|
||
this tension [#hugunin]_; none have been really satisfactory.
|
||
Currently, most numerical Python code uses ``*`` for elementwise
|
||
multiplication, and function/method syntax for matrix multiplication;
|
||
however, this leads to ugly and unreadable code in common
|
||
circumstances. The problem is bad enough that significant amounts of
|
||
code continue to use the opposite convention (which has the virtue of
|
||
producing ugly and unreadable code in *different* circumstances), and
|
||
this API fragmentation across codebases then creates yet more
|
||
problems. There does not seem to be any *good* solution to the
|
||
problem of designing a numerical API within current Python syntax --
|
||
only a landscape of options that are bad in different ways. The
|
||
minimal change to Python syntax which is sufficient to resolve these
|
||
problems is the addition of a single new infix operator for matrix
|
||
multiplication.
|
||
|
||
Matrix multiplication has a singular combination of features which
|
||
distinguish it from other binary operations, which together provide a
|
||
uniquely compelling case for the addition of a dedicated infix
|
||
operator:
|
||
|
||
* Just as for the existing numerical operators, there exists a vast
|
||
body of prior art supporting the use of infix notation for matrix
|
||
multiplication across all fields of mathematics, science, and
|
||
engineering; ``@`` harmoniously fills a hole in Python's existing
|
||
operator system.
|
||
|
||
* ``@`` greatly clarifies real-world code.
|
||
|
||
* ``@`` provides a smoother onramp for less experienced users, who are
|
||
particularly harmed by hard-to-read code and API fragmentation.
|
||
|
||
* ``@`` benefits a substantial and growing portion of the Python user
|
||
community.
|
||
|
||
* ``@`` will be used frequently -- in fact, evidence suggests it may
|
||
be used more frequently than ``//`` or the bitwise operators.
|
||
|
||
* ``@`` allows the Python numerical community to reduce fragmentation,
|
||
and finally standardize on a single consensus duck type for all
|
||
numerical array objects.
|
||
|
||
|
||
Background: What's wrong with the status quo?
|
||
---------------------------------------------
|
||
|
||
When we crunch numbers on a computer, we usually have lots and lots of
|
||
numbers to deal with. Trying to deal with them one at a time is
|
||
cumbersome and slow -- especially when using an interpreted language.
|
||
Instead, we want the ability to write down simple operations that
|
||
apply to large collections of numbers all at once. The *n-dimensional
|
||
array* is the basic object that all popular numeric computing
|
||
environments use to make this possible. Python has several libraries
|
||
that provide such arrays, with numpy being at present the most
|
||
prominent.
|
||
|
||
When working with n-dimensional arrays, there are two different ways
|
||
we might want to define multiplication. One is elementwise
|
||
multiplication::
|
||
|
||
[[1, 2], [[11, 12], [[1 * 11, 2 * 12],
|
||
[3, 4]] x [13, 14]] = [3 * 13, 4 * 14]]
|
||
|
||
and the other is `matrix multiplication`_:
|
||
|
||
.. _matrix multiplication: https://en.wikipedia.org/wiki/Matrix_multiplication
|
||
|
||
::
|
||
|
||
[[1, 2], [[11, 12], [[1 * 11 + 2 * 13, 1 * 12 + 2 * 14],
|
||
[3, 4]] x [13, 14]] = [3 * 11 + 4 * 13, 3 * 12 + 4 * 14]]
|
||
|
||
Elementwise multiplication is useful because it lets us easily and
|
||
quickly perform many multiplications on a large collection of values,
|
||
without writing a slow and cumbersome ``for`` loop. And this works as
|
||
part of a very general schema: when using the array objects provided
|
||
by numpy or other numerical libraries, all Python operators work
|
||
elementwise on arrays of all dimensionalities. The result is that one
|
||
can write functions using straightforward code like ``a * b + c / d``,
|
||
treating the variables as if they were simple values, but then
|
||
immediately use this function to efficiently perform this calculation
|
||
on large collections of values, while keeping them organized using
|
||
whatever arbitrarily complex array layout works best for the problem
|
||
at hand.
|
||
|
||
Matrix multiplication is more of a special case. It's only defined on
|
||
2d arrays (also known as "matrices"), and multiplication is the only
|
||
operation that has an important "matrix" version -- "matrix addition"
|
||
is the same as elementwise addition; there is no such thing as "matrix
|
||
bitwise-or" or "matrix floordiv"; "matrix division" and "matrix
|
||
to-the-power-of" can be defined but are not very useful, etc.
|
||
However, matrix multiplication is still used very heavily across all
|
||
numerical application areas; mathematically, it's one of the most
|
||
fundamental operations there is.
|
||
|
||
Because Python syntax currently allows for only a single
|
||
multiplication operator ``*``, libraries providing array-like objects
|
||
must decide: either use ``*`` for elementwise multiplication, or use
|
||
``*`` for matrix multiplication. And, unfortunately, it turns out
|
||
that when doing general-purpose number crunching, both operations are
|
||
used frequently, and there are major advantages to using infix rather
|
||
than function call syntax in both cases. Thus it is not at all clear
|
||
which convention is optimal, or even acceptable; often it varies on a
|
||
case-by-case basis.
|
||
|
||
Nonetheless, network effects mean that it is very important that we
|
||
pick *just one* convention. In numpy, for example, it is technically
|
||
possible to switch between the conventions, because numpy provides two
|
||
different types with different ``__mul__`` methods. For
|
||
``numpy.ndarray`` objects, ``*`` performs elementwise multiplication,
|
||
and matrix multiplication must use a function call (``numpy.dot``).
|
||
For ``numpy.matrix`` objects, ``*`` performs matrix multiplication,
|
||
and elementwise multiplication requires function syntax. Writing code
|
||
using ``numpy.ndarray`` works fine. Writing code using
|
||
``numpy.matrix`` also works fine. But trouble begins as soon as we
|
||
try to integrate these two pieces of code together. Code that expects
|
||
an ``ndarray`` and gets a ``matrix``, or vice-versa, may crash or
|
||
return incorrect results. Keeping track of which functions expect
|
||
which types as inputs, and return which types as outputs, and then
|
||
converting back and forth all the time, is incredibly cumbersome and
|
||
impossible to get right at any scale. Functions that defensively try
|
||
to handle both types as input and DTRT, find themselves floundering
|
||
into a swamp of ``isinstance`` and ``if`` statements.
|
||
|
||
PEP 238 split ``/`` into two operators: ``/`` and ``//``. Imagine the
|
||
chaos that would have resulted if it had instead split ``int`` into
|
||
two types: ``classic_int``, whose ``__div__`` implemented floor
|
||
division, and ``new_int``, whose ``__div__`` implemented true
|
||
division. This, in a more limited way, is the situation that Python
|
||
number-crunchers currently find themselves in.
|
||
|
||
In practice, the vast majority of projects have settled on the
|
||
convention of using ``*`` for elementwise multiplication, and function
|
||
call syntax for matrix multiplication (e.g., using ``numpy.ndarray``
|
||
instead of ``numpy.matrix``). This reduces the problems caused by API
|
||
fragmentation, but it doesn't eliminate them. The strong desire to
|
||
use infix notation for matrix multiplication has caused a number of
|
||
specialized array libraries to continue to use the opposing convention
|
||
(e.g., scipy.sparse, pyoperators, pyviennacl) despite the problems
|
||
this causes, and ``numpy.matrix`` itself still gets used in
|
||
introductory programming courses, often appears in StackOverflow
|
||
answers, and so forth. Well-written libraries thus must continue to
|
||
be prepared to deal with both types of objects, and, of course, are
|
||
also stuck using unpleasant funcall syntax for matrix multiplication.
|
||
After nearly two decades of trying, the numerical community has still
|
||
not found any way to resolve these problems within the constraints of
|
||
current Python syntax (see `Rejected alternatives to adding a new
|
||
operator`_ below).
|
||
|
||
This PEP proposes the minimum effective change to Python syntax that
|
||
will allow us to drain this swamp. It splits ``*`` into two
|
||
operators, just as was done for ``/``: ``*`` for elementwise
|
||
multiplication, and ``@`` for matrix multiplication. (Why not the
|
||
reverse? Because this way is compatible with the existing consensus,
|
||
and because it gives us a consistent rule that all the built-in
|
||
numeric operators also apply in an elementwise manner to arrays; the
|
||
reverse convention would lead to more special cases.)
|
||
|
||
So that's why matrix multiplication doesn't and can't just use ``*``.
|
||
Now, in the rest of this section, we'll explain why it nonetheless
|
||
meets the high bar for adding a new operator.
|
||
|
||
|
||
Why should matrix multiplication be infix?
|
||
------------------------------------------
|
||
|
||
Right now, most numerical code in Python uses syntax like
|
||
``numpy.dot(a, b)`` or ``a.dot(b)`` to perform matrix multiplication.
|
||
This obviously works, so why do people make such a fuss about it, even
|
||
to the point of creating API fragmentation and compatibility swamps?
|
||
|
||
Matrix multiplication shares two features with ordinary arithmetic
|
||
operations like addition and multiplication on numbers: (a) it is used
|
||
very heavily in numerical programs -- often multiple times per line of
|
||
code -- and (b) it has an ancient and universally adopted tradition of
|
||
being written using infix syntax. This is because, for typical
|
||
formulas, this notation is dramatically more readable than any
|
||
function call syntax. Here's an example to demonstrate:
|
||
|
||
One of the most useful tools for testing a statistical hypothesis is
|
||
the linear hypothesis test for OLS regression models. It doesn't
|
||
really matter what all those words I just said mean; if we find
|
||
ourselves having to implement this thing, what we'll do is look up
|
||
some textbook or paper on it, and encounter many mathematical formulas
|
||
that look like:
|
||
|
||
.. math::
|
||
|
||
S = (H \beta - r)^T (H V H^T)^{-1} (H \beta - r)
|
||
|
||
Here the various variables are all vectors or matrices (details for
|
||
the curious: [#lht]_).
|
||
|
||
Now we need to write code to perform this calculation. In current
|
||
numpy, matrix multiplication can be performed using either the
|
||
function or method call syntax. Neither provides a particularly
|
||
readable translation of the formula::
|
||
|
||
import numpy as np
|
||
from numpy.linalg import inv, solve
|
||
|
||
# Using dot function:
|
||
S = np.dot((np.dot(H, beta) - r).T,
|
||
np.dot(inv(np.dot(np.dot(H, V), H.T)), np.dot(H, beta) - r))
|
||
|
||
# Using dot method:
|
||
S = (H.dot(beta) - r).T.dot(inv(H.dot(V).dot(H.T))).dot(H.dot(beta) - r)
|
||
|
||
With the ``@`` operator, the direct translation of the above formula
|
||
becomes::
|
||
|
||
S = (H @ beta - r).T @ inv(H @ V @ H.T) @ (H @ beta - r)
|
||
|
||
Notice that there is now a transparent, 1-to-1 mapping between the
|
||
symbols in the original formula and the code that implements it.
|
||
|
||
Of course, an experienced programmer will probably notice that this is
|
||
not the best way to compute this expression. The repeated computation
|
||
of :math:`H \beta - r` should perhaps be factored out; and,
|
||
expressions of the form ``dot(inv(A), B)`` should almost always be
|
||
replaced by the more numerically stable ``solve(A, B)``. When using
|
||
``@``, performing these two refactorings gives us::
|
||
|
||
# Version 1 (as above)
|
||
S = (H @ beta - r).T @ inv(H @ V @ H.T) @ (H @ beta - r)
|
||
|
||
# Version 2
|
||
trans_coef = H @ beta - r
|
||
S = trans_coef.T @ inv(H @ V @ H.T) @ trans_coef
|
||
|
||
# Version 3
|
||
S = trans_coef.T @ solve(H @ V @ H.T, trans_coef)
|
||
|
||
Notice that when comparing between each pair of steps, it's very easy
|
||
to see exactly what was changed. If we apply the equivalent
|
||
transformations to the code using the .dot method, then the changes
|
||
are much harder to read out or verify for correctness::
|
||
|
||
# Version 1 (as above)
|
||
S = (H.dot(beta) - r).T.dot(inv(H.dot(V).dot(H.T))).dot(H.dot(beta) - r)
|
||
|
||
# Version 2
|
||
trans_coef = H.dot(beta) - r
|
||
S = trans_coef.T.dot(inv(H.dot(V).dot(H.T))).dot(trans_coef)
|
||
|
||
# Version 3
|
||
S = trans_coef.T.dot(solve(H.dot(V).dot(H.T)), trans_coef)
|
||
|
||
Readability counts! The statements using ``@`` are shorter, contain
|
||
more whitespace, can be directly and easily compared both to each
|
||
other and to the textbook formula, and contain only meaningful
|
||
parentheses. This last point is particularly important for
|
||
readability: when using function-call syntax, the required parentheses
|
||
on every operation create visual clutter that makes it very difficult
|
||
to parse out the overall structure of the formula by eye, even for a
|
||
relatively simple formula like this one. Eyes are terrible at parsing
|
||
non-regular languages. I made and caught many errors while trying to
|
||
write out the 'dot' formulas above. I know they still contain at
|
||
least one error, maybe more. (Exercise: find it. Or them.) The
|
||
``@`` examples, by contrast, are not only correct, they're obviously
|
||
correct at a glance.
|
||
|
||
If we are even more sophisticated programmers, and writing code that
|
||
we expect to be reused, then considerations of speed or numerical
|
||
accuracy might lead us to prefer some particular order of evaluation.
|
||
Because ``@`` makes it possible to omit irrelevant parentheses, we can
|
||
be certain that if we *do* write something like ``(H @ V) @ H.T``,
|
||
then our readers will know that the parentheses must have been added
|
||
intentionally to accomplish some meaningful purpose. In the ``dot``
|
||
examples, it's impossible to know which nesting decisions are
|
||
important, and which are arbitrary.
|
||
|
||
Infix ``@`` dramatically improves matrix code usability at all stages
|
||
of programmer interaction.
|
||
|
||
|
||
Transparent syntax is especially crucial for non-expert programmers
|
||
-------------------------------------------------------------------
|
||
|
||
A large proportion of scientific code is written by people who are
|
||
experts in their domain, but are not experts in programming. And
|
||
there are many university courses run each year with titles like "Data
|
||
analysis for social scientists" which assume no programming
|
||
background, and teach some combination of mathematical techniques,
|
||
introduction to programming, and the use of programming to implement
|
||
these mathematical techniques, all within a 10-15 week period. These
|
||
courses are more and more often being taught in Python rather than
|
||
special-purpose languages like R or Matlab.
|
||
|
||
For these kinds of users, whose programming knowledge is fragile, the
|
||
existence of a transparent mapping between formulas and code often
|
||
means the difference between succeeding and failing to write that code
|
||
at all. This is so important that such classes often use the
|
||
``numpy.matrix`` type which defines ``*`` to mean matrix
|
||
multiplication, even though this type is buggy and heavily
|
||
disrecommended by the rest of the numpy community for the
|
||
fragmentation that it causes. This pedagogical use case is, in fact,
|
||
the *only* reason ``numpy.matrix`` remains a supported part of numpy.
|
||
Adding ``@`` will benefit both beginning and advanced users with
|
||
better syntax; and furthermore, it will allow both groups to
|
||
standardize on the same notation from the start, providing a smoother
|
||
on-ramp to expertise.
|
||
|
||
|
||
But isn't matrix multiplication a pretty niche requirement?
|
||
-----------------------------------------------------------
|
||
|
||
The world is full of continuous data, and computers are increasingly
|
||
called upon to work with it in sophisticated ways. Arrays are the
|
||
lingua franca of finance, machine learning, 3d graphics, computer
|
||
vision, robotics, operations research, econometrics, meteorology,
|
||
computational linguistics, recommendation systems, neuroscience,
|
||
astronomy, bioinformatics (including genetics, cancer research, drug
|
||
discovery, etc.), physics engines, quantum mechanics, geophysics,
|
||
network analysis, and many other application areas. In most or all of
|
||
these areas, Python is rapidly becoming a dominant player, in large
|
||
part because of its ability to elegantly mix traditional discrete data
|
||
structures (hash tables, strings, etc.) on an equal footing with
|
||
modern numerical data types and algorithms.
|
||
|
||
We all live in our own little sub-communities, so some Python users
|
||
may be surprised to realize the sheer extent to which Python is used
|
||
for number crunching -- especially since much of this particular
|
||
sub-community's activity occurs outside of traditional Python/FOSS
|
||
channels. So, to give some rough idea of just how many numerical
|
||
Python programmers are actually out there, here are two numbers: In
|
||
2013, there were 7 international conferences organized specifically on
|
||
numerical Python [#scipy-conf]_ [#pydata-conf]_. At PyCon 2014, ~20%
|
||
of the tutorials appear to involve the use of matrices
|
||
[#pycon-tutorials]_.
|
||
|
||
To quantify this further, we used Github's "search" function to look
|
||
at what modules are actually imported across a wide range of
|
||
real-world code (i.e., all the code on Github). We checked for
|
||
imports of several popular stdlib modules, a variety of numerically
|
||
oriented modules, and various other extremely high-profile modules
|
||
like django and lxml (the latter of which is the #1 most downloaded
|
||
package on PyPI). Starred lines indicate packages which export array-
|
||
or matrix-like objects which will adopt ``@`` if this PEP is
|
||
approved::
|
||
|
||
Count of Python source files on Github matching given search terms
|
||
(as of 2014-04-10, ~21:00 UTC)
|
||
================ ========== =============== ======= ===========
|
||
module "import X" "from X import" total total/numpy
|
||
================ ========== =============== ======= ===========
|
||
sys 2374638 63301 2437939 5.85
|
||
os 1971515 37571 2009086 4.82
|
||
re 1294651 8358 1303009 3.12
|
||
numpy ************** 337916 ********** 79065 * 416981 ******* 1.00
|
||
warnings 298195 73150 371345 0.89
|
||
subprocess 281290 63644 344934 0.83
|
||
django 62795 219302 282097 0.68
|
||
math 200084 81903 281987 0.68
|
||
threading 212302 45423 257725 0.62
|
||
pickle+cPickle 215349 22672 238021 0.57
|
||
matplotlib 119054 27859 146913 0.35
|
||
sqlalchemy 29842 82850 112692 0.27
|
||
pylab *************** 36754 ********** 41063 ** 77817 ******* 0.19
|
||
scipy *************** 40829 ********** 28263 ** 69092 ******* 0.17
|
||
lxml 19026 38061 57087 0.14
|
||
zlib 40486 6623 47109 0.11
|
||
multiprocessing 25247 19850 45097 0.11
|
||
requests 30896 560 31456 0.08
|
||
jinja2 8057 24047 32104 0.08
|
||
twisted 13858 6404 20262 0.05
|
||
gevent 11309 8529 19838 0.05
|
||
pandas ************** 14923 *********** 4005 ** 18928 ******* 0.05
|
||
sympy 2779 9537 12316 0.03
|
||
theano *************** 3654 *********** 1828 *** 5482 ******* 0.01
|
||
================ ========== =============== ======= ===========
|
||
|
||
These numbers should be taken with several grains of salt (see
|
||
footnote for discussion: [#github-details]_), but, to the extent they
|
||
can be trusted, they suggest that ``numpy`` might be the single
|
||
most-imported non-stdlib module in the entire Pythonverse; it's even
|
||
more-imported than such stdlib stalwarts as ``subprocess``, ``math``,
|
||
``pickle``, and ``threading``. And numpy users represent only a
|
||
subset of the broader numerical community that will benefit from the
|
||
``@`` operator. Matrices may once have been a niche data type
|
||
restricted to Fortran programs running in university labs and military
|
||
clusters, but those days are long gone. Number crunching is a
|
||
mainstream part of modern Python usage.
|
||
|
||
In addition, there is some precedence for adding an infix operator to
|
||
handle a more-specialized arithmetic operation: the floor division
|
||
operator ``//``, like the bitwise operators, is very useful under
|
||
certain circumstances when performing exact calculations on discrete
|
||
values. But it seems likely that there are many Python programmers
|
||
who have never had reason to use ``//`` (or, for that matter, the
|
||
bitwise operators). ``@`` is no more niche than ``//``.
|
||
|
||
|
||
So ``@`` is good for matrix formulas, but how common are those really?
|
||
----------------------------------------------------------------------
|
||
|
||
We've seen that ``@`` makes matrix formulas dramatically easier to
|
||
work with for both experts and non-experts, that matrix formulas
|
||
appear in many important applications, and that numerical libraries
|
||
like numpy are used by a substantial proportion of Python's user base.
|
||
But numerical libraries aren't just about matrix formulas, and being
|
||
important doesn't necessarily mean taking up a lot of code: if matrix
|
||
formulas only occurred in one or two places in the average
|
||
numerically-oriented project, then it still wouldn't be worth adding a
|
||
new operator. So how common is matrix multiplication, really?
|
||
|
||
When the going gets tough, the tough get empirical. To get a rough
|
||
estimate of how useful the ``@`` operator will be, the table below
|
||
shows the rate at which different Python operators are actually used
|
||
in the stdlib, and also in two high-profile numerical packages -- the
|
||
scikit-learn machine learning library, and the nipy neuroimaging
|
||
library -- normalized by source lines of code (SLOC). Rows are sorted
|
||
by the 'combined' column, which pools all three code bases together.
|
||
The combined column is thus strongly weighted towards the stdlib,
|
||
which is much larger than both projects put together (stdlib: 411575
|
||
SLOC, scikit-learn: 50924 SLOC, nipy: 37078 SLOC). [#sloc-details]_
|
||
|
||
The ``dot`` row (marked ``******``) counts how common matrix multiply
|
||
operations are in each codebase.
|
||
|
||
::
|
||
|
||
==== ====== ============ ==== ========
|
||
op stdlib scikit-learn nipy combined
|
||
==== ====== ============ ==== ========
|
||
= 2969 5536 4932 3376 / 10,000 SLOC
|
||
- 218 444 496 261
|
||
+ 224 201 348 231
|
||
== 177 248 334 196
|
||
* 156 284 465 192
|
||
% 121 114 107 119
|
||
** 59 111 118 68
|
||
!= 40 56 74 44
|
||
/ 18 121 183 41
|
||
> 29 70 110 39
|
||
+= 34 61 67 39
|
||
< 32 62 76 38
|
||
>= 19 17 17 18
|
||
<= 18 27 12 18
|
||
dot ***** 0 ********** 99 ** 74 ****** 16
|
||
| 18 1 2 15
|
||
& 14 0 6 12
|
||
<< 10 1 1 8
|
||
// 9 9 1 8
|
||
-= 5 21 14 8
|
||
*= 2 19 22 5
|
||
/= 0 23 16 4
|
||
>> 4 0 0 3
|
||
^ 3 0 0 3
|
||
~ 2 4 5 2
|
||
|= 3 0 0 2
|
||
&= 1 0 0 1
|
||
//= 1 0 0 1
|
||
^= 1 0 0 0
|
||
**= 0 2 0 0
|
||
%= 0 0 0 0
|
||
<<= 0 0 0 0
|
||
>>= 0 0 0 0
|
||
==== ====== ============ ==== ========
|
||
|
||
These two numerical packages alone contain ~780 uses of matrix
|
||
multiplication. Within these packages, matrix multiplication is used
|
||
more heavily than most comparison operators (``<`` ``!=`` ``<=``
|
||
``>=``). Even when we dilute these counts by including the stdlib
|
||
into our comparisons, matrix multiplication is still used more often
|
||
in total than any of the bitwise operators, and 2x as often as ``//``.
|
||
This is true even though the stdlib, which contains a fair amount of
|
||
integer arithmetic and no matrix operations, makes up more than 80% of
|
||
the combined code base.
|
||
|
||
By coincidence, the numeric libraries make up approximately the same
|
||
proportion of the 'combined' codebase as numeric tutorials make up of
|
||
PyCon 2014's tutorial schedule, which suggests that the 'combined'
|
||
column may not be *wildly* unrepresentative of new Python code in
|
||
general. While it's impossible to know for certain, from this data it
|
||
seems entirely possible that across all Python code currently being
|
||
written, matrix multiplication is already used more often than ``//``
|
||
and the bitwise operations.
|
||
|
||
|
||
But isn't it weird to add an operator with no stdlib uses?
|
||
----------------------------------------------------------
|
||
|
||
It's certainly unusual (though extended slicing existed for some time
|
||
builtin types gained support for it, ``Ellipsis`` is still unused
|
||
within the stdlib, etc.). But the important thing is whether a change
|
||
will benefit users, not where the software is being downloaded from.
|
||
It's clear from the above that ``@`` will be used, and used heavily.
|
||
And this PEP provides the critical piece that will allow the Python
|
||
numerical community to finally reach consensus on a standard duck type
|
||
for all array-like objects, which is a necessary precondition to ever
|
||
adding a numerical array type to the stdlib.
|
||
|
||
|
||
Compatibility considerations
|
||
============================
|
||
|
||
Currently, the only legal use of the ``@`` token in Python code is at
|
||
statement beginning in decorators. The new operators are both infix;
|
||
the one place they can never occur is at statement beginning.
|
||
Therefore, no existing code will be broken by the addition of these
|
||
operators, and there is no possible parsing ambiguity between
|
||
decorator-@ and the new operators.
|
||
|
||
Another important kind of compatibility is the mental cost paid by
|
||
users to update their understanding of the Python language after this
|
||
change, particularly for users who do not work with matrices and thus
|
||
do not benefit. Here again, ``@`` has minimal impact: even
|
||
comprehensive tutorials and references will only need to add a
|
||
sentence or two to fully document this PEP's changes for a
|
||
non-numerical audience.
|
||
|
||
|
||
Intended usage details
|
||
======================
|
||
|
||
This section is informative, rather than normative -- it documents the
|
||
consensus of a number of libraries that provide array- or matrix-like
|
||
objects on how ``@`` will be implemented.
|
||
|
||
This section uses the numpy terminology for describing arbitrary
|
||
multidimensional arrays of data, because it is a superset of all other
|
||
commonly used models. In this model, the *shape* of any array is
|
||
represented by a tuple of integers. Because matrices are
|
||
two-dimensional, they have len(shape) == 2, while 1d vectors have
|
||
len(shape) == 1, and scalars have shape == (), i.e., they are "0
|
||
dimensional". Any array contains prod(shape) total entries. Notice
|
||
that `prod(()) == 1`_ (for the same reason that sum(()) == 0); scalars
|
||
are just an ordinary kind of array, not a special case. Notice also
|
||
that we distinguish between a single scalar value (shape == (),
|
||
analogous to ``1``), a vector containing only a single entry (shape ==
|
||
(1,), analogous to ``[1]``), a matrix containing only a single entry
|
||
(shape == (1, 1), analogous to ``[[1]]``), etc., so the dimensionality
|
||
of any array is always well-defined. Other libraries with more
|
||
restricted representations (e.g., those that support 2d arrays only)
|
||
might implement only a subset of the functionality described here.
|
||
|
||
.. _prod(()) == 1: https://en.wikipedia.org/wiki/Empty_product
|
||
|
||
Semantics
|
||
---------
|
||
|
||
The recommended semantics for ``@`` for different inputs are:
|
||
|
||
* 2d inputs are conventional matrices, and so the semantics are
|
||
obvious: we apply conventional matrix multiplication. If we write
|
||
``arr(2, 3)`` to represent an arbitrary 2x3 array, then ``arr(2, 3)
|
||
@ arr(3, 4)`` returns an array with shape (2, 4).
|
||
|
||
* 1d vector inputs are promoted to 2d by prepending or appending a '1'
|
||
to the shape, the operation is performed, and then the added
|
||
dimension is removed from the output. The 1 is always added on the
|
||
"outside" of the shape: prepended for left arguments, and appended
|
||
for right arguments. The result is that matrix @ vector and vector
|
||
@ matrix are both legal (assuming compatible shapes), and both
|
||
return 1d vectors; vector @ vector returns a scalar. This is
|
||
clearer with examples.
|
||
|
||
* ``arr(2, 3) @ arr(3, 1)`` is a regular matrix product, and returns
|
||
an array with shape (2, 1), i.e., a column vector.
|
||
|
||
* ``arr(2, 3) @ arr(3)`` performs the same computation as the
|
||
previous (i.e., treats the 1d vector as a matrix containing a
|
||
single *column*, shape = (3, 1)), but returns the result with
|
||
shape (2,), i.e., a 1d vector.
|
||
|
||
* ``arr(1, 3) @ arr(3, 2)`` is a regular matrix product, and returns
|
||
an array with shape (1, 2), i.e., a row vector.
|
||
|
||
* ``arr(3) @ arr(3, 2)`` performs the same computation as the
|
||
previous (i.e., treats the 1d vector as a matrix containing a
|
||
single *row*, shape = (1, 3)), but returns the result with shape
|
||
(2,), i.e., a 1d vector.
|
||
|
||
* ``arr(1, 3) @ arr(3, 1)`` is a regular matrix product, and returns
|
||
an array with shape (1, 1), i.e., a single value in matrix form.
|
||
|
||
* ``arr(3) @ arr(3)`` performs the same computation as the
|
||
previous, but returns the result with shape (), i.e., a single
|
||
scalar value, not in matrix form. So this is the standard inner
|
||
product on vectors.
|
||
|
||
An infelicity of this definition for 1d vectors is that it makes
|
||
``@`` non-associative in some cases (``(Mat1 @ vec) @ Mat2`` !=
|
||
``Mat1 @ (vec @ Mat2)``). But this seems to be a case where
|
||
practicality beats purity: non-associativity only arises for strange
|
||
expressions that would never be written in practice; if they are
|
||
written anyway then there is a consistent rule for understanding
|
||
what will happen (``Mat1 @ vec @ Mat2`` is parsed as ``(Mat1 @ vec)
|
||
@ Mat2``, just like ``a - b - c``); and, not supporting 1d vectors
|
||
would rule out many important use cases that do arise very commonly
|
||
in practice. No-one wants to explain to new users why to solve the
|
||
simplest linear system in the obvious way, they have to type
|
||
``(inv(A) @ b[:, np.newaxis]).flatten()`` instead of ``inv(A) @ b``,
|
||
or perform an ordinary least-squares regression by typing
|
||
``solve(X.T @ X, X @ y[:, np.newaxis]).flatten()`` instead of
|
||
``solve(X.T @ X, X @ y)``. No-one wants to type ``(a[np.newaxis, :]
|
||
@ b[:, np.newaxis])[0, 0]`` instead of ``a @ b`` every time they
|
||
compute an inner product, or ``(a[np.newaxis, :] @ Mat @ b[:,
|
||
np.newaxis])[0, 0]`` for general quadratic forms instead of ``a @
|
||
Mat @ b``. In addition, sage and sympy (see below) use these
|
||
non-associative semantics with an infix matrix multiplication
|
||
operator (they use ``*``), and they report that they haven't
|
||
experienced any problems caused by it.
|
||
|
||
* For inputs with more than 2 dimensions, we treat the last two
|
||
dimensions as being the dimensions of the matrices to multiply, and
|
||
'broadcast' across the other dimensions. This provides a convenient
|
||
way to quickly compute many matrix products in a single operation.
|
||
For example, ``arr(10, 2, 3) @ arr(10, 3, 4)`` performs 10 separate
|
||
matrix multiplies, each of which multiplies a 2x3 and a 3x4 matrix
|
||
to produce a 2x4 matrix, and then returns the 10 resulting matrices
|
||
together in an array with shape (10, 2, 4). The intuition here is
|
||
that we treat these 3d arrays of numbers as if they were 1d arrays
|
||
*of matrices*, and then apply matrix multiplication in an
|
||
elementwise manner, where now each 'element' is a whole matrix.
|
||
Note that broadcasting is not limited to perfectly aligned arrays;
|
||
in more complicated cases, it allows several simple but powerful
|
||
tricks for controlling how arrays are aligned with each other; see
|
||
[#broadcasting]_ for details. (In particular, it turns out that
|
||
when broadcasting is taken into account, the standard scalar *
|
||
matrix product is a special case of the elementwise multiplication
|
||
operator ``*``.)
|
||
|
||
If one operand is >2d, and another operand is 1d, then the above
|
||
rules apply unchanged, with 1d->2d promotion performed before
|
||
broadcasting. E.g., ``arr(10, 2, 3) @ arr(3)`` first promotes to
|
||
``arr(10, 2, 3) @ arr(3, 1)``, then broadcasts the right argument to
|
||
create the aligned operation ``arr(10, 2, 3) @ arr(10, 3, 1)``,
|
||
multiplies to get an array with shape (10, 2, 1), and finally
|
||
removes the added dimension, returning an array with shape (10, 2).
|
||
Similarly, ``arr(2) @ arr(10, 2, 3)`` produces an intermediate array
|
||
with shape (10, 1, 3), and a final array with shape (10, 3).
|
||
|
||
* 0d (scalar) inputs raise an error. Scalar * matrix multiplication
|
||
is a mathematically and algorithmically distinct operation from
|
||
matrix @ matrix multiplication, and is already covered by the
|
||
elementwise ``*`` operator. Allowing scalar @ matrix would thus
|
||
both require an unnecessary special case, and violate TOOWTDI.
|
||
|
||
|
||
Adoption
|
||
--------
|
||
|
||
We group existing Python projects which provide array- or matrix-like
|
||
types based on what API they currently use for elementwise and matrix
|
||
multiplication.
|
||
|
||
**Projects which currently use * for elementwise multiplication, and
|
||
function/method calls for matrix multiplication:**
|
||
|
||
The developers of the following projects have expressed an intention
|
||
to implement ``@`` on their array-like types using the above
|
||
semantics:
|
||
|
||
* numpy
|
||
* pandas
|
||
* blaze
|
||
* theano
|
||
|
||
The following projects have been alerted to the existence of the PEP,
|
||
but it's not yet known what they plan to do if it's accepted. We
|
||
don't anticipate that they'll have any objections, though, since
|
||
everything proposed here is consistent with how they already do
|
||
things:
|
||
|
||
* pycuda
|
||
* panda3d
|
||
|
||
**Projects which currently use * for matrix multiplication, and
|
||
function/method calls for elementwise multiplication:**
|
||
|
||
The following projects have expressed an intention, if this PEP is
|
||
accepted, to migrate from their current API to the elementwise-``*``,
|
||
matmul-``@`` convention (i.e., this is a list of projects whose API
|
||
fragmentation will probably be eliminated if this PEP is accepted):
|
||
|
||
* numpy (``numpy.matrix``)
|
||
* scipy.sparse
|
||
* pyoperators
|
||
* pyviennacl
|
||
|
||
The following projects have been alerted to the existence of the PEP,
|
||
but it's not known what they plan to do if it's accepted (i.e., this
|
||
is a list of projects whose API fragmentation may or may not be
|
||
eliminated if this PEP is accepted):
|
||
|
||
* cvxopt
|
||
|
||
**Projects which currently use * for matrix multiplication, and which
|
||
don't really care about elementwise multiplication of matrices:**
|
||
|
||
There are several projects which implement matrix types, but from a
|
||
very different perspective than the numerical libraries discussed
|
||
above. These projects focus on computational methods for analyzing
|
||
matrices in the sense of abstract mathematical objects (i.e., linear
|
||
maps over free modules over rings), rather than as big bags full of
|
||
numbers that need crunching. And it turns out that from the abstract
|
||
math point of view, there isn't much use for elementwise operations in
|
||
the first place; as discussed in the Background section above,
|
||
elementwise operations are motivated by the bag-of-numbers approach.
|
||
So these projects don't encounter the basic problem that this PEP
|
||
exists to address, making it mostly irrelevant to them; while they
|
||
appear superficially similar to projects like numpy, they're actually
|
||
doing something quite different. They use ``*`` for matrix
|
||
multiplication (and for group actions, and so forth), and if this PEP
|
||
is accepted, their expressed intention is to continue doing so, while
|
||
perhaps adding ``@`` as an alias. These projects include:
|
||
|
||
* sympy
|
||
* sage
|
||
|
||
|
||
Implementation details
|
||
======================
|
||
|
||
New functions ``operator.matmul`` and ``operator.__matmul__`` are
|
||
added to the standard library, with the usual semantics.
|
||
|
||
A corresponding function ``PyObject* PyObject_MatrixMultiply(PyObject
|
||
*o1, PyObject *o2)`` is added to the C API.
|
||
|
||
A new AST node is added named ``MatMult``, along with a new token
|
||
``ATEQUAL`` and new bytecode opcodes ``BINARY_MATRIX_MULTIPLY`` and
|
||
``INPLACE_MATRIX_MULTIPLY``.
|
||
|
||
Two new type slots are added; whether this is to ``PyNumberMethods``
|
||
or a new ``PyMatrixMethods`` struct remains to be determined.
|
||
|
||
|
||
Rationale for specification details
|
||
===================================
|
||
|
||
Choice of operator
|
||
------------------
|
||
|
||
Why ``@`` instead of some other spelling? There isn't any consensus
|
||
across other programming languages about how this operator should be
|
||
named [#matmul-other-langs]_; here we discuss the various options.
|
||
|
||
Restricting ourselves only to symbols present on US English keyboards,
|
||
the punctuation characters that don't already have a meaning in Python
|
||
expression context are: ``@``, backtick, ``$``, ``!``, and ``?``. Of
|
||
these options, ``@`` is clearly the best; ``!`` and ``?`` are already
|
||
heavily freighted with inapplicable meanings in the programming
|
||
context, backtick has been banned from Python by BDFL pronouncement
|
||
(see PEP 3099), and ``$`` is uglier, even more dissimilar to ``*`` and
|
||
:math:`\cdot`, and has Perl/PHP baggage. ``$`` is probably the
|
||
second-best option of these, though.
|
||
|
||
Symbols which are not present on US English keyboards start at a
|
||
significant disadvantage (having to spend 5 minutes at the beginning
|
||
of every numeric Python tutorial just going over keyboard layouts is
|
||
not a hassle anyone really wants). Plus, even if we somehow overcame
|
||
the typing problem, it's not clear there are any that are actually
|
||
better than ``@``. Some options that have been suggested include:
|
||
|
||
* U+00D7 MULTIPLICATION SIGN: ``A × B``
|
||
* U+22C5 DOT OPERATOR: ``A ⋅ B``
|
||
* U+2297 CIRCLED TIMES: ``A ⊗ B``
|
||
* U+00B0 DEGREE: ``A ° B``
|
||
|
||
What we need, though, is an operator that means "matrix
|
||
multiplication, as opposed to scalar/elementwise multiplication".
|
||
There is no conventional symbol with this meaning in either
|
||
programming or mathematics, where these operations are usually
|
||
distinguished by context. (And U+2297 CIRCLED TIMES is actually used
|
||
conventionally to mean exactly the wrong things: elementwise
|
||
multiplication -- the "Hadamard product" -- or outer product, rather
|
||
than matrix/inner product like our operator). ``@`` at least has the
|
||
virtue that it *looks* like a funny non-commutative operator; a naive
|
||
user who knows maths but not programming couldn't look at ``A * B``
|
||
versus ``A × B``, or ``A * B`` versus ``A ⋅ B``, or ``A * B`` versus
|
||
``A ° B`` and guess which one is the usual multiplication, and which
|
||
one is the special case.
|
||
|
||
Finally, there is the option of using multi-character tokens. Some
|
||
options:
|
||
|
||
* Matlab and Julia use a ``.*`` operator. Aside from being visually
|
||
confusable with ``*``, this would be a terrible choice for us
|
||
because in Matlab and Julia, ``*`` means matrix multiplication and
|
||
``.*`` means elementwise multiplication, so using ``.*`` for matrix
|
||
multiplication would make us exactly backwards from what Matlab and
|
||
Julia users expect.
|
||
|
||
* APL apparently used ``+.×``, which by combining a multi-character
|
||
token, confusing attribute-access-like . syntax, and a unicode
|
||
character, ranks somewhere below U+2603 SNOWMAN on our candidate
|
||
list. If we like the idea of combining addition and multiplication
|
||
operators as being evocative of how matrix multiplication actually
|
||
works, then something like ``+*`` could be used -- though this may
|
||
be too easy to confuse with ``*+``, which is just multiplication
|
||
combined with the unary ``+`` operator.
|
||
|
||
* PEP 211 suggested ``~*``. This has the downside that it sort of
|
||
suggests that there is a unary ``*`` operator that is being combined
|
||
with unary ``~``, but it could work.
|
||
|
||
* R uses ``%*%`` for matrix multiplication. In R this forms part of a
|
||
general extensible infix system in which all tokens of the form
|
||
``%foo%`` are user-defined binary operators. We could steal the
|
||
token without stealing the system.
|
||
|
||
* Some other plausible candidates that have been suggested: ``><`` (=
|
||
ascii drawing of the multiplication sign ×); the footnote operator
|
||
``[*]`` or ``|*|`` (but when used in context, the use of vertical
|
||
grouping symbols tends to recreate the nested parentheses visual
|
||
clutter that was noted as one of the major downsides of the function
|
||
syntax we're trying to get away from); ``^*``.
|
||
|
||
So, it doesn't matter much, but ``@`` seems as good or better than any
|
||
of the alternatives:
|
||
|
||
* It's a friendly character that Pythoneers are already used to typing
|
||
in decorators, but the decorator usage and the math expression
|
||
usage are sufficiently dissimilar that it would be hard to confuse
|
||
them in practice.
|
||
|
||
* It's widely accessible across keyboard layouts (and thanks to its
|
||
use in email addresses, this is true even of weird keyboards like
|
||
those in phones).
|
||
|
||
* It's round like ``*`` and :math:`\cdot`.
|
||
|
||
* The mATrices mnemonic is cute.
|
||
|
||
* The swirly shape is reminiscent of the simultaneous sweeps over rows
|
||
and columns that define matrix multiplication
|
||
|
||
* Its asymmetry is evocative of its non-commutative nature.
|
||
|
||
* Whatever, we have to pick something.
|
||
|
||
|
||
Precedence and associativity
|
||
----------------------------
|
||
|
||
There was a long discussion [#associativity-discussions]_ about
|
||
whether ``@`` should be right- or left-associative (or even something
|
||
more exotic [#group-associativity]_). Almost all Python operators are
|
||
left-associative, so following this convention would be the simplest
|
||
approach, but there were two arguments that suggested matrix
|
||
multiplication might be worth making right-associative as a special
|
||
case:
|
||
|
||
First, matrix multiplication has a tight conceptual association with
|
||
function application/composition, so many mathematically sophisticated
|
||
users have an intuition that an expression like :math:`R S x` proceeds
|
||
from right-to-left, with first :math:`S` transforming the vector
|
||
:math:`x`, and then :math:`R` transforming the result. This isn't
|
||
universally agreed (and not all number-crunchers are steeped in the
|
||
pure-math conceptual framework that motivates this intuition
|
||
[#oil-industry-versus-right-associativity]_), but at the least this
|
||
intuition is more common than for other operations like :math:`2 \cdot
|
||
3 \cdot 4` which everyone reads as going from left-to-right.
|
||
|
||
Second, if expressions like ``Mat @ Mat @ vec`` appear often in code,
|
||
then programs will run faster (and efficiency-minded programmers will
|
||
be able to use fewer parentheses) if this is evaluated as ``Mat @ (Mat
|
||
@ vec)`` then if it is evaluated like ``(Mat @ Mat) @ vec``.
|
||
|
||
However, weighing against these arguments are the following:
|
||
|
||
Regarding the efficiency argument, empirically, we were unable to find
|
||
any evidence that ``Mat @ Mat @ vec`` type expressions actually
|
||
dominate in real-life code. Parsing a number of large projects that
|
||
use numpy, we found that when forced by numpy's current funcall syntax
|
||
to choose an order of operations for nested calls to ``dot``, people
|
||
actually use left-associative nesting slightly *more* often than
|
||
right-associative nesting [#numpy-associativity-counts]_. And anyway,
|
||
writing parentheses isn't so bad -- if an efficiency-minded programmer
|
||
is going to take the trouble to think through the best way to evaluate
|
||
some expression, they probably *should* write down the parentheses
|
||
regardless of whether they're needed, just to make it obvious to the
|
||
next reader that they order of operations matter.
|
||
|
||
In addition, it turns out that other languages, including those with
|
||
much more of a focus on linear algebra, overwhelmingly make their
|
||
matmul operators left-associative. Specifically, the ``@`` equivalent
|
||
is left-associative in R, Matlab, Julia, IDL, and Gauss. The only
|
||
exceptions we found are Mathematica, in which ``a @ b @ c`` would be
|
||
parsed non-associatively as ``dot(a, b, c)``, and APL, in which all
|
||
operators are right-associative. There do not seem to exist any
|
||
languages that make ``@`` right-associative and ``*``
|
||
left-associative. And these decisions don't seem to be controversial
|
||
-- I've never seen anyone complaining about this particular aspect of
|
||
any of these other languages, and the left-associativity of ``*``
|
||
doesn't seem to bother users of the existing Python libraries that use
|
||
``*`` for matrix multiplication. So, at the least we can conclude from
|
||
this that making ``@`` left-associative will certainly not cause any
|
||
disasters. Making ``@`` right-associative, OTOH, would be exploring
|
||
new and uncertain ground.
|
||
|
||
And another advantage of left-associativity is that it is much easier
|
||
to learn and remember that ``@`` acts like ``*``, than it is to
|
||
remember first that ``@`` is unlike other Python operators by being
|
||
right-associative, and then on top of this, also have to remember
|
||
whether it is more tightly or more loosely binding than
|
||
``*``. (Right-associativity forces us to choose a precedence, and
|
||
intuitions were about equally split on which precedence made more
|
||
sense. So this suggests that no matter which choice we made, no-one
|
||
would be able to guess or remember it.)
|
||
|
||
On net, therefore, the general consensus of the numerical community is
|
||
that while matrix multiplication is something of a special case, it's
|
||
not special enough to break the rules, and ``@`` should parse like
|
||
``*`` does.
|
||
|
||
|
||
(Non)-Definitions for built-in types
|
||
------------------------------------
|
||
|
||
No ``__matmul__`` or ``__matpow__`` are defined for builtin numeric
|
||
types (``float``, ``int``, etc.) or for the ``numbers.Number``
|
||
hierarchy, because these types represent scalars, and the consensus
|
||
semantics for ``@`` are that it should raise an error on scalars.
|
||
|
||
We do not -- for now -- define a ``__matmul__`` method on the standard
|
||
``memoryview`` or ``array.array`` objects, for several reasons. Of
|
||
course this could be added if someone wants it, but these types would
|
||
require quite a bit of additional work beyond ``__matmul__`` before
|
||
they could be used for numeric work -- e.g., they have no way to do
|
||
addition or scalar multiplication either! -- and adding such
|
||
functionality is beyond the scope of this PEP. In addition, providing
|
||
a quality implementation of matrix multiplication is highly
|
||
non-trivial. Naive nested loop implementations are very slow and
|
||
shipping such an implementation in CPython would just create a trap
|
||
for users. But the alternative -- providing a modern, competitive
|
||
matrix multiply -- would require that CPython link to a BLAS library,
|
||
which brings a set of new complications. In particular, several
|
||
popular BLAS libraries (including the one that ships by default on
|
||
OS X) currently break the use of ``multiprocessing`` [#blas-fork]_.
|
||
Together, these considerations mean that the cost/benefit of adding
|
||
``__matmul__`` to these types just isn't there, so for now we'll
|
||
continue to delegate these problems to numpy and friends, and defer a
|
||
more systematic solution to a future proposal.
|
||
|
||
There are also non-numeric Python builtins which define ``__mul__``
|
||
(``str``, ``list``, ...). We do not define ``__matmul__`` for these
|
||
types either, because why would we even do that.
|
||
|
||
|
||
Non-definition of matrix power
|
||
------------------------------
|
||
|
||
Earlier versions of this PEP also proposed a matrix power operator,
|
||
``@@``, analogous to ``**``. But on further consideration, it was
|
||
decided that the utility of this was sufficiently unclear that it
|
||
would be better to leave it out for now, and only revisit the issue if
|
||
-- once we have more experience with ``@`` -- it turns out that ``@@``
|
||
is truly missed. [#atat-discussion]_
|
||
|
||
|
||
Rejected alternatives to adding a new operator
|
||
==============================================
|
||
|
||
Over the past few decades, the Python numeric community has explored a
|
||
variety of ways to resolve the tension between matrix and elementwise
|
||
multiplication operations. PEP 211 and PEP 225, both proposed in 2000
|
||
and last seriously discussed in 2008 [#threads-2008]_, were early
|
||
attempts to add new operators to solve this problem, but suffered from
|
||
serious flaws; in particular, at that time the Python numerical
|
||
community had not yet reached consensus on the proper API for array
|
||
objects, or on what operators might be needed or useful (e.g., PEP 225
|
||
proposes 6 new operators with unspecified semantics). Experience
|
||
since then has now led to consensus that the best solution, for both
|
||
numeric Python and core Python, is to add a single infix operator for
|
||
matrix multiply (together with the other new operators this implies
|
||
like ``@=``).
|
||
|
||
We review some of the rejected alternatives here.
|
||
|
||
**Use a second type that defines __mul__ as matrix multiplication:**
|
||
As discussed above (`Background: What's wrong with the status quo?`_),
|
||
this has been tried this for many years via the ``numpy.matrix`` type
|
||
(and its predecessors in Numeric and numarray). The result is a
|
||
strong consensus among both numpy developers and developers of
|
||
downstream packages that ``numpy.matrix`` should essentially never be
|
||
used, because of the problems caused by having conflicting duck types
|
||
for arrays. (Of course one could then argue we should *only* define
|
||
``__mul__`` to be matrix multiplication, but then we'd have the same
|
||
problem with elementwise multiplication.) There have been several
|
||
pushes to remove ``numpy.matrix`` entirely; the only counter-arguments
|
||
have come from educators who find that its problems are outweighed by
|
||
the need to provide a simple and clear mapping between mathematical
|
||
notation and code for novices (see `Transparent syntax is especially
|
||
crucial for non-expert programmers`_). But, of course, starting out
|
||
newbies with a dispreferred syntax and then expecting them to
|
||
transition later causes its own problems. The two-type solution is
|
||
worse than the disease.
|
||
|
||
**Add lots of new operators, or add a new generic syntax for defining
|
||
infix operators:** In addition to being generally un-Pythonic and
|
||
repeatedly rejected by BDFL fiat, this would be using a sledgehammer
|
||
to smash a fly. The scientific python community has consensus that
|
||
adding one operator for matrix multiplication is enough to fix the one
|
||
otherwise unfixable pain point. (In retrospect, we all think PEP 225
|
||
was a bad idea too -- or at least far more complex than it needed to
|
||
be.)
|
||
|
||
**Add a new @ (or whatever) operator that has some other meaning in
|
||
general Python, and then overload it in numeric code:** This was the
|
||
approach taken by PEP 211, which proposed defining ``@`` to be the
|
||
equivalent of ``itertools.product``. The problem with this is that
|
||
when taken on its own terms, it's pretty clear that
|
||
``itertools.product`` doesn't actually need a dedicated operator. It
|
||
hasn't even been deemed worth of a builtin. (During discussions of
|
||
this PEP, a similar suggestion was made to define ``@`` as a general
|
||
purpose function composition operator, and this suffers from the same
|
||
problem; ``functools.compose`` isn't even useful enough to exist.)
|
||
Matrix multiplication has a uniquely strong rationale for inclusion as
|
||
an infix operator. There almost certainly don't exist any other
|
||
binary operations that will ever justify adding any other infix
|
||
operators to Python.
|
||
|
||
**Add a .dot method to array types so as to allow "pseudo-infix"
|
||
A.dot(B) syntax:** This has been in numpy for some years, and in many
|
||
cases it's better than dot(A, B). But it's still much less readable
|
||
than real infix notation, and in particular still suffers from an
|
||
extreme overabundance of parentheses. See `Why should matrix
|
||
multiplication be infix?`_ above.
|
||
|
||
**Use a 'with' block to toggle the meaning of * within a single code
|
||
block**: E.g., numpy could define a special context object so that
|
||
we'd have::
|
||
|
||
c = a * b # element-wise multiplication
|
||
with numpy.mul_as_dot:
|
||
c = a * b # matrix multiplication
|
||
|
||
However, this has two serious problems: first, it requires that every
|
||
array-like type's ``__mul__`` method know how to check some global
|
||
state (``numpy.mul_is_currently_dot`` or whatever). This is fine if
|
||
``a`` and ``b`` are numpy objects, but the world contains many
|
||
non-numpy array-like objects. So this either requires non-local
|
||
coupling -- every numpy competitor library has to import numpy and
|
||
then check ``numpy.mul_is_currently_dot`` on every operation -- or
|
||
else it breaks duck-typing, with the above code doing radically
|
||
different things depending on whether ``a`` and ``b`` are numpy
|
||
objects or some other sort of object. Second, and worse, ``with``
|
||
blocks are dynamically scoped, not lexically scoped; i.e., any
|
||
function that gets called inside the ``with`` block will suddenly find
|
||
itself executing inside the mul_as_dot world, and crash and burn
|
||
horribly -- if you're lucky. So this is a construct that could only
|
||
be used safely in rather limited cases (no function calls), and which
|
||
would make it very easy to shoot yourself in the foot without warning.
|
||
|
||
**Use a language preprocessor that adds extra numerically-oriented
|
||
operators and perhaps other syntax:** (As per recent BDFL suggestion:
|
||
[#preprocessor]_) This suggestion seems based on the idea that
|
||
numerical code needs a wide variety of syntax additions. In fact,
|
||
given ``@``, most numerical users don't need any other operators or
|
||
syntax; it solves the one really painful problem that cannot be solved
|
||
by other means, and that causes painful reverberations through the
|
||
larger ecosystem. Defining a new language (presumably with its own
|
||
parser which would have to be kept in sync with Python's, etc.), just
|
||
to support a single binary operator, is neither practical nor
|
||
desirable. In the numerical context, Python's competition is
|
||
special-purpose numerical languages (Matlab, R, IDL, etc.). Compared
|
||
to these, Python's killer feature is exactly that one can mix
|
||
specialized numerical code with code for XML parsing, web page
|
||
generation, database access, network programming, GUI libraries, and
|
||
so forth, and we also gain major benefits from the huge variety of
|
||
tutorials, reference material, introductory classes, etc., which use
|
||
Python. Fragmenting "numerical Python" from "real Python" would be a
|
||
major source of confusion. A major motivation for this PEP is to
|
||
*reduce* fragmentation. Having to set up a preprocessor would be an
|
||
especially prohibitive complication for unsophisticated users. And we
|
||
use Python because we like Python! We don't want
|
||
almost-but-not-quite-Python.
|
||
|
||
**Use overloading hacks to define a "new infix operator" like *dot*,
|
||
as in a well-known Python recipe:** (See: [#infix-hack]_) Beautiful is
|
||
better than ugly. This is... not beautiful. And not Pythonic. And
|
||
especially unfriendly to beginners, who are just trying to wrap their
|
||
heads around the idea that there's a coherent underlying system behind
|
||
these magic incantations that they're learning, when along comes an
|
||
evil hack like this that violates that system, creates bizarre error
|
||
messages when accidentally misused, and whose underlying mechanisms
|
||
can't be understood without deep knowledge of how object oriented
|
||
systems work.
|
||
|
||
**Use a special "facade" type to support syntax like arr.M * arr:**
|
||
This is very similar to the previous proposal, in that the ``.M``
|
||
attribute would basically return the same object as ``arr *dot` would,
|
||
and thus suffers the same objections about 'magicalness'. This
|
||
approach also has some non-obvious complexities: for example, while
|
||
``arr.M * arr`` must return an array, ``arr.M * arr.M`` and ``arr *
|
||
arr.M`` must return facade objects, or else ``arr.M * arr.M * arr``
|
||
and ``arr * arr.M * arr`` will not work. But this means that facade
|
||
objects must be able to recognize both other array objects and other
|
||
facade objects (which creates additional complexity for writing
|
||
interoperating array types from different libraries who must now
|
||
recognize both each other's array types and their facade types). It
|
||
also creates pitfalls for users who may easily type ``arr * arr.M`` or
|
||
``arr.M * arr.M`` and expect to get back an array object; instead,
|
||
they will get a mysterious object that throws errors when they attempt
|
||
to use it. Basically with this approach users must be careful to
|
||
think of ``.M*`` as an indivisible unit that acts as an infix operator
|
||
-- and as infix-operator-like token strings go, at least ``*dot*``
|
||
is prettier looking (look at its cute little ears!).
|
||
|
||
|
||
Discussions of this PEP
|
||
=======================
|
||
|
||
Collected here for reference:
|
||
|
||
* Github pull request containing much of the original discussion and
|
||
drafting: https://github.com/numpy/numpy/pull/4351
|
||
|
||
* sympy mailing list discussions of an early draft:
|
||
|
||
* https://groups.google.com/forum/#!topic/sympy/22w9ONLa7qo
|
||
* https://groups.google.com/forum/#!topic/sympy/4tGlBGTggZY
|
||
|
||
* sage-devel mailing list discussions of an early draft:
|
||
https://groups.google.com/forum/#!topic/sage-devel/YxEktGu8DeM
|
||
|
||
* 13-Mar-2014 python-ideas thread:
|
||
https://mail.python.org/pipermail/python-ideas/2014-March/027053.html
|
||
|
||
* numpy-discussion thread on whether to keep ``@@``:
|
||
http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069448.html
|
||
|
||
* numpy-discussion threads on precedence/associativity of ``@``:
|
||
* http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069444.html
|
||
* http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069605.html
|
||
|
||
|
||
References
|
||
==========
|
||
|
||
.. [#preprocessor] From a comment by GvR on a G+ post by GvR; the
|
||
comment itself does not seem to be directly linkable: https://plus.google.com/115212051037621986145/posts/hZVVtJ9bK3u
|
||
.. [#infix-hack] http://code.activestate.com/recipes/384122-infix-operators/
|
||
http://www.sagemath.org/doc/reference/misc/sage/misc/decorators.html#sage.misc.decorators.infix_operator
|
||
.. [#scipy-conf] http://conference.scipy.org/past.html
|
||
.. [#pydata-conf] http://pydata.org/events/
|
||
.. [#lht] In this formula, :math:`\beta` is a vector or matrix of
|
||
regression coefficients, :math:`V` is the estimated
|
||
variance/covariance matrix for these coefficients, and we want to
|
||
test the null hypothesis that :math:`H\beta = r`; a large :math:`S`
|
||
then indicates that this hypothesis is unlikely to be true. For
|
||
example, in an analysis of human height, the vector :math:`\beta`
|
||
might contain one value which was the average height of the
|
||
measured men, and another value which was the average height of the
|
||
measured women, and then setting :math:`H = [1, -1], r = 0` would
|
||
let us test whether men and women are the same height on
|
||
average. Compare to eq. 2.139 in
|
||
http://sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore/tutorials/xegbohtmlnode17.html
|
||
|
||
Example code is adapted from https://github.com/rerpy/rerpy/blob/0d274f85e14c3b1625acb22aed1efa85d122ecb7/rerpy/incremental_ls.py#L202
|
||
|
||
.. [#pycon-tutorials] Out of the 36 tutorials scheduled for PyCon 2014
|
||
(https://us.pycon.org/2014/schedule/tutorials/), we guess that the
|
||
8 below will almost certainly deal with matrices:
|
||
|
||
* Dynamics and control with Python
|
||
|
||
* Exploring machine learning with Scikit-learn
|
||
|
||
* How to formulate a (science) problem and analyze it using Python
|
||
code
|
||
|
||
* Diving deeper into Machine Learning with Scikit-learn
|
||
|
||
* Data Wrangling for Kaggle Data Science Competitions – An etude
|
||
|
||
* Hands-on with Pydata: how to build a minimal recommendation
|
||
engine.
|
||
|
||
* Python for Social Scientists
|
||
|
||
* Bayesian statistics made simple
|
||
|
||
In addition, the following tutorials could easily involve matrices:
|
||
|
||
* Introduction to game programming
|
||
|
||
* mrjob: Snakes on a Hadoop *("We'll introduce some data science
|
||
concepts, such as user-user similarity, and show how to calculate
|
||
these metrics...")*
|
||
|
||
* Mining Social Web APIs with IPython Notebook
|
||
|
||
* Beyond Defaults: Creating Polished Visualizations Using Matplotlib
|
||
|
||
This gives an estimated range of 8 to 12 / 36 = 22% to 33% of
|
||
tutorials dealing with matrices; saying ~20% then gives us some
|
||
wiggle room in case our estimates are high.
|
||
|
||
.. [#sloc-details] SLOCs were defined as physical lines which contain
|
||
at least one token that is not a COMMENT, NEWLINE, ENCODING,
|
||
INDENT, or DEDENT. Counts were made by using ``tokenize`` module
|
||
from Python 3.2.3 to examine the tokens in all files ending ``.py``
|
||
underneath some directory. Only tokens which occur at least once
|
||
in the source trees are included in the table. The counting script
|
||
is available `in the PEP repository
|
||
<http://hg.python.org/peps/file/tip/pep-0465/scan-ops.py>`_.
|
||
|
||
Matrix multiply counts were estimated by counting how often certain
|
||
tokens which are used as matrix multiply function names occurred in
|
||
each package. This creates a small number of false positives for
|
||
scikit-learn, because we also count instances of the wrappers
|
||
around ``dot`` that this package uses, and so there are a few dozen
|
||
tokens which actually occur in ``import`` or ``def`` statements.
|
||
|
||
All counts were made using the latest development version of each
|
||
project as of 21 Feb 2014.
|
||
|
||
'stdlib' is the contents of the Lib/ directory in commit
|
||
d6aa3fa646e2 to the cpython hg repository, and treats the following
|
||
tokens as indicating matrix multiply: n/a.
|
||
|
||
'scikit-learn' is the contents of the sklearn/ directory in commit
|
||
69b71623273ccfc1181ea83d8fb9e05ae96f57c7 to the scikit-learn
|
||
repository (https://github.com/scikit-learn/scikit-learn), and
|
||
treats the following tokens as indicating matrix multiply: ``dot``,
|
||
``fast_dot``, ``safe_sparse_dot``.
|
||
|
||
'nipy' is the contents of the nipy/ directory in commit
|
||
5419911e99546401b5a13bd8ccc3ad97f0d31037 to the nipy repository
|
||
(https://github.com/nipy/nipy/), and treats the following tokens as
|
||
indicating matrix multiply: ``dot``.
|
||
|
||
.. [#blas-fork] BLAS libraries have a habit of secretly spawning
|
||
threads, even when used from single-threaded programs. And threads
|
||
play very poorly with ``fork()``; the usual symptom is that
|
||
attempting to perform linear algebra in a child process causes an
|
||
immediate deadlock.
|
||
|
||
.. [#threads-2008] http://fperez.org/py4science/numpy-pep225/numpy-pep225.html
|
||
|
||
.. [#broadcasting] http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html
|
||
|
||
.. [#matmul-other-langs] http://mail.scipy.org/pipermail/scipy-user/2014-February/035499.html
|
||
|
||
.. [#github-details] Counts were produced by manually entering the
|
||
string ``"import foo"`` or ``"from foo import"`` (with quotes) into
|
||
the Github code search page, e.g.:
|
||
https://github.com/search?q=%22import+numpy%22&ref=simplesearch&type=Code
|
||
on 2014-04-10 at ~21:00 UTC. The reported values are the numbers
|
||
given in the "Languages" box on the lower-left corner, next to
|
||
"Python". This also causes some undercounting (e.g., leaving out
|
||
Cython code, and possibly one should also count HTML docs and so
|
||
forth), but these effects are negligible (e.g., only ~1% of numpy
|
||
usage appears to occur in Cython code, and probably even less for
|
||
the other modules listed). The use of this box is crucial,
|
||
however, because these counts appear to be stable, while the
|
||
"overall" counts listed at the top of the page ("We've found ___
|
||
code results") are highly variable even for a single search --
|
||
simply reloading the page can cause this number to vary by a factor
|
||
of 2 (!!). (They do seem to settle down if one reloads the page
|
||
repeatedly, but nonetheless this is spooky enough that it seemed
|
||
better to avoid these numbers.)
|
||
|
||
These numbers should of course be taken with multiple grains of
|
||
salt; it's not clear how representative Github is of Python code in
|
||
general, and limitations of the search tool make it impossible to
|
||
get precise counts. AFAIK this is the best data set currently
|
||
available, but it'd be nice if it were better. In particular:
|
||
|
||
* Lines like ``import sys, os`` will only be counted in the ``sys``
|
||
row.
|
||
|
||
* A file containing both ``import X`` and ``from X import`` will be
|
||
counted twice
|
||
|
||
* Imports of the form ``from X.foo import ...`` are missed. We
|
||
could catch these by instead searching for "from X", but this is
|
||
a common phrase in English prose, so we'd end up with false
|
||
positives from comments, strings, etc. For many of the modules
|
||
considered this shouldn't matter too much -- for example, the
|
||
stdlib modules have flat namespaces -- but it might especially
|
||
lead to undercounting of django, scipy, and twisted.
|
||
|
||
Also, it's possible there exist other non-stdlib modules we didn't
|
||
think to test that are even more-imported than numpy -- though we
|
||
tried quite a few of the obvious suspects. If you find one, let us
|
||
know! The modules tested here were chosen based on a combination
|
||
of intuition and the top-100 list at pypi-ranking.info.
|
||
|
||
Fortunately, it doesn't really matter if it turns out that numpy
|
||
is, say, merely the *third* most-imported non-stdlib module, since
|
||
the point is just that numeric programming is a common and
|
||
mainstream activity.
|
||
|
||
Finally, we should point out the obvious: whether a package is
|
||
import**ed** is rather different from whether it's import**ant**.
|
||
No-one's claiming numpy is "the most important package" or anything
|
||
like that. Certainly more packages depend on distutils, e.g., then
|
||
depend on numpy -- and far fewer source files import distutils than
|
||
import numpy. But this is fine for our present purposes. Most
|
||
source files don't import distutils because most source files don't
|
||
care how they're distributed, so long as they are; these source
|
||
files thus don't care about details of how distutils' API works.
|
||
This PEP is in some sense about changing how numpy's and related
|
||
packages' APIs work, so the relevant metric is to look at source
|
||
files that are choosing to directly interact with that API, which
|
||
is sort of like what we get by looking at import statements.
|
||
|
||
.. [#hugunin] The first such proposal occurs in Jim Hugunin's very
|
||
first email to the matrix SIG in 1995, which lays out the first
|
||
draft of what became Numeric. He suggests using ``*`` for
|
||
elementwise multiplication, and ``%`` for matrix multiplication:
|
||
https://mail.python.org/pipermail/matrix-sig/1995-August/000002.html
|
||
|
||
.. [#atat-discussion] http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069502.html
|
||
|
||
.. [#associativity-discussions]
|
||
http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069444.html
|
||
http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069605.html
|
||
|
||
.. [#oil-industry-versus-right-associativity]
|
||
http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069610.html
|
||
|
||
.. [#numpy-associativity-counts]
|
||
http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069578.html
|
||
|
||
.. [#group-associativity]
|
||
http://mail.scipy.org/pipermail/numpy-discussion/2014-March/069530.html
|
||
|
||
|
||
Copyright
|
||
=========
|
||
|
||
This document has been placed in the public domain.
|