How can I write bounds of parameters by using basinhopping?

Question

I have difficulty with writing the bounds of parameters in basinhopping.

(x0)=(a, b, c )

a = (0, 100)

b = (0, 0.100)

c = (0, 10)

from scipy.optimize import basinhopping

minimizer_kwargs = { "method": "Nelder-Mead" }

min = basinhopping(rmse, x0, minimizer_kwargs=minimizer_kwargs, T=0.5, stepsize=0.1, niter=200, disp=True)

Answer 1

There are multiple approaches for this, each potentially behaving differently (common in non-convex global-optimization). The best approach always takes a-priori information about the optimization-problem into consideration!

The most robust general approach (and in my opinion the best) would be a combination of:

• A: inner level: bound-constrained local-search
• B: outer level: bound-constrained steps

The original author of this optimizer says, relying only on A (as done in both other answers as of now) might fail!

Code:

import numpy as np
from scipy.optimize import basinhopping


""" Example problem
    https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.optimize.basinhopping.html
"""
def func2d(x):
    f = np.cos(14.5 * x[0] - 0.3) + (x[1] + 0.2) * x[1] + (x[0] + 0.2) * x[0]
    df = np.zeros(2)
    df[0] = -14.5 * np.sin(14.5 * x[0] - 0.3) + 2. * x[0] + 0.2
    df[1] = 2. * x[1] + 0.2
    return f, df

""" Example bounds """
bx0 = (-0.175, 1.)
bx1 = (-0.09, 1.)
bounds = [bx0, bx1]

""" Solve without bounds """
minimizer_kwargs = {"method":"L-BFGS-B", "jac":True}
x0 = [1.0, 1.0]
ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs, niter=200)
print(ret.message)
print("unconstrained minimum: x = [%.4f, %.4f], f(x0) = %.4f" % (ret.x[0], ret.x[1],ret.fun))


""" Custom step-function """
class RandomDisplacementBounds(object):
    """random displacement with bounds:  see: https://stackoverflow.com/a/21967888/2320035
        Modified! (dropped acceptance-rejection sampling for a more specialized approach)
    """
    def __init__(self, xmin, xmax, stepsize=0.5):
        self.xmin = xmin
        self.xmax = xmax
        self.stepsize = stepsize

    def __call__(self, x):
        """take a random step but ensure the new position is within the bounds """
        min_step = np.maximum(self.xmin - x, -self.stepsize)
        max_step = np.minimum(self.xmax - x, self.stepsize)

        random_step = np.random.uniform(low=min_step, high=max_step, size=x.shape)
        xnew = x + random_step

        return xnew

bounded_step = RandomDisplacementBounds(np.array([b[0] for b in bounds]), np.array([b[1] for b in bounds]))

""" Custom optimizer """
minimizer_kwargs = {"method":"L-BFGS-B", "jac":True, "bounds": bounds}

""" Solve with bounds """
x0 = [1.0, 1.0]
ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs, niter=200, take_step=bounded_step)
print(ret.message)
print("constrained minimum: x = [%.4f, %.4f], f(x0) = %.4f" % (ret.x[0], ret.x[1],ret.fun))

Output:

['requested number of basinhopping iterations completed successfully']
unconstrained minimum: x = [-0.1951, -0.1000], f(x0) = -1.0109
['requested number of basinhopping iterations completed successfully']
constrained minimum: x = [-0.1750, -0.0900], f(x0) = -0.9684

Answer 2

You should use "bounds" parameter in minimizer_kwargs which is passing to scipy.optimize.minimize(). Here is an example:

bnds = ((1, 100), (1, 100), (1,100))# your bounds

def rmse(X):# your function
    a,b,c = X[0],X[1],X[2]
    return a**2+b**2+c**2

x0 = [10., 10., 10.]
minimizer_kwargs = { "method": "L-BFGS-B","bounds":bnds }
ret = basinhopping(rmse, x0, minimizer_kwargs=minimizer_kwargs,niter=10)
print("global minimum: a = %.4f, b = %.4f c = %.4f | f(x0) = %.4f" % (ret.x[0], ret.x[1], ret.x[2], ret.fun))

The result is

global minimum: a = 1.0000, b = 1.0000 c = 1.0000 | f(x0) = 3.0000

Without boundaries it is (clear) 0,0,0

Answer 3

I know this question is old (the accepted answer is from about 5 years ago) and I wish I could comment to it rather than provide a new "answer" (I lack the reputation to do that, unfortunately), but I would like to ask for some clarification on the accepted answer provided by @sascha on Nov 1, 2017.

Specifically, by looking at the code of the function __call__ of class RandomDisplacementBounds(object), I don't understand why x_new is defined as xnew = x + random_step rather than simply as xnew = random_step.

Probably, I am missing out on something and/or misunderstanding the stated problem and context. My reasoning is that if xmin and xmax passed to the function __init__ are indeed the lower and upper bounds of the vector x and random_step is a point randomly drawn from a uniform distribution between xmin and xmax, then what is called random_step is indeed the new point for the local optimization routine. Isn't this true? I would greatly appreciate any help with this.