基于FPGA的LSTM加速器设计（MNIST数据集为例）

摘要

本文以MNIST手写数字识别任务为例，使用FPGA搭建了一个LSTM网络加速器，并选取MNIST数据集中的10张图片，通过vivado软件进行仿真验证。实验结果表明，本文设计的基于FPGA的LSTM网络加速器可以完成图片分类任务，其准确率为80%(20张图片，4张分类错误)。本文主要分为四部分，第一章为LSTM硬件加速器的原理介绍，第二章为软件部分的程序设计思路，第三章为FPGA硬件部分的设计思路。本文所设计的LSTM硬件加速器的完整的工程文件已上传，并在文末对工程文件进行了简单的介绍。

一、基于FPGA的LSTM加速器设计原理

基于FPGA的LSTM加速器设计流程如下所示：

• 1. 在Pytorch框架下搭建LSTM网络模型，并用GPU进行模型训练；
• 2. 对训练好的LSTM网络模型进行剪枝压缩，并通过重训练恢复精度；
• 3. 使用分段线性函数替代非线性激活函数，并通过重训练恢复精度；
• 4. 对模型的权值进行量化，并导入到FPGA的ROM资源中；
• 5. 在FPGA上搭建LSTM硬件加速器；
• 6. 对搭建的硬件加速器进行验证评估；

本章主要对各流程中的原理部分进行介绍，本章尽可能的解释了FPGA实现LSTM硬件加速器所需的原理知识，详细的介绍可以参考文末给出的参考资料。

1. 长短期神经网络（Long Short Term Memory，LSTM)原理

传统循环神经网络（Recuurent Nerual Network，RNN）通过将上一时间步的输出$h_t$作为输入的一部分，和当前输入$x_t$一起作为输入信息输入到网络中，从而能够捕获序列信号的特性。然而，传统RNN存在梯度消失和梯度下降的问题，
而LSTM通过引入记忆细胞机制缓解了传统循环神经网络（Recuurent Nerual Network，RNN）梯度消失和梯度爆炸的问题。
LSTM的表达式如下所示：

其中$i_t$为输入门，取值范围为（0~1），表示更新记忆细胞的程度。$f_t$为遗忘门，取值范围为（0~1），表示上一时间步的记忆细胞$c_{t-1}$的剔除程度。${\tilde{c}}_t$表示待更新入记忆细胞的信息。$o_t$为输出门，决定记忆细胞与输出信息的关系。LSTM通过门控单元，能及时的剔除记忆细胞中的无用信息，并及时准确地更新信息，从而能缓解传统RNN的梯度消失和梯度爆炸的问题，但也因此引入了大量的参数，导致其难以直接运行在存储、计算资源受限的平台，例如FPGA。

2. 对LSTM网络模型中权值矩阵的剪枝（Top-k剪枝）

由于神经网络具有很强的鲁棒性，即使被大幅度的压缩，也能保证其准确率。$ELST{M}^{[1]}$中提出了一种top-$k$剪枝方案。该剪枝方案将权值矩阵的每个相邻的c个权值分为一组，每组只保留前k个绝对值最大的非零权值，其余权值均设为0。

c=8,k=2的top-$k$剪枝示意图如上所示。存储top-$k$剪枝方案压缩后的权值矩阵时，只需要3bits（$lo{g}_{2}c=3$）即可表示该非零权值的位置信息。

3. 对LSTM网络模型中权值的量化

除了通过剪枝压缩模型外，还可以通过量化来压缩模型。由于GPU训练的模型为32位浮点数，而FPGA处理数据大多是定点数，因此需要对模型进行量化。常用的量化方式用两种，一种是线性量化，一种是对数量化。

3.1 对数量化原理

对数量化的表达式如下所示：

其中m,f分别代表量化后的整数位位数和小数位位数。这种量化方式量化后的数为2的幂次方，如0011就表示$2^{-3}$，而与2的幂次方的乘法运算，可以用移位运算替代，如$a\ast 2^{-b}=a>>>b$。

如图所示，蓝色线表示函数 $y=log{Q}_{m,-f}(x)$，橘色线条表示函数$y=x$。由图可知，对数量化对数值较大的数，量化产生的误差较大。因此，只被用于量化LSTM网络的权值矩阵参数。

3.2 线性量化原理

线性量化的表达式如下所示：

其中m,f分别代表量化后的整数位位数和小数位位数。

如图所示，蓝色线表示函数 $y=Q_{m,f}(x)$，橘色线条表示函数$y=x$。线性量化的误差较稳定，因此输入、输出、中间运算结果均采用线性量化。

4. 非线性激活函数近似法

LSTM中的$\sigma$和$\tanh$函数，均为带有$e^x$的指数运算函数，FPGA难以直接实现这种复杂函数。常用的解决方案是用分段线性函数替代这两个非线性激活函数。如下所示：

为了便于FPGA的实现，本文采用Ptanh函数替代$tanh$，用$Hsigm$替代$\sigma$。

二、Pytorch框架下的LSTM实现手写数字图片分类任务（MNIST数据集）

本章主要对软件部分的设计进行介绍，并对关键代码进行解释，完整代码参考工程文件。

1. 分类器模型的搭建

本文搭建的LSTM网络模型包含一个输入维度为28，隐藏层维度28的单层单向的LSTM层，一个输入维度28，输出维度10的全连接层（Fully Connect，FC），以及一个4bit输出的分类器（输出最大值的位置信息，0~9）。

class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.rnn = LSTM(28, 28, num_layers=1,
                         bias=True, return_sequences=True, grad_clip=None)
        self.fc = nn.Linear(28, 10)
    def forward(self, x):
        zeros = Variable(torch.zeros(x.size()[0], 28, device=torch.device("cuda")))
        initial_states = [(zeros, zeros)]
        x = x.squeeze(1)
        hidden, _ = self.rnn(x, initial_states) #LSTM层
        x = self.fc(hidden[:, -1, :]) #全连接层        
        return x
net = Net()
outputs = net(images)
 _, predicted = torch.max(outputs.data, 1) #分类器

由于官方未提供（也可能是我没找到）LSTM对数量化的API，因此本文选用的LSTM为自定义的程序文件（详见工程文件）。

2. 对分类器的剪枝压缩（Top-k剪枝）

训练好的网络模型难以直接运行在资源受限的平台，因此需要剪枝压缩。top-$k$剪枝压缩代码如下所式：

def topk(para, k):
    c = torch.zeros(para.size()[0], para.size()[1],dtype = torch.int) #初始化一个和权值矩阵相同大小的掩膜矩阵
    l = int(para.size()[1]/7) #将每行的每7个权值分为一组，l为分组的数量
    parameter = torch.abs(para)  #将权值矩阵取绝对值
    _, b = torch.topk(parameter[:,:7], k, 1, largest = True) #b为0~6之间的k个数，表示该组最大的前k个权值的位置
    for i in range(1,l):
        _, b1 = torch.topk(parameter[:,i*7:(i+1)*7], k, 1, largest = True) #遍历每一组最大的前k个值的位置
        b1 = b1 + i * 7  #得到每一行中保留的权值位置信息的绝对位值
        b = torch.cat((b,b1),dim=1) #将每一段拼接起来

    for j in range(c.size()[0]):
        c[j, b[j, :]] = 1 #将c中，b中位置信息的对应的位置，置1（保留），其他部分为0
    return c

由于信息维度为28，因此本文选取每7个权值为一组，每组保留k个权值。本文选取的k为4，通过函数topk生成掩膜矩阵，再根据pytorch提供的API进行剪枝训练。

c1 = topk(net.rnn.cell0.weight_ix, 4) #根据W^{i}生成掩膜矩阵

class FooBarPruningMethod1(prune.BasePruningMethod): #调用剪枝的API，更改掩膜矩阵
    """Prune every other entry in a tensor
    """
    PRUNING_TYPE = 'unstructured'

    def compute_mask(self, t, default_mask):
        mask = c1
        return mask
        
FooBarPruningMethod1.apply(net.rnn.cell0, 'weight_ix') #对权值矩阵W^{i}进行剪枝

权值矩阵$W^i$的剪枝程序设计如上所示，以此类推，对LSTM其他权值矩阵，及FC层的权值矩阵进行剪枝。

3. 非线性激活函数近似法的程序设计

Ptanh和Hsigm函数如下所示，本文将剪枝后的网络中的非线性激活函数用这两个分段线性函数进行替代，并通过重训练恢复模型精度。

def Ptanh(inn):
       we1 = inn < -2.5
       we2 = (inn >= -2.5) & (inn < -0.5)
       we3 = (inn >= -0.5) & (inn < 0.5)
       we4 = (inn >= 0.5) & (inn < 2.5)
       we5 = inn >= 2.5
       out1 = -1
       out2 = 0.25 * inn - 0.375
       out3 = inn
       out4 = 0.25 * inn + 0.375
       out5 = 1
       out = out1 * we1 + out2 * we2 + out3 * we3 + out4 * we4 + out5 * we5
       return out

def Hsigm(inn)
      out = torch.clip(0.25 * inn + 0.5, 0, 1, out=None)
      return out

4. 对数量化的程序设计

def log2_Q(a):
   a = a.to('cpu')
   b = a.detach().numpy()
   e = np.sign(b)
   b = np.clip(np.round(np.log2(np.fabs(b))+0.4),-7,0) #得到最接近原始a的2的幂次方，不改变a的其他属性，因此只使用data属性
   b = np.power(2,b) * e
   a.data = torch.from_numpy(b).data
   return a

对数量化的函数如上所示，权值被量化为4bits数，详见[1]，然而量化后的权值矩阵不能再训练，因此本文先量化的候选记忆细胞的权值矩阵（$W^c$、$U^c$），之后对其他剩余的权值进行重新训练；再量化输出门和遗忘门的权值矩阵（$W^o$、$U^o$、$W^f$、$U^f$），并对其他剩余的权值进行训练，最后再量化输入门的权值矩阵（$W^i$、$U^i$）。

5. 线性量化的程序设计

def Q(a):             #输入训练参数，输出量化后的训练参数
    a = a.to('cpu')
    b = a.detach().numpy() #由于a是训练参数，requires_grad为True，因此不能直接用numpy函数操作，需转换
    b = np.clip(b,-0.875,0.875) #0.875是1 - (1/2)^3
    b = np.round(b * 8 + 0.5) / 8
    a.data = torch.from_numpy(b).data     #得到最接近原始a的定点数
    return a

权值被量化为4bits数，详见[1]，本文对FC层的权值矩阵进行线性量化，将其量化为4bits数。

6. 压缩后模型的数据导入FPGA的程序介绍

由于FPGA中存储和运算为二进制补码形式，因此需要将权值转换为补码形式，转换程序如下所示：

def p_d2b(n, m, f): #将一个10进制正数转换为一个2进制数，保留m位整数，f位小数，首位符号位
    b = []
    x = 2
    n = n * np.power(2, f)
    n = int(n)
    while True:
        s = n // x
        y = n % x
        b = b + [y]
        if s == 0:
            break
        n = s
    b.reverse()
    if(len(b) > (m+f)):
        for i in range(m+f):
            b[i] = 1
            b = b[:m+f]
    elif(len(b) < (m+f)):
        for i in range(m+f-len(b)):
            b.insert(0,0)
    b.insert(0,0)
    a = [str(i) for i in b ]
    return a

def n_d2b(n, m, f): #求一个10进制负数转换为一个2进制补码形式，保留m位整数，f位小数，首位符号位
    n = -1 * n
    b = p_d2b(n, m, f)
    b[0] = '1'
    flag = 1
    for i in range(len(b)-1,0,-1):
        if b[i]== '1' and flag == 1:
            b[i] = '1'
            flag = 0
        elif b[i] == '0' and flag == 1:
            b[i] = '0'
            flag = 1
        elif b[i] == '0':
            b[i] = '1'
        else:
            b[i] = '0'
    a = [str(i) for i in b ]
    return a

def d2b(n, m, f): #求一个数n的补码，保留m位整数，n位小数，首位符号位
    if n < 0:
        c = n_d2b(n, m, f)
    else:
        c = p_d2b(n, m, f)
    return c

对数量化后的数的存储与线性量化不同，如0.125（$2^{-3}$），应该存储为0011，表示右移3位即可完成乘法操作。对数量化后的数据的二进制生成函数如下所示：

def logd2b(n, f):
    if n > 0:
        n = np.log2(n)
        n = np.floor(-1 * n + 0.5)
        if n >= np.power(2,f):
            n = np.power(2,f)
        a = p_d2b(n,f,0)
    else:
        n = -1 * n
        n = np.log2(n)
        n = np.floor(-1 * n + 0.5)
        if n >= np.power(2, f):
            n = np.power(2, f)
        a = p_d2b(n, f, 0)
        a[0] = '1'
    return a

通过以上方式生成FPGA的ROM可加载的coe文件，将GPU训练出的网络模型参数导入到FPGA中。

def output_file_log(weight,name):
    name = 'coe/'+name
    f1 = open(str(name)+"_data.coe","a")
    f2 = open(str(name)+"_index.coe","a")
    data =';\nmemory_initialization_radix = 2;\nmemory_initialization_vector='
    f1.writelines(data)
    f2.writelines(data)
    para = weight.numpy()
    for i in range(para.shape[0]):
        f1.writelines('\n')
        f2.writelines('\n')
        for j in range(para.shape[1]):
            if para[i,j]!=0:
                data = logd2b(para[i,j], 3)
                index = p_d2b(j%7,3,0)[1:] #topk剪枝后的非零权值需要3bit表示其在分组中的位置信息
                f1.writelines(data)
                f2.writelines(index)
    f1.writelines(';')
    f1.close()
    f2.writelines(';')
    f2.close()

def output_file_Q(weight,name):
    name = 'coe/' + name
    f1 = open(name + "_data.coe", "a")
    f2 = open(name + "_index.coe", "a")
    data = ';\nmemory_initialization_radix = 2;\nmemory_initialization_vector='
    f1.writelines(data)
    f2.writelines(data)
    para = weight.numpy()
    for i in range(para.shape[0]):
        f1.writelines('\n')
        f2.writelines('\n')
        for j in range(para.shape[1]):
            if para[i, j] != 0:
                data = d2b(para[i, j], 2,13)
                index = p_d2b(j % 7, 3, 0)[1:]#topk剪枝后的非零权值需要3bit表示其在分组中的位置信息
                f1.writelines(data)
                f2.writelines(index)
    f1.writelines(';')
    f1.close()
    f2.writelines(';')
    f2.close()

三、Xilinx FPGA上的LSTM硬件加速器设计

LSTM整体框架如图所示：
其中LSTM层的结构如下所示：

输入数据经过S-P-X转换为多维向量，和非零权值的位置信息一起送入KMUX单元进行筛选，筛选后的非零权值和输入信息在MVMs模块中完成矩阵乘加运算，并在EWU模块中完成激活，点乘等运算，计算出的记忆细胞的值存储在FIFO-C中，输出信息$h_t$则存储在S-P-H中，以作为下一个时间步的输入。全连接层和LSTM的架构大体相似。
输入‘4’的图片，仿真结果如图所示，结果正确。

硬件部分没有太多要说的，需要注意的有如下几点：

• 1. 位宽问题：如果整数位位宽设置太少，会出现溢出，小数位位宽较少，则会截断一些信息，从而导致误差，因此中间结果的位宽，需要根据实际情况进行调整；
• 2. 补码的移位运算：本文将LSTM的权值矩阵量化为2的幂次方，从而可以用移位运算来替代乘法运算，可有一个问题，比如负数补码的移位，直接使用有符号的移位运算，出来是错的。比如一个1符号位，3整数位的数1111（-1），右移4位后是1111（-1，带符号的右移，补充符号位），而这并不是我们想要的结果，我们想要的是-1右移4位是-0.0625，截断后位0，而不是依然-1，而原码的移位则运算正确。因此，本文LSTM的乘法运算（用移位运算替代），先将输入转换为原码形式（C2t模块），进行移位操作后，再转换为补码形式（t2c模块）。
• 3. 其余的比较复杂的就是时序了，需要慢慢的捋，完整工程文件已上传，没有积分的可以评论区留邮箱发。

四、工程文件介绍

完整工程文件：基于FPGA的LSTM加速器设计（MNIST数据集为例）

• python (软件代码)
• • 1-MNIST-LSTM.py (初始网络模型)
• • 2-MNIST-LSTM.py (top-k剪枝)
• • 3-MNIST-LSTM-topk-linear.py （分段线性函数替代非线性函数）
• • 4.5.6.7.8.9.10.11的py文件为量化程序
• • 12-output-weight.py （权值导出函数）
• • 13-out-img.py（图片数据导出程序）
• FPGA （硬件代码）
• • 其中MNIST为顶层模块，Test为仿真验证程序

参考资料

1. ELSTM（topk剪枝）