最小平方法

最小平方法（英語：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）基本呈線性關係。散點圖如下：

以下圖表列出了各戰艦的數據，隨後步驟是採用最小平方法確定兩變量間的線性關係。

i	t_i	y_i	t_i*	y_i*	t_iy_i	t_it_i	y_iy_i
編號	長度 (m)	寬度 (m)	t_i - t	y_i - y
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".

最小平方法

目次