how to find global minimum in python optimization with bounds?

Question

I have a Python function with 64 variables, and I tried to optimise it using L-BFGS-B method in the minimise function, however this method have quite a strong dependence on the initial guess, and failed to find the global minimum.

But I liked its ability to set bounds for the variables. Is there a way/function to find the global minimum while having boundaries for the variables ?

Answer 1

This can be done with scipy.optimize.basinhopping. Basinhopping is a function designed to find the global minimum of an objective function. It does repeated minimizations using the function scipy.optimize.minimize and takes a random step in coordinate space after each minimization. Basinhopping can still respect bounds by using one of the minimizers that implement bounds (e.g. L-BFGS-B). Here is some code that shows how to do this

# an example function with multiple minima
def f(x): return x.dot(x) + sin(np.linalg.norm(x) * np.pi)

# the starting point
x0 = [10., 10.]

# the bounds
xmin = [1., 1.]
xmax = [11., 11.]

# rewrite the bounds in the way required by L-BFGS-B
bounds = [(low, high) for low, high in zip(xmin, xmax)]

# use method L-BFGS-B because the problem is smooth and bounded
minimizer_kwargs = dict(method="L-BFGS-B", bounds=bounds)
res = basinhopping(f, x0, minimizer_kwargs=minimizer_kwargs)
print res

The above code will work for a simple case, but you can still end up in a forbidden region if basinhopping random displacement routine takes you there. Luckily that can be overridden by passing a custom step taking routine using the keyword take_step

class RandomDisplacementBounds(object):
    """random displacement with bounds"""
    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"""
        while True:
            # this could be done in a much more clever way, but it will work for example purposes
            xnew = x + np.random.uniform(-self.stepsize, self.stepsize, np.shape(x))
            if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
                break
        return xnew

# define the new step taking routine and pass it to basinhopping
take_step = RandomDisplacementBounds(xmin, xmax)
result = basinhopping(f, x0, niter=100, minimizer_kwargs=minimizer_kwargs,
                      take_step=take_step)
print result

Answer 2

Some common-sense suggestions for debugging and visualizing any optimizer on your function:

Are your objective function and your constraints reasonable ?
If the objective function is a sum say f() + g(), print those separately for all the x in "fx-opt.nptxt" (below); if f() is 99 % of the sum and g() 1 %, investigate.

Constraints: how many of the components x_i in xfinal are stuck at bounds, x_i <= lo_i or >= hi_i ?

---

How bumpy is your function on a global scale ?
Run with several random startpoints, and save the results to analyze / plot:

title = "%s  n %d  ntermhess %d  nsample %d  seed %d" % (  # all params!
    __file__, n, ntermhess, nsample, seed )
print title
...
np.random.seed(seed)  # for reproducible runs
np.set_printoptions( threshold=100, edgeitems=10, linewidth=100,
        formatter = dict( float = lambda x: "%.3g" % x ))  # float arrays %.3g

lo, hi = bounds.T  # vecs of numbers or +- np.inf
print "lo:", lo
print "hi:", hi

fx = []  # accumulate all the final f, x
for jsample in range(nsample):
        # x0 uniformly random in box lo .. hi --
    x0 = lo + np.random.uniform( size=n ) * (hi - lo)

    x, f, d = fmin_l_bfgs_b( func, x0, approx_grad=1,
                m=ntermhess, factr=factr, pgtol=pgtol )
    print "f: %g  x: %s  x0: %s" % (f, x, x0)
    fx.append( np.r_[ f, x ])

fx = np.array(fx)  # nsample rows, 1 + dim cols
np.savetxt( "fx-opt.nptxt", fx, fmt="%8.3g", header=title )  # to analyze / plot

ffinal = fx[:,0]
xfinal = fx[:,1:]
print "final f values, sorted:", np.sort(ffinal)
jbest = ffinal.argmin()
print "best x:", xfinal[jbest]

If some of the ffinal values look reasonably good, try more random startpoints near those -- that's surely better than pure random.

If the x s are curves, or anything real, plot the best few x0 and xfinal.
(A rule of thumb is nsample ~ 5*d or 10*d in d dimensions. Too slow, too many ? Reduce maxiter / maxeval, reduce ftol -- you don't need ftol 1e-6 for exploration like this.)

If you want reproducible results, then you must list ALL relevant parameters in the title and in derived files and plots. Otherwise, you'll be asking "where did this come from ??"

---

How bumpy is your function on epsilon scale ~ 10^-6 ?
Methods that approximate a gradient sometimes return their last estimate, but if not:

from scipy.optimize._numdiff import approx_derivative  # 3-point, much better than
## from scipy.optimize import approx_fprime
for eps in [1e-3, 1e-6]:
    grad = approx_fprime( x, func, epsilon=eps )
    print "approx_fprime eps %g: %s" % (eps, grad)

If however the gradient estimate is poor / bumpy before the optimizer quit, you won't see that. Then you have to save all the intermediate [f, x, approx_fprime] to watch them too; easy in python -- ask if that's not clear.

In some problem areas it's common to back up and restart from a purported xmin. For example, if you're lost on a country road, first find a major road, then restart from there.

---

Summary:
don't expect any black-box optimizer to work on a function that's large-scale bumpy, or epsilon-scale bumpy, or both.
Invest in test scaffolding, and in ways to see what the optimizer is doing.