Smith form affine (#13)

* adding cleaning up of torch

* adding jacobian structure and jacobian only value computation in interface

* making analytic center of ellipsoid work

* much more progress on everything, sparse jacobians and smith form

---------

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

Author: Transurgeon

cvxpy/atoms/affine/binary_operators.py

@@ -41,7 +41,6 @@
 class BinaryOperator(AffAtom):
     """
     Base class for expressions involving binary operators. (other than addition)
-
     """
 
     OP_NAME = 'BINARY_OP'
@@ -172,37 +171,40 @@ def is_decr(self, idx) -> bool:
         return self.args[1-idx].is_nonpos()
 
     def _grad(self, values):
-        """Gives the (sub/super)gradient of the atom w.r.t. each argument.
-
-        Matrix expressions are vectorized, so the gradient is a matrix.
-
+        """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)
+        
         Args:
-            values: A list of numeric values for the arguments.
-
+            values: A list of numeric values for the arguments [X, Y].
+        
         Returns:
-            A list of SciPy CSC sparse matrices or None.
+            A list of SciPy CSC sparse matrices [DX, DY].
         """
+        # Handle constant cases
         if self.args[0].is_constant() or self.args[1].is_constant():
             return super(MulExpression, self)._grad(values)
-
-        # TODO(akshayka): Verify that the following code is correct for
-        # non-affine arguments.
+        
         X = values[0]
         Y = values[1]
-
-        DX_rows = self.args[0].size
-        cols = self.args[0].size
-
-        # DX = [diag(Y11), diag(Y12), ...]
-        #      [diag(Y21), diag(Y22), ...]
-        #      [   ...        ...     ...]
-        DX = sp.dok_array((DX_rows, cols))
-        for k in range(self.args[0].shape[0]):
-            DX[k::self.args[0].shape[0], k::self.args[0].shape[0]] = Y
-        DX = sp.csc_array(DX)
-        cols = 1 if len(self.args[1].shape) == 1 else self.args[1].shape[1]
-        DY = sp.block_diag([X.T for k in range(cols)], 'csc')
-
+        
+        # Get dimensions
+        m, n = self.args[0].shape if len(self.args[0].shape) == 2 else (self.args[0].size, 1)
+        n2, p = self.args[1].shape if len(self.args[1].shape) == 2 else (self.args[1].size, 1)
+        
+        # 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(Y) = I_p ⊗ X
+        # This is a (m*p) × (n*p) matrix
+        DY = sp.kron(sp.eye(p, format='csc'), X, format='csc')
+        
         return [DX, DY]
 
     def graph_implementation(

cvxpy/atoms/atom.py

@@ -465,6 +465,10 @@ def _domain(self) -> List['Constraint']:
         """
         # Default is no constraints.
         return []
+    
+    def point_in_domain(self) -> np.ndarray:
+        """default point in domain of zero"""
+        return np.zeros(self.shape)
 
     @staticmethod
     def numpy_numeric(numeric_func):

cvxpy/atoms/elementwise/log.py

@@ -95,3 +95,8 @@ def _domain(self) -> List[Constraint]:
         """Returns constraints describing the domain of the node.
         """
         return [self.args[0] >= 0]
+
+    def point_in_domain(self) -> np.ndarray:
+        """Returns a point in the domain of the node.
+        """
+        return np.ones(self.size)

cvxpy/atoms/elementwise/trig.py

@@ -75,10 +75,13 @@ def _domain(self) -> List[Constraint]:
         """
         return []
 
-    def _grad(self) -> List[Constraint]:
+    def _grad(self, values) -> List[Constraint]:
         """Returns the gradient of the node.
         """
-        return []
+        rows = self.args[0].size
+        cols = self.size
+        grad_vals = np.cos(values[0])
+        return [sin.elemwise_grad_to_diag(grad_vals, rows, cols)]
 
 
 class cos(Elementwise):
@@ -135,10 +138,13 @@ def _domain(self) -> List[Constraint]:
         """
         return []
 
-    def _grad(self) -> List[Constraint]:
+    def _grad(self, values) -> List[Constraint]:
         """Returns the gradient of the node.
         """
-        return []
+        rows = self.args[0].size
+        cols = self.size
+        grad_vals = -np.sin(values[0])
+        return [cos.elemwise_grad_to_diag(grad_vals, rows, cols)]
 
 
 class tan(Elementwise):
@@ -195,8 +201,10 @@ def _domain(self) -> List[Constraint]:
         """
         return []
 
-    def _grad(self) -> List[Constraint]:
+    def _grad(self, values) -> List[Constraint]:
         """Returns the gradient of the node.
         """
-        return []
-    
+        rows = self.args[0].size
+        cols = self.size
+        grad_vals = 1/np.cos(values[0])**2
+        return [tan.elemwise_grad_to_diag(grad_vals, rows, cols)]

cvxpy/reductions/expr2smooth/canonicalizers/__init__.py

@@ -14,10 +14,12 @@
 limitations under the License.
 """
 from cvxpy.atoms import maximum
+from cvxpy.atoms.elementwise.log import log
 from cvxpy.atoms.elementwise.minimum import minimum
 from cvxpy.atoms.elementwise.power import power
 from cvxpy.atoms.pnorm import Pnorm
 from cvxpy.atoms.elementwise.abs import abs
+from cvxpy.reductions.expr2smooth.canonicalizers.log_canon import log_canon
 from cvxpy.reductions.expr2smooth.canonicalizers.minimum_canon import minimum_canon
 from cvxpy.reductions.expr2smooth.canonicalizers.abs_canon import abs_canon
 from cvxpy.reductions.expr2smooth.canonicalizers.pnorm_canon import pnorm_canon
@@ -28,8 +30,8 @@
     abs: abs_canon,
     maximum : maximum_canon,
     minimum: minimum_canon,
-    # log: log_canon,
+    log: log_canon,
     power: power_canon,
-    #Pnorm : pnorm_canon,
+    Pnorm : pnorm_canon,
     # inv: inv_pos_canon,
 }

cvxpy/reductions/expr2smooth/canonicalizers/log_canon.py

@@ -0,0 +1,23 @@
+"""
+Copyright 2025 CVXPY developers
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+
+from cvxpy.expressions.variable import Variable
+
+
+def log_canon(expr, args):
+    t = Variable(args[0].size)
+    t.value = expr.point_in_domain()
+    return expr.copy([t]), [t==args[0]]

cvxpy/reductions/expr2smooth/canonicalizers/power_canon.py

@@ -39,8 +39,11 @@ def power_canon(expr, args):
                 t.value = np.power(np.abs(x.value), p)
             return t, [t**(1/p) == x, t >= 0]
         elif p > 1:
-            return x**p, []
-        else:  # p < 0
-            raise ValueError(
-                "Power canonicalization does not support negative powers."
-            )
+            t = Variable(args[0].shape)
+            if args[0].value is not None:
+                t.value = np.power(args[0].value, p)
+            else:
+                t.value = expr.point_in_domain()
+            return expr.copy([t]), [t==args[0]]
+        else:
+            pass

cvxpy/reductions/expr2smooth/expr2smooth.py

@@ -18,6 +18,7 @@
 
 from cvxpy import problems
 from cvxpy.expressions.expression import Expression
+from cvxpy.expressions.variable import Variable
 from cvxpy.problems.objective import Minimize
 from cvxpy.reductions.canonicalization import Canonicalization
 from cvxpy.reductions.expr2smooth.canonicalizers import CANON_METHODS as smooth_canon_methods
@@ -76,7 +77,6 @@ def canonicalize_tree(self, expr, affine_above: bool) -> Tuple[Expression, list]
         -------
         A tuple of the canonicalized expression and generated constraints.
         """
-        # TODO don't copy affine expressions?
         affine_atom = type(expr) not in self.smooth_canon_methods
         canon_args = []
         constrs = []
@@ -108,56 +108,11 @@ def canonicalize_expr(self, expr, args, affine_above: bool) -> Tuple[Expression,
 
         if type(expr) in self.smooth_canon_methods:
             return self.smooth_canon_methods[type(expr)](expr, args)
-
+        """
+        elif hasattr(expr, "curvature") and expr.curvature == "AFFINE":
+            if type(expr) is Variable:
+                return expr, []
+            t = Variable(expr.shape)
+            return t, [t == expr.copy(args)]
+        """
         return expr.copy(args), []
-
-"""
-def example_max():
-    # Define variables
-    x = cp.Variable(1)
-    y = cp.Variable(1)
-    
-    objective = cp.Minimize(-cp.maximum(x,y))
-    
-    constraints = [x - 14 == 0, y - 6 == 0]
-    
-    problem = cp.Problem(objective, constraints)
-    return problem
-
-def example_sqrt():
-    # Define variables
-    x = cp.Variable(1)
-    
-    objective = cp.Minimize(cp.sqrt(x))
-    
-    constraints = [x - 4 == 0]
-    
-    problem = cp.Problem(objective, constraints)
-    return problem
-
-def example_pnorm_even():
-    # Define variables
-    x = cp.Variable(2)
-    
-    objective = cp.Minimize(cp.pnorm(x, p=2))
-    
-    constraints = [x[0] - 3 == 0, x[1] - 4 == 0]
-    
-    problem = cp.Problem(objective, constraints)
-    return problem
-
-def example_pnorm_odd():
-    # Define variables
-    x = cp.Variable(2)
-    
-    objective = cp.Minimize(cp.pnorm(x, p=3))
-    
-    constraints = [x[0] - 3 == 0, x[1] - 4 == 0]
-    
-    problem = cp.Problem(objective, constraints)
-    return problem
-
-prob = example_sqrt()
-new_problem, inverse = Expr2smooth(prob).apply(prob)
-print(new_problem)
-"""