機器學習之梯度降低法與線性迴歸

梯度降低法不是一個機器學習算法,而是一種基於搜索的最優化方法,用於最小化一個效用函數。python

簡單理解梯度降低法

假設存在一個只有一個參數 $\theta$ 的損失函數 $J$,想找到最小極值處的 $\theta$,如圖所示:git

clipboard.png

藉助於損失函數 $J$ 在 $\theta$ 處的切線,能夠直觀的反映出損失函數 $J$ 在 $\theta$ 處的導數大小;導數的大小表明着 $\theta$ 變化時 $J$ 相應的變化。github

同時導數也能夠表明 $J$ 增大的方向,若是將導數與 $-\eta$ 相乘,即 $-\eta\frac{dJ}{d\theta}$ 表明了 $J$ 減少的方向。算法

$\eta$ 表示學習率,是梯度降低法的一個超參數,其取值影響最優解的速度。過小會減慢收斂學習速度,太大可能致使不收斂。

若是 $J$ 中存在多處導數爲零的狀況,即存在一個全局最優解和多個局部最優解,此時能夠屢次使用梯度降低法,每次隨機化一個初始點。dom

對於有多個 $\theta$ 的 $J$ 相似,即找出全局最優解處的這些 $\theta$ 的值。機器學習

模擬梯度降低法

首先,模擬一個損失曲線 $J$:函數

import numpy as np

plot_x = np.linspace(-1, 6, 141)
plot_y = (plot_x - 2.5) ** 2 - 1

做圖表示以下:學習

clipboard.png

定義函數 dJ() 用於求 $J$ 在 $\theta$ 處的導數:優化

def dJ(theta):
    return 2 * (theta - 2.5)

函數 J() 用於求 $J$ 在 $\theta$ 處的大小:spa

def J(theta):
    return (theta - 2.5) ** 2 - 1

接着使用梯度降低法,首先給 $\theta$ 賦一個初值 0 及學習率 $\eta$ 爲 0.1,接着在循環裏進行屢次梯度降低。每次循環都要求得 $J$ 在 $\theta$ 處的導數值 $gradient$,而且 $\theta$ 嚮導數的負方向移動,即:$\theta=\theta-\eta*gradient$。

因爲計算機計算浮點數存在偏差,對於求得的 $\theta$ 可能不能恰好等於 0,所以設定一個精度值(epsilon = 1e-8),若是新的 $\theta$ 對應的損失函數的值與上一次 $\theta$ 對應的損失函數的值的差值知足精度要求,就表示找到了要找的 $\theta$。程序以下:

theta = 0.0
eta = 0.1
epsilon = 1e-8

while True:
    gradient = dJ(theta)
    last_theta = theta
    
    theta = theta - eta * gradient
    if (abs(J(theta) - J(last_theta)) < epsilon):
        break

運行程序求得的 $\theta$ 爲:2.499891109642585。

對於 $\theta$ 的取值變化,能夠用圖片表示,以下(紅色的點):

clipboard.png

對於學習率 $\eta$,這裏取值 0.1 是沒有問題的,但若是取值 1.1 程序運行就會報錯:

OverflowError                             Traceback (most recent call last)
<ipython-input-6-5bd2217401d5> in <module>
      8 
      9     theta = theta - eta * gradient
---> 10     if (abs(J(theta) - J(last_theta)) < epsilon):
     11         break
     12 

<ipython-input-5-bd41aa589cfc> in J(theta)
      1 def J(theta):
----> 2     return (theta - 2.5) ** 2 - 1

OverflowError: (34, 'Result too large')

這是由於學習率過大會致使 J(theta) 愈來愈大。爲了使程序不會報錯,修改 J() 方法:

def J(theta):
    try:
        return (theta - 2.5) ** 2 - 1
    except:
        return float('inf')

注意,當無窮減無窮時,結果時 nan 而不是 0,此時 if (abs(J(theta) - J(last_theta)) < epsilon) 將永遠沒法觸發而使得程序進入死循環。爲了解決這個問題,增長一個新的超參數 n_iters,表示可以執行循環的最大次數。

好比使 n_iters=10,$\theta$ 取值變化如圖:

clipboard.png

最後,把梯度降低法封裝到方法中:

def gradient_descent(initial_theta, eta, n_iters=1e4, epsilon=1e-8):
    theta = initial_theta
    i_ters = 0

    while i_ters < n_iters:
        gradient = dJ(theta)
        last_theta = theta

        theta = theta - eta * gradient
        if (abs(J(theta) - J(last_theta)) < epsilon):
            break
        
        i_ters += 1
    
    return theta

多元線性迴歸中的梯度降低法

原理

多元線性迴歸的損失函數爲:

$$ J=\sum_{i=1}^{m}(y^{(i)} - \hat{y}^{(i)})^2 $$

其中:$\hat{y}^{(i)} = \theta_{0} + \theta_{1}X_{1}^{(i)} + \theta_{2}X_{2}^{(i)} + ... + \theta_{n}X_{n}^{(i)}$ 。

對 $J$ 求導爲:

$$ \nabla J=(\frac{\partial J}{\partial \theta_0},\frac{\partial J}{\partial \theta_1},...,\frac{\partial J}{\partial \theta_n}) $$

其中:$\frac{\partial J}{\partial \theta_i}$ 爲偏導數,與導數的求法同樣。

對 $\nabla J$ 進一步計算:

$$ \nabla J(\theta) = \begin{pmatrix} \frac{\partial J}{\partial \theta_0} \\\ \frac{\partial J}{\partial \theta_1} \\\ \frac{\partial J}{\partial \theta_2} \\\ \cdots \\\ \frac{\partial J}{\partial \theta_n} \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)}\theta)·(-1) \\\ \sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)}\theta)·(-X_1^{(i)}) \\\ \sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)}\theta)·(-X_2^{(i)}) \\\ \cdots \\\ \sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)}\theta)·(-X_n^{(i)}) \end{pmatrix} = 2·\begin{pmatrix} \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)}) \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_1^{(i)} \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_2^{(i)} \\\ \cdots \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_n^{(i)} \end{pmatrix} $$

其中:$X_b = \begin{pmatrix} 1 & X_1^{(1)} & X_2^{(1)} & \cdots & X_n^{(1)} \\\ 1 & X_1^{(2)} & X_2^{(2)} & \cdots & X_n^{(2)} \\\ \cdots & & & & \cdots \\\ 1 & X_1^{(m)} & X_2^{(m)} & \cdots & X_n^{(m)} \end{pmatrix}$

這個結果是與樣本數量 m 相關的,爲了使結果與 m 無關,對這個梯度除以 m,即:

$$ \nabla J(\theta) = \frac{2}{m}·\begin{pmatrix} \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)}) \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_1^{(i)} \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_2^{(i)} \\\ \cdots \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_n^{(i)} \end{pmatrix} $$

此時,目標函數就成了使 $\frac{1}{m}\sum_{i=1}^{m}(y^{(i)} - \hat{y}^{(i)})^2$ 儘量小,即均方偏差儘量小:

$$ J(\theta) = MSE(y, \hat{y}) $$

使用梯度降低法訓練模型

首先模擬訓練數據集:

import numpy as np

x = np.random.random(size=100)
y = x * 3.0 + 4.0 + np.random.normal(size=100)

X = x.reshape(-1, 1)

定義函數 J() 計算損失函數的值:

def J(theta, X_b, y):
    try:
        return np.sum((y - X_b.dot(theta)) ** 2) / len(X_b)
    except:
        return float('inf')

函數 dJ() 對 $\theta$ 求導數:

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)

注意:對 $J$ 求導更好的方式是進行向量化處理,即 $\nabla J(\theta) = \frac{2}{m}·X_b^T·(X_b\theta-y)$,dJ() 改寫爲:

def dJ(theta, X_b, y):
    return X_b.T.dot(X_b.dot(theta) - y) * 2 / len(X_b)

梯度降低的過程爲:

def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
    theta = initial_theta
    i_ters = 0

    while i_ters < 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_ters += 1
        
    return theta

執行程序找出最優的 $\theta$ 以下:

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$ 結果爲:

array([4.0269033, 3.0043078])

將梯度降低法封裝到線性迴歸算法的模型訓練方法 fit_gd() 中:

class LinearRegression:
    # other codes here
    
    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(X_b)
            except:
                return float('inf')

        def dJ(theta, X_b, y):
            return X_b.T.dot(X_b.dot(theta) - y) * 2 /len(X_b)

        def gradient_descent(X_b, y, initial_theta, eta, n_iters=n_iters, epsilon=1e-8):
            theta = initial_theta
            i_ters = 0

            while i_ters < 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_ters += 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)
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self
注意:

在真實的數據集 (X, y) 中,X 總體不在一個規模上會影響梯度的結果,而梯度的結果再乘以 $\eta$ 獲得步長就太大或過小,從而致使訓練出的模型可能不好。所以在使用梯度降低法以前,最好進行數據歸一化。

隨機梯度降低法

$\nabla J(\theta) = \frac{2}{m}·\begin{pmatrix} \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)}) \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_1^{(i)} \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_2^{(i)} \\\ \cdots \\\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_n^{(i)} \end{pmatrix}$ 中每一項都要對全部樣本進行計算,所以這種梯度降低法稱爲批量梯度降低法(Batch Gradient Descent)。若是 m 很是大,使用批量梯度降低法計算梯度的計算量就會很是大。

改進的方法是對每一項計算時再也不計算全部的樣本而是選取其中一個樣本進行計算,即:

$$ 2·\begin{pmatrix} (X_b^{(i)}\theta - y^{(i)})·X_0^{(i)} \\\ (X_b^{(i)}\theta - y^{(i)})·X_1^{(i)} \\\ (X_b^{(i)}\theta - y^{(i)})·X_2^{(i)} \\\ \cdots \\\ (X_b^{(i)}\theta - y^{(i)})·X_n^{(i)} \end{pmatrix} = 2·(X_b^{(i)})^T·(X_b^{(i)}\theta-y^{(i)}) $$

這樣的方式就是隨機梯度降低法(Stochastic Gradient Descent),此時搜索路徑如圖所示:

clipboard.png

隨機梯度降低法不能保證梯度降低的方向就是損失函數減少最快的方向(有時會是增大的方向),但總體上是會到達最小值附近的(知足必定的精度)。

同時在隨機梯度降低法中學習率 $\eta$ 的取值是逐漸遞減的,爲了防止固定取值的學習率使得梯度降低到達最優解附近時繼續跳出這個範圍。一個合理的取值方式爲:

$$ \eta = \frac{t0}{i\_iter + t1} $$

其中:$t0、t1$ 爲超參數。

定義函數 dJ_sgd() 對應批量梯度降低法中對損失函數求導的過程,此時傳入函數的就再也不是全部樣本 $(X_b, y)$ 了,而是其中一個樣本 $(X_b^{(i)}, y^{(i)})$:

def dJ_sgd(theta, X_b_i, y_i):
    return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2.

函數 sgd() 中定義了一個 learning_rate() 方法用來計算學習率,傳入參數爲當前迭代次數。

由於要隨機的選取樣本中的一個,但又要將全部樣本看一遍,因此咱們將這個樣本集打亂造成一個新的樣本集 $(X_{b,new}^{(i)}, y_{new}^{(i)})$,同時指定參數 n_iters 表示將這個樣本集看幾遍:

def sgd(X_b, y, initial_theta, n_iters, t0, t1):

    def learning_rate(t):
        return t0 / (t + t1)

    theta = initial_theta
    m = len(X_b)

    for cur_iter in range(n_iters):
        indexes = np.random.permutation(m)
        X_b_new = X_b[indexes]
        y_new = y[indexes]

        for i in range(m):
            grandient = dJ_sgd(theta, X_b_new[i], y_new[i])
            theta = theta - learning_rate(cur_iter * m + i) * grandient

            return theta

將隨機梯度降低法封裝到線性迴歸算法的模型訓練方法 fit_sgd() 中:

class LinearRegression:
    # other codes here
    
    def fit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50):
        def dJ_sgd(theta, X_b_i, y_i):
            return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2.

        def sgd(X_b, y, initial_theta, n_iters, t0, t1):
            
            def learning_rate(t):
                return t0 / (t + t1)
            
            theta = initial_theta
            m = len(X_b)

            for cur_iter in range(n_iters):
                indexes = np.random.permutation(m)
                X_b_new = X_b[indexes]
                y_new = y[indexes]

                for i in range(m):
                    grandient = dJ_sgd(theta, X_b_new[i], y_new[i])
                    theta = theta - learning_rate(cur_iter * m + i) * grandient
                
            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])

        self._theta = sgd(X_b, y_train, initial_theta, n_iters, t0, t1)
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

在 Scikit Learn 的 linear_model 模塊中提供了一個使用隨機梯度降低法的迴歸算法 SGDRegressor:

from sklearn.linear_model import SGDRegressor

源碼地址

Github | ML-Algorithms-Action

相關文章
相關標籤/搜索