線性回歸中的梯度下降法（實現以及向量化並進行數據歸一化）
多元線性回歸中的梯度下降法
我們試一下應用在多元線性回歸中，對於線性回歸的問題，前面的基本都是使每一次模型預測出的結果和數據所對應的真值的差的平方的和為損失函數，對於參數來說，有n+1個元素，這種情況下，我們就需要變換式子
這實際上就是求對應的梯度值，梯度本身也是一個向量，含有n+1個元素，即對每一個參數進行一次偏導數，這樣一來，梯度就代表了方向，對應J增大最快的方向，區別在於這里變成了向量操作
原理終了，這樣我們的目標就變成了使下面的式子盡可能的小
梯度就是J對每一個維度上的參數進行求導，導數怎么求呢，我們除了對截距求導的時候是要在其基礎上乘上-1以外
其余的n個參數都是按照鏈式求導來計算
我們將其排成一個列向量，最后整理就可以得到
很明顯的是，我們不難發現其中每一項都是若干數量的樣本進行計算以后的求和結果，這就說明了樣本數量是影響着梯度的的，也就是說樣本數量越大，求出來的梯度中，元素相應的也就越大，這是不太合理的，應該和m無關才比較好，這樣我們再為其除以m，這樣我們整體的式子就變成了
這很明顯，J實際上就是MSE的值
有的時候會對J取不一樣的公式
實際上出來的式子的結果差別不大，但是如果我們沒有1/m的話，我們數據過大就會導致一些問題，這就說明了，當我們使用梯度下降法來進行求函數最小值，有時候我們需要對目標函數進行特殊的設計，雖然理論上沒問題，但是實際上數據過大可能影響效率，這種一般被稱為批量梯度下降法
實現：
（在notebook中）
首先使用我們自己設計的虛假的例子來實現梯度下降法的方式
import numpy as np
import matplotlib.pyplot as plt
生成一組一維數組x，size為100，對於y，使用線性的方式生成，在系數截距的上加上一個噪音
x = 2 * np.random.random(size=100)
y = x * 3. + 4. + np.random.normal(size=100)
為了代碼容易的擴展到多維，用X來代替x，設置成一百行一列這個結構，然后看一下shape
X = x.reshape(-1,1)
X.shape
y.shape
輸出如下
然后簡單的繪制出圖像
plt.scatter(x,y)
圖像如下
為了保證隨機性的可重復性，可以在隨機生成前用np.random.seed（）的方式設置一個種子，我這里就沒有設置出來
下面使用梯度下降法訓練
首先我們先設置一個計算損失函數J的值，我們要存進去theta，X_b這個矩陣以及y，相應的直接返回真值與預測值的差值的平方再除以樣本數這個結果，設置一個異常檢測來防止數據過大導致出現問題
def J(theta, X_b, y):
try:
return np.sum((y - X_b.dot(theta))**2) / len(X_b)
except:
return float('inf')
然后我們來設置dJ來作為計算其導數的函數，參數與J一樣，求導方面，設置最后導數為res，其空間為theta的量，然后我們對第一個元素進行特殊處理，對X_b.dot(theta) - y在進行求和，剩下的n項使用循環就可以得出，對X_b.dot(theta) - y和X_b[:,i]每一個樣本取出第i個特征的這兩個向量進行點乘，然后返回乘上2/m的結果
def dJ(theta,X_b,y):
res = np.empty(len(theta))
res[0] = np.sum(X_b.dot(theta) - y)
for i in range(1,len(theta)):
res[i] = (X_b.dot(theta) - y).dot(X_b[:,i])
return res * 2 / len(X_b)
盡管例子是簡單的線性回歸，但是這個方法對於多元來說也是通用的，然后我們就可以開始梯度下降的過程，我們需要傳進X_b和y，相應的計算也是要更改一下，其他的和之前一樣不變，相應的，將之前設置用來跟蹤theta的刪除掉即可，最后返回theta值
def gradient_descent(X_b, y, initial_theta,eta,n_iters = 1e4,epsilon=1e-8):
theta = initial_theta
i_iter = 0
while i_iter < n_iters:
gradient = dJ(theta,X_b,y)
last_theta = theta
theta = theta - eta * gradient
if(abs(J(theta,X_b,y) - J(last_theta,X_b,y)) < epsilon):
break
i_iter += 1
return theta
下面我們就可以開始調用了，我們先為原先的X添加一列，每一個元素都為1，相應的y值不變，我們需要一個初始化theta，將initial_theta設置成元素數量為特征數的數量加1的向量，元素數量即為X_b的列數，對eta取0.01
X_b = np.hstack([np.ones((len(X),1)),X])
initial_theta = np.zeros(X_b.shape[1])
eta = 0.01
然后我們進行調用
theta = gradient_descent(X_b, y, initial_theta, eta)
查看一下theta值是多少
結果如下，其對應的是截距和斜率，雖然有些差距，但是可以發現還是成功地訓練了
我們在LinearRegression中新加一個函數fit_gd，將上面的代碼封裝進取，包括了計算損失函數的函數，計算梯度的函數以及梯度下降法的函數
def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
"""根據訓練數據集使用梯度算法訓練模型"""
assert X_train.shape[0] == y_train.shape[0], \
"the size of X_train must be equal to the size of y_train"
def J(theta, X_b, y):
try:
return np.sum((y - X_b.dot(theta)) ** 2) / len(y)
except:
return float('inf')
def dJ(theta, X_b, y):
res = np.empty(len(theta))
res[0] = np.sum(X_b.dot(theta) - y)
for i in range(1, len(theta)):
res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
return res * 2 / len(X_b)
def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
theta = initial_theta
cur_iter = 0
while cur_iter < n_iters:
gradient = dJ(theta, X_b, y)
last_theta = theta
theta = theta - eta * gradient
if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
break
cur_iter += 1
return theta
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
initial_theta = np.zeros(X_b.shape[1])
self._theta = gradient_descent(X_b, y_train, initial_theta, eta, n_iters)
self.interception_ = self._theta[0]
self.coef_ = self._theta[1:]
return self
然后我們試一下使用我們自己的代碼
from LinearRegression import LinearRegression
實例化對象並調用一下梯度下降法進行訓練
lin_reg = LinearRegression()
lin_reg.fit_gd(X,y)
使用lin_reg.coef_可以看到此時系數結果
通過lin_reg.intercept_可以看到此時截距結果
可以發現和之前是一樣的
梯度下降法的向量化和數據標准化
向量化的處理主要集中在梯度的過程中，我們之前是一項一項的將對應的元素給求出來，那么我們嘗試將其轉換成矩陣的運算，這樣算起來就會簡單很多，我們對第一項乘上一個恆等於1的X0，這樣J就會變成
那么我們開始進行向量化，我們將這看作為兩個向量進行點乘的形式，我們可以對梯度里面的式子理解成一個矩陣的運算，我們將前一個向量給單獨提出來並將m項給拆開
我們可以將后面的向量寫成這個樣式，其就是X_b這個矩陣
這就相當於是我們對這個矩陣中每一個對應的列進行點乘，其實際上就是原來式子中的計算方式，即對第0列進行點乘，得到的就是以前的式子的第0項，以此類推，知道點乘到第n列，這樣我們就有了n+1項，這樣寫其實就是將前一個向量寫作成一個1m的向量，矩陣則是mn+1的矩陣
將求之前的梯度過程轉換成了兩個矩陣的乘法運算，最終就可以得到一個1*n+1的向量，而我們原來要求的梯度中也有n+1個元素，這兩個中的元素是相等一一對應的，這樣我們就得到了新的式子（由於之后使用的時候都是列向量，但是我們上面的分析中其是一個行向量，所以要進行一下轉置）
雖然在numpy中的計算是不區分行列向量的，但是為了直觀以及好看，我們將其也進行一下轉換，對最后的結果整體進行一下轉置，很好理解，不再贅述，新的式子即為
這個列向量就是我們要求的梯度
這里我們只要將對梯度的求解部分dJ進行一下改變，將其變成
def dJ(theta, X_b, y):
return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)
那么我們使用真實的數據進行訓練（老波士頓房價了）
（notebook中）
老樣子引用加限制范圍再進行分割，分割為測試數據集以及訓練數據集，隨機種子為666
import numpy as np
from sklearn import datasets
boston = datasets.load_boston()
X = boston.data
y = boston.target
X = X[y < 50.0]
y = y[y < 50.0]
from model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,seed=666)
首先先用正規方程式來訓練一個線性回歸模型，同時打印出准確度
from LinearRegression import LinearRegression
lin_reg1 = LinearRegression()
%time lin_reg1.fit_normal(X_train,y_train)
lin_reg1.score(X_test,y_test)
結果如下
然后我們使用梯度下降法，首先先進行實例化，調用fit_gd並將訓練集傳進來
lin_reg2 = LinearRegression()
lin_reg2.fit_gd(X_train,y_train)
可以發現，又出現了一些報錯
我們使用lin_reg2.coef_來看一下系數，可以發現結果都是nan
原因是我們將其放在了一個真實的數據集中，並沒有觀察這個數據的情況就進行了操作
我們可以使用X_train[:10,:]來看看前十行的數據
我們可以發現，在其中每一個特征的數據規模是不一樣的，面對這樣的一個數據，我們最后求出來的結果很有可能也是非常大的，我們使用的默認的eta為0.01，其出來的結果還是發散的，因此我們降低eta，設置eta為0.000001
lin_reg2.fit_gd(X_train,y_train,eta=0.000001)
輸出結果可以看出來，此時是沒有報錯的
我們可以看看這個情況下的值是多少
lin_reg2.score(X_test,y_test)
結果為
通過這個值我們可以發現，我們還沒有達到這個的最小值，很有可能是因為現在的eta太小了，導致行進的步長很小，那么我們使用更多的循環次數可能就可以了
我們設置成循環次數為一百萬次
%time lin_reg2.fit_gd(X_train,y_train,eta=0.000001,n_iters=1e6)
結果為
此時我們用lin_reg2.score(X_test,y_test)來看一下這個最后的值是多少
很顯然還是沒有達到對應的最小值，雖然我們可以使用更多的循環來使其達到最小值，但是很顯然，這太耗時了，因此我們在使用梯度下降前，使用數據歸一化就可以很好地解決這個情況
數據歸一化
先加載上StandardScaler這個類，然后就是實例化加fit操作就可以了
from sklearn.preprocessing import StandardScaler
standardScaler = StandardScaler()
standardScaler.fit(X_train)
可以得到
然后我們就可以得到標准化以后的結果X_train_standard
X_train_standard = standardScaler.transform(X_train)
這樣我們就對原始數據進行了歸一化處理，同樣的進行fit，使用默認值並進行計時
lin_reg3 = LinearRegression()
%time lin_reg3.fit_gd(X_train_standard,y_train)
得出的輸出結果
可以看出來這是非常快就訓練完成了，由於我們之前已經進行了數據歸一化，同樣的也要對test進行數據歸一化
X_test_standard = standardScaler.transform(X_test)
lin_reg3.score(X_test_standard,y_test)
得到的輸出結果為
這是我們已經找到了這個損失函數的最小值，同時這個速度是非常快的
我們通過設計一個虛擬的測試用例來說明梯度下降法的優勢在哪里
我們設置樣本數為1000個，相應的特征值為5000，為什么樣本數會少於特征值呢，這是因為我們在計算梯度的時候，其實是需要所有的樣本都來參與計算的，當樣本量比較大的時候，我們計算也是會比較慢的，后面可以進行隨機梯度來改進，不過這里就不變了
我們設置一個使用隨機化的正態分布的基於這個規模的量，這里由於之前歸一化過，這里就不需要歸一化了，在設置一個初始化的隨機生成的theta，然后我們就可以生成一個y了，將矩陣點乘上系數再加上true_theta第0位元素再加上一個噪音即可得到
m =1000
n = 5000
big_X = np.random.normal(size=(m,n))
true_theta = np.random.uniform(0.0,100.0,size=n+1)
big_y = big_X.dot(true_theta[1:]) + true_theta[0] + np.random.normal(0.,10.,size=m)
基於上述的數據，我們先實例化后在使用正規方程來解析
big_reg1 = LinearRegression()
%time big_reg1.fit_normal(big_X,big_y)
結果如下
然后我們再實例化一個對象，使用梯度下降法來解析
big_reg2 = LinearRegression()
%time big_reg2.fit_gd(big_X,big_y)
結果為
可以看出來，我們使用梯度下降法的方式來訓練是比正規方程的耗時要短的，加大數據集的話，優勢會更加的明顯
