Maximum function evaluation with NLOPT in Python

Question

I am having an issue with the implementation of NLOPT in Python. My objective is to minimize a somewhat complicated Maximum Likelihood function.

My function is called mle and there are 6 parameters to estimate.

Finding the gradient to this MLE is not trivial, so I decided to turn to a numerical gradient function:

def numgrad(f, x, step=1e-6):
    """numgrad(f: function, x: num array, step: num) -> num array

    Numerically estimates the gradient of a function f which takes an array as
    its argument.
    """
    ary = len(x)

    curr = x * sp.ones((ary, ary))
    next = curr + sp.identity(ary) * step

    delta = sp.apply_along_axis(f, 1, next) - sp.apply_along_axis(f, 1, curr)

    return delta / step

Then my implementation of NLOPT goes like this:

def myfunc(x, grad):
    if grad.size > 0:
        grad = numgrad(mle, [x[0], x[1], x[2], x[3], x[4], x[5]], step=1e-14)
    return mle([x[0], x[1], x[2], x[3], x[4], x[5]])

opt = nlopt.opt(nlopt.LD_SLSQP, 6)
opt.set_lower_bounds([mmin, smin, ming, bmin, vmin, pmin]) #min bound for each of the param.
opt.set_upper_bounds([mmax, smax, maxg, bmax, vmax, pmax])
opt.set_min_objective(myfunc)
opt.set_xtol_rel(1e-15)
opt.maxeval=10000
x = opt.optimize([x1, x2, x3, x4, x5, x6])
minf = opt.last_optimum_value()
print "optimum at ", x[0], x[1], x[2], x[3], x[4], x[5]
print "minimum value = ", minf
print "result code = ", opt.last_optimize_result()

Now the issue is this .... the minimization process goes wayyy tooo fast. In matlab, it takes approx 1 hour and here in Python 12 seconds ... I don't get the same results in Matlab using fmincon.

My feeling is that the code does not recognize the opt.set_xtol_rel(1e-15) and opt.maxeval=10000 statements because even if I increase the number ... no change in the time process...

Or the problem is elsewhere... what am I doing wrong?

Answer 1

You should essentially never estimate the gradient numerically.

You say the gradient is difficult to estimate. In general, if it is at all possible to get the gradient exactly you should do so, and use an appropriate algorithm (NLopt has several, the one you're using should be fine).

However, if you cannot get the gradient exactly, NLopt features several derivative free algorithm, which you could use instead and expect to get better results. I think this is probably the easiest solution to your problem, and would give you better results.

The speed difference can easily be real however. Differences in optimization algorithm and the fact that python is generally faster than Matlab could explain the difference easily.

TLDR: Change from nlopt.LD_SLSQP to nlopt.LN_BOBYQA .

Hope this helps.

Answer 2

To pass gradients from Python to NLopt, you have to

grad[:] = gradient()  # or np.copyto

NOT

grad = gradient()  # wrong, a new pointer

The NLopt Python doc has a big warning on this; unfortunately it's not in the tutorial example you used.

By the way, on random startpoints, see NLopt MLSL
"a sequence of local optimizations (using some other local optimization algorithm) from random startpoints ... clustering heuristic ...". This theme has many many variations.
Also, fwiw, this plot shows 3 NLopt methods from 10 random startpoints on a dozen academic test functions in 8d.