[Ready for review] Initialization order (#114)

dance858 · web-flow · commit d147964d6c72 · 2025-12-13T11:22:03.000-05:00
* initialization

* clean up

* sparse matrix variable fix (thanks William)

* huber comment

* removed assertion

* minor

* fixed to WZZ's comments
diff --git a/cvxpy/problems/problem.py b/cvxpy/problems/problem.py
@@ -1246,6 +1246,7 @@ def _solve(self,
 
             # standard solve
             if best_of == 1:
+                self.set_NLP_initial_point()
                 canon_problem, inverse_data = nlp_chain.apply(problem=self)
                 solution = nlp_chain.solver.solve_via_data(canon_problem, warm_start,
                                                             verbose, solver_opts=kwargs)
@@ -1646,52 +1647,95 @@ def unpack_results(self, solution, chain: SolvingChain, inverse_data) -> None:
         self._solver_stats = SolverStats.from_dict(self._solution.attr,
                                          chain.solver.name())
 
+    
+    def set_NLP_initial_point(self) -> dict:
+        """ Constructs an initial point for the optimization problem. If no
+        initial value is specified, look at the bounds. If both lb and ub are 
+        specified, we initialize the variables to be their midpoints. If only
+        one of them is specified, we initialize the variable one unit from 
+        the bound. If none of them is specified, we initialize it to zero.
+        """
+        for var in self.variables():
+            if var.value is not None:
+                continue
+
+            bounds = var.bounds
+            is_nonneg = var.is_nonneg()
+            is_nonpos = var.is_nonpos()
+            
+            if bounds is None:
+                if is_nonneg:
+                    x0 = np.ones(var.shape)
+                elif is_nonpos:
+                    x0 = -np.ones(var.shape)
+                else:
+                    x0 = np.zeros(var.shape)
+
+                var.save_value(x0)
+            else:
+                lb, ub = bounds
+
+                if is_nonneg:
+                    lb = np.maximum(lb, 0)
+                elif is_nonpos:
+                    ub = np.maximum(ub, 0)
+
+                lb_finite = np.isfinite(lb)
+                ub_finite = np.isfinite(ub)
+                # Replace infs with zero for arithmetic
+                lb0 = np.where(lb_finite, lb, 0.0)
+                ub0 = np.where(ub_finite, ub, 0.0)
+                # Midpoint if both finite, one from bound if only one finite, zero if none
+                init = (lb_finite * ub_finite * 0.5 * (lb0 + ub0) +
+                        lb_finite * (~ub_finite) * (lb0 + 1.0) +
+                        (~lb_finite) * ub_finite * (ub0 - 1.0))
+                var.save_value(init)
+
+
     def set_random_NLP_initial_point(self, run) -> dict:
         """ Generates a random initial point for DNLP problems.
         A variable is initialized randomly in the following cases:
-        1. the initial value specified by the user is None and 
-          'sample_bounds' is set for that variable.
+        1. 'sample_bounds' is set for that variable.
         2. the initial value specified by the user is None,
            'sample_bounds' is not set for that variable, but the
            variable has both finite lower and upper bounds.
-    
-           In other words, a variable that has already been 
-           initialized by the user will not be changed, even if 
-           sample_bounds is set for that variable.
         """
 
-        # store user-specified initial values
+        # store user-specified initial values for variables that do 
+        # not have sample bounds assigned
         if run == 0:
             self._user_initials = {}
             for var in self.variables():
-                self._user_initials[var.id] = var.value
-
+                if var.sample_bounds is not None:
+                    self._user_initials[var.id] = None
+                else:
+                    self._user_initials[var.id] = var.value
+        
         for var in self.variables():
             
             # skip variables with user-specified initial value
+            # (note that any variable with sample bounds set will have
+            #  _user_initials[var.id] == None)
             if self._user_initials[var.id] is not None:
-                # reset to user-specified initial value
-                # from last solve
+                # reset to user-specified initial value from last solve
                 var.value = self._user_initials[var.id]
                 continue
             else:
                 # reset to None from last solve
                 var.value = None 
             
             # set sample_bounds to variable bounds if sample_bounds is None
-            # and variable has finite bounds
+            # and variable bounds (possibly infinite) are set
             if var.sample_bounds is None and var.bounds is not None:
-                low, high = var.bounds 
-                if np.all(np.isfinite(low)) and np.all(np.isfinite(high)):
-                   var.sample_bounds = var.bounds
-
+                var.sample_bounds = var.bounds
+                
             # sample initial value if sample_bounds is set
             if var.sample_bounds is not None:
                 low, high = var.sample_bounds
                 if not np.all(np.isfinite(low)) or not np.all(np.isfinite(high)):
                     raise ValueError(
-                        "Variable %s has non-finite sample_bounds %s."
-                        " Cannot generate random initial point."
+                        "Variable %s has non-finite sample_bounds %s. Cannot generate"
+                        " random initial point. Either add sample bounds or set the value. "
                         % (var.name(), var.sample_bounds)
                     )
 
diff --git a/cvxpy/reductions/cvx_attr2constr.py b/cvxpy/reductions/cvx_attr2constr.py
@@ -156,6 +156,8 @@ def apply(self, problem):
                 elif var.attributes['sparsity']:
                     n = len(var.sparse_idx[0])
                     sparse_var = Variable(n, var_id=var.id, **new_attr)
+                    if var.value_sparse is not None:
+                        sparse_var.value = var.value_sparse.data
                     sparse_var.set_variable_of_provenance(var)
                     id2new_var[var.id] = sparse_var
                     row_idx = np.ravel_multi_index(var.sparse_idx, var.shape, order='F')
diff --git a/cvxpy/reductions/dnlp2smooth/canonicalizers/geo_mean_canon.py b/cvxpy/reductions/dnlp2smooth/canonicalizers/geo_mean_canon.py
@@ -30,8 +30,8 @@ def geo_mean_canon(expr, args):
     """
     t = Variable(expr.shape, nonneg=True)
 
-    if args[0].value is not None and args[0].value > MIN_INIT:
-        t.value = expr.numeric(args[0].value)
+    if args[0].value is not None:
+        t.value = np.max((expr.numeric(args[0].value), MIN_INIT))
     else:
         t.value = np.ones(expr.shape)
 
diff --git a/cvxpy/reductions/dnlp2smooth/canonicalizers/huber_canon.py b/cvxpy/reductions/dnlp2smooth/canonicalizers/huber_canon.py
@@ -17,6 +17,8 @@
 # This file is identical to the dcp2cone version, except that it uses the 
 # dnlp2smooth canonicalizer for the power atom. This is necessary.
 
+import numpy as np
+
 from cvxpy.atoms.elementwise.abs import abs
 from cvxpy.atoms.elementwise.power import power
 from cvxpy.expressions.variable import Variable
@@ -31,6 +33,14 @@ def huber_canon(expr, args):
     n = Variable(shape)
     s = Variable(shape)
 
+    if x.value is None:
+        x.value = np.zeros(x.shape)
+
+    # this choice of initial value follows from how the smooth epigraph
+    # form of the huber function is constructed
+    n.value = np.minimum(np.abs(x.value), M.value) * np.sign(x.value)
+    s.value = x.value - n.value
+
     # n**2 + 2*M*|s|
     power_expr = power(n, 2)
     n2, constr_sq = power_canon(power_expr, power_expr.args)
@@ -41,4 +51,5 @@ def huber_canon(expr, args):
     # x == s + n
     constraints = constr_sq + constr_abs
     constraints.append(x == s + n)
+
     return obj, constraints
diff --git a/cvxpy/reductions/dnlp2smooth/canonicalizers/pnorm_canon.py b/cvxpy/reductions/dnlp2smooth/canonicalizers/pnorm_canon.py
@@ -24,6 +24,10 @@ def pnorm_canon(expr, args):
     p = expr.p
     shape = expr.shape
     t = Variable(shape, nonneg=True)
+
+    # expression will always have a value here in DNLP
+    t.value = expr.value
+
     # we canonicalize 2-norm as follows:
     # ||x||_2 <= t  <=>  quad_over_lin(x, t) <= t
     if p == 2:
diff --git a/cvxpy/reductions/dnlp2smooth/canonicalizers/power_canon.py b/cvxpy/reductions/dnlp2smooth/canonicalizers/power_canon.py
@@ -20,6 +20,7 @@
 from cvxpy.expressions.constants import Constant
 from cvxpy.expressions.variable import Variable
 
+MIN_INIT = 1e-4
 
 def power_canon(expr, args):
     x = args[0]
@@ -45,7 +46,7 @@ def power_canon(expr, args):
         t = Variable(shape, nonneg=True)
 
         if x.value is not None:
-            t.value = x.value
+            t.value = np.maximum(x.value, MIN_INIT)
         else:
             t.value = expr.point_in_domain()
 
diff --git a/cvxpy/reductions/eliminate_pwl/canonicalizers/abs_canon.py b/cvxpy/reductions/eliminate_pwl/canonicalizers/abs_canon.py
@@ -14,11 +14,18 @@
 limitations under the License.
 """
 
+
 from cvxpy.expressions.variable import Variable
 
 
 def abs_canon(expr, args):
     x = args[0]
     t = Variable(expr.shape)
     constraints = [t >= x, t >= -x]
+
+    # for DNLP we must initialize the new variable (DNLP guarantees that 
+    # x.value will be set when this function is called)
+    if expr.value is not None:
+        t.value = expr.value
+
     return t, constraints
diff --git a/cvxpy/reductions/eliminate_pwl/canonicalizers/max_canon.py b/cvxpy/reductions/eliminate_pwl/canonicalizers/max_canon.py
@@ -35,4 +35,10 @@ def max_canon(expr, args):
         promoted_t = reshape(t, (x.shape[0], 1), order='F') @ Constant(np.ones((1, x.shape[1])))
 
     constraints = [x <= promoted_t]
+
+    # for DNLP we must initialize the new variable (DNLP guarantees that 
+    # x.value will be set when this function is called)
+    if expr.value is not None:
+        t.value = expr.value
+        
     return t, constraints
diff --git a/cvxpy/reductions/eliminate_pwl/canonicalizers/maximum_canon.py b/cvxpy/reductions/eliminate_pwl/canonicalizers/maximum_canon.py
@@ -21,11 +21,17 @@
 def maximum_canon(expr, args):
     shape = expr.shape
     t = Variable(shape)
-    
+
+    # for DNLP we must initialize the new variable (DNLP guarantees that 
+    # x.value will be set when this function is called)
+    if expr.value is not None:
+        t.value = expr.value
+
     if expr.is_nonneg():
         t = nonneg_wrap(t)
     if expr.is_nonpos():
         t = nonpos_wrap(t)
     
     constraints = [t >= elem for elem in args]
+
     return t, constraints
diff --git a/cvxpy/reductions/eliminate_pwl/canonicalizers/norm_inf_canon.py b/cvxpy/reductions/eliminate_pwl/canonicalizers/norm_inf_canon.py
@@ -34,4 +34,9 @@ def norm_inf_canon(expr, args):
     else:  # shape = (m, 1)
         promoted_t = reshape(t, (x.shape[0], 1), order='F') @ Constant(np.ones((1, x.shape[1])))
 
+    # for DNLP we must initialize the new variable (DNLP guarantees that 
+    # x.value will be set when this function is called)
+    if expr.value is not None:
+        t.value = expr.value
+
     return t, [x <= promoted_t, x + promoted_t >= 0]
diff --git a/cvxpy/reductions/eliminate_pwl/canonicalizers/sum_largest_canon.py b/cvxpy/reductions/eliminate_pwl/canonicalizers/sum_largest_canon.py
@@ -14,6 +14,8 @@
 limitations under the License.
 """
 
+import numpy as np
+
 from cvxpy.atoms.affine.sum import sum
 from cvxpy.expressions.variable import Variable
 
@@ -29,4 +31,17 @@ def sum_largest_canon(expr, args):
     q = Variable()
     obj = sum(t) + k*q
     constraints = [x <= t + q, t >= 0]
+
+    # for DNLP we must initialize the new variable (DNLP guarantees that 
+    # x.value will be set when this function is called). The initialization 
+    # below is motivated by the optimal solution of min sum(t) + kq subject
+    # to x <= t + q, 0 <= t
+    if x.value is not None:
+        sorted_indices = np.argsort(x.value)
+        idx_of_smallest = sorted_indices[:k]
+        idx_of_largest = sorted_indices[-k:]
+        q.value = np.max(x.value[idx_of_smallest])
+        t.value = np.zeros_like(x.value)
+        t.value[idx_of_largest] = x.value[idx_of_largest] - q.value
+        
     return obj, constraints
diff --git a/cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py b/cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py
@@ -169,6 +169,12 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
 
         _, info = nlp.solve(data["x0"])
 
+        if oracles.iterations == 0:
+            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
         return info
diff --git a/cvxpy/reductions/solvers/nlp_solvers/nlp_solver.py b/cvxpy/reductions/solvers/nlp_solvers/nlp_solver.py
@@ -144,39 +144,17 @@ def get_variable_bounds(self):
         self.lb = np.array(var_lower)
         self.ub = np.array(var_upper)
     
+
     def construct_initial_point(self):
-        """
-        Constructs an initial point for the optimization problem.
-        If no initial value is specified, look at the bounds.
-        If both lb and ub are specified, we initialize the
-        variables to be their midpoints. If only one of them
-        is specified, we initialize the variable one unit
-        from the bound. If none of them is specified, we
-        initialize it to zero.
-        """
-        initial_values = []
-        offset = 0
-        lbs = self.lb
-        ubs = self.ub
-        for var in self.problem.variables():
-            if var.value is not None:
-                initial_values.append(np.atleast_1d(var.value).flatten(order='F'))
-            else:
-                lb = lbs[offset:offset + var.size]
-                ub = ubs[offset:offset + var.size]
-                lb_finite = np.isfinite(lb)
-                ub_finite = np.isfinite(ub)
-                # Replace infs with zero for arithmetic
-                lb0 = np.where(lb_finite, lb, 0.0)
-                ub0 = np.where(ub_finite, ub, 0.0)
-                # Midpoint if both finite, one from bound if only one finite, zero if none
-                init = (lb_finite * ub_finite * 0.5 * (lb0 + ub0) +
-                        lb_finite * (~ub_finite) * (lb0 + 1.0) +
-                        (~lb_finite) * ub_finite * (ub0 - 1.0))
-                initial_values.append(init)
-            offset += var.size
-        self.x0 = np.concatenate(initial_values, axis=0)
+        """ Loop through all variables and collect the intial point."""
+        x0 = []
+        for var in self.main_var:
+            if var.value is None:
+                raise ValueError("Variable %s has no value. This is a bug and should be reported."
+                                  % var.name())
 
+            x0.append(np.atleast_1d(var.value).flatten(order='F'))
+        self.x0 = np.concatenate(x0, axis=0)
 
 class Oracles():
     def __init__(self, problem, initial_point, num_constraints):
diff --git a/cvxpy/tests/NLP_tests/test_best_of.py b/cvxpy/tests/NLP_tests/test_best_of.py
diff --git a/cvxpy/tests/NLP_tests/test_matmul.py b/cvxpy/tests/NLP_tests/test_matmul.py
diff --git a/cvxpy/tests/NLP_tests/test_nlp_solvers.py b/cvxpy/tests/NLP_tests/test_nlp_solvers.py
diff --git a/cvxpy/tests/NLP_tests/test_problem.py b/cvxpy/tests/NLP_tests/test_problem.py
