[SciPy-user] Another leastsq Jacobian bug

[Lin Shao]

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 at 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 at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>