Basin hopping bounded optimisation

Viewed 37

I don't understand how to incorporate bounds into scipy's basinhopping routine.

I've got the following script based on scipy's tutorial to test basinhopping with bounds in one dimension.

from scipy.optimize import basinhopping
import matplotlib.pyplot as plt

class MyBounds:
    def __init__(self, xmax, xmin ):
        self.xmax = np.array(xmax)
        self.xmin = np.array(xmin)
    def __call__(self, **kwargs):
        x = kwargs["x_new"]
        tmax = bool(np.all(x <= self.xmax))
        tmin = bool(np.all(x >= self.xmin))
        return tmax and tmin

my_bounds = MyBounds([-4], [2])

func = lambda x: np.cos(14.5 * x - 0.3) + (x + 0.2) * x
x0=[1.]

minimizer_kwargs = {"method": "BFGS"}
ret = basinhopping(func, x0, minimizer_kwargs=minimizer_kwargs,
                   niter=2000, accept_test=my_bounds)
print("global minimum: x = %.4f, f(x) = %.4f" % (ret.x, ret.fun))

bounds = np.linspace(-5, 5, 100)
vals = [func(x) for x in bounds]

plt.plot(bounds, vals)

Without the accept_test = my_bounds the output is global minimum: x = -0.1951, f(x) = -1.0009. With this line included, it outputs global minimum: x = 1.0926, f(x) = 0.4259. This doesn't make sense - -0.1951 is included in [-4, 2] so shouldn't it be the output when the bounds are applied?

Also, changing the bounds to [4, 0] outputs global minimum: x = 1.0926, f(x) = 0.4259. This isn't even in the bounds.

What's the right way to handle this in scipy?

1 Answers

Basinhopping invokes Metropolis Hasting to generate step proposals in an iterative fashion. This is stochastic and happens under the hood of the optimizer. However, this is most likely where your answer is hidden. You can uncover the sequence of proposals (and, if they were accepted or not) with the function

def debug(x, f, accepted):
    print("f(x)={}, \tIs accepted: {}".format(round(f, 3), int(accepted)))

and modifying the line

ret = basinhopping(func, x0, minimizer_kwargs=minimizer_kwargs,
                   niter=2000, accept_test=my_bounds, callback=debug)

in your code.

This prints the status at each step of the iterative procedure. You will be able to see where steps are not accepted (keeping the optimizer away from the minimum you expect.)

However, x = -0.1951 is in the interval [-4, 2]. The same is true for x = 1.0926 and [0, 4]. Unfortunately, the myBounds.__init__(self, xmax, xmin) has a counterintuitive order of parameters listing the large one first and the small one last.

Related