LM(Levenberg–Marquardt)算法原理及其python自定义实现

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]]).T


tao = 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,5
    y = 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.000001
    params2[param_index,0] -= 0.000001
    data_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_output
    return 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 params
    k = 0#set the init iter count is 0
    #calculating the init residual
    residual = cal_residual(params,input_data,output_data)
    #calculating the init Jocobian matrix
    Jacobian = cal_Jacobian(params,input_data)
    
    A = Jacobian.T.dot(Jacobian)#calculating the init A
    g = Jacobian.T.dot(residual)#calculating the init gradient g
    stop = (np.linalg.norm(g, ord=np.inf) <= threshold_stop)#set the init stop
    u = get_init_u(A,tao)#set the init u
    v = 2#set the init v=2
    
    while((not stop) and (k<num_iter)):
        k+=1
        while(1):
            Hessian_LM = A + u*np.eye(num_params)#calculating Hessian matrix in LM
            step = np.linalg.inv(Hessian_LM).dot(g)#calculating the update step
            if(np.linalg.norm(step) <= threshold_step):
                stop = True
            else:
                new_params = params + step#update params using step
                new_residual = cal_residual(new_params,input_data,output_data)#get new residual using new params
                rou = (np.linalg.norm(residual)**2 - np.linalg.norm(new_residual)**2) / (step.T.dot(u*step+g))
                if rou > 0:
                    params = new_params
                    residual = new_residual
                    residual_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 params
                    A = Jacobian.T.dot(Jacobian)#recalculating A
                    g = Jacobian.T.dot(residual)#recalculating gradient g
                    stop = (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 = 2
                else:
                    u = u*v
                    v = 2*v
            if(rou > 0 or stop):
                break;
        
    return params
  


        
        
def main():
    #set the true params for generate_data() function
    params = np.zeros((2,1))
    params[0,0]=10.0
    params[1,0]=0.8
    num_data = 100# set the data number
    data_input,data_output = generate_data(params,num_data)#generate data as requested
    
    #set the init params for LM algorithm 
    params[0,0]=6.0
    params[1,0]=0.3

    #using LM algorithm estimate params
    num_iter=100    # the number of iteration
    est_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.1
    plt.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()

运行结果:

在这里插入图片描述在这里插入图片描述

  • 81
    点赞
  • 434
    收藏
    觉得还不错? 一键收藏
  • 20
    评论
LM算法,全称为Levenberg-Marquard算法,它可用于解决非线性最小二乘问题,多用于曲线拟合等场合。 LM算法实现并不算难,它的关键是用模型函数 f 对待估参数向量 p 在其邻域内做线性近似,忽略掉二阶以上的导数项,从而转化为线性最小二乘问题,它具有收敛速度快等优点。LM算法属于一种“信赖域法”——所谓的信赖域法,此处稍微解释一下:在最优化算法中,都是要求一个函数的极小值,每一步迭代中,都要求目标函数值是下降的,而信赖域法,顾名思义,就是从初始点开始,先假设一个可以信赖的最大位移 s ,然后在以当前点为中心,以 s 为半径的区域内,通过寻找目标函数的一个近似函数(二次的)的最优点,来求解得到真正的位移。在得到了位移之后,再计算目标函数值,如果其使目标函数值的下降满足了一定条件,那么就说明这个位移是可靠的,则继续按此规则迭代计算下去;如果其不能使目标函数值的下降满足一定的条件,则应减小信赖域的范围,再重新求解。 事实上,你从所有可以找到的资料里看到的LM算法的说明,都可以找到类似于“如果目标函数值增大,则调整某系数再继续求解;如果目标函数值减小,则调整某系数再继续求解”的迭代过程,这种过程与上面所说的信赖域法是非常相似的,所以说LM算法是一种信赖域法。
Levenberg-Marquardt算法是非线性最小二乘问题的一种最优化方法,通常用于曲线拟合或参数估计问题。下面是该算法的一个简单实现示例: 假设我们有一组数据点(x,y),并且我们想要拟合一个函数f(x)来最小化均方误差。 首先,我们定义一个函数f(x,params)来表示所要拟合的函数,其中params代表待估参数。 下面是Levenberg-Marquardt算法的代码: ```python def levenberg_marquardt(x, y, param0, f, tol=1e-6, max_iter=1000): # x, y: 输入数据 # param0: 待估参数的初始值 # f: 所要拟合的函数 # tol: 收敛阈值 # max_iter: 最大迭代次数 # 返回:估计的参数值 # 初始化参数 params = np.array(param0) n = len(params) # 计算jacobian矩阵 def jacobian(x, params): h = 1e-6 J = np.zeros((len(x), n)) for i in range(n): p = params.copy() p[i] += h J[:, i] = (f(x, p) - f(x, params)) / h return J # 初始化lambda值和误差矩阵 lamda = 0.01 err = y - f(x, params) # 迭代 for i in range(max_iter): J = jacobian(x, params) # 计算增量 A = np.dot(J.T, J) + lamda * np.eye(n) b = np.dot(J.T, err) dp = np.linalg.solve(A, b) new_params = params + dp # 计算新的误差 new_err = y - f(x, new_params) if np.sum(new_err**2) < np.sum(err**2): params = new_params err = new_err lamda /= 10 else: lamda *= 10 # 判断是否达到收敛 if np.max(np.abs(dp)) < tol: return params # 达到最大迭代次数,返回估计的参数值 return params ``` 以上代码实现了Levenberg-Marquardt算法的主要步骤:计算jacobian矩阵、初始化lambda值和误差矩阵、计算增量并更新参数,以及判断是否达到收敛。 在实际使用中,可以根据具体的问题调整迭代次数、lambda值和收敛阈值等参数,以达到更好的拟合效果。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 20
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值