remove unnecessary point in dom + clean up C_problem + clean up oracle (#168)

* remove unnecessary point in dom + clean up C_problem + clean up oracle

* fix william's comments

Author: dance858

cvxpy/atoms/affine/binary_operators.py

@@ -560,9 +560,6 @@ def is_decr(self, idx) -> bool:
         else:
             return self.args[0].is_nonneg()
 
-    def point_in_domain(self):
-        return np.ones(self.args[1].shape)
-
     def graph_implementation(
         self, arg_objs, shape: Tuple[int, ...], data=None
     ) -> Tuple[lo.LinOp, List[Constraint]]:

cvxpy/atoms/atom.py

@@ -550,10 +550,6 @@ 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):
         """Wraps an atom's numeric function that requires numpy ndarrays as input.

cvxpy/atoms/elementwise/log.py

@@ -109,5 +109,4 @@ def _domain(self) -> List[Constraint]:
     def point_in_domain(self) -> np.ndarray:
         """Returns a point in the domain of the node.
         """
-        dim = (1, ) if self.size == 1 else self.shape
-        return np.ones(dim)
+        return np.ones(self.shape)

cvxpy/atoms/elementwise/logistic.py

@@ -82,6 +82,3 @@ def _grad(self, values):
         cols = self.size
         grad_vals = np.exp(values[0] - np.logaddexp(0, values[0]))
         return [logistic.elemwise_grad_to_diag(grad_vals, rows, cols)]
-
-    def point_in_domain(self):
-        return np.zeros(self.shape)

cvxpy/atoms/elementwise/power.py

@@ -400,11 +400,6 @@ def _domain(self) -> List[Constraint]:
         else:
             return []
 
-    def point_in_domain(self) -> np.ndarray:
-        """Returns a point in the domain of the node.
-        """
-        return np.ones(self.shape)
-
     def get_data(self):
         return [self._p_orig, self.max_denom]
 

cvxpy/reductions/solvers/nlp_solvers/copt_nlpif.py

@@ -122,13 +122,11 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         use_hessian = True
 
         if solver_cache is None:
-            oracles = Oracles(bounds.new_problem, bounds.x0, len(bounds.cl),
-                            verbose=verbose, use_hessian=use_hessian)
+            oracles = Oracles(bounds.new_problem, verbose=verbose, use_hessian=use_hessian)
         elif 'oracles' in solver_cache:
             oracles = solver_cache['oracles']
         else:
-            oracles = Oracles(bounds.new_problem, bounds.x0, len(bounds.cl),
-                            verbose=verbose, use_hessian=use_hessian)
+            oracles = Oracles(bounds.new_problem, verbose=verbose, use_hessian=use_hessian)
             solver_cache['oracles'] = oracles
 
         class COPTNlpCallbackCVXPY(copt.NlpCallbackBase):

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

@@ -16,7 +16,6 @@
 """
 
 import numpy as np
-from scipy import sparse
 
 import cvxpy as cp
 
@@ -42,31 +41,21 @@ def __init__(self, cvxpy_problem: cp.Problem, verbose: bool = True):
         c_obj = convert_expr(cvxpy_problem.objective.expr, var_dict, n_vars)
         c_constraints = [convert_expr(c.expr, var_dict, n_vars) for c in cvxpy_problem.constraints]
         self._capsule = _diffengine.make_problem(c_obj, c_constraints, verbose)
-        self._jacobian_allocated = False
-        self._hessian_allocated = False
 
-    def init_jacobian(self):
-        """Initialize Jacobian structures only. Must be called before jacobian()."""
-        _diffengine.problem_init_jacobian(self._capsule)
-        self._jacobian_allocated = True
-    
     def init_jacobian_coo(self):
-        """Initialize Jacobian COO structures only.
+        """Fill sparsity for the constraint Jacobian in COO format.
 
-        Must be called before get_jacobian_sparsity_coo().
+        Must be called once before get_jacobian_sparsity_coo() or eval_jacobian_vals().
         """
         _diffengine.problem_init_jacobian_coo(self._capsule)
-        self._jacobian_allocated = True
-
-    def init_hessian(self):
-        """Initialize Hessian structures only. Must be called before hessian()."""
-        _diffengine.problem_init_hessian(self._capsule)
-        self._hessian_allocated = True
 
     def init_hessian_coo_lower_tri(self):
-        """Initialize Hessian COO structures only. Must be called before get_hessian()."""
+        """Fill sparsity for the Lagrangian Hessian (lower triangle, COO).
+
+        Must be called once before get_problem_hessian_sparsity_coo() or
+        eval_hessian_vals_coo_lower_tri().
+        """
         _diffengine.problem_init_hessian_coo_lower_triangular(self._capsule)
-        self._hessian_allocated = True
 
     def objective_forward(self, u: np.ndarray) -> float:
         """Evaluate objective. Returns obj_value float."""
@@ -79,63 +68,43 @@ def constraint_forward(self, u: np.ndarray) -> np.ndarray:
     def gradient(self) -> np.ndarray:
         """Compute gradient of objective. Call objective_forward first. Returns gradient array."""
         return _diffengine.problem_gradient(self._capsule)
-
-    def jacobian(self) -> sparse.csr_matrix:
-        """Compute constraint Jacobian. Call constraint_forward first."""
-        data, indices, indptr, shape = _diffengine.problem_jacobian(self._capsule)
-        return sparse.csr_matrix((data, indices, indptr), shape=shape)
 
     def get_jacobian_sparsity_coo(self) -> tuple[np.ndarray, np.ndarray]:
-        """Get Jacobian sparsity pattern as COO. This function does not evaluate the jacobian."""
+        """Return the sparsity pattern (row, col) of the constraint Jacobian.
+
+        Does not evaluate the Jacobian; only returns structural nonzero indices.
+        Call init_jacobian_coo() first.
+        """
         rows, cols, unused_shape = _diffengine.get_jacobian_sparsity_coo(self._capsule)
         return rows, cols
 
     def eval_jacobian_vals(self) -> np.ndarray:
-        """Evaluate Jacobian values only.
+        """Evaluate the constraint Jacobian and return its nonzero values.
 
-        Call constraint_forward first. Returns jacobian values array.
+        The values correspond to the sparsity pattern from get_jacobian_sparsity_coo().
+        Call constraint_forward() first to set the evaluation point.
         """
         return _diffengine.problem_eval_jacobian_vals(self._capsule)
 
-    def get_jacobian(self) -> sparse.csr_matrix:
-        """Get constraint Jacobian. This function does not evaluate the jacobian. """
-        data, indices, indptr, shape = _diffengine.get_jacobian(self._capsule)
-        return sparse.csr_matrix((data, indices, indptr), shape=shape)
-
     def get_problem_hessian_sparsity_coo(self) -> tuple[np.ndarray, np.ndarray]:
-        """Get Hessian sparsity pattern as COO. This function does not evaluate the hessian."""
+        """Return the sparsity pattern (row, col) of the lower-triangular Lagrangian Hessian.
+
+        Does not evaluate the Hessian; only returns structural nonzero indices.
+        Call init_hessian_coo_lower_tri() first.
+        """
         rows, cols, unused_shape = _diffengine.get_problem_hessian_sparsity_coo(self._capsule)
         return rows, cols
-    
+
     def eval_hessian_vals_coo_lower_tri(
         self, obj_factor: float, lagrange: np.ndarray
     ) -> np.ndarray:
-        """Evaluate Hessian values only for lower triangular part.
-
-        Call objective_forward and constraint_forward first.
-        """
-        return _diffengine.problem_eval_hessian_vals_coo(self._capsule, obj_factor, lagrange)
+        """Evaluate the lower-triangular Lagrangian Hessian and return its nonzero values.
 
-    def hessian(self, obj_factor: float, lagrange: np.ndarray) -> sparse.csr_matrix:
-        """Compute Lagrangian Hessian.
+        Computes obj_factor * hess_f + sum(lagrange[i] * hess_gi), where f is the objective
+        and gi are the constraints. The values correspond to the sparsity pattern from
+        get_problem_hessian_sparsity_coo(). Only the lower triangle is returned.
 
-        Computes: obj_factor * H_obj + sum(lagrange_i * H_constraint_i)
-
-        Call objective_forward and constraint_forward before this.
-
-        Args:
-            obj_factor: Weight for objective Hessian
-            lagrange: Array of Lagrange multipliers (length = total_constraint_size)
-
-        Returns:
-            scipy CSR matrix of shape (n_vars, n_vars)
+        Call objective_forward() and constraint_forward() first to set the evaluation point.
         """
-        data, indices, indptr, shape = _diffengine.problem_hessian(
-            self._capsule, obj_factor, lagrange
-        )
-        return sparse.csr_matrix((data, indices, indptr), shape=shape)
+        return _diffengine.problem_eval_hessian_vals_coo(self._capsule, obj_factor, lagrange)
 
-    def get_hessian(self) -> sparse.csr_matrix:
-        """Get Lagrangian Hessian. This function does not evaluate the hessian."""
-        data, indices, indptr, shape = _diffengine.get_hessian(self._capsule)
-        return sparse.csr_matrix((data, indices, indptr), shape=shape)

cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py

@@ -76,12 +76,14 @@ def invert(self, solution, inverse_data):
         """
         Returns the solution to the original problem given the inverse_data.
         """
-        attr = {}
+        attr = {
+            s.NUM_ITERS: solution.get('num_iters'),
+        }
         status = self.STATUS_MAP[solution['status']]
         # the info object does not contain all the attributes we want
         # see https://github.com/mechmotum/cyipopt/issues/17
         # attr[s.SOLVE_TIME] = solution.solve_time
-        attr[s.NUM_ITERS] = solution['iterations']
+        #attr[s.NUM_ITERS] = solution['iterations']
         # more detailed statistics here when available
         # attr[s.EXTRA_STATS] = solution.extra.FOO
         if 'all_objs_from_best_of' in solution:
@@ -150,13 +152,11 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         use_hessian = (hessian_approx == 'exact')
 
         if solver_cache is None:
-            oracles = Oracles(bounds.new_problem, bounds.x0, len(bounds.cl),
-                            verbose=verbose, use_hessian=use_hessian)
+            oracles = Oracles(bounds.new_problem, verbose=verbose, use_hessian=use_hessian)
         elif 'oracles' in solver_cache:
             oracles = solver_cache['oracles']
         else:
-            oracles = Oracles(bounds.new_problem, bounds.x0, len(bounds.cl),
-                            verbose=verbose, use_hessian=use_hessian)
+            oracles = Oracles(bounds.new_problem, verbose=verbose, use_hessian=use_hessian)
             solver_cache['oracles'] = oracles
 
         nlp = cyipopt.Problem(
@@ -186,16 +186,19 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         for option_name, option_value in default_options.items():
             nlp.add_option(option_name, option_value)
 
+        # ipopt will evaluate the gradient of the Lagrangian at the initial point to decide 
+        # without doing the forward pass for the objective and constraints, so we need to do
+        # a forward pass here to fill in any necessary values for the derivative evaluation.
+        oracles.objective(data["x0"])
+        oracles.constraints(data["x0"])
+    
         _, info = nlp.solve(data["x0"])
 
-        if oracles.iterations == 0 and info['status'] == s.OPTIMAL:
-            print("Warning: IPOPT returned after 0 iterations. This may indicate that\n"
-                  "the initial point passed to Ipopt is a stationary point, and it is\n"
-                  "quite unlikely that the initial point is also a local minimizer. \n"
-                  "Perturb the initial point and try again.")
-
-        # add number of iterations to info dict from oracles
-        info['iterations'] = oracles.iterations
+        # cyipopt does currently not expose the number of iterations, see
+        # https://github.com/mechmotum/cyipopt/issues/17. We set it to "Not available" for now,
+        # but we should update this when the information becomes available.
+        info['num_iters'] = "Not available"
+        
         return info
 
     def cite(self, data):

cvxpy/reductions/solvers/nlp_solvers/knitro_nlpif.py

@@ -177,13 +177,11 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         use_hessian = (hessopt == 1)
 
         if solver_cache is None:
-            oracles = Oracles(bounds.new_problem, bounds.x0, len(bounds.cl),
-                            verbose=verbose, use_hessian=use_hessian)
+            oracles = Oracles(bounds.new_problem, verbose=verbose, use_hessian=use_hessian)
         elif 'oracles' in solver_cache:
             oracles = solver_cache['oracles']
         else:
-            oracles = Oracles(bounds.new_problem, bounds.x0, len(bounds.cl),
-                            verbose=verbose, use_hessian=use_hessian)
+            oracles = Oracles(bounds.new_problem, verbose=verbose, use_hessian=use_hessian)
             solver_cache['oracles'] = oracles
 
         # Extract data from the data dictionary

cvxpy/reductions/solvers/nlp_solvers/nlp_solver.py

@@ -174,8 +174,6 @@ class Oracles:
     def __init__(
         self,
         problem: Problem,
-        initial_point: np.ndarray,
-        num_constraints: int,
         verbose: bool = True,
         use_hessian: bool = True,
     ) -> None:
@@ -191,40 +189,24 @@ def __init__(
         if use_hessian:
             self.c_problem.init_hessian_coo_lower_tri()
 
-        self.initial_point = initial_point
-        self.num_constraints = num_constraints
-        self.iterations = 0
-
         # Cached sparsity structures
         self._jac_structure: tuple[np.ndarray, np.ndarray] | None = None
         self._hess_structure: tuple[np.ndarray, np.ndarray] | None = None
-        self.constraints_forward_passed = False
-        self.objective_forward_passed = False
 
     def objective(self, x: np.ndarray) -> float:
         """Returns the scalar value of the objective given x."""
-        self.objective_forward_passed = True
         return self.c_problem.objective_forward(x)
 
     def gradient(self, x: np.ndarray) -> np.ndarray:
         """Returns the gradient of the objective with respect to x."""
-
-        if not self.objective_forward_passed:
-            self.objective(x)
-
         return self.c_problem.gradient()
 
     def constraints(self, x: np.ndarray) -> np.ndarray:
         """Returns the constraint values."""
-        self.constraints_forward_passed = True
         return self.c_problem.constraint_forward(x)
 
     def jacobian(self, x: np.ndarray) -> np.ndarray:
         """Returns the Jacobian values in COO format at the sparsity structure. """
-
-        if not self.constraints_forward_passed:
-            self.constraints(x)
-
         return self.c_problem.eval_jacobian_vals()
 
     def jacobianstructure(self) -> tuple[np.ndarray, np.ndarray]:
@@ -239,46 +221,23 @@ def jacobianstructure(self) -> tuple[np.ndarray, np.ndarray]:
     def hessian(self, x: np.ndarray, duals: np.ndarray, obj_factor: float) -> np.ndarray:
         """Returns the lower triangular Hessian values in COO format. """
         if not self.use_hessian:
-            # Shouldn't be called when using quasi-Newton, but return empty array
-            return np.array([])
-        
-        if not self.objective_forward_passed:
-            self.objective(x)
-        if not self.constraints_forward_passed:
-            self.constraints(x)
-
+            raise ValueError("Hessian oracle called but use_hessian is False. "
+                             "This is a bug and should be reported.")
+       
         return self.c_problem.eval_hessian_vals_coo_lower_tri(obj_factor, duals)
 
     def hessianstructure(self) -> tuple[np.ndarray, np.ndarray]:
         """Returns the COO sparsity structure of the lower part of the Hessian.
            The returned rows are ascending, and within each row the columns are 
            ascending."""
         if not self.use_hessian:
-            # Return empty structure when using quasi-Newton approximation
-            return (np.array([], dtype=np.int32), np.array([], dtype=np.int32))
-
+            # IPOPT calls this function even when hessian_approximation='limited-memory', 
+            # so return empty structure
+            return (np.array([]), np.array([]))  
+         
         if self._hess_structure is not None:
             return self._hess_structure
 
         rows, cols = self.c_problem.get_problem_hessian_sparsity_coo()
         self._hess_structure = (rows, cols)
         return self._hess_structure
-        
-    def intermediate(
-        self,
-        alg_mod: int,
-        iter_count: int,
-        obj_value: float,
-        inf_pr: float,
-        inf_du: float,
-        mu: float,
-        d_norm: float,
-        regularization_size: float,
-        alpha_du: float,
-        alpha_pr: float,
-        ls_trials: int,
-    ) -> None:
-        """Prints information at every Ipopt iteration."""
-        self.iterations = iter_count
-        self.objective_forward_passed = False
-        self.constraints_forward_passed = False