非线性最小二乘法拟合函数-3

在前面两个博客中，推导得到了Levenberg-Marquardt算法（简称LM算法）的迭代公式，这里将讲述如何使用Fortran编写一个简单的LM算法。

此处介绍待定系数只有一个的情况。

1 程序设计

1.1 Levenberg_Maquardt_Fit

• 简介：LM算法子例程，输入离散数据，迭代拟合待定系数

• 传递参数：

  • n 离散数据的长度
  • x 离散数据的自变量
  • y 离散数据的因变量
  • a 待定系数

1.2 myfunc

• 简介：待拟合函数的数值计算子例程。为了具有普适性，假定待拟合函数的解析式未知，函数结果只能由此子例程数值计算得到。
• 传递参数：
  • n 离散数据的长度
  • x 离散数据的自变量
  • a 待定系数的当前值
  • fx 待拟合函数的数值计算结果

1.3 Calculate_Jacobian

• 简介：计算待拟合函数的Jacobian矩阵的子例程。为了具有普适性，使用差分法计算一阶导数。
• 传递参数：
  • n 离散数据的长度
  • x 离散数据的自变量
  • y 待拟合函数的数值计算结果
  • a 待定系数的当前值
  • J 待拟合函数的Jacobian矩阵

2 源代码

2.1 Levenberg_Maquardt_Fit

• 此处迭代次数上限设置为30次，可根据需求自行更改
• 对于自变量有多个的情况，可将输入x改为x1,x2,...，满足自变量个数即可

SUBROUTINE Levenberg_Maquardt_Fit(n, x, y, a)
    IMPLICIT NONE

    ! 输入参数
    INTEGER, INTENT(IN) :: n ! 数据点个数
    REAL, DIMENSION(n), INTENT(IN) :: x, y ! 数据点
    ! REAL, DIMENSION(3), INTENT(INOUT) :: a ! 待求的拟合系数
    REAL, INTENT(INOUT) :: a ! 待求的拟合系数

    ! 定义常量
    INTEGER, PARAMETER :: m = 1 ! 待求的系数个数
    REAL, PARAMETER :: eps = 1.0E-6 ! 收敛阈值

    ! 定义变量
    REAL :: da(m), r(n), J(n,m), H(m,m) ! 拟合系数、残差、雅可比矩阵、Hessian矩阵
    REAL :: lambda, alpha(m) ! 调节因子、步长
    INTEGER :: i, iter ! 循环计数器

    REAL :: fx(n)

    ! 初始化调节因子
    lambda = 0.001

    ! 开始迭代
    iter = 0
    DO WHILE(iter < 30) ! 迭代次数上限为10000
      iter = iter + 1

      ! 计算被拟合函数值
      CALL myfunc(n,x,a,fx)

      ! 计算残差向量
      r = y - fx

      ! 计算雅可比矩阵
      CALL Calculate_Jacobian(n, x, fx, a, J)

      ! 计算Hessian矩阵
      H = MATMUL(TRANSPOSE(J), J)

      ! 计算梯度向量
      da = MATMUL(TRANSPOSE(J), r)

      ! 计算搜索方向
      DO i = 1, m
        H(i,i) = H(i,i) + lambda
      END DO
    !   CALL SOLVE_LINEAR_SYSTEM(m, H, da, alpha)

      alpha = da(1)/H(1,1)

      ! 计算新的拟合系数
      a = a + alpha(1)

      ! 更新调节因子
      IF (NORM2(r) < eps) THEN
        EXIT
      ELSE IF (NORM2(r) < NORM2(r - MATMUL(J, alpha))/2) THEN
        lambda = lambda/10.
      ELSE
        lambda = lambda*10.
      END IF

    END DO

END SUBROUTINE

2.2 myfunc

• 此处使用了一个非常简单的一次函数作为示例，使用时需要将fx = a*x修改为实际的待拟合函数的形式。若待拟合函数是一个子例程，则在此处使用CALL调用
• 若在Levenberg_Maquardt_Fit子例程中修改了自变量个数，此处需要对应修改

SUBROUTINE myfunc(n,x,a,fx)
    IMPLICIT NONE
    ! 输入参数
    INTEGER, INTENT(IN) :: n
    REAL, INTENT(IN) :: x(n) ! 自变量
    REAL, INTENT(IN) :: a ! 拟合系数
    REAL, INTENT(OUT) :: fx(n)

    fx = a*x
END SUBROUTINE

2.3 Calculate_Jacobian

• 本子例程使用了差分法计算一阶导数，因此几乎适用于所有的待拟合函数
• 若在前面修改了自变量的个数，此处需对应修改

SUBROUTINE Calculate_Jacobian(n, x, y, a, J)
    IMPLICIT NONE

    ! 输入参数
    INTEGER, INTENT(IN) :: n ! 数据点个数
    REAL, DIMENSION(n), INTENT(IN) :: x, y ! 数据点
    REAL, INTENT(IN) :: a ! 拟合系数

    ! 输出参数
    REAL, DIMENSION(n,1), INTENT(OUT) :: J ! 雅可比矩阵

    ! 定义变量
    INTEGER :: i ! 循环计数器
    REAL :: dela, y1(n)

    ! 步进
    dela = 0.01

    ! 计算 a+dela 对应的 y1
    CALL myfunc(n,x,a+dela,y1)

    ! 计算雅可比
    DO i = 1, n
        J(i,1) = (y1(i) - y(i))/dela
    END DO
END SUBROUTINE

3 测试

• 编写程序主体，设定y约为x的两倍，调用LM子例程计算待定系数a

program name
    implicit none
    INTEGER n 
    real :: x(5),y(5),a 

    x = (/1.,2.,3.,4.,5./)
    y = (/2.1,4.1,5.9,8.1,9.9/)
    ! y = 2*x**2 + 

    n = 5
    a = 1.0

    call Levenberg_Maquardt_Fit(n, x, y, a)

    print *,a

end program name

• 使用gfortran编译生成可执行文件，运行可执行文件后，在终端打印出

   1.99818182

• 离散数据的因变量约为自变量的两倍，因此上述计算结果符合预期