最小二乘法求解直線擬合問題

本文瀏覽次數

发布日期

2022年1月4日

摘要
本文以直綫擬合問題為例,介紹了最小二乘法的計算方法。文章基於SVD分解和拉格朗日乘數法,給出了直綫擬合問題的閉式解。

1801年1月1日,意大利天文學家Giuseppe Piazzi發現了穀神星。在连续观测了几十天後,穀神星湮沒在太陽的眩光之中。爲了再次找到它,24嵗的高斯憑藉現在被稱爲最小二乘法的方法,預測了穀神星的軌道。1801年12月31日,天文學家在高斯預測的軌道附近再度找回了穀神星。这就是最小二乘法在历史上首次登场的故事。

最小二乘(least squares)的含義是最小化誤差的平方和(該術語受到日文影響。其中的“二乘”即爲“平方”)。給定一系列樣本點,以最小二乘為目標,求模型參數的優化方法叫做最小二乘法。其中有綫性解的方法叫做綫性最小二乘法

下面的小程序給出了一個例子。你可以通過鼠標點擊在方框中添加二維點。如果方框中有兩個以上的點,程序就會應用綫性最小二乘法,找到使誤差平方和最小的直綫。這個例子中,樣本點是我們手動添加的二維坐標,模型參數是直綫的係數。下文將會詳細介紹程序的計算原理。

瀏覽器不支持!!

直綫擬合問題

上面展示的是一個求解直綫擬合問題的程序。我們先形式化地描述這個問題,然後給出數學求解步驟。

問題:已知二維平面上的有一系列點 \((x_k, y_k)\)\(k=1,2,3,\cdots\)要求用最小二乘法擬合直線。

解: 設直線方程爲\(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)=0 \]
方程滿足\(A\v{x}=0\)的形式。其中\(A\)\(k\times 3\)的矩陣,\(x=(a, b, c)^T\).

所謂最小二乘法,即是求解\(x\)使得 \[||Ax||^2=(Ax)^T Ax=x^T A^T A x\] 的值最小。

這個問題顯然有一個解\(x=0\),但這個解沒有意義,因爲它不能構成直線的係數,人們把這個解叫做“平凡解(trivial solution)”。

爲了排除平凡解,引入一個約束條件\[x^Tx=1\]

問題變爲約束條件下求最小值的問題,可以運用拉格朗日乘數法。令 \[L(x)=x^TA^TAx-\lambda(x^Tx-1)\]

\(L(x)\)最小時,\[{\partial L(x)\over \partial x} = 0,\] 因此 \[A^TAx - \lambda x = 0\]

因此\(\lambda\)\(A^TA\)的特徵值,\(x\)\(A^T A\)對應\(\lambda\)的一個單位特徵向量,記\(x=\boldsymbol{e_\lambda}\)。代入得\(L(\boldsymbol{e_\lambda})=\lambda \boldsymbol{e_\lambda}^T \boldsymbol{e_\lambda}=\lambda\)\(\lambda\)最小時,\(L(x)\)取得最小值。換言之,\(A^TA\)的最小特徵值對應的單位特徵向量即是所求的解。

求矩陣的特徵特徵向量有幾種方法。注意這裏的\(A^T A\)是實對稱矩陣,可以用Jacobi eigenvalue algorithm求特徵值。另一種常見的方法則是使用SVD分解。

用SVD分解求解最小二乘問題

根據SVD分解, \[A=UDV^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] \] 注意到\(A^TA=(UDV^T)^T(UDV^T)=VDU^TUDV^T=V\Lambda V^T\) 其中,\(\Lambda\)的第i行的非0元素\(\lambda_i=\sigma_i^2\) 兩邊右乘\(V\)得到 \[(A^TA)V=V\Lambda\] 即是說,\(V\)的每一個列向量\(\boldsymbol{v_j}\)都是\(A^TA\)的特徵向量,對應的特徵值是\(\lambda_j\).

總結

為了解決文中所述的帶平凡解的最小二乘問題,本文先引入約束條件\(x^Tx=1\),然後運用拉格朗日乘數法來求解最小值。實際上,還有其它種形式的約束條件,例如,可以約定\(b=1\),但這樣一來, 直線方程從一般式退化成斜截式,就無法表示平行於y軸的直線了。儘管有缺陷,這種形式的約束條件也還是很常用。

若是採用斜截式,經過求導等數學推導步驟,一樣可以求出直線的參數。但不論採用何種約束條件,這兩種方法都是以最小化誤差的平方和為目的的,都可以叫做最小二乘法。

最小二乘法不僅可以用於直線的擬合,也可以推廣到圓、橢圓、二項式、平面……等各類模型的擬合。

參考資料

By @執迷 in
Tags : #最小二乘法,