Hi there, 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. --lin
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
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
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
## 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
Lin Shao wrote:
## 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
you forgot about the chain rule . Christian
No I didn't. The simplest refute to your answer is that your J[1] is -2*params[0]*(xx*params[1]-params[1]^2). How could there be a term with params[1]'s quadratic in the derivative if the original function is a quadratic function? On 6/25/07, Christian K <ckkart@hoc.net> wrote:
Lin Shao wrote:
## 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
you forgot about the chain rule .
Christian
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
Lin Shao wrote:
No I didn't. The simplest refute to your answer is that your J[1] is -2*params[0]*(xx*params[1]-params[1]^2). How could there be a term with params[1]'s quadratic in the derivative if the original function is a quadratic function?
On 6/25/07, Christian K <ckkart@hoc.net> wrote:
Lin Shao wrote:
## 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 you forgot about the chain rule .
You're right. Sorry. Christian
participants (3)
-
Christian K -
Lin Shao -
Stefan van der Walt