18

In python 3.5, the @ operator was introduced for matrix multiplication, following PEP465. This is implemented e.g. in numpy as the matmul operator.

However, as proposed by the PEP, the numpy operator throws an exception when called with a scalar operand:

>>> import numpy as np
>>> np.array([[1,2],[3,4]]) @ np.array([[1,2],[3,4]])    # works
array([[ 7, 10],
       [15, 22]])
>>> 1 @ 2                                                # doesn't work
Traceback (most recent call last):
  File "<input>", line 1, in <module>
TypeError: unsupported operand type(s) for @: 'int' and 'int'

This is a real turnoff for me, since I'm implementing numerical signal processing algorithms that should work for both scalars and matrices. The equations for both cases are mathematically exactly equivalent, which is no surprise, since "1-D x 1-D matrix multiplication" is equivalent to scalar multiplication. The current state however forces me to write duplicate code in order to handle both cases correctly.

So, given that the current state is not satisfactory, is there any reasonable way I can make the @ operator work for scalars? I thought about adding a custom __matmul__(self, other) method to scalar data types, but this seems like a lot of hassle considering the number of involved internal data types. Could I change the implementation of the __matmul__ method for numpy array data types to not throw an exception for 1x1 array operands?

And, on a sidenote, which is the rationale behind this design decision? Off the top of my head, I cannot think of any compelling reasons not to implement that operator for scalars as well.

Eike P.
  • 3,333
  • 1
  • 27
  • 38
  • 1
    how about `[1] @ [2]` ? Scalars have already `*` so why to duplicate it. – furas Dec 16 '16 at 14:21
  • 9
    It sounds like the real problem is that your code sometimes returns scalars and sometimes returns matrices. Why not refactor so your code returns 1 x 1 matrices instead of scalars? Or write a quick function that takes a matrix or scalar and returns that matrix or a 1x1 matrix with the scalar in it – Patrick Haugh Dec 16 '16 at 14:23
  • Why can't you use a try - except routine? – Jalo Dec 16 '16 at 14:23
  • 4
    As to why, from the PEP you linked: `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.` – Patrick Haugh Dec 16 '16 at 14:25
  • @furas `[1] @ [2]` doesn't work. Interestingly, `np.array([1]) @ np.array([2])` does (which I didn't notice before), while `np.array(1) @ np.array(2)` does not, although `type()` of both expressions yields `np.ndarray`. – Eike P. Dec 16 '16 at 14:29
  • @PatrickHaugh Nice hint, as I just commented, `np.array([1]) @ np.array([2])` does work, which I didn't notice before. – Eike P. Dec 16 '16 at 14:31
  • @PatrickHaugh and furas - I saw that note in the PEP, but I think I stated my reasons for critisizing it in the question. Actually, is it really an algorithmically distinct operation? Mathematically, it is equivalent. – Eike P. Dec 16 '16 at 14:34
  • @Jalo I could, but the whole point of my question was to make code for signal processing applications as concise as possible, to stay close to the underlying mathematical algorithms. (Which btw is in the spirit of PEP465.) – Eike P. Dec 16 '16 at 14:46
  • 2
    when I wrote `[1] @ [2]` I was thinking about something like `np.array([1]) @ np.array([2])` but I should describe this :) – furas Dec 16 '16 at 15:02
  • 1
    Perhaps you could use `atleast_1d`, e.g. `np.atleast_1d(5) @ np.atleast_1d(6)` works just fine. – Alex Riley Dec 16 '16 at 15:08
  • You'd pretty much have to recompile both Python and NumPy if you wanted to add `@` support to scalars. – user2357112 Dec 29 '16 at 05:28

1 Answers1

8

As ajcr suggested, you can work around this issue by forcing some minimal dimensionality on objects being multiplied. There are two reasonable options: atleast_1d and atleast_2d which have different results in regard to the type being returned by @: a scalar versus a 1-by-1 2D array. 

x = 3
y = 5
z = np.atleast_1d(x) @ np.atleast_1d(y)   # returns 15 
z = np.atleast_2d(x) @ np.atleast_2d(y)   # returns array([[15]])

However:

  • Using atleast_2d will lead to an error if x and y are 1D-arrays that would otherwise be multiplied normally
  • Using atleast_1d will result in the product that is either a scalar or a matrix, and you don't know which.
  • Both of these are more verbose than np.dot(x, y) which would handle all of those cases.

Also, the atleast_1d version suffers from the same flaw that would also be shared by having scalar @ scalar = scalar: you don't know what can be done with the output. Will z.T or z.shape throw an error? These work for 1-by-1 matrices but not for scalars. In the setting of Python, one simply cannot ignore the distinction between scalars and 1-by-1 arrays without also giving up all the methods and properties that the latter have.