一、前言
这个LSTM系列是在学习时间序列预测过程中的一些学习笔记，包含理论分析和源码实现两部分。本质属于进阶内容，因此神经网络的基础内容不做过多讲解，想学习基础，可看之前的神经网络入门系列文章：
本系列重心放在解析LSTM算法逻辑、前向和反向传播数学原理、推导过程、以及LSTM模型的源码实现上。
本文详细讲解了LSTM源码的实现过程，以及数据在LSTM网络中流转的全过程，尽量做到每一行代码都讲解清楚，即是自己对知识做总结，也方便大家学习。本文是建立在前两篇文章的基础上，很多数学表达式不在重新推导，详细过程可查阅本系列第一篇和第二篇文章。
二、模型结构及训练数据说明
2.1.结构说明
本文构建的LSTM模型结构图如下所示：
从上图可以看出，模型核心部分包含三层LSTM结构，LSTM后接一个全连接层FNN，最后是一个softmax层，用于将输出结果映射成分布律，损失函数对应的是二元交叉熵函数。
上面的结构很明显可以看出是一个解决分类问题的模型，下一节我们来构建一个可用于分类的训练数据集。
2.2.训练数据构建
LSTM的强项是解决时间序列预测问题，但这里为了后续便于分析代码，我们构造一个相对简单的数字序列用于预测。
假设有两个小于等于50的随机数字，将这两个数字求和，如果两数之和大于60则输出1，如果两数之和小于60则输出0。示例如下：
34 23 => 0
45 34 => 1
34 33 => 1
10 13 => 0
11 24 => 0
44 46 => 1
此时我们就有了一个基本的数据集，根据设定好的参数，将生成的数据集划分为训练集和测试集，输入X.shape(32,1,2)，输出y.shape(32,),代码实现如下：
# 求和结果分类，x1+x2>60
def ClassifyData(self):
xArray = []
yArray = []
for _ in range(Params.TRAINING_EXAMPLES + Params.TESTING_EXAMPLES):
num1 = np.random.randint(0, 50)
num2 = np.random.randint(0, 50)
sum = num1 + num2
xArray.append([num1, num2])
if sum >= 60:
yArray.append(1)
else:
yArray.append(0)
# 监督学习数据 n*[X1, X2] -> n*[y] <=> X.shape(sample, 1 , 2) -> Y.shape(sample, 1, 1)
trainX = np.array(xArray[:Params.TRAINING_EXAMPLES]).reshape(Params.TRAINING_EXAMPLES, 1, 2)
trainY = np.array(yArray[:Params.TRAINING_EXAMPLES])
testX = np.array(xArray[Params.TRAINING_EXAMPLES:]).reshape(Params.TESTING_EXAMPLES, 1, 2)
testY = np.array(yArray[Params.TRAINING_EXAMPLES:])
return trainX, trainY,testX,testY三、网络结构参数设置
模型的参数设置如下，训练epoch为30，LSTM中隐藏节点数为30个，学习率为0.01，每个sample的batch大小为32，训练数据10000条，测试数据1000条，LSTM层数为3层。
EPOCH_NUM = 30 # EPOCH
MINI_BATCH_SIZE = 32 # batch_size
ITERATION = 1 # 每batch训练轮数
LEARNING_RATE = 0.01 # LSTM
VAL_FREQ = 5 # val per how many batches
# LOG_FREQ = 10 # log per how many batches
LOG_FREQ = 1 # log per how many batches
HIDDEN_SIZE = 30 # LSTM中隐藏节点的个数,每个时间节点上的隐藏节点的个数，是w的维度.
# RNN/LSTM/GRU每个层次的的时间节点个数，有输入数据的元素个数确定。
NUM_LAYERS = 2 # RNN/LSTM的层数。
# 设置缺省数值类型
DTYPE_DEFAULT = np.float32
INIT_W = 0.01 # 权重矩阵初始化参数
DROPOUT_R_RATE = 1 # dropout比率
TIMESTEPS = 1 # 循环神经网络的训练序列长度。
PRED_STEPS = TIMESTEPS # 预测序列长度
TRAINING_STEPS = 10000 # 训练轮数。
TRAINING_EXAMPLES = 10000 # 训练数据个数。
TESTING_EXAMPLES = 1000 # 测试数据个数。
SAMPLE_GAP = 0.01 # 采样间隔。
VALIDATION_CAPACITY = TESTING_EXAMPLES-TIMESTEPS # 验证集大小
TYPE_K = 2 # 分类类别
# 持久化开关
TRACE_FLAG = False
# loss曲线开关
SHOW_LOSS_CURVE = True
# Optimizer params
BETA1 = 0.9
BETA2 = 0.999
EPS = 1e-8
EPS2 = 1e-10
REG_PARA = 0.5 # 正则化乘数
LAMDA = 1e-4 # 正则化系数lamda
INIT_RNG=1e-4对应模型结构如下：
四、Lstm前向传播
4.1.初始化参数矩阵
随机初始化权重矩阵
和
，以及偏置项
。
：
- 对应的是“状态=>控制门”的权重矩阵，三层LSTM包含3个矩阵。
- 每层隐藏层为30个神经元，所以单个控制门对应的权重矩阵为(30,30)。
- 一个LSTM单元包含4个控制门，所以三层LSTM对应的权重分别是：(30,120)，(30,120)，(30,120)。
：
- 对应的是“输入=>控制门”的权重矩阵，三层LSTM包含3个矩阵
- 首层的权重矩阵，输入特征为2，隐藏层为30个神经元。所以单个控制门对应的权重矩阵为(2,30)，第二层和第三层LSTM对应的输入权重矩阵是(30,30)。
- 一个LSTM单元包含4个控制门，所以三层LSTM对应的权重分别是：(2,120)，(30,120)，(30,120)。
：
- 偏置项，同样包含4个控制门，b.shape(120,)
矩阵初始化代码实现如下，这里的初始化使用了矩阵奇异值分解的方式实现的。奇异值分解产生的矩阵有减少存储空间和方便计算的特点。
# 函数：np.linalg.svd(a,full_matrices=1,compute_uv=1)。
# 参数：
# a是一个形如(M,N)矩阵
# full_matrices的取值是为0或者1，默认值为1，这时u的大小为(M,M)，v的大小为(N,N) 。
# 否则u的大小为(M,K)，v的大小为(K,N) ，K=min(M,N)。
# compute_uv的取值是为0或者1，默认值为1，表示计算u,s,v。为0的时候只计算s。
# 返回值：
# 总共有三个返回值u,s,v
# u大小为(M,M)，s大小为(M,N)，v大小为(N,N)。
# -------
# A = u*s*v
# 其中s是对矩阵a的奇异值分解。s除了对角元素不为0，其他元素都为0，并且对角元素从大到小排列。
# s中有n个奇异值，一般排在后面的比较接近0，所以仅保留比较大的r个奇异值。
# 矩阵的奇异值分解可将一个大矩阵分解成三个小矩阵，减少了存储空间同时也便于计算
@staticmethod
def initOrthogonal(shape,initRng,dType):
reShape = (shape[0], np.prod(shape[1:]))
# 在区间范围内按reShape形状取样
x = np.random.uniform(-1 * initRng, initRng, reShape).astype(dType)
# x = np.random.normal(-1 * initRng, initRng, reShape).astype(dType)
# x = np.random.normal(0, 1, reShape).astype(dType)
# 矩阵奇异值分解
u,_,vt= np.linalg.svd(x,full_matrices =False)
w = u if u.shape == reShape else vt
w = w.reshape(shape)
return w4.2.LSTM的层间循环
LSTM的前向传播中，首层因为接收的参数和2层3层不同，因此对于首层我们总是要单独处理。
函数包含两个for循环，第一个for循环负责对三个LSTM层进行循环计算，第二个for循环负责对每一层LSTM中包含的
个时间步进行循环计算。 对于时间序列LSTM模型来说，在循环计算每个LSTM层之前都要对输入的状态参数
和
进行初始化操作，因此可知他的长序列预测能力集中在一个sample的
个时间步之间。而每个sample之间的状态是不能互相传递的。 在一个LSTM层内部计算过程中，共要计算T个时间步，当前时间步产生的
和
会传到下一个时间步中，作为下一个时间步的
和
输入。还有一点要注意的是，
个sample的数据会同时传入并行计算。比如
个sample包含有
个时间步，那么在计算第一个时间步
时，会将
个sample的第一个时间步都传入LSTM模型中进行计算。然后输出
个
和
传入下一个时间步进行计算。 本层LSTM计算完成后，最终会输出一个新状态
和
，
会在计算下个LSTM层时被舍弃掉，而Ht则会转换成下一层LSTM的新输入
，以此往复循环向后传递。在每一层LSTM计算过程中，都要将中间计算产生的参数记录到cache中，这些参数在反向传播时会用到。
#############################################################################
# 多层LSTM 多时间步 前向传播算法，此处需注意xh值在不同层对应的参数意义
# Input
# - x: 训练集输入数据 (N, T, D)
# - xh: 每层LSTM间传递的参数，首层xh为训练集输入x，后续xh为上层LSTM每个时间步的输出状态h(N, T, D)
# - h0: 首层LSTM传入的h状态 shape (N, H)
# - c0: 首层LSTM传入的c状态 shape(N, H)
# - h: 后续层LSTM传入的h状态 shape(N, T, H)
# - c: 后续层LSTM传入的c状态 shape(N, T, H)
# - Wx: x到f,i,g,o的权重矩阵 shape(D, 4H)，首层(2,120)，2,3层(30,120)
# - Wh: h到f,i,g,o的权重矩阵 shape(H, 4H)，3层都是(30,120)
# - b: 偏置项 shape(4H,)
# Returns a tuple of:
# - h: 每个时间步输出的状态 shape(N, T, H)
# - cache: 每个时间步产生的关于f,i,g,o门的中间参数
#############################################################################
def lstm_forward(self, x):
h, cache = None, None
# x.shape(32, 1, 2)
N, T, D = x.shape
# 根据权重矩阵中偏置项b的shape来获取Hidden层的节点数
H = int(self.lstmParams[0]['b'].shape[0] / 4) # 取整
# 首次计算时只存在输入x，不存在h，所以传入xh值为输入x,xh表示上个时间步的状态h
# 1. 首次输入xh=x.shape(32, 1, 2)
xh = x
'''
当前循环负责：3个LSTM层之间的参数传递，每层之间传参时，h,c,cache全都重新初始化；
其中xh表示每个LSTM层的输入，可以为x/h，h[]，c[]记录每个时间步生成的两个状态，cache存储f,i,g,o门产生的中间参数；
首层的xh为输入的训练数据x，后续层的xh为上一层LSTM在每个时间步所生成的状态h[,,]，二者shape可能会有差异
'''
for layer in range(self.layersNum):
# 每个LSTM层首次计算时要初始化当前h和c为0矩阵，类似reset_states作用，(N, T, H)=(32, 1, 30)
h = np.zeros((N, T, H))
c = np.zeros((N, T, H))
# h0,c0作为本层首个时间步的初始化参数，(N, H)=shape(32, 30)
h0 = np.zeros((N, H))
c0 = np.zeros((N, H))
cache = []
'''
当前循环负责：对N个sample中 每一个sample内部的 每个timesteps间的参数传递。每轮循环将一个batch中包含的N个sample同时传入，并行计算，最后也会同时输出N组预测结果；
这里的T对应的是每个sample中所包含的时间步timesteps个数，程序按照每个时间步来进行循环前向传播计算；
每次计算时将上一时间步输出的状态h[:,t-1,:], c[:,t-1,:]作为当前时间步的输入；当前时间步输出的状态存入h[:,t,:],c[:,t,:]中，作为后续的输入
'''
for t in range(T):
# (h0,c0).shape = (h[:,t-1,:],c[:,t-1,:]).shape = (32, 30)
# 每轮出参的h,cshape相同,h(32, 1, 30),c(32, 1, 30)，此例子时间步T为1，所以只进行一轮循环计算
h[:, t, :], c[:, t, :], tmp_cache = self.lstm_step_forward(xh[:, t, :],
h[:, t - 1, :] if t > 0 else h0,
c[:, t - 1, :] if t > 0 else c0,
self.lstmParams[layer]['Wx'], self.lstmParams[layer]['Wh'], self.lstmParams[layer]['b'])
cache.append(tmp_cache)
# 计算完当前LSTM层所有时间步，将每个时间步生成的h,c集合代入下一层的LSTM进行跨层运算
# 2.xh(32, 1, 30) xh为上个LSTM层每个时间步的状态集合h(32, 1, 30)
# 3.xh(32, 1, 30) xh为上个LSTM层每个时间步的状态集合h(32, 1, 30)
xh = h
##############################################################################
# END OF YOUR CODE #
##############################################################################
self.lstmParams[layer]['h'] = h
self.lstmParams[layer]['c'] = c
self.lstmParams[layer]['cache'] = cache
# 最终将最后一个LSTM层的每个时间步生成的h值返回，作为新输入传给FNN全连接层
return h4.3.LSTM的步间循环
LSTM内部控制门的功能函数实现如下，函数内部实现了4个控制门的计算逻辑，这个函数完成的是每一个时间步之间的计算逻辑。输入包含三个参数：
、
、
，其中
和
分别与权重矩阵
和
先进行仿射变换。
值得注意的是原先四个矩阵合并在一起，现在计算的时候需要拆分成4部分分别计算。单最红计算完后4个控制门的参数矩阵仍然是合并在一起的。
#############################################################################
# LSTM单元控制门内部算法实现
# Inputs:
# - x: 输入为训练数据x，或上一层LSTM生成的每个时间步的状态参数h集合, shape(N, D)
# - prev_h: 本层LSTM上个时间步生成的状态h, shape(N, H)
# - prev_c: 本层LSTM上个时间步生成的状态c, shape(N, H)
# - Wx: 输入xh-figo门权重矩阵, shape(D, 4H)
# - Wh: 隐藏状态h-figo门权重矩阵, shape(H, 4H)
# - b: 偏置项, shape(4H,)
# Returns a tuple of:
# - next_h: 当前时间步计算的h.shape(N, H)，传给下个时间步
# - next_c: 当前时间步计算的c.shape(N, H)，传给下个时间步
# - cache: 反向传播需要用到的f,i,g,o门参数，组成的数据集合
#############################################################################
def lstm_step_forward(self, x, prev_h, prev_c, Wx, Wh, b):
"""
LSTM计算技巧说明 ：
i_t = σ(W_xi*x_t + W_hi*h_(t-1) + b_i)
f_t = σ(W_xf*x_t + W_hf*h_(t-1) + b_f)
o_t = σ(W_xo*x_t + W_ho*h_(t-1) + b_o)
c^_t = tanh(W_xc*x_t + W_hc*h_(t-1) + b_c)
// g_t = tanh(W_ig*x_t + b_ig + W_hg*h_(t-1) + b_hg)
// c_t = f_t ⊙ c_(t-1) + i_t ⊙ g_t
#此处说明LSTM如何解决梯度消失原因，c_(t-1)表示过去信息，c^_t表示当前信息，此时c_t和c_(t-1)是线性关系而不再是乘积关系
c_t = f_t ⊙ c_(t-1) + i_t ⊙ c^_t
h_t = o_t ⊙ tanh(c_t)
通过前4个表达式可以看出，其实是x和h与f,i,g,o门对应的权重矩阵Wx和Wh进行了相同的矩阵运算，只是使用的权重矩阵不同，
所以我们可以构建一个4倍大小的W，将f,i,g,o门对应的4个W矩阵拼接起来，计算之后再将4个矩阵分别分离出来
这样可以减少计算量
"""
next_h, next_c, cache = None, None, None
# prev_h.shape(32, 30)
H = prev_h.shape[1]
# 合并之后的i,f,o,g在这里可以统一计算
# 1.matmul作矩阵乘法(32, 120)=(32, 2)⊙(2, 120) + (32, 30)⊙(30, 120) + (120,)
# 2.matmul作矩阵乘法(32, 120)=(32, 30)⊙(30, 120) + (32, 30)⊙(30, 120) + (120,)
# 3.matmul作矩阵乘法(32, 120)=(32, 30)⊙(30, 120) + (32, 30)⊙(30, 120) + (120,)
z = Tools.matmul(x, Wx) + Tools.matmul(prev_h, Wh) + b
# 之前将i,f,o,g四个矩阵合并了，这里将z(32,120)拆分4块进行计算，i,f,o,g的shape都是(32, 30)
# 计算方式见注释“计算技巧部分” of shape(N,H)
i = Tools.sigmoid(z[:, : H])
f = Tools.sigmoid(z[:, H : 2*H])
o = Tools.sigmoid(z[:,2*H : 3*H])
g = np.tanh( z[:,3*H : ])
# next_c(32, 30) = (32, 30)*(32, 30) + (32, 30)*(32, 30)
next_c = f * prev_c + i * g
# next_h(32, 30)
next_h = o * np.tanh(next_c)
# i,f,o,g门产生的中间参数
cache = (x, prev_h, prev_c, Wx, Wh, i, f, o, g, next_c)
##############################################################################
# END OF YOUR CODE #
##############################################################################
return next_h, next_c, cache五、FNN前向传播
LSTM层计算完成后，输出一个
矩阵传递到全连接层FNN中，全连接层只是先将张量降维:(32,1,30)=>(32,30)，然后又进行了一个简单的仿射变换，最终全连接层输出的矩阵是(32,2)
# 全连接层的前向传播,激活后再输出
def fp(self, input):
# 全连接层首先对输入进行拉伸变形处理，相当于Flatten()的功能，将(32,10,30)->(32,300)
self.shapeOfOriIn = input.shape
self.inputReshaped = input if self.needReshape is False else input.reshape(input.shape[0],-1)
# 先将输入矩阵与全连接层权重矩阵相乘，再进行激活函数运算
self.out = self.activator.activate(Tools.matmul(self.inputReshaped, self.w) + self.b)
return self.out六、损失函数用法及代码实现
上述内容可知，最终全连接层FNN会输出一个shape(32, 2)的预测结果矩阵，这个矩阵中对应了对数值预测的两种结果，即0或1。程序则要对这个预测结果进行评估检验，计算当前预测结果与真实值之间的误差大小，并将误差值反向传播给LSTM网络，使其修正参数，然后进行新一轮的学习。
1.在计算预测值和真实值之间的误差时，先用softmax函数对输出数值进行概率转换，然后获取最大概率所对应的数字就是预测值，最后来判断预测值是否正确。关于softmax为什么能概率转换可参考文章：
神经网络中的softmax层为何可以解决分类问题——softmax前世今生系列(3)导读：softmax的前世今生系列是作者在学习NLP神经网络时，以softmax层为何能对文本进行分类、预测等问题为入手点，顺藤摸瓜进行的一系列研究学习。其中包含：1.softmax函数的正推原理，softmax的代数和几何意义，softmax为什么能用作分类预测，softmax链式求导的过程。2.从数学的角度上研究了神经网络为什么能通过反向传播来训练网络的原理。3.结合信息熵理论... 2.有了预测值和真实值，可以使用二元交叉熵来计算预测结果和真实值之间的损失（误差）值，再将损失值均摊到32个预测数字上。就可以让我们直观的看到每次误差的变化。
3.二元交叉熵计算出的误差方便我们量化观察，而向模型中反向传播的误差，则是softmax输出的误差矩阵。对输出的误差矩阵(32, 2)中每个元素除以32，作为误差矩阵输出。后续就可以利用这个误差矩阵来实现反向传播的算法了。原理可以参考下面这篇文章：
BP神经网络中交叉熵作为损失函数的原理——softmax前世今生系列(4)导读：softmax的前世今生系列是作者在学习NLP神经网络时，以softmax层为何能对文本进行分类、预测等问题为入手点，顺藤摸瓜进行的一系列研究学习。其中包含：1.softmax函数的正推原理，softmax的代数和几何意义，softmax为什么能用作分类预测，softmax链式求导的过程。2.从数学的角度上研究了神经网络为什么能通过反向传播来训练网络的原理。3.结合信息熵理论...
"""
二元交叉熵损失函数
"""
class SoftmaxCrossEntropyLoss:
@staticmethod
def loss(y,y_, n):
y_argmax = np.argmax(y, axis=1)
softmax_y = Tools.softmax(y)
acc = np.mean(y_argmax == y_)
# loss
corect_logprobs = Tools.crossEntropy(softmax_y, y_)
data_loss = np.sum(corect_logprobs) / n
# delta
softmax_y[range(n), y_] -= 1
delta = softmax_y / n
return data_loss, delta, acc, y_argmaxsoftmax函数的实现代码
# 输出层结果转换为标准化概率分布，
# 入参为原始线性模型输出y ，N*K矩阵，
# 输出矩阵规格不变
@staticmethod
def softmax(y):
# 对每一行：所有元素减去该行的最大的元素,避免exp溢出,得到1*N矩阵,
max_y = np.max(y, axis=1)
# 极大值重构为N * 1 数组
max_y.shape = (-1, 1)
# 每列都减去该列最大值
y1 = y - max_y
# 计算exp
exp_y = np.exp(y1)
# 按行求和，得1*N 累加和数组
sigma_y = np.sum(exp_y, axis=1)
# 累加和reshape为N*1 数组
sigma_y.shape = (-1, 1)
# 计算softmax得到N*K矩阵
softmax_y = exp_y / sigma_y
return softmax_y七、总结
上述就是LSTM神经网络前向传播的核心代码实现部分。文章比较详细的介绍了输入数据在模型中流转的整个过程，在每一行代码的注释中都尽量写清了数据的变化形态。
本来准备一篇文章把前向传播和反向传播一起写，但全文超过3万字，考虑到内容太多会导致可读性变差，所以反向传播部分放在下篇文章中继续讲解。
参考文献：
GitHub - ljpzzz/machinelearning: My blogs and code for machine learning.
https://github.com/NLP-LOVE/ML-NLP
