[SciPy-User] optimize.leastsq does not converge with full Jacobian?

[Sebastian Walter]

is it ok with you if I use your code in the unit test of the automatic
differentiation tool pyadolc?
(http://github.com/b45ch1/pyadolc)

regards,
Sebastian



On Thu, Oct 1, 2009 at 9:12 AM, Sturla Molden <sturla at molden.no> wrote:
> I was trying to fit some Michaelis-Menten data from Bates and Watts
> (1988) puromycin experiment to test non-linear regression in SciPy. I
> get some very strange results that I don't understand.
>
> First the non-linear model to be fitted is:
>
> y = Vmax * x / (x + Km)
>
> Not using Jacobian works:
>
> import numpy as np
> import scipy
> from scipy.linalg import qr, solve
> from scipy.optimize import leastsq
>
>
> # Bates and Watts (1988) puromycin data
>
> data = ((0.02, 47, 76),
>        (0.06, 97, 107),
>        (0.11, 123, 139),
>        (0.22, 152, 159),
>        (0.56, 191, 201),
>        (1.10, 200, 207))
>
> data = np.array(data)
>
> X = data[:,0:1].repeat(2,axis=1).flatten()
> Y = data[:,1:].flatten()
>
> # initial fit from Linewaver-Burk plot
> y = Y**-1
> x = np.vstack((np.ones(X.shape),X**-1)).T
> q,r = qr(x, econ=True)
> b = solve(r, (np.mat(y) * np.mat(q)).T).ravel()
>
> # Michaelis-Menten fit from Lineweaver-Burk
> Vmax = 1.0 / b[0]
> Km = b[1] * Vmax
>
> # refit with Levenberg-Marquardt method
>
> def michaelis_menten(t, x):
>    Vmax, Km = t
>    return Vmax*x/(x + Km)
>
> def residuals(t, x, y):
>    return y - michaelis_menten(t, x)
>
> (Vmax,Km),ierr = leastsq(residuals, (Vmax,Km), args=(X,Y))
>
> This gives Vmax,Km = 2.12683559e+02, 6.41209954e-02, which is the
> "correct" answer.
>
> However, when I use the full Jacobian,
>
> def jacobian(t, x, y):
>    j = np.zeros((2,x.shape[0]))
>    Vmax, Km = t
>    j[0,:] = x/(Km + x)
>    j[1,:] = -Vmax*x/((Km + x)**2)
>    return j
>
> (Vmax,Km),ierr = leastsq(residuals, (Vmax,Km), args=(X,Y),
> Dfun=jacobian, col_deriv=1)
>
> I always get the start value of (Vmax,Km) returned, which is (195.8027,
> 0.0484065).
>
> What on earth is going on?
>
> Is there a bug in SciPy or am I being incredibly stupid?
>
>
> Sturla Molden
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>