Levenberg-Marquardt (LM) 算法进行非线性拟合

1. LM算法

LM算法是一种非线性最小二乘优化算法，用于求解非线性最小化问题。LM主要用于解决具有误差函数的非线性最小二乘问题，其中误差函数是参数的非线性函数，需要通过调整参数使误差函数最小化。算法的基本思想是通过迭代的方式逐步调整参数，使得误差函数在参数空间中逐渐收敛到最小值。在每一次迭代中，算法通过求解一个线性方程组来更新参数。这个线性方程组由误差函数的雅可比矩阵和参数更新量构成。

LM算法的优点在于它能够快速收敛到局部最小值，并且对于初始参数的选择不太敏感。此外，算法还能够处理参数个数多于观测数据个数的问题，并且对于存在噪声的数据也比较鲁棒。

2. 调包实现

如图1所示，调用scipy.optimize的least_squares函数实现对测试函数 $exp(-ax^2-by^2)$ 的拟合结果。目标参数为 $[0.5, 0.5]$ ，初始参数设置为 $[1.0, 1.0]$ ，经过22次迭代，由于观测值暂未添加噪声，所以最终拟合参数与目标参数完全一致。

在这里插入图片描述

Fig. 1. 三维目标拟合: $\exp(-ax^2-by^2)$

3. LM算法实现

使用Python对LM做了简单实现，并对测试函数 $exp(ax^2+bx+c)$ 进行拟合，观测值添加高斯噪声。目标参数为 $[1.0, 2.0, 3.0]$ ，初始参数设置为 $[3.0, 9.0, 6.0]$ ，经过41次迭代，拟合参数为 $[2.0, 0.6, 3.5]$ ，MSE损失小于0.000001，符合拟合误差要求。图2绘制了第12（蓝），13（黄），15（绿）次迭代结果以及最终拟合结果（红）。

在这里插入图片描述

Fig. 2. 二维目标拟合: $\exp(ax^2+bx+c)$

# 部分函数代码：

def Func(abc,iput):   # 需要拟合的函数，abc是包含三个参数的一个矩阵[[a],[b],[c]]
    a = abc[0,0]
    b = abc[1,0]
    c = abc[2,0]
    return np.exp(a*iput**2+b*iput+c)
 
def Deriv(abc,iput,n):  # 对函数求偏导
    x1 = abc.copy()
    x2 = abc.copy()
    x1[n,0] -= 0.000001
    x2[n,0] += 0.000001
    p1 = Func(x1,iput)
    p2 = Func(x2,iput)
    d = (p2-p1)*1.0/(0.000002)
    return d
         
xk_l = []  # 用来存放每次迭代的结果
while conve:
       
    mse,mse_tmp = 0,0
    step += 1  
    fx = Func(xk,h) - y
    mse += sum(fx**2)
    for j in range(3): 
        J[:,j] = Deriv(xk,h,j) # 数值求导                                                    
    mse /= n  # 范围约束
 
    H = J.T*J + u*np.eye(3)   # 3*3
    dx = -H.I * J.T*fx        # 
    xk_tmp = xk.copy()
    xk_tmp += dx
    fx_tmp =  Func(xk_tmp,h) - y  
    mse_tmp = sum(fx_tmp[:,0]**2)
    mse_tmp /= n
    #判断是否下降
    q = float((mse - mse_tmp)/((0.5*dx.T*(u*dx - J.T*fx))[0,0]))
    if q > 0:
        s = 1.0/3.0
        v = 2
        mse = mse_tmp
        xk = xk_tmp
        temp = 1 - pow(2*q-1,3)
 
        if s > temp:
            u = u*s
        else:
            u = u*temp
    else:
        u = u*v
        v = 2*v
        xk = xk_tmp
 
    print ("step = %d,abs(mse-lase_mse) = %.8f" %(step,abs(mse-lase_mse)))  
    if abs(mse-lase_mse)<0.000001:
        break
       
    lase_mse = mse  # 记录上一个 mse 的位置
    conve -= 1
    
    xk_l.append(xk)