# [SciPy-User] scipy.optimize.leastsq question

Sebastian Walter sebastian.walter@gmail....
Tue Jun 29 03:59:21 CDT 2010

```On Tue, Jun 29, 2010 at 9:52 AM,  <josef.pktd@gmail.com> wrote:
> On Tue, Jun 29, 2010 at 3:46 AM, Sebastian Walter
> <sebastian.walter@gmail.com> wrote:
>> On Tue, Jun 29, 2010 at 9:28 AM, Ralph Kube <ralphkube@googlemail.com> wrote:
>>> Thank you, I found the error this way. The Jacobian is indeed very
>>> hard to compute, and the leastsq routine computes a zero Jacobian.
>>> The albedo function I want to minimize does not change for values
>>> of h \approx 10**-8, but on scales h \approx 10**-3.
>>> I now use the fmin function and working with other functions that do not
>>> require any information about the derivative. They seem more appropriate
>>> to my problem.
>>
>> 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)
>> 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)
>>>>>
>>>>> 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
>
```