最小二乘法

最小二乘法（英语：least squares method），又称最小平方法，是一种数学优化建模方法。它通过最小化误差的平方和寻找数据的最佳函数匹配。

使用二次函数拟合一组数据点的结果

利用最小二乘法可以简便的求得未知的数据，并使得求得的数据与实际数据之间误差的平方和为最小。

“最小二乘法”是对线性方程组，即方程个数比未知数更多的方程组，以回归分析求得近似解的标准方法。在这整个解决方案中，最小二乘法演算为每一方程式的结果中，将残差平方和的总和最小化。

最重要的应用是在曲线拟合上。最小二乘所涵义的最佳拟合，即残差（残差为：观测值与模型提供的拟合值之间的差距）平方总和的最小化。当问题在自变量（x变量）有重大不确定性时，那么使用简易回归和最小二乘法会发生问题；在这种情况下，须另外考虑变量-误差-拟合模型所需的方法，而不是最小二乘法。

最小二乘问题分为两种：线性或普通的最小二乘法，和非线性的最小二乘法，取决于在所有未知数中的残差是否为线性。线性的最小二乘问题发生在统计回归分析中；它有一个封闭形式的解决方案。非线性的问题通常经由迭代细致化来解决；在每次迭代中，系统由线性近似，因此在这两种情况下核心演算是相同的。

最小二乘法所得出的多项式，即以拟合曲线的函数来描述自变量与预计因变量的方差关系。

当观测值来自指数族且满足轻度条件时，最小二乘估计和最大似然估计是相同的。最小二乘法也能从动差法得出。

以下讨论大多是以线性函数形式来表示，但对于更广泛的函数族，最小二乘法也是有效和实用的。此外，迭代地将局部的二次近似应用于或然性（借由费希尔信息），最小二乘法可用于拟合广义线性模型。

最小二乘法通常归功于高斯（Carl Friedrich Gauss，1795），但最小二乘法是由阿德里安-马里·勒让德（Adrien-Marie Legendre）首先发表的。

历史背景

最小二乘法发展于天文学和大地测量学领域，科学家和数学家尝试为大航海探索时期的海洋航行挑战提供解决方案。准确描述天体的行为是船舰在大海洋上航行的关键，水手不能再依靠陆上目标导航作航行。

这个方法是在十八世纪期间一些进步的集大成：

• 不同观测值的组合是真实值的最佳估计；多次观测会减少误差而不是增加，也许在1722年由Roger Cotes首先阐明。

• 在相同条件下采取的不同观察结果，与只尝试记录一次最精确的观察结果是对立的。这个方法被称为平均值方法。托马斯·马耶尔（Tobias Mayer）在1750年研究月球的天平动时，特别使用这种方法，而拉普拉斯（Pierre-Simon Laplace）在1788年他的工作成果中以此解释木星和土星的运动差异。

• 在不同条件下进行的不同观测值组合。该方法被称为最小绝对偏差法，出现在Roger Joseph Boscovich在1757年他对地球形体的著名作品，而拉普拉斯在1799年也表示了同样的问题。

• 评定对误差达到最小的解决方案标准，拉普拉斯指明了误差的概率密度的数学形式，并定义了误差最小化的估计方法。为此，拉普拉斯使用了一双边对称的指数分布，现在称为拉普拉斯分布作为误差分布的模型，并将绝对偏差之和作为估计误差。他认为这是他最简单的假设，他期待得出算术平均值而成为最佳的估计。可相反地，他的估计是后验中位数。

最小二乘法

高斯

1801年，意大利天文学家朱塞普·皮亚齐发现了第一颗小行星谷神星。经过40天的追踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。当年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥伯斯根据高斯计算出来的轨道重新发现了谷神星。

高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中，而法国科学家勒让德于1806年独立发现“最小二乘法”，但因不为世人所知而默默无闻。两人曾为谁最早创立最小二乘法原理发生争执。

1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，见高斯-马尔可夫定理。

示例

数据点（红色）、使用最小二乘法求得的最佳解（蓝色）、误差（绿色）。

某次实验得到了四个数据点 $(x,y)$：$(1,6)$、$(2,5)$、$(3,7)$、$(4,10)$（右图红色的点）。我们希望找出一条和这四个点最匹配的直线 $y=\beta _{2}x+\beta _{1}$，即找出在某种“最佳情况”下能够大致符合如下超定线性方程组的 $\beta _{1}$ 和 $\beta _{2}$：

${\begin{alignedat}{4}\beta _{1}+1\beta _{2}&&\;=\;&&6&\\\beta _{1}+2\beta _{2}&&\;=\;&&5&\\\beta _{1}+3\beta _{2}&&\;=\;&&7&\\\beta _{1}+4\beta _{2}&&\;=\;&&10&\\\end{alignedat}}$

最小二乘法采用的方法是尽量使得等号两边差的平方最小，也就是找出这个函数的最小值：

${\begin{aligned}S(\beta _{1},\beta _{2})=&\left[6-(\beta _{1}+1\beta _{2})\right]^{2}+\left[5-(\beta _{1}+2\beta _{2})\right]^{2}\\&+\left[7-(\beta _{1}+3\beta _{2})\right]^{2}+\left[10-(\beta _{1}+4\beta _{2})\right]^{2}.\\\end{aligned}}$

最小值可以通过对 $S(\beta _{1},\beta _{2})$ 分别求 $\beta _{1}$ 和 $\beta _{2}$ 的偏导数，然后使他们等于零得到。

${\frac {\partial S}{\partial \beta _{1}}}=0=8\beta _{1}+20\beta _{2}-56$

${\frac {\partial S}{\partial \beta _{2}}}=0=20\beta _{1}+60\beta _{2}-154.$

如此就得到了一个只有两个未知数的方程组，很容易就可以解出：

$\beta _{1}=3.5$

$\beta _{2}=1.4$

也就是说直线 $y=3.5+1.4x$ 是最佳的。

方法

人们对由某一变量$t$ 或多个变量$t_{1}$……$t_{n}$ 构成的相关变量$y$感兴趣。如弹簧的形变与所用的力相关，一个企业的盈利与其营业额，投资收益和原始资本有关。为了得到这些变量同$y$之间的关系，便用不相关变量去构建$y$，使用如下函数模型

$y_{m}=f(t_{1},\dots ,t_{q};b_{1},\dots ,b_{p})$,

$q$个独立变量或$p$个系数去拟合。

通常人们将一个可能的、对不相关变量t的构成都无困难的函数类型称作函数模型（如抛物线函数或指数函数）。参数b是为了使所选择的函数模型同观测值y相匹配。（如在测量弹簧形变时，必须将所用的力与弹簧的膨胀系数联系起来）。其目标是合适地选择参数，使函数模型最好的拟合观测值。一般情况下，观测值远多于所选择的参数。

其次的问题是怎样判断不同拟合的质量。高斯和勒让德的方法是，假设测量误差的平均值为0。令每一个测量误差对应一个变量并与其它测量误差不相关（随机无关）。人们假设，在测量误差中绝对不含系统误差，它们应该是纯偶然误差(有固定的方差)，围绕真值波动。除此之外，测量误差符合正态分布，这保证了偏差值在最后的结果y上忽略不计。

确定拟合的标准应该被重视，并小心选择，较大误差的测量值应被赋予较小的权。并建立如下规则：被选择的参数，应该使算出的函数曲线与观测值之差的平方和最小。用函数表示为：

$\min _{\vec {b}}{\sum _{i=1}^{n}(y_{m}-y_{i})^{2}}.$

用欧几里得度量表达为：

$\min _{\vec {b}}\|{\vec {y}}_{m}({\vec {b}})-{\vec {y}}\|_{2}^{2}\ .$

又因为$\|{\vec {y}}_{m}({\vec {b}})-{\vec {y}}\|_{2}$ ≥0,

所以也可以表示为$\min _{\vec {b}}\|{\vec {y}}_{m}({\vec {b}})-{\vec {y}}\|_{2}\ .$

最小化问题的精度，依赖于所选择的函数模型。

线性函数模型

典型的一类函数模型是线性函数模型。最简单的线性式是$y=b_{0}+b_{1}t$，写成矩阵式，为

$\min _{b_{0},b_{1}}\left\|{\begin{pmatrix}1&t_{1}\\\vdots &\vdots \\1&t_{n}\end{pmatrix}}{\begin{pmatrix}b_{0}\\b_{1}\end{pmatrix}}-{\begin{pmatrix}y_{1}\\\vdots \\y_{n}\end{pmatrix}}\right\|_{2}=\min _{b}\|Ab-Y\|_{2}.$

直接给出该式的参数解：

$b_{1}={\frac {\sum _{i=1}^{n}t_{i}y_{i}-n\cdot {\bar {t}}{\bar {y}}}{\sum _{i=1}^{n}t_{i}^{2}-n\cdot ({\bar {t}})^{2}}}$ 和 $b_{0}={\bar {y}}-b_{1}{\bar {t}}$

其中${\bar {t}}={\frac {1}{n}}\sum _{i=1}^{n}t_{i}$，为t值的算术平均值。也可解得如下形式：

$b_{1}={\frac {\sum _{i=1}^{n}(t_{i}-{\bar {t}})(y_{i}-{\bar {y}})}{\sum _{i=1}^{n}(t_{i}-{\bar {t}})^{2}}}$

简单线性模型 y = b₀ + b₁t 的例子

随机选定10艘战舰，并分析它们的长度与宽度，寻找它们长度与宽度之间的关系。由下面的描点图可以直观地看出，一艘战舰的长度（t）与宽度（y）基本呈线性关系。散点图如下：

以下图表列出了各战舰的数据，随后步骤是采用最小二乘法确定两变量间的线性关系。

编号	长度 (m)	宽度 (m)	ti - t	yi - y
i	ti	yi	ti*	yi*	ti*yi*	ti*ti*	yi*yi*
1	208	21.6	40.2	3.19	128.238	1616.04	10.1761
2	152	15.5	-15.8	-2.91	45.978	249.64	8.4681
3	113	10.4	-54.8	-8.01	438.948	3003.04	64.1601
4	227	31.0	59.2	12.59	745.328	3504.64	158.5081
5	137	13.0	-30.8	-5.41	166.628	948.64	29.2681
6	238	32.4	70.2	13.99	982.098	4928.04	195.7201
7	178	19.0	10.2	0.59	6.018	104.04	0.3481
8	104	10.4	-63.8	-8.01	511.038	4070.44	64.1601
9	191	19.0	23.2	0.59	13.688	538.24	0.3481
10	130	11.8	-37.8	-6.61	249.858	1428.84	43.6921
总和（Σ）	1678	184.1	0.0	0.00	3287.820	20391.60	574.8490

仿照上面给出的例子

${\bar {t}}={\frac {\sum _{i=1}^{n}t_{i}}{n}}={\frac {1678}{10}}=167{.}8$ 并得到相应的${\bar {y}}=18{.}41$.

然后确定b₁

$b_{1}={\frac {\sum _{i=1}^{n}(t_{i}-{\bar {t}})(y_{i}-{\bar {y}})}{\sum _{i=1}^{n}(t_{i}-{\bar {t}})^{2}}}$

$={\frac {3287{.}820}{20391{.}60}}=0{.}1612\;,$

可以看出，战舰的长度每变化1m，相对应的宽度便要变化16cm。并由下式得到常数项b₀：

$b_{0}={\bar {y}}-b_{1}{\bar {t}}=18{.}41-0{.}1612\cdot 167{.}8=-8{.}6394\;,$

在这里随机理论不加阐述。可以看出点的拟合非常好，长度和宽度的相关性大约为96.03％。 利用Matlab得到拟合直线：

一般线性情况

若含有更多不相关模型变量$t_{1},...,t_{q}$，可如组成线性函数的形式

$y(t_{1},\dots ,t_{q};b_{0},b_{1},\dots ,b_{q})=b_{0}+b_{1}t_{1}+\cdots +b_{q}t_{q}$

即线性方程组

${\begin{matrix}b_{0}+b_{1}t_{11}+\cdots +b_{j}t_{1j}+\cdots +b_{q}t_{1q}=y_{1}\\b_{0}+b_{1}t_{21}+\cdots +b_{j}t_{2j}+\cdots +b_{q}t_{2q}=y_{2}\\\vdots \\b_{0}+b_{1}t_{i1}+\cdots +b_{j}t_{ij}+\cdots +b_{q}t_{iq}=y_{i}\\\vdots \\b_{0}+b_{1}t_{n1}+\cdots +b_{j}t_{nj}+\cdots +b_{q}t_{nq}=y_{n}\end{matrix}}$

通常人们将t_ij记作数据矩阵 A，参数b_j记做参数向量b，观测值y_i记作Y，则线性方程组又可写成：

${\begin{pmatrix}1&t_{11}&\cdots &t_{1j}\cdots &t_{1q}\\1&t_{21}&\cdots &t_{2j}\cdots &t_{2q}\\\vdots \\1&t_{i1}&\cdots &t_{ij}\cdots &t_{iq}\\\vdots \\1&t_{n1}&\cdots &t_{nj}\cdots &t_{nq}\end{pmatrix}}\cdot {\begin{pmatrix}b_{0}\\b_{1}\\b_{2}\\\vdots \\b_{j}\\\vdots \\b_{q}\end{pmatrix}}={\begin{pmatrix}y_{1}\\y_{2}\\\vdots \\y_{i}\\\vdots \\y_{n}\end{pmatrix}}$ 即 $Ab=Y$

上述方程运用最小二乘法导出为线性平方差计算的形式为：

$\min _{b}\|Ab-Y\|_{2}$。

最小二乘法的解

$\min _{b}\left\|{\boldsymbol {Ab}}-{\boldsymbol {Y}}\right\|_{2},{\boldsymbol {A}}\in \mathbf {C} ^{n\times m},{\boldsymbol {Y}}\in \mathbf {C} ^{n}$

的特解为A的广义逆矩阵与Y的乘积，这同时也是二范数极小的解，其通解为特解加上A的零空间。证明如下：

先将Y拆成A的值域及其正交补两部分

${\boldsymbol {Y}}={\boldsymbol {Y}}_{1}+{\boldsymbol {Y}}_{2}$

${\boldsymbol {Y}}_{1}={\boldsymbol {A}}{\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}\in R\left({\boldsymbol {A}}\right)$

${\boldsymbol {Y}}_{2}=\left({\boldsymbol {I}}-{\boldsymbol {A}}{\boldsymbol {A}}^{\dagger }\right){\boldsymbol {Y}}\in R\left({\boldsymbol {A}}\right)^{\bot }$

所以${\boldsymbol {Ab}}-{\boldsymbol {Y}}_{1}\in R\left({\boldsymbol {A}}\right)$，可得

$\left\|{\boldsymbol {Ab}}-{\boldsymbol {Y}}\right\|^{2}=\left\|{\boldsymbol {Ab}}-{\boldsymbol {Y}}_{1}+\left(-{\boldsymbol {Y}}_{2}\right)\right\|^{2}=\left\|{\boldsymbol {Ab}}-{\boldsymbol {Y}}_{1}\right\|^{2}+\left\|{\boldsymbol {Y}}_{2}\right\|^{2}$

故当且仅当${\boldsymbol {b}}$是${\boldsymbol {Ab}}={\boldsymbol {Y}}_{1}={\boldsymbol {A}}{\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}$解时，${\boldsymbol {b}}$即为最小二乘解，即${\boldsymbol {b}}={\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}={\left({{{\mathbf {A} }^{H}}{\mathbf {A} }}\right)^{-1}}{{\mathbf {A} }^{H}}{\mathbf {Y} }$。

又因为

$N\left({\boldsymbol {A}}\right)=N\left({\boldsymbol {A}}^{\dagger }{\boldsymbol {A}}\right)=R\left({\boldsymbol {I}}-{\boldsymbol {A}}^{\dagger }{\boldsymbol {A}}\right)=\left\{\left({\boldsymbol {I}}-{\boldsymbol {A}}^{\dagger }{\boldsymbol {A}}\right){\boldsymbol {h}}:{\boldsymbol {h}}\in \mathbf {C} ^{n}\right\}$

故${\boldsymbol {Ab}}={\boldsymbol {A}}{\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}$的通解为

${\boldsymbol {b}}={\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}+\left({\boldsymbol {I}}-{\boldsymbol {A}}^{\dagger }{\boldsymbol {A}}\right){\boldsymbol {h}}:{\boldsymbol {h}}\in \mathbf {C} ^{n}$

因为

${\begin{aligned}\left\|{\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}\right\|^{2}&<\left\|{\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}\right\|^{2}+\left\|\left({\boldsymbol {I}}-{\boldsymbol {A}}^{\dagger }{\boldsymbol {A}}\right){\boldsymbol {h}}\right\|^{2}\\&=\left\|{\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}+\left({\boldsymbol {I}}-{\boldsymbol {A}}^{\dagger }{\boldsymbol {A}}\right){\boldsymbol {h}}\right\|^{2},\left({\boldsymbol {I}}-{\boldsymbol {A}}^{\dagger }{\boldsymbol {A}}\right){\boldsymbol {h}}\neq {\boldsymbol {0}}\\\end{aligned}}$

所以${\boldsymbol {A}}^{\dagger }{\boldsymbol {Y}}$又是二范数极小的最小二乘解。

参考文献

书籍

• Wang Guorong; Wei Yimin; Qiao SanZheng. Equation Solving Generalized Inverses. Generalized Inverses: Theory and Computations. Beijing: Science Press. 2004: 6. ISBN 7-03-012437-5 （英语）. 
• Mansfield Merriman. "A List of Writings Relating to the Method of Least Squares". 

外部链接

• 3种统计学点估计的理论推演:动差法,最小二乘法,最大似然估计法 （页面存档备份，存于互联网档案馆）
• http://www.physics.csbsju.edu/stats/least_squares.html （页面存档备份，存于互联网档案馆）
• https://web.archive.org/web/20200218162523/http://www.zunzun.com/
• http://www.orbitals.com/self/least/least.htm （页面存档备份，存于互联网档案馆）
• 最小二乘法. PlanetMath.