[SciPy-User] Optimization problem using fmin_slsqp

[Matt Newville]

Hi Qing,

On Mon, Apr 7, 2014 at 3:56 PM, Qing Yu <yuqingpsy at gmail.com> wrote:
> Hi Matt,
>
> Thanks so much for your help! I've now changed the residual function as
> below:
>
> def residual(a,x,y):
>         return sum((y-f(x,a))**2)
>
> the script seems to be working after the change, but I have another problem
> now:
>
> because I've set the upper and lower bounds of the parameters in fmin_slsqp
> (and that is why I use fmin_slsqp instead of leastsq):
> bounds=[(None,None),(0,pi),(0,None),(None,None)], if I keep the bounds, I
> get the error:
>             Iteration limit exceeded    (Exit mode 9)
>             Current function value: nan
>             Iterations: 101
>             Function evaluations: 1606
>             Gradient evaluations: 101
> I tried changing iter to 10000, the error still existed.
> If I delete the bound option, I can get the successful output. However, the
> fitted values are quite different from what I get using lsqcurvefit in
> matlab. Moreover, if I plot the data, the model fit is indeed worse using
> fmin_slsqp. I'm not sure if the difference is caused by the missing bounds.
> But the outputs from fmin_slsqp did not go out of bounds anyway. Does anyone
> have any idea why I could not set the bounds in fmin_slsqp?
>
> Thanks so much for any help!
>

I don't know enough about the internals of fmin_slsqp() to know for
sure what's going wrong.


You might find the lmfit-py (http://lmfit.github.io/lmfit-py/) package
useful.   This builds on scipy.optimize to allow upper/lower bounds to
be used with any of the scalar minimizers or with least-squares.  It
also allows you to easily change between fitting methods, and several
other bells and whistles.

Translating your problem into lmfit would look something like this:

############
import numpy as np
import pylab

from lmfit import minimize, Parameters, fit_report

pi = np.pi

def f(x, scale, phase, expon, offset):
    return scale*np.exp((expon**2)*np.cos(x-phase)-1) + offset

def resid(params, x, y):
    scale = params['scale'].value
    phase = params['phase'].value
    expon = params['expon'].value
    offset = params['offset'].value
    return f(x, scale, phase, expon, offset) - y

x = np.array([0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6])
y = np.array([0.17,0.36,0.61,0.38,0.17,0.16])

pars = Parameters()
pars.add('scale', value=1, vary=True)
pars.add('phase', value=1, vary=True, min=0, max=pi)
pars.add('expon', value=1, vary=True)
pars.add('offset', value=0, vary=True)

# first fit with the Nelder-Mead method (fmin, scalar)
result = minimize(resid, pars, args=(x,y), method='Nelder')

# then use these values as initial guess to least--squares
# which provides estimate of uncertainties
result = minimize(resid, pars, args=(x,y), method='leastsq')

print fit_report(result.params)

fit = y + result.residual
pylab.plot(x, y, 'r+')
pylab.plot(x, fit, 'b-')
pylab.show()

############

which would print out the results:

[Variables]]
     expon:      2.371419 +/- 0.07606183 (3.21%) initial =  2.371424
     offset:     0.1487316 +/- 0.008012265 (5.39%) initial =  0.1487296
     phase:      1.062886 +/- 0.01148847 (1.08%) initial =  1.062866
     scale:      0.004551185 +/- 0.001651263 (36.28%) initial =  0.004551088
[[Correlations]] (unreported correlations are <  0.100)
    C(expon, scale)              = -0.998
    C(offset, scale)             = -0.719
    C(expon, offset)             =  0.684


(initial values here are those from fmin(), and leastsq() agrees with
the results).   Note that parameters can have bounds (specified by
min, max), and that one can use fmin() or leastsq() without modifying
the objective --   for scalar minimizers it automatically converts
arrays returned by the objective function to a sum of squares.

Hope that helps,

--Matt