通过雅可比矩阵进行 lmfit 最小二乘法

Question

我正在使用 lmfit python 包进行非线性优化（网址：http://lmfit.github.io/lmfit-py/）。我想知道在使用最小二乘拟合方法时是否可以通过雅可比函数？如果是的话，是否可以为我提供一个最小的例子？

谢谢！ 卡布罗

P.S.：代码如下：

# 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：我在 github 上找到了这个有价值的例子：https://github.com/lmfit/lmfit-py/blob/master/examples/example_derivfunc.py

注意事项：在设法将 Jacobian 函数传递给 lmft.leastsq 之后，我意识到，在我的测试用例中，由 lmfit 返回的优化解决方案不再收敛到真正的解决方案。但是，当使用实际的 scipy.optimize.leastq 函数（由 lmfit 调用）时，一切正常，即返回的解决方案收敛，也包括适合的雅可比行列式。我并不是说 lmfit.leastsq 在为它提供雅可比函数时不能正常工作，但我建议谨慎对待这种情况。到目前为止，我还没有时间深入研究造成这种情况的原因。

Answer 1

您可以传入一个函数来计算雅可比矩阵作为Dfun 参数到Minimizer.leastsq 方法：http://lmfit.github.io/lmfit-py/fitting.html#leastsq

默认情况下，为Dfun 传入的函数应返回一个数组，该数组具有与参数相同的行数，并且每一行对应于每个被拟合的参数的导数。确保使用Parameters 对象指定参数，以便以正确的顺序处理参数。我相信这是必要的 IIRC，尽管它可能无关紧要。