Fix conv and add convolve (cvxpy#2047)

* make conv more like np.convolve

* more tests

* make backwards compatible

* add convolve, deprecate conv

* added convolve, deprecated conv

* update docs

cvxpy/atoms/__init__.py

@@ -18,7 +18,7 @@
                                                  scalar_product,)
 from cvxpy.atoms.affine.bmat import bmat
 from cvxpy.atoms.affine.conj import conj
-from cvxpy.atoms.affine.conv import conv
+from cvxpy.atoms.affine.conv import conv, convolve
 from cvxpy.atoms.affine.cumsum import cumsum
 from cvxpy.atoms.affine.diag import diag
 from cvxpy.atoms.affine.diff import diff

cvxpy/atoms/affine/conv.py

@@ -14,6 +14,7 @@
 limitations under the License.
 """
 
+import warnings
 from typing import List, Tuple
 
 import numpy as np
@@ -44,19 +45,21 @@ class conv(AffAtom):
     rh_expr : Expression
         A 1D vector or a 2D column vector.
     """
-    # TODO work with right hand constant.
-    # TODO(akshayka): make DGP-compatible
 
     def __init__(self, lh_expr, rh_expr) -> None:
+        warnings.warn("conv is deprecated. Use convolve instead.", DeprecationWarning)
         super(conv, self).__init__(lh_expr, rh_expr)
 
     @AffAtom.numpy_numeric
     def numeric(self, values):
         """Convolve the two values.
         """
-        # Convert values to 1D.
-        values = list(map(intf.from_2D_to_1D, values))
-        return np.convolve(values[0], values[1])
+        flat_values = list(map(intf.from_2D_to_1D, values))
+        output = np.convolve(flat_values[0], flat_values[1])
+        if values[0].ndim == 2 or values[1].ndim == 2:
+            return output[:, None]
+        else:
+            return output
 
     def validate_arguments(self) -> None:
         """Checks that both arguments are vectors, and the first is constant.
@@ -69,9 +72,94 @@ def validate_arguments(self) -> None:
     def shape_from_args(self) -> Tuple[int, int]:
         """The sum of the argument dimensions - 1.
         """
-        lh_length = self.args[0].shape[0]
-        rh_length = self.args[1].shape[0]
-        return (lh_length + rh_length - 1, 1)
+        lh_length = self.args[0].size
+        rh_length = self.args[1].size
+        output_length = lh_length + rh_length - 1
+        if self.args[0].ndim == 2 or self.args[1].ndim == 2:
+            return (output_length, 1)
+        else:
+            return (output_length,)
+
+    def sign_from_args(self) -> Tuple[bool, bool]:
+        """Same as times.
+        """
+        return u.sign.mul_sign(self.args[0], self.args[1])
+
+    def is_incr(self, idx) -> bool:
+        """Is the composition non-decreasing in argument idx?
+        """
+        return self.args[0].is_nonneg()
+
+    def is_decr(self, idx) -> bool:
+        """Is the composition non-increasing in argument idx?
+        """
+        return self.args[0].is_nonpos()
+
+    def graph_implementation(
+        self, arg_objs, shape: Tuple[int, ...], data=None
+    ) -> Tuple[lo.LinOp, List[Constraint]]:
+        """Convolve two vectors.
+
+        Parameters
+        ----------
+        arg_objs : list
+            LinExpr for each argument.
+        shape : tuple
+            The shape of the resulting expression.
+        data :
+            Additional data required by the atom.
+
+        Returns
+        -------
+        tuple
+            (LinOp for objective, list of constraints)
+        """
+        return (lu.conv(arg_objs[0], arg_objs[1], shape), [])
+
+
+class convolve(AffAtom):
+    """ 1D discrete convolution of two vectors.
+
+    The discrete convolution :math:`c` of vectors :math:`a` and :math:`b` of
+    lengths :math:`n` and :math:`m`, respectively, is a length-:math:`(n+m-1)`
+    vector where
+
+    .. math::
+
+        c_k = \\sum_{i+j=k} a_ib_j, \\quad k=0, \\ldots, n+m-2.
+
+    Matches numpy.convolve
+
+    Parameters
+    ----------
+    lh_expr : Constant
+        A constant scalar or 1D vector.
+    rh_expr : Expression
+        A scalar or 1D vector.
+    """
+    # TODO work with right hand constant.
+    # TODO(akshayka): make DGP-compatible
+
+    @AffAtom.numpy_numeric
+    def numeric(self, values):
+        """Convolve the two values.
+        """
+        return np.convolve(values[0], values[1])
+
+    def validate_arguments(self) -> None:
+        """Checks that both arguments are vectors, and the first is constant.
+        """
+        if not self.args[0].ndim <= 1 or not self.args[1].ndim <= 1:
+            raise ValueError("The arguments to conv must be scalar or 1D.")
+        if not self.args[0].is_constant():
+            raise ValueError("The first argument to conv must be constant.")
+
+    def shape_from_args(self) -> Tuple[int, int]:
+        """The sum of the argument dimensions - 1.
+        """
+        lh_length = self.args[0].size
+        rh_length = self.args[1].size
+        return (lh_length + rh_length - 1,)
 
     def sign_from_args(self) -> Tuple[bool, bool]:
         """Same as times.

cvxpy/tests/test_atoms.py

@@ -861,6 +861,24 @@ def test_kron_expr(self) -> None:
         self.assertEqual(str(cm.exception),
                          "At least one argument to kron must be constant.")
 
+    def test_convolve(self) -> None:
+        """Test the convolve atom.
+        """
+        a = np.ones((3,))
+        b = Parameter(2, nonneg=True)
+        expr = cp.convolve(a, b)
+        assert expr.is_nonneg()
+        self.assertEqual(expr.shape, (4,))
+        b = Parameter(2, nonpos=True)
+        expr = cp.convolve(a, b)
+        assert expr.is_nonpos()
+        with self.assertRaises(Exception) as cm:
+            cp.convolve(self.x, -1)
+        self.assertEqual(str(cm.exception),
+                         "The first argument to conv must be constant.")
+        with pytest.raises(ValueError, match="scalar or 1D"):
+            cp.convolve([[0, 1], [0, 1]], self.x)
+
     def test_partial_optimize_dcp(self) -> None:
         """Test DCP properties of partial optimize.
         """

cvxpy/tests/test_convolution.py

@@ -15,6 +15,7 @@
 """
 
 import numpy as np
+import pytest
 
 import cvxpy as cvx
 import cvxpy.problems.iterative as iterative
@@ -31,29 +32,95 @@ def test_1D_conv(self) -> None:
         """
         n = 3
         x = cvx.Variable(n)
-        f = [1, 2, 3]
-        g = [0, 1, 0.5]
-        f_conv_g = [0., 1., 2.5,  4., 1.5]
-        expr = cvx.conv(f, g)
+        f = np.array([1, 2, 3])
+        g = np.array([0, 1, 0.5])
+        f_conv_g = np.array([0., 1., 2.5,  4., 1.5])
+        with pytest.warns(DeprecationWarning, match="Use convolve"):
+            expr = cvx.conv(f, g)
         assert expr.is_constant()
-        self.assertEqual(expr.shape, (5, 1))
+        self.assertEqual(expr.shape, (5,))
+        self.assertEqual(expr.shape, expr.value.shape)
         self.assertItemsAlmostEqual(expr.value, f_conv_g)
 
         expr = cvx.conv(f, x)
         assert expr.is_affine()
+        self.assertEqual(expr.shape, (5,))
+        # Matrix stuffing.
+        prob = cvx.Problem(cvx.Minimize(cvx.norm(expr, 1)),
+                           [x == g])
+        result = prob.solve(solver=cvx.SCS)
+        self.assertAlmostEqual(result, sum(f_conv_g), places=3)
+        self.assertItemsAlmostEqual(expr.value, f_conv_g)
+
+        # Test other shape configurations.
+        expr = cvx.conv(2, g)
+        self.assertEqual(expr.shape, (3,))
+        self.assertEqual(expr.shape, expr.value.shape)
+        self.assertItemsAlmostEqual(expr.value, 2 * g)
+
+        expr = cvx.conv(f, 2)
+        self.assertEqual(expr.shape, (3,))
+        self.assertEqual(expr.shape, expr.value.shape)
+        self.assertItemsAlmostEqual(expr.value, 2 * f)
+
+        expr = cvx.conv(f, g[:, None])
+        self.assertEqual(expr.shape, (5, 1))
+        self.assertEqual(expr.shape, expr.value.shape)
+        self.assertItemsAlmostEqual(expr.value, f_conv_g)
+
+        expr = cvx.conv(f[:, None], g)
+        self.assertEqual(expr.shape, (5, 1))
+        self.assertEqual(expr.shape, expr.value.shape)
+        self.assertItemsAlmostEqual(expr.value, f_conv_g)
+
+        expr = cvx.conv(f[:, None], g[:, None])
         self.assertEqual(expr.shape, (5, 1))
+        self.assertEqual(expr.shape, expr.value.shape)
+        self.assertItemsAlmostEqual(expr.value, f_conv_g)
+
+    def test_convolve(self) -> None:
+        """Test convolve.
+        """
+        n = 3
+        x = cvx.Variable(n)
+        f = np.array([1, 2, 3])
+        g = np.array([0, 1, 0.5])
+        f_conv_g = np.array([0., 1., 2.5,  4., 1.5])
+        expr = cvx.convolve(f, g)
+        assert expr.is_constant()
+        self.assertEqual(expr.shape, (5,))
+        self.assertEqual(expr.shape, expr.value.shape)
+        self.assertItemsAlmostEqual(expr.value, f_conv_g)
+
+        expr = cvx.convolve(f, x)
+        assert expr.is_affine()
+        self.assertEqual(expr.shape, (5,))
         # Matrix stuffing.
         prob = cvx.Problem(cvx.Minimize(cvx.norm(expr, 1)),
                            [x == g])
         result = prob.solve(solver=cvx.SCS)
         self.assertAlmostEqual(result, sum(f_conv_g), places=3)
         self.assertItemsAlmostEqual(expr.value, f_conv_g)
 
-        # # Expression trees.
-        # prob = Problem(Minimize(norm(expr, 1)))
-        # self.prob_mat_vs_mul_funcs(prob)
-        # result = prob.solve(solver=SCS, expr_tree=True, verbose=True)
-        # self.assertAlmostEqual(result, 0, places=1)
+        # Test other shape configurations.
+        expr = cvx.convolve(2, g)
+        self.assertEqual(expr.shape, (3,))
+        self.assertEqual(expr.shape, expr.value.shape)
+        self.assertItemsAlmostEqual(expr.value, 2 * g)
+
+        expr = cvx.convolve(f, 2)
+        self.assertEqual(expr.shape, (3,))
+        self.assertEqual(expr.shape, expr.value.shape)
+        self.assertItemsAlmostEqual(expr.value, 2 * f)
+
+        with pytest.raises(ValueError, match="must be scalar or 1D"):
+            expr = cvx.convolve(f, g[:, None])
+
+        with pytest.raises(ValueError, match="must be scalar or 1D"):
+            expr = cvx.convolve(f[:, None], g)
+
+        with pytest.raises(ValueError, match="must be scalar or 1D"):
+            expr = cvx.convolve(f[:, None], g[:, None])
 
     def prob_mat_vs_mul_funcs(self, prob) -> None:
         data, dims = prob.get_problem_data(solver=cvx.SCS)
@@ -101,24 +168,37 @@ def mat_from_func(self, func, rows, cols):
     def test_conv_prob(self) -> None:
         """Test a problem with convolution.
         """
-        import numpy as np
         N = 5
+        # Test conv.
         y = np.random.randn(N, 1)
         h = np.random.randn(2, 1)
-        x = cvx.Variable(N)
+        x = cvx.Variable((N, 1))
         v = cvx.conv(h, x)
         obj = cvx.Minimize(cvx.sum(cvx.multiply(y, v[0:N])))
-        print(cvx.Problem(obj, []).solve(solver=cvx.ECOS))
+        prob = cvx.Problem(obj, [])
+        prob.solve(solver=cvx.ECOS)
+        assert prob.status is cvx.UNBOUNDED
+
+        # Test convolve.
+        y = np.random.randn(N)
+        h = np.random.randn(2)
+        x = cvx.Variable((N))
+        v = cvx.convolve(h, x)
+        obj = cvx.Minimize(cvx.sum(cvx.multiply(y, v[0:N])))
+        prob = cvx.Problem(obj, [])
+        prob.solve(solver=cvx.ECOS)
+        assert prob.status is cvx.UNBOUNDED
 
     def test_0D_conv(self) -> None:
         """Convolution with 0D input.
         """
-        x = cvx.Variable((1,))  # or cvx.Variable((1,1))
-        problem = cvx.Problem(
-            cvx.Minimize(
-                cvx.max(cvx.conv([1.], cvx.multiply(1., x)))
-            ),
-            [x >= 0]
-        )
-        problem.solve(cvx.ECOS)
-        assert problem.status == cvx.OPTIMAL
+        for func in [cvx.conv, cvx.convolve]:
+            x = cvx.Variable((1,))  # or cvx.Variable((1,1))
+            problem = cvx.Problem(
+                cvx.Minimize(
+                    cvx.max(func(1., cvx.multiply(1., x)))
+                ),
+                [x >= 0]
+            )
+            problem.solve(cvx.ECOS)
+            assert problem.status == cvx.OPTIMAL

doc/source/tutorial/functions/index.rst

@@ -770,7 +770,7 @@ Expression and a negative Expression) then the returned Expression will have unk
      - |affine| affine
      - |incr| incr.
 
-   * - :ref:`conv(c, x) <conv>`
+   * - :ref:`convolve(c, x) <convolve>`
 
        :math:`c\in\mathbf{R}^m`
      - :math:`c*x`
@@ -880,7 +880,7 @@ The input to :math:`\texttt{bmat}` is a list of lists of CVXPY expressions.
 It constructs a block matrix.
 The elements of each inner list are stacked horizontally and then the resulting block matrices are stacked vertically.
 
-The output :math:`y = \mathbf{conv}(c, x)` has size :math:`n+m-1` and is defined as
+The output :math:`y = \mathbf{convolve}(c, x)` has size :math:`n+m-1` and is defined as
 :math:`y_k =\sum_{j=0}^{k} c[j]x[k-j]`.
 
 The output :math:`y = \mathbf{vec}(X)` is the matrix :math:`X` flattened in column-major order into a vector.