Merge branch 'hessian' of github.com:cvxgrp/cvxpy-ipopt into hessian

Author: dance858

cvxpy/reductions/expr2smooth/canonicalizers/abs_canon.py

@@ -16,13 +16,9 @@
 
 import numpy as np
 
-
-from cvxpy.expressions.variable import Variable
-from cvxpy.atoms.elementwise.power import power
-
 # TODO (DCED): ask William if this the multiplication we want to use
 from cvxpy.atoms.affine.binary_operators import multiply
-from cvxpy.reductions.expr2smooth.canonicalizers.power_canon import power_canon
+from cvxpy.expressions.variable import Variable
 
 #def abs_canon(expr, args):
 #    shape = expr.shape
@@ -41,13 +37,13 @@ def abs_canon(expr, args):
     t1 = Variable(shape, bounds = [0, None])
     y = Variable(shape, bounds = [-1.01, 1.01])
     if args[0].value is not None:
-        #t1.value = np.sqrt(expr.value**2)
         t1.value = np.abs(args[0].value)
         y.value = np.sign(args[0].value)
 
     t1.value = np.ones(shape)
     y.value = np.zeros(shape)
 
-    # TODO (DCED): check how multiply is canonicalized. We don't want to introduce a new variable for
-    # y inside multiply. But args[0] should potentially be canonicalized further?
+    # TODO (DCED): check how multiply is canonicalized. We don't want to introduce a 
+    # new variable for y inside multiply. But args[0] should potentially be canonicalized 
+    # further?
     return t1, [y ** 2 == np.ones(shape), t1 == multiply(y, args[0])]

cvxpy/reductions/expr2smooth/canonicalizers/entr_canon.py

@@ -14,9 +14,11 @@
 limitations under the License.
 """
 
-from cvxpy.expressions.variable import Variable
 import numpy as np
 
+from cvxpy.expressions.variable import Variable
+
+
 def entr_canon(expr, args):
     t = Variable(args[0].shape, bounds=[0, None])
     if args[0].value is not None and np.all(args[0].value >= 1):

cvxpy/reductions/expr2smooth/canonicalizers/kl_div_canon.py

@@ -14,9 +14,11 @@
 limitations under the License.
 """
 
-from cvxpy.expressions.variable import Variable
 import numpy as np
 
+from cvxpy.expressions.variable import Variable
+
+
 def kl_div_canon(expr, args):
     constraints = []
 

cvxpy/reductions/expr2smooth/canonicalizers/log_canon.py

@@ -14,10 +14,11 @@
 limitations under the License.
 """
 
-from cvxpy.expressions.variable import Variable
-from cvxpy.expressions.constants import Constant
 import numpy as np
-from cvxpy.atoms.elementwise.exp import exp
+
+from cvxpy.expressions.constants import Constant
+from cvxpy.expressions.variable import Variable
+
 
 def collect_constant_and_variable(expr, constants, variable):
     if isinstance(expr, Constant):

cvxpy/reductions/expr2smooth/canonicalizers/quad_over_lin.py

@@ -14,7 +14,6 @@
 limitations under the License.
 """
 
-from cvxpy.expressions.variable import Variable
 
 def quad_over_lin_canon(expr, args):
     assert(False)

cvxpy/reductions/expr2smooth/canonicalizers/rel_entr_canon.py

@@ -14,9 +14,11 @@
 limitations under the License.
 """
 
-from cvxpy.expressions.variable import Variable
 import numpy as np
 
+from cvxpy.expressions.variable import Variable
+
+
 def rel_entr_canon(expr, args):
     constraints = []
 

cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py

@@ -93,7 +93,7 @@ def invert(self, solution, inverse_data):
         # 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
 
@@ -135,28 +135,44 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         import cyipopt
         bounds = self.Bounds(data["problem"])
         x0 = self.construct_initial_point(bounds)
+        # Create oracles object
+        oracles = self.Oracles(bounds.new_problem, x0)
         nlp = cyipopt.Problem(
         n=len(x0),
         m=len(bounds.cl),
-        problem_obj=self.Oracles(bounds.new_problem, x0),
+        problem_obj=oracles,
         lb=bounds.lb,
         ub=bounds.ub,
         cl=bounds.cl,
         cu=bounds.cu,
         )
-        nlp.add_option('mu_strategy', 'adaptive')
-        nlp.add_option('tol', 1e-7)
+        # Set default IPOPT options, but use solver_opts if provided
+        default_options = {
+            'mu_strategy': 'adaptive',
+            'tol': 1e-7,
+            'bound_relax_factor': 0.0,
+            'hessian_approximation': 'limited-memory',
+            'derivative_test': 'first-order',
+            'least_square_init_duals': 'yes'
+        }
+        
         #nlp.add_option('honor_original_bounds', 'yes')
-        nlp.add_option('bound_relax_factor', 0.0)
-        nlp.add_option('hessian_approximation', "exact")
-        nlp.add_option('derivative_test', 'second-order')
-        #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)
+
+        # Update defaults with user-provided options
+        if solver_opts:
+            default_options.update(solver_opts)
+        if not verbose:
+            default_options['print_level'] = 3
+        # Apply all options to the nlp object
+        for option_name, option_value in default_options.items():
+            nlp.add_option(option_name, option_value)
+
         _, info = nlp.solve(x0)
+        # add number of iterations to info dict from oracles
+        info['iterations'] = oracles.iterations
         return info
 
     def cite(self, data):
@@ -396,6 +412,12 @@ def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu,
             """Prints information at every Ipopt iteration."""
             self.iterations = iter_count
 
+        def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu,
+                        d_norm, regularization_size, alpha_du, alpha_pr,
+                        ls_trials):
+            """Prints information at every Ipopt iteration."""
+            self.iterations = iter_count
+
     class Bounds():
         def __init__(self, problem):
             self.problem = problem

cvxpy/sandbox/bound_investigation_problems/ML_Gaussian_mean.py

@@ -0,0 +1,50 @@
+import pdb
+from math import pi
+
+import numpy as np
+import numpy.linalg as LA
+
+import cvxpy as cp
+from cvxpy import log, square
+
+np.random.seed(1234)
+n = 2
+data = 10 * np.random.randn(n)
+sigma_opt = (1 / np.sqrt(n)) * LA.norm(data - np.mean(data))
+mu_opt = np.mean(data)
+mu = cp.Variable((1, ), name="mu")
+TO_RUN = 2
+
+# for n = 200, the first one doesn't work if we use cp.sum(cp.square(data-mu)) but it works 
+# with sum_of_squares
+
+# how is the prod canoncalized? maybe that's the issue. Hmm, or the chain rule! Start in the
+# opt solution and see what happens
+
+if TO_RUN == 1:
+    # here we wont induce that sigma is nonnegative so it can be useful to mention it
+    sigma = cp.Variable((1, ), nonneg=True, name="sigma")
+    obj = (n / 2) * log(2*pi*square(sigma)) + (1 / (2 * square(sigma))) * cp.sum(cp.square(data-mu))
+    constraints = []
+elif TO_RUN == 2:
+    # here we will induce that sigma2 is nonnegative so no need to mention it
+    sigma2 = cp.Variable((1, ), name="sigma2")
+    obj = (n / 2) * log(2*pi*sigma2) + (1 / (2 * sigma2)) * cp.sum(cp.square(data-mu))
+    constraints = []
+    sigma = cp.sqrt(sigma2)
+elif TO_RUN == 3:
+    sigma = cp.Variable((1, ))
+    #sigma.value = np.array([1 * sigma_opt])
+    #mu.value = np.array([1 * mu_opt])
+    #t = cp.Variable((n, ))
+    #v = cp.Variable((1, ), bounds=[0, None])
+    obj = n  * log(np.sqrt(2*pi)*sigma) + (1 / (2 * square(sigma))) * cp.sum(cp.square(data-mu))
+    #obj = n  * log(np.sqrt(2*pi)*sigma) + (1 / (2 * square(sigma))) * cp.sum_squares(data-mu)
+    constraints = []
+
+problem = cp.Problem(cp.Minimize(obj), constraints)
+problem.solve(solver=cp.IPOPT, nlp=True)
+
+print("mu difference: ", mu.value - np.mean(data))
+print("sigma difference: ", sigma.value - sigma_opt)
+pdb.set_trace()

cvxpy/sandbox/bound_investigation_problems/ML_Gaussian_mean_stress_test.py

@@ -0,0 +1,64 @@
+from math import pi
+
+import numpy as np
+import numpy.linalg as LA
+
+import cvxpy as cp
+from cvxpy import log, square
+
+np.random.seed(1234)
+TOL = 1e-3
+METHODS = [1, 2, 3, 4, 5]
+all_n = np.arange(2, 100, 5)
+scaling_factors = [1]
+mu = cp.Variable((1, ), name="mu")
+
+for n in all_n:
+    for factor in scaling_factors:
+        data = factor*np.random.randn(n)
+        sigma_opt = (1 / np.sqrt(n)) * LA.norm(data - np.mean(data))
+        mu_opt = np.mean(data)
+        for method in METHODS:
+            mu.value = None 
+            print("Method, n, scale factor: ", method, n, factor)
+            if method == 1:
+                # here we wont deduce that sigma is nonnegative so it can be useful to mention it
+                sigma = cp.Variable((1, ), nonneg=True)
+                obj = (n / 2) * log(2*pi*square(sigma)) + \
+                      (1 / (2 * square(sigma))) * cp.sum(cp.square(data-mu))
+                constraints = []
+            elif method == 2:
+                # here we will deduce that sigma2 is nonnegative so no need to mention it
+                sigma2 = cp.Variable((1, ), name="Sigma2")
+                obj = (n / 2) * log( 2 * pi * sigma2) + \
+                      (1 / (2 * sigma2)) * cp.sum(cp.square(data-mu))
+                constraints = []
+                sigma = cp.sqrt(sigma2)
+            elif method == 3:
+                # here we will deduce that sigma is nonnegative so no need to mention it
+                sigma = cp.Variable((1, ), name="Sigma")
+                obj = n  * log(np.sqrt(2*pi)*sigma) + \
+                      (1 / (2 * square(sigma))) * cp.sum(cp.square(data-mu))
+                constraints = []
+            elif method == 4:
+                # here we will deduce that sigma is nonnegative so no need to mention it
+                sigma2 = cp.Variable((1, ), name="Sigma2")
+                obj = (n / 2) * log(sigma2 * 2 * pi * -1 * -1) + \
+                      (1 / (2 * sigma2)) * cp.sum(cp.square(data-mu))
+                constraints = []
+                sigma = cp.sqrt(sigma2)
+            elif method == 5:
+                # here we will deduce that sigma is nonnegative so no need to mention it
+                sigma = cp.Variable((1, ), name="Sigma")
+                obj = n  * log(np.sqrt(2*pi)*sigma * -1 * -1 * 2 * 0.5) + \
+                      (1 / (2 * square(sigma))) * cp.sum(cp.square(data-mu))
+                constraints = []
+
+            problem = cp.Problem(cp.Minimize(obj), constraints)
+            problem.solve(solver=cp.IPOPT, nlp=True)
+            print("sigma.value: ", sigma.value)
+            print("sigma_opt: ", sigma_opt)
+            assert(np.abs(sigma.value - sigma_opt) / np.max([1, np.abs(sigma_opt)]) <= TOL)
+            assert(np.abs(mu.value - mu_opt) / np.max([1, np.abs(mu_opt)]) <= TOL)
+
+

cvxpy/sandbox/bound_investigation_problems/ML_Gaussian_zero_mean.py

@@ -0,0 +1,34 @@
+from math import pi
+
+import numpy as np
+import numpy.linalg as LA
+
+import cvxpy as cp
+from cvxpy import log, square
+
+np.random.seed(1234)
+n = 10
+data = np.random.randn(n)
+sigma_opt = (1 / np.sqrt(n)) * LA.norm(data)
+res = LA.norm(data) ** 2
+
+TO_RUN = 2
+
+if TO_RUN == 1:
+    sigma = cp.Variable((1, ))
+    obj = (n / 2) * log(2*pi*square(sigma)) + (1 / (2 * square(sigma))) * res
+    constraints = []
+elif TO_RUN == 2:
+    sigma2 = cp.Variable((1, ))
+    obj = (n / 2) * log(2*pi*sigma2) + (1 / (2 * sigma2)) * res
+    constraints = []
+    sigma = cp.sqrt(sigma2)
+elif TO_RUN == 3:
+    sigma = cp.Variable((1, ))
+    obj = n  * log(np.sqrt(2*pi)*sigma) + (1 / (2 * square(sigma))) * res
+    constraints = []
+
+problem = cp.Problem(cp.Minimize(obj), constraints)
+problem.solve(solver=cp.IPOPT, nlp=True)
+print("difference sigma:", sigma.value - sigma_opt)
+