ENH: Implement matmul function

[abalkin]

The @-operator PEP (gh-4351) describes two operations that would be a valuable addition to numpy independently of the main purpose of the PEP.

Most of the functionality is already available in various places (np.dot, np.linalg.matrix_power, etc.) but no single function conforms to the PEP.

[njsmith]

The @@ in the current draft PEP is the same as np.linalg.matrix_power,
except that np.linalg.matrix_power currently does not handle >2d inputs --
it needs to be gufuncified (or just fixed by hand).

The @ operator in the current PEP draft is what we want np.dot to be
eventually, regardless of the PEP's fate, so it would indeed be good to get
an implementation into numpy somewhere -- we'll need something to direct
people to during the np.dot deprecation cycle.

On Sun, Mar 9, 2014 at 3:53 AM, abalkin [email protected] wrote:

The @-operator PEP (gh-4351 #4351)
describes two operations that would be a valuable addition to numpy
independently of the main purpose of the PEP.

Most of the functionality is already available in various places (np.dot,
np.linalg.matrix_power, etc.) but no single function conforms to the PEP.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/4464
.

Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org

[pv]

Related: gh-4397

[abalkin]

[@asmeurer]

Re boolean arrays: are Numpy booleans different from Python booleans?
In Python. True and False act exactly like 1 and 0 except for their printing.

They are very different:

from numpy import *
bool_(1) + bool_(1)
True

See also #4105.

[asmeurer]

I see. That's really a good thing, because the logical operations ~, &, |, and so on act on booleans, not bitwise (e.g., in Python ~True gives something that has a boolean value of True). We do the same thing in SymPy, see http://docs.sympy.org/latest/modules/logic.html#sympy.logic.boolalg.BooleanTrue (except our BooleanTrue doesn't allow addition; we would use a distinct FF(2) object for that).

[abalkin]

@asmeurer - so what should matmul do for two boolean matrices?

[asmeurer]

Just the regular algorithm. It looks like NumPy considers booleans to be the elements of the finite field Z_2, so true + true = true, true*true = true, and so on.

[abalkin]

It's more like Z_2 semi-ring. In the Z_2 field 1 + 1 = 0.

[asmeurer]

Oh you're right. I didn't notice that inconsistency. I that case, I have no idea why they did that, and really no idea what matmul should do, other than just give whatever you would get with the normal matrix multiplication formula.

[abalkin]

Since final PEP 465 does not include matpow, I am changing the title of this issue to focus on matmul.

[abalkin]

Here is the preliminary list of classes that should implement __matmul__:

1. numpy.ndarray (C)
2. numpy.matrixlib.defmatrix.matrix (Python)
3. numpy.ma.MaskedArray (Python)

Any additions?

[njsmith]

The ideal implementation I think would be:

1. add support for "optional axes" in gufunc signatures, as per this
   proposal:
   http://mail.scipy.org/pipermail/numpy-discussion/2013-July/067217.html

2. use this to make a matrix multiply gufunc which just has the PEP dot
   semantics directly, without needing messy wrapper stuff.

3. wire up nb_matmul to call this gufunc.

This would involve messing with two potentially messy pieces of code
(gufunc signature interpretation, and the special dtype-based dot
dispatch), but it's probably worth wading in and attempting the proper
solution... Maybe it will turn out to not be that bad.
On 8 Apr 2014 18:27, "abalkin" [email protected] wrote:

Here is the preliminary list of classes that should implement matmul:

1. numpy.ndarray (C)
2. numpy.matrixlib.defmatrix.matrix (Python)
3. numpy.ma.MaskedArray (Python)

Any additions?

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/4464#issuecomment-39877055
.

[nicktimko]

Python 3.5.0a3 has added the @ binary operator to the syntax, is there some NumPy branch available to test with? Naively trying to subclass from np.array fails because np.ndarray is actually the type and it gets really messy to recast everything (I don't think I can monkey-patch the class...)

[charris]

I had a toy implementation and it worked fine with 3.5 as far as I tested. I'm working on an 'official' version now for Numpy 1.10.

[abalkin]

@nicktimko - if you would like to contribute pure python code to this effort, try implementing __matmul__ for masked arrays. (See numpy.ma subpackage.)

[nicktimko]

What's the procedure for writing tests that involve new syntax? Just don't use it and manually call the method, i.e. a.__matmul__(b) instead of a @ b?

[pv]

Check Python version, and use eval("a @ b")?

[njsmith]

Both! (Because we'd like the tests to actually exercise the code now, even
before we have 3.5 on our CI servers, but we'd also like to test that the
syntax actually works.)
On Apr 7, 2015 9:15 AM, "Pauli Virtanen" [email protected] wrote:

Check Python version, and use eval("a @ b")?

—
Reply to this email directly or view it on GitHub
#4464 (comment).

[charris]

We would want to test both older an newer versions of Python. For the older
versions we would do something like the operator.matmul. Not sure what
the best way to test the new syntax is, maybe have a test install that is
version dependent?

On Tue, Apr 7, 2015 at 10:09 AM, Nick Timkovich [email protected]
wrote:

What's the procedure for writing tests that involve new syntax? Just don't
use it and manually call the method, i.e. a.matmul(b) instead of a @ b
?

—
Reply to this email directly or view it on GitHub
#4464 (comment).

[seberg]

Not sure if it is practical, but I guess you could also avoid writing "@" completly by doing something like this (not quite sure I got all the magic python does right):

try:
    from operator import matmul
except ImportError:
    def matmul(a, b):
        # subtypes have priority:
        reversed_first = issubtype(type(b), type(a))
        if reversed_first:
            res = b.__rmatmul__(a)
        else:
            res = a.__matmul__(b)

        if res is not NotImplemented:
            return res
        # Try the opposite way around, since NotImplemented was returned:
        if reversed_first:
             ret = a.__matmul__(b)
        else:
            ret = b.__rmatmul__(a)

        if ret is not NotImplemented:
            return ret
        return False

Edit: had forgotten to return False ;).

[nicktimko]

Can you capture NotImplemented as a return value from an unimplemented __matmul__? I thought only the comparisons (e.g. __lt__) give those, and for other magic operator methods it's just a proper AttributeError.

[seberg]

Oh, yeah....

[asmeurer]

Operators seem to return NotImplemented for me

In [1]: (1).__add__([])
Out[1]: NotImplemented

[argriffing]

Operators seem to return NotImplemented for me

Maybe I'm not following correctly, but I thought the problem was something like this:

>>> from operator import matmul
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: cannot import name matmul
>>> (1).__add__([])
NotImplemented
>>> (1).__matmul__([])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'int' object has no attribute '__matmul__'

[asmeurer]

Ah, I see. That is even in Python 3.5.

[nicktimko]

Would it work to directly dig out the method for some tests (a_array.__matmul__(b_array), et al.), and do the syntax based on something like:

matrix_op_syntax = True
try:
    compile('x @ x', '<string>', 'eval')
except SyntaxError:
    matrix_op_syntax = False

...

def test_syntax_whatever(self):
    if not matrix_op_syntax:
        raise nose.SkipTest('Matrix operator syntax not supported')
    c = eval('a @ b')

[njsmith]

Sure, that would work.
On Apr 9, 2015 5:17 PM, "Nick Timkovich" [email protected] wrote:

Would it work to directly dig out the method for some tests (
a_array.matmul(b_array), et al.), and do the syntax based on
something like:

matrix_op_syntax = True
try:
compile('x @ x', '', 'eval')
except SyntaxError:
matrix_op_syntax = False

...

def test_syntax_whatever(self):
if not matrix_op_syntax:
raise nose.SkipTest('Matrix operator syntax not supported')
c = eval('a @ b')

—
Reply to this email directly or view it on GitHub
#4464 (comment).

[argriffing]

This is just a thought, but what about using a module that is not imported unless you know the Python version can deal with a @ b? The idea would be that Pythons that don't understand the syntax wouldn't even look at the file, so you could avoid wrapping the tests in eval.

[charris]

On Mon, Apr 6, 2015 at 4:16 PM, Nick Timkovich [email protected]
wrote:

Python 3.5.0a3 has added the @ binary operator to the syntax, is there
some NumPy branch available to test with? Naively trying to subclass from
np.array fails because np.ndarray is actually the type and it gets really
messy to recast everything (I don't think I can monkey-patch the class...)

BTW, the class I used to experiment with is an attachment to this post
http://article.gmane.org/gmane.comp.python.numeric.general/58500/match=matmul
.

Chuck

[rgommers]

@charris maybe put it up as a gist for convenience?

[charris]

@nicktimko I've implemented __matmul__ and __rmatmul__ in #5878. The methods seem to work, but they do not get called by the @ operator. Perhaps that is because ndarray is a type defined in C.

[shoyer]

This was fixed by #5878