\[ \def \v#1{{\boldsymbol{#1}}} \def \m#1{{\mathbf{#1}}} \def \argmax#1{\underset{#1}{\operatorname{argmax}}} \]
1801年1月1日,意大利天文學家Giuseppe Piazzi發現了穀神星。在连续观测了几十天後,穀神星湮沒在太陽的眩光之中。爲了再次找到它,24嵗的高斯憑藉現在被稱爲最小二乘法的方法,預測了穀神星的軌道。1801年12月31日,天文學家在高斯預測的軌道附近再度找回了穀神星。这就是最小二乘法在历史上首次登场的故事。
最小二乘(least squares)的含義是最小化誤差的平方和(該術語受到日文影響。其中的“二乘”即爲“平方”)。給定一系列樣本點,以最小二乘為目標,求模型參數的優化方法叫做最小二乘法。其中有綫性解的方法叫做綫性最小二乘法。
下面的小程序給出了一個例子。你可以通過鼠標點擊在方框中添加二維點。如果方框中有兩個以上的點,程序就會應用綫性最小二乘法,找到使誤差平方和最小的直綫。這個例子中,樣本點是我們手動添加的二維坐標,模型參數是直綫的係數。下文將會詳細介紹程序的計算原理。
1 最小二乘法求解直綫擬合問題
上面展示的是一個求解直綫擬合問題的程序。我們先形式化地描述這個問題,然後給出數學求解步驟。
問題:已知二維平面上的有一系列點 \((x_k, y_k)\),\(k=1,2,3,\cdots\)要求用最小二乘法擬合直線。
針對這個問題,根據設置的目標函數不同,解法也不同。下面將給出幾種不同的目標函數,説明其特點,然後給出對應的解法。
1.1 斜截式
例子 1 在這個例子中,我們設直綫方程為\(y=ax+b\),其中\(a\)為斜率,\(b\)為截距。這種類型的直綫方程被稱爲斜截式。
根據方程可列出如下式子: \[ \left(\begin{matrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \end{matrix} \right) \left(\begin{matrix} a \\ b \\ \end{matrix}\right) = \left(\begin{matrix} y_1\\ y_2\\ \vdots \end{matrix}\right) \] 該方程滿足\(\mathbf A \v x = \v b\)的形式,其中\(\mathbf A\in \mathbb R^{k\times 2}\),\(\v x=(a,b)^T\), \(\v b \in \mathbb R^{k\times 1}\).
解:
我們的目標為最小化 \[ \begin{aligned} L(\v x) &= ||\mathbf A\v x - \v b||^2 \\ &= \v x^T \mathbf A^T \mathbf A \v x - 2 \v x ^T \mathbf A^T \v b +\v b^T \v b \end{aligned} \] 令其導數為0 \[ \frac{\partial (L(\v x))}{\partial x} = 0 \] 得 \[ \begin{aligned} &\mathbf A^T \mathbf A \v x - \mathbf A^T \v b = 0 \\ \Rightarrow &\mathbf A^T \mathbf A \v x = \mathbf A^T \v b\\ \Rightarrow &\v x= (\mathbf A^T \mathbf A)^{-1} \mathbf A^T \v b \end{aligned} \]
這樣就得到了該問題的最優解。
1.2 一般式
例子 2 設直線方程爲\(ax+by+c=0\),則有 \[\left(\begin{matrix}
x_1 & y_1 & 1\\
x_2 & y_2 & 1\\
\vdots & \vdots & \vdots\end{matrix}\right)
\left(\begin{matrix}
a\\b\\c
\end{matrix}\right)=\v 0
\]
方程滿足\(\mathbf A\v{x}=\v 0\)的形式。其中\(\mathbf A\)是\(k\times 3\)的矩陣,\(\v x=(a, b, c)^T\).
解:
與開頭的方法不同,此方法采用一般式來表示直綫,因此能夠正常表示垂直於\(x\)軸的直綫。
根據給定的目標,求解\(\v x\)使得 \[||\mathbf A\v x||^2=(\mathbf A \v x)^T\mathbf A\v x=\v x^T \mathbf A^T \mathbf A \v x\] 的值最小。
這個問題顯然有一個解\(\v x=\v 0\),但這個解沒有意義,因爲它不能構成直線的係數,人們把這個解叫做“平凡解(trivial solution)”。
爲了排除平凡解,引入一個約束條件: \[\v x^T \v x=1\]
問題變爲約束條件下求最小值的問題,可以運用拉格朗日乘數法。令 \[L(\v x)=\v x^T\mathbf A^T\mathbf A\v x-\lambda(\v x^T\v x-1)\]
\(L(\v x)\)最小時,\[{\partial L(\v x)\over \partial \v x} = \v 0,\] 因此 \[\mathbf A^T\mathbf A\v x - \lambda \v x = \v 0\]
因此\(\lambda\)是\(\mathbf A^T\mathbf A\)的特徵值,\(\v x\)是\(\mathbf A^T \mathbf A\)對應\(\lambda\)的一個單位特徵向量,記\(\v x=\boldsymbol{e_\lambda}\)。代入得\(L(\boldsymbol{e_\lambda})=\lambda \boldsymbol{e_\lambda}^T \boldsymbol{e_\lambda}=\lambda\)。\(\lambda\)最小時,\(L(\v x)\)取得最小值。換言之,\(\mathbf A^T\mathbf A\)的最小特徵值對應的單位特徵向量即是所求的解。
求矩陣的特徵特徵向量有幾種方法。注意這裏的\(\mathbf A^T\mathbf A\)是實對稱矩陣,可以用Jacobi eigenvalue algorithm求特徵值。另一種常見的方法則是使用SVD分解。
1.3 用SVD分解求解最小二乘問題
根據SVD分解, \[\m A=\m U\m D\m V^T= \left[\begin{matrix} \boldsymbol{u_0} | \cdots | \boldsymbol{u_{p-1}} \end{matrix}\right] \left[\begin{matrix} \boldsymbol{\sigma_0} & & \\ & \ddots & \\ & & \boldsymbol{\sigma_{p-1}} \end{matrix}\right] \left[\begin{matrix} \boldsymbol{v_0^T}\\ \hline \vdots\\ \hline \boldsymbol{v_{p-1}^T} \end{matrix}\right] \] 注意到\(\m A^T\m A=(\m U\m D\m V^T)^T(\m U\m D\m V^T)=\m V\m D\m U^T\m U\m D\m V^T=\m V\m \Lambda \m V^T\) 其中,\(\m \Lambda\)的第i行的非0元素\(\lambda_i=\sigma_i^2\) 兩邊右乘\(\m V\)得到 \[(\m A^T\m A)\m V=\m V\m \Lambda\] 即是說,\(\m V\)的每一個列向量\(\boldsymbol{v_j}\)都是\(\m A^T\m A\)的特徵向量,對應的特徵值是\(\lambda_j\).
1.4 局限性
前文中我們分斜截式和一般式兩種情況討論了直綫擬合問題,並指出一般式具有能夠擬合任意方向直綫的優點。但還應當注意到兩種方法都存在有共同的局限性:其并非以最小化點到直綫的垂直距離為目標,即最小化 \[ \sum_i {|ax_i + by_i + c|\over \sqrt{a^2 + b^2}} \]
後來提出的“正交回歸”方法解決了這一問題。但正交回歸已經超出了本文要討論的範疇,這裏就不再繼續展開。
2 總結
為了解決文中所述的帶平凡解的最小二乘問題,本文先使用斜截式建立直綫方程並給出解,然後指出了斜截式的局限性。
而後本文嘗試了使用一般式表示直綫,並引入約束條件\(\v x^T \v x=1\),然後運用拉格朗日乘數法來求解最小值。
最後,文章討論了此類方法的局限性。
雖然本文圍繞著直綫擬合問題展開,但最小二乘法不僅可以用於直線的擬合,也可以推廣到圓、橢圓、二項式、平面……等各類模型,是一個非常實用的數學工具。