Implemented projection-constrained initial directions

pybobyqa/controller.py

@@ -268,7 +268,121 @@ def initialise_random_directions(self, number_of_samples, num_directions, params
 
         return None
 
-    def trust_region_step(self):
+    def initialise_general_constraints(self, number_of_samples, num_directions, params):
+        # Initialisation when we have general projection constraints, as well as bound constraints
+        # Basic idea: start with bound-respecting coordinate directions, just like initialise_coordinate_directions
+        # Then, use geometry-improving procedures to make all these points feasible wrt projection constraints
+        # Only after this is all done, then evaluate the initial points
+        # (after the initial coordinate directions given, add points to model with nan objective value)
+        if self.do_logging:
+            module_logger.debug("Initialising to satisfy general projection constraints")
+        # self.model already has x0 evaluated (and x0 feasible), so only need to initialise the other points
+        # num_directions = params("growing.ndirs_initial")
+        assert self.model.num_pts <= (self.n() + 1) * (self.n() + 2) // 2, "prelim: must have npt <= (n+1)(n+2)/2"
+        assert 1 <= num_directions < self.model.num_pts, "Initialisation: must have 1 <= ndirs_initial < npt"
+
+        at_lower_boundary = (self.model.sl > -0.01 * self.delta)  # sl = xl - x0, should be -ve, actually < -rhobeg
+        at_upper_boundary = (self.model.su < 0.01 * self.delta)  # su = xu - x0, should be +ve, actually > rhobeg
+
+        # **** STEP 1 ****
+        # Build initial coordinate-based directions
+        xpts_added = np.zeros((num_directions + 1, self.n()))
+        for k in range(1, num_directions + 1):
+            # k = 0 --> base point (xpt = 0)  [ not here]
+            # k = 1, ..., 2n --> coordinate directions [1,...,n and n+1,...,2n]
+            # k = 2n+1, ..., (n+1)(n+2)/2 --> off-diagonal directions
+            if 1 <= k < self.n() + 1:  # first step along coord directions
+                dirn = k - 1  # direction to move in (0,...,n-1)
+                stepa = self.delta if not at_upper_boundary[dirn] else -self.delta
+                stepb = None
+                xpts_added[k, dirn] = stepa
+
+            elif self.n() + 1 <= k < 2 * self.n() + 1:  # second step along coord directions
+                dirn = k - self.n() - 1  # direction to move in (0,...,n-1)
+                stepa = xpts_added[k - self.n(), dirn]
+                stepb = -self.delta
+                if at_lower_boundary[dirn]:
+                    stepb = min(2.0 * self.delta, self.model.su[dirn])  # su = xu - x0, should be +ve
+                if at_upper_boundary[dirn]:
+                    stepb = max(-2.0 * self.delta, self.model.sl[dirn])  # sl = xl - x0, should be -ve
+                xpts_added[k, dirn] = stepb
+
+            else:  # k = 2n+1, ..., (n+1)(n+2)/2
+                # p = (k - 1) % n + 1  # cycles through (1,...,n), starting at 2n+1 --> 1
+                # l = (k - 2 * n - 1) / n + 1  # (1,...,1, 2, ..., 2, etc.) where each number appears n times
+                # q = (p + l if p + l <= n else p + l - n)
+                stepa = None
+                stepb = None
+                itemp = (k - self.n() - 1) // self.n()
+                q = k - itemp * self.n() - self.n()
+                p = q + itemp
+                if p > self.n():
+                    p, q = q, p - self.n()  # does swap correctly in Python
+
+                xpts_added[k, p - 1] = xpts_added[p, p - 1]
+                xpts_added[k, q - 1] = xpts_added[q, q - 1]
+
+            # Add this new point to the model with a nan function value for now, since it is likely to be moved
+            # to ensure the projection constraints are respected
+            x = self.model.as_absolute_coordinates(xpts_added[k, :])
+            self.model.change_point(k, x - self.model.xbase, np.nan)  # expect step, not absolute x
+
+            # (no need to switch the k-th and (k-N)-th points, since we can't compare function values at this stage)
+
+        # **** STEP 2 ****
+        # Now go through all interpolation points and check if they are feasible wrt projection constraints
+        # If not, do a geometry-improving step (again with nan objective value)
+        for k in range(1, num_directions + 1):
+            x = self.model.as_absolute_coordinates(xpts_added[k, :])
+            this_x_feasible = True
+            for i, proj in enumerate(self.model.projections):
+                if np.linalg.norm(x - proj(x)) > params("projections.feasible_tol"):
+                    this_x_feasible = False
+                    break  # quit loop, only need to find one bad constraint
+
+            # Calculate Lagrange polynomial and update this point
+            if not this_x_feasible:
+                try:
+                    c, g, H = self.model.lagrange_polynomial(k)  # based at xopt
+                    if np.any(np.isinf(g)) or np.any(np.isnan(g)) or np.any(np.isinf(H)) or np.any(np.isnan(H)):
+                        raise ValueError
+                    xnew = ctrsbox_geometry(self.model.xbase, self.model.xopt(), c, g, H, self.model.sl, self.model.su,
+                                            self.model.projections, self.delta,
+                                            dykstra_max_iters=params("projections.dykstra.max_iters"),
+                                            dykstra_tol=params("projections.feasible_tol"),
+                                            gtol=params("projections.pgd_tol"))
+                except LA.LinAlgError:
+                    exit_info = ExitInformation(EXIT_LINALG_ERROR, "Singular matrix encountered in initialization geometry step")
+                    return exit_info  # didn't fix geometry - return & quit
+                except ValueError:
+                    # A ValueError may be raised if gopt or H have nan/inf values (issue #23)
+                    # Ideally this should be picked up earlier in self.model.lagrange_polynomial(...)
+                    exit_info = ExitInformation(EXIT_LINALG_ERROR, "Error when calculating initialization geometry-improving step")
+                    return exit_info  # didn't fix geometry - return & quit
+
+                # Update point and fill with np.nan again, just for now
+                self.model.change_point(k, xnew, np.nan)  # expect step, not absolute x
+
+        # **** STEP 3 ****
+        # Now we are ready to evaluate our initial points
+        for k in range(1, num_directions + 1):
+            x = self.model.xpt(k, abs_coordinates=True)
+            f_list, num_samples_run, exit_info = self.evaluate_objective(x, number_of_samples, params)
+
+            # Handle exit conditions (f < min obj value or maxfun reached)
+            if exit_info is not None:
+                if num_samples_run > 0:
+                    self.model.save_point(x, np.mean(f_list[:num_samples_run]), num_samples_run, x_in_abs_coords=True)
+                return exit_info  # didn't fix geometry - return & quit
+
+            # Otherwise, add new results
+            self.model.change_point(k, x - self.model.xbase, f_list[0])  # expect step, not absolute x
+            for i in range(1, num_samples_run):
+                self.model.add_new_sample(k, f_extra=f_list[i])
+
+        return None   # return & continue
+
+    def trust_region_step(self, params):
         # Build model for full least squares objectives
         gopt, H = self.model.build_full_model()
         if np.any(np.isinf(gopt)) or np.any(np.isnan(gopt)) or np.any(np.isinf(H)) or np.any(np.isnan(H)):
@@ -281,7 +395,11 @@ def trust_region_step(self):
                 if self.model.projections is None:
                     d, gnew, crvmin = trsbox(self.model.xopt(), gopt, H, self.model.sl, self.model.su, self.delta)
                 else:
-                    d, gnew, crvmin = ctrsbox(self.model.xbase, self.model.xopt(), gopt, H, self.model.sl, self.model.su, self.model.projections, self.delta)
+                    d, gnew, crvmin = ctrsbox(self.model.xbase, self.model.xopt(), gopt, H, self.model.sl, self.model.su,
+                                              self.model.projections, self.delta,
+                                              dykstra_max_iters=params("projections.dykstra.max_iters"),
+                                              dykstra_tol=params("projections.feasible_tol"),
+                                              gtol=params("projections.pgd_tol"))
             except ValueError:
                 # A ValueError may be raised if gopt or H have nan/inf values (issue #14)
                 # Although this should be picked up earlier, in this situation just return a zero
@@ -307,7 +425,11 @@ def geometry_step(self, knew, adelt, number_of_samples, params):
             if self.model.projections is None:
                 xnew = trsbox_geometry(self.model.xopt(), c, g, H, self.model.sl, self.model.su, adelt)
             else:
-                xnew = ctrsbox_geometry(self.model.xbase, self.model.xopt(), c, g, H, self.model.sl, self.model.su, self.model.projections, adelt)
+                xnew = ctrsbox_geometry(self.model.xbase, self.model.xopt(), c, g, H, self.model.sl, self.model.su,
+                                        self.model.projections, adelt,
+                                        dykstra_max_iters=params("projections.dykstra.max_iters"),
+                                        dykstra_tol=params("projections.feasible_tol"),
+                                        gtol=params("projections.pgd_tol"))
         except ValueError:
             # A ValueError may be raised if gopt or H have nan/inf values (issue #23)
             # Ideally this should be picked up earlier in self.model.lagrange_polynomial(...)

pybobyqa/params.py

@@ -86,6 +86,12 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False, seek_global_minimum=F
         self.params["restarts.auto_detect.min_chg_model_slope"] = 1.5e-2
         self.params["restarts.auto_detect.min_correl"] = 0.1
 
+        # Projections
+        self.params["projections.dykstra.d_tol"] = 1e-10
+        self.params["projections.dykstra.max_iters"] = 100
+        self.params["projections.feasible_tol"] = 1e-10
+        self.params["projections.pgd_tol"] = 1e-8
+
         self.params_changed = {}
         for p in self.params:
             self.params_changed[p] = False
@@ -193,6 +199,14 @@ def param_type(self, key, npt):
             type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
         elif key == "restarts.auto_detect.min_correl":
             type_str, nonetype_ok, lower, upper = 'float', False, 0.0, 1.0
+        elif key == "projections.dykstra.d_tol":
+            type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
+        elif key == "projections.dykstra.max_iters":
+            type_str, nonetype_ok, lower, upper = 'int', False, 1, None
+        elif key == "projections.feasible_tol":
+            type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
+        elif key == "projections.pgd_tol":
+            type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
         else:
             assert False, "ParameterList.param_type() has unknown key: %s" % key
         return type_str, nonetype_ok, lower, upper

pybobyqa/solver.py

@@ -231,7 +231,7 @@ def solve_main(objfun, x0, args, xl, xu, projections, npt, rhobeg, rhoend, maxfu
 
 
         # Trust region step
-        d, gopt, H, gnew, crvmin = control.trust_region_step()
+        d, gopt, H, gnew, crvmin = control.trust_region_step(params)
         if do_logging:
             module_logger.debug("Trust region step is d = " + str(d))
         xnew = control.model.xopt() + d
@@ -813,6 +813,20 @@ def solve(objfun, x0, args=(), bounds=None, projections=None, npt=None, rhobeg=N
         warnings.warn("x0 above upper bound, adjusting", RuntimeWarning)
     x0[idx] = xu[idx]
 
+    # Enforce feasibility of x0 with respect to projection constraints
+    if projections is not None:
+        x0_feasible = True
+        for i, proj in enumerate(projections):
+            if np.linalg.norm(x0 - proj(x0)) > params("projections.feasible_tol"):
+                x0_feasible = False
+                warnings.warn("x0 not feasible with respect to projections[%g], adjusting" % i, RuntimeWarning)
+                break  # quit loop, only need to find one bad constraint
+        if not x0_feasible:
+            proj_box = lambda w: pbox(w, xl, xu)
+            P = list(projections)  # make a copy of the projections list
+            P.append(proj_box)
+            x0 = dykstra(P, x0)
+
     # Call main solver (first time)
     diagnostic_info = DiagnosticInfo()
     nruns = 0