最小平方法

維基百科,自由的百科全書
前往: 導覽搜尋
線性代數
\mathbf{A} = \begin{bmatrix}
1 & 2 \\
3 & 4 \end{bmatrix}
向量 · 矩陣  · 行列式  · 線性空間

最小平方法(又稱最小平方法)是一種數學優化技術。它通過最小化誤差的平方和尋找數據的最佳函數匹配。

利用最小平方法可以簡便地求得未知的數據,並使得這些求得的數據與實際數據之間誤差的平方和為最小。

最小平方法還可用於曲線擬合

其他一些優化問題也可通過最小化能量或最大化用最小平方法來表達。

示例[編輯]

數據點(紅色)、使用最小平方法求得的最佳解(藍色)、誤差(綠色)。

某次實驗得到了四個數據點 (x, y)(1, 6)(2, 5)(3, 7)(4, 10)(右圖中紅色的點)。我們希望找出一條和這四個點最匹配的直線 y=\beta_1+\beta_2 x,即找出在某種「最佳情況」下能夠大致符合如下超定線性方程組的 \beta_1\beta_2

\begin{alignat}{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{alignat}

最小平方法採用的手段是盡量使得等號兩邊的方差最小,也就是找出這個函數的最小值:

\begin{align}
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{align}

最小值可以通過對 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 是最佳的。

簡介[編輯]

高斯

1801年,義大利天文學家朱賽普·皮亞齊發現了第一顆小行星穀神星。經過40天的跟蹤觀測後,由於穀神星運行至太陽背後,使得皮亞齊失去了穀神星的位置。隨後全世界的科學家利用皮亞齊的觀測數據開始尋找穀神星,但是根據大多數人計算的結果來尋找穀神星都沒有結果。時年24歲的高斯也計算了穀神星的軌道。奧地利天文學家海因里希·奧爾伯斯根據高斯計算出來的軌道重新發現了穀神星。

高斯使用的最小平方法的方法發表於1809年他的著作《天體運動論》中,而法國科學家勒讓德於1806年獨立發現「最小平方法」,但因不為時人所知而默默無聞。兩人曾為誰最早創立最小平方法原理發生爭執。

1829年,高斯提供了最小平方法的優化效果強於其他方法的證明,見高斯-馬爾可夫定理

方法[編輯]

人們對由某一變數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} \ .

最小化問題的精度,依賴於所選擇的函數模型。

線性函數模型[編輯]

典型的一類函數模型是線性函數模型。最簡單的線性式是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_iy_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 = b0 + b1t 的例子[編輯]

隨機選定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 18.41 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.

然後確定b1

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。並由下式得到常數項b0

b_0 = \bar y - b_1 \bar t = 18{.}41 - 0{.}1612 \cdot 167{.}8 = -8{.}6394
\;,

在這裡隨機理論不加闡述。可以看出點的擬合非常好,長度寬度的相關性大約為96.03%。 利用Matlab得到擬合直線: 最小平方法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}

通常人們將tij記作數據矩陣 A,參數bj記做參數向量b,觀測值yi記作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{align}
\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{align}

所以\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 (English). 

外部連結[編輯]