Using Py-BOBYQA
This section describes the main interface to Py-BOBYQA and how to use it.
Nonlinear Minimization
Py-BOBYQA is designed to solve the local optimization problem
where the bound constraints \(a \leq x \leq b\) are optional. The upper and lower bounds on the variables are non-relaxable (i.e. Py-BOBYQA will never ask to evaluate a point outside the bounds). The objective function \(f(x)\) is usually nonlinear and nonquadratic. If you know your objective is linear or quadratic, you should consider a solver designed for such functions (see here for details).
Py-BOBYQA iteratively constructs an interpolation-based model for the objective, and determines a step using a trust-region framework. For an in-depth technical description of the algorithm see the paper [CFMR2018], and for the global optimization heuristic, see [CRO2022].
How to use Py-BOBYQA
The main interface to Py-BOBYQA is via the function solve
soln = pybobyqa.solve(objfun, x0)
The input objfun is a Python function which takes an input \(x\in\mathbb{R}^n\) and returns the objective value \(f(x)\in\mathbb{R}\). The input of objfun must be one-dimensional NumPy arrays (i.e. with x.shape == (n,)) and the output must be a single Float.
The input x0 is the starting point for the solver, and (where possible) should be set to be the best available estimate of the true solution \(x_{min}\in\mathbb{R}^n\). It should be specified as a one-dimensional NumPy array (i.e. with x0.shape == (n,)).
As Py-BOBYQA is a local solver, providing different values for x0 may cause it to return different solutions, with possibly different objective values.
The output of pybobyqa.solve is an object containing:
soln.x- an estimate of the solution, \(x_{min}\in\mathbb{R}^n\), a one-dimensional NumPy array.soln.f- the objective value at the calculated solution, \(f(x_{min})\), a Float.soln.gradient- an estimate of the gradient vector of first derivatives of the objective, \(g_i \approx \partial f(x_{min})/\partial x_i\), a NumPy array of length \(n\).soln.hessian- an estimate of the Hessian matrix of second derivatives of the objective, \(H_{i,j} \approx \partial^2 f(x_{min})/\partial x_i \partial x_j\), a NumPy array of size \(n\times n\).soln.nf- the number of evaluations ofobjfunthat the algorithm needed, an Integer.soln.nx- the number of points \(x\) at whichobjfunwas evaluated, an Integer. This may be different tosoln.nfif sample averaging is used.soln.nruns- the number of runs performed by Py-BOBYQA (more than 1 if using multiple restarts), an Integer.soln.flag- an exit flag, which can take one of several values (listed below), an Integer.soln.msg- a description of why the algorithm finished, a String.soln.diagnostic_info- a table of diagnostic information showing the progress of the solver, a Pandas DataFrame.
The possible values of soln.flag are defined by the following variables:
soln.EXIT_SUCCESS- Py-BOBYQA terminated successfully (the objective value or trust region radius are sufficiently small).soln.EXIT_MAXFUN_WARNING- maximum allowed objective evaluations reached. This is the most likely return value when using multiple restarts.soln.EXIT_SLOW_WARNING- maximum number of slow iterations reached.soln.EXIT_FALSE_SUCCESS_WARNING- Py-BOBYQA reached the maximum number of restarts which decreased the objective, but to a worse value than was found in a previous run.soln.EXIT_INPUT_ERROR- error in the inputs.soln.EXIT_TR_INCREASE_ERROR- error occurred when solving the trust region subproblem.soln.EXIT_LINALG_ERROR- linear algebra error, e.g. the interpolation points produced a singular linear system.
These variables are defined in the soln object, so can be accessed with, for example
if soln.flag == soln.EXIT_SUCCESS: print("Success!")
Optional Arguments
The solve function has several optional arguments which the user may provide:
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, do_logging=True, print_progress=False)
These arguments are:
args- a tuple of extra arguments passed to the objective function.bounds- a tuple(lower, upper)with the vectors \(a\) and \(b\) of lower and upper bounds on \(x\) (default is \(a_i=-10^{20}\) and \(b_i=10^{20}\)). To set bounds for eitherlowerorupper, but not both, pass a tuple(lower, None)or(None, upper).npt- the number of interpolation points to use (default is \(2n+1\) for a problem withlen(x0)=nifobjfun_has_noise=False, otherwise it is set to \((n+1)(n+2)/2\)). Py-BOBYQA requires \(n+1 \leq npt \leq (n+1)(n+2)/2\). Larger values are particularly useful for noisy problems.rhobeg- the initial value of the trust region radius (default is 0.1 ifscaling_within_bounds=True, otherwise \(0.1\max(\|x_0\|_{\infty}, 1)\)).rhoend- minimum allowed value of trust region radius, which determines when a successful termination occurs (default is \(10^{-8}\)).maxfun- the maximum number of objective evaluations the algorithm may request (default is \(\min(100(n+1),1000)\)).nsamples- a Python functionnsamples(delta, rho, iter, nrestarts)which returns the number of times to evaluateobjfunat 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 (delta), the lower bound on the trust region radius (rho), how many iterations the algorithm has been running for (iter), and how many restarts have been performed (nrestarts). Default is no averaging (i.e.nsamples(delta, rho, iter, nrestarts)=1).user_params- a Python dictionary{'param1': val1, 'param2':val2, ...}of optional parameters. A full list of available options is given in the next section Advanced Usage.objfun_has_noise- a flag to indicate whether or notobjfunhas stochastic noise; i.e. will callingobjfun(x)multiple times at the same value ofxgive 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 inuser_params.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, all of which can be overridden by the values provided inuser_params. IfTrue, both upper and lower bounds must be set. Note that Py-BOBYQA only implements a heuristic method, so there are no guarantees it will find a global minimum. However, by using this flag, it is more likely to escape local minima if there are better values nearby. The method used is a multiple restart mechanism, where we repeatedly re-initialize Py-BOBYQA from the best point found so far, but where we use a larger trust reigon radius each time (note: this is different to more common multi-start approach to global optimization).scaling_within_bounds- a flag to indicate whether the algorithm should internally shift and scale the entries ofxso that the bounds become \(0 \leq x \leq 1\). This is useful is you are settingboundsand the bounds have different orders of magnitude. Ifscaling_within_bounds=True, the values ofrhobegandrhoendapply to the shifted variables.do_logging- a flag to indicate whether logging output should be produced. This is not automatically visible unless you use the Python logging module (see below for simple usage).print_progress- a flag to indicate whether to print a per-iteration progress log to terminal.
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.
The scaling_within_bounds flag is designed to provide an easy way to achieve this, if you have set the bounds lower and upper.
A Simple Example
Suppose we wish to minimize the Rosenbrock test function:
This function has exactly one local minimum \(f(x_{min})=0\) at \(x_{min}=(1,1)\). A commonly-used starting point for testing purposes is \(x_0=(-1.2,1)\). The following script shows how to solve this problem using Py-BOBYQA:
# Py-BOBYQA example: minimize the Rosenbrock function from __future__ import print_function import numpy as np import pybobyqa # Define the objective function def rosenbrock(x): return 100.0 * (x[1] - x[0] ** 2) ** 2 + (1.0 - x[0]) ** 2 # Define the starting point x0 = np.array([-1.2, 1.0]) # Call Py-BOBYQA soln = pybobyqa.solve(rosenbrock, x0) # Display output print(soln)
Note that Py-BOBYQA is a randomized algorithm: in its first phase, it builds an internal approximation to the objective function by sampling it along random directions. In the code above, we set NumPy’s random seed for reproducibility over multiple runs, but this is not required. The output of this script, showing that Py-BOBYQA finds the correct solution, is
****** Py-BOBYQA Results ****** Solution xmin = [1. 1.] Objective value f(xmin) = 1.013856052e-20 Needed 151 objective evaluations (at 151 points) Approximate gradient = [ 2.35772499e-08 -1.07598803e-08] Approximate Hessian = [[ 802.00799968 -400.04089119] [-400.04089119 199.99228723]] Exit flag = 0 Success: rho has reached rhoend ******************************
This and all following problems can be found in the examples directory on the Py-BOBYQA Github page.
Adding Bounds and More Output
We can extend the above script to add constraints. To do this, we can add the lines
# Define bound constraints (lower <= x <= upper) lower = np.array([-10.0, -10.0]) upper = np.array([0.9, 0.85]) # Call Py-BOBYQA (with bounds) soln = pybobyqa.solve(rosenbrock, x0, bounds=(lower,upper))
Py-BOBYQA correctly finds the solution to the constrained problem:
****** Py-BOBYQA Results ****** Solution xmin = [0.9 0.81] Objective value f(xmin) = 0.01 Needed 146 objective evaluations (at 146 points) Approximate gradient = [-2.00000006e-01 -4.74578563e-09] Approximate Hessian = [[ 649.66398033 -361.03094781] [-361.03094781 199.94223213]] Exit flag = 0 Success: rho has reached rhoend ******************************
However, we also get a warning that our starting point was outside of the bounds:
RuntimeWarning: x0 above upper bound, adjusting
Py-BOBYQA automatically fixes this, and moves \(x_0\) to a point within the bounds, in this case \(x_0=(-1.2,0.85)\).
We can also get Py-BOBYQA to print out more detailed information about its progress using the logging module. To do this, we need to add the following lines:
import logging logging.basicConfig(level=logging.INFO, format='%(message)s') # ... (call pybobyqa.solve)
And we can now see each evaluation of objfun:
Function eval 1 at point 1 has f = 39.65 at x = [-1.2 0.85] Initialising (coordinate directions) Function eval 2 at point 2 has f = 14.337296 at x = [-1.08 0.85] Function eval 3 at point 3 has f = 55.25 at x = [-1.2 0.73] ... Function eval 145 at point 145 has f = 0.0100000013172792 at x = [0.89999999 0.80999999] Function eval 146 at point 146 has f = 0.00999999999999993 at x = [0.9 0.81] Did a total of 1 run(s)
If we wanted to save this output to a file, we could replace the above call to logging.basicConfig() with
logging.basicConfig(filename="myfile.log", level=logging.INFO, format='%(message)s', filemode='w')
If you have logging for some parts of your code and you want to deactivate all Py-BOBYQA logging, you can use the optional argument do_logging=False in pybobyqa.solve().
An alternative option available is to get Py-BOBYQA to print to terminal progress information every iteration, by setting the optional argument print_progress=True in pybobyqa.solve(). If we do this for the above example, we get
Run Iter Obj Grad Delta rho Evals 1 1 1.43e+01 1.74e+02 1.20e-01 1.20e-01 5 1 2 5.57e+00 1.20e+02 3.66e-01 1.20e-01 6 1 3 5.57e+00 1.20e+02 6.00e-02 1.20e-02 6 ... 1 132 1.00e-02 2.00e-01 1.50e-08 1.00e-08 144 1 133 1.00e-02 2.00e-01 1.50e-08 1.00e-08 145
Example: Noisy Objective Evaluation
As described in Overview, derivative-free algorithms such as Py-BOBYQA are particularly useful when objfun has noise. Let’s modify the previous example to include random noise in our objective evaluation, and compare it to a derivative-based solver:
# Py-BOBYQA example: minimize the noisy Rosenbrock function from __future__ import print_function import numpy as np import pybobyqa # Define the objective function def rosenbrock(x): return 100.0 * (x[1] - x[0] ** 2) ** 2 + (1.0 - x[0]) ** 2 # Modified objective function: add 1% Gaussian noise def rosenbrock_noisy(x): return rosenbrock(x) * (1.0 + 1e-2 * np.random.normal(size=(1,))[0]) # Define the starting point x0 = np.array([-1.2, 1.0]) # Set random seed (for reproducibility) np.random.seed(0) print("Demonstrate noise in function evaluation:") for i in range(5): print("objfun(x0) = %g" % rosenbrock_noisy(x0)) print("") # Call Py-BOBYQA soln = pybobyqa.solve(rosenbrock_noisy, x0) # Display output print(soln) # Compare with a derivative-based solver import scipy.optimize as opt soln = opt.minimize(rosenbrock_noisy, x0) print("") print("** SciPy results **") print("Solution xmin = %s" % str(soln.x)) print("Objective value f(xmin) = %.10g" % (soln.fun)) print("Needed %g objective evaluations" % soln.nfev) print("Exit flag = %g" % soln.status) print(soln.message)
The output of this is:
Demonstrate noise in function evaluation: objfun(x0) = 24.6269 objfun(x0) = 24.2968 objfun(x0) = 24.4369 objfun(x0) = 24.7423 objfun(x0) = 24.6519 ****** Py-BOBYQA Results ****** Solution xmin = [-1.04327395 1.09935385] Objective value f(xmin) = 4.080030471 Needed 42 objective evaluations (at 42 points) Approximate gradient = [-3786376.5065785 5876675.51763198] Approximate Hessian = [[ 1.32881117e+14 -2.68241358e+14] [-2.68241358e+14 6.09785319e+14]] Exit flag = 0 Success: rho has reached rhoend ****************************** ** SciPy results ** Solution xmin = [-1.20013817 0.99992915] Objective value f(xmin) = 23.86371763 Needed 80 objective evaluations Exit flag = 2 Desired error not necessarily achieved due to precision loss.
Although Py-BOBYQA does not find the true solution (and it cannot produce a good estimate of the objective gradient and Hessian), it still gives a reasonable decrease in the objective. However SciPy’s derivative-based solver, which has no trouble solving the noise-free problem, is unable to make any progress.
As noted above, Py-BOBYQA has an input parameter objfun_has_noise to indicate if objfun has noise in it, which it does in this case. Therefore we can call Py-BOBYQA with
soln = pybobyqa.solve(rosenbrock_noisy, x0, objfun_has_noise=True)
This time, we find the true solution, and better estimates of the gradient and Hessian:
****** Py-BOBYQA Results ****** Solution xmin = [1. 1.] Objective value f(xmin) = 1.237351799e-19 Needed 300 objective evaluations (at 300 points) Did a total of 5 runs Approximate gradient = [-2.17072738e-07 9.62304351e-08] Approximate Hessian = [[ 809.56521044 -400.33737779] [-400.33737779 198.36487985]] Exit flag = 1 Warning (max evals): Objective has been called MAXFUN times ******************************
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.
Note that Py-BOBYQA only implements a heuristic method, so there are no guarantees it will find a global minimum. However, by using the seek_global_minimum flag, it is more likely to escape local minima if there are better values nearby.
# 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:
First run - search for local minimum only ****** Py-BOBYQA Results ****** Solution xmin = [11.41277902 -0.89680525] Objective value f(xmin) = 48.98425368 Needed 143 objective evaluations (at 143 points) Approximate gradient = [-1.64941396e-07 9.69795781e-07] Approximate Hessian = [[ 7.74717421 -104.51102613] [-104.51102613 1135.76500421]] 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) = 3.659891409e-17 Needed 500 objective evaluations (at 500 points) Did a total of 5 runs Approximate gradient = [ 8.70038835e-10 -4.64918043e-07] Approximate Hessian = [[ 4.28883646 64.16836253] [ 64.16836253 3722.93109385]] Exit flag = 1 Warning (max evals): Objective has been called MAXFUN times ******************************
As we can see, the seek_global_minimum flag helped Py-BOBYQA escape the local minimum from the first run, and find the global minimum. More details are given in [CRO2022].
References
Coralia Cartis, Jan Fiala, Benjamin Marteau and Lindon Roberts, Improving the Flexibility and Robustness of Model-Based Derivative-Free Optimization Solvers, ACM Transactions on Mathematical Software, 45:3 (2019), pp. 32:1-32:41 [preprint]
Coralia Cartis, Lindon Roberts and Oliver Sheridan-Methven, Escaping local minima with derivative-free methods: a numerical investigation, Optimization, 71:8 (2022), pp. 2343-2373. [arXiv preprint: 1812.11343]