機器學習基礎——推導線性回歸公式

在之前的文章當中，我們介紹過了簡單的朴素貝葉斯分類模型，介紹過最小二乘法，所以這期文章我們順水推舟，來講講線性回歸模型。

線性回歸的本質其實是一種統計學當中的回歸分析方法，考察的是自變量和因變量之間的線性關聯。后來也許是建模的過程和模型訓練的方式和機器學習的理念比較接近，所以近年來，這個模型被歸入到了機器學習的領域當中。然而，不管它屬於哪個領域，整個模型的思想並沒有發生變化。我們只要有所了解即可。

模型概念

線性回歸的定義非常簡單，它最簡單的形式其實就是一元一次方程組。比如，我們有如下式子：

\[y=wx+b \]

我們知道若干的x和y，要求w和b。解的方法很簡單，我們通過消元法，就可以很容易求出來w和b。

我們針對以上的式子做兩個變形，第一個變形是我們的自變量x不再是一個單值，而是一個m * n的矩陣。m表示樣本數，n表示特征數，我們寫成X。X矩陣的每一行是一個n維的行向量，它代表一個樣本。它的系數W也不再是一個值，而是一個n * 1的列向量，它的每一維代表一個樣本當中這一維的權重。我們把上面的公式寫成矩陣相乘的形式：

\[Y=XW+b \]

式子里的Y、X和W分別是m * 1, m * n和n * 1的矩陣。

這里有兩點要注意，第一點是這里的b我們可以當做是一個浮點數的參數，但是實際上它也是一個m * 1的矩陣（列向量）。但即使我們用的是浮點數也沒關系，因為在我們實現模型的時候，numpy或者TensorFlow或者是其他的框架會自動地使用廣播將它轉化成向量來做加法。

第二點是這里的X寫在了W的前面，這也是為了矩陣乘法計算方便。當然我們也可以將X和W都轉置，寫成WX，但這樣得到的結果是一個1 * m的行向量，如果要和Y進行比較，那么還需要再進行一次轉置。所以為了簡便，我們對調了X和W的順序。所以大家不要覺得疑惑，明明是WX+b怎么寫出來就成了XW+b了。

我們把式子列出來之后，目標就很明確了，就是要通過計算求到一個W和b使得式子成立。但是在此之前，我們需要先明確一點：在實際的工程應用場景當中，是不可能找到W和b使得XW+b恰好和Y完全相等的。因為真實的場景當中數據都存在誤差，所以精確的解是不存在的，我們只能退而求其次，追求盡可能精確的解。

最小二乘法與均方差

在之前的文章當中我們介紹過最小二乘法，遺忘的同學可以點擊下方鏈接回顧一下。

概率統計——期望、方差與最小二乘法

在機器學習的過程當中，模型不是直接達到最佳的，而是通過一步一步的迭代，效果逐漸提高，最終收斂不再劇烈變化。我們明白了這個過程，就能理解，在學習的過程當中，我們需要一個量化的指標來衡量模型當前學習到的能力。就好像學生在上學的時候需要考試來測試學生的能力一樣，我們也需要一個指標來測試模型的能力。

對於回歸模型而言，預測的目標是一個具體的值。顯然這個預測值和真實值越接近越好。我們假設預測值是\(\hat{y}\)，真實值是y，顯然應該是\(|y-\hat{y}|\)越小越好。

但是絕對值的計算非常麻煩，也不方便求導，所以我們通常會將它平方，即：\((y-\hat{y})^2\)最小。對於m個樣本而言，我們希望它們的平方和盡量小：\(\sum_{i=1}^m(y_i - \hat{y_i})^2\)。

這個式子和我們之前介紹的方差非常相似，只不過在方差當中減的是期望值，而在這里我們減的是真實值。所以這個平方差也有一個類似的名稱，叫做均方差。

方差反應的是變量在期望值附近的震盪程度，同樣的，均方差反應的是模型預測值距離真實值的震盪程度。

尋找最佳的參數來使得均方差盡量小，就是最小二乘法。

推導過程

到這里，我們已經搞清楚了模型的公式，也寫出了優化的目標，已經非常接近了，只剩下最后一步，就是優化這個目標的方法。

如果我們觀察一下均方差，我們把它寫全：\((Y-(XW+b))^2\)，我們將W視作變量的話，這其實是一個廣義的二次函數。二次函數怎么求最小值？當然是求導了。

在求導之前，我們先對均方差做一個簡單的變形：我們想辦法把b處理掉，讓式子盡可能簡潔。

首先，我們在X當中增加一列1，也就是將X變成m * (n+1)的矩陣，它的第一列是常數1，新的矩陣寫成\(X_{new}\)

同樣，我們在W中也增加一行，它的第一行寫成b，我們將新的矩陣寫成\(\theta\)，我們可以得到：

\[XW+b = X_{new}\theta=\left[ \begin{matrix} 1 & x_{11} & x_{12} & \cdots & x_{1n}\\ 1 & x_{21} & x_{22} & \cdots & x_{2n}\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_{m1} & x_{m2} & \cdots & x_{mn} \end{matrix} \right] \cdot \left[ \begin{matrix} b & w_1 & w_2 \cdots w_n \end{matrix} \right]\]

之后，我們對均方差進行變形：

\[J(\theta) = \frac{1}{2m}\sum_{i=1}^m(x_i\cdot \theta - y_i)^2 \]

這里的m是樣本的數量，是一個常數，我們對均方差乘上這個系數並不會影響\(\theta\)的取值。這個\(J(\theta)\)就是我們常說的模型的損失函數，這里它代表的意義是所有樣本均方差的平均數的1/2。

這里的損失其實是誤差的意思，損失函數的值越小，說明模型的誤差越小，和真實結果越接近。

我們計算\(J(\theta)\)對\(\theta\)的導數：

\[\frac{\partial J(\theta)}{\partial \theta}=\frac{1}{m}X^T(X\theta - Y) \]

我們令導數等於0，由於m是常數，可以消掉，得到：

\[\begin{aligned} X^T(X\theta - Y) &= 0 \\ X^TX\theta &= X^TY \\ \theta &= (X^T \cdot X)^{-1}X^TY \end{aligned} \]

上面的推導過程初看可能覺得復雜，但實際上並不困難。只是一個簡單的求偏導的推導，我們就可以寫出最優的\(\theta\)的取值。

從這個公式來看並不難計算，實際上是否真的是這么簡單呢？我們試着用代碼來實驗一下。

代碼實驗

為了簡單起見，我們針對最簡單的場景：樣本只有一個特征，我們首先先試着自己生產一批數據：

import numpy as np
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

import matplotlib.pyplot as plt
plt.scatter(X, y)

我們先生成一百個0~2范圍內的x，然后y= 3x + 4，為了模擬真實存在誤差的場景，我們再對y加上一個0~1范圍內的誤差。

我們把上面的點通過plt畫出來可以得到這樣一張圖：

接着，我們就可以通過上面的公式直接計算出\(\theta\)的取值了：

def get_theta(x, y):
    m = len(y)
    # x中新增一列常數1
    x = np.c_[np.ones(m).T, x]
    # 通過公式計算theta
    theta = np.dot(np.dot(np.linalg.inv(np.dot(x.T, x)), x.T), y)
    return theta

我們傳入x和y得到\(\theta\)，打印出來看，會發現和我們設置的非常接近：

最后，我們把模型擬合的結果和真實樣本的分布都畫在一張圖上：

# 我們畫出模型x在0到2區間內的值
X_new = np.array([[0],[2]])
# 新增一列常數1的結果
X_new_b = np.c_[np.ones((2, 1)), X_new]
# 預測的端點值
y_predict = X_new_b.dot(theta)

# 畫出模型擬合的結果
plt.plot(X_new, y_predict,"r-")
# 畫出原來的樣本
plt.scatter(X,y)
plt.show()

得到的結果如下：

從結果上來看模型的效果非常棒，和我們的逾期非常吻合，並且實現的代碼實在是太簡單了，只有短短幾行。

但實際上，有一點我必須要澄清，雖然上面的代碼只有幾行，非常方便，但是在實際使用線性回歸的場景當中，我們並不會直接通過公式來計算\(\theta\)，而是會使用其他的方法迭代來優化。

那么，我們為什么不直接計算，而要繞一圈用其他方法呢？

原因也很簡單，第一個原因是我們計算\(\theta\)的公式當中用到了逆矩陣的操作。在之前線性代數的文章當中我們曾經說過，只有滿秩矩陣才有逆矩陣。如果\(X^T \cdot X\)是奇異矩陣，那么它是沒有逆矩陣的，自然這個公式也用不了了

當然這個問題並不是不能解決的，\(X^T \cdot X\)是奇異矩陣的條件是矩陣\(X\)當中存在一行或者一列全為0。我們通過特征預處理，是可以避免這樣的事情發生的。所以樣本無法計算逆矩陣只是原因之一，並不是最關鍵的問題。

最關鍵的問題是復雜度，雖然我們看起來上面核心的代碼只有一行，但實際上由於我們用到了逆矩陣的計算，它背后的開銷非常大，\(X^T \cdot X\)的結果是一個n * n的矩陣，這里的n是特征的維度。這樣一個矩陣計算逆矩陣的復雜度大概在\(n^{2.4}\)到\(n^3\)之間。隨着我們使用特征的增加，整個過程的耗時以指數級增長，並且很多時候我們的樣本數量也非常大，我們計算矩陣乘法也會有很大的開銷。

正是因為以上這些原因，所以通常我們並不會使用直接通過公式計算的方法來求模型的參數。

那么問題來了，如果我們不通過公式直接計算，還有其他方法求解嗎？

答案當然是有的，由於篇幅限制，我們放到下一篇文章當中。

今天的文章就到這里，掃描下方二維碼，關注我的微信公眾號，獲取更多文章，你們的支持是我最大的動力。