address Parth's reviews

Author: Transurgeon

cvxpy/reductions/solvers/nlp_solvers/diff_engine/converters.py

@@ -31,11 +31,18 @@
 
 
 def convert_matmul(expr, children, var_dict, n_vars, param_dict):
-    """Convert matrix multiplication A @ f(x), f(x) @ A, or X @ Y."""
+    """Convert matrix multiplication A @ f(x), f(x) @ A, or X @ Y.
+
+    numpy convention: a 1D left operand acts as a row vector, a 1D right
+    operand acts as a column vector. The C engine normalizes (n,) to (1, n),
+    so 1D right operands need an explicit reshape to (n, 1).
+    """
     left_arg, right_arg = expr.args
 
     if left_arg.is_constant():
         A = left_arg.value
+        if A.ndim == 1:
+            A = A.reshape(1, -1)
         param_node = children[0] if left_arg.parameters() else None
         if sparse.issparse(A):
             return make_sparse_left_matmul(param_node, children[1], A)
@@ -51,17 +58,16 @@ def convert_matmul(expr, children, var_dict, n_vars, param_dict):
         return make_dense_right_matmul(param_node, children[0], A)
 
     else:
-        return _diffengine.make_matmul(children[0], children[1])
+        right_node = children[1]
+        if len(right_arg.shape) == 1:
+            right_node = _diffengine.make_reshape(right_node, right_arg.shape[0], 1)
+        return _diffengine.make_matmul(children[0], right_node)
 
 # TODO we should support sparse elementwise multiply at some point.
 def convert_multiply(expr, children, var_dict, n_vars, param_dict):
     """Convert elementwise multiplication."""
     left_arg, right_arg = expr.args
 
-    # TODO: would be nice to catch promote here so we correctly create a
-    # a scalar multiply. What is even the convention with promoting a parameter?
-    # This is a very deep question.
-
     if left_arg.is_constant():
         if left_arg.size == 1:
             return _diffengine.make_param_scalar_mult(children[0], children[1])

cvxpy/reductions/solvers/nlp_solvers/diff_engine/helpers.py

@@ -29,16 +29,20 @@ def normalize_shape(shape):
 def to_dense_float(value):
     """Convert a value to a dense float64 numpy array."""
     if sparse.issparse(value):
-        value = value.todense()
+        value = value.toarray()
     return np.asarray(value, dtype=np.float64)
 
 
 def chain_add(children):
-    """Chain multiple children with binary adds: a + b + c -> add(add(a, b), c)."""
-    result = children[0]
-    for child in children[1:]:
-        result = _diffengine.make_add(result, child)
-    return result
+    """Combine children with a balanced binary tree of adds.
+
+    Tree depth is ceil(log2(N)) instead of N-1, which keeps the AD graph shallow
+    when summing many terms.
+    """
+    if len(children) == 1:
+        return children[0]
+    mid = len(children) // 2
+    return _diffengine.make_add(chain_add(children[:mid]), chain_add(children[mid:]))
 
 
 def make_sparse_left_matmul(param_node, child, A):
@@ -53,8 +57,6 @@ def make_sparse_left_matmul(param_node, child, A):
 
 
 def make_dense_left_matmul(param_node, child, A):
-    if A.ndim == 1:
-        A = A.reshape(1, -1)
     m, n = A.shape
     return _diffengine.make_left_matmul(
         param_node, child, 'dense', A.flatten(order='C'), m, n)

cvxpy/reductions/solvers/nlp_solvers/uno_nlpif.py

@@ -197,9 +197,7 @@ def constraints_callback(x, constraint_values):
             # Define Jacobian callback
             # Signature: jacobian(x, jacobian_values) -> void
             def jacobian_callback(x, jacobian_values):
-                jac_vals = oracles.jacobian(x)
-                # Flatten in case it's returned as a 2D memoryview
-                jacobian_values[:] = np.asarray(jac_vals).flatten()
+                jacobian_values[:] = oracles.jacobian(x)
 
             # set_constraints(number_constraints, constraint_functions,
             #     constraints_lower_bounds, constraints_upper_bounds,
@@ -224,9 +222,7 @@ def jacobian_callback(x, jacobian_values):
         # Define Lagrangian Hessian callback
         # Signature: hessian(x, objective_multiplier, multipliers, hessian_values)
         def hessian_callback(x, objective_multiplier, multipliers, hessian_values):
-            hess_vals = oracles.hessian(x, multipliers, objective_multiplier)
-            # Flatten in case it's returned as a 2D array
-            hessian_values[:] = np.asarray(hess_vals).flatten()
+            hessian_values[:] = oracles.hessian(x, multipliers, objective_multiplier)
 
         # set_lagrangian_hessian(number_hessian_nonzeros, hessian_triangular_part,
         #     hessian_row_indices, hessian_column_indices, lagrangian_hessian)

cvxpy/tests/nlp_tests/test_affine_atoms.py

@@ -147,3 +147,22 @@ def test_convolve_alias(self):
 
         checker = DerivativeChecker(problem)
         checker.run_and_assert()
+
+    def test_matmul_1d_dot_product(self):
+        """x.T @ x for a 1D variable. Pins both the transpose 1D no-op and
+        the var-var matmul handling of a 1D right operand (numpy column-vector
+        convention)."""
+        np.random.seed(0)
+        n = 5
+        x = cp.Variable(n, bounds=[-1, 1], name='x')
+        x.value = np.random.rand(n)
+        assert (cp.sin(x).T @ cp.sin(x)).shape == ()
+
+        obj = cp.sin(x).T @ cp.sin(x)
+        problem = cp.Problem(cp.Minimize(obj))
+
+        problem.solve(solver=cp.IPOPT, nlp=True, verbose=False)
+        assert(problem.status == cp.OPTIMAL)
+
+        checker = DerivativeChecker(problem)
+        checker.run_and_assert()

pyproject.toml

@@ -85,7 +85,7 @@ dependencies = [
     "scipy >= 1.13.0",
     "highspy >= 1.11.0",
     "qdldl >= 0.1.7.post0",
-    "sparsediffpy >= 0.2.2",
+    "sparsediffpy >= 0.2.2, < 0.3.0",
 ]
 requires-python = ">=3.11"
 urls = {Homepage = "https://github.com/cvxpy/cvxpy"}