【作者主页】Francek Chen
【专栏介绍】$\lceil$Python机器学习$\rfloor$ 机器学习是一门人工智能的分支学科，通过算法和模型让计算机从数据中学习，进行模型训练和优化，做出预测、分类和决策支持。Python成为机器学习的首选语言，依赖于强大的开源库如Scikit-learn、TensorFlow和PyTorch。本专栏介绍机器学习的相关算法以及基于Python的算法实现。
【GitCode】专栏资源保存在我的GitCode仓库：https://gitcode.com/Morse_Chen/Python_machine_learning。

---

  在前面两篇文章中，我们分别介绍了神经网络的基础概念和最简单的MLP，以及适用于图像处理的CNN。从中我们可以意识到，不同结构的神经网络具有不同的特点，在不同任务上具有自己的优势。例如MLP复杂度低、训练简单、适用范围广，适合解决普通任务或作为大型网络的小模块；CNN可以捕捉到输入中不同尺度的关联信息，适合从图像中提取特征。而对于具有序列特征的数据，例如一年内随时间变化的温度、一篇文章中的文字等，它们具有明显的前后关联。然而这些关联的数据在序列中出现的位置可能间隔非常远，例如文章在开头和结尾描写了同一个事物，如果用CNN来提取这些关联的话，其卷积核的大小需要和序列的长度相匹配。当数据序列较长时，这种做法会大大增加网络复杂度和训练难度。因此，我们需要引入一种新的网络结构，使其能够充分利用数据的序列性质，从前到后分析数据、提取关联。这就是本文要介绍的循环神经网络（recurrent neural networks，RNN）。

一、循环神经网络的基本原理

  我们先从最简单的模型开始考虑。对于不存在序列关系的数据，我们采用一个两层的MLP来拟合它，如图1(a)所示，输入样本为$\mathbf{x}$，经过第一个权重为${\mathbf{W}}_i$和${\mathbf{b}}_i$的隐层得到中间向量 $\mathbf{h}={\mathbf{f}}_h({\mathbf{W}}_i\mathbf{x}+{\mathbf{b}}_i)$，再经过权重为${\mathbf{W}}_o$和${\mathbf{b}}_o$的隐层得到输出 $\mathbf{y}=f_o({\mathbf{W}}_o\mathbf{h}+{\mathbf{b}}_o)$，其中$f_h$和$f_o$为激活函数。这是一个标准的MLP的预测流程。

图1 从MLP到RNN

  假设数据集中的数据分别是在时刻1和时刻2采集到的，并且我们知道时刻2的结果与时刻1有关。这时，由于两个时刻的数据产生了依赖关系，如果我们用相同的模型权重来进行预测而忽略其关联，预测的准确度就会降低。为了利用上额外的关联信息，我们将MLP的结构拓展一下，如图1(b)所示，第二个MLP的中间向量与一般的MLP不同。在计算时刻2的中间向量${\mathbf{h}}_2$时，我们将时刻1的中间向量${\mathbf{h}}_1$也纳入进来，得到 ${\mathbf{h}}_2=f_h({\mathbf{W}}_h{\mathbf{h}}_1+{\mathbf{W}}_i{\mathbf{x}}_2+{\mathbf{b}}_i)$，再将${\mathbf{h}}_2$传给第二个隐层，计算出输出 $y_2=f_o({\mathbf{W}}_o{\mathbf{h}}_2+{\mathbf{b}}_o)$。这样，我们就在时刻2的预测中用到了时刻1的信息。如果将这种思想进一步扩展，如图1(c)所示，我们可以将MLP沿着序列不断扩展下去，中间的每个MLP都将上一时刻的中间向量${\mathbf{h}}_{t-1}$与当前的输入${\mathbf{x}}_t$组合得到中间向量，再进行后续处理。同时，由于序列中每一位置之间又存在对称性，为了减小网络的复杂度，每一MLP前后的权重与中间组合的权重可以共用，不随序列位置变化。因此，这样重复的网络结构可以用图2中的循环来表示，称为循环神经网络。

图2 RNN的循环表示

  RNN的输入与输出并不一定要像上面展示的一样，在每一时刻都有一个输入样本和一个预测输出。根据任务的不同，RNN的输入输出对应可以有多种形式。图3展示了一些不同对应形式的RNN结构，从左到右依次是一对多、多对一、同步多对多和异步多对多，它们都有合适的任务场景。例如，如果我们要根据一个关键词生成一句话，以词语作为最小单元，那么RNN的输入只有一个，而生成的句子需要有连贯的含义和语义，因此可以利用RNN在每一时刻输出一个词，从前到后连成完整的句子。这样的任务就更适合采用一对多的结构。再比如，常见的时间序列预测任务需要我们根据一段时间中收集的数据，预测接下来一定时间内数据的情况。这时，我们就可以用异步多对多的结构，先分析样本的规律和特征，再生成紧接着样本所在时间之后的结果。

图3 适用与不同任务的RNN结构

  当我们训练RNN时，由于每一时刻的中间向量都会组合上一时刻的中间向量，如果把时刻$t$的中间向量全部展开，就得到
$$\begin{array}{rl}{\mathbf{h}}_t& =f_h({\mathbf{W}}_h{\mathbf{h}}_{t-1}+{\mathbf{W}}_i{\mathbf{x}}_t+{\mathbf{b}}_i)\\ & =f_h({\mathbf{W}}_h{f}_h({\mathbf{W}}_h{\mathbf{h}}_{t-2}+{\mathbf{W}}_i{\mathbf{x}}_{t-1}+{\mathbf{b}}_i)+{\mathbf{W}}_i{\mathbf{x}}_t+{\mathbf{b}}_i)\\ & =\cdots \\ & =f_h({\mathbf{W}}_h{f}_h(\cdots {\mathbf{W}}_h{f}_h({\mathbf{W}}_h({\mathbf{W}}_i{\mathbf{x}}_1+{\mathbf{b}}_i)+{\mathbf{W}}_i{\mathbf{x}}_2+{\mathbf{b}}_i)\cdots )+{\mathbf{W}}_i{\mathbf{x}}_t+{\mathbf{b}}_i)\end{array}$$

  如果在时刻$t$存在输出，我们可计算时刻$t$的损失函数，并使用梯度回传方法优化参数。然而，随着反向传播的步数增加，RNN有可能会出现梯度消失或梯度爆炸的现象。为了详细解释这一现象，我们考虑时刻$t$的损失${\mathcal{L}}_t$关于参数${\mathbf{W}}_i$的导数。根据求导的链式法则，我们可以计算如下：
$$\begin{array}{rl}\frac{\partial {\mathcal{L}}_t}{\partial {\mathbf{W}}_i}& =\frac{\partial {\mathcal{L}}_t}{\partial y_t}\frac{\partial y_t}{\partial {\mathbf{W}}_i}\\ & =\frac{\partial {\mathcal{L}}_t}{\partial y_t}\frac{\partial y_t}{\partial {\mathbf{h}}_t}\frac{\mathrm{d}{\mathbf{h}}_t}{\mathrm{d}{\mathbf{W}}_i}\\ & =\frac{\partial {\mathcal{L}}_t}{\partial y_t}\frac{\partial y_t}{\partial {\mathbf{h}}_t}\left(\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{W}}_i}+\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{h}}_{t-1}}\frac{\mathrm{d}{\mathbf{h}}_{t-1}}{\mathrm{d}{\mathbf{W}}_i}\right)\\ & =\cdots \\ & =\frac{\partial {\mathcal{L}}_t}{\partial y_t}\frac{\partial y_t}{\partial {\mathbf{h}}_t}\left(\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{W}}_i}+\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{h}}_{t-1}}\frac{\partial {\mathbf{h}}_{t-1}}{\partial {\mathbf{W}}_i}+\cdots +\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{h}}_{t-1}}\frac{\partial {\mathbf{h}}_{t-1}}{\partial {\mathbf{h}}_{t-2}}\cdots \frac{\partial {\mathbf{h}}_2}{\partial {\mathbf{h}}_1}\frac{\partial {\mathbf{h}}_1}{\partial {\mathbf{W}}_i}\right)\\ & =\frac{\partial {\mathcal{L}}_t}{\partial y_t}\frac{\partial y_t}{\partial {\mathbf{h}}_t}\left(\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{W}}_i}+\sum _{j=1}^{t-1}\left(\prod _{k=j+1}^t\frac{\partial {\mathbf{h}}_k}{\partial {\mathbf{h}}_{k-1}}\right)\frac{\partial {\mathbf{h}}_j}{\partial {\mathbf{W}}_i}\right)\end{array}$$ 将 $\begin{array}{r}\frac{\partial {\mathbf{h}}_k}{\partial {\mathbf{h}}_{k-1}}=f_h^{\prime }{\mathbf{W}}_h\end{array}$ 和 $\begin{array}{r}\frac{\partial {\mathbf{h}}_k}{\partial {\mathbf{W}}_i}={\mathbf{x}}_k\end{array}$ 代入，就得到$$\frac{\partial {\mathcal{L}}_t}{\partial {\mathbf{W}}_i}=\frac{\partial {\mathcal{L}}_t}{\partial y_t}\frac{\partial y_t}{\partial {\mathbf{h}}_t}\left({\mathbf{x}}_t+\sum _{j=1}^{t-1}\left(\prod _{k=j+1}^t{f}_h^{\prime }{\mathbf{W}}_h\right){\mathbf{x}}_j\right)$$

  观察上式可以发现，梯度中会出现一些$f_h^{\prime }{\mathbf{W}}_h$的连乘项。如果 $f_h^{\prime }{\mathbf{W}}_h<1$，当时刻$t$与时刻$j$距离较远时，该连乘的值就会趋近于0，因此由时刻$t$的损失函数计算出的梯度在回传时会逐渐消失；反之，如果 $f_h^{\prime }{\mathbf{W}}_h>1$，该连乘会趋于无穷大，梯度在回传时会出现发散的现象。我们将这两种情况分别称为梯度消失和梯度爆炸。无论出现哪种情况，网络的参数都无法正常更新，模型的性能也会大打折扣。当出现梯度消失时，时刻$t$的梯度只能影响时刻$t$之前的少数几步，而无法影响到较远的位置。换句话说，距离时刻$t$较远的信息已经丢失，模型很难捕捉到序列中的长期关联。而当出现梯度爆炸时，网络的梯度会迅速发散，出现数值溢出等错误。

  为了防止上述现象发生，最简单的做法是对梯度进行裁剪，为梯度设置上限和下限，当梯度过大或过小时，直接用上下限来代替梯度的值。但是，这种做法在复杂情况下仍然会导致信息丢失，通常只作为一种辅助手段。我们还可以选用合适的激活函数$f_h$并调整网络参数${\mathbf{W}}_h$初始化的值，使得乘积$f_h^{\prime }{\mathbf{W}}_h$始终稳定在1附近。但是，随着网络参数不断更新，${\mathbf{W}}_h$总会变化，要始终控制它们的乘积比较困难。因此，我们可以将网络中关联起相邻两步的$f_h$和${\mathbf{W}}_h$扩展成一个小的网络，通过设计其结构来达到稳定梯度的目的。

二、门控循环单元

  本节，我们就来介绍一种较为简单的设计——门控循环单元（gated recurrent unit，GRU）。为了解决梯度消失与梯度爆炸的问题，GRU在普通RNN的设计上进行改进，通过门控单元来调整${\mathbf{h}}_t$和${\mathbf{h}}_{t-1}$的关系。我们不妨将输入${\mathbf{x}}_t$理解为外部输入的信息，${\mathbf{h}}_t$理解为网络记住的信息，它从时刻1的${\mathbf{h}}_1$开始向后传递。然而，由于模型本身复杂度的限制，模型并不需要、也无法将所有时刻的信息都保留下来。因此，在由上一时刻的信息${\mathbf{h}}_{t-1}$计算${\mathbf{h}}_t$时，必须有选择地进行遗忘。同时，在时刻$t$有新的信息${\mathbf{x}}_t$输入进网络，我们需要在过去的信息${\mathbf{h}}_{t-1}$与新信息${\mathbf{x}}_t$之间做到平衡。

  图4展示了GRU单元的内部结构，GRU设置的门控单元共有两个，分别称为更新门和重置门。每个门控单元输出一个数值或向量，由上一时刻的信息${\mathbf{h}}_{t-1}$和当前时刻的输入${\mathbf{x}}_t$组合计算得到
$$\begin{array}{r}{\mathbf{z}}_t=\sigma ({\mathbf{W}}_z{\mathbf{x}}_t+{\mathbf{U}}_z{\mathbf{h}}_{t-1}+{\mathbf{b}}_z)\\ {\mathbf{r}}_t=\sigma ({\mathbf{W}}_r{\mathbf{x}}_t+{\mathbf{U}}_r{\mathbf{h}}_{t-1}+{\mathbf{b}}_r)\end{array}$$ 其中，${\mathbf{z}}_t$是更新单元，${\mathbf{r}}_t$是重置单元，${\mathbf{W}}_z$、${\mathbf{W}}_r$、${\mathbf{U}}_z$、${\mathbf{U}}_r$、${\mathbf{b}}_z$和${\mathbf{b}}_r$都是网络的参数，$\sigma$是逻辑斯谛函数，从而门控单元的值都在$(0,1)$区间内。

图4 GRU结构示意

  虽然这两个单元的计算方式完全相同，但是接下来它们会发挥不同的作用。利用重置单元${\mathbf{r}}_t$，我们对过去的信息${\mathbf{h}}_{t-1}$进行选择性遗忘：$${\mathbf{h}}_{t-1}^{\prime }={\mathbf{r}}_t\odot {\mathbf{h}}_{t-1}$$ 其中，$\odot$称为阿达马积（Hadamard product），表示向量或矩阵的逐元素相乘。例如，形状均为 $m\times n$ 的矩阵$\mathbf{A}$和$\mathbf{B}$的阿达玛积为
$$\mathbf{A}\odot \mathbf{B}=\left(\begin{array}{cccc}{a}_{11}{b}_{11}& a_{12}{b}_{12}& \cdots & a_{1n}{b}_{1n}\\ a_{21}{b}_{21}& a_{22}{b}_{22}& \cdots & a_{2n}{b}_{2n}\\ ⋮& ⋮& & ⋮\\ a_{m1}{b}_{m1}& a_{m2}{b}_{m2}& \cdots & a_{mn}{b}_{mn}\end{array}\right)$$ 当${\mathbf{r}}_t$某一维度的值接近0时，网络就更倾向于遗忘${\mathbf{h}}_{t-1}$的相应维度；反之，当${\mathbf{r}}_t$某一维度的值接近1时，网络更倾向于保留${\mathbf{h}}_{t-1}$的相应维度。之后，我们再将重置过的${\mathbf{h}}_{t-1}^{\prime }$与${\mathbf{x}}_t$组合，得到${\hat{\mathbf{h}}}_t$：$${\hat{\mathbf{h}}}_t=\tanh ({\mathbf{W}}_h{\mathbf{x}}_t+{\mathbf{U}}_h{\mathbf{h}}_{t-1}^{\prime }+{\mathbf{b}}_h)$$ 这里得到的${\hat{\mathbf{h}}}_t$混合了当前的${\mathbf{x}}_t$与部分过去的信息${\mathbf{h}}_{t-1}^{\prime }$，并由$\tanh$函数映射到了$(-1,1)$范围内。观察上式与普通RNN的更新方式 ${\mathbf{h}}_t=f_h({\mathbf{W}}_i{\mathbf{x}}_i+{\mathbf{W}}_h{\mathbf{h}}_{t-1}+{\mathbf{b}}_t)$，可以看出，普通的RNN相当于令重置单元${\mathbf{r}}_t$的所有维度都为1，从而保留了所有过去的信息；而 ${\mathbf{r}}_t=0$ 会消除所有过去的信息，使得RNN退化为与过去无关的单个MLP。可以通过这样的对比体会重置单元的意义。

  最后，我们要决定${\mathbf{h}}_t$是要更倾向于旧的信息${\mathbf{h}}_{t-1}$，还是旧信息与新输入${\mathbf{x}}_t$的混合${\hat{\mathbf{h}}}_t$。利用更新单元${\mathbf{z}}_t$，我们令$${\mathbf{h}}_t={\mathbf{z}}_t\odot {\mathbf{h}}_{t-1}+(1-{\mathbf{z}}_t)\odot {\hat{\mathbf{h}}}_t$$

  在上式中，如果更新单元${\mathbf{z}}_t$接近$1$，我们将保留更多的旧信息${\mathbf{h}}_{t-1}$，而忽略${\mathbf{x}}_t$的影响；反之，如果${\mathbf{z}}_t$接近$0$，我们将让旧信息与新信息混合，保留${\hat{\mathbf{h}}}_{t-1}$。需要注意，重置单元和更新单元的作用并不相同，两者不能合为一个单元。简单来说，重置单元控制旧信息保留的比例，而更新单元同时控制旧信息和新输入的比例。虽然理论上我们可以用类似 ${\mathbf{h}}_t={\mathbf{z}}_t\odot f_h({\mathbf{U}}_h{\mathbf{h}}_{t-1}+{\mathbf{b}}_h)+(1-{\mathbf{z}}_t)\odot f_x({\mathbf{W}}_x{\mathbf{x}}_t+{\mathbf{b}}_x)$ 这样的式子，仅用一个更新单元来计算${\mathbf{h}}_t$，但是其灵活性将大打折扣。

  为什么GRU的设计可以缓解梯度爆炸与梯度消失问题呢？上文我们已经提到，导致梯度问题的最大因素是$\begin{array}{r}\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{h}}_{t-1}}\end{array}$的连乘。在GRU中，我们可以通过调整门控单元${\mathbf{r}}_t$与${\mathbf{z}}_t$的值，使该梯度始终保持稳定。以文本分析为例，假如某一事物在一段话的开头和结尾出现，为了让模型保留它们之间的关联，我们只需要将重置单元${\mathbf{r}}_t$的值减小、更新单元${\mathbf{z}}_t$的值增大，就可以使网络在间隔很多时间步之后，仍然保留最初的记忆信息。最极端的情况下，如果令 ${\mathbf{z}}_2,\cdots ,{\mathbf{z}}_{t-1}=1$，那么从时刻2到时刻 $t-1$ 的所有输入都将被忽略，可以直接得到 ${\mathbf{h}}_t={\mathbf{h}}_1$。这样，梯度的连乘为$$\prod _{k=2}^t\frac{\partial {\mathbf{h}}_k}{\partial {\mathbf{h}}_{k-1}}=\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{h}}_1}=\mathbf{I}$$

  虽然门控单元的值也是由网络训练得到的，但是门控单元的引入使得GRU可以自我调节梯度。也就是说，如果${\mathbf{h}}_1$非常重要，那么门控单元会让${\mathbf{h}}_1$保留下来，其梯度较大；如果${\mathbf{h}}_1$重要性不高，随着时间推移被遗忘，那么其梯度即使消失也不会产生什么问题。因此，GRU几乎不会发生普通RNN的梯度爆炸或梯度消失现象。

三、动手实现GRU

  本节我们使用PyTorch库中的工具来实现GRU模型，完成简单的时间序列预测任务。时间序列预测任务是指根据一段连续时间内采集的数据、分析其变化规律、预测接下来数据走向的任务。如果当前数据与历史数据存在依赖关系，或者有随时间有一定的规律性，该任务就很适合用RNN求解。本节中，我们生成了一条经过一定处理的正弦曲线作为数据集，存储在sindata_1000.csv中。该曲线包含1000个数据点。其中前800个点作为训练集，后200个点作为测试集。由于本任务是时序预测任务，我们在划分训练集和测试集时无须将其打乱。我们首先导入必要的库和数据集，并将数据集的图像绘制出来。

import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
import torch
import torch.nn as nn# 导入数据集
data = np.loadtxt('sindata_1000.csv', delimiter=',')
num_data = len(data)
split = int(0.8 * num_data)
print(f'数据集大小：{num_data}')
# 数据集可视化
plt.figure()
plt.scatter(np.arange(split), data[:split], color='blue', s=10, label='training set')
plt.scatter(np.arange(split, num_data), data[split:], color='none', edgecolor='orange', s=10, label='test set')
plt.xlabel('X axis')
plt.ylabel('Y axis')
plt.legend()
plt.show()
# 分割数据集
train_data = np.array(data[:split])
test_data = np.array(data[split:])

  在训练RNN模型时，虽然我们可以把每个时间步$t$单独输入，得到模型的预测值${\hat{y}}_t$，但是这样无法体现出数据的序列相关性质。因此，我们通常会把一段时间序列 ${\mathbf{x}}_t,\cdots ,{\mathbf{x}}_{t+k}$整体作为输入，PyTorch中的GRU模块输出这段序列对应的中间变量 ${\mathbf{h}}_t,\cdots ,{\mathbf{h}}_{t+k}$。下面的实现中，我们每次输入 ${\mathbf{x}}_t,\cdots ,{\mathbf{x}}_{t+k}$ 的时间序列，预测输入向后错一步 ${\mathbf{x}}_{t+1},\cdots ,{\mathbf{x}}_{t+k+1}$ 的数据。参照图4的结构可以发现，GRU模型只输出中间变量。如果要得到我们最后的输出，还需要将这些中间变量经过自定义的其他网络。这一点和CNN里卷积层负责提取特征、MLP负责根据特征完成特定任务的做法非常相似。因此，我们在GRU之后拼接一个全连接层，通过中间变量序列 ${\mathbf{h}}_{t+1},\cdots ,{\mathbf{h}}_{t+k+1}$ 来预测未来的数据分布。

# 输入序列长度
seq_len = 20
# 处理训练数据，把切分序列后多余的部分去掉
train_num = len(train_data) // (seq_len + 1) * (seq_len + 1)
train_data = np.array(train_data[:train_num]).reshape(-1, seq_len + 1, 1)
np.random.seed(0)
torch.manual_seed(0)x_train = train_data[:, :seq_len] # 形状为(num_data, seq_len, input_size)
y_train = train_data[:, 1: seq_len + 1]
print(f'训练序列数：{len(x_train)}')# 转为PyTorch张量
x_train = torch.from_numpy(x_train).to(torch.float32)
y_train = torch.from_numpy(y_train).to(torch.float32)
x_test = torch.from_numpy(test_data[:-1]).to(torch.float32)
y_test = torch.from_numpy(test_data[1:]).to(torch.float32)

  考虑到GRU的模型结构较为复杂，我们直接使用在PyTorch库中封装好的GRU模型。我们只需要为该模型提供两个参数，第一个参数input_size表示输入$\mathbf{x}$的维度，第二个参数hidden_size表示中间向量$\mathbf{h}$的维度，其余参数我们保持默认值。在前向传播时，GRU接受序列$\mathbf{x}$和初始的中间变量$\mathbf{h}$。如果最开始我们不知道中间变量的值，GRU会自动将其初始化为全零。前向传播的输出是out和hidden，前者是整个时间序列上中间变量的值，而后者只包含是最后一步。out[-1]和hidden在GRU内部的层数不同时会有区别，但本节只使用单层网络，因此不详细展开。感兴趣的可以参考PyTorch的官方文档。我们将out作为最后全连接层的输入，得到预测值，再把预测值和hidden返回。hidden将作为下一次前向传播的初始中间变量。

class GRU(nn.Module):# 包含PyTorch的GRU和拼接的MLPdef __init__(self, input_size, output_size, hidden_size):super().__init__()# GRU模块self.gru = nn.GRU(input_size=input_size, hidden_size=hidden_size) # 将中间变量映射到预测输出的MLPself.linear = nn.Linear(hidden_size, output_size)def forward(self, x, hidden):# 前向传播# x的维度为(batch_size, seq_len, input_size)# GRU模块接受的输入为(seq_len, batch_size, input_size)# 因此需要对x进行变换# transpose函数可以交换x的坐标轴# out的维度是(seq_len, batch_size, hidden_size)out, hidden = self.gru(torch.transpose(x, 0, 1), hidden) # 取序列最后的中间变量输入给全连接层out = self.linear(out.view(-1, hidden_size))return out, hidden

  接下来，我们设置超参数并实例化GRU。在训练之前，我们还要强调时序模型在测试时与普通模型的区别。GRU在测试时，我们将输入的时间序列长度降为1，即只输入${\mathbf{x}}_t$，让GRU预测 $t+1$ 时刻的值。之后，不像普通的任务那样把所有测试数据都给模型，而是让GRU将自己预测的${\hat{\mathbf{x}}}_{t+1}$作为输入，再预测 $t+2$ 时刻的值，循环往复。这样的测试方式对模型在时序上的建模能力有相当高的要求，否则就会很快因为预测值的误差累积，和真实值偏差很大。

# 超参数
input_size = 1 # 输入维度
output_size = 1 # 输出维度
hidden_size = 16 # 中间变量维度
learning_rate = 5e-4# 初始化网络
gru = GRU(input_size, output_size, hidden_size)
gru_optim = torch.optim.Adam(gru.parameters(), lr=learning_rate)# GRU测试函数，x和hidden分别是初始的输入和中间变量
def test_gru(gru, x, hidden, pred_steps):pred = []inp = x.view(-1, input_size)for i in range(pred_steps):gru_pred, hidden = gru(inp, hidden)pred.append(gru_pred.detach())inp = gru_predreturn torch.concat(pred).reshape(-1)

  作为对比，我们用相同的数据同步训练一个3层的MLP模型。该MLP将同样将 ${\mathbf{x}}_t,\cdots ,{\mathbf{x}}_{t+k}$ 的数据拼接在一起作为输入，此时$k$被理解为输入的批量大小，并输出 ${\mathbf{x}}_{t+1},\cdots ,{\mathbf{x}}_{t+k+1}$ 的预测值，与GRU保持一致。在测试时，MLP同样只接受测试集第一个时间步的数据，以和GRU相同的方式进行自循环预测。

# MLP的超参数
hidden_1 = 32
hidden_2 = 16
mlp = nn.Sequential(nn.Linear(input_size, hidden_1),nn.ReLU(),nn.Linear(hidden_1, hidden_2),nn.ReLU(),nn.Linear(hidden_2, output_size)
)
mlp_optim = torch.optim.Adam(mlp.parameters(), lr=learning_rate)# MLP测试函数，相比于GRU少了中间变量
def test_mlp(mlp, x, pred_steps):pred = []inp = x.view(-1, input_size)for i in range(pred_steps):mlp_pred = mlp(inp)pred.append(mlp_pred.detach())inp = mlp_predreturn torch.concat(pred).reshape(-1)

  我们用完全相同的数据训练GRU和MLP。由于已经有了序列长度，我们不再设置SGD的批量大小，直接将每个训练样本单独输入模型进行优化。

max_epoch = 150
criterion = nn.functional.mse_loss
hidden = None # GRU的中间变量# 训练损失
gru_losses = []
mlp_losses = []
gru_test_losses = []
mlp_test_losses = []
# 开始训练
with tqdm(range(max_epoch)) as pbar:for epoch in pbar:st = 0gru_loss = 0.0mlp_loss = 0.0# 随机梯度下降for X, y in zip(x_train, y_train):# 更新GRU模型# 我们不需要通过梯度回传更新中间变量# 因此将其从有梯度的部分分离出来if hidden is not None:hidden.detach_()gru_pred, hidden = gru(X[None, ...], hidden)gru_train_loss = criterion(gru_pred.view(y.shape), y)gru_optim.zero_grad()gru_train_loss.backward()gru_optim.step()gru_loss += gru_train_loss.item()# 更新MLP模型# 需要对输入的维度进行调整，变成(seq_len, input_size)的形式mlp_pred = mlp(X.view(-1, input_size))mlp_train_loss = criterion(mlp_pred.view(y.shape), y)mlp_optim.zero_grad()mlp_train_loss.backward()mlp_optim.step()mlp_loss += mlp_train_loss.item()gru_loss /= len(x_train)mlp_loss /= len(x_train)gru_losses.append(gru_loss)mlp_losses.append(mlp_loss)# 训练和测试时的中间变量序列长度不同，训练时为seq_len，测试时为1gru_pred = test_gru(gru, x_test[0], hidden[:, -1], len(y_test))mlp_pred = test_mlp(mlp, x_test[0], len(y_test))gru_test_loss = criterion(gru_pred, y_test).item()mlp_test_loss = criterion(mlp_pred, y_test).item()gru_test_losses.append(gru_test_loss)mlp_test_losses.append(mlp_test_loss)pbar.set_postfix({'Epoch': epoch,'GRU loss': f'{gru_loss:.4f}','MLP loss': f'{mlp_loss:.4f}','GRU test loss': f'{gru_test_loss:.4f}','MLP test loss': f'{mlp_test_loss:.4f}'})

  最后，我们在测试集上对比GRU和MLP模型的效果并绘制出来。图中包含了原始数据的训练集和测试集的曲线，可以看出，GRU的预测基本符合测试集的变化规律，而MLP很快就因为缺乏足够的时序信息与测试集偏离。

# 最终测试结果
gru_preds = test_gru(gru, x_test[0], hidden[:, -1], len(y_test)).numpy()
mlp_preds = test_mlp(mlp, x_test[0], len(y_test)).numpy()plt.figure(figsize=(13, 5))# 绘制训练曲线
plt.subplot(121)
x_plot = np.arange(len(gru_losses)) + 1
plt.plot(x_plot, gru_losses, color='blue', label='GRU training loss')
plt.plot(x_plot, mlp_losses, color='red', ls='-.', label='MLP training loss')
plt.plot(x_plot, gru_test_losses, color='blue', ls='--', label='GRU test loss')
plt.plot(x_plot, mlp_test_losses, color='red', ls=':', label='MLP test loss')
plt.xlabel('Training step')
plt.ylabel('Loss')
plt.legend(loc='lower left')# 绘制真实数据与模型预测值的图像
plt.subplot(122)
plt.scatter(np.arange(split), data[:split], color='blue', s=10, label='training set')
plt.scatter(np.arange(split, num_data), data[split:], color='none', edgecolor='orange', s=10, label='test set')
plt.scatter(np.arange(split, num_data - 1), mlp_preds, color='violet', marker='x', alpha=0.4, s=20, label='MLP preds')
plt.scatter(np.arange(split, num_data - 1), gru_preds, color='green', marker='*', alpha=0.4, s=20, label='GRU preds')
plt.legend(loc='lower left')
plt.savefig('output_20_0.png')
plt.savefig('output_20_0.pdf')
plt.show()

附：以上文中的数据集及相关资源下载地址：
链接：https://pan.quark.cn/s/b485bdc0e8eb
提取码：NAn2