[SciPy-user] Another leastsq Jacobian bug

[Lin Shao]

> > ## 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]  ?

I think I was right. Think about what's the derivative of -x^2 -- it's -2x


>
> 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.

Ok, thanks for the tip.

> > ## 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
>

As said before, it doesn't seg fault everywhere. It happens on my
Linux 2.6.18, i686, same scipy and numpy version as yours.


> > ## 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.
>
I think this last bug is most important, because based on my
understanding of the docstring, the difference between col_deriv=0 and
1 is just whether the Jacobian is transposed or not. But apparently
there's more.


--lin