[SciPy-User] scipy.optimize.leastsq question
Joshua Holbrook
josh.holbrook at gmail.com
Tue Jul 6 14:19:38 EDT 2010
On Tue, Jul 6, 2010 at 10:02 AM, ms <devicerandom at gmail.com> wrote:
> 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 at 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/<http://www.tau.ac.il/%7Ekineret/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 at scipy.org
> >>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
> >>>>>>
> >>>>> _______________________________________________
> >>>>> SciPy-User mailing list
> >>>>> SciPy-User at scipy.org
> >>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
> >>>>
> >>>> --
> >>>>
> >>>> Cheers, Ralph
> >>>> _______________________________________________
> >>>> SciPy-User mailing list
> >>>> SciPy-User at scipy.org
> >>>> http://mail.scipy.org/mailman/listinfo/scipy-user
> >>>>
> >>> _______________________________________________
> >>> SciPy-User mailing list
> >>> SciPy-User at scipy.org
> >>> http://mail.scipy.org/mailman/listinfo/scipy-user
> >>>
> >> _______________________________________________
> >> SciPy-User mailing list
> >> SciPy-User at scipy.org
> >> http://mail.scipy.org/mailman/listinfo/scipy-user
> >>
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User at scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
m. said:
> 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.
I say this not as someone intimately familiar with scipy.optimize, but as
someone who has implemented a least squares-ish algorithm himself.
You are almost certainly correct in that leastsq calculates an approximate
Jacobian using a finite difference method on its own. However, if you can
symbolically differentiate your promblem without too much heartache, then
supplying an exact Jacobian is probably preferrable due to higher precision
and less function evaluation (f(x) and f(x+h), differenced and normalized,
vs. simply f'(x)).
On the other hand: When I implemented my algorithm (nearly two years ago),
my equations were pretty nasty. My derivatives just happened to be much much
worse (as can be seen at
http://modzer0.cs.uaf.edu/~jesusabdullah/gradients.html, at least for a
little while), and at the time sympy honestly wasn't production-ready. So, I
ended up using a finite difference method to calculate them (I believe I
used scipy's derivative function), with which I did have to tweak step
sizes.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100706/e25af8e2/attachment.html>
More information about the SciPy-User
mailing list