梯度降低法不是一個機器學習算法,而是一種基於搜索的最優化方法,用於最小化一個效用函數。python
假設存在一個只有一個參數 $\theta$ 的損失函數 $J$,想找到最小極值處的 $\theta$,如圖所示:git
藉助於損失函數 $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
做圖表示以下:學習
定義函數 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$ 的取值變化,能夠用圖片表示,以下(紅色的點):
對於學習率 $\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$ 取值變化如圖:
最後,把梯度降低法封裝到方法中:
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),此時搜索路徑如圖所示:
隨機梯度降低法不能保證梯度降低的方向就是損失函數減少最快的方向(有時會是增大的方向),但總體上是會到達最小值附近的(知足必定的精度)。
同時在隨機梯度降低法中學習率 $\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