Cvxpy: Cannot use matmul on 1x1 matrix

[wkschwartz]

In cvxpy 1.0.6 and Python 3.7, I would expect the following example to work

import cvxpy, numpy
v = cvxpy.Variable((1,1))
col_scalar = numpy.array([[2]])
assert v.shape == col_scalar.shape == col_scalar.T.shape
print(repr(col_scalar.T @ col_scalar))
print(v.is_scalar())
print(repr(col_scalar.T @ v))

Instead I get

array([[4]])
True
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-99-86b127272bd5> in <module>()
      5 print(repr(col_scalar.T @ col_scalar))
      6 print(v.is_scalar())
----> 7 print(repr(col_scalar.T @ v))

~/Documents/wdp2/jupyter-venv/lib/python3.7/site-packages/cvxpy/expressions/expression.py in cast_op(self, other)
     47         """
     48         other = self.cast_to_const(other)
---> 49         return binary_op(self, other)
     50     return cast_op
     51 

~/Documents/wdp2/jupyter-venv/lib/python3.7/site-packages/cvxpy/expressions/expression.py in __rmatmul__(self, other)
    436         """Expression : Called for matrix @ Expression.
    437         """
--> 438         return other.__matmul__(self)
    439 
    440     def __neg__(self):

~/Documents/wdp2/jupyter-venv/lib/python3.7/site-packages/cvxpy/expressions/expression.py in cast_op(self, other)
     47         """
     48         other = self.cast_to_const(other)
---> 49         return binary_op(self, other)
     50     return cast_op
     51 

~/Documents/wdp2/jupyter-venv/lib/python3.7/site-packages/cvxpy/expressions/expression.py in __matmul__(self, other)
    395         """
    396         if self.is_scalar() or other.is_scalar():
--> 397             raise ValueError("Scalar operands are not allowed, use '*' instead")
    398         return self.__mul__(other)
    399 

ValueError: Scalar operands are not allowed, use '*' instead

Notice that the numpy 1x1 matrix multiplication worked and that the 1x1 Variable object says that it's a scalar.

Because of this I have to special case for when there's one variable versus more than one.

[rileyjmurray]

I agree that this behavior should be fixed. The recommended semantics of the @ operator are given in https://legacy.python.org/dev/peps/pep-0465/#semantics.

Right now Expression.is_scalar() simply returns all( d == 1 for d in self.shape ). As a result, Expressions of shape (1,1), (1,), and () are all considered scalars in cvxpy. Meanwhile the recommended semantics of the @ operator use the term "scalar" only to mean things with .shape == ().
So @wkschwartz it's easy enough to fix your particular example. We would just need to change the check on line 396 of expression.py to check self.shape == () or other.shape == () rather than using cvxpy's is_scalar() function. I have tested this change on my fork of cvxpy and all unittests still pass.

Unfortunately the above one-line change will not completely resolve the discrepancies between cvxpy's implementation of @ and the recommended semantics (which have been adopted by numpy, pandas, blaze, and theano, among others). I haven't determined the full scope of the differences, but the following example

import cvxpy
import numpy as np
a = np.array([1])
x = cvxpy.Variable(shape=(1,)
expr = a @ x

has expr.shape == (1,), instead of expr.shape == () as per the recommended semantics.

[wkschwartz]

@rileyjmurray Do you have a branch or pull request where you've made this change? (I just poked around your fork and didn't see anything.)

[rileyjmurray]

@wkschwartz I've lost track of where I made these changes. (I have two computers I work on. From time to time I have some uncommitted changes on one machine and I don't remember what they are, so I just reset to the most recent commit.)

Luckily, my comment was detailed enough that it shouldn't be hard for someone to recreate the changes. The main reason why I didn't create a PR is because I consider this issue unresolved until cvxpy completely matches the recommended semantics of the @ operator (I don't know what changes will be necessary for that to happen).

[rileyjmurray]

Follow-up:

The example below shows that when using the matmul operator between an n-by-n matrix and an n-by-1 vector, cvxpy drops the trailing 1 dimension of the result. This is a meaningful departure from the recommended semantics of the @ operator, and one that will affect cvxpy users (who might inadvertently create a malformed constraint expr1 == expr2 when expr1.shape == (n, 1) and expr2.shape == (n,)).

import numpy as np
import cvxpy
A = np.random.randn(4,4)
z = cvxpy.Variable(4, 1)
A @ z
>>> Expression(AFFINE, UNKNOWN, (4,))

[wkschwartz]

Can all of the shape semantics of @ be fixed in https://github.com/cvxgrp/cvxpy/blob/16ad9adad944d4cb34275e2f7dd0e27788b47b18/cvxpy/expressions/expression.py#L390-L396 by .reshapeing the result, or does this have to be fixed in https://github.com/cvxgrp/cvxpy/blob/16ad9adad944d4cb34275e2f7dd0e27788b47b18/cvxpy/expressions/expression.py#L378-L388?

[wkschwartz]

Actually, maybe it's best fixed here:
https://github.com/cvxgrp/cvxpy/blob/16ad9adad944d4cb34275e2f7dd0e27788b47b18/cvxpy/atoms/affine/binary_operators.py#L110-L113
and here:
https://github.com/cvxgrp/cvxpy/blob/16ad9adad944d4cb34275e2f7dd0e27788b47b18/cvxpy/atoms/affine/binary_operators.py#L274-L277

[SteveDiamond]

This is a bad situation. I hadn't realized there were so many failed cases. I'll take a look when I have time, and at least make @ cover the test cases here.

[wkschwartz]

The "Semantics" section of PEP 465 may be a source of some test cases

[SteveDiamond]

This should be fixed now. These semantics are more painful than they appear, once they start interacting with the rest of the system. Maybe MATLAB is right to have everything as a matrix.

[SteveDiamond]

Also @rileyjmurray in you're example above you accidentally defined a Variable of shape (4,) with name 1. This should probably raise a warning.

[SteveDiamond]

I made naming a variable an integer an error.