[SciPy-User] scipy.optimize.leastsq question
ms
devicerandom@gmail....
Tue Jul 6 13:02:22 CDT 2010
On 29/06/10 09:59, Sebastian Walter wrote:
>>> Only use derivative free optimization methods if your problem is not continuous.
>>> If your problem is differentiable, you should compute the Jacobian
>>> yourself, e.g. with
>>>
>>> def myJacobian(x):
>>> h = 10**-3
>>> # do finite differences approximation
>>> return ....
>>>
>>> and provide the Jacobian to
>>> scipy.optimize.leastsq(..., Dfun = myJacobian)
Uh, I am real newbie in this field, but I expected that the Jacobian was
needed if there was an analytical expression for the derivatives; I
thought the leastsq routine calculated the finite difference
approximation by itself otherwise. So I never bothered providing an
"approximate" Jacobian. Or maybe I do not get what do you mean by finite
difference.
Could someone provide some insight on this?
thanks!
m.
>>> This should work much better/reliable/faster than any of the alternatives.
>>
>> Maybe increasing the step length in the options to leastsq also works:
>>
>> epsfcn – A suitable step length for the forward-difference
>> approximation of the Jacobian (for Dfun=None).
>>
>> I don't think I have tried for leastsq, but for some fmin it works
>> much better with larger step length for the finite difference
>> approximation.
>
> choosing the right "step length" h is an art that I don't know much about.
> But apparently one rule of thumb is to use
>
> h = abs(x)* sqrt(numpy.finfo(float).eps)
> to compute
> f'(x) = (f(x+h) - f(x))/h
>
> i.e. if one has x = [1,10**-3, 10**4] one would have to scale h with
> 1, 10**-3 and 10**4.
>
> Regarding epsfcn: I find the documentation of leastsq a "little" confusing.
>
> epsfcn -- A suitable step length for the forward-difference
> approximation of the Jacobian (for Dfun=None). If
> epsfcn is less than the machine precision, it is assumed
> that the relative errors in the functions are of
> the order of the machine precision.
>
> In particular I don't quite get what is meant by "relative errors in
> the functions". Which "functions" does it refer to?
>
>
> Sebastian
>
>>
>> Josef
>>
>>
>>
>>>
>>> Also, using Algorithmic Differentiation to compute the Jacobian would
>>> probably help in terms of robustness and convergence speed of leastsq.
>>>
>>> Sebastian
>>>
>>>
>>>
>>>
>>>
>>>>
>>>> Cheers, Ralph
>>>>
>>>> Den 28.06.10 17.13, skrev Sebastian Walter:
>>>>> there may be others who have more experience with scipy.optimize.leastsq.
>>>>>
>>>>>> From the mathematical point of view you should be certain that your
>>>>> function is continuously differentiable or at least
>>>>> (Lipschitz-)continuous.
>>>>> This is because scipy.optimize.leastsq uses the Levenberg-Marquardt
>>>>> algorithm which requires the Jacobian J(x) = dF/dx.
>>>>>
>>>>> You do not provide an analytic Jacobian for scipy.optimize.leastsq.
>>>>> That means that scipy.optimize.leastsq uses some finite differences
>>>>> approximation to approximate the Jacobian J(x).
>>>>> It can happen that this finite differences approximation is so poor
>>>>> that no descent direction for the residual can be found.
>>>>>
>>>>> So the first thing I would check is if the Jacobian J(x) makes sense.
>>>>> You should be able to extract it from
>>>>> scipy.optimize.leastsq's output infodict['fjac'].
>>>>>
>>>>> Then I'd check if
>>>>> F(x + h*v) - F(x)/h, for h \approx 10**-8
>>>>>
>>>>> gives the same vector as dot(J(x),v)
>>>>> if this doesn't match at all, then your Jacobian is wrong resp. your
>>>>> function is not continuously differentiable.
>>>>>
>>>>> Hope this helps a little,
>>>>> Sebastian
>>>>>
>>>>>
>>>>>
>>>>> On Mon, Jun 28, 2010 at 2:36 PM, Ralph Kube<ralphkube@googlemail.com> wrote:
>>>>>> Hello people,
>>>>>> I am having a problem using the leastsq routine. My goal is to
>>>>>> determine three parameters r_i, r_s and ppw so that the residuals
>>>>>> to a model function a(r_i, r_s, ppw) to a measurement are minimal.
>>>>>> When I call the leastsq routine with a good guess of starting values, it
>>>>>> iterates 6 times without changing the vales of the initial parameters
>>>>>> and then exits without an error.
>>>>>> The function a is very complicated and expensive to evaluate. Some
>>>>>> evaluation is done by using the subprocess module of python. Can this
>>>>>> pose a problem for the leastsq routine?
>>>>>>
>>>>>>
>>>>>> This is in the main routine:
>>>>>>
>>>>>> import numpy as N
>>>>>>
>>>>>> for t_idx, t in enumerate(time_var):
>>>>>>
>>>>>> r_i = 300.
>>>>>> r_s = 1.0
>>>>>> ppw=1e-6
>>>>>> sza = 70.
>>>>>> wl = N.arange(300., 3001., 1.)
>>>>>>
>>>>>> albedo_true = compute_albedo(r_i, r_s, ppw, sza, wl)
>>>>>> # This emulates the measurement data
>>>>>> albedo_meas = albedo_true + 0.01*N.random.randn(len(wl))
>>>>>>
>>>>>> print 'Optimizing albedo'
>>>>>> p0 = [2.*r_i, 1.4*r_s, 4.*ppw]
>>>>>> plsq2 = leastsq(albedo_residual, p0, args=(albedo_meas, sza,
>>>>>> wl))
>>>>>> print '... done: ', plsq2[0][0], plsq2[0][1], plsq2[0][2]
>>>>>> albedo_model = compute_albedo(plsq2[0][0], plsq2[0][1], plsq2[0][2],
>>>>>> sza, wl)
>>>>>>
>>>>>> The residual function:
>>>>>> def albedo_residual(p, y, sza, wvl):
>>>>>> r_i, r_s, ppw = p
>>>>>> albedo = compute_albedo(r_i, r_s, ppw, sza, wvl)
>>>>>> err = albedo - y
>>>>>> print 'Albedo for r_i = %4.0f, r_s = %4.2f, ppw = %3.2e \
>>>>>> Residual squared: %5f' % (r_i, r_s, ppw, N.sum(err**2))
>>>>>>
>>>>>> return err
>>>>>>
>>>>>> The definition of the function a(r_i, r_s, ppw)
>>>>>> def compute_albedo(radius_ice, radius_soot, ppw, sza, wvl):
>>>>>>
>>>>>> The output is:
>>>>>> Optimizing albedo
>>>>>> Albedo for r_i = 600, r_s = 1.40, ppw = 4.00e-06 Residual squared:
>>>>>> 0.973819
>>>>>> Albedo for r_i = 600, r_s = 1.40, ppw = 4.00e-06 Residual squared:
>>>>>> 0.973819
>>>>>> Albedo for r_i = 600, r_s = 1.40, ppw = 4.00e-06 Residual squared:
>>>>>> 0.973819
>>>>>> Albedo for r_i = 600, r_s = 1.40, ppw = 4.00e-06 Residual squared:
>>>>>> 0.973819
>>>>>> Albedo for r_i = 600, r_s = 1.40, ppw = 4.00e-06 Residual squared:
>>>>>> 0.973819
>>>>>> Albedo for r_i = 600, r_s = 1.40, ppw = 4.00e-06 Residual squared:
>>>>>> 0.973819
>>>>>> ... done: 600.0 1.4 4e-06
>>>>>>
>>>>>> To check for errors, I implemented the example code from
>>>>>> http://www.tau.ac.il/~kineret/amit/scipy_tutorial/ in my code and it
>>>>>> runs successfully.
>>>>>>
>>>>>> I would be glad for any suggestion.
>>>>>>
>>>>>>
>>>>>> Cheers, Ralph
>>>>>> _______________________________________________
>>>>>> SciPy-User mailing list
>>>>>> SciPy-User@scipy.org
>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>
>>>>> _______________________________________________
>>>>> SciPy-User mailing list
>>>>> SciPy-User@scipy.org
>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>
>>>> --
>>>>
>>>> Cheers, Ralph
>>>> _______________________________________________
>>>> SciPy-User mailing list
>>>> SciPy-User@scipy.org
>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-User
mailing list