DL中的優化演算法與實現

04-18

這篇文章的內容主要是參考沐神的mxnet/gluon視頻中，Aston Zhang講解的優化部分。

文章可能有些長，小夥伴們可以耐心讀完，收穫一定會很大。。。

通過這麼長時間的學習，我們應該對於通過深度學習解決問題的大體流程有個宏觀的概念了吧？

設計模型
構造loss function
通過優化演算法，找到一個某種意義上的optim值

其實在找optim的值過程，就是我們通常所說的調參的過程。

同一個model，一樣的loss function，為什麼不同的人對於同一個數據集會得到不同的結果？這其實就是調參的魅力了。。。

當然，調參並不是盲目的瞎調，也不是靠著運氣亂調參，其背後也是遵循著一定的數學原理的。（本人是調參新手，如果說的不對，請各位評論區拍磚指教）

通過前面關於深度學習的介紹和沐神的視頻教程，我們已經接觸到了很多優化演算法。比如說，在訓練模型的時候，不斷迭代參數以最小化損失函數。

在大多數的ML問題中，我們往往都是這樣做的：

定義一個損失函數 loss function
通過優化演算法來最小化這個loss function
找到使得loss function最大或者最小的optim解作為model的parameter值

其實，這個優化過程中更多是針對訓練集上進行的，而實際的ML問題求解過程中，我們更需要的是對於測試集上的表現來衡量，比如前面講過的各種正則化，weight decay等手段來應對過擬合現象。

在這裡，主要關注loss function在訓練集上的得到的結果。又把loss function稱作目標函數。

現在，再來看看求解優化問題的challenge，我們知道，絕大數深度學習中的目標函數都很複雜。因此，很多優化問題並不存在解析解，所以，我們就需要通過基於數值方法的優化演算法來找到目標函數的近似解。而尋找近似解的過程，我們就需要不斷迭代更新解的數值，從而找到那個在某種意義上最optim的解。

從梯度下降說起

有關梯度下降的具體概念這裡就不詳細展開了。

如果不知道的小夥伴，可以去我的模式識別和機器學習筆記專欄中學習。

首先假設我們要學習的參數是 $x$
那麼， $x := x-eta frac{partial Loss}{partial x}$
這裡的 $x和Loss$ 都是一維實數集上的一個映射

那麼，我們需要回憶下，大學學過的泰勒展開。

因自百度百科中的概念：

數學中，泰勒公式是一個用函數在某點的信息描述其附近取值的公式。如果函數足夠平滑的話，在已知函數在某一點的各階導數值的情況之下，泰勒公式可以用這些導數值做係數構建一個多項式來近似函數在這一點的鄰域中的值。泰勒公式還給出了這個多項式和實際的函數值之間的偏差。

通過數學公式表達出來就是：

$f(x+epsilon) approx f(x)+f(x)epsilon+frac{f(x)}{2!}epsilon^2+frac{f(x)}{3!}epsilon^3+...R_n(x)$

上述公式中的 $epsilon$ 往往被替換成為 $-eta f(x)$

那麼，此時上述的公式就可以被寫成

$f(x-eta f(x)) approx f(x)-eta f(x) f』(x)+O(epsilon^2)$

因為，在定義中， $epsilon$ 往往就是一個非常小的值，那麼，其2次方，3次方，一直到n次方，那就幾乎等於0了。所以，我們把從第二次方項開始的所有結果的和定義為 $O(epsilon^2)$ ，這一項被稱作無窮小項。

那麼此時的式子就被簡化成為了： $f(x-eta f(x)) approx f(x)-eta f(x) f』(x)$

注意到 $[f(x)]^2$ 是一個大於等於0的數
$eta$ 也是一個大於0的數，那麼 $eta f(x) f』(x)$ 也是一個大於等於0的數
$x-eta f(x)$ 可以類比到開始我們講過的公式 $x := x-eta frac{partial Loss}{partial x}$
那麼，此時式子又可以被簡化為： $f(x-eta f(x)) leq f(x)$

其實，上面這個被簡化的式子就是在做一件事情：不斷的修改函數自變數的取值，使得其能夠滿足左邊的數值儘可能都小於右邊的值，也就是朝著 $f(x)$ 減少的方向一步一步的走。

（怎麼感覺自己說了句廢話）

我們還記著 $epsilon$ 的定義嗎？它是一個非常小的數值且接近0，可以通過 $-eta f(x)$ 來替換。

當 $f(x)$ 固定的時候，如果 $eta$ 是一個非常大的數值的時候，那麼此時的 $epsilon$ 就不能滿足非常小的定義。那麼，這個時候 $O(epsilon^2)$ 這一項就不能被drop掉，如果不能被drop掉，就不能滿足泰勒展開的基本定義。這也就對應了我們在訓練的過程中，如果把learning_rate（此處的 $eta$ 對應到learning_rate）調的太大，就出現了NaN的情況。
當 $eta$ 是一個固定的很小值的時候，如果 $x$ 的初始值選擇了一個很大的數值，比如說，這裡 $f(x)=x^2$ ，那麼此時的 $f(x) = 2x$ , $[f(x)]^2=4x^2$ ，這一項的數值也會變得非常大。導致其不能滿足 $epsilon$ 不是一個非常小的數值的條件。此時訓練肯定也會出問題。

總結：這裡就發現了導致model難以訓練的兩個問題：

學習率的大小選擇
初始化參數值的選擇

其實，這僅僅是從一維空間稍微為小夥伴們介紹了下「調參和其背後數學原理」的聯繫。。。

學習率

上述梯度下降演算法中的 $eta$ （取正數）叫做學習率或步長。

我們現在就來討論下，學習率過大和過小會帶來什麼問題。

當我們學習率太小的時候：

如上圖所示，當我們的初始值是從最左邊的這個值開始的時候，隨著 $x := x-eta frac{partial Loss}{partial x}$ 不斷迭代，使得 $x$ 不斷的增加，我們發現，其一步一步的走向了最低點，如果learning_rate設置的太小，那麼走的步數就會相當多，以至於很長時間，loss都沒辦法收斂。這就是所說的undershoot問題。

試想：如果以相同的學習率走了很多步後，終於走到了上圖中的最低點附近，那麼是不是很容易就能走到最低點了呢？

答案：不行的。因為如果我們不對learning_rate進行一個自適應的decay，比如說，隨著epoch的增多，learning_rate不斷的減小。是很難保證其最終能走到這個真正意義上的最低點的。

2. 當我們的學習率設置的過大

此時，就會產生一個走「之」字型的現象。

假設此時的初始值點是圖中左下角的點，那麼根據 $x := x-eta frac{partial Loss}{partial x}$ ，第一步就走到圖中右下角的這個點，然後不斷這樣迭代，有可能最後整個結果都沒辦法收斂下來。。。

這也就是在我們調參的過程中，設定太大的learning_rate會導致NaN的情況。

隨機梯度下降

然而，當訓練數據集很大的時候，梯度下降法可能會難以使用。

接下來，我們從數學原理方面來解釋下為什麼數據量變大的時候，往往不去採用梯度下降法

先來看看我們優化問題的目標函數 $f(x) = frac{1}{n}sum_{i=1}^nf_i(x)$
其中 $f_i(x)$ 表示的是第i個樣本所帶來的損失
可以觀察到，梯度下降每次進行迭代的開銷隨著n的增長呈線性增長。
因此，如果我們的訓練樣本個數非常多，那麼一次迭代的開銷也將非常大

針對這個問題，我們就引入了隨機梯度下降(sgd)方法，該方法從所有訓練樣本中隨機均勻採樣i，然後計算 $frac{partial f_i(x)}{partial x}$ 。

實際上，隨機梯度 $frac{partial f_i(x)}{partial x}$ 是對梯度 $frac{partial f(x)}{partial x}$ 的一個無偏估計，也就是說

$E[frac{partial f_i(x)}{partial x}] = frac{1}{n}sum_{i=1}^nfrac{partial f_i(x)}{partial x} = frac{partial f(x)}{partial x}$

Mini-Batch的隨機梯度下降

雖然已經有了梯度下降和隨機梯度下降，在實際的訓練過程中，我們還是更傾向於使用帶有mini-batch的sgd。

它就是說，隨機均勻採樣一個由訓練數據構成的小的batch。然後，通過這個batch上的所有樣本點的損失得到最終的損失，更新x的公式都和前面一樣，這裡就不寫了。。。

同樣，也可以通過數學來證明mini-batch的sgd對於原始梯度來說，也是無偏估計，證明方法同上。

演算法實現

其實，我們只要實現一個mini-batch的sgd就行了。

當batch_size == 1的時候，就是sgd
當batch_size == 整個訓練集大小的時候，就是梯度下降

def sgd( params, lr, batch_size ): for param in params: param[:] = param-lr*param.grad/batch_size

實驗中，我們以最簡單的線性回歸為例，完整代碼如下：

from mxnet import gluonfrom mxnet.gluon import nnfrom mxnet import ndarray as ndfrom mxnet import autogradimport mxnet as mxfrom time import timefrom mxnet import symimport matplotlib.pyplot as pltimport numpy as npimport randommx.random.seed(1)random.seed(1)def sgd( params, lr, batch_size ): for param in params: param[:] = param-lr*param.grad/batch_size# 生成數據集num_inputs = 2num_examples = 1000 true_w = [2,-3.4]true_b = 4.2X = nd.random_normal(scale = 1,shape = (num_examples,num_inputs))y = true_w[0]*X[:,0] + true_w[1]*X[:,1] + true_by += 0.01*nd.random_normal(scale = 1, shape = y.shape)dataset = gluon.data.ArrayDataset(X,y)def data_iter(batch_size): idx = list(range(num_examples)) random.shuffle(idx) for batch_i, i in enumerate(range(0,num_examples,batch_size)): j = nd.array(idx[i:min(i+batch_size,num_examples)]) yield batch_i,X.take(j),y.take(j)def init_params(): w = nd.random_normal(scale = 1,shape = (num_inputs,1)) b = nd.zeros(shape = (1,)) params = [w,b] for param in params: param.attach_grad() return paramsdef net(X,w,b): return nd.dot(X,w)+bdef square_loss(yhat,y): return (yhat-y.reshape(yhat.shape))**2/2def train(batch_size, lr, epochs, period): assert period >= batch_size and period % batch_size == 0 w, b = init_params() total_loss = [np.mean(square_loss(net(X, w, b), y).asnumpy())] # 注意epoch從1開始計數。 for epoch in range(1, epochs + 1): # 學習率自我衰減。 if epoch > 2: lr *= 0.1 for batch_i, data, label in data_iter(batch_size): with autograd.record(): output = net(data, w, b) loss = square_loss(output, label) loss.backward() sgd([w, b], lr, batch_size) if batch_i * batch_size % period == 0: total_loss.append( np.mean(square_loss(net(X, w, b), y).asnumpy())) print("Batch size %d, Learning rate %f, Epoch %d, loss %.4e" % (batch_size, lr, epoch, total_loss[-1])) print(w:, np.reshape(w.asnumpy(), (1, -1)), b:, b.asnumpy()[0], ) x_axis = np.linspace(0, epochs, len(total_loss), endpoint=True) plt.semilogy(x_axis, total_loss) plt.xlabel(epoch) plt.ylabel(loss) plt.show()

關於上述代碼中參數的說明，當epoch大於2的時候，學習率會有一個自我衰減的過程，也就是lr = lr *0.1。

period參數：每次採用到與period相同數目的數據點後，記錄當前目標函數值用於作圖。

比如，當batch_size = 10, period = 10，那麼，每次迭代一個batch後，都會記錄loss的值用於作圖。

調參

現在，我們來分析下不同參數對於loss變化的影響。

當batch_size = 1時，該訓練方式就是sgd，當前lr的情況下，loss function在前期快速減少，當epoch > 2的時候，lr會出現自我衰減，loss function下降後基本收斂，最終學習到的parameter和真實parameter相當。

train(batch_size = 1,lr = 0.1, epochs = 5, period = 10)

當batch_size = 1000時，由於訓練數據集包含1000個樣本，此時訓練使用的是標註的梯度下降演算法。loss function還是在前兩個epoch的時候下降較快。當epoch > 2的時候，lr自我衰減。loss function下降緩慢。最終學習到的parameter和真實parameter相當

當batch_size = 10時，由於訓練樣本中含有1000個樣本，此時訓練使用mini-batch的sgd來進行。最終學到的parameter與真實parameter相當

還是針對mini-batch的sgd，當我們把lr調大，那麼loss就會出現了nan的情況，也就是前面所分析的由於lr太大，導致loss不斷上升。

train(batch_size = 10, lr = 3, epochs = 5, period = 10)

還是針對mini-batch的sgd，當我們把lr調的非常小，那麼此時，已經經過了兩個epoch，loss下降還是非常緩慢的，這就是前面分析的lr太小，導致收斂過慢

動量法 momentum

動量法在前面的學習中應該有所接觸，但是理解的不深。通過Aston Zhang的講解，我對於為什麼要發明這個方法，以及這個方法所能帶來的好處有了更進一步的認識。

先來看下面的圖：

如果我們是從最左邊的眼睛看，那麼會看到一個變化非常「陡峭」的部分
如果我們是從最下面的眼睛看，那麼會看到一個相對「平緩」的部分

那麼，結合前面講過的sgd中lr的選擇問題

如果lr選擇的過大，就會出現overshot問題，模型loss出現nan，黑色線條部分
如果lr選擇的過小，就會出現undershot問題，迭代多步仍不收斂，紅色線條部分

那麼，我們就需要找到一個處於紅色部分和黑色部分的折中方式，使得其前進的方向朝著我們的optim值（紅色的點）能夠更快，更準的移動。

基於這個問題，就出現了動量法

知道了問題產生的背景，那麼我們接下里定義一些符號化的表示。

$v := gamma v+eta frac{partial f}{partial x}$

$x := x - v$

注意到，當 $gamma=0$ 的時候，那麼就退化成為了前面講過的梯度下降方法

現在，我們在對上面有關 $v$ 的公式進行改寫：

$v := gamma v + (1-gamma) frac{eta}{1-gamma} frac{partial f}{partial x}$

我們在引入一個叫做EMA（Estimation Moving Average）的東西，可以把上面這個有關 $v$ 的式子以另外一種更加優美的方式來表示

$y^{(t)} = gamma y^{(t+1)}+(1-gamma)x^{(t)}$

其實這兩種表示在某種意義上是等價的。。。

關於EMA，Aston Zhang通過虛擬貨幣交易市場中的例子來解釋的

就是說，在一個時時刻刻都在發生頻繁交易的市場中，我們想儘可能的通過一條光滑的曲線來模擬這種頻繁的抖動和變化。。。類似於moving average的思想。

此時，我們通過對式子中的參數賦予具體的值來達到更加具體化的學習。

$gamma = 0.95， t=20$

$egin{align*} y^{(20)} & = 0.05x^{(20)}+0.95y^{(19)} \ & = 0.05x^{(20)}+0.95*0.05x^{(19)}+0.95^2y^{(18)} \ & = 0.05x^{(20)}+0.95*0.05x^{(19)}+0.05*0.95^2x^{(18)}+0.95^3y^{(17)} \ & = ... end{align*}$

通過這樣一步一步的迭代下去。我們可以發現：

隨著迭代的進行， $x$ 的指數越大，那麼前面的係數就越小
$x$ 的指數越接近20的項，其前面的係數就越大

那麼，再來把上面涉及到的數學公式放在一起看：

$v := gamma v+eta frac{partial f}{partial x}$

$v := gamma v + (1-gamma) frac{eta}{1-gamma} frac{partial f}{partial x}$

$y^{(t)} = gamma y^{(t+1)}+(1-gamma)x^{(t)}$

就相當於對 $eta frac{partial f}{partial x}$ 變成了原來的 $frac{1}{1-gamma}$ 的倍數。

通過圖來形象的解釋就是：

其實，就是說，對於：

出現overshot的問題的時候，momentum方法可以對其進行一個裁剪。使得其在方差過大的走偏情況下，能夠盡量削弱這種情況，得到儘可能小的方差
出現undershot的問題的時候，momentum方法可以對其進行一個加速，使得這種趨勢能夠加強。

總結一下：加入momentum能夠加速達到optim的值。

現在用代碼實現下：

import mxnet as mxfrom mxnet import autogradfrom mxnet import gluonfrom mxnet import ndarray as ndimport numpy as npimport randommx.random.seed(1)random.seed(1)# 生成數據集。num_inputs = 2num_examples = 1000true_w = [2, -3.4]true_b = 4.2X = nd.random_normal(scale=1, shape=(num_examples, num_inputs))y = true_w[0] * X[:, 0] + true_w[1] * X[:, 1] + true_by += .01 * nd.random_normal(scale=1, shape=y.shape)dataset = gluon.data.ArrayDataset(X, y)net = gluon.nn.Sequential()net.add(gluon.nn.Dense(1))square_loss = gluon.loss.L2Loss()%matplotlib inlineimport matplotlib as mplmpl.rcParams[figure.dpi]= 120import matplotlib.pyplot as pltdef train(batch_size, lr, mom, epochs, period): assert period >= batch_size and period % batch_size == 0 net.collect_params().initialize(mx.init.Normal(sigma=1), force_reinit=True) # 動量法。 trainer = gluon.Trainer(net.collect_params(), sgd, {learning_rate: lr, momentum: mom}) data_iter = gluon.data.DataLoader(dataset, batch_size, shuffle=True) total_loss = [np.mean(square_loss(net(X), y).asnumpy())] for epoch in range(1, epochs + 1): # 重設學習率。 if epoch > 2: trainer.set_learning_rate(trainer.learning_rate * 0.1) for batch_i, (data, label) in enumerate(data_iter): with autograd.record(): output = net(data) loss = square_loss(output, label) loss.backward() trainer.step(batch_size) if batch_i * batch_size % period == 0: total_loss.append(np.mean(square_loss(net(X), y).asnumpy())) print("Batch size %d, Learning rate %f, Epoch %d, loss %.4e" % (batch_size, trainer.learning_rate, epoch, total_loss[-1])) print(w:, np.reshape(net[0].weight.data().asnumpy(), (1, -1)), b:, net[0].bias.data().asnumpy()[0], ) x_axis = np.linspace(0, epochs, len(total_loss), endpoint=True) plt.semilogy(x_axis, total_loss) plt.xlabel(epoch) plt.ylabel(loss) plt.show()

通過momentum的方法，我們先將mom/ $gamma$ 設置為0.9，此時正常，並且epoch>2後快速收斂了

train(batch_size=10, lr=0.2, mom=0.9, epochs=5, period=10)

再把 $gamma$ 設置的更大，0.99，此時梯度應該變為100倍，已經訓練飛了

Adagrad

在前面講過的這些優化演算法中，基本都是使用同一個learning_rate來更新所有的參數。

舉個二元函數的例子， $f(x_1,x_2)$ ，假設學習率為 $eta$ ，那麼參數的更新過程就是：

$x_1 := x_1 -eta frac{partial f}{partial x_1}$

$x_2 := x_2 -eta frac{partial f}{partial x_2}$

那麼，Adagrad要做的，就是對於不同的parameter，使用不同的learning_rate進行更新，並且其在迭代的過程中，能夠不斷自我調整learning_rate。

Adagrad演算法具體是這樣操作的：

使用一個梯度按元素平方的累加變數 $s:=s+g odot g$
其中 $g$ 就是通過mini-batch的計算得到的梯度
然後通過下面的式子對模型中每個參數的學習率通過按照元素重新調整 $g:=frac{eta}{sqrt{s+epsilon}}odot g$
其中 $eta$ 是初始學習率， $epsilon$ 是為了維持數值穩定性而添加的元素，防止分母除以0
然後再通過 $x:=x-g$ 對相應的parameter進行了更新

Adagrad演算法的核心思想：我們注意到 $s$ 其實是一個累加項的過程，

如果loss function相對於某一個parameter的偏導數一直都很大，那麼就讓他下降的快一點。
如果loss function相對於某一個parameter的偏導數一直都比較小，那麼就讓他下降的慢一點。

關於Adagrad的實現，我們只要將上述數學公式翻譯成python代碼即可：

def adagrad( params, sqrs, lr, batch_size): eps_stable = 1e-7 for param, sqr in zip(params,sqrs): g = param.grad/batch_size sqr[:] += nd.square(g) div = lr*g/nd.sqrt(sqr+eps_stable) param[:] -= div

整個程序的實現代碼和上述相似，這裡就不寫了，具體可以看看沐神的gluon教程。

RMSProp

在前面剛剛講過的Adagrad中，每個參數都有一個適應自己的learning_rate去更新，但是，當學習率在迭代早起降得比較快且這個時候的解依然比較不理想的時候，那麼有可能在就找不到一個更加理想的解了，也就是early stopping到了一個並不是我們認為最optim的點。

所以，RMSProp是對adagrad的一個改進。。。

其實，看到RMSProp的第一個式子，我就相當了Aston將的EMA的例子。看下面的公式

$s:=gamma s+(1-gamma ) g odot g$

相比adagrad演算法，RMSProp增加了一個衰減係數來控制對歷史信息的獲取多少。

剩下的公式和adagrad是一樣的。。。

寫成代碼：

def rmsprop(params, sqrs, lr, gamma, batch_size): eps_stable = 1e-8 for param, sqr in zip(params, sqrs): g = param.grad / batch_size sqr[:] = gamma * sqr + (1. - gamma) * nd.square(g) div = lr * g / nd.sqrt(sqr + eps_stable) param[:] -= div

當我們把這個係數 $gamma$ 設置為0.9，得到下面的曲線

再將 $gamma$ 設置的大一些，比如0.999，發現後期的曲線就比較平滑了

Adam

Adam演算法其實是前面講過的momentum方法和RMSProp方法的組合。

它使用了一個動量變數 $v$ 和一個RMSProp中梯度安裝元素平方的指數加權移動平均變數 $s$

$t:=t+1$ 每次的迭代
$v:=eta_1v+(1-eta_1)v$
$s:=eta_2s+(1-eta_2)g odot g$

在Adam演算法里，為了減輕 $v$ 和 $s$ 被初始化為0，在迭代初期對於計算指數加權移動平均的影響，進行了如下的修正：

$hat{v}:=frac{v}{1-eta^t_1}$
$hat{s}:=frac{s}{1-eta^t_2}$

在教程中講到， $eta_1=0.9,eta_2=0.999$ ，並且他們都在 $(0,1)$ 之間。當我們隨著迭代的進行， $t$ 比較大的時候，其對於 $v和s$ 將不會有太大的影響。

接下來通過：

$g:=frac{eta hat{v}}{sqrt{hat{s}+epsilon}}$
$x:=x-g$

Adam的實現：

def adam(params, vs, sqrs, lr, batch_size, t): beta1 = 0.9 beta2 = 0.999 eps_stable = 1e-8 for param, v, sqr in zip(params, vs, sqrs): g = param.grad / batch_size v[:] = beta1 * v + (1. - beta1) * g sqr[:] = beta2 * sqr + (1. - beta2) * nd.square(g) v_bias_corr = v / (1. - beta1 ** t) sqr_bias_corr = sqr / (1. - beta2 ** t) div = lr * v_bias_corr / (nd.sqrt(sqr_bias_corr) + eps_stable) param[:] = param - div

當學習率為0.1的時候，曲線圖如下

Adadelta

前面，我們已經介紹了有關momentum，RMSProp，Adagrad，Adam演算法，這些演算法有一個共性就是都帶有學習率，這一部分介紹的Adadelta演算法是沒有learning_rate這個參數的。

其計算過程和RMSProp一樣：

首先 $s:= ho s+(1- ho)g odot g$
然後計算需要更新的parameter的變化量： $g = frac{sqrt{Delta x+epsilon}}{sqrt{s+epsilon}}odot g$
這裡相比前面的演算法出現了 $Delta x$ ，初始化為零張量。
且其計算表達式為 $Delta x := ho Delta x + (1- ho)g odot g$
$x:=x-g$
$ho$ 的取值一般在0.9~0.999範圍內，當然也可以根據情況來調

代碼實現：

def adadelta(params, sqrs, deltas, rho, batch_size): eps_stable = 1e-5 for param, sqr, delta in zip(params, sqrs, deltas): g = param.grad / batch_size sqr[:] = rho * sqr + (1. - rho) * nd.square(g) cur_delta = nd.sqrt(delta + eps_stable) / nd.sqrt(sqr + eps_stable) * g delta[:] = rho * delta + (1. - rho) * cur_delta * cur_delta param[:] -= cur_delta

當 $ho = 0.999$ ，batch_size = 10時候的loss曲線圖

可以發現，整個過程中，沒有出現有關learning_rate的任何信息，同樣做到了很好地優化。。

有關優化的演算法，大體上就按照Aston zhang的講解介紹這麼多，希望大家在理解了基本的概念以及每一個優化演算法背後的原理後，在使用gluon的時候，就能「自信」的在trainer中設置自己想要的優化演算法了。