From 7173cab4925053d0590226455b6fcd5a1aa3fc05 Mon Sep 17 00:00:00 2001 From: Lindon Roberts Date: Mon, 9 Sep 2024 12:13:42 +1000 Subject: [PATCH] Implemented projection-constrained initial directions --- pybobyqa/controller.py | 128 ++++++++++++++++++++++++++++++++++++++++- pybobyqa/params.py | 14 +++++ pybobyqa/solver.py | 16 +++++- 3 files changed, 154 insertions(+), 4 deletions(-) diff --git a/pybobyqa/controller.py b/pybobyqa/controller.py index 1761e1a..40e5aeb 100644 --- a/pybobyqa/controller.py +++ b/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(...) diff --git a/pybobyqa/params.py b/pybobyqa/params.py index 506e435..2f7e9f5 100644 --- a/pybobyqa/params.py +++ b/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 diff --git a/pybobyqa/solver.py b/pybobyqa/solver.py index fa76eb2..7c96677 100644 --- a/pybobyqa/solver.py +++ b/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