1 什么是“最小二乘”

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

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

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

因此，最好能够建立尽可能多的方程(包含多对y_i和x_i)，从而获得更好、更可靠的统计结果，以更加准确地估计参数a。

假设得到一组观测值y_i及其相应的条件x_i，数据拟合的目标就是要利用最小二乘方法，将观测值与基于模型函数计算出的数值之间的残差平方和最小化，其优化模型为：

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

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

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

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

2 求解最小化问题的一般性算法

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

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

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

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

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

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

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

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

3 需要注意的问题

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

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

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

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

4 模型未知时的曲线拟合问题

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

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

下面举出三个例子。

4.1 例1

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

下图描绘了其拟合结果：

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

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

4.2 例2

4.3 例3

5 模型已知时的曲线拟合问题

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

5.1 常数拟合

5.2 直线拟合

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

5.3 余弦函数拟合

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

6 拟合非线性模型函数

6.1 误差曲面的近似

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




其中：

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

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

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

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

6.2 Gauss-Newton方法

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

当利用最小二乘方法进行数据拟合时，需要将观测值和模型函数之间的误差平方和最小化，即寻求超曲面χ²(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)中的乘法项y_i-f(x_i|a)可以近似看成是实验条件x_i下的观测误差，只要参数向量a接近真实值，该误差就会或正或负，如果乘法项y_i-f(x_i|a)与f(x_i|a)的二阶导数不相关，那么当对i进行累加求和时，这些项就会相互抵消。

若将矩阵Q忽略可得：

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

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

6.3 梯度下降方法

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

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

于是：

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

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

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

6.4 Levenberg-Marquardt方法

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

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

上式与式(6.25)的梯度下降方法相对应。

反之，非常小的μ值则会使迭代过程与Gauss-Newton方法相似，这有助于最终收敛到全局极小值点。

7 参考书籍

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

END

一、最小二乘法的思路

最小二乘法（Least Squares）是回归分析中的一种标准方法，它是当方程数量大于未知数个数时，利用数据点构建的方程组，对未知参数进行一种近似估计的方法。之所以叫做“最小二乘”，是因为利用的优化项是由所有数据点与模型观测点残差的平方和构成的，通过极小化残差的平方和，达到一种从整体上最“接近”实际观测数据的模型参数。

二、最小二乘法的公式推导

以最简单的线性最小二乘法拟合为例。设最终需要得到的拟合直线公式为

其中，，为待估计的参数。

设我们拥有的一组数据的观测点对为，一共有对数据点，则将这些数据点分别代入方程，可以得到一共个方程构成的方程组，通过极小化实际观测数据与模型观测数据的残差平方和：

将分别对参数求偏导，并令偏导等于零，可得

化简后得到

在这里，由于都是对下标从求和，所以将求和上小标均省略。

对参数求偏导，可得

化简后得到

利用方程组成的关于未知参数的二元一次方程组，联立求解，可得

其中上的横线代表对观测数据点求平均值。

至此，我们推导出了线性最小二乘法对于最简单两个直线参数估计情形的计算公式。

普通最小二乘法简介

总体最小二乘法简介