Merge pull request #25 from numericalalgorithmsgroup/dev_existing_evaluations

Use existing evaluations

Author: lindonroberts

dfols/__init__.py

@@ -39,9 +39,11 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 
 # DFO-LS version
-__version__ = '1.5.4'
+__version__ = '1.6'
 
 # Main solver & exit flags
 from .solver import *
 __all__ = ['solve', 'OptimResults']
 
+from .evaluations_database import *
+__all__ += ['EvaluationDatabase']

dfols/controller.py

@@ -414,6 +414,48 @@ def initialise_random_directions(self, number_of_samples, num_directions, params
 
         return None
 
+    def initialise_from_database(self, eval_database, number_of_samples, params):
+        # Here, eval_database has at least one entry, and the base index has already been used
+        # to evaluate (x0,r0), which has already been added to self.model
+        # Now, find exactly n feasible perturbations (either from database or new evals) and add them to the model
+        base_idx, perturbation_idx, new_perturbations = eval_database.select_starting_evals(self.delta,
+                                                                                            xl=self.model.xbase + self.model.sl,
+                                                                                            xu=self.model.xbase + self.model.su,
+                                                                                            projections=self.model.projections,
+                                                                                            tol=params("database.new_direction_tol"),
+                                                                                            dykstra_max_iters=params("dykstra.max_iters"),
+                                                                                            dykstra_tol=params("dykstra.d_tol"))
+
+        # Add suitable pre-existing evaluations
+        for i, idx in enumerate(perturbation_idx):
+            module_logger.info("Adding pre-existing evaluation %g to initial model" % idx)
+            x, rx = eval_database.get_eval(idx)
+            self.model.change_point(i + 1, x - self.model.xbase, rx, -idx)  # use eval_num = -idx
+
+        if new_perturbations is not None:
+            num_perturbations = new_perturbations.shape[0]
+            module_logger.debug("Adding %g new evaluations to initial model" % num_perturbations)
+            for i in range(num_perturbations):
+                new_point = (eval_database.get_x(base_idx) - self.model.xbase) + new_perturbations[i,:]  # new_perturbations[i,:] has length <= self.delta
+
+                # Evaluate objective
+                x = self.model.as_absolute_coordinates(new_point)
+                rvec_list, obj_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(rvec_list[:num_samples_run, :], axis=0), num_samples_run,
+                                              self.nx, x_in_abs_coords=True)
+                    return exit_info  # return & quit
+
+                # Otherwise, add new results (increments model.npt_so_far)
+                self.model.change_point(len(perturbation_idx) + 1 + i, x - self.model.xbase, rvec_list[0, :], self.nx)  # expect step, not absolute x
+                for j in range(1, num_samples_run):
+                    self.model.add_new_sample(len(perturbation_idx) + 1 + i, rvec_extra=rvec_list[j, :])
+
+        return None
+
     def add_new_direction_while_growing(self, number_of_samples, params, min_num_steps=0):
         num_steps = max(params('growing.num_new_dirns_each_iter'), min_num_steps)
         step_length = params('growing.delta_scale_new_dirns') * self.delta

dfols/evaluations_database.py

@@ -0,0 +1,208 @@
+"""
+Class to create/store database of existing evaluations, and routines to select
+existing evaluations to build an initial linear model
+"""
+import logging
+import numpy as np
+
+from .util import apply_scaling, dykstra
+from .trust_region import ctrsbox_geometry, trsbox_geometry
+
+__all__ = ['EvaluationDatabase']
+
+module_logger = logging.getLogger(__name__)
+
+
+# Class to store set of evaluations (x, rx)
+class EvaluationDatabase(object):
+    def __init__(self, eval_list=None, starting_eval=None):
+        # eval_list is a list of tuples (x, rx)
+        self._evals = []
+        if eval_list is not None:
+            for e in eval_list:
+                self._evals.append(e)
+
+        # Which evaluation index should be the starting point of the optimization?
+        self.starting_eval = None
+        if starting_eval is not None and 0 <= starting_eval <= len(self._evals):
+            self.starting_eval = starting_eval
+
+    def __len__(self):
+        return len(self._evals)
+
+    def append(self, x, rx, make_starting_eval=False):
+        self._evals.append((x, rx))
+        if make_starting_eval:
+            self.starting_eval = len(self) - 1
+
+    def set_starting_eval(self, index):
+        if 0 <= index < len(self):
+            self.starting_eval = index
+        else:
+            raise IndexError("Invalid index %g given current set of %g evaluations" % (index, len(self)))
+
+    def get_starting_eval_idx(self):
+        if len(self) == 0:
+            raise RuntimeError("No evaluations available, no suitable starting evaluation ")
+        elif self.starting_eval is None:
+            module_logger.warning("Starting evaluation index not set, using most recently appended evaluation")
+            self.starting_eval = len(self) - 1
+
+        return self.starting_eval
+
+    def get_eval(self, index):
+        # Return (x, rx) for given index
+        if 0 <= index < len(self):
+            return self._evals[index][0], self._evals[index][1]
+        else:
+            raise IndexError("Invalid index %g given current set of %g evaluations" % (index, len(self)))
+
+    def get_x(self, index):
+        return self.get_eval(index)[0]
+
+    def get_rx(self, index):
+        return self.get_eval(index)[1]
+
+    def apply_scaling(self, scaling_changes):
+        # Adjust all input x values based on scaling
+        if scaling_changes is not None:
+            for i in range(len(self)):
+                x, rx = self._evals[i]
+                self._evals[i] = (apply_scaling(x, scaling_changes), rx)
+        return
+
+    def select_starting_evals(self, delta, xl=None, xu=None, projections=[], tol=1e-8,
+                              dykstra_max_iters=100, dykstra_tol=1e-10):
+        # Given a database 'evals' with prescribed starting index, and initial trust-region radius delta > 0
+        # determine a subset of the database to use
+
+        # The bounds xl <= x <= xu and projection list are used to determine where to evaluate any new points
+        # (ensuring they are feasible)
+
+        if delta <= 0.0:
+            raise RuntimeError("delta must be strictly positive")
+        if len(self) == 0:
+            raise RuntimeError("Need at least one evaluation to select starting evaluations")
+
+        base_idx = self.get_starting_eval_idx()
+        xbase = self.get_x(self.get_starting_eval_idx())
+        n = len(xbase)
+        module_logger.debug("Selecting starting evaluations from existing database")
+        module_logger.debug("Have %g evaluations to choose from" % len(self))
+        module_logger.debug("Using base index %g" % base_idx)
+
+        # For linear interpolation, we will use the matrix
+        # M = [[1, 0], [0, L]] where L has rows (xi-xbase)/delta
+        # So, just build a large matrix Lfull with everything
+        n_perturbations = len(self) - 1
+        Lfull = np.zeros((n_perturbations, n))
+        row_idx = 0
+        for i in range(n_perturbations + 1):
+            if i == base_idx:
+                continue
+            Lfull[row_idx, :] = (self.get_x(i) - xbase) / delta  # Lfull[i,:] = (xi-xbase) / delta
+            row_idx += 1
+
+        xdist = np.linalg.norm(Lfull, axis=1)  # xdist[i] = ||Lfull[i,:]|| = ||xi-xbase|| / delta
+        # module_logger.debug("xdist =", xdist)
+
+        # We ideally want xdist ~ 1, so reweight these distances based on that (large xdist_reweighted --> xdist ~ 1 --> good)
+        xdist_reweighted = 1.0 / np.maximum(xdist, 1.0 / xdist)
+        # module_logger.debug("xdist_reweighted =", xdist_reweighted)
+
+        if n_perturbations == 0:
+            module_logger.debug("Only one evaluation available, just selecting that")
+            return base_idx, [], delta * np.eye(n)
+
+        # Now, find as many good perturbations as we can
+        # Good = not too far from xbase (relative to delta) and sufficiently linearly independent
+        #        from other selected perturbations (i.e. Lfull[perturbation_idx,:] well-conditioned
+        #        and len(perturbation_idx) <= n
+        perturbation_idx = []  # what point indices to use as perturbations
+
+        for iter in range(min(n_perturbations, n)):
+            # Add one more good perturbation, if available
+            # Note: can only add at most the number of available perturbations, or n perturbations, whichever is smaller
+            if iter == 0:
+                # First perturbation: every direction is equally good, so pick the point closest to the
+                # trust-region boundary
+                idx = int(np.argmax(xdist_reweighted))
+                module_logger.debug("Adding index %g with ||xi-xbase|| / delta = %g" % (idx if idx < base_idx else idx+1, xdist[idx]))
+                perturbation_idx.append(idx)
+            else:
+                Q, R = np.linalg.qr(Lfull[perturbation_idx, :].T, mode='reduced')
+                # module_logger.debug("Current perturbation_idx =", perturbation_idx)
+                L_rem = Lfull @ (np.eye(n) - Q @ Q.T)  # part of (xi-xbase)/delta orthogonal to current perturbations
+                # rem_size = fraction of original length ||xi-xbase||/delta that is orthogonal to current perturbations
+                # all entries are in [0,1], and is zero for already selected perturbations
+                rem_size = np.linalg.norm(L_rem, axis=1) / xdist
+                rem_size[perturbation_idx] = 0  # ensure this holds exactly
+                # module_logger.debug("rem_size =", rem_size)
+                # module_logger.debug("rem_size * xdist_reweighted =", rem_size * xdist_reweighted)
+
+                # We want a point with large rem_size and xdist ~ 1 (i.e. xdist_reweighted large)
+                idx = int(np.argmax(rem_size * xdist_reweighted))
+                if rem_size[idx] * xdist_reweighted[idx] > tol:
+                    # This ensures new perturbation is sufficiently linearly independent of existing perturbations
+                    # (and also ensures idx hasn't already been chosen)
+                    module_logger.debug("Adding index %g" % (idx if idx < base_idx else idx+1))
+                    perturbation_idx.append(idx)
+                else:
+                    module_logger.debug("No more linearly independent directions, quitting")
+                    break
+
+        # Find new linearly independent directions
+        if len(perturbation_idx) < n:
+            module_logger.debug("Selecting %g new linearly independent directions" % (n - len(perturbation_idx)))
+            Q, _ = np.linalg.qr(Lfull[perturbation_idx, :].T, mode='complete')
+            new_perturbations = delta * Q[:, len(perturbation_idx):].T
+
+            # Make perturbations feasible w.r.t. xl <= x <= xu and projections
+            # Note: if len(projections) > 0, then the projection list *already* includes bounds
+            # Don't need to make pre-existing evaluations feasible, since we already have r(x) for these
+
+            # Start construction of interpolation matrix for later
+            L = np.zeros((n, n), dtype=float)
+            L[:len(perturbation_idx), :] = Lfull[perturbation_idx, :]
+            L[len(perturbation_idx):, :] = new_perturbations / delta
+
+            # Since we already have a full set of linearly independent directions,
+            # we do this by moving each infeasible perturbation to a geometry-improving location
+            for i in range(new_perturbations.shape[0]):
+                xnew = xbase + new_perturbations[i, :]
+                # Check feasibility
+                if len(projections) == 0:
+                    # Bounds only
+                    feasible = np.all(xnew >= xl) and np.all(xnew <= xu)
+                else:
+                    # Projections
+                    xnew_C = dykstra(projections, xnew, max_iter=dykstra_max_iters, tol=dykstra_tol)
+                    feasible = np.linalg.norm(xnew - xnew_C) < dykstra_tol
+
+                if feasible:
+                    # Skip feasible points, nothing to do
+                    continue
+
+                # If infeasible, build Lagrange polynomial and move to geometry-improving location in B(xbase,delta)
+                # which will automatically be feasible
+                module_logger.debug("Moving default %g-th new perturbation to ensure feasibility" % i)
+                c = 0.0  # Lagrange polynomial centered at xbase
+                ei = np.zeros((n,), dtype=float)
+                ei[len(perturbation_idx) + i] = 1.0
+                g = np.linalg.solve(L, ei) / delta  # divide by delta because L is scaled by 1/delta
+                if len(projections) == 0:
+                    new_perturbations[i, :] = trsbox_geometry(xbase, c, g, xl, xu, delta)
+                else:
+                    new_perturbations[i, :] = ctrsbox_geometry(xbase, c, g, projections, delta)
+
+                # Update L after replacement
+                L[len(perturbation_idx) + i, :] = new_perturbations[i,:] / delta
+        else:
+            module_logger.debug("Full set of directions found, no need for new evaluations")
+            new_perturbations = None
+
+        # perturbation_idx in [0, ..., n_perturbations-1], reset to be actual indices
+        for i in range(len(perturbation_idx)):
+            if perturbation_idx[i] >= base_idx:
+                perturbation_idx[i] += 1
+        return base_idx, perturbation_idx, new_perturbations

dfols/params.py

@@ -122,6 +122,9 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
         self.params["func_tol.tr_step"] = 1-1e-1
         self.params["func_tol.max_iters"] = 500
         self.params["sfista.max_iters_scaling"] = 2.0
+
+        # Evaluation database
+        self.params["database.new_direction_tol"] = 1e-8
 
         self.params_changed = {}
         for p in self.params:
@@ -284,6 +287,8 @@ def param_type(self, key, npt):
             type_str, nonetype_ok, lower, upper = 'int', False, 0, None
         elif key == "sfista.max_iters_scaling":
             type_str, nonetype_ok, lower, upper = 'float', False, 1.0, None
+        elif key == "database.new_direction_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