在《神经网络的梯度推导与代码验证》之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快速实现前向传播是什么样的，代码挺短的：
y_true = np.array([[[0.3, 0.5, 0.2],
[0.2, 0.3, 0.5],
[0.5, 0.2, 0.3]]]).astype(np.float32) # --------inputs---------
inputs = np.random.random([1, 3, 4]).astype(np.float32)
init_state = [tf.constant(np.random.random((inputs.shape[0], 2)).astype(np.float32))]
# --------vanilla rnn---------
rnn = tf.keras.layers.RNN(tf.keras.layers.SimpleRNNCell(2), return_sequences=True)
h_seq = rnn(inputs=inputs, initial_state=init_state)
# --------fnn-------------
dense = tf.keras.layers.Dense(3)
output_dense = dense(h_seq)
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官方给出的关于这部分的教程（自动微分）。
with tf.GradientTape(persistent=True) as t2:
# -------rnn analysis by steps---------
# 下面是手动演算的rnn展开，和上面的whole_sequence_output是对得上的 # ------time stpe 1--------
z1 = tf.matmul(inputs[:, 0, :], rnn.weights[0]) + tf.matmul(init_state[0], rnn.weights[1]) + rnn.weights[2]
t2.watch(z1)
h1 = tf.math.tanh(z1)
t2.watch(h1)
out1 = dense(h1)
t2.watch(out1)
a1 = tf.math.softmax(out1)
t2.watch(a1)
# ------time stpe 2--------
z2 = tf.matmul(inputs[:, 1, :], rnn.weights[0]) + tf.matmul(h1, rnn.weights[1]) + rnn.weights[2]
t2.watch(z2)
h2 = tf.math.tanh(z2)
t2.watch(h2)
out2 = dense(h2)
t2.watch(out2)
a2 = tf.math.softmax(out2)
t2.watch(a2)
# ------time stpe 3--------
z3 = tf.matmul(inputs[:, 2, :], rnn.weights[0]) + tf.matmul(h2, rnn.weights[1]) + rnn.weights[2]
t2.watch(z3)
h3 = tf.math.tanh(z3)
t2.watch(h3)
out3 = dense(h3)
t2.watch(out3)
a3 = tf.math.softmax(out3)
t2.watch(a3) # -------loss---------
my_seqout = tf.stack([a1, a2, a3], axis=1)
my_loss = get_crossentropy(y_pred=my_seqout, y_true=y_true)
# -------L(t)---------
loss_1 = get_crossentropy(y_pred=my_seqout[:, 0, :], y_true=y_true[:, 0, :])
loss_2 = get_crossentropy(y_pred=my_seqout[:, 1, :], y_true=y_true[:, 1, :])
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)}$
根据上面两条式子，我们写出相应代码：
# ---------delta_t与delta_(t+1)------------
delta_3 = t2.gradient(my_loss, h3) # = t2.gradient(loss_3, h3)
my_delta_3 = tf.matmul((a3 - y_true[:, 2, :]), tf.transpose(V))
delta_2 = t2.gradient(my_loss, h2) # = t2.gradient(loss_2, h2) + t2.gradient(loss_3, h2)
delta_1 = t2.gradient(my_loss, h1)
# 已知delta_3，求 my_delta2
my_delta_2 = tf.matmul((a2 - y_true[:, 1, :]), tf.transpose(V)) + \
tf.matmul(tf.matmul(delta_3, np.diag(tf.squeeze(1 - h3**2))), tf.transpose(rnn.weights[1]))
# 已知delta_2，求my_delta1
my_delta_1 = tf.matmul((a1 - y_true[:, 0, :]), tf.transpose(V)) + \
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）
