[SciPy-User] beginner's question regarding optimize.fmin_l_bfgs_b

Fri Oct 15 05:05:05 EDT 2010

-----Original Message-----
From: scipy-user-bounces at scipy.org [mailto:scipy-user-bounces at scipy.org] On Behalf Of Sebastian Haase
Sent: 15. oktober 2010 10:31
To: SciPy Users List
Subject: Re: [SciPy-User] beginner's question regarding optimize.fmin_l_bfgs_b

On Thu, Oct 14, 2010 at 11:58 PM, Skipper Seabold <jsseabold at gmail.com> wrote:
> On Thu, Oct 14, 2010 at 5:32 PM, Skipper Seabold <jsseabold at gmail.com> wrote:
>> On Thu, Oct 14, 2010 at 5:21 PM, Skipper Seabold <jsseabold at gmail.com> wrote:
>>> On Tue, Oct 12, 2010 at 9:10 AM, Tveraa, Torkild <Torkild.Tveraa at nina.no> wrote:
>>>> Dear All,
>>>>
>>>> I have been able to use the optimize.leastsq - module to minimize a given function (see below), but since my data is sparse I have convergence problems and would ideally be able to put bounds on the parameters. If I have understood this correctly this can be done with the optimize.fmin_l_bfgs_b - module, but I am unable to figure out how to do this. Some helps & hints would be most appreciated :-)
>>>>
>>>>        Cheers,
>>>>        Torkild
>>>>
>>>> -------------------------------------------------------
>>>> import numpy
>>>> import pylab
>>>> from scipy import *
>>>> from scipy import optimize
>>>>
>>>> ## This is y-data:
>>>> y_data = (([0.2867, 0.1171, -0.0087, 0.1326, 0.2415, 0.2878, 0.3133, 0.3701, 0.3996, 0.3728, 0.3551, 0.3587, 0.1408, 0.0416, 0.0708, 0.1142, 0, 0, 0]))
>>>>
>>>> ## This is x-data:
>>>> t = (([67, 88, 104, 127, 138, 160, 169, 188, 196, 215, 240, 247, 271, 278, 303, 305, 321, 337, 353]))
>>>>
>>>> ## This is the equation:
>>>> fitfunc = lambda p, x:    p[0] + (p[1] -p[0]) * ((1/(1+exp(-p[2]*(t-p[3])))) + (1/(1+exp(p[4]*(t-p[5])))) -1)
>>>>
>>>> ##
>>>> errfunc = lambda p, x, y: fitfunc(p,x) -y
>>>>
>>>> guess = [0, max(y_data), 0.1, 140, -0.1, 270]
>>>>
>>>> bounds = [(-0.2, 0.1),(0.1,0.97), (0.05,0.8), (120,190), (-0.8, -0.05), (200,300) ]
>>>>
>>>> ## This seems to work ok:
>>>> p2,success = optimize.leastsq(errfunc, guess, args=(t, y_data),full_output=0)
>>>> print 'Estimates from leastsq \n', p2,success
>>>>
>>>>
>>>> ## But this does not:
>>>> best, val, d = optimize.fmin_l_bfgs_b(errfunc, guess, bounds=bounds, args=(t, y_data), iprint=2)
>>>
>>> The minimization routines, I believe, in fmin expect a function that
>>> maps from to a scalar.  So you need to tell fmin_l_bfgs that you want
>>> to minimize the sum of squared errors, optimze.leastsq assumes this.
>>> So just define one more function that sums the squared errors and
>>> minimize it
>>>
>>> errfuncsumsq = lambda p, x, y: np.sum(errfunc(p,x,y)**2)
>>>
>>> Now, run it without bounds to make sure we get the same thing
>>>
>>> boundsnone = [(None,None)]*6
>>>
>>> Notice that you also have to tell fmin_l_bfgs_b to approximate the
>>> gradient or else it assumes that your objective function also returns
>>> its gradient
>>>
>>> best, val, d = optimize.fmin_l_bfgs_b(errfuncsum, guess,
>>> approx_grad=True, bounds=boundsnone, args=(t, y_data), iprint=2)
>>>
>>> p2
>>> array([  6.79548883e-02,   3.68922503e-01,   7.55565728e-02,
>>>         1.41378227e+02,   2.91307814e+00,   2.70608242e+02])
>>>
>>> best
>>> array([  6.79585333e-02,  -2.33026316e-01,  -7.55409880e-02,
>>>         1.41388265e+02,  -1.36069434e+00,   2.70160779e+02])
>>>
>>
>> I just realized that these don't come up with the same thing.  I don't
>> have an answer for why yet.
>>
>> Skipper
>>
>
> Oh,
>
> ret = optimize.leastsq(errfunc, guess, args=(t,y_data))
>
> ret2 = optimize.fmin_l_bfgs_b(errfuncsumsq, guess, approx_grad=True,
> bounds=boundsnone, args=(t, y_data), iprint=2)
>
> fitfunc(ret[0],t)
> array([ 0.0690421 ,  0.0731951 ,  0.08481868,  0.14388978,  0.199337  ,
>        0.30971974,  0.33570587,  0.3602918 ,  0.36414477,  0.36777158,
>        0.36874788,  0.36881958,  0.14080121,  0.06794499,  0.06795339,
>        0.0679536 ,  0.0679545 ,  0.06795477,  0.06795485])
>
> fitfunc(ret2[0],t)
> array([ 0.06904625,  0.07319943,  0.0848205 ,  0.14386744,  0.19929593,
>        0.30968735,  0.3356897 ,  0.36029973,  0.3641578 ,  0.36779021,
>        0.36876834,  0.3688402 ,  0.14077023,  0.06795562,  0.06795703,
>        0.06795724,  0.06795815,  0.06795842,  0.0679585 ])
>
> errfuncsumsq(ret[0], t, y_data)
> 0.079297668259408899
>
> errfuncsumsq(ret2[0], t, y_data)
> 0.079298042836826454
>

Very nice example !
However, is it correct, that *within* the given bounds the result is
just a constant ?
>>> ret3 = optimize.fmin_l_bfgs_b(errfuncsumsq, guess, approx_grad=True, bounds=bounds, args=(t, y_data), iprint=2)
>>> ret3
([  1.00000000e-01   1.00000000e-01   5.00000000e-02   1.39979309e+02
  -5.00000000e-02   2.70003604e+02], 0.55408092, {'warnflag': 0,
'task': 'CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL', 'grad':
array([-4.41837799,  1.03037829,  0.        ,  0.        ,  0.
,  0.        ]), 'funcalls': 2})
>>>
>>> errfuncsumsq(ret3[0], t, y_data)
0.55408092
>>>
>>>
>>> bounds
[(-0.2, 0.1), (0.1, 0.97), (0.05, 0.8), (120, 190), (-0.8, -0.05), (200, 300)]
>>> fitfunc(ret3[0],t)
[ 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
  0.1  0.1  0.1  0.1]
>>>

(Note:  fitfunc above used wronge 't' instead of 'x', correct is:
>>> fitfunc = lambda p, x:    p[0] + (p[1] -p[0]) * ((1/(1+exp(-p[2]*(x-p[3])))) + (1/(1+exp(p[4]*(x-p[5])))) -1)
)

 Thanks,
Sebastian Haase
_______________________________________________
SciPy-User mailing list
SciPy-User at scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

Jepp thanks. - For me it was a very, very helpful lesson.

I realize that equation is here somewhat taken out of the air, sorry. The equation is slightly modified from Beck et al 2006. "Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. Remote Sensing of Environment 100:321-334". So the first parameter (p[0]) estimates the minimum NDVI/EVI. Parameter 2 estimates maximum NDVI. Parameter 3 & 5  estimates the rate of increase in EVI in spring and the rate of decrease in NDVI in fall. Parameter 4 & 6 estimates the inflection points in spring and fall (cf. Fig 3 in Beck et al). 

Thanks to you both,

Torkild