[SciPy-User] Optimization problem using fmin_slsqp

[Qing Yu]

Hi Matt,

Thanks so much! I've actually tried the 'minimize' function in
scipy.optimize. There are several constrained optimization methods in
'minimize' and among them 'SLSQP' seems to be equivalent to the method used
in fmin_slsqp; I got similar bad fitting results using these two functions.
Meanwhile, using the 'L-BFGS-B' method in 'minimize', I can get similar
good fitting results as 'leastsq'. But I'm not sure how these two
algorithms differ.

I then tried the 'lmfit' method you provided, and it seems that this method
is exactly what I want, though I've no idea how it adds bounds to
'leastsq'. Thanks again for your help!

Best,
Qing


On Tue, Apr 8, 2014 at 12:31 AM, Matt Newville
<newville at cars.uchicago.edu>wrote:

> 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
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20140408/5dc43b7d/attachment.html>