从数学推导至程序实现，对非线性最小二乘拟合原理及应用进行简单的介绍，希望能够帮助到有需要的人。
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 结束迭代的条件
- 残差的范数小于某个值
- 超出最大循环次数
评论 (0)