背景

上一篇文章用最小二乘法（即公式法）求出了線性迴歸的引數 theta ；本篇程式碼介紹用梯度下降法求極小值。

原理

實在不知道怎麼描述啊，okay，從 山頂走向山腳有 n 條路，問題來了：

• 捷徑？最快的那條路徑。
• 以多大的步伐走比較合適？比較走的太快，容易走過頭。
  用計算機模擬實現，我們要容易想到 x = x - deta

求解過程模擬

已知：y = (x-2.5)^2 -1 ，求 y 的最小值。
當然極小值（這裡也是最小值）一眼就能看出是 -1 ，但對於 三次曲線等更復雜的函式，怎麼求極小值呢？這裡我們模擬梯度下降法搜尋拋物線極小值的過程：

import numpy as np
import matplotlib.pyplot as plt

plot_x = np.linspace(-1.,6.,141) #  資料準備
plot_y = (plot_x-2.5)**2 - 1 # 生成一個目標二次函式

 

def J(theta):  # 即待求最小值的二次函式；模擬損失函式
    try:
        return (theta-2.5)**2 -1
    except:
        return float('inf')
    
def dJ(theta): # 二次函式的導數
    return 2*(theta-2.5)

theta_history = []
def gradient_descent(initial_theta, eta, n_iters = 1e4, epsilon=1e-8):
    theta = initial_theta
    i_iter = 0
    theta_history.append(initial_theta)

    while i_iter < n_iters: # 保證方法能退出？eta值過大，可能無法收斂
        gradient = dJ(theta)
        last_theta = theta
        theta = theta - eta * gradient
        theta_history.append(theta)
        if(abs(J(theta) - J(last_theta)) < epsilon):  
            break           
        i_iter += 1
        
    return

def plot_theta_history():
    plt.plot(plot_x, J(plot_x))
    plt.plot(np.array(theta_history), J(np.array(theta_history)), color="r", marker='+')
    plt.show()
 

測試結果：

eta = 0.8 
theta_history = []
gradient_descent(0,eta)
plot_theta_history()

eta = 0.3
theta_history = []
gradient_descent(0,eta)
plot_theta_history()

eta = 1.1
theta_history = []
gradient_descent(0, eta, n_iters=10)
plot_theta_history()

輸出：

• 紅色折線點描述 theta 的取值，對應min(J)的搜尋路徑
• eta ：步長，對應下降速度
• n_iters ：保證 gradient_descent 方法可以退出？eta 過大可能導致演算法不收斂，一直執行
• epsilon ：計算機取值浮點誤差

線上性迴歸中的應用

先推導一下求偏導數的公式：

import numpy as np
class LinearRegression:
    def __init__(self):
        """初始化Linear Regression模型"""
        self.coef_ = None
        self.intercept_ = None
        self._theta = None
        
    def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
        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):
            return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(y)

        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.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    def predict(self, X_predict):
        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        return X_b.dot(self._theta)

    def score(self, X_test, y_test):
        """根據測試資料集 X_test 和 y_test 確定當前模型的準確度"""
        y_predict = self.predict(X_test)
        return r2_score(y_test, y_predict)

import numpy as np
from math import sqrt

def mean_squared_error(y_true, y_predict):
    """計算y_true和y_predict之間的MSE"""
    assert len(y_true) == len(y_predict), \
        "the size of y_true must be equal to the size of y_predict"

    return np.sum((y_true - y_predict)**2) / len(y_true)

def r2_score(y_true, y_predict):
    """計算y_true和y_predict之間的R Square"""
    return 1 - mean_squared_error(y_true, y_predict)/np.var(y_true)

測 試：

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

boston = datasets.load_boston()
X = boston.data
y = boston.target

X = X[y < 50.0]
y = y[y < 50.0]

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)

standardScaler = StandardScaler()
standardScaler.fit(X_train)
X_train_standard = standardScaler.transform(X_train)
X_test_standard = standardScaler.transform(X_test)

lin_reg = LinearRegression()
lin_reg.fit_gd(X_train_standard, y_train)
lin_reg.score(X_test_standard, y_test)  # 輸出：0.80087954605863654

注： 使用梯度下降法求最小值，一般要將資料歸一化。

批量梯度下降 VS 隨機梯度下降法

粗糙的理解一下，Batch Gradient Descent (bgd) 正如上面所做的那樣，每走一步，都需要考慮所有的樣本;而Stochastic Gradient Descent (sgd) 每走一步，只考慮其中一個樣本，So ：

1. SGD 速度比 BGD快（迭代次數少）
2. BGD一定能得到一個區域性最優解（線上性迴歸模型中一定是得到一個全域性最優解），SGD由於隨機性的存在,可能導致最終結果比BGD的差
3. SGD有可能跳出某些小的區域性最優解，所以不會比BGD差
4. BGD 能保證每次 J 的值都向最快的方向減少，SGD 可能存在使 J 值增加的步法
   注：SGD為了保證收斂，步驟 eta 值需要漸漸減少
   
   其中 i_iters 為迴圈的次數，為了保證第1次 與第100迴圈時，eta值相差不那麼大，額外增加兩個平衡引數 a與b,經驗值：a = 5 ; b = 50

總 結

梯度下降法是一種基於搜尋的最優化方法（不是一個機器學習演算法），主要用於最小化一個損失函式；
使用梯度下降法前，為了提高收斂（求出極小值）的速度，需要對資料進行歸一化；
梯度上升法，使用者最大化一個效用函式