LM(Levenberg–Marquardt)算法原理

LM算法作为非线性优化的“标准”方法，算法的数学原理有很多优秀的参考资料。我看过这些参考资料之后，觉得再重新写一遍已经是无力且多余的事情了。我简单说明一下这些参考资料，然后贴上自己的手写笔记。
参考资料：
1.《Methods for non-linear least squares problems》这本书将非线性最小二乘问题的优化方法讲了一大通，非常值得一看。因为LM算法是从Gauss-Newton方法演进来的，而Gauss-Newton方法又是从Newton方法演进过来的，所以追根溯源应该从Newton法开始看起。而比Newton法更简洁的就是最速下降法了，这本书将所有的非线性优化问题讲了个底朝天，聪明人仔细读一读不吃亏。
2.A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar
如果说上面那本书是准备好给搞理论看的版本的话，那这篇文章一定就是准备好给工程师看的了，文章对LM算法的实现给出了很好的讲解，工程师读一下，醍醐灌顶就可以写代码了。
3.[blog]原理及C++实现：Levenberg–Marquardt算法学习
4.[blog]原理及matlab实现：Levenberg-Marquardt
5.[blog]另一篇python实现：Python 算例实现Levenberg-Marquardt算法
6.我的笔记，你可以放心略过的部分:)D

LM算法python实现

实现步骤：

在LM算法原理中提到的参考资料提供了一些算法实现的伪代码，但是他们略有不同，主要的不同点是在公式表述以及u、v的更新比率上有小的差异。
我运行过他们的部分代码，发现优化效果也能够快速收敛，并不影响实际效果。
我按照文章A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar中的步骤，重新写了python代码，代码实现步骤如下：

代码：

1.我随机产生了100个input_data,设定正确的参数a和b，然后按照我要拟合的公式a×np.exp(b×input_data)加上一些高斯噪声计算出了100个对应的output_data, 作为观察。
2.初始化参数a和b，使之不要与真实值太离谱
3.用LM算法对其优化拟合，画出拟合曲线和迭代误差曲线。

'''
#Implement LM algorithm only using basic python
#Author:Leo Ma
#For csmath2019 assignment4,ZheJiang University
#Date:2019.04.28
'''
import numpy as np
import matplotlib.pyplot as plt #input data, whose shape is (num_data,1)
#data_input=np.array([[0.25, 0.5, 1, 1.5, 2, 3, 4, 6, 8]]).T
#data_output=np.array([[19.21, 18.15, 15.36, 14.10, 12.89, 9.32, 7.45, 5.24, 3.01]]).Ttao = 10**-3
threshold_stop = 10**-15
threshold_step = 10**-15
threshold_residual = 10**-15
residual_memory = []#construct a user function
def my_Func(params,input_data):a = params[0,0]b = params[1,0]#c = params[2,0]#d = params[3,0]return a*np.exp(b*input_data)#return a*np.sin(b*input_data[:,0])+c*np.cos(d*input_data[:,1])#generating the input_data and output_data,whose shape both is (num_data,1)
def generate_data(params,num_data):x = np.array(np.linspace(0,10,num_data)).reshape(num_data,1)       # 产生包含噪声的数据mid,sigma = 0,5y = my_Func(params,x) + np.random.normal(mid, sigma, num_data).reshape(num_data,1)return x,y#calculating the derive of pointed parameter,whose shape is (num_data,1)
def cal_deriv(params,input_data,param_index):params1 = params.copy()params2 = params.copy()params1[param_index,0] += 0.000001params2[param_index,0] -= 0.000001data_est_output1 = my_Func(params1,input_data)data_est_output2 = my_Func(params2,input_data)return (data_est_output1 - data_est_output2) / 0.000002#calculating jacobian matrix,whose shape is (num_data,num_params)
def cal_Jacobian(params,input_data):num_params = np.shape(params)[0]num_data = np.shape(input_data)[0]J = np.zeros((num_data,num_params))for i in range(0,num_params):J[:,i] = list(cal_deriv(params,input_data,i))return J#calculating residual, whose shape is (num_data,1)
def cal_residual(params,input_data,output_data):data_est_output = my_Func(params,input_data)residual = output_data - data_est_outputreturn residual'''
#calculating Hessian matrix, whose shape is (num_params,num_params)
def cal_Hessian_LM(Jacobian,u,num_params):H = Jacobian.T.dot(Jacobian) + u*np.eye(num_params)return H#calculating g, whose shape is (num_params,1)
def cal_g(Jacobian,residual):g = Jacobian.T.dot(residual)return g#calculating s,whose shape is (num_params,1)
def cal_step(Hessian_LM,g):s = Hessian_LM.I.dot(g)return s'''#get the init u, using equation u=tao*max(Aii)
def get_init_u(A,tao):m = np.shape(A)[0]Aii = []for i in range(0,m):Aii.append(A[i,i])u = tao*max(Aii)return u#LM algorithm
def LM(num_iter,params,input_data,output_data):num_params = np.shape(params)[0]#the number of paramsk = 0#set the init iter count is 0#calculating the init residualresidual = cal_residual(params,input_data,output_data)#calculating the init Jocobian matrixJacobian = cal_Jacobian(params,input_data)A = Jacobian.T.dot(Jacobian)#calculating the init Ag = Jacobian.T.dot(residual)#calculating the init gradient gstop = (np.linalg.norm(g, ord=np.inf) <= threshold_stop)#set the init stopu = get_init_u(A,tao)#set the init uv = 2#set the init v=2while((not stop) and (k<num_iter)):k+=1while(1):Hessian_LM = A + u*np.eye(num_params)#calculating Hessian matrix in LMstep = np.linalg.inv(Hessian_LM).dot(g)#calculating the update stepif(np.linalg.norm(step) <= threshold_step):stop = Trueelse:new_params = params + step#update params using stepnew_residual = cal_residual(new_params,input_data,output_data)#get new residual using new paramsrou = (np.linalg.norm(residual)**2 - np.linalg.norm(new_residual)**2) / (step.T.dot(u*step+g))if rou > 0:params = new_paramsresidual = new_residualresidual_memory.append(np.linalg.norm(residual)**2)#print (np.linalg.norm(new_residual)**2)Jacobian = cal_Jacobian(params,input_data)#recalculating Jacobian matrix with new paramsA = Jacobian.T.dot(Jacobian)#recalculating Ag = Jacobian.T.dot(residual)#recalculating gradient gstop = (np.linalg.norm(g, ord=np.inf) <= threshold_stop) or (np.linalg.norm(residual)**2 <= threshold_residual)u = u*max(1/3,1-(2*rou-1)**3)v = 2else:u = u*vv = 2*vif(rou > 0 or stop):break;return paramsdef main():#set the true params for generate_data() functionparams = np.zeros((2,1))params[0,0]=10.0params[1,0]=0.8num_data = 100# set the data numberdata_input,data_output = generate_data(params,num_data)#generate data as requested#set the init params for LM algorithm params[0,0]=6.0params[1,0]=0.3#using LM algorithm estimate paramsnum_iter=100    # the number of iterationest_params = LM(num_iter,params,data_input,data_output)print(est_params)a_est=est_params[0,0]b_est=est_params[1,0]#老子画个图看看状况plt.scatter(data_input, data_output, color='b')x = np.arange(0, 100) * 0.1 #生成0-10的共100个数据，然后设置间距为0.1plt.plot(x,a_est*np.exp(b_est*x),'r',lw=1.0)plt.xlabel("2018.06.13")plt.savefig("result_LM.png")plt.show()plt.plot(residual_memory)plt.xlabel("2018.06.13")plt.savefig("error-iter.png")plt.show()if __name__ == '__main__':main()

运行结果：