《神经网络的梯度推导与代码验证》之vanilla RNN前向和反向传播的代码验证
在《神经网络的梯度推导与代码验证》之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就行
import tensorflow as tf import numpy as np np.random.seed(0)
然后是定义交叉熵损失函数：
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）
