Bounds investigation (#29)

* div canon stuff

* some test problems in sandbox

* preparing for merge

* some progress

* weird value check

* stress test ML zero mean passes

* very subtle implementation details. MLE stress test (with and without zero mean) work

* minor

* fixed dimension error inside div-canon and added Sharpe ratio as a test

* fixing some linting issues

* more fixing imports

* fixing final imports, giving up

---------

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

Author: dance858

cvxpy/atoms/elementwise/power.py

@@ -395,6 +395,11 @@ def _domain(self) -> List[Constraint]:
             return [self.args[0] >= 0]
         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/expr2smooth/canonicalizers/__init__.py

@@ -35,5 +35,5 @@
     log: log_canon,
     power: power_canon,
     Pnorm : pnorm_canon,
-    DivExpression: div_canon,       
+    DivExpression: div_canon
 }

cvxpy/reductions/expr2smooth/canonicalizers/div_canon.py

@@ -14,23 +14,33 @@
 limitations under the License.
 """
 
+import numpy as np
+
 from cvxpy.expressions.variable import Variable
 
+
 # We canonicalize div(x, y) as z * y = x.
 def div_canon(expr, args):
-    # TODO: potential bounds here?
-    # TODO: wrong dimension for z here
-    z = Variable(args[1].shape) 
+    dim = (1, ) if args[0].shape == () else args[0].shape
+    sgn_z = args[0].sign
+
+    if sgn_z == 'NONNEGATIVE':
+        z = Variable(dim, bounds=[0, None])
+    elif sgn_z == 'NONPOSITIVE':
+        z = Variable(dim, bounds=[None, 0])
+    else:
+        z = Variable(dim)
+    
     y = Variable(args[1].shape, bounds=[0, None])
 
-    if args[1].value is not None and args[1].value != 0.0:
+    # DCED: perhaps initialize further away from boundary?
+    if args[1].value is not None and np.all(args[1].value != 0.0):
         y.value = args[1].value
     else:
         y.value = expr.point_in_domain()
 
     if args[0].value is not None:
-        z.value = args[0].value / y.value 
+        z.value = np.atleast_1d(args[0].value) / y.value 
     else:
         z.value = expr.point_in_domain()
-    # TODO: should we also include y >= 0 here?
-    return z, [z * y == args[0], y == args[1]]#], #y >= 0, z >= 0]
+    return z, [z * y == args[0], y == args[1]]

cvxpy/reductions/expr2smooth/canonicalizers/log_canon.py

@@ -14,13 +14,60 @@
 limitations under the License.
 """
 
+import numpy as np
+
+from cvxpy.expressions.constants import Constant
 from cvxpy.expressions.variable import Variable
 
 
+def collect_constant_and_variable(expr, constants, variable):
+    if isinstance(expr, Constant):
+        constants.append(expr)
+    elif isinstance(expr, Variable):
+        variable.append(expr)
+    elif hasattr(expr, "args"):
+        for subexpr in expr.args:
+            collect_constant_and_variable(subexpr, constants, variable)
+
+    assert(len(variable) <= 1)
+
+# DCED: Without this lower bound the stress test for ML Gaussian non-zero mean fails.
+# Perhaps this should be a parameter exposed to the user?
+LOWER_BOUND = 1e-5
+
 def log_canon(expr, args):
-    t = Variable(args[0].size)
-    if args[0].value is not None:
+    t = Variable(args[0].size, bounds=[LOWER_BOUND, None], name='t')
+
+    # DCED: if args[0] is a * x for a constant scalar or vector 'a' 
+    # and a vector variable 'x', we want to add bounds to x if x
+    # does not have any bounds. We also want to initialize x far 
+    # from its bounds. 
+    constants, variable = [], []
+    collect_constant_and_variable(args[0], constants, variable) 
+    a = 1
+    is_special_case = True
+    for constant in constants:
+        if len(constant.shape) == 2:
+            is_special_case = False
+            break
+        else:
+            a *= constant.value
+
+    if not is_special_case:
+        if args[0].value is not None and np.all(args[0].value > 0):
+            t.value = args[0].value
+        else:
+            t.value = expr.point_in_domain()
+    else:  
+        if variable[0].value is None:
+            variable[0].value = expr.point_in_domain() * np.sign(a)
+        
+        lbs = -np.inf * np.ones(args[0].size)
+        ubs = np.inf * np.ones(args[0].size)
+        lbs[a > 0] = 0
+        ubs[a < 0] = 0
+        variable[0].bounds = [lbs, ubs]        
+        assert(args[0].value is not None and np.all(args[0].value > 0.0))
         t.value = args[0].value
-    else:
-        t.value = expr.point_in_domain()
+
     return expr.copy([t]), [t==args[0]]

cvxpy/reductions/expr2smooth/canonicalizers/pnorm_canon.py

@@ -14,8 +14,9 @@
 limitations under the License.
 """
 
-from cvxpy.expressions.variable import Variable
 from cvxpy.atoms.affine.sum import sum
+from cvxpy.expressions.variable import Variable
+
 
 def pnorm_canon(expr, args):
     x = args[0]

cvxpy/reductions/expr2smooth/canonicalizers/power_canon.py

@@ -23,8 +23,7 @@
 def power_canon(expr, args):
     x = args[0]
     p = expr.p_rational
-    #w = expr.w
-
+    
     if p == 1:
         return x, []
 
@@ -35,16 +34,22 @@ def power_canon(expr, args):
     else:
         t = Variable(shape)
         if 0 < p < 1:
-            if x.value is not None:
-                # TODO: check if this initialization is correct
-                t.value = np.power(np.abs(x.value), p)
-            return t, [t**(1/p) == x, t >= 0]
+            raise NotImplementedError('This power is not yet supported.')
+        
+            # DCED: This is not a smooth implementation so we raise an error for now
+            #if x.value is not None:
+            #    t.value = np.power(np.abs(x.value), p)
+            #return t, [t**(1/p) == x, t >= 0]
         elif p > 1:
             t = Variable(args[0].shape)
-            if args[0].value is not None:
+
+            # DCED: for Gaussian ML it works better if we include the second
+            #       condition here
+            if args[0].value is not None and np.all(args[0].value >= 1):
                 t.value = args[0].value
             else:
                 t.value = expr.point_in_domain()
+
             return expr.copy([t]), [t==args[0]]
         else:
-            pass
+            raise NotImplementedError('This power is not yet supported.')

cvxpy/reductions/expr2smooth/expr2smooth.py

@@ -18,7 +18,6 @@
 
 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

cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py

@@ -16,7 +16,6 @@
 
 import numpy as np
 import scipy.sparse as sp
-from time import time
 
 import cvxpy.settings as s
 from cvxpy.constraints import (
@@ -134,14 +133,7 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         """
         import cyipopt
         bounds = self.Bounds(data["problem"])
-        initial_values = []
-        for var in bounds.main_var:
-            if var.value is not None:
-                initial_values.append(var.value.flatten(order='F'))
-            else:
-                # If no initial value, use zero
-                initial_values.append(np.zeros(var.size))
-        x0 = np.concatenate(initial_values, axis=0)
+        x0 = self.construct_initial_point(bounds)
         nlp = cyipopt.Problem(
         n=len(x0),
         m=len(bounds.cl),
@@ -153,9 +145,14 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         )
         nlp.add_option('mu_strategy', 'adaptive')
         nlp.add_option('tol', 1e-7)
-        nlp.add_option('honor_original_bounds', 'yes')
+        #nlp.add_option('honor_original_bounds', 'yes')
+        nlp.add_option('bound_relax_factor', 0.0)
         nlp.add_option('hessian_approximation', "limited-memory")
         nlp.add_option('derivative_test', 'first-order')
+        nlp.add_option('least_square_init_duals', 'yes')
+        #nlp.add_option('constr_mult_init_max', 1e10) 
+        #nlp.add_option('derivative_test_perturbation', 1e-5)
+        #nlp.add_option('point_perturbation_radius', 0.1)
         _, info = nlp.solve(x0)
         return info
 
@@ -168,6 +165,43 @@ def cite(self, data):
             Data generated via an apply call.
         """
         return CITATION_DICT["IPOPT"]
+    
+    # TODO (DCED): Ask WZ where to put this.
+    def construct_initial_point(self, bounds):
+            initial_values = []
+            offset = 0
+            lbs = bounds.lb 
+            ubs = bounds.ub
+            for var in bounds.main_var:
+                if var.value is not None:
+                    initial_values.append(var.value.flatten(order='F'))
+                else:
+                    # 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.
+                    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
+            x0 = np.concatenate(initial_values, axis=0)
+            return x0
 
     class Oracles():
         def __init__(self, problem, inital_point):
@@ -176,7 +210,7 @@ def __init__(self, problem, inital_point):
             self.initial_point = inital_point
             for var in self.problem.variables():
                 self.main_var.append(var)
-        
+
         def objective(self, x):
             """Returns the scalar value of the objective given x."""
             # Set the variable value
@@ -190,6 +224,8 @@ def objective(self, x):
             return obj_value
 
         def gradient(self, x):
+            #import pdb 
+            #pdb.set_trace()
             """Returns the gradient of the objective with respect to x."""
             # compute the gradient using _grad
             offset = 0
@@ -214,6 +250,8 @@ def constraints(self, x):
             """Returns the constraint values."""
             # Set the variable value
             offset = 0
+            #import pdb 
+            #pdb.set_trace()
             for var in self.main_var:
                 size = var.size
                 var.value = x[offset:offset+size].reshape(var.shape, order='F')
@@ -222,12 +260,13 @@ def constraints(self, x):
             # Evaluate all constraints
             constraint_values = []
             for constraint in self.problem.constraints:
-                constraint_values.append(np.asarray(constraint.args[0].value).flatten())
+                constraint_values.append(np.asarray(constraint.args[0].value).flatten(order='F'))
             return np.concatenate(constraint_values)
 
         def jacobian(self, x):
             """Returns only the non-zero values of the Jacobian."""
-            
+            #import pdb 
+            #pdb.set_trace()
             # Set variable values
             offset = 0
             for var in self.main_var:
@@ -245,7 +284,10 @@ def jacobian(self, x):
                         jacobian = grad_dict[var].T
                         if sp.issparse(jacobian):
                             jacobian = sp.dok_matrix(jacobian)
-                            data = np.array([jacobian.get((r, c), 0) for r, c in zip(rows, cols)])        
+                            data = np.array([
+                                jacobian.get((r, c), 0)
+                                for r, c in zip(rows, cols)
+                            ])
                             values.append(np.atleast_1d(data))
                         else:
                             values.append(np.atleast_1d(jacobian))
@@ -328,12 +370,29 @@ def get_variable_bounds(self):
             for var in self.main_var:
                 size = var.size
                 if var.bounds:
-                    var_lower.extend(var.bounds[0].flatten(order='F'))
-                    var_upper.extend(var.bounds[1].flatten(order='F'))
+                    lb = var.bounds[0].flatten(order='F')
+                    ub = var.bounds[1].flatten(order='F')
+
+                    if var.is_nonneg():
+                        lb = np.maximum(lb, 0)
+                    
+                    if var.is_nonpos():
+                        ub = np.minimum(ub, 0)
+                    
+                    var_lower.extend(lb)
+                    var_upper.extend(ub)
                 else:
-                    # No bounds specified, use infinite bounds
-                    var_lower.extend([-np.inf] * size)
-                    var_upper.extend([np.inf] * size)
+                    # No bounds specified, use infinite bounds or bounds
+                    # set by the nonnegative or nonpositive attribute
+                    if var.is_nonneg():
+                        var_lower.extend([0.0] * size)
+                    else:
+                        var_lower.extend([-np.inf] * size)
+                    
+                    if var.is_nonpos():
+                        var_upper.extend([0.0] * size)
+                    else:
+                        var_upper.extend([np.inf] * size)
 
             self.lb = np.array(var_lower)
             self.ub = np.array(var_upper)

cvxpy/sandbox/control_of_car_matrix.py

@@ -1,6 +1,7 @@
-import cvxpy as cp
-import numpy as np
 import matplotlib.pyplot as plt
+import numpy as np
+
+import cvxpy as cp
 
 
 def solve_car_control(x_final, L=0.1, N=2, h=0.1, gamma=10):
@@ -91,7 +92,7 @@ def plot_trajectory(x_opt, u_opt, L, h, title="Car Trajectory"):
 
     # Plot car position and orientation at several time steps
     car_length = L
-    car_width = L * 0.6
+    # car_width = L * 0.6
 
     # Select time steps to show car outline (every 5th step)
     steps_to_show = range(0, len(x_opt), 5)
@@ -153,7 +154,11 @@ def plot_trajectory(x_opt, u_opt, L, h, title="Car Trajectory"):
 
             if x_opt is not None and u_opt is not None:
                 print("Optimization successful!")
-                print(f"Final position: p1={x_opt[-1, 0]:.3f}, p2={x_opt[-1, 1]:.3f}, theta={x_opt[-1, 2]:.3f}")
+                print(
+                    f"Final position: p1={x_opt[-1, 0]:.3f}, "
+                    f"p2={x_opt[-1, 1]:.3f}, "
+                    f"theta={x_opt[-1, 2]:.3f}"
+                )
 
                 # Plot the trajectory
                 fig, ax = plot_trajectory(x_opt, u_opt, L=0.1, h=0.1, title=description)

cvxpy/sandbox/entr_max.py

@@ -1,6 +1,7 @@
-import cvxpy as cp
 import numpy as np
 
+import cvxpy as cp
+
 k = 5
 x = cp.Variable((k, k))
 x.value = np.ones((k, k)) / k**2