作者：chen_h
微訊號 & QQ：862251340
微信公眾號：coderpai
我的部落格：請點選這裡
這篇教程是翻譯Peter Roelants寫的迴圈神經網路教程，作者已經授權翻譯，這是原文。
該教程將介紹如何實現一個迴圈神經網路（RNN），一共包含兩部分。你可以在以下連結找到完整內容。
這篇教程中的程式碼是由 Python 2 IPython Notebook產生的，在教程的最後，我會給出全部程式碼的連結，幫助學習。神經網路中有關矩陣的運算我們採用NumPy來構建，畫圖使用Matplotlib來構建。如果你來沒有安裝這些軟體，那麼我強烈建議你使用Anaconda Python來安裝，這個軟體包中包含了執行這個教程的所有軟體包，非常方便使用。
迴圈神經網路
本教程主要包含三部分：
一個非常簡單的迴圈神經網路（RNN）
基於時序的反向傳播（BPTT）
彈性優化演算法
迴圈神經網路是一種可以解決序列資料的模型。在時序模型上面，這種迴圈關係可以定義成如下式子：
其中，Sk表示在時間k時刻的狀態，Xk是在時序k時刻的輸入資料，Wrec和Wx都是神經網路的連結權重。如果簡單的理解，可以把RNN理解成是一個帶反饋迴路的狀態模型。由於迴圈關係和延時處理，時序狀態被加入了模型之中。這個延時操作賦予了模型記憶力，因為它能記住模型前面一個狀態。
神經網路最後的輸出結果Yk是在時間k時刻計算出來的，即是通過前面一個或者多個狀態Sk，....，Sk+j計算出來的。
接下來，我們就可以通過輸入的資料Xk和前一步的狀態S(k-1)，來計算當前的狀態S(k)，或者通過輸入的資料Xk和前一步的狀態S(k)來預測下一步的狀態S(k+1)。
這篇教程會說明迴圈神經網路和一般的前饋神經網路沒有很大的不同，但是在訓練的方式上面可能會有一些不同。
線性迴圈神經網路
這部分教程我們來設計一個簡單的RNN模型，這個模型的輸入是一個二進位制的資料流，任務是去計算這個二進位制的資料流中存在幾個1。
在這個教程中，我們設計的RNN模型中的狀態只有一維，在每個時間點上，輸入資料也是一維的，最後輸出的結果就是序列狀態的最後一個狀態，即y = S(k)。我們將RNN模型進行展開，就可以得到下圖的模型。注意，展開的模型可以看做是一個 (n+1) 層的神經網路，每一層使用相同的連結權重Wrec和Wx。
雖然實現和訓練這個模型是一件非常有意思的事情，但是我們可以很容易得到，當W(rec) = W(x) = 1時，模型是最優的。
我們先匯入教程需要的軟體包
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import LogNorm複製程式碼定義資料集
輸入資料集 X 一共有20組資料，每組資料的長度是10，即每組資料的時間狀態步長是10。輸入資料是由均勻的隨機分佈產生的，取值 0 或者 1 。
輸出結果是輸入的二進位制資料流中存在幾個1，也就是把序列的每一位都加起來求和的結果。
# Create dataset
nb_of_samples = 20
sequence_len = 10
# Create the sequences
X = np.zeros((nb_of_samples, sequence_len))
for row_idx in range(nb_of_samples):
X[row_idx,:] = np.around(np.random.rand(sequence_len)).astype(int)
# Create the targets for each sequence
t = np.sum(X, axis=1)複製程式碼通過基於時序的反向傳播（BPTT）演算法進行訓練
訓練RNN的一個典型演算法是BPTT（backpropagation through time）演算法。通過名字，你也能發現這是一個基於BP的演算法。
如果你很瞭解常規的BP演算法，那麼BPTT演算法和常規的BP演算法沒有很大的不同。唯一的不同是，RNN需要每一個特定的時間步驟中，將每個神經元進行展開處理而已。展開圖已經在教程的最前面進行了說明。展開後，模型就和規則的神經網路模型很像了。唯一不同是，RNN有多個輸入源（前一個時間步驟的輸入狀態和當前的輸入資料）和每一層中的連結矩陣（ W(rec)和W(x) ）都是一樣的。
正向傳播計算RNN的輸出結果
正向傳播的時候，我們會把RNN展開進行處理，這樣就可以按照規則的神經網路進行處理了。RNN模型最後的輸出結果將會被使用在損失函式的計算中，用於訓練網路。（其實這些都和常規的多層神經網路一樣。）
當我們將RNN進行展開計算時，在不同的時間點上面，其實迴圈關係是相同的，我們將這個相同的迴圈關係在 update_state 函式中實現了。
forward_states函式通過 for 迴圈，將update_state函式應用到每一個時間點上面。如果我們將這些步驟都向量化，那麼就可以進行平行計算了。跟常規神經網路一樣，我們需要給權重進行初始化。在這個教程中，我們將權重初始化為0。
最後，我們通過累加所以輸入資料的誤差進行計算均方誤差函式（MSE）來得到損失函式 ξ 。在程式中，我們使用 cost 函式來實現。
# Define the forward step functions
def update_state(xk, sk, wx, wRec):
"""
Compute state k from the previous state (sk) and current input (xk),
by use of the input weights (wx) and recursive weights (wRec).
"""
return xk * wx + sk * wRec
def forward_states(X, wx, wRec):
"""
Unfold the network and compute all state activations given the input X,
and input weights (wx) and recursive weights (wRec).
Return the state activations in a matrix, the last column S[:,-1] contains the
final activations.
"""
# Initialise the matrix that holds all states for all input sequences.
# The initial state s0 is set to 0.
S = np.zeros((X.shape[0], X.shape[1]+1))
# Use the recurrence relation defined by update_state to update the
# states trough time.
for k in range(0, X.shape[1]):
# S[k] = S[k-1] * wRec + X[k] * wx
S[:,k+1] = update_state(X[:,k], S[:,k], wx, wRec)
return S
def cost(y, t):
"""
Return the MSE between the targets t and the outputs y.
"""
return ((t - y)**2).sum() / nb_of_samples複製程式碼反向傳播的梯度計算
在進行反向傳播過程之前，我們需要先計算誤差的對於輸出結果的梯度∂ξ/∂y，函式 output_gradient 實現了這個梯度計算過程。這個梯度將會被通過反向傳播演算法一層一層的向前傳播，函式 backward_gradient 實現了這個計算過程。具體的數學推導如下所示：
梯度最開始的計算公式為：
其中，n 表示RNN展開之後的時間步長。需要注意的是，引數 Wrec 擔當著反向傳遞誤差的角色。
損失函式對於權重的梯度是通過累加每一層中的梯度得到的。具體數學公式如下：
def output_gradient(y, t):
"""
Compute the gradient of the MSE cost function with respect to the output y.
"""
return 2.0 * (y - t) / nb_of_samples
def backward_gradient(X, S, grad_out, wRec):
"""
Backpropagate the gradient computed at the output (grad_out) through the network.
Accumulate the parameter gradients for wX and wRec by for each layer by addition.
Return the parameter gradients as a tuple, and the gradients at the output of each layer.
"""
# Initialise the array that stores the gradients of the cost with respect to the states.
grad_over_time = np.zeros((X.shape[0], X.shape[1]+1))
grad_over_time[:,-1] = grad_out
# Set the gradient accumulations to 0
wx_grad = 0
wRec_grad = 0
for k in range(X.shape[1], 0, -1):
# Compute the parameter gradients and accumulate the results.
wx_grad += np.sum(grad_over_time[:,k] * X[:,k-1])
wRec_grad += np.sum(grad_over_time[:,k] * S[:,k-1])
# Compute the gradient at the output of the previous layer
grad_over_time[:,k-1] = grad_over_time[:,k] * wRec
return (wx_grad, wRec_grad), grad_over_time複製程式碼梯度檢查
對於RNN，我們也需要對其進行梯度檢查，具體的檢查方法可以參考在常規多層神經網路中的梯度檢查。如果在反向傳播中的梯度計算正確，那麼這個梯度值應該和數值計算出來的梯度值應該是相同的。
# Perform gradient checking
# Set the weight parameters used during gradient checking
params = [1.2, 1.2] # [wx, wRec]
# Set the small change to compute the numerical gradient
eps = 1e-7
# Compute the backprop gradients
S = forward_states(X, params[0], params[1])
grad_out = output_gradient(S[:,-1], t)
backprop_grads, grad_over_time = backward_gradient(X, S, grad_out, params[1])
# Compute the numerical gradient for each parameter in the layer
for p_idx, _ in enumerate(params):
grad_backprop = backprop_grads[p_idx]
# + eps
params[p_idx] += eps
plus_cost = cost(forward_states(X, params[0], params[1])[:,-1], t)
# - eps
params[p_idx] -= 2 * eps
min_cost = cost(forward_states(X, params[0], params[1])[:,-1], t)
# reset param value
params[p_idx] += eps
# calculate numerical gradient
grad_num = (plus_cost - min_cost) / (2*eps)
# Raise error if the numerical grade is not close to the backprop gradient
if not np.isclose(grad_num, grad_backprop):
raise ValueError('Numerical gradient of {:.6f} is not close to the backpropagation gradient of {:.6f}!'.format(float(grad_num), float(grad_backprop)))
print('No gradient errors found')複製程式碼No gradient errors found
引數更新
由於不穩定的梯度，RNN是非常難訓練的。這也使得一般對於梯度的優化演算法，比如梯度下降，都不能使得RNN找到一個好的區域性最小值。
我們在下面的兩張圖中說明了RNN梯度的不穩定性。第一張圖表示，當我們給定 w(x) 和 w(rec) 時得到的損失表面圖。圖中帶顏色標記的地方，是我們取了幾個值做的實驗結果。從圖中，我們可以發現，當誤差表面的值接近於0時，w(x) = w(rec) = 1。但是當 |w(rec)| > 1時，誤差表面的值增加的非常迅速。
第二張圖我們通過幾組資料模擬了梯度的不穩定性，這個隨著時間步長而不穩定的梯度的形式和等比數列的形式很像，具體數學公式如下：
在狀態S(k)時的梯度，反向傳播m步得到的狀態S(k-m)可以被寫成：
在我們簡單的線性模型中，如果 |w(rec)| > 1，那麼梯度是一個指數爆炸的增長。如果 |w(rec)| < 1，那麼梯度將會消失。
關於指數暴漲，在第二張圖中，當我們取 w(x) =1, w(rec) = 2時，在圖中顯示梯度是指數爆炸增長的，當我們取 w(x) =1, w(rec) = -2時，正負徘徊指數增長，為什麼會出現徘徊？是因為我們把引數 w(rec) 取成了負數。這個指數爆炸說明了，模型的訓練對引數 w(rec) 是非常敏感的。
關於梯度消失，在第二張圖中，當我們取 w(x) = 1, w(rec) = 0.5和 w(x) = 1, w(rec) = -0.5時，那麼梯度將會指數下降，直至消失。這個梯度消失表示模型不能長時間的訓練，因為最後梯度將會消失。
如果 w(rec) = 0 時，梯度馬上變成了0。當 w(rec) = 1時，梯度隨著時間不變。
在下一部分，我們將說明怎麼去優化一個不穩定的誤差函式。
# Define plotting functions
# Define points to annotate (wx, wRec, color)
points = [(2,1,'r'), (1,2,'b'), (1,-2,'g'), (1,0,'c'), (1,0.5,'m'), (1,-0.5,'y')]
def get_cost_surface(w1_low, w1_high, w2_low, w2_high, nb_of_ws, cost_func):
"""Define a vector of weights for which we want to plot the cost."""
w1 = np.linspace(w1_low, w1_high, num=nb_of_ws) # Weight 1
w2 = np.linspace(w2_low, w2_high, num=nb_of_ws) # Weight 2
ws1, ws2 = np.meshgrid(w1, w2) # Generate grid
cost_ws = np.zeros((nb_of_ws, nb_of_ws)) # Initialize cost matrix
# Fill the cost matrix for each combination of weights
for i in range(nb_of_ws):
for j in range(nb_of_ws):
cost_ws[i,j] = cost_func(ws1[i,j], ws2[i,j])
return ws1, ws2, cost_ws
def plot_surface(ax, ws1, ws2, cost_ws):
"""Plot the cost in function of the weights."""
surf = ax.contourf(ws1, ws2, cost_ws, levels=np.logspace(-0.2, 8, 30), cmap=cm.pink, norm=LogNorm())
ax.set_xlabel('$w_{in}$', fontsize=15)
ax.set_ylabel('$w_{rec}$', fontsize=15)
return surf
def plot_points(ax, points):
"""Plot the annotation points on the given axis."""
for wx, wRec, c in points:
ax.plot(wx, wRec, c+'o', linewidth=2)
def get_cost_surface_figure(cost_func, points):
"""Plot the cost surfaces together with the annotated points."""
# Plot figures
fig = plt.figure(figsize=(10, 4))
# Plot overview of cost function
ax_1 = fig.add_subplot(1,2,1)
ws1_1, ws2_1, cost_ws_1 = get_cost_surface(-3, 3, -3, 3, 100, cost_func)
surf_1 = plot_surface(ax_1, ws1_1, ws2_1, cost_ws_1 + 1)
plot_points(ax_1, points)
ax_1.set_xlim(-3, 3)
ax_1.set_ylim(-3, 3)
# Plot zoom of cost function
ax_2 = fig.add_subplot(1,2,2)
ws1_2, ws2_2, cost_ws_2 = get_cost_surface(0, 2, 0, 2, 100, cost_func)
surf_2 = plot_surface(ax_2, ws1_2, ws2_2, cost_ws_2 + 1)
plot_points(ax_2, points)
ax_2.set_xlim(0, 2)
ax_2.set_ylim(0, 2)
# Show the colorbar
fig.subplots_adjust(right=0.8)
cax = fig.add_axes([0.85, 0.12, 0.03, 0.78])
cbar = fig.colorbar(surf_1, ticks=np.logspace(0, 8, 9), cax=cax)
cbar.ax.set_ylabel('$\\xi$', fontsize=15, rotation=0, labelpad=20)
cbar.set_ticklabels(['{:.0e}'.format(i) for i in np.logspace(0, 8, 9)])
fig.suptitle('Cost surface', fontsize=15)
return fig
def plot_gradient_over_time(points, get_grad_over_time):
"""Plot the gradients of the annotated point and how the evolve over time."""
fig = plt.figure(figsize=(6.5, 4))
ax = plt.subplot(111)
# Plot points
for wx, wRec, c in points:
grad_over_time = get_grad_over_time(wx, wRec)
x = np.arange(-grad_over_time.shape[1]+1, 1, 1)
plt.plot(x, np.sum(grad_over_time, axis=0), c+'-', label='({0}, {1})'.format(wx, wRec), linewidth=1, markersize=8)
plt.xlim(0, -grad_over_time.shape[1]+1)
# Set up plot axis
plt.xticks(x)
plt.yscale('symlog')
plt.yticks([10**8, 10**6, 10**4, 10**2, 0, -10**2, -10**4, -10**6, -10**8])
plt.xlabel('timestep k', fontsize=12)
plt.ylabel('$\\frac{\\partial \\xi}{\\partial S_{k}}$', fontsize=20, rotation=0)
plt.grid()
plt.title('Unstability of gradient in backward propagation.\n(backpropagate from left to right)')
# Set legend
leg = plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False, numpoints=1)
leg.set_title('$(w_x, w_{rec})$', prop={'size':15})
def get_grad_over_time(wx, wRec):
"""Helper func to only get the gradient over time from wx and wRec."""
S = forward_states(X, wx, wRec)
grad_out = output_gradient(S[:,-1], t).sum()
_, grad_over_time = backward_gradient(X, S, grad_out, wRec)
return grad_over_time複製程式碼# Plot cost surface and gradients
# Get and plot the cost surface figure with markers
fig = get_cost_surface_figure(lambda w1, w2: cost(forward_states(X, w1, w2)[:,-1] , t), points)
# Get the plots of the gradients changing by backpropagating.
plot_gradient_over_time(points, get_grad_over_time)
# Show figures
plt.show()複製程式碼
彈性優化演算法
在上面的部分，我們已經介紹了RNN的梯度是非常不穩定的，所以梯度在損失表面的跳躍度是非常大的，也就是說優化程式可能將最優值帶到離真實最優值很遠的地方，如下圖：
根據在我們神經網路裡面的基礎教程，梯度下降法更新引數的公式如下：
其中，W(i) 表示在第 i 次迭代時 W 的值，μ 是學習率。
在訓練過程中，當我們取 w(x) = 1 和 w(rec) = 2時，誤差表面上的藍色點的梯度值將達到 10^7。儘管我們把學習率取的非常小，比如0.000001(1e-6)，但是引數 W 也將離開原來的距離 10 個單位，在我們的模型中，這將會導致災難性的結果。一個解決方案是我們再降低學習率的值，但是這樣做將導致，當梯度很小時，更新的點將保持在原地不動。
對於這個問題，研究者們已經找到了很多的方法來解決不穩定的梯度，比如Gradient clipping，Hessian-Free Optimization，Momentum。
我們可以使用一些優化演算法來處理這個不穩定梯度，以此來減小梯度的敏感度。其中一個技術就是使用彈性反向傳播（Rprop）。彈性反向傳播演算法和之前教程中的動量演算法非常相似，但是這裡只是用在梯度上面，用來更新引數。Rprop演算法描述如下：
一般情況下，模型的超引數被設定為 η^+ = 1.2 和 η^- = 0.5 。如果我們將這個Rprop演算法和之前的動量演算法進行對比的話，我們可以發現：當梯度的符合不改變時，我們將增加 20% 的權重；當梯度的符合改變時，我們將減小 50% 的權重。注意，Rprop演算法的更新值 Δ 類似於動量中的速度引數。不同點是Rprop演算法的值只是反映了動量中的速度的值，不包括方向。方向是由當前梯度的方向來決定的。
在這個教程中，我們迭代這個Rprop演算法 500 次。下圖中的藍色點就是在誤差表面的更新值。注意圖中，儘管權重引數開始的位置是在一個很高的誤差值和一個很高的梯度位置，但是在我們的迭代最後，Rprop演算法還是將最優值鎖定在座標 (1, 1) 左右。
# Define Rprop optimisation function
def update_rprop(X, t, W, W_prev_sign, W_delta, eta_p, eta_n):
"""
Update Rprop values in one iteration.
X: input data.
t: targets.
W: Current weight parameters.
W_prev_sign: Previous sign of the W gradient.
W_delta: Rprop update values (Delta).
eta_p, eta_n: Rprop hyperparameters.
"""
# Perform forward and backward pass to get the gradients
S = forward_states(X, W[0], W[1])
grad_out = output_gradient(S[:,-1], t)
W_grads, _ = backward_gradient(X, S, grad_out, W[1])
W_sign = np.sign(W_grads) # Sign of new gradient
# Update the Delta (update value) for each weight parameter seperately
for i, _ in enumerate(W):
if W_sign[i] == W_prev_sign[i]:
W_delta[i] *= eta_p
else:
W_delta[i] *= eta_n
return W_delta, W_sign複製程式碼# Perform Rprop optimisation
# Set hyperparameters
eta_p = 1.2
eta_n = 0.5
# Set initial parameters
W = [-1.5, 2] # [wx, wRec]
W_delta = [0.001, 0.001] # Update values (Delta) for W
W_sign = [0, 0] # Previous sign of W
ls_of_ws = [(W[0], W[1])] # List of weights to plot
# Iterate over 500 iterations
for i in range(500):
# Get the update values and sign of the last gradient
W_delta, W_sign = update_rprop(X, t, W, W_sign, W_delta, eta_p, eta_n)
# Update each weight parameter seperately
for i, _ in enumerate(W):
W[i] -= W_sign[i] * W_delta[i]
ls_of_ws.append((W[0], W[1])) # Add weights to list to plot
print('Final weights are: wx = {0}, wRec = {1}'.format(W[0], W[1]))複製程式碼Final weights are: wx = 1.00135554721, wRec = 0.999674473785
# Plot the cost surface with the weights over the iterations.
# Define plot function
def plot_optimisation(ls_of_ws, cost_func):
"""Plot the optimisation iterations on the cost surface."""
ws1, ws2 = zip(*ls_of_ws)
# Plot figures
fig = plt.figure(figsize=(10, 4))
# Plot overview of cost function
ax_1 = fig.add_subplot(1,2,1)
ws1_1, ws2_1, cost_ws_1 = get_cost_surface(-3, 3, -3, 3, 100, cost_func)
surf_1 = plot_surface(ax_1, ws1_1, ws2_1, cost_ws_1 + 1)
ax_1.plot(ws1, ws2, 'b.')
ax_1.set_xlim([-3,3])
ax_1.set_ylim([-3,3])
# Plot zoom of cost function
ax_2 = fig.add_subplot(1,2,2)
ws1_2, ws2_2, cost_ws_2 = get_cost_surface(0, 2, 0, 2, 100, cost_func)
surf_2 = plot_surface(ax_2, ws1_2, ws2_2, cost_ws_2 + 1)
ax_2.set_xlim([0,2])
ax_2.set_ylim([0,2])
surf_2 = plot_surface(ax_2, ws1_2, ws2_2, cost_ws_2)
ax_2.plot(ws1, ws2, 'b.')
# Show the colorbar
fig.subplots_adjust(right=0.8)
cax = fig.add_axes([0.85, 0.12, 0.03, 0.78])
cbar = fig.colorbar(surf_1, ticks=np.logspace(0, 8, 9), cax=cax)
cbar.ax.set_ylabel('$\\xi$', fontsize=15)
cbar.set_ticklabels(['{:.0e}'.format(i) for i in np.logspace(0, 8, 9)])
plt.suptitle('Cost surface', fontsize=15)
plt.show()
# Plot the optimisation
plot_optimisation(ls_of_ws, lambda w1, w2: cost(forward_states(X, w1, w2)[:,-1] , t))
plt.show()複製程式碼
測試模型
最後我們編寫測試程式碼。從程式碼的執行中，我們能發現目標值和真實值非常的相近。如果我們取模型輸出值的最靠近的整數，那麼預測值的輸出將更加完美。
test_inpt = np.asmatrix([[0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1]])
print test_inpt
test_outpt = forward_states(test_inpt, W[0], W[1])[:,-1]
print 'Target output: {:d} vs Model output: {:.2f}'.format(test_inpt.sum(), test_outpt[0])複製程式碼Target output: 5 vs Model output: 4.99
CoderPai 是一個專注於演算法實戰的平臺，從基礎的演算法到人工智慧演算法都有設計。如果你對演算法實戰感興趣，請快快關注我們吧。加入AI實戰微信群，AI實戰QQ群，ACM演算法微信群，ACM演算法QQ群。詳情請關注 “CoderPai” 微訊號（coderpai）。