技術標籤：DeepLearning學習python深度學習機器學習

動量法

目標函式有關自變數的梯度代表了目標函式在自變數當前位置下降最快的方向，因此，梯度下降也叫作最陡下降（steepest descent）。在每次迭代中，梯度下降根據自變數當前位置，沿著當前位置的梯度更新自變數。然而，如果自變數的迭代方向僅僅取決於自變數當前位置，這可能會帶來一些問題。

梯度下降的問題

考慮一個目標函式：

• 輸入為二維向量 x = [ x 1 , x 2 ] ⊤ \boldsymbol{x} = [x_1, x_2]^\top x=[x1​,x2​]⊤
• 輸出為標量
• f ( x ) = 0.1 x 1 2 + 2 x 2 2 f(\boldsymbol{x})=0.1x_1^2+2x_2^2
  f(x)=0.1x12​+2x22​。

這裡將 x 1 2 x_1^2 x12​係數從 1 1 1減小到了 0.1 0.1 0.1。下面實現基於這個目標函式的梯度下降，並演示使用學習率為 0.4 0.4 0.4時自變數的迭代軌跡。

%matplotlib inline
import sys
sys.path.append("..") 
from matplotlib import pyplot as plt
import torch

eta = 0.4 # 學習率

def f_2d(x1, x2):
    return 0.1 * x1 ** 2 + 2 * x2 ** 2

def gd_2d
 
(x1, x2, s1, s2):
    return (x1 - eta * 0.2 * x1, x2 - eta * 4 * x2, 0, 0)
    
def show_trace_2d(f, results):  
    plt.plot(*zip(*results), '-o', color='#ff7f0e')
    x1, x2 = np.meshgrid(np.arange(-5.5, 1.0, 0.1), np.arange(-3.0, 1.0, 0.1))
    plt.contour(x1, x2, f(x1, x2), colors='#1f77b4')
    plt.xlabel(
 
'x1')
    plt.ylabel('x2')

def train_2d(trainer):  
    x1, x2, s1, s2 = -5, -2, 0, 0  # s1和s2是自變數狀態
    results = [(x1, x2)]
    for i in range(20):
        x1, x2, s1, s2 = trainer(x1, x2, s1, s2)
        results.append((x1, x2))
    print('epoch %d, x1 %f, x2 %f' % (i + 1, x1, x2))
    return results
    
show_trace_2d(f_2d, train_2d(gd_2d))

輸出：

epoch 20, x1 -0.943467, x2 -0.000073

這幅圖可以理解為 f ( x ) f(\boldsymbol{x}) f(x)在 x 1 x_1 x1​、 x 2 x_2 x2​平面上的投影，可以看到

• 同一位置上，目標函式在豎直方向（ x 2 x_2 x2​軸方向）比在水平方向（ x 1 x_1 x1​軸方向）的斜率的絕對值更大。
• 因此，給定學習率，梯度下降迭代自變數時會使自變數在豎直方向比在水平方向移動幅度更大。
• 我們需要一個較小的學習率從而避免自變數在豎直方向上越過目標函式最優解。然而，這會造成自變數在水平方向上朝最優解移動變慢。

下面我們試著將學習率調得稍大一點，此時自變數在豎直方向不斷越過最優解並逐漸發散。

eta = 0.6
show_trace_2d(f_2d, train_2d(gd_2d))

動量法

動量法的提出是為了解決梯度下降的上述問題。

• 設時間步 t t t的自變數為 x t \boldsymbol{x}_t xt​
• 學習率為 η t \eta_t ηt​
• 在時間步 0 0 0，動量法建立速度變數 v 0 \boldsymbol{v}_0 v0​，並將其元素初始化成0。在時間步 t > 0 t>0 t>0，動量法對每次迭代的步驟做如下修改：

v t ← γ v t − 1 + η t g t , x t ← x t − 1 − v t , \begin{aligned} \boldsymbol{v}_t &\leftarrow \gamma \boldsymbol{v}_{t-1} + \eta_t \boldsymbol{g}_t, \\ \boldsymbol{x}_t &\leftarrow \boldsymbol{x}_{t-1} - \boldsymbol{v}_t, \end{aligned} vt​xt​​←γvt−1​+ηt​gt​,←xt−1​−vt​,​

其中，動量超引數 γ \gamma γ滿足 0 ≤ γ < 1 0 \leq \gamma < 1 0≤γ<1。當 γ = 0 \gamma=0 γ=0時，動量法等價於小批量隨機梯度下降。

在解釋動量法的數學原理前，讓我們先從實驗中觀察梯度下降在使用動量法後的迭代軌跡。

def momentum_2d(x1, x2, v1, v2):
    v1 = gamma * v1 + eta * 0.2 * x1
    v2 = gamma * v2 + eta * 4 * x2
    return x1 - v1, x2 - v2, v1, v2

eta, gamma = 0.4, 0.5
show_trace_2d(f_2d, train_2d(momentum_2d))

可以看到，使用較小的學習率 η = 0.4 \eta=0.4 η=0.4和動量超引數 γ = 0.5 \gamma=0.5 γ=0.5時

• 動量法在豎直方向上的移動更加平滑
• 且在水平方向上更快逼近最優解

下面使用較大的學習率 η = 0.6 \eta=0.6 η=0.6，此時自變數也不再發散。

eta = 0.6
show_trace_2d(f_2d, train_2d(momentum_2d))

指數加權移動平均

為了從數學上理解動量法，需要先解釋一下指數加權移動平均（exponentially weighted moving average）。

• 給定超引數 0 ≤ γ < 1 0 \leq \gamma < 1 0≤γ<1
• 當前時間步 t t t的變數 y t y_t yt​
• 上一時間步 t − 1 t-1 t−1的變數 y t − 1 y_{t-1} yt−1​
• 當前時間步另一變數 x t x_t xt​

當前時間步 t t t的變數 y t y_t yt​是上一時間步 t − 1 t-1 t−1的變數 y t − 1 y_{t-1} yt−1​與當前時間步另一變數 x t x_t xt​的線性組合：

y t = γ y t − 1 + ( 1 − γ ) x t . y_t = \gamma y_{t-1} + (1-\gamma) x_t. yt​=γyt−1​+(1−γ)xt​.

我們可以對 y t y_t yt​展開：

y t = ( 1 − γ ) x t + γ y t − 1 = ( 1 − γ ) x t + ( 1 − γ ) ⋅ γ x t − 1 + γ 2 y t − 2 = ( 1 − γ ) x t + ( 1 − γ ) ⋅ γ x t − 1 + ( 1 − γ ) ⋅ γ 2 x t − 2 + γ 3 y t − 3 … \begin{aligned} y_t &= (1-\gamma) x_t + \gamma y_{t-1}\\ &= (1-\gamma)x_t + (1-\gamma) \cdot \gamma x_{t-1} + \gamma^2y_{t-2}\\ &= (1-\gamma)x_t + (1-\gamma) \cdot \gamma x_{t-1} + (1-\gamma) \cdot \gamma^2x_{t-2} + \gamma^3y_{t-3}\ &\ldots \end{aligned} yt​​=(1−γ)xt​+γyt−1​=(1−γ)xt​+(1−γ)⋅γxt−1​+γ2yt−2​=(1−γ)xt​+(1−γ)⋅γxt−1​+(1−γ)⋅γ2xt−2​+γ3yt−3​​…​

令 n = 1 / ( 1 − γ ) n = 1/(1-\gamma) n=1/(1−γ)，那麼 ( 1 − 1 / n ) n = γ 1 / ( 1 − γ ) \left(1-1/n\right)^n = \gamma^{1/(1-\gamma)} (1−1/n)n=γ1/(1−γ)。因為

lim ⁡ n → ∞ ( 1 − 1 n ) n = exp ⁡ ( − 1 ) ≈ 0.3679 , \lim_{n \rightarrow \infty} \left(1-\frac{1}{n}\right)^n = \exp(-1) \approx 0.3679, n→∞lim​(1−n1​)n=exp(−1)≈0.3679,

所以當 γ → 1 \gamma \rightarrow 1 γ→1時， γ 1 / ( 1 − γ ) = exp ⁡ ( − 1 ) \gamma^{1/(1-\gamma)}=\exp(-1) γ1/(1−γ)=exp(−1)，如 0.9 5 20 ≈ exp ⁡ ( − 1 ) 0.95^{20} \approx \exp(-1) 0.9520≈exp(−1)。如果把 exp ⁡ ( − 1 ) \exp(-1) exp(−1)當作一個比較小的數，我們可以在近似中忽略所有含 γ 1 / ( 1 − γ ) \gamma^{1/(1-\gamma)} γ1/(1−γ)和比 γ 1 / ( 1 − γ ) \gamma^{1/(1-\gamma)} γ1/(1−γ)更高階的係數的項。例如，當 γ = 0.95 \gamma=0.95 γ=0.95時，

y t ≈ 0.05 ∑ i = 0 19 0.9 5 i x t − i . y_t \approx 0.05 \sum_{i=0}^{19} 0.95^i x_{t-i}. yt​≈0.05i=0∑19​0.95ixt−i​.

因此，在實際中，我們常常將 y t y_t yt​看作是對最近 1 / ( 1 − γ ) 1/(1-\gamma) 1/(1−γ)個時間步的 x t x_t xt​值的加權平均。例如，

• 當 γ = 0.95 \gamma = 0.95 γ=0.95時， y t y_t yt​可以被看作對最近20個時間步的 x t x_t xt​值的加權平均
• 當 γ = 0.9 \gamma = 0.9 γ=0.9時， y t y_t yt​可以看作是對最近10個時間步的 x t x_t xt​值的加權平均
• 而且，離當前時間步 t t t越近的 x t x_t xt​值獲得的權重越大（越接近1）

由指數加權移動平均理解動量法

現在，我們對動量法的速度變數做變形：

v t ← γ v t − 1 + ( 1 − γ ) ( η t 1 − γ g t ) . \boldsymbol{v}_t \leftarrow \gamma \boldsymbol{v}{t-1} + (1 - \gamma) \left(\frac{\eta_t}{1 - \gamma} \boldsymbol{g}_t\right). vt​←γvt−1+(1−γ)(1−γηt​​gt​).

由指數加權移動平均的形式可得，速度變數 v t \boldsymbol{v}_t vt​實際上對序列 { η t − i g t − i / ( 1 − γ ) : i = 0 , … , 1 / ( 1 − γ ) − 1 } \{{\eta_{t-i}\boldsymbol{g}_{t-i} /(1-\gamma):i=0,\ldots,1/(1-\gamma)-1}\} {ηt−i​gt−i​/(1−γ):i=0,…,1/(1−γ)−1}做了指數加權移動平均。

• 換句話說，相比於小批量隨機梯度下降，動量法在每個時間步的自變數更新量近似於將最近 1 / ( 1 − γ ) 1/(1-\gamma) 1/(1−γ)個時間步的普通更新量（即學習率乘以梯度）做了指數加權移動平均後再除以 1 − γ 1-\gamma 1−γ

所以，在動量法中，自變數在各個方向上的移動幅度不僅取決當前梯度，還取決於過去的各個梯度在各個方向上是否一致。

之前的優化問題中，所有梯度在水平方向上為正（向右），而在豎直方向上時正（向上）時負（向下）。這樣，我們就可以使用較大的學習率，從而使自變數向最優解更快移動。

實現動量法

相對於小批量隨機梯度下降，動量法需要對每一個自變數維護一個同它一樣形狀的速度變數，且超引數裡多了動量超引數。實現中，將速度變數用更廣義的狀態變數states表示。

def get_data_ch7():  
    data = np.genfromtxt('../../data/airfoil_self_noise.dat', delimiter='\t')
    data = (data - data.mean(axis=0)) / data.std(axis=0)
    return torch.tensor(data[:1500, :-1], dtype=torch.float32), \
        torch.tensor(data[:1500, -1], dtype=torch.float32) # 前1500個樣本(每個樣本5個特徵)

def train_ch7(optimizer_fn, states, hyperparams, features, labels,
              batch_size=10, num_epochs=2):
    # 初始化模型
    net, loss = linreg, squared_loss
    
    w = torch.nn.Parameter(torch.tensor(np.random.normal(0, 0.01, size=(features.shape[1], 1)), dtype=torch.float32),
                           requires_grad=True)
    b = torch.nn.Parameter(torch.zeros(1, dtype=torch.float32), requires_grad=True)

    def eval_loss():
        return loss(net(features, w, b), labels).mean().item()

    ls = [eval_loss()]
    data_iter = torch.utils.data.DataLoader(
        torch.utils.data.TensorDataset(features, labels), batch_size, shuffle=True)

features, labels = get_data_ch7()

def init_momentum_states():
    v_w = torch.zeros((features.shape[1], 1), dtype=torch.float32)
    v_b = torch.zeros(1, dtype=torch.float32)
    return (v_w, v_b)

def sgd_momentum(params, states, hyperparams):
    for p, v in zip(params, states):
        v.data = hyperparams['momentum'] * v.data + hyperparams['lr'] * p.grad.data
        p.data -= v.data

先將動量超引數momentum設0.5，這時可以看成是特殊的小批量隨機梯度下降：

• 其小批量隨機梯度為最近2個時間步的2倍小批量梯度的加權和。

train_ch7(sgd_momentum, init_momentum_states(),
              {'lr': 0.02, 'momentum': 0.5}, features, labels)

將動量超引數momentum增大到0.9，這時依然可以看成是特殊的小批量隨機梯度下降：

• 其小批量隨機梯度為最近10個時間步的10倍小批量梯度的加權和。我們先保持學習率0.02不變

目標函式值在後期迭代過程中的變化不夠平滑。直覺上，10倍小批量梯度比2倍小批量梯度大了5倍，試著將學習率減小到原來的1/5。此時目標函式值在下降了一段時間後變化更加平滑。

train_ch7(sgd_momentum, init_momentum_states(),
              {'lr': 0.004, 'momentum': 0.9}, features, labels)

• 動量法使用了指數加權移動平均的思想。它將過去時間步的梯度做了加權平均，且權重按時間步指數衰減。
• 動量法使得相鄰時間步的自變數更新在方向上更加一致。