完全解读：用最小二乘法求模型参数

要想确定N个未知参数，就得至少有N个线性独立的方程。其中：1、如果方程个数低于未知数个数，就称该问题为“欠定问题”；2、如果线性独立的方程个数大于未知数个数，则称该问题为“超定问题”。

最小二乘方法处理第二类问题(即超定问题)。

由于实验观测值yi会受到随机误差εi的影响，即有：

因此，最好能够建立尽可能多的方程(包含多对yi和xi)，从而获得更好、更可靠的统计结果，以更加准确地估计参数a。假设得到一组观测值yi及其相应的条件xi，数据拟合的目标就是要利用最小二乘方法，将观测值与基于模型函数计算出的数值之间的残差平方和最小化，其优化模型为：

数据拟合问题对应于一个最小化问题。

χ2的数值取决于所选择的模型函数f(xi|a)以及观测值yi中的随机误差εi。除了已经引入的变量以外，式(2.2)中还包含权值wi，每个权值应能反映相应观测值的可靠性。测量结果yi的不确定度越大，其权值wi就应该越低。比如，如果某个观测值看起来像是个异常值，那么其权值就应该非常接近或者等于零。但在大多数情况下，观测值的可靠性是无法预先获得的，所以有必要对权值进行估计。当对观测值进行加权时，一般称其为广义最小二乘方法或者是加权最小二乘方法，与之相反的是，如果每个观测值的权值都相等(即wi=w=常数)，就称其为常规最小二乘方法。

不同的参数向量a对于χ2的影响可以通过误差曲面来直观显现。求解χ2的最小值过程就是在误差曲面上寻求其全局最小值点。如果模型参数的个数超过两个，就将其对应的误差曲面称为误差超曲面。

非线性模型函数可以利用梯度下降法来求解，该方法通常要从某个起始点出发，通过迭代逐步到达全局最小值点。此时需要参数a的一个初始值，并且该值最好尽可能接近全局最小值点。

线性问题可以看成是非线性问题的一个子集，能够从多方面简化最小二乘拟合的求解过程。

还需要一个包含全部残差的列向量，残差是指观测值yi与利用当前参数计算出的模型函数值之差，如下式所示：

在计算刚开始的时候需要给出参数a的初始值。另外，对角加权矩阵W：

矩阵W是求解最小化问题的第3个组成部分(可选择也可不选)。该矩阵包含了对每个观测数据所赋予的权值，权值反映了观测值对于确定模型参数的重要性或是其自身的可靠性。如果所有观测值都具有相同影响，就可以将该矩阵设为单位矩阵，否则就应该选择合理的权值。

(本节，即第2节的余下内容只列出公式，不讲细节原理。在后面的【6.2 Gauss-Newton方法】一节中会给出详细推导过程)

综合上面各个要素可以计算下式：

此外，通过限制迭代次数的上限值可以避免陷入无限循环。

该正规矩阵显然是关于其主对角线对称的。

上述非线性最小二乘问题的求解过程仍然存在一些问题值得警惕。

1、第一个问题是如何合理给出参数a的初始值。显然，初始值应该尽可能接近全局最小值点，从而避免陷入局部极小值点，并且初始值需要使得迭代过程具有快速收敛性。为了实现这一目标，不妨利用一个能反映模型参数和观测值之间的简化关系式，即通过一些近似假设得到参数a的粗略估计，又或者是，在参数空间中利用网格搜索得到一个有希望的初始值。

2、参数的取值范围是第二个值得考虑的问题。如果参数的取值范围是已知的，就应该利用这一先验知识以监督迭代过程，从而避免其徒劳地收敛到非常大的值，又或者是收敛到无限小的值。可以修正误差(超)曲面的梯度，以使迭代过程不会进入到临界(错误)区域。

3、第三个潜在问题涉及到矩阵求逆运算。矩阵N=JTWJ中非常大的数值在其逆矩阵N-1中会变得很小。在很多应用中需要衡量该数值，以使其能与观测数据在同一量级上进行综合，否则的话，计算机的有限数值精度会产生难以预期的影响。

在理想的数据拟合条件下，实验条件与观测结果之间的函数关系是已知的，最小二乘拟合的目标就是要确定模型函数中的未知参数。

然而，在很多实际场景中，这个函数关系式可能是未知的，也就是说模型函数的解析形式无法获知。

下面举出三个例子。

针对5个采样信号，其条件值xi和采样值yi如下所示：

下图描绘了其拟合结果：

在上图中，模型函数的多项式阶数从0(M=1)增加至4(M=5)。零阶多项式用于确定观测值的均值，一阶多项式用于拟合直线，其他类推。从图中可以清楚发现，曲线越接近数据点，模型参数的个数就越多，但此时观测数据点之间的泛化性能就会不可靠。比如，当M=5时，模型函数会在前两个采样点之间发生振荡，这看起来不像是一个很充分的拟合。如果M相对于N的值过大，就称为“过拟合”。

仅仅将初等函数的叠加作为模型函数无法取得很好的预测性能。

下面将针对几个模型函数给出一些计算实例，以使得各种模型的物理含义和功能更加清晰。

直线拟合问题的参数个数从一个增加至两个，相应的模型函数为：

注意：在阅读完前面5.1、5.2两小节之后，读者对于最小二乘问题的求解过程应该可以明晰，因此下面仅给出雅可比矩阵。

最小二乘方法可以看成是更一般性的最大似然方法的一个特例。

其中：

由于P|a需要最大化，因此参数向量a的最优值应该使得χ2(a)取最小值。式(6.4)新增的系数“1/2”是为了便于后续处理，其对于拟合过程和最终结果并没有实质影响。系数1/σ2i起到加权的作用，并且依据每个观测值的可靠性或重要性来分配权值。

上面的数学推导是在观测误差服从正态(高斯)分布的假设条件下进行的。在单次测量可以获得多个观测值的情况下，需要利用观测向量y代替标量y，相应地，式(6.3)中的概率应该根据多维高斯概率密度函数来计算。

只要模型函数关于其中一个模型参数是非线性的，那么拟合该函数就需要迭代运算。向量a中包含的M个参数可以扩展成一个M维空间，该空间中的每一个点由χ2的大小来刻画，空间中所有点形成的实体称为超曲面。拟合算法是从某个初始值出发，试图寻求χ2(a)的全局最小值点。

下面分析如何利用该近似表达式确定超曲面上的迭代路径，即如何在当前迭代值的基础上寻求更优的参数值。

Gauss-Newton方法是一种在一维或者多维参数空间中寻求极小值点的方法。

当利用最小二乘方法进行数据拟合时，需要将观测值和模型函数之间的误差平方和最小化，即寻求超曲面χ2(a)上的极小值点。这可以通过目标函数的二阶泰勒级数展开来实现。

式(6.14)中的变量则需要被替换。首先，问题的自变量为a，于是有：

式(6.15)中的H称为Hessian矩阵，其中的元素为：

该矩阵中的元素可以利用导数乘法规则(uv)’=u’v+uv’来推导，根据式(6.17)，可得：

观察g和H中的元素(见式(6.17)和式(6.19))，可得：

J，r和W分别为：

请注意下面的内容。

矩阵Q包含了模型函数f(x|a)的二阶导数项，而该项在很多教科书中是被忽略的。有两个原因：

1、第一个原因是，对于线性问题，该项就等于0。

2、第二个原因是，式(6.20)中的乘法项yi-f(xi|a)可以近似看成是实验条件xi下的观测误差，只要参数向量a接近真实值，该误差就会或正或负，如果乘法项yi-f(xi|a)与f(xi|a)的二阶导数不相关，那么当对i进行累加求和时，这些项就会相互抵消。

若将矩阵Q忽略可得：

上述最小化过程称为Gauss-Newton方法，其主要优势在于能够快速收敛到邻近的极小值点。但是，能否成功收敛到极小值点取决于目标函数的曲率。若不做出一些针对性的处理，该算法可能会陷入鞍点，或者沿错误方向进行迭代。

不用担心，6.3节会提出梯度下降法。

尽管Gauss-Newton方法十分有效，但是当参数向量a=(a1,a2,…,aj,…,aM)T的初始值不能充分接近真实值时，该方法有可能达不到期望的极小值点。因为此时利用有限阶泰勒级数近似会对优化过程产生不利影响，比如，导致太大的参数调整步长或者是错误的迭代方向。

在这种情况下，更适合采用梯度下降方法，以保证严格沿“下坡”方向迭代。

于是：

上式(6.24)利用梯度向量的范数做归一化，因子参数γ决定了步长幅度，而向量g则仅仅控制迭代方向。

如果参数空间中的当前迭代值与最终的解相差甚远，梯度下降方法还是相当稳健的，只是有时难以估计γ的数值，通常它的值要保持相对较小。然而，当迭代过程接近目标函数的极小值点时，梯度下降方法就将会变得不太有利。如果没有对||g||做归一化(见式(6.24))，收敛过程将会变得十分缓慢，因为在这个区域内的梯度会很小。另一方面，鞍点也是阻止梯度下降方法成功收敛的另一不利因素。当按照式(6.24)归一化时，如果迭代已经越过极小值点，并且算法出现振荡现象，就必须要降低γ的数值。

上面一段就是为什么在极小值点附近使用Gauss-Newton方法代替梯度下降方法的原因。

Levenberg-Marquardt方法综合利用了梯度下降方法和Gauss-Newton方法两者的优势，它引入了一个阻尼因子。利用阻尼因子可将式(6.22)修改为如下形式：

阻尼因子μ取正数，以确保Δa是个下降方向。当μ取值较大时，所引入的项“μⅠ”将起主导作用，此时的解近似为：

上式与式(6.25)的梯度下降方法相对应。反之，非常小的μ值则会使迭代过程与Gauss-Newton方法相似，这有助于最终收敛到全局极小值点。

书名: 《数据拟合与不确定度 加权最小二乘及其推广的实用指南》 作者: (德)汤露·舒茨(Tilo Strutz)

END