title: pass Jacobian for lmfit least squares
body: <p>I am using the lmfit python package for non-linear optimisation (url: <a href="http://lmfit.github.io/lmfit-py/" rel="nofollow">http://lmfit.github.io/lmfit-py/</a>). 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 ?</p>
  
  <p>Thank you!
  kabrrrp</p>
  
  <p>P.S.: code is as follow:</p>
  
  <pre><code># 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) &gt;&gt;&gt;&gt; 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 &lt; 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)
  </code></pre>
  
  <p>P.P.S: I have found this valuable example at github: <a href="https://github.com/lmfit/lmfit-py/blob/master/examples/example_derivfunc.py" rel="nofollow">https://github.com/lmfit/lmfit-py/blob/master/examples/example_derivfunc.py</a></p>
  
  <p>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.</p>