在《神经网络的梯度推导与代码验证》之LSTM的前向传播和反向梯度推导 中，我们学习了LSTM的前向传播和反向梯度求导，但知识仍停留在纸面。本篇章将基于深度学习框架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快速实现前向传播是什么样的：
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 = np.random.random([1, 3, 4]).astype(np.float32) # 初始化c和h状态： 第一个是h状态，第二个是c init_state = [tf.constant(np.random.random((inputs.shape[0], 2)).astype(np.float32)), tf.constant(np.random.random((inputs.shape[0], 2)).astype(np.float32))] # 定义LSTM layer lstm_cell = tf.keras.layers.LSTMCell(2) lstm = tf.keras.layers.RNN(lstm_cell, return_sequences=True, return_state=True) lstm_seq, last_h, last_c = lstm(inputs=inputs, initial_state=init_state) # fnn layer dense = tf.keras.layers.Dense(3) # 最后得到y output_seq = tf.math.softmax(dense(lstm_seq))
我们随便造一条样本(inputs, y_ture)，跟vanilla RNN的代码验证一样，输入inputs是一条步长为3的有4维特征的数据；标签数据步长也为3，每个步长上是长度为3的概率向量。
inputs.shape Out[3]: (1, 3, 4) y_true.shape Out[4]: (1, 3, 3)
至于初始状态init_state，区别于vanilla RNN，LSTM的内部状态不仅有h，还有c状态。
init_state Out[5]: [<tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[0.56804454, 0.92559665]], dtype=float32)>, <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[0.07103606, 0.0871293 ]], dtype=float32)>]
它们的dim都是2，意味着接下来定义LSTM layer的时候，units等于2，也就是第10行代码。
代码12行输出了inputs经过LSTM layer的h状态序列，并返回了最后一个time step的h和c 状态：
lstm_seq Out[6]: <tf.Tensor: shape=(1, 3, 2), dtype=float32, numpy= array([[[0.0734415 , 0.0337011 ], [0.12446897, 0.04911498], [0.08435698, 0.1005574 ]]], dtype=float32)> last_h Out[7]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[0.08435698, 0.1005574 ]], dtype=float32)> last_c Out[8]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[0.4424333 , 0.25549415]], dtype=float32)>
显然lstm_seq[0, -1, :] 就是last_h，这点在上面这段代码得到了验证。
lstm_seq经过全连接最后得到shape和y_true一样的输出序列：
output_seq Out[9]: <tf.Tensor: shape=(1, 3, 3), dtype=float32, numpy= array([[[0.33866003, 0.33178926, 0.32955068], [0.34256846, 0.3297772 , 0.3276543 ], [0.3379959 , 0.3387031 , 0.32330096]]], dtype=float32)>
通过调用上面几行代码，我们似乎能够实现上图的操作。但尚未确定代码lstm = tf.keras.layers.RNN(lstm_cell, return_sequences=True, return_state=True)是否按正真按照LSTM的前传公式进行计算的。
从下面tensorflow的LSTM源码上看其实不太好对应上前向传播公式，因为kernel和bias实际上是包含了四个门的kernel和bias的大tensor，如果想要对照着前传公式并配合源码一起消化的话，可以参考这篇博文的源码解读。
def step(cell_inputs, cell_states): """Step function that will be used by Keras RNN backend.""" h_tm1 = cell_states[0] # previous memory state c_tm1 = cell_states[1] # previous carry state z = K.dot(cell_inputs, kernel) z += K.dot(h_tm1, recurrent_kernel) z = K.bias_add(z, bias) z0, z1, z2, z3 = array_ops.split(z, 4, axis=1) i = nn.sigmoid(z0) f = nn.sigmoid(z1) c = f * c_tm1 + i * nn.tanh(z2) o = nn.sigmoid(z3) h = o * nn.tanh(c) return h, [h, c]
但我们最好还是动手实操一下，这样印象可以更加深刻，所以接下来我们手动按照LSTM的前传公式实现一遍前向传播看跟tensorflow给出的输出结果是否一致。
下面这段代码看似很长，但实际上只是反复做同样的事情而已，它实现了LSTM在时间上展开3步的前向操作，并计算了loss，所以真正的代码量你可以认为大约只有1/3这样。
这里 tf.GradientTape(persistent=True) ，t.watch()是用于后面计算变量的导数用的，不太熟悉的可参考tensorflow官方给出的关于这部分的教程（自动微分）。
1 # 提取lstm的变量 2 kernel = lstm.weights[0] 3 recurrent_kernel = lstm.weights[1] 4 bias = lstm.weights[2] 5 6 with tf.GradientTape(persistent=True) as t: 7 h0, c0 = init_state 8 # --------stpe1----------- 9 z1 = tf.matmul(inputs[:, 0, :], kernel) 10 z1 += tf.matmul(init_state[0], recurrent_kernel) 11 z1 += bias 12 13 i1 = tf.math.sigmoid(z1[:, 0:2]) 14 f1 = tf.math.sigmoid(z1[:, 2:4]) 15 c1 = f1 * init_state[1] + i1 * tf.math.tanh(z1[:, 4:6]) 16 t.watch(c1) 17 o1 = tf.math.sigmoid(z1[:, 6:8]) 18 19 h1 = o1 * tf.math.tanh(c1) 20 t.watch(h1) 21 out1 = dense(h1) 22 t.watch(out1) 23 a1 = tf.math.softmax(out1) 24 t.watch(a1) 25 # --------stpe2------------- 26 z2 = tf.matmul(inputs[:, 1, :], kernel) 27 z2 += tf.matmul(h1, recurrent_kernel) 28 z2 += bias 29 30 i2 = tf.math.sigmoid(z2[:, 0:2]) 31 f2 = tf.math.sigmoid(z2[:, 2:4]) 32 c2 = f2 * c1 + i2 * tf.math.tanh(z2[:, 4:6]) 33 t.watch(c2) 34 o2 = tf.math.sigmoid(z2[:, 6:8]) 35 36 h2 = o2 * tf.math.tanh(c2) 37 t.watch(h2) 38 out2 = dense(h2) 39 t.watch(out2) 40 a2 = tf.math.softmax(out2) 41 t.watch(a2) 42 # ---------step3--------------- 43 z3 = tf.matmul(inputs[:, 2, :], kernel) 44 z3 += tf.matmul(h2, recurrent_kernel) 45 z3 += bias 46 47 i3 = tf.math.sigmoid(z3[:, 0:2]) 48 f3 = tf.math.sigmoid(z3[:, 2:4]) 49 c3 = f3 * c2 + i3 * tf.math.tanh(z3[:, 4:6]) 50 t.watch(c3) 51 o3 = tf.math.sigmoid(z3[:, 6:8]) 52 53 h3 = o3 * tf.math.tanh(c3) 54 t.watch(c3) 55 out3 = dense(h3) 56 t.watch(out3) 57 a3 = tf.math.softmax(out3) 58 t.watch(a3) 59 # ---------loss---------------- 60 my_seqout = tf.stack([a1, a2, a3], axis=1) 61 my_loss = get_crossentropy(y_pred=my_seqout, y_true=y_true) 62 # -------L(t)--------- 63 loss_1 = get_crossentropy(y_pred=my_seqout[:, 0, :], y_true=y_true[:, 0, :]) 64 loss_2 = get_crossentropy(y_pred=my_seqout[:, 1, :], y_true=y_true[:, 1, :]) 65 loss_3 = get_crossentropy(y_pred=my_seqout[:, 2, :], y_true=y_true[:, 2, :])
为方便结合公式理解，下面是LSTM前传的核心公式：
接下来解释上述代码，首先是看看lstm.weights都有什么：
lstm.weights Out[11]: [<tf.Variable 'rnn/lstm_cell/kernel:0' shape=(4, 8) dtype=float32, numpy= array([[-0.20611948, -0.13917124, -0.4464248 , -0.5220002 , -0.07079256, -0.11765444, -0.6348312 , -0.09944093], [-0.36945656, 0.49826902, -0.38441902, 0.5155092 , 0.43278998, -0.6451189 , 0.25279564, -0.5977526 ], [-0.16671354, -0.34911466, -0.36100882, -0.441833 , 0.22657168, 0.6603723 , -0.64222896, 0.39119774], [ 0.50059 , -0.46351027, 0.21986896, 0.5533599 , 0.34416074, 0.46849674, -0.6191046 , -0.6988822 ]], dtype=float32)>, <tf.Variable 'rnn/lstm_cell/recurrent_kernel:0' shape=(2, 8) dtype=float32, numpy= array([[ 0.0145992 , -0.7155518 , -0.02687099, -0.21511525, -0.5779975 , 0.2309799 , -0.06609378, -0.22130255], [-0.06615384, -0.62732273, -0.27503067, 0.35276616, 0.4759796 , -0.19590497, -0.18337685, 0.3216232 ]], dtype=float32)>, <tf.Variable 'rnn/lstm_cell/bias:0' shape=(8,) dtype=float32, numpy=array([0., 0., 1., 1., 0., 0., 0., 0.], dtype=float32)>]
跟vanilla RNN类似，也是有3类变量。但注意这里3类变量的shape中都有一个维度等于8，因为这是4种门的kernel，recurrent_kernel和bias拼在一起的结果，这样方便进行批量矩阵乘法。
回到代码，代码9~11实现的就是上面前传公式组中的：
$\boldsymbol{W}_{f}\boldsymbol{h}^{(t - 1)} + \boldsymbol{U}_{f}\boldsymbol{x}^{(t)} + \boldsymbol{b}_{f}$
$\boldsymbol{W}_{i}\boldsymbol{h}^{(t - 1)} + \boldsymbol{U}_{i}\boldsymbol{x}^{(t)} + \boldsymbol{b}_{i}$
$\boldsymbol{W}_{a}\boldsymbol{h}^{(t - 1)} + \boldsymbol{U}_{a}\boldsymbol{x}^{(t)} + \boldsymbol{b}_{a}$
$\boldsymbol{W}_{o}\boldsymbol{h}^{(t - 1)} + \boldsymbol{U}_{o}\boldsymbol{x}^{(t - 1)} + \boldsymbol{b}_{o}$
只不过这四组运算的结果都合一块变成了z1，如果要专门提取某个门的z结果，需要对z1进行索引操作。例如i1 = tf.math.sigmoid(z1[:, 0:2])在做的就是：
$\boldsymbol{i}^{(t)} = \sigma\left( {\boldsymbol{W}_{i}\boldsymbol{h}^{(t - 1)} + \boldsymbol{U}_{i}\boldsymbol{x}^{(t)} + \boldsymbol{b}_{i}} \right)$
其他3个门的计算以此类推。
c1 = f1 * init_state[1] + i1 * tf.math.tanh(z1[:, 4:6])实现的是c状态的递推，即：
$\boldsymbol{C}^{(t)} = \boldsymbol{C}^{(t - 1)}\bigodot\boldsymbol{f}^{(t)} + \boldsymbol{a}^{(t)}\bigodot\boldsymbol{i}^{(t)}$
h1 = o1 * tf.math.tanh(c1)则实现了上述公式组的h状态的计算，即：
$\boldsymbol{h}^{(t)} = \boldsymbol{o}^{(t)}\bigodot tanh\left( \boldsymbol{C}^{(t)} \right)$
21+23行代码则实现了LSTM的输出计算，即：
${\hat{\boldsymbol{y}}}^{(t)} = \sigma\left( {\boldsymbol{V}\boldsymbol{h}^{(t)} + \boldsymbol{c}} \right)$
其中desne layer的kernel和bias分别对应$\boldsymbol{V}$和$\boldsymbol{c}$
为突出重点这里不会去验证前面定义好了的dense layer是否真的按照FNN的前向传播公式在做，这部分验证可参考FNN（DNN）前向和反向传播过程的代码验证。
至此，我们完成了c和h状态在时间步上的递进，也完成了当前time step的输出，剩下的代码就是继续循环再进行多两次。注意在计算time step2的c和h状态时，递推公式用到的就是上一个time step的c和h状态init_state了。
最后我们看看3个time step上的前向输出跟tensorflow给出rnn layer的输出output_seq 是否一致：
output_seq Out[12]: <tf.Tensor: shape=(1, 3, 3), dtype=float32, numpy= array([[[0.33866003, 0.33178926, 0.32955068], [0.34256846, 0.3297772 , 0.3276543 ], [0.3379959 , 0.3387031 , 0.32330096]]], dtype=float32)> my_seqout Out[13]: <tf.Tensor: shape=(1, 3, 3), dtype=float32, numpy= array([[[0.33866003, 0.33178926, 0.32955068], [0.34256846, 0.3297772 , 0.3276543 ], [0.3379959 , 0.3387031 , 0.32330096]]], dtype=float32)>
看来并没有问题。
--------反向传播验证---------
我们先将4个门的recurrent_kernel，即前传公式中组中的$\boldsymbol{W}_{i}$，$\boldsymbol{W}_{f}$，$\boldsymbol{W}_{a}$和$\boldsymbol{W}_{o}$分别提取出来：
w_i = recurrent_kernel[:, 0:2] w_f = recurrent_kernel[:, 2:4] w_a = recurrent_kernel[:, 4:6] w_o = recurrent_kernel[:, 6:8]
因为是BPTT，所以首先我们验证$\boldsymbol{\delta}_{h}^{(T)}$和$\boldsymbol{\delta}_{C}^{(T)}$，根据公式，它们满足：
$\boldsymbol{\delta}_{h}^{(T)} = \boldsymbol{V}^{T}\left( {{\hat{\boldsymbol{y}}}^{(T)} - \boldsymbol{y}^{(T)}} \right)$
$\boldsymbol{\delta}_{C}^{(T)} = \left( \frac{\partial\boldsymbol{h}^{(T)}}{\partial\boldsymbol{C}^{(T)}} \right)^{T}\frac{\partial L}{\partial\boldsymbol{h}^{(T)}} = \boldsymbol{\delta}_{h}^{(T)}\bigodot\boldsymbol{o}^{(T)}\bigodot{tanh}^{'}\left( \boldsymbol{C}^{(T)} \right)$
下面是对比结果，其中delta_hT = t.gradient(my_loss, h3)表示这是通过tensorflow自动微分工具求得的$\boldsymbol{\delta}_{h}^{(T)}$，而带my_前缀的则是根据LSTM的前向传播和反向梯度推导 中的公式（就是上面的式子）手动实现的结果。后续的符号同样沿用这样的命名规则。
（.transpose()的作用和意义见FNN（DNN）前向和反向传播过程的代码验证 给出的解释，这里不再赘述）
# --------to validate delta_hT--------- delta_hT = t.gradient(my_loss, h3) my_delta_hT = tf.matmul((a3 - y_true[:, 2, :]), tf.transpose(dense.kernel)) # --------to validate delta_cT--------- delta_cT = t.gradient(my_loss, c3) my_delta_cT = delta_hT * o3 * (1 - tf.math.tanh(c3)**2)
delta_hT Out[14]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[-0.07147076, 0.0525395 ]], dtype=float32)> my_delta_hT Out[15]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[-0.07147078, 0.05253952]], dtype=float32)> delta_cT Out[16]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[-0.01199877, 0.01980529]], dtype=float32)> my_delta_cT Out[17]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[-0.01199877, 0.01980529]], dtype=float32)>
验证无误，公式是正确的。
接下来是验证两个重要的递推公式：
第一个递推公式，从$\boldsymbol{\delta}_{C}^{(t + 1)}$和$\boldsymbol{\delta}_{h}^{(t + 1)}$推得$\boldsymbol{\delta}_{h}^{(t)}$：
$\boldsymbol{\delta}_{h}^{(t)} = \frac{\partial l\left( t \right)}{\partial\boldsymbol{h}^{(t)}} + \left( \frac{\partial\boldsymbol{C}^{(t + 1)}}{\partial\boldsymbol{h}^{(t)}} \right)^{T}\boldsymbol{\delta}_{C}^{(t + 1)} + \left( {\frac{\partial\boldsymbol{h}^{(t + 1)}}{\partial\boldsymbol{o}^{(t)}}\frac{\partial\boldsymbol{o}^{(t + 1)}}{\partial\boldsymbol{h}^{(t)}}} \right)^{T}\boldsymbol{\delta}_{h}^{(t + 1)}$
其中，
$\frac{\partial\boldsymbol{C}^{(t + 1)}}{\partial\boldsymbol{h}^{(t)}} = diag\left( {\boldsymbol{C}^{(t)}\bigodot\boldsymbol{f}^{({t + 1})}\bigodot\left( {1 - \boldsymbol{f}^{({t + 1})}} \right)} \right)\boldsymbol{W}_{f} + diag\left( {\boldsymbol{a}^{({t + 1})}\bigodot\boldsymbol{i}^{({t + 1})}\bigodot\left( {1 - \boldsymbol{i}^{({t + 1})}} \right)} \right)\boldsymbol{W}_{i} + diag\left( {\boldsymbol{i}^{({t + 1})}\bigodot\left( {1 - {\boldsymbol{a}^{({t + 1})}}^{2}} \right)} \right)\boldsymbol{W}_{a}$
$\frac{\partial\boldsymbol{h}^{(t + 1)}}{\partial\boldsymbol{o}^{(t)}}\frac{\partial\boldsymbol{o}^{(t + 1)}}{\partial\boldsymbol{h}^{(t)}} = diag\left( {tanh\left( \boldsymbol{C}^{({t + 1})} \right)\bigodot\boldsymbol{o}^{({t + 1})}\bigodot\left( {1 - \boldsymbol{o}^{({t + 1})}} \right)} \right)$
验证代码如下：
# -----------from delta_c2 and delta_h2 to delta_h1---------- a_2 = tf.math.tanh(z2[:, 4:6]) delta_h2 = t.gradient(my_loss, h2) delta_c2 = t.gradient(my_loss, c2) dc2_dh1 = tf.matmul(w_f, np.diag(tf.squeeze(c1 * f2 * (1 - f2)))) + \ tf.matmul(w_i, np.diag(tf.squeeze(a_2 * i2 * (1 - i2)))) + \ tf.matmul(w_a, np.diag(tf.squeeze(i2 * (1 - a_2**2)))) delta_h1 = t.gradient(my_loss, h1) my_delta_h1 = tf.matmul((a1 - y_true[:, 0, :]), tf.transpose(dense.kernel)) + \ tf.matmul(delta_c2, tf.transpose(dc2_dh1)) + \ tf.matmul(delta_h2, tf.transpose(tf.matmul(w_o, np.diag(tf.squeeze(tf.math.tanh(c2) * o2 * (1 - o2))))))
delta_h1 Out[18]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[ 0.0437731 , -0.09992676]], dtype=float32)> my_delta_h1 Out[19]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[ 0.0437731 , -0.09992676]], dtype=float32)>
没有问题，递推公式正确。
再是另外一个递推公式， 从$\boldsymbol{\delta}_{h}^{(t)}$和$\boldsymbol{\delta}_{C}^{(t + 1)}$推得$\boldsymbol{\delta}_{C}^{(t)}$：
$\boldsymbol{\delta}_{C}^{(t)} = \left( \frac{\partial\boldsymbol{h}^{(t)}}{\partial\boldsymbol{c}^{(t)}} \right)^{T}\boldsymbol{\delta}_{h}^{(t)} + \left( \frac{\partial\boldsymbol{c}^{(t + 1)}}{\partial\boldsymbol{c}^{(t)}} \right)^{T}\boldsymbol{\delta}_{C}^{(t + 1)} = \boldsymbol{o}^{(t)} \odot \left( {1 - {tanh}^{2}\left( \boldsymbol{C}^{(t)} \right)} \right)^{2}{\odot \boldsymbol{\delta}}_{h}^{(t)} + \boldsymbol{f}^{({t + 1})} \odot \boldsymbol{\delta}_{C}^{(t + 1)}$
验证代码如下：
# --------------from delta_h1 and delta_c2 to delta_c1--------- delta_c1 = t.gradient(my_loss, c1) my_delta_c1 = o1 * (1 - tf.math.tanh(c1)**2) * delta_h1 + f2 * delta_c2
delta_c1 Out[20]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[ 0.01076444, -0.01738559]], dtype=float32)> my_delta_c1 Out[21]: <tf.Tensor: shape=(1, 2), dtype=float32, numpy=array([[ 0.01076444, -0.01738559]], dtype=float32)>
没问题+1。
最后是验证$\frac{\partial L}{\partial\boldsymbol{W}_{f}}$，根据公式，它满足：
$\frac{\partial L}{\partial\boldsymbol{W}_{f}} = {\sum\limits_{t = 1}^{T}\left\lbrack {\boldsymbol{\delta}_{C}^{(t)} \odot \boldsymbol{C}^{(t - 1)} \odot \boldsymbol{f}^{(t)} \odot \left\lbrack {1 - \boldsymbol{f}^{(t)}} \right\rbrack} \right\rbrack}\left( \boldsymbol{h}^{(t - 1)} \right)^{T}$
代码验证如下：
# -------------to validate dl_dW_f--------------- dl_dW_f = t.gradient(my_loss, recurrent_kernel)[:, 2:4] my_dl_dW_f = tf.matmul(tf.transpose(h2), (delta_cT * c2 * f3 * (1 - f3))) + \ tf.matmul(tf.transpose(h1), (delta_c2 * c1 * f2 * (1 - f2))) + \ tf.matmul(tf.transpose(h0), (delta_c1 * c0 * f1 * (1 - f1)))
dl_dW_f Out[22]: <tf.Tensor: shape=(2, 2), dtype=float32, numpy= array([[-6.2376050e-05, -2.1518332e-05], [ 1.0945055e-04, -1.8339025e-04]], dtype=float32)> my_dl_dW_f Out[23]: <tf.Tensor: shape=(2, 2), dtype=float32, numpy= array([[-6.2376050e-05, -2.1518332e-05], [ 1.0945055e-04, -1.8339025e-04]], dtype=float32)>
没问题+1。
至此，LSTM的所有前向和部分反向传播公式已验证完，剩下的反向传播公式例如$\frac{\partial L}{\partial\boldsymbol{W}_{i}}$和$\frac{\partial L}{\partial\boldsymbol{W}_{o}}$等，都是类似这样的操作罢了。
如果本文对您有所帮助的话，还请点下“推荐”让它能帮到更多的人，谢谢。
（欢迎转载，转载请注明出处。欢迎留言或沟通交流： [email protected]）