pass Jacobian for lmfit least squares

Question

I am using the lmfit python package for non-linear optimisation (url: http://lmfit.github.io/lmfit-py/). I would like to know if it is possible to pass the Jacobian function when using the least squares fitting method ? If yes, would it be possible to provide me with a minimal example ?

Thank you! kabrrrp

P.S.: code is as follow:

# f(t,g,p) = dg_dt(t,g,p) = R*(c^h/(c^h+K^h))-l*g
# returns rhs of an ODE (dg_dt)
def hill_1g1c(t, g_in, p_in):
    R = p_in['R'].value
    K = p_in['K'].value
    h = p_in['h'].value
    l = p_in['l'].value

    dg_dt = R*((c_int(t)**h)/((c_int(t)**h)+(K**h))) - l*g_in
    return dg_dt

# f_deriv(t,g,p)
# is intended to return the derivatives of f with respect to 4 parameter
def hill_1g1c_jac(p_in, y_in, dt, t_max, g_data, g_err):
    t=1
    R = p_in['R'].value
    K = p_in['K'].value
    h = p_in['h'].value
    l = p_in['l'].value
    dg_dR = (c_int(t)**h) / (c_int(t)**h + K**h) - l * 1
    dg_dK = (-1 * R * c_int(t)**h * h * k**(h-1)) / ((c_int(t)**h + K**h)**2) - l * 1
    dg_dh = (-1 * R * c_int(t)**h * k**h * (np.log(k) - np.log(c_int(t)))) / ((c_int(t)**h + K**h)**2) - l * 1
    dg_dl = -y_in - l * 1

    return np.array([dg_dR, dg_dK, dg_dh, dg_dl])

# y = ODE_solve(y0,p,dt, t_max) >>>> wrapper around ode.integrate, returns array of g
def ODE_solve(y0, p_in, dt, t_max):
    t = [0]
    y = [y0]
    r.set_initial_value(y0, t=0.0)
    r.set_f_params(p_in)
    while r.successful() and r.t < t_max:
         r.integrate(r.t+dt)
         t.append(r.t)
         y.append(r.y)
    return np.array(y)

# weighted least squares, objective function to be minimised
def ODE_wres(p_in, y0, dt, t_max, g_data, g_err):
    g_extended = ODE_solve(y0, p_in, dt, t_max)
    g_model = g_extended[-25:]/g_extended[-25]
    weighted_residuals = (g_data - g_model)/(g_err + 0.00000001)
    return weighted_residuals


# specs for inegrate.ode
y0 = 1
t0 = 0
r = integrate.ode(hill_1g1c).set_integrator('vode', with_jacobian=False)
t_stim = 15
t_max = t_stim + 24
t_plus = 5
dt = 1
t_extended = np.linspace(0,t_max+t_plus,t_max+t_plus+1)

# set history of all inputs to 1
c_history = [1 for val in range(t_stim)]

# data (is read in from text file)
g_data = Y_data[:,i]
# error in g
g_err = Y_error[:,i]
# input c
c_data = Y_data[:,k]
# interpolation of c, contains history (=1) and future (=endval)
c_future = [c_data[-1] for val in range(t_plus)]
c_extended = np.hstack((c_history, c_data, c_future))
c_int = interp1d(t_extended, c_extended, kind='linear')

# initial parameter vector
R_ini = random.uniform(0.01, 500.0)
K_ini = random.uniform(0.01, 20.0)
h_ini = random.uniform(-4.0, 4.0)
l_ini = random.uniform(0.07, 7.0)

p_ini = Parameters()
p_ini.add('R', value= R_ini, min=0.01, max=500)
p_ini.add('K', value= K_ini, min=0.01, max=20)
p_ini.add('h', value= h_ini, min=-4, max=4)
p_ini.add('l', value= l_ini, min=0.07, max=7.0)

res_ini = ODE_wres(p_ini, y0, dt, t_max, g_data, g_err)
chisqr_ini = np.sum(res_ini**2)

# optimise
lmsol = Minimizer(ODE_wres, p_ini, fcn_args=(y0, dt, t_max, g_data, g_err))
lmsol.leastsq(Dfun=hill_1g1c_jac, col_deriv=True)

P.P.S: I have found this valuable example at github: https://github.com/lmfit/lmfit-py/blob/master/examples/example_derivfunc.py

NOTE OF CAUTION: After managing to pass the Jacobian function to lmft.leastsq, I realized, that in my test case the optimised solution as returned by lmfit was no longer converging to the true solution. However, when using the actual scipy.optimize.leastq function (which is called by lmfit) everything worked fine, that is the returned solution converged also including the Jacobian to fit. I am not saying that lmfit.leastsq doesn't work properly when supplying it with the Jacobian function, but I recommend to treat this case with some caution. So far I have not found the time to look deeper into the cause of this.

Answer 1

You can pass in a function to compute the Jacobian as the Dfun argument to the Minimizer.leastsq method: http://lmfit.github.io/lmfit-py/fitting.html#leastsq

By default the function passed in for Dfun should return an array with the same number of rows as parameters, and each row the derivative with respect to each parameter being fitted. Make sure to specify the parameters with a Parameters object so that the parameters are treated in the correct order. I believe this this necessary IIRC though it might not matter.