摘要：線性循環神經網絡這部分教程我們來設計一個簡單的模型，這個模型的輸入是一個二進制的數據流，任務是去計算這個二進制的數據流中存在幾個。

作者：chen_h
微信號 & QQ：862251340
微信公眾號：coderpai
簡書地址：https://www.jianshu.com/p/160...

這篇教程是翻譯Peter Roelants寫的循環神經網絡教程，作者已經授權翻譯，這是原文。

該教程將介紹如何實現一個循環神經網絡（RNN），一共包含兩部分。你可以在以下鏈接找到完整內容。

（一）線性循環神經網絡（RNN）

（二）非線性循環神經網絡（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)

訓練RNN的一個典型算法是BPTT（backpropagation through time）算法。通過名字，你也能發現這是一個基于BP的算法。

如果你很了解常規的BP算法，那么BPTT算法和常規的BP算法沒有很大的不同。唯一的不同是，RNN需要每一個特定的時間步驟中，將每個神經元進行展開處理而已。展開圖已經在教程的最前面進行了說明。展開后，模型就和規則的神經網絡模型很像了。唯一不同是，RNN有多個輸入源（前一個時間步驟的輸入狀態和當前的輸入數據）和每一層中的鏈接矩陣（ W(rec)和W(x) ）都是一樣的。

正向傳播的時候，我們會把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.
(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

完整代碼，點擊這里

作者：chen_h
微信號 & QQ：862251340
簡書地址：https://www.jianshu.com/p/160...

CoderPai 是一個專注于算法實戰的平臺，從基礎的算法到人工智能算法都有設計。如果你對算法實戰感興趣，請快快關注我們吧。加入AI實戰微信群，AI實戰QQ群，ACM算法微信群，ACM算法QQ群。長按或者掃描如下二維碼，關注 “CoderPai” 微信號（coderpai）