在《神經網路的梯度推導與程式碼驗證》之vanilla RNN的前向傳播和反向梯度推導中，我們學習了vanilla RNN的前向傳播和反向梯度求導，但知識仍停留在紙面。本篇章將基於深度學習框架tensorflow驗證我們所得結論的準確性，以便將抽象的數學符號和實際資料結合起來，將知識固化。更多相關內容請見《神經網路的梯度推導與程式碼驗證》系列介紹。
提醒：
- 後續會反覆出現$\boldsymbol{\delta}^{l}$這個（類）符號，它的定義為$\boldsymbol{\delta}^{l} = \frac{\partial l}{\partial\boldsymbol{z}^{\boldsymbol{l}}}$，即loss $l$對$\boldsymbol{z}^{\boldsymbol{l}}$的導數
- 其中$\boldsymbol{z}^{\boldsymbol{l}}$表示第$l$層（DNN，CNN，RNN或其他例如max pooling層等）未經過啟用函式的輸出。
- $\boldsymbol{a}^{\boldsymbol{l}}$則表示$\boldsymbol{z}^{\boldsymbol{l}}$經過啟用函式後的輸出。
這些符號會貫穿整個系列，還請留意。
需要用到的庫有tensorflow和numpy，其中tensorflow其實版本>=2.0.0就行
def get_crossentropy(y_pred, y_true): return -tf.reduce_sum(y_true * tf.math.log(y_pred))
然後是定義交叉熵損失函式：
def get_crossentropy(y_pred, y_true): return -tf.reduce_sum(y_true * tf.math.log(y_pred))
--------前向傳播驗證---------
下面開始實現前向傳播：
我們先來看如果拿tensorflow快速實現前向傳播是什麼樣的，程式碼挺短的：
1 y_true = np.array([[[0.3, 0.5, 0.2], 2 [0.2, 0.3, 0.5], 3 [0.5, 0.2, 0.3]]]).astype(np.float32) 4 5 # --------inputs--------- 6 inputs = np.random.random([1, 3, 4]).astype(np.float32) 7 init_state = [tf.constant(np.random.random((inputs.shape[0], 2)).astype(np.float32))] 8 # --------vanilla rnn--------- 9 rnn = tf.keras.layers.RNN(tf.keras.layers.SimpleRNNCell(2), return_sequences=True) 10 h_seq = rnn(inputs=inputs, initial_state=init_state) 11 # --------fnn------------- 12 dense = tf.keras.layers.Dense(3) 13 output_dense = dense(h_seq) 14 output_seq = tf.math.softmax(output_dense)
先看樣本(inputs, y_ture)，輸入inputs是一條步長為3的有4維特徵的資料；標籤資料同樣也只有一條，步長也為3，每個步長上是長度為3的概率向量。
inputs.shape Out[5]: (1, 3, 4) y_true.shape Out[6]: (1, 3, 3)
再看我們的init_state，它就是在計算$\boldsymbol{h}^{(1)}$的時候顯然我們需要的$\boldsymbol{h}^{(0)}$，關於$\boldsymbol{h}^{(t)}$的計算公式見下面這個式子：
$\boldsymbol{h}^{(t)} = \sigma\left( \boldsymbol{z}^{(t)} \right) = \sigma\left( {\boldsymbol{U}\boldsymbol{x}^{(t)} + \boldsymbol{W}\boldsymbol{h}^{(t - 1)} + \boldsymbol{b}} \right)$
接著就是建立一個vanilla RNN 單元 tf.keras.layers.SimpleRNNCell(2)，它的輸出是2維的。外層還要包一個tf.keras.layers.RNN(tf.keras.layers.SimpleRNNCell(2), return_sequences=True)，這樣才算是實現了vanilla RNN層。
我們看看vanilla rnn的輸出h_seq：
h_seq Out[8]: <tf.Tensor: shape=(1, 3, 2), dtype=float32, numpy= array([[[ 0.85244656, 0.01066787], [ 0.3758163 , 0.21824013], [ 0.8450084 , -0.6304215 ]]], dtype=float32)>
h_seq = rnn(inputs=inputs, initial_state=init_state) 這句話在做的就是下圖右邊的那部分的操作：
- init_state對應著$\boldsymbol{h}^{(...)}$，i
- inputs[0, 0, :]對應著$\boldsymbol{x}^{(t-1)}$，inputs[0, 1, :]對應著$\boldsymbol{x}^{(t)}$，inputs[0, 2, :]對應著$\boldsymbol{x}^{(t+1)}$
- h_seq[0, 0, :]對應著$\boldsymbol{h}^{(t-1)}$，h_seq[0, 1, :]對應著$\boldsymbol{h}^{(t)}$，h_seq[0, 2, :]對應著$\boldsymbol{h}^{(t+1)}$
h狀態經過FNN之後得到下面的output_seq：
output_seq Out[11]: <tf.Tensor: shape=(1, 3, 3), dtype=float32, numpy= array([[[0.43937117, 0.37462497, 0.18600382], [0.32436833, 0.3793582 , 0.29627347], [0.6205071 , 0.2681237 , 0.11136916]]], dtype=float32)>
其中，output_seq[0, 0, :]對應著$\boldsymbol{o}^{(t-1)}$，output_seq[0, 1, :]對應著$\boldsymbol{o}^{(t)}$，output_seq[0, 2, :]對應著$\boldsymbol{o}^{(t+1)}$
通過呼叫上面幾行程式碼，我們似乎能夠實現上圖的操作。但這還不夠細緻入微，因為目前充其量只是稍微驗證了下輸入和輸出的shape而已，我們尚未確定程式碼 rnn = tf.keras.layers.RNN(tf.keras.layers.SimpleRNNCell(2), return_sequences=True)是否按正真按照vanilla RNN的前向傳播和反向梯度推導中給出的前傳公式進行的。所以接下來我們手動按照vanilla RNN的前傳公式實現一遍前向傳播看跟tensorflow給出的輸出結果是否一致。
下面這段程式碼看似很長，但實際上只是反覆做同樣的事情而已，它實現了一個vanilla RNN在時間上展開3步的前向操作。
這裡 tf.GradientTape(persistent=True) ，t.watch()是用於後面計算變數的導數用的，不太熟悉的可參考tensorflow官方給出的關於這部分的教程（自動微分）。
1 with tf.GradientTape(persistent=True) as t2: 2 # -------rnn analysis by steps--------- 3 # 下面是手動演算的rnn展開，和上面的whole_sequence_output是對得上的 4 5 # ------time stpe 1-------- 6 z1 = tf.matmul(inputs[:, 0, :], rnn.weights[0]) + tf.matmul(init_state[0], rnn.weights[1]) + rnn.weights[2] 7 t2.watch(z1) 8 h1 = tf.math.tanh(z1) 9 t2.watch(h1) 10 out1 = dense(h1) 11 t2.watch(out1) 12 a1 = tf.math.softmax(out1) 13 t2.watch(a1) 14 # ------time stpe 2-------- 15 z2 = tf.matmul(inputs[:, 1, :], rnn.weights[0]) + tf.matmul(h1, rnn.weights[1]) + rnn.weights[2] 16 t2.watch(z2) 17 h2 = tf.math.tanh(z2) 18 t2.watch(h2) 19 out2 = dense(h2) 20 t2.watch(out2) 21 a2 = tf.math.softmax(out2) 22 t2.watch(a2) 23 # ------time stpe 3-------- 24 z3 = tf.matmul(inputs[:, 2, :], rnn.weights[0]) + tf.matmul(h2, rnn.weights[1]) + rnn.weights[2] 25 t2.watch(z3) 26 h3 = tf.math.tanh(z3) 27 t2.watch(h3) 28 out3 = dense(h3) 29 t2.watch(out3) 30 a3 = tf.math.softmax(out3) 31 t2.watch(a3) 32 33 # -------loss--------- 34 my_seqout = tf.stack([a1, a2, a3], axis=1) 35 my_loss = get_crossentropy(y_pred=my_seqout, y_true=y_true) 36 # -------L(t)--------- 37 loss_1 = get_crossentropy(y_pred=my_seqout[:, 0, :], y_true=y_true[:, 0, :]) 38 loss_2 = get_crossentropy(y_pred=my_seqout[:, 1, :], y_true=y_true[:, 1, :]) 39 loss_3 = get_crossentropy(y_pred=my_seqout[:, 2, :], y_true=y_true[:, 2, :])
為方便結合公式理解，下面是vanilla RNN前傳的核心公式：
$\boldsymbol{h}^{(t)} = \sigma\left( \boldsymbol{z}^{(t)} \right) = \sigma\left( {\boldsymbol{U}\boldsymbol{x}^{(t)} + \boldsymbol{W}\boldsymbol{h}^{(t - 1)} + \boldsymbol{b}} \right)$
$\boldsymbol{o}^{(t)} = \boldsymbol{V}\boldsymbol{h}^{(t - 1)} + \boldsymbol{c}$
${\hat{\boldsymbol{y}}}^{(t)} = \sigma\left( \boldsymbol{o}^{(t)} \right)$
我們先看看先前我們定義的rnn layer的weights：
rnn.weights Out[12]: [<tf.Variable 'rnn/simple_rnn_cell/kernel:0' shape=(4, 2) dtype=float32, numpy= array([[ 0.8315637 , -0.6328325 ], [-0.6249914 , -0.02345729], [ 0.3253579 , -0.66497874], [ 0.5830157 , -0.03220487]], dtype=float32)>, <tf.Variable 'rnn/simple_rnn_cell/recurrent_kernel:0' shape=(2, 2) dtype=float32, numpy= array([[-0.26513684, 0.96421075], [ 0.96421075, 0.2651369 ]], dtype=float32)>, <tf.Variable 'rnn/simple_rnn_cell/bias:0' shape=(2,) dtype=float32, numpy=array([0., 0.], dtype=float32)>]
可以看到有3類變數，shape=(4, 2)的kernel 對應著上述公式的$\boldsymbol{U}$；shape=(2, 2)的recurrent_kernel對應著$\boldsymbol{W}$；shape=(2, )的bias對應著$\boldsymbol{b}$
所以第6+8行程式碼就是在實現h狀態遞推公式（因為bias等於0，所以沒在程式碼中體現出來）。
接下來的程式碼10+12則是實現了上面公式組的後兩行。其中$\boldsymbol{V}$和$\boldsymbol{c}$分別對應下面的kernel和bias：
dense.weights Out[14]: [<tf.Variable 'dense/kernel:0' shape=(2, 3) dtype=float32, numpy= array([[ 0.16630626, -0.03400385, -0.8589585 ], [-1.0909375 , -0.02843404, 0.2594756 ]], dtype=float32)>, <tf.Variable 'dense/bias:0' shape=(3,) dtype=float32, numpy=array([0., 0., 0.], dtype=float32)>]
為突出重點這裡不會去驗證前面定義好了的dense layer是否真的按照FNN的前向傳播公式在做，這部分驗證可參考FNN（DNN）前向和反向傳播過程的程式碼驗證。
至此，我們完成了h狀態在時間步上的遞進，也完成了當前time step的輸出，剩下的程式碼就是繼續迴圈再進行多兩次。注意在計算time step2的h狀態時，h的遞推公式用到的就是上一個time step的h狀態h1而不是init_state[0]了
最後我們看看3個time step上的前向輸出跟tensorflow給出rnn layer的輸出output_seq 是否一致：
output_seq Out[18]: <tf.Tensor: shape=(1, 3, 3), dtype=float32, numpy= array([[[0.43937117, 0.37462497, 0.18600382], [0.32436833, 0.3793582 , 0.29627347], [0.6205071 , 0.2681237 , 0.11136916]]], dtype=float32)> my_seqout Out[19]: <tf.Tensor: shape=(1, 3, 3), dtype=float32, numpy= array([[[0.43937117, 0.37462497, 0.18600382], [0.32436833, 0.3793582 , 0.29627347], [0.6205071 , 0.2681237 , 0.11136916]]], dtype=float32)>
看來並沒有問題。
--------反向傳播驗證---------
先看$\frac{\partial L}{\partial\boldsymbol{V}}$，按照公式它應當滿足：
$\frac{\partial L}{\partial\boldsymbol{V}} = {\sum\limits_{t = 1}^{T}\frac{\partial L^{(t)}}{\partial\boldsymbol{V}}} = {\sum\limits_{t = 1}^{T}\left( {{\hat{\boldsymbol{y}}}^{(t)} - \boldsymbol{y}^{(t)}} \right)}\left( \boldsymbol{h}^{(t)} \right)^{T}$
下面是對比結果，其中dl_dV = t2.gradient(my_loss, dense.kernel) 表示這是通過tensorflow自動微分工具求得的$\frac{\partial L}{\partial\boldsymbol{V}}$，而帶my_字首的則是根據vanilla RNN的前向傳播和反向梯度推導 中的公式（就是上面這個式子）手動實現的結果。後續的符號同樣沿用這樣的命名規則。
（.transpose()的作用和意義見FNN（DNN）前向和反向傳播過程的程式碼驗證 給出的解釋，這裡不再贅述）
dl_dV = t2.gradient(my_loss, dense.kernel) my_dl_dV = tf.matmul(tf.transpose(h1), (a1 - y_true[:, 0, :])) + \ tf.matmul(tf.transpose(h2), (a2 - y_true[:, 1, :])) + \ tf.matmul(tf.transpose(h3), (a3 - y_true[:, 2, :])) dl_dV Out[20]: <tf.Tensor: shape=(2, 3), dtype=float32, numpy= array([[ 0.26737562, -0.01948635, -0.2478894 ], [-0.04734133, -0.02696498, 0.07430632]], dtype=float32)> my_dl_dV Out[21]: <tf.Tensor: shape=(2, 3), dtype=float32, numpy= array([[ 0.26737565, -0.01948633, -0.2478894 ], [-0.04734133, -0.02696498, 0.07430633]], dtype=float32)>
看來並沒有問題。
然後是$\frac{\partial L}{\partial\boldsymbol{c}}$，根據公式，它滿足：
$\frac{\partial L}{\partial\boldsymbol{c}} = {\sum\limits_{t = 1}^{T}\frac{\partial L^{(t)}}{\partial\boldsymbol{c}}} = {\sum\limits_{t = 1}^{T}{{\hat{\boldsymbol{y}}}^{(t)} - \boldsymbol{y}^{(t)}}}$
# ---------dl_dc------------ dl_dc = t2.gradient(my_loss, dense.bias) my_dl_dc = (a1 - y_true[:, 0, :]) + (a2 - y_true[:, 1, :]) + (a3 - y_true[:, 2, :]) dl_dc Out[22]: <tf.Tensor: shape=(3,), dtype=float32, numpy=array([ 0.3842466 , 0.02210681, -0.40635356], dtype=float32)> my_dl_dc Out[23]: <tf.Tensor: shape=(1, 3), dtype=float32, numpy=array([[ 0.3842466 , 0.02210684, -0.40635356]], dtype=float32)>
也沒有問題。
接下來是$\boldsymbol{\delta}^{(t)}$，求出它的目的是方便後面進一步求loss關於$\boldsymbol{U},\boldsymbol{W},\boldsymbol{b}$的導數。
根據公式，$\boldsymbol{\delta}^{(t)}$的逆推公式如下：
- $t = T$時，$\boldsymbol{\delta}^{(T)} = \boldsymbol{V}^{T}\left( {{\hat{\boldsymbol{y}}}^{(T)} - \boldsymbol{y}^{(T)}} \right)$
- $t < T$時，$\boldsymbol{\delta}^{(t)} = \boldsymbol{V}^{T}\left( {{\hat{\boldsymbol{y}}}^{(t)} - \boldsymbol{y}^{(t)}} \right) + \boldsymbol{W}^{T}diag\left( {\sigma^{'}\left( \boldsymbol{h}^{(t + 1)} \right)} \right)\boldsymbol{\delta}^{(t + 1)}$
根據上面兩條式子，我們寫出相應程式碼：
1 # ---------delta_t與delta_(t+1)------------ 2 delta_3 = t2.gradient(my_loss, h3) # = t2.gradient(loss_3, h3) 3 my_delta_3 = tf.matmul((a3 - y_true[:, 2, :]), tf.transpose(V)) 4 delta_2 = t2.gradient(my_loss, h2) # = t2.gradient(loss_2, h2) + t2.gradient(loss_3, h2) 5 delta_1 = t2.gradient(my_loss, h1) 6 # 已知delta_3，求 my_delta2 7 my_delta_2 = tf.matmul((a2 - y_true[:, 1, :]), tf.transpose(V)) + \ 8 tf.matmul(tf.matmul(delta_3, np.diag(tf.squeeze(1 - h3**2))), tf.transpose(rnn.weights[1])) 9 # 已知delta_2，求my_delta1 10 my_delta_1 = tf.matmul((a1 - y_true[:, 0, :]), tf.transpose(V)) + \ 11 tf.matmul(tf.matmul(delta_2, np.diag(tf.squeeze(1 - h2**2))), tf.transpose(rnn.weights[1]))
為提高效率，這裡直接看$\boldsymbol{\delta}^{(1)}$的對比結果，因為它的計算要用到$\boldsymbol{\delta}^{(2)}$和$\boldsymbol{\delta}^{(3)}$，如果它沒問題那麼剩下兩個應該也就沒問題。
delta_1 Out[24]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[-0.13369548, -0.13435042]], dtype=float32)> my_delta_1 Out[25]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[-0.13369548, -0.13435042]], dtype=float32)>
結果確實是一致的。
接著我們看$\frac{\partial L}{\partial\boldsymbol{W}}$，它滿足：
$\frac{\partial L}{\partial\boldsymbol{W}} = {\sum\limits_{t = 1}^{T}{diag\left( {\sigma^{'}\left( \boldsymbol{h}^{(t)} \right)} \right)\boldsymbol{\delta}^{(t)}\left( \boldsymbol{h}^{(t - 1)} \right)^{T}}}$
# ---------dl_dW----------- dl_dW = t2.gradient(my_loss, rnn.weights[1]) # tanh'(x) = (1 - tanh(x)**2) my_dl_dW = tf.matmul(tf.transpose(h2), (delta_3 * (1 - h3**2))) + \ tf.matmul(tf.transpose(h1), (delta_2 * (1 - h2**2))) + \ tf.matmul(tf.transpose(init_state[0]), (delta_1 * (1 - h1**2)))
dl_dW Out[28]: <tf.Tensor: shape=(2, 2), dtype=float32, numpy= array([[ 0.05229464, -0.2559137 ], [-0.02193429, -0.15005064]], dtype=float32)> my_dl_dW Out[29]: <tf.Tensor: shape=(2, 2), dtype=float32, numpy= array([[ 0.05229464, -0.2559137 ], [-0.02193429, -0.15005064]], dtype=float32)>
同樣沒有問題。
繼續看到$\frac{\partial L}{\partial\boldsymbol{b}}$，它應當滿足：
$\frac{\partial L}{\partial\boldsymbol{b}} = {\sum\limits_{t = 1}^{T}{diag\left( {\sigma^{'}\left( \boldsymbol{h}^{(t)} \right)} \right)\boldsymbol{\delta}^{(t)}}}$
# --------dl_db------------ dl_db = t2.gradient(my_loss, rnn.weights[2]) my_dl_db = tf.matmul(delta_3, np.diag(tf.squeeze(1 - h3**2))) + \ tf.matmul(delta_2, np.diag(tf.squeeze(1 - h2**2))) + \ tf.matmul(delta_1, np.diag(tf.squeeze(1 - h1**2)))
dl_db Out[30]: <tf.Tensor: shape=(2,), dtype=float32, numpy=array([ 0.07789478, -0.40646493], dtype=float32)> my_dl_db Out[31]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[ 0.07789478, -0.40646493]], dtype=float32)>
沒有問題+1。
最後是$\frac{\partial L}{\partial\boldsymbol{U}}$，它滿足：
$\frac{\partial L}{\partial\boldsymbol{U}} = {\sum\limits_{t = 1}^{T}{diag\left( {\sigma^{'}\left( \boldsymbol{h}^{(t)} \right)} \right)\boldsymbol{\delta}^{(t)}\left( \boldsymbol{x}^{(t)} \right)^{T}}}$
# --------dl_dU----------- dl_dU = t2.gradient(my_loss, rnn.weights[0]) # tanh'(x) = (1 - tanh(x)**2) my_dl_dU = tf.matmul(tf.transpose(inputs[:, 2, :]), (delta_3 * (1 - h3**2))) + \ tf.matmul(tf.transpose(inputs[:, 1, :]), (delta_2 * (1 - h2**2))) + \ tf.matmul(tf.transpose(inputs[:, 0, :]), (delta_1 * (1 - h1**2)))
dl_dU Out[32]: <tf.Tensor: shape=(4, 2), dtype=float32, numpy= array([[ 0.05618405, -0.24834853], [ 0.03428898, -0.24300456], [ 0.04625289, -0.23896447], [ 0.06348854, -0.27600297]], dtype=float32)> my_dl_dU Out[33]: <tf.Tensor: shape=(4, 2), dtype=float32, numpy= array([[ 0.05618405, -0.24834856], [ 0.03428898, -0.24300456], [ 0.04625289, -0.23896447], [ 0.06348854, -0.27600297]], dtype=float32)>
沒有問題+1。
至此，vanilla RNN的所有前向和反向傳播公式已驗證完。
（歡迎轉載，轉載請註明出處。歡迎留言或溝通交流： lxwalyw@gmail.com）