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