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]) 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 ## 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 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 Sorry it wasn't true what I said in my first report. The leastsq() result when col_deriv is not 1 is not an exactly the same as the initial guess. But one does have to agree there's a bug somewhere. --lin On 6/23/07, Stefan van der Walt <stefan@sun.ac.za> wrote:
Hi Lin
On Fri, Jun 22, 2007 at 07:35:05PM -0700, Lin Shao wrote:
It seems like when calling leastsq(), if the Jacobian matrix is organized as rows of partial derivatives (w.r.t. each variable), then no optimization is done at all -- the return value is the same as the initial guess. It only works when the matrix is columns of derivatives and set the parameter col_deriv to 1.
It would be helpful if you could provide two short snippets of code to illustrate the problems you mention.
Regards Stéfan _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user