New global optimization functionality, bugfix for trust region radius when scaling bounds

Author: lindonroberts

docs/advanced.rst

@@ -4,7 +4,7 @@ This section describes different optional user parameters available in Py-BOBYQA
 
 In the last section (:doc:`userguide`), we introduced :code:`pybobyqa.solve()`, which has the optional input :code:`user_params`. This is a Python dictionary of user parameters. We will now go through the settings which can be changed in this way. More details are available in the paper [CFMR2018]_.
 
-The default values, used if no override is given, in some cases vary depending on whether :code:`objfun` has stochastic noise; that is, whether evaluating :code:`objfun(x)` several times at the same :code:`x` gives the same result or not. Whether or not this is the case is determined by the :code:`objfun_has_noise` input to :code:`pybobyqa.solve()` (and not by inspecting :code:`objfun`, for instance).
+The default values, used if no override is given, in some cases vary depending on whether :code:`objfun` has stochastic noise; that is, whether evaluating :code:`objfun(x)` several times at the same :code:`x` gives the same result or not. Whether or not this is the case is determined by the :code:`objfun_has_noise` input to :code:`pybobyqa.solve()` (and not by inspecting :code:`objfun`, for instance). Similarly, the default values depend on the input flag :code:`seek_global_minimum`, i.e. if a global minimum is desired.
 
 General Algorithm Parameters
 ----------------------------
@@ -59,8 +59,10 @@ Interpolation Management
 
 Multiple Restarts
 -----------------
-* :code:`restarts.use_restarts` - Whether to do restarts when :math:`\rho_k` reaches :math:`\rho_{end}`, or (optionally) when all points are within noise level of :math:`f(x_k)`. Default is :code:`False` for smooth problems or :code:`True` for noisy problems. 
+* :code:`restarts.use_restarts` - Whether to do restarts when :math:`\rho_k` reaches :math:`\rho_{end}`, or (optionally) when all points are within noise level of :math:`f(x_k)`. Default is :code:`False` for smooth problems or :code:`True` for noisy problems or when seeking a global minimum. 
 * :code:`restarts.max_unsuccessful_restarts` - Maximum number of consecutive unsuccessful restarts allowed (i.e.~restarts which did not reduce the objective further). Default is 10. 
+* :code:`restarts.max_unsuccessful_restarts_total` - Maximum number of total unsuccessful restarts allowed. Default is 20 when seeking a global minimum, otherwise it is :code:`maxfun` (i.e.~not restricted).
+* :code:`restarts.rhobeg_scale_after_unsuccessful_restart` - Factor to increase :math:`\rho_{beg}` by after unsuccessful restarts. Default is 1.1 when seeking a global minimum, otherwise it is 1.
 * :code:`restarts.rhoend_scale` - Factor to reduce :math:`\rho_{end}` by with each restart. Default is 1. 
 * :code:`restarts.use_soft_restarts` - Whether to use soft or hard restarts. Default is :code:`True`. 
 * :code:`restarts.soft.num_geom_steps` - For soft restarts, the number of points to move. Default is 3. 

docs/history.rst

@@ -15,3 +15,9 @@ Version 1.0.2 (20 Jun 2018)
 * Extra optional input :code:`args` which passes through arguments for :code:`objfun` (pull request from `logangrado <https://github.com/logangrado>`_).
 * Bug fixes: default parameters for reduced initialization cost regime, returning correct value from safety steps, retrieving dependencies during installation.
 
+Version 1.1a0 (13 Jul 2018)
+---------------------------
+* Extra parameters to control the trust region radius over multiple restarts, designed for global optimization.
+* New input flag :code:`seek_global_minimum` to set sensible default parameters for global optimization. New example script to demonstrate this functionality.
+* Bug fix: default trust region radius when scaling variables within bounds.
+

docs/index.rst

@@ -25,6 +25,8 @@ Full details of the Py-BOBYQA algorithm are given in our paper: C. Cartis, J. Fi
 
 If you are interested in solving least-squares minimization problems, you may wish to try `DFO-LS <https://github.com/numericalalgorithmsgroup/dfols>`_, which has the same features as Py-BOBYQA (plus some more), and exploits the least-squares problem structure, so performs better on such problems.
 
+**New feature: global optimization (July 2018)!** Py-BOBYQA now has a heuristic for global optimization (see :doc:`userguide` for details). As this is a heuristic, there are no guarantees it will find a global minimum, but it is more likely to escape local minima if there are better values nearby.
+
 Py-BOBYQA is released under the GNU General Public License. Please `contact NAG <http://www.nag.com/content/worldwide-contact-information>`_ for alternative licensing.
 
 .. toctree::

docs/userguide.rst

@@ -68,19 +68,21 @@ The :code:`solve` function has several optional arguments which the user may pro
       pybobyqa.solve(objfun, x0, args=(), bounds=None, npt=None,
 		  rhobeg=None, rhoend=1e-8, maxfun=None, nsamples=None, 
                   user_params=None, objfun_has_noise=False, 
+                  seek_global_minimum=False, 
                   scaling_within_bounds=False)
 
 These arguments are:
 
 * :code:`args` - a tuple of extra arguments passed to the objective function.
 * :code:`bounds` - a tuple :code:`(lower, upper)` with the vectors :math:`a` and :math:`b` of lower and upper bounds on :math:`x` (default is :math:`a_i=-10^{20}` and :math:`b_i=10^{20}`). To set bounds for either :code:`lower` or :code:`upper`, but not both, pass a tuple :code:`(lower, None)` or :code:`(None, upper)`.
 * :code:`npt` - the number of interpolation points to use (default is :code:`2*len(x0)+1`). Py-BOBYQA requires :code:`n+1 <= npt <= (n+1)*(n+2)/2` for a problem with :code:`len(x0)=n`. Larger values are particularly useful for noisy problems.
-* :code:`rhobeg` - the initial value of the trust region radius (default is :math:`0.1\max(\|x_0\|_{\infty}, 1)`).
+* :code:`rhobeg` - the initial value of the trust region radius (default is 0.1 if :code:`scaling_within_bounds=True`, otherwise :math:`0.1\max(\|x_0\|_{\infty}, 1)`).
 * :code:`rhoend` - minimum allowed value of trust region radius, which determines when a successful termination occurs (default is :math:`10^{-8}`).
 * :code:`maxfun` - the maximum number of objective evaluations the algorithm may request (default is :math:`\min(100(n+1),1000)`).
 * :code:`nsamples` - a Python function :code:`nsamples(delta, rho, iter, nrestarts)` which returns the number of times to evaluate :code:`objfun` at a given point. This is only applicable for objectives with stochastic noise, when averaging multiple evaluations at the same point produces a more accurate value. The input parameters are the trust region radius (:code:`delta`), the lower bound on the trust region radius (:code:`rho`), how many iterations the algorithm has been running for (:code:`iter`), and how many restarts have been performed (:code:`nrestarts`). Default is no averaging (i.e. :code:`nsamples(delta, rho, iter, nrestarts)=1`).
 * :code:`user_params` - a Python dictionary :code:`{'param1': val1, 'param2':val2, ...}` of optional parameters. A full list of available options is given in the next section :doc:`advanced`.
 * :code:`objfun_has_noise` - a flag to indicate whether or not :code:`objfun` has stochastic noise; i.e. will calling :code:`objfun(x)` multiple times at the same value of :code:`x` give different results? This is used to set some sensible default parameters (including using multiple restarts), all of which can be overridden by the values provided in :code:`user_params`.
+* :code:`seek_global_minimum` - a flag to indicate whether to search for a global minimum, rather than a local minimum. This is used to set some sensible default parameters (including using multiple restarts), all of which can be overridden by the values provided in :code:`user_params`. If :code:`True`, both upper and lower bounds must be set.
 * :code:`scaling_within_bounds` - a flag to indicate whether the algorithm should internally shift and scale the entries of :code:`x` so that the bounds become :math:`0 \leq x \leq 1`. This is useful is you are setting :code:`bounds` and the bounds have different orders of magnitude. If :code:`scaling_within_bounds=True`, the values of :code:`rhobeg` and :code:`rhoend` apply to the *shifted* variables.
 
 In general when using optimization software, it is good practice to scale your variables so that moving each by a given amount has approximately the same impact on the objective function.
@@ -305,6 +307,87 @@ This time, we find the true solution, and better estimates of the gradient and H
       ******************************
 
 
+Example: Global Optimization
+----------------------------
+The following example shows how to use the global optimization features of Py-BOBYQA. Here, we try to minimize the Freudenstein and Roth function (problem 2 in J.J. Moré, B.S. Garbow, B.S. and K.E. Hillstrom, Testing Unconstrained Optimization Software, *ACM Trans. Math. Software* 7:1 (1981), 17-41). This function has two local minima, one of which is global.
+
+  .. code-block:: python
+  
+      # Py-BOBYQA example: globally minimize the Freudenstein and Roth function
+      from __future__ import print_function
+      import numpy as np
+      import pybobyqa
+      
+      # Define the objective function
+      # This function has a local minimum f = 48.98 
+      # at x = np.array([11.41, -0.8968])
+      # and a global minimum f = 0 at x = np.array([5.0, 4.0])
+      def freudenstein_roth(x):
+          r1 = -13.0 + x[0] + ((5.0 - x[1]) * x[1] - 2.0) * x[1]
+          r2 = -29.0 + x[0] + ((1.0 + x[1]) * x[1] - 14.0) * x[1]
+          return r1 ** 2 + r2 ** 2
+      
+      # Define the starting point
+      x0 = np.array([5.0, -20.0])
+      
+      # Define bounds (required for global optimization)
+      lower = np.array([-30.0, -30.0])
+      upper = np.array([30.0, 30.0])
+      
+      # Set random seed (for reproducibility)
+      np.random.seed(0)
+      
+      print("First run - search for local minimum only")
+      print("")
+      soln = pybobyqa.solve(freudenstein_roth, x0, maxfun=500, 
+                            bounds=(lower, upper))
+      print(soln)
+      
+      print("")
+      print("")
+      
+      print("Second run - search for global minimum")
+      print("")
+      soln = pybobyqa.solve(freudenstein_roth, x0, maxfun=500, 
+                            bounds=(lower, upper), 
+                            seek_global_minimum=True)
+      print(soln)
+
+The output of this is:
+
+  .. code-block:: none
+  
+      First run - search for local minimum only
+      
+      ****** Py-BOBYQA Results ******
+      Solution xmin = [ 11.41277906  -0.89680525]
+      Objective value f(xmin) = 48.98425368
+      Needed 203 objective evaluations (at 203 points)
+      Approximate gradient = [ -1.61348180e-06  -3.61662651e-07]
+      Approximate Hessian = [[ 132.10265455  -45.5426821 ]
+       [ -45.5426821   976.15808779]]
+      Exit flag = 0
+      Success: rho has reached rhoend
+      ******************************
+      
+      
+      
+      Second run - search for global minimum
+      
+      ****** Py-BOBYQA Results ******
+      Solution xmin = [ 5.  4.]
+      Objective value f(xmin) = 9.734692105e-19
+      Needed 500 objective evaluations (at 500 points)
+      Did a total of 4 runs
+      Approximate gradient = [  4.28964221e-08   4.58344260e-07]
+      Approximate Hessian = [[    4.06992486    61.15006935]
+       [   61.15006935  3728.06826545]]
+      Exit flag = 1
+      Warning (max evals): Objective has been called MAXFUN times
+      ******************************
+
+As we can see, the :code:`seek_global_minimum` flag helped Py-BOBYQA escape the local minimum from the first run, and find the global minimum.
+
 References
 ----------
 

examples/global_opt_demo.py

@@ -0,0 +1,35 @@
+# Py-BOBYQA example: globally minimize the Freudenstein and Roth function
+from __future__ import print_function
+import numpy as np
+import pybobyqa
+
+# Define the objective function
+# This function has a local minimum f = 48.98 at x = np.array([11.41, -0.8968])
+# and a global minimum f = 0 at x = np.array([5.0, 4.0])
+def freudenstein_roth(x):
+    r1 = -13.0 + x[0] + ((5.0 - x[1]) * x[1] - 2.0) * x[1]
+    r2 = -29.0 + x[0] + ((1.0 + x[1]) * x[1] - 14.0) * x[1]
+    return r1 ** 2 + r2 ** 2
+
+# Define the starting point
+x0 = np.array([5.0, -20.0])
+
+# Define bounds (required for global optimization)
+lower = np.array([-30.0, -30.0])
+upper = np.array([30.0, 30.0])
+
+# Set random seed (for reproducibility)
+np.random.seed(0)
+
+print("First run - search for local minimum only")
+print("")
+soln = pybobyqa.solve(freudenstein_roth, x0, maxfun=500, bounds=(lower, upper))
+print(soln)
+
+print("")
+print("")
+
+print("Second run - search for global minimum")
+print("")
+soln = pybobyqa.solve(freudenstein_roth, x0, maxfun=500, bounds=(lower, upper), seek_global_minimum=True)
+print(soln)

pybobyqa/controller.py

@@ -119,6 +119,7 @@ def __init__(self, objfun, x0, args, f0, f0_nsamples, xl, xu, npt, rhobeg, rhoen
         self.last_successful_iter = 0  # when ||d|| >= rho
         # For counting the number of soft restarts
         self.last_successful_run = 0
+        self.total_unsuccessful_restarts = 0
         self.last_run_fopt = f0
         self.scaling_changes = scaling_changes
 
@@ -465,16 +466,22 @@ def soft_restart(self, number_of_samples, nruns_so_far, params, x_in_abs_coords_
         # A successful run is one where we reduced fopt
         if self.model.fopt() < self.last_run_fopt:
             self.last_successful_run = nruns_so_far
+        else:
+            # Unsuccessful run
+            self.total_unsuccessful_restarts += 1
+            self.rhobeg = self.rhobeg * params("restarts.rhobeg_scale_after_unsuccessful_restart")
         self.last_run_fopt = self.model.fopt()
 
         ok_to_do_restart = (nruns_so_far - self.last_successful_run < params("restarts.max_unsuccessful_restarts")) and \
-                           (self.nf < self.maxfun)
+                           (self.nf < self.maxfun) and self.total_unsuccessful_restarts < params("restarts.max_unsuccessful_restarts_total")
 
         if not ok_to_do_restart:
             # last outputs are (exit_flag, exit_str, return_to_new_tr_iteration)
             exit_info = ExitInformation(EXIT_MAXFUN_WARNING, "Objective has been called MAXFUN times")
             if nruns_so_far - self.last_successful_run >= params("restarts.max_unsuccessful_restarts"):
-                exit_info = ExitInformation(EXIT_SUCCESS, "Reached maximum number of unsuccessful restarts")
+                exit_info = ExitInformation(EXIT_SUCCESS, "Reached maximum number of consecutive unsuccessful restarts")
+            elif self.total_unsuccessful_restarts >= params("restarts.max_unsuccessful_restarts_total"):
+                exit_info = ExitInformation(EXIT_SUCCESS, "Reached maximum total number of unsuccessful restarts")
             return exit_info
 
         # Save current best point and the one the user wanted to save too

pybobyqa/params.py

@@ -32,7 +32,7 @@
 
 
 class ParameterList(object):
-    def __init__(self, n, npt, maxfun, objfun_has_noise=False):
+    def __init__(self, n, npt, maxfun, objfun_has_noise=False, seek_global_minimum=False):
         self.params = {}
         # Rounding error constant (for shifting base)
         self.params["general.rounding_error_constant"] = 0.1  # 0.1 in DFBOLS, 1e-3 in BOBYQA
@@ -71,8 +71,10 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
         self.params["noise.multiplicative_noise_level"] = None
         self.params["noise.additive_noise_level"] = None
         # Restarts
-        self.params["restarts.use_restarts"] = True if objfun_has_noise else False
+        self.params["restarts.use_restarts"] = True if (objfun_has_noise or seek_global_minimum) else False
         self.params["restarts.max_unsuccessful_restarts"] = 10
+        self.params["restarts.max_unsuccessful_restarts_total"] = 20 if seek_global_minimum else maxfun  # i.e. not used (in line with previous behaviour)
+        self.params["restarts.rhobeg_scale_after_unsuccessful_restart"] = 1.1 if seek_global_minimum else 1.0
         self.params["restarts.rhoend_scale"] = 1.0  # how much to decrease rhoend by after each restart
         self.params["restarts.use_soft_restarts"] = True
         self.params["restarts.soft.num_geom_steps"] = 3
@@ -161,6 +163,10 @@ def param_type(self, key, npt):
             type_str, nonetype_ok, lower, upper = 'bool', False, None, None
         elif key == "restarts.max_unsuccessful_restarts":
             type_str, nonetype_ok, lower, upper = 'int', False, 0, None
+        elif key == "restarts.max_unsuccessful_restarts_total":
+            type_str, nonetype_ok, lower, upper = 'int', False, 0, None
+        elif key == "restarts.rhobeg_scale_after_unsuccessful_restart":
+            type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
         elif key == "restarts.rhoend_scale":
             type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
         elif key == "restarts.use_soft_restarts":