Remove Python _jacobian and _hess_vec oracle methods (#141)

These methods are now provided by the C-based diff_engine_core library.
The NLP solver Oracles class uses c_problem from diff_engine, not the
Python methods on atoms.

Removed from base classes:
- atom.py: jacobian(), hess_vec(), _jacobian(), _hess_vec(),
  _verify_jacobian_args(), _verify_hess_vec_args()
- variable.py: jacobian(), hess_vec()
- constant.py: jacobian(), hess_vec()

Removed from 13 affine atoms, 10 elementwise atoms, and 4 non-elementwise
atoms (pnorm, prod, quad_form, quad_over_lin).

Deleted test files:
- cvxpy/tests/nlp_tests/jacobian_tests/ (10 files)
- cvxpy/tests/nlp_tests/hess_tests/ (10 files)

Co-authored-by: Claude Opus 4.5 <noreply@anthropic.com>

Author: Transurgeon

cvxpy/atoms/affine/add_expr.py

@@ -18,7 +18,6 @@
 from typing import Any, Iterable, List, Tuple
 
 import numpy as np
-from scipy.sparse import coo_matrix
 
 import cvxpy.lin_ops.lin_op as lo
 import cvxpy.lin_ops.lin_utils as lu
@@ -156,81 +155,3 @@ def graph_implementation(
             if arg.shape != shape and lu.is_scalar(arg):
                 arg_objs[i] = lu.promote(arg, shape)
         return (lu.sum_expr(arg_objs), [])
-
-    def _verify_hess_vec_args(self):
-        return True
-
-    def _hess_vec(self, vec):
-        """
-        Computes the merged Hessian-vector product dictionary for all arguments.
-        If a key appears in several, their values are summed.
-        """
-        hess_dict = {}
-        keys_require_summing = []
-        
-        for arg in self.args:
-            if not arg.is_affine():
-                arg_hess = arg.hess_vec(vec)
-                for k, v in arg_hess.items():
-                    if k in hess_dict:
-                        hess_dict[k][0].extend(v[0])
-                        hess_dict[k][1].extend(v[1])
-                        hess_dict[k][2].extend(v[2])
-                        keys_require_summing.append(k)
-                    else:
-                        hess_dict[k] = ([], [], [])
-                        hess_dict[k][0].extend(np.atleast_1d(v[0]))
-                        hess_dict[k][1].extend(np.atleast_1d(v[1]))
-                        hess_dict[k][2].extend(np.atleast_1d(v[2]))
-
-        # sum duplicates
-        for key in set(keys_require_summing): 
-            rows, cols, vals = hess_dict[key]
-            shape = (key[0].size, key[0].size)
-            coo = coo_matrix((vals, (rows, cols)), shape=shape)
-            coo.sum_duplicates()
-            hess_dict[key] = (coo.row, coo.col, coo.data)
-
-        # convert lists to arrays
-        for k, v in hess_dict.items():
-            rows, cols, vals = v
-            hess_dict[k] = (np.array(rows), np.array(cols), np.array(vals))
-
-        return hess_dict
-
-    def _verify_jacobian_args(self):
-        return True
-
-    def _jacobian(self):
-        jacobian_dict = {}
-        keys_require_summing = []
-
-        for arg in self.args:
-            if not arg.is_constant():
-                arg_jac = arg.jacobian()
-            
-                for k, v in arg_jac.items():
-                    if k in jacobian_dict:
-                        jacobian_dict[k][0].extend(v[0])
-                        jacobian_dict[k][1].extend(v[1])
-                        jacobian_dict[k][2].extend(v[2])
-                        keys_require_summing.append(k)
-                    else:
-                        jacobian_dict[k] = ([], [], [])
-                        jacobian_dict[k][0].extend(np.atleast_1d(v[0]))
-                        jacobian_dict[k][1].extend(np.atleast_1d(v[1]))
-                        jacobian_dict[k][2].extend(np.atleast_1d(v[2]))
-
-        # sum duplicates
-        for key in set(keys_require_summing): 
-            rows, cols, vals = jacobian_dict[key]
-            coo = coo_matrix((vals, (rows, cols)), shape=(self.size, key.size))
-            coo.sum_duplicates()
-            jacobian_dict[key] = (coo.row, coo.col, coo.data)
-
-        # convert lists to arrays
-        for k, v in jacobian_dict.items():
-            rows, cols, vals = v
-            jacobian_dict[k] = (np.array(rows), np.array(cols), np.array(vals))
-
-        return jacobian_dict

cvxpy/atoms/affine/affine_atom.py

@@ -169,7 +169,4 @@ def _grad(self, values) -> List[Any]:
                 grad_list.append(grad_matrix[var_start:var_end, :])
                 var_start = var_end
 
-        return grad_list
-
-    def _verify_jacobian_args(self):
-        return True
+        return grad_list

cvxpy/atoms/affine/binary_operators.py

@@ -28,7 +28,6 @@
 from cvxpy.atoms.affine.affine_atom import AffAtom
 from cvxpy.atoms.affine.broadcast_to import broadcast_to
 from cvxpy.atoms.affine.conj import conj
-from cvxpy.atoms.affine.promote import Promote
 from cvxpy.atoms.affine.reshape import deep_flatten, reshape
 from cvxpy.atoms.affine.sum import sum as cvxpy_sum
 from cvxpy.constraints.constraint import Constraint
@@ -38,7 +37,6 @@
     is_param_free,
 )
 from cvxpy.expressions.expression import Expression
-from cvxpy.expressions.variable import Variable
 from cvxpy.utilities import bounds as bounds_utils
 
 
@@ -311,137 +309,6 @@ def _grad(self, values):
         DY = sp.kron(sp.eye_array(n), X, format='csc').T
 
         return [DX, DY]
-    
-    def _verify_hess_vec_args(self):
-        X = self.args[0]
-        Y = self.args[1]
-
-        # if X is an atom, Y must be constant
-        if not isinstance(X, Variable) and not X.is_constant():
-            if not Y.is_constant():
-                return False
-        
-        # if Y is an atom, X must be constant
-        if not isinstance(Y, Variable) and not Y.is_constant():
-            if not X.is_constant():
-                return False
-
-        # if both are variables, check that they are not the same variable
-        if isinstance(X, Variable) and isinstance(Y, Variable):
-            if X.id == Y.id:
-                return False
-
-        return True
-    
-    def _hess_vec(self, vec):
-        X = self.args[0]
-        Y = self.args[1]
-
-        m, n = self.get_dimensions(X)
-        _, p = self.get_dimensions(Y)
-
-        if X.is_constant():
-            B = X.value.T @ np.reshape(vec, (m, p), order='F')
-            hess_dict = Y.hess_vec(B.flatten(order='F'))
-            return hess_dict 
-            
-        if Y.is_constant():
-            B = np.reshape(vec, (m, p), order='F') @ Y.value.T
-            hess_dict = X.hess_vec(B.flatten(order='F'))
-            return hess_dict
-
-        # here both are variables by themselves so we only get a cross term
-        rows = np.tile(np.arange(m * n), p)
-        cols = np.repeat(np.arange(n * p), m)
-        vals = vec[(cols // n) * m + (rows % m)]
-        return {(X, Y): (rows, cols, vals), (Y, X): (cols, rows, vals)}
-            
-    def _verify_jacobian_args(self):
-        X = self.args[0]
-        Y = self.args[1]
-
-        X_vars = X.variables()
-        Y_vars = Y.variables()
-
-        # no variable can appear in both arguments
-        for x_var in X_vars:
-            for y_var in Y_vars:
-                if x_var.id == y_var.id:
-                    return False
-        
-        return True
-
-    def get_dimensions(self, X):
-        """Get the dimensions of X as (rows, cols).   
-        """
-        if len(X.shape) == 0:
-            return (1, 1)
-        elif len(X.shape) == 1:
-            return (X.shape[0], 1)
-        else:
-            return X.shape
-
-    def _jacobian(self):
-        """
-        The atom is phi(X, Y) = X @ Y. It is vectorized as
-        z = vec(phi(X, Y)) = (I ⊗ X) vec(Y) = (Y.T ⊗ I) vec(X).
-        Let x = vec(X) and y = vec(Y). Then the Jacobian is given by
-        dz/dx = kron(Y.T, I) and dz/dy = kron(I, X).
-        """
-
-        X = self.args[0]
-        Y = self.args[1]
-        
-        m, _ = self.get_dimensions(X)
-        _, p = self.get_dimensions(Y)
-
-        dx_dict = {}
-        dy_dict = {}
-
-        if not X.is_constant():
-            dx = sp.kron(Y.value.T, sp.eye(m), format='csr')
-            
-            if not isinstance(X, Variable):
-                X_jac_dict = X.jacobian()
-                for var in X_jac_dict:
-                    rows, cols, vals = X_jac_dict[var]
-                    X_jac = sp.coo_array((vals, (rows, cols)), 
-                                         shape=(dx.shape[1], var.size)).tocsc()
-                    X_jac = (dx @ X_jac).tocoo()
-                    X_jac_dict[var] = (X_jac.row, X_jac.col, X_jac.data)
-                
-                dx_dict = X_jac_dict 
-            else: 
-                dx = dx.tocoo()
-                dx_dict = {X: (dx.row, dx.col, dx.data)}
-            
-            
-        if not Y.is_constant():
-            dy = sp.kron(sp.eye(p), X.value, format='csr')
-
-            if not isinstance(Y, Variable):
-                Y_jac_dict = Y.jacobian()
-                for var in Y_jac_dict:
-                    rows, cols, vals = Y_jac_dict[var]
-                    Y_jac = sp.coo_array((vals, (rows, cols)),
-                                          shape=(dy.shape[1], var.size)).tocsc()
-                    Y_jac = (dy @ Y_jac).tocoo()
-                    Y_jac_dict[var] = (Y_jac.row, Y_jac.col, Y_jac.data)
-                
-                dy_dict = Y_jac_dict 
-            else: 
-                dy = dy.tocoo()
-                dy_dict = {Y: (dy.row, dy.col, dy.data)}
-
-        if X.is_constant() and not Y.is_constant():
-            return dy_dict
-        
-        if not X.is_constant() and Y.is_constant():
-            return dx_dict
-
-        # merge the two dictionaries together
-        dx_dict.update(dy_dict)
-        return dx_dict
 
     def graph_implementation(
         self, arg_objs, shape: Tuple[int, ...], data=None
@@ -574,114 +441,6 @@ def _grad(self, values):
 
         return [DX, DY]
 
-    def _verify_hess_vec_args(self):
-        x = self.args[0]
-        y = self.args[1]
-        if x.size != y.size:
-            return False
-
-        if x.is_constant() and y.is_constant():
-            return False
-
-        # one of the following must be true:
-        # 1. both arguments are variables
-        # 2. one argument is a constant
-        # 3. one argument is a Promote of a variable and the other is a variable
-        both_are_variables = isinstance(x, Variable) and isinstance(y, Variable)
-        one_is_constant = x.is_constant() or y.is_constant()
-        x_is_promote = type(x) == Promote and isinstance(y, Variable)
-        y_is_promote = type(y) == Promote and isinstance(x, Variable)
-
-        if not (both_are_variables or one_is_constant or x_is_promote or y_is_promote):
-            return False
-        
-        if both_are_variables and x.id == y.id:
-            return False
-
-        return True
-
-    def _hess_vec(self, vec):
-        x = self.args[0]
-        y = self.args[1]
-        
-        # constant * atom
-        if x.is_constant(): 
-            y_hess_vec = y.hess_vec(x.value.flatten(order='F') * vec)
-            return y_hess_vec
-        
-        # atom * constant
-        if y.is_constant():
-            x_hess_vec = x.hess_vec(y.value.flatten(order='F') * vec)
-            return x_hess_vec
-
-        # x * y with x a scalar variable, y a vector variable
-        if not isinstance(x, Variable) and x.is_affine():
-            assert(type(x) == Promote)
-            x_var = x.args[0] # here x is a Promote because of how we canonicalize
-            zeros_x = np.zeros(x_var.size, dtype=int)
-            cols = np.arange(y.size,  dtype=int)
-            return {(x_var, y): (zeros_x, cols, vec),
-                    (y, x_var): (cols, zeros_x, vec)}
-        
-        # x * y with x a vector variable, y a scalar
-        if not isinstance(y, Variable) and y.is_affine():
-            assert(type(y) == Promote)
-            y_var = y.args[0] # here y is a Promote because of how we canonicalize
-            zeros_y = np.zeros(y_var.size, dtype=int)
-            cols = np.arange(x.size,  dtype=int)
-            return {(x, y_var): (cols, zeros_y, vec),
-                    (y_var, x): (zeros_y, cols, vec)}
-        
-        # if we arrive here both arguments are variables of the same size
-        rows = np.arange(x.size,  dtype=int)
-        cols = np.arange(x.size,  dtype=int)
-        return {(x, y): (rows, cols, vec), (y, x): (rows, cols, vec)}
-
-    def _verify_jacobian_args(self):
-        return self._verify_hess_vec_args()
-    
-
-    def _jacobian(self):
-        x = self.args[0]
-        y = self.args[1]
-
-        if x.is_constant():
-            dy = y.jacobian()
-            for k in dy:
-                rows, cols, vals = dy[k]
-                # this is equivalent to forming the matrix defined
-                # rows, cols, vals and scaling each row i by y.value[i]
-                dy[k] = (rows, cols, np.atleast_1d(x.value).flatten(order='F')[rows] * vals)
-            return dy
-
-        if y.is_constant():
-            dx = x.jacobian()
-            for k in dx:
-                rows, cols, vals = dx[k]
-                dx[k] = (rows, cols, np.atleast_1d(y.value).flatten(order='F')[rows] * vals)
-            return dx
-
-        if not isinstance(x, Variable) and x.is_affine():
-            assert(type(x) == Promote)
-            x_var = x.args[0] # here x is a Promote because of how we canonicalize
-            idxs = np.arange(y.size,  dtype=int)
-            return {(x_var): (idxs, np.zeros(y.size, dtype=int), y.value),
-                    (y): (idxs, idxs, x.value)}
-        
-        # x * y with x a vector variable, y a scalar
-        if not isinstance(y, Variable) and y.is_affine():
-            assert(type(y) == Promote)
-            y_var = y.args[0] # here y is a Promote because of how we canonicalize
-            idxs = np.arange(x.size, dtype=int)
-            return {(x): (idxs, idxs, y.value),
-                    (y_var): (idxs, np.zeros(x.size, dtype=int), x.value)}
-        
-        # here both are variables
-        idxs = np.arange(x.size, dtype=int)
-        jacobian_dict = {x: (idxs, idxs, y.value.flatten(order='F')),
-                        y: (idxs, idxs, x.value.flatten(order='F'))}
-        return jacobian_dict
-
     def graph_implementation(
         self, arg_objs, shape: Tuple[int, ...], data=None
     ) -> Tuple[lo.LinOp, List[Constraint]]:
@@ -799,14 +558,6 @@ def is_decr(self, idx) -> bool:
     def point_in_domain(self):
         return np.ones(self.args[1].shape)
 
-    def _verify_hess_vec_args(self):
-        raise RuntimeError("The _verify_hess_vec_args method of"
-                           " the division atom should never be called.")
-
-    def _hess_vec(self, vec):
-        raise RuntimeError("The hess_vec method of the division atom should never "
-                           "be called.")
-
     def graph_implementation(
         self, arg_objs, shape: Tuple[int, ...], data=None
     ) -> Tuple[lo.LinOp, List[Constraint]]: