1
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I will try your suggestion right away. Thank you $\endgroup$ – Plug4 Aug 16 '14 at 23:19
  • $\begingroup$ It turns out that I have the best solution with the following algo: Nelder-Mead Simplex but yet my minimization of my MLE results are not as good as the one that I get in MATLAB ... beside opt.set_xtol_rel(.) and opt.maxeval= what are the other parameters that I can set to make the algo search longer for better results ...... ? $\endgroup$ – Plug4 Aug 16 '14 at 23:41
  • $\begingroup$ AHHH... but wait ... maybe this was not the best initial random variables to begin with... I have to try multiple initial random variables. $\endgroup$ – Plug4 Aug 16 '14 at 23:41
  • $\begingroup$ Yep that was the issue. Thanks for the help! $\endgroup$ – Plug4 Aug 17 '14 at 0:30
  • $\begingroup$ @LKlevin: For this problem, I agree with you. If the number of decision variables is say, less than 10-30, and the gradient cannot be calculated exactly, a derivative-free algorithm is a good choice. For many more variables than that, it is worth attempting to estimate carefully the gradient, if possible. $\endgroup$ – Geoff Oxberry Oct 3 '14 at 20:54
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.