非线性最小二乘法拟合函数-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 离散数据的因变量约为自变量的两倍，因此上述计算结果符合预期