一看就會(廢)的最小二乘法的推導

一、 預備知識：方程組解的存在性及引入

  最小二乘法可以用來做函式的擬合或者求函式極值。在機器學習的迴歸模型中，我們經常使用最小二乘法。我們先舉一個小例子來走進最小二乘法。

  某次實驗得到了四個資料點\((x,y):(1,6)、(2,5)、(3,7)、(4,10)\) (下圖中紅色的點)。我們希望找出一條與這四個點最匹配的直線 \(y = \theta_{1} + \theta_{2}x\) ，即找出在某種"最佳情況"下能夠大致符合如下超定線性方程組的 \(\theta_{1}\) 和 \(\theta_{2}\) ，我們把四個點代入該直線方程可得:

\[\theta_{1} + 1 \theta_{2} = 6\\ \theta_{1}+2\theta_{2}=5\\ \theta_{1}+3\theta_{2}=7\\ \theta_{1}+4\theta_{2}=10 \]

我們要求的是\(\theta_{1}\)和\(\theta_{2}\) 兩個變數，但是這裡列出了四個方程組，我們是無法求解的。我們現在以向量空間的角度來解釋為何無解，以及最小二乘法如何處理這種無解的情況。

\[Ax = b\\ \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \\ \end{bmatrix} _{4\times2} \begin{bmatrix} \theta1\\ \theta2 \end{bmatrix}_{2\times1} = \begin{bmatrix} 6\\ 5\\ 7\\ 10 \end{bmatrix}_{4\times1} \]

我們將四個方程組成的方程組寫成矩陣形式。矩陣A代表係數，\(x\)即待求的引數，\(b\)是每個方程對應的值。

從線性代數的角度來看，要判斷\(Ax=b\)是否有解可以從向量空間角度來看。

這裡先給出向量空間的性質：

• 向量空間要求取其任意取兩個列向量\(v\)和\(w\)，\(v+w\)或者\(cv\)(c是一個常數)都要屬於該向量空間，並且任意列向量數乘的排列組合\(cv+dw+ek\)(c,d,e表示常數，v,w,k表示任意的列向量)也要屬於該向量空間。若滿足此條件，則可說明這是一個向量空間。

• 最大的向量空間是某個維度的全部空間，比如四維空間就是就是向量空間\(R_4\)的最大向量空間，最小向量空間是零空間，即由0組成的向量空間。

要看一個方程組有沒有解，我們要看矩陣A的列向量所有的排列組合是否能夠表示向量b的值。矩陣A有四個方程，即矩陣A由四個分量組成，所以它的向量空間記為\(R^4\),而列向量有兩個，列向量的所有排列組合組成的空間是列空間(column space)，列空間是\(R^4\)的子空間。對於矩陣\(A_{m\times n}\)，如果\(m>n\)，就像這個例子\(m=4,n=2\)，二維的列向量排列組合所得到的列向量空間撐死也只能充滿整個二維空間(如果有列的線性無關項就更少了)，而對於四維的向量空間\(R^4\)，這個例子中列向量空間是遠遠達不到的。那麼，當向量\(b\)的值在列向量空間之外時，就無法求解反之則有解。所以說，只要向量\(b\)是A的列向量的排列組合就能使\(x\)有解。

順便再從另外一個角度來思考，為什麼對於\(Ax=0\)來說是一定有解的呢？很簡單的啦，你發現列向量的排列組合肯定是能夠得到0向量的，只需要讓列向量和0進行數乘即可。所以0向量\(b\)一定是能由列向量的排列組合得到，因此它一定有解。

注意：\(m>n\)不一定就代表無解，有時候矩陣可能會有行的線性無關項，某行線性無關就是說某一行可以通過另外兩行進行乘、加操作後得到。\(m>n\)的方程組可以叫做超定方程組。而曲線擬合的最小二乘法要解決的問題，實際上就是以上超定方程組最小二乘解的問題。

我們把A的咧空間抽象成一個超平面，b可以抽象成一個向量。當對於無解的超定方程組時(即b不在超平面上)，我們可以在超平面上找一個最接近的向量b，假設這個最接近的向量是Ax*(其中x*就是最優的引數解)，那麼，向量b和向量Ax*之間的距離應該是最小的。所以可以通過最小化b和Ax*之間的距離來求出x*。

二、 最小二乘法的原理及解釋

  最小二乘法的原理還是很簡單的，形式如下：

\[目標函式 = \sum (觀測值 - 理論值)^2 \]

  觀測值(樣本值)就是我們的多組樣本，理論值(預測值)就是我們假設的擬合函式對應的值。目標函式也叫做損失函式，如上所說，我們的目標是得到使損失函式最小化時擬合函式的模型(即此時對應的引數)。舉一個最簡單的一元線性迴歸的例子，比如我們有m個樣本，每個樣本只有一個特徵：

\[(x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),...(x^{(m)},y^{(m)}) \]

我們選定一個用於擬合曲線的函式為：

\[h_{\theta}(x) = \theta_{0} + \theta_{1}x \]

我們的目標函式一般形式為：

\[J(\theta_{0},\theta_{1}) = \sum_{i=1}^m(y^{(i)}-h_{\theta}(x))^2\\=\sum_{i=1}^m(y^{(i)}-\theta_{0}-\theta_{1}x^{(i)})^2 \]

我們就是要求出當函式\(J（\theta_0,\theta_1）\)最小時的\(\theta_0\)和\(\theta1\)。

下面這張圖是第一部分的四個資料點(紅色)，以及用最小二乘法擬合的曲線(藍色的直線)。通過影像形象的表示就是最小二乘法要根據預先選擇的擬合函式(\(y=\theta_0+\theta_1 x\))找到一條曲線，讓所有的綠色線段的距離加起來的值最小。

下面我們會通過代數法和矩陣法來進行求解過程的推導。

三、用代數法求解最小二乘法

  代數解法就是分別對引數\(\theta_0,\theta_1\)分別求偏導數，並令其等於0，得到一個關於\(\theta_0\)和\(\theta_1\)的二元方程組。求解這個方程組就可以得到\(\theta_0\)和\(\theta_1\)的值，求導後得到以下方程組:

\[\frac{\partial E(\theta_0,\theta_1)}{ \partial \ \theta_0} =\sum_i^m(y^{(i)}-\theta_0 -\theta_1x^{(i)}) = 0 \tag{1} \]

\[\frac{\partial E(\theta_0,\theta_1)}{ \partial \ \theta_1} =\sum_i^m(y^{(i)}-\theta_0-\theta_1x^{(i)})x^{(i)} = 0 \tag{2} \]

(1)和(2)組成了一個二元一次方程組，可以求出\(\theta_0\)和\(\theta_1\)的值:

\[b = \theta_0 = \frac{1}{m} \sum_{i=1}^m(y^{(i)}-\theta_1x^{(i)})\\ \]

\[w = \frac{\sum_{i=1}^my_{i}(x^{(i)}- \overline x)}{\sum_{i=1}^m(x^{(i)})^2-\frac{1}{m}(\sum_{i=1}^mx^{(i)})^2} \tag{其中$\overline x = \frac{1}{m}\sum_{i=1}^mx^{(i)}$為$x$的均值} \]

這個方法很容易推導到多個樣本特徵的非線性擬合。原理和上面的一樣，都是用損失函式對各個引數求導為0，然後求解方程組得到引數值。

四、用矩陣法求解最小二乘法

  因為矩陣運算可以取代迴圈(也就是將其向量化)，可以讓矩陣法比代數法要更加的簡潔。

對於一個資料集，樣本由\(n\)個類似的，我們用最小二乘法來對\(w\)和\(b\)來進行估計。為了方便討論，我們把\(w\)和\(b\)進行合併乘一個向量 \(\hat w = (w_i;b)\)；（這裡為了比較好分辨，就不要再用\(\theta_0 \theta_1\)來表示了，用\(w,b\)）把。把資料集表示成一個\(M \times N+1\)大小的矩陣\(X\)，其中每行對應一條資料，該行的前N列對應於資料的N個屬性值（特徵值），最後一個元素恆置為1，目的是當\(X\hat w\)時令1與\(\hat w\)中的b相乘，則所得的每行都是$x_1^Tw_1 +b $的形式。

\[X = \left[ \begin{matrix} x_{11} & x_{12} & \cdots & x_{1n} & 1 \\ x_{21} & x_{22} & \cdots & x_{2n} & 1 \\ \vdots & \vdots & \ddots & \vdots &\vdots\\ x_{m1} & x_{m2} & \cdots &x_{mn} & 1\\ \end{matrix} \right] = \left[ \begin{matrix} x_1^T & 1 \\ x_2^T & 1 \\ \vdots &\vdots\\ x_m^T & 1\\ \end{matrix} \right] \]

再把標記(即資料集的\(y\)值)也寫成向量的形式\(y = (y_1;y_2;y_3;...;y_m)\),則有：

\[\hat w^* = \begin{equation} \mathop{\arg\min}_{\hat w} \end{equation} (y - X \hat w )^T(y - X\hat w). \]

令\(E_{\hat w} = (y - X \hat w)^T(y - X \hat w)\)，對\(\hat w\)求導得到：

\[\frac{\partial E_{\hat w}}{\partial \hat w} = 2 X^T(X\hat w - y) \]

令上面的式子為0即可得到\(\hat w\)最優解的閉式解。當\(X^TX\)為滿秩矩陣時，令其等於0可以得到，

\[\hat y^* = (X^TX)^{-1}X^Ty \]

其中\((X^TX)^{-1}\)是矩陣\((X^TX)\)的逆矩陣，因為只有滿秩矩陣才可逆。令\(\hat x_i = (x_i;1)\)，則最終學得的多元線性迴歸模型為 \(f(\hat x_i) = \hat x_i^T(X^TX)^{-1}X^Ty\)。

但是，由於現實生活中\(X^TX\)往往不是滿秩矩陣，例如有時候會遇到大量變數，數目甚至超過了樣本數目，導致\(X\)的列數多於行數，此時明顯不是滿秩。從線性方程組角度來思考就是，當因變數過多的時候(即方程個數小於未知數個數的方程組)，對於這種方程組我們叫做欠定方程組，該方程組的解不是唯一的，可能有無窮多組解。我們可以引入正則化來解決。這裡就先不具體去闡述什麼是正則化了。

五、 總結

此部分摘自 劉建平Pinard 部落格園對於最小二乘法的總結。

  從上面可以看出，最小二乘法適用簡潔高效，比梯度下降這樣的迭代法似乎方便很多。但是這裡我們就聊聊最小二乘法的侷限性。

  首先，最小二乘法需要計算\(X^TX\)的逆矩陣，有可能它的逆矩陣不存在，這樣就沒有辦法直接用最小二乘法了，此時梯度下降法仍然可以使用。當然，我們可以通過對樣本資料進行整理，去掉冗餘特徵。讓\(X^TX\)可逆(滿秩，行列式值不為0)，然後繼續使用最小二乘法。

  第二，當樣本特徵n非常的大的時候，計算\(X^TX\)的逆矩陣是一個非常耗時的工作（\(n\times n\)的矩陣求逆），甚至不可行。此時以梯度下降為代表的迭代法仍然可以使用。那這個\(n\)到底多大就不適合最小二乘法呢？如果你沒有很多的分散式大資料計算資源，建議超過10000個特徵就用迭代法吧。或者通過主成分分析降低特徵的維度後再用最小二乘法。

  第三，如果擬合函式不是線性的，這時無法使用最小二乘法，需要通過一些技巧轉化為線性才能使用，此時梯度下降仍然可以用。

  第四，講一些特殊情況。當樣本量m很少，小於特徵數n的時候，這時擬合方程是欠定的，常用的優化方法都無法去擬合資料。當樣本量m等於特徵數n的時候，用方程組求解就可以了。當m大於n時，擬合方程是超定的，也就是我們常用與最小二乘法的場景了。