Fix matmul gradient issue (#14)

* fix matmul grad

* fix matmul and add sparsity in blocks

---------

Co-authored-by: William Zijie Zhang <[email protected]>

Author: Transurgeon

cvxpy/atoms/affine/binary_operators.py

@@ -174,8 +174,8 @@ def _grad(self, values):
         """Compute the gradient of matrix multiplication w.r.t. each argument.
         
         For Z = X @ Y, this computes the Jacobian matrices:
-        - ∂vec(Z)/∂vec(X) = Y.T ⊗ I_m  where X is (m, n)
-        - ∂vec(Z)/∂vec(Y) = I_p ⊗ X    where Y is (n, p)
+        - ∂vec(Z)/∂vec(X) = (Y.T ⊗ I_m).T  where X is (m, n)
+        - ∂vec(Z)/∂vec(Y) = (I_p ⊗ X).T    where Y is (n, p)
         
         Args:
             values: A list of numeric values for the arguments [X, Y].
@@ -197,13 +197,13 @@ def _grad(self, values):
         # Verify dimension compatibility
         assert n == n2, f"Inner dimensions must match for multiplication: {n} != {n2}"
 
-        # Compute ∂vec(Z)/∂vec(X) = Y.T ⊗ I_m
-        # This is a (m*p) × (m*n) matrix
-        DX = sp.kron(Y.T, sp.eye(m, format='csc'), format='csc')
+        # Compute ∂vec(Z)/∂vec(X) = (Y.T ⊗ I_m).T
+        # This is a (m*n) × (m*p) matrix 
+        DX = sp.kron(Y.T, sp.eye(m, format='csc'), format='csc').T
 
-        # Compute ∂vec(Z)/∂vec(Y) = I_p ⊗ X
-        # This is a (m*p) × (n*p) matrix
-        DY = sp.kron(sp.eye(p, format='csc'), X, format='csc')
+        # Compute ∂vec(Z)/∂vec(Y) = (I_p ⊗ X).T
+        # This is a (n*p) × (m*p) matrix
+        DY = sp.kron(sp.eye(p, format='csc'), X, format='csc').T
 
         return [DX, DY]
 

cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py

@@ -239,7 +239,9 @@ def jacobian(self, x):
                     if var in grad_dict:
                         jacobian = grad_dict[var].T
                         if sp.issparse(jacobian):
-                            jacobian = jacobian.toarray().flatten(order='F')
+                            jacobian = jacobian.tocoo()
+                            jacobian = jacobian.data
+                            #jacobian = jacobian.toarray().flatten(order='F')
                             values.append(np.atleast_1d(jacobian))
                         else:
                             values.append(np.atleast_1d(jacobian))
@@ -264,9 +266,9 @@ def jacobianstructure(self):
                     if var in grad_dict:
                         jacobian = grad_dict[var].T
                         if sp.issparse(jacobian):
-                            row_indices, col_indices = np.indices(jacobian.shape)
-                            rows.extend(row_indices.flatten(order='F') + row_offset)
-                            cols.extend(col_indices.flatten(order='F') + col_offset)
+                            jacobian = jacobian.tocoo()
+                            rows.extend(jacobian.row + row_offset)
+                            cols.extend(jacobian.col + col_offset)
                         else:
                             rows.extend(np.ones(jacobian.size)*row_offset)
                             cols.extend(np.arange(col_offset, col_offset + var.size))

cvxpy/sandbox/control_of_car.py

@@ -4,7 +4,7 @@
 import cvxpy as cp
 
 
-def solve_car_control(x_final, L=0.1, N=50, h=0.1, gamma=10):
+def solve_car_control(x_final, L=0.1, N=5, h=0.1, gamma=10):
     """
     Solve the nonlinear optimal control problem for car trajectory planning.
     

cvxpy/tests/test_nlp_solvers.py

@@ -194,11 +194,23 @@ def test_acopf(self):
             theta >= theta_lb,
             theta <= theta_ub
         ]
-        P_balance = cp.multiply(v, (G * cp.cos(theta_col - theta_col.T) + B * cp.sin(theta_col - theta_col.T)) @ v)
+        P_balance = cp.multiply(
+            v,
+            (
+                G * cp.cos(theta_col - theta_col.T)
+                + B * cp.sin(theta_col - theta_col.T)
+            ) @ v
+        )
         constraints.append(P == P_balance)
 
         # Reactive power balance
-        Q_balance = cp.multiply(v, (G * cp.sin(theta_col - theta_col.T) - B * cp.cos(theta_col - theta_col.T)) @ v)
+        Q_balance = cp.multiply(
+            v,
+            (
+                G * cp.sin(theta_col - theta_col.T)
+                - B * cp.cos(theta_col - theta_col.T)
+            ) @ v
+        )
         constraints.append(Q == Q_balance)
 
         # Objective: minimize reactive power at buses 1 and 3 (indices 0 and 2)