问题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矩阵变为一个标量。
评论 (0)