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

非线性最小二乘法拟合函数-2 问题1：假设有$n$个离散点$(x_i,y_i)$，存在一个系数未知的被拟合函数$f(x,a)$，如何使用非线性最小二乘法拟合得到最优的待定系数$a$。 问题2：假设有$n$个离散点$(x_i,y_i)$，存在一个系数未知的被拟合函数$f(x,a_1,a_2,\dots,a_m)$，如何使用非线性最小二乘法拟合得到最优的待定系数组$a_1,a_2,\dots,a_m$。 1 数学推导 上一个博客讲解了问题1的数学推导，这里关注待定系数有多个的问题2。 1.1 问题数学化 对于问题2，第$i$个数据的残差表达式如下 $$ r_i (a_1,a_2,\dots,a_m) = y_i - f(x_i,a_1,a_2,\dots,a_m) = y_i - f_i(a_1,a_2,\dots,a_m) $$ 上式中，$x_i,y_i$均为已知量，为书写方便，将$f(x_i,a_1,a_2,\dots,a_m)$简写为$f_i(a_1,a_2,\dots,a_m)$。 使用残差的平方和作为损失函数，如下 $$ L(a_1,a_2,\dots,a_m) = \sum_{i=1}^n r_i^2(a_1,a_2,\dots,a_m) = \sum_{i=1}^n[y_i - f_i(a_1,a_2,\dots,a_m)]^2 $$ 则问题2可表示为求解损失函数$L(a_1,a_2,\dots,a_m)$最小值对应的待定系数组$a_1,a_2,\dots,a_m$。同样的，将最小值问题近似为极小值问题。为书写方便，使用向量$\mathbf{a}$表示$a_1,a_2,\dots,a_m$，则上式修改为 $$ L(\mathbf{a}) = \sum_{i=1}^n [y_i - f_i(\mathbf{a})]^2 $$ 1.2 求函数极值点 1.2.1 迭代法递推公式 求解函数的极值点问题，即为求解一阶导数的零点问题，求函数$L(\mathbf{a})$对$\mathbf{a}$的一阶导数，如下 $$ \frac{d L(\mathbf{a})}{d \mathbf{a}} = -2 \sum_{i=1}^n [y_i - f_i(\mathbf{a})] \frac{d f_i(\mathbf{a})}{d \mathbf{a}} = -2 \sum_{i=1}^n [y_i - f_i(\mathbf{a})] J_i(\mathbf{a}) $$ 其中$J_i(\mathbf{a}) = d f_i(\mathbf{a}) / d \mathbf{a}$。令 $$ g(\mathbf{a}) = -\frac{1}{2} \frac{d L(\mathbf{a})}{d \mathbf{a}} = \sum_{i=1}^n J_i(\mathbf{a}) [y_i - f_i(\mathbf{a})] $$ 至此，将求解函数$L(\mathbf{a})$的极值点问题转化为求解函数$g(\mathbf{a})$的零点问题。对$g(\mathbf{a})$在$\mathbf{a}=\mathbf{a_0}$处进行泰勒展开，其中$\mathbf{a_0}$表示待定系数组$\mathbf{a}$的初始值，只保留一阶项有 $$ g(\mathbf{a}) \approx \sum_{i=1}^n J_i(\mathbf{a_0}) [y_i - f_i(\mathbf{a_0}) - J_i(\mathbf{a_0}) (\mathbf{a} - \mathbf{a_0})] = \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) - \mathbf{H} (\mathbf{a_0}) \Delta \mathbf{a} $$ 其中$\mathbf{H} (\mathbf{a}) = \mathbf{J^T} \mathbf{J}$为Hession矩阵，$\mathbf{J}(\mathbf{a})$为Jacobian矩阵，表达式如下 $$ \mathbf{J} (\mathbf{a}) = \left[ \begin{matrix} \frac{\partial J_1}{\partial a_1} & \frac{\partial J_1}{\partial a_2} & \cdots & \frac{\partial J_1}{\partial a_m} \\ \frac{\partial J_2}{\partial a_1} & \frac{\partial J_2}{\partial a_2} & \cdots & \frac{\partial J_2}{\partial a_m} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial J_n}{\partial a_1} & \frac{\partial J_n}{\partial a_2} & \cdots & \frac{\partial J_n}{\partial a_m} \end{matrix} \right] $$ 令$g(\mathbf{a})=0$可得 $$ \mathbf{H} (\mathbf{a_0}) \Delta \mathbf{a} = \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) $$ 上式即为高斯-牛顿法迭代公式。引入阻尼因子$\mu$得到Levenberg-Marquardt算法迭代公式如下 $$ [ \mathbf{H} (\mathbf{a_0}) + \mu \mathbf{I} ] \Delta \mathbf{a} = \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) $$ 其中$\mathbf{I}$为单位矩阵。整理后可得 $$ \Delta \mathbf{a} = [ \mathbf{H} (\mathbf{a_0}) + \mu \mathbf{I} ]^{-1} \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) $$ 进而可以得到递推公式为 $$ \mathbf{a_{k+1}} = \mathbf{a_{k}} + [ \mathbf{H} (\mathbf{a_0}) + \mu \mathbf{I} ]^{-1} \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) $$ 至此，得到了非线性最小二乘拟合的Levenberg-Marquardt算法迭代公式。 上述迭代公式得到的是待定系数组。 1.2.2 Jocabian矩阵 对于可直接求导得到一阶导数解析表达式的函数，直接代入相应的解析表达式即可。对于只能数值求解的函数，可使用差分近似求解，即 $$ J(\mathbf{a}) \approx \frac{f(\mathbf{a}+\delta \mathbf{a}) - f(\mathbf{a})}{\delta \mathbf{a}} $$ 其中$\delta \mathbf{a}$表示一个小量。至此，即可使用Levenberg-Marquardt算法迭代公式进行迭代求解。 1.2.3 结束迭代的条件 残差的范数小于某个值 超出最大循环次数 2 两个问题的联系 在本文中，Jacobian矩阵的大小为n行m列，Hession矩阵的大小为m行m列。其中，n为离散点的个数，m为待定系数的个数。当$m=1$时，问题2退化为问题1，即Jacobian矩阵变为n行1列的向量，Hession矩阵变为一个标量。

非线性最小二乘法拟合函数-1 从数学推导至程序实现，对非线性最小二乘拟合原理及应用进行简单的介绍，希望能够帮助到有需要的人。 0 问题引入 问题1：假设有$n$个离散点$(x_i,y_i)$，存在一个系数未知的被拟合函数$f(x,a)$，如何使用非线性最小二乘法拟合得到最优的待定系数$a$。 问题2：假设有$n$个离散点$(x_i,y_i)$，存在一个系数未知的被拟合函数$f(x,a_1,a_2,\dots,a_m)$，如何使用非线性最小二乘法拟合得到最优的待定系数组$a_1,a_2,\dots,a_m$。 问题2是问题1的拓展，由一个待定系数变为多个待定系数。之后的介绍和推导将由浅入深，从问题1开始讲起。 1 非线性最小二乘法 概念：非线性最小二乘法是一种优化方法，它通过最小化非线性函数与实验数据之间的残差平方和来确定未知参数的最优解。 非线性函数：参数未知的非线性被拟合函数，常用的有多项式函数、指数函数、幂函数、高斯函数等。 残差：是指在拟合过程中，模型预测值与实际观测值之间的差异。 对残差的要求： 独立同分布：残差之间应该是独立的且来自同一个概率分布。一般指高斯分布，如果不是，可能会出现偏差、误差较大的问题； 零均值：残差的平均值应该为0； 常数方差：残差的方差应该为常数。 求解方法：通常使用迭代优化算法来寻找最优解，如高斯-牛顿法、Levenberg-Marquardt算法等。 2 迭代法 概念：迭代法也称辗转法，是一种不断用变量的旧值递推新值的过程。迭代法是用计算机解决问题的一种基本方法，它利用计算机运算速度快、适合做重复性操作的特点，让计算机对一组指令(或一定步骤)进行重复执行，在每次执行这组指令(或这些步骤)时，都从变量的原值推出它的一个新值，迭代法又分为精确迭代和近似迭代。 百度百科——迭代法 下面介绍两种使用迭代法的场景。 2.1 求函数的零点 对函数$f(x)$在$x=x_0$处进行泰勒展开，只保留到一阶项，即 $$ f(x) \approx f(x_0) + f'(x_0) (x - x_0) $$ 令上式等于零，若函数$f(x)$的一阶导数存在且不为零，可得 $$ x = x_0 - f(x_0)/f'(x_0) $$ 进而可以得到递推公式 $$ x_{k+1} = x_k - f(x_k)/f'(x_k) $$ 设增量$\Delta = x_{k+1} - x_k$，则 $$ \Delta = - f(x_k)/f'(x_k) $$ 根据递推公式即可完成迭代求解，本方法为牛顿迭代法。 2.2 求函数的极值点 已知极值点处的一阶导数为零，即问题可修改为求函数$f'(x)$的零点。利用2.1节的方法，对$f'(x)$进行泰勒展开，可得 $$ f'(x) \approx f'(x_0) + f''(x_0) (x - x_0) $$ 同样的，可得增量表达式为 $$ \Delta = - f'(x_k)/f''(x_k) $$ ## 3 数学推导 由浅入深，先讲解待定系数只有一个的问题1。 3.1 问题数学化 问题1：假设有$n$个离散点$(x_i,y_i)$，存在一个系数未知的被拟合函数$f(x,a)$，如何使用非线性最小二乘法拟合得到最优的待定系数$a$。 对于问题1，第$i$个数据的残差表达式如下 $$ r_i (a) = y_i - f(x_i,a) = y_i - f_i(a) $$ 上式中，$x_i,y_i$均为已知量，为书写方便，将$f(x_i,a)$简写为$f_i(a)$。 根据第1节的概念，使用残差的平方和作为损失函数，表达式如下 $$ L(a) = \sum_{i=1}^n r_i^2(a) = \sum_{i=1}^n[y_i - f_i(a)]^2 $$ 则问题1可表示为求解损失函数$L(a)$最小值对应的待定系数$a$。最小值的问题通常不容易求解，因此可将其近似为求解损失函数$L(a)$极小值对应的待定系数$a$。若最终求得的极小值不是最小值，则为局部最优解。 3.2 求函数L(a)的极值点 3.2.1 迭代法递推公式 求解函数的极值点问题，即为求解一阶导数的零点问题，求函数$L(a)$对$a$的一阶导数，如下 $$ \frac{d L(a)}{d a} = -2 \sum_{i=1}^n [y_i - f_i(a)] \frac{d f_i(a)}{d a} $$ 令$J_i (a) = d f_i (a) / d a$，则有 $$ g(a) = -\frac{1}{2} \frac{d L(a)}{d a} = \sum_{i=1}^n [y_i - f_i(a)] J_i (a) $$ 至此，将求解函数$L(a)$的极值点问题转化为求解函数$g(a)$的零点问题。对$g(a)$在$a=a_0$处进行泰勒展开，并且只保留一阶项，有 $$ g(a) \approx \sum_{i=1}^n J_i (a_0) [y_i - f_i(a_0) - J_i(a_0) (a - a_0)] = \sum_{i=1}^n J_i (a_0) r_i(a_0) - H(a_0) \Delta a $$ 其中$H(a_0) = \sum_{i=1}^n J_i^2(a_0)$，$\Delta a = a - a_0$。令$g(a)=0$可得 $$ H(a_0) \Delta a = \sum_{i=1}^n J_i (a_0) r_i(a_0) $$ 上式即为高斯-牛顿法迭代公式。引入阻尼因子$\mu$得到Levenberg-Marquardt算法迭代公式如下 $$ [H(a_0) + \mu] \Delta a = \sum_{i=1}^n J_i (a_0) r_i(a_0) $$ 整理后可得 $$ \Delta a = [H(a_0) + \mu]^{-1} \sum_{i=1}^n J_i(a_0) r_i(a_0) $$ 进而可以得到递推公式为 $$ a_{k+1} = a_k + [H(a_k) + \mu]^{-1} \sum_{i=1}^n J_i(a_k) r_i(a_k) $$ 至此，得到了非线性最小二乘拟合的Levenberg-Marquardt算法迭代公式。 阻尼因子：是用来平衡算法的收敛速度和稳定性的一个参数。 当阻尼因子较小时，迭代算法会更快地收敛，但也会更容易陷入局部最小值； 当阻尼因子较大时，算法会更稳定，但会收敛得更慢。 3.2.2 一阶导数 对于可直接求导得到一阶导数解析表达式的函数，直接代入相应的解析表达式即可。对于只能数值求解的函数，可使用差分近似求解，即 $$ J(a) \approx \frac{f(a+\delta a) - f(a)}{\delta a} $$ 其中$\delta a$表示一个小量。至此，即可使用Levenberg-Marquardt算法迭代公式进行迭代求解。 3.2.3 结束迭代的条件 残差的范数小于某个值 超出最大循环次数