算法學習筆記——最小二乘法的回歸方程求解

最小二乘法的回歸方程求解

最近短暫告別大數據，開始進入到了算法學習的領域，這時才真的意識到學海無涯啊，數學領域充滿了無限的魅力和樂趣，可以說更甚於計算機帶給本人的樂趣，由於最近正好看到線性代數，因此，今天我們就來好好整理一下機器學習領域中的一個非常重要的算法——最小二乘法，那么，廢話不多說，我們直接開始吧 !

1. 最小二乘法介紹

1.1 舉例

現實生活中，我們經常會觀察到這樣一類現象，比如說某個男的，情商很高，身高180，家里很有錢，有房，有車，是個現充，結果就是他有好幾個女朋友，那么從一個觀測者的角度來看，該男性具備好多個特征(比如EQ值較高，身高較高，有錢對應的布爾值是True等等)，輸出結果就是女友的個數；這只是一條記錄，那么，當我們將觀測的樣本數擴大到很多個時，每個個體作為輸入，而輸出就是每個個體的女朋友數量；於是在冥冥之中，我們就能感覺到一個男性擁有的女友數量應該和上述特征之間存在着某種必然的聯系。然后可以這樣理解，決定一個男性可以交到女友數量的因素有很多，那么，在那么多的因素之中，肯定有幾項因素比較重要，有幾項相對不那么重要，我們暫時將每個因素的重要程度用一個數值來表示，可以近似理解為權重，然后將每個權重和因素的數值相乘相加，最后再加上一個常數項，那么這個式子就可以理解為一個回歸方程。

1.2 SSE,SST和SSR

有了上述的基礎，我們就可以做這樣一件事，預先設定好一個方程(先簡單一點，假設該方程只有一個自變量)：y = ax + b，a和b是我們要求出來的；那么，我們可不可以這樣理解，每輸入一個x，即能通過這個計算式輸出一個結果y，如果輸出的y和真實的y偏差是最小的，那么不就能說明這個方程擬合的是最佳的了嗎？順着這個思路，原問題就可以演變成一個求解當a和b各為多少時能使得這個偏差值最小的求最優化問題了，或者說我們的目標就是求使得SSE最小的a和b的值。於是乎，SSE這個概念由此得來：SSE的英文全稱為Sum Squared Error，誤差平方和，它在數值上等於預測值減去真實值的平方和，計算公式為：

 

既然都已經價紹SSE了，我們就把另外兩個成員也一並介紹了吧，先是SST，Sum Squared Total，總平方和，它在數值上等於真實值減去真實值的均值的平方和，注意不是方差！計算公式為：

再是SSR，Sum Squared Regression，回歸平方和，它在數值上等於預測值減去真實值均值的平方和，計算公式為：

三者在數值上存在這樣的關系：

1.3 R-Square及其他評估指標

經過上述公式的介紹，其實可以大致感覺到，SSE值越小則說明回歸方程擬合得越好，但是到底要多小，小到什么程度才算比較好呢。比如，我們計算出來SSE等於10000，足夠大了，但是我們並不能說這個方程擬合得不好，因為有可能，我們面對的數據量是100000000，相對這個數據量，10000的SSE絕對不能說不好。因此，我們在制定指標衡量模型擬合好壞時就必須將原始數據的因素考慮進去，於是乎，出現了我們常說的r平方，即R-Square指標，計算公式為：

有上述公式可知，SST是由原始數據得出的，因此可以看成是一個定值，SSE越小，模型擬合度越高，R-Square就越接近於1，反之，SSE越大，R-Square越接近於0，則說明模型解釋性越差，由於R-Sqare在數值上是介於0和1之間的，因此我們可以很明確的知道模型的擬合得到底是好還是不好。

當然，除了R-Square之外，還有其他好幾個指標也能夠評價一個回歸方程擬合的好壞，這里也簡單地介紹一下：

MSE(Mean Squared Error)：均方誤差

可以理解為就是SSE除了一個樣本數，計算公式為：

 

RMSE(Root Mean Squared Error)：均方根誤差

就是給均方誤差開了一個根號，計算公式為：

MAE(平均絕對誤差)：平均絕對誤差

計算的不再是平方和，而是絕對值之和，有點類似於曼哈頓距離和歐式距離之間的區別，計算公式為：

2. 二維平面內求解最小二乘法

2.1 假設條件及問題描述

設在一個直角坐標平面內，最優直線方程為：，為簡化后續書寫步驟，之后會一律表示為 ，將該直角坐標平面上的一系列的點，代入到SSE的計算公式：，現求解a和b的值，使得SSE能取到最小值。

2.2 計算過程詳解

2.2.1 SSE圖像形狀

由SSE的計算式可知，該式可以表示的一般形式為：，這樣的函數可以由基本函數進行一系列的拉伸平移等等的變換得到，研究該圖像的形狀，可知在三維空間中它是一個開口向上的拋物面，是一個凸函數，因此SSE必定可以取到最小值

2.2.2 求導，聯立方程組並求解

下一步就比較簡單了，對於這個方程，我們分別對a和b求偏導數，然后再聯立方程組將a和b求解出來即可

將SSE表示成關於a的表達式(即將b看作常量)，可得：

將SSE表示成關於b的表達式(即將a看作常量)，可得：

將SSE分別對a和b求偏導數，得：

求解該方程，最終得到：

注意到該式分子分母都有1/n²，將1/n² = 1/n * 1/n進行分配並表示成期望的形式，最終可將a的表達式變形為：

，同理b的表達式也可以這樣變形，最終求得：

以上，二維平面內的回歸方程的表達式求解步驟宣告完成！

3. 高維空間內求解最小二乘法

現在開始對原問題進行深入，現實情況錯綜復雜，如果需要用一個回歸方程來解釋某種現象，自變量的數量往往會有很多個，那這樣，為求解出回歸方程所需要用到的點組成的空間就必定是一個高維空間了，此時，需要使用到線性代數中的矩陣運算，現在開始求解這個方程的各項系數組成的矩陣

3.1 假設條件及問題描述

現將原問題擴展到p維空間，樣本數量為n，將這些樣本全部進行聯立，並設置截距項組成一個多元方程組：

用線性代數中矩陣的思想將上式表示成矩陣相乘的形式，可得：

我們令矩陣，矩陣，矩陣

原式就變成了  ，現在，我們開始求解SSE，可得：

注意：上述求解過程利用到了矩陣轉置公式 

最終問題變為：求解出矩陣A，使得！

3.2 計算過程詳解

我們首先將SSE展開，可得：

對該式進行適當簡單整理，可得：

觀察上式，注意到該式的第二以及第三項經過矩陣運算過后的結果均為一個數字，該數字正好等於 ，因此上式就可以進一步簡化為：

對SSE求導，並令它的導數為0，可得(此處用到了矩陣求導方面的知識點，本人也還未對其進行深入的研究，對此有興趣的朋友可以深入了解一下)：

以上，矩陣A的求解過程宣告完成！

3.3 注意事項

由於A式包含逆矩陣，因此在使用該算法之前必須先計算XTX的行列式的值並判斷其是否滿秩，若XTX中的行或列存在線性相關，就無法運用該算法，這時，就需要用到Ridge或是Lasso算法來優化該模型，該知識點會在之后的博客中介紹。

4. 逐步回歸算法

接下來介紹逐步回歸算法，之前我們講到的算法是根據公式推導直接得到最終結果的，而逐步回歸算法則是通過反復地迭代一步一步將結果優化，最終得到一個最優結果，其基本思路如下：

(1). 原方程組可以表示成：

(2). 現對到這 p + 1 個元素進行迭代，每次迭代時，保持其他元素不變，而只變化當前元素，變化方式為預先設置一個lambda值，將當前元素加上這個lambda以及減去這個lambda各計算一次SSE，取兩個SSE中的較小值，然后求得所有 p + 1 個元素的SSE的最小值，改變該元素的值即可，以上完成了一次迭代

 

(3). 將上述過程循環進行很多次，不斷地迭代之后能夠得到最優值

該算法具有一個一步步迭代的思想，由此可以引申到另一個機器學習領域中重要的算法：梯度下降法，這會在接下來的博客中進行總結。

5. 回歸算法的代碼實現

將上述算法使用python代碼實現寫成包如下：

import numpy as np
import pandas as pd

class regression:
    
    def __init__(self, data, intercept = True):
        self.X = np.mat(data.iloc[:,:-1].values)
        self.Y = np.mat(data.iloc[:,-1].values).T
        self.data = data
        self.intercept = intercept
        
    def basicfit(self):
        #This function aims to realize the basic linear regression algorithm
        if self.intercept == False:
            XTX = self.X.T * self.X
            if np.linalg.det(XTX) != 0:
                return np.ravel(XTX.I * self.X.T * self.Y)
            return
        ones = np.array([1] * self.data.shape[0])
        self.X = np.mat(np.column_stack((ones, self.X)))
        XTX = self.X.T * self.X
        if np.linalg.det(XTX) == 0:
            print("linear dependent columns exist!!!")
            return
        return np.ravel(XTX.I * self.X.T * self.Y)
    
    def cal_SSE(self, theta):
        #This function aims to calculate the SSE value from a given theta
        Y_pred = self.X * np.mat(theta).T
        return np.power(self.Y - Y_pred, 2).sum()
    
    def stepfit(self, lam, iternum):
        #This function aims to realize the step linear regression algorithm
        theta = np.zeros(self.X.shape[1])
        for _ in range(iternum):
            L = []
            for i in range(len(theta)):
                v_plus = theta[i] + lam
                v_minus = theta[i] - lam
                theta[i] = v_plus
                SSE_plus = self.cal_SSE(theta)
                theta[i] = v_minus
                SSE_minus = self.cal_SSE(theta)
                L_tmp = [(SSE_plus, v_plus, i), (SSE_minus, v_minus, i)]
                L_tmp.sort()
                L.append(L_tmp[0])
            L.sort()
            theta[L[0][2]] = L[0][1]
        return thetadef evaluate(self, coef):
        #This function aims to calculate SSE and R2 in order to evaluate whether the model is accurate or not
        Y_pred = self.X * np.mat(coef).T
        SSE = np.power(np.ravel(self.Y) - np.ravel(Y_pred), 2).sum()
        SST = np.power(np.ravel(self.Y) - np.ravel(self.Y).mean(), 2).sum()
        R2 = 1 - SSE / SST
        return SSE, R2