[SciPy-user] Another leastsq Jacobian bug

Christian K ckkart@hoc....
Mon Jun 25 05:32:46 CDT 2007


Lin Shao wrote:
> Try this:
> 
> import numpy as N
> import scipy.optimize as O
> 
> ## Create my data to fit:
> x=N.arange(-10,10,dtype=N.float32)
> y=(x-1.234)**2+3.456
> ## I want to fit y=p0*(x-p1)^2+p2
> 
> ## Now I define my objective
> def obj_func(params, xx, yy, mode='col'):
>     return params[0] * (xx-params[1])**2 + params[2] - yy
> ## 'mode' is dummy here, because I want to use it for my Jocobian
> 
> ## Now define my Jacobian
> def Jacobian(params,xx,yy,mode='col'):
>     J = N.empty((len(params),xx.size))
>     J[0] = (xx-params[1])**2
>     J[1] = -2*params[0]*(xx-params[1])

shouldn't that be -2*params[0]*(xx-params[1])*params[1]  ?

>     J[2] = 1
>     if mode=='col':
>         return J
>     elif mode=='row:
>         return J.transpose()
>     else:
>         raise ValueError, "Unkown mode %s in Jacobian()" % mode
> ## Hopefully I did my calculus correctly
> 
> ## First, try without Jocobian given
> b=O.leastsq(obj_func, (1.,1.,2.), args=(x, y, 'col'))
> print b[0]
> ## the result is [ 1.  1.  2.]
> ## obviously a failure, but no warning is returned

Not really a failure - this is numerics. The step length for finding the numeric
derivative is too small (epsfcn keyword arg). Try with epsfcn = 1e-10.
Have look at the return message of leastsq, which is b[3] in case full_output=1.

> 
> ## Second, try Jacobian with col_deriv=1
> b=O.leastsq(obj_func, (1.,1.,2.), args=(x, y, 'col'), Dfun=Jacobian,
> full_output=1, col_deriv=1)
> ## full_output has to be 1 because of a seg fault bug I reported earlier

does not segfault here. unix, python2.5,
scipy 0.5.3.dev3062, numpy 1.0.2.dev3484

> print b[0]
> ## the result is [ 1.00000004  1.23400008  3.45599962]
> ## Good Job!
> 
> ## Last, try Jacobian with col_deriv=0
> b=O.leastsq(obj_func, (1.,1.,2.), args=(x, y, 'row'), Dfun=Jacobian,
> full_output=1, col_deriv=0)
> print b[0]
> ## the result is [ 1.02507777  1.222815    2.2651552 ]
> ## completely different result and wrong

I've no idea why this fails, though.

Christian



More information about the SciPy-user mailing list