高阶矩阵并行求逆的FPGA设计与实现

伍松，毛宇河

（1.广西科技大学 机械与交通学院，广西 柳州 545006；
2.广西汽车零部件与整车技术重点实验室（广西科技大学），广西 柳州 545006）

在现代数字信号处理中，许多先进的技术方法都涉及高阶矩阵运算，例如在声全息领域里的波叠加法近场声全息技术[1]中或在对碳纤维复合材料的转动轴进行模态分析时都会用到传递矩阵或传递矩阵的逆变换过程[2]，这些均与高阶矩阵运算息息相关。矩阵求逆过程亦是如此，其在数据分析、信号处理、系统理论、多元分析、现代控制等领域都有重要的使用价值，可见矩阵求逆的应用范围十分广泛。

但由于求逆过程自身的繁琐性和计算复杂性，因此大多数求逆运算都过于依赖通用型计算机中的软件开发平台（如MAT⁃LAB、LabVEIW等），并且在需要进行矩阵求逆等相关运算时，还需先掌握开发平台的配套语言和使用方式，增加了研究人员的设计难度。

另外矩阵求逆的传统方法大多都是通过最小均方算法或最小二乘法求解近似值，这样的求解过程往往精度低，收敛慢，难以达到设计要求[3]。同时在硬件实现方面，矩阵求逆也会受到资源及传输复杂程度的限制，相关的研究大都是针对特殊矩阵或者小规模矩阵进行实现的，例如，文献[4]是针对三角矩阵求逆进行硬件实现的过程，而针对高阶任意矩阵或大规模矩阵的相关文献却较少，因此对于众多计算密集型问题的处理中，研究高阶矩阵及大规模求逆的硬件实现过程十分具有挑战性[5]。

即针对高阶矩阵求逆问题，采用并行算法，基于现场可编程门阵列（Field Programmmable Gate Array，FPGA）的高速并行处理能力，采用流水线法以及LU分解求逆方法来完成基于FPGA[6]的8阶满秩方阵快速求逆处理器的设计与实现。这种快速实现高阶矩阵求逆过程，不仅解决了传统方法收敛性差、精度低的问题，而且还解决了矩阵求逆自身的计算复杂性以及高阶矩阵难以求逆的问题，同时也避免了求逆运算对于通用性计算机的过度依赖，能为很多使用到矩阵求逆过程的研究领域提供便利，提高运算效率，具有广泛的应用价值，对工程应用也有一定的帮助。

LU分解是矩阵求逆过程中使用较为广泛的一种分解求解办法。其基本原理是先将n阶方阵A分解成一下三角矩阵L同一上三角矩阵U的乘积形式后，利用三角矩阵易求逆的原则，分别求出下三角矩阵L和上三角矩阵U的逆形式，再通过矩阵相乘得到最终的逆矩阵A−1。即：

该方法可较大程度地减少求逆过程的运算量和复杂程度，同时能解决最小均方、最小二乘等方法求近似解精度低收敛性差的问题。

2.1 矩阵的LU分解

设n阶方阵的各阶顺序主子式均不为零，则A的LU分解则能仅包含一种结果。若以矩阵乘法原理为支撑，可推导得出LU分解的迭代结果算法[7]，即：

由上面的结果式可知，矩阵的LU分解是一种循环迭代过程，U矩阵和L矩阵中各结果元素的求得过程均是使用递推迭代，过程类似。但由式（4）、式（5）分析可知，L矩阵实际上要比U矩阵少一步，因此在此处的设计上要注意小心防止出错。

2.2 上三角及下三角矩阵求逆

假设下三角矩阵L的逆矩阵为R，上三角矩阵U的逆矩阵为S，虽然主要计算8阶矩阵的逆结果，但对于n阶上下三角逆矩阵而言，其结果形式可规律呈现。例如，n阶下三角逆矩阵R的一般形式为：

上三角逆矩阵S的结果式，其表达式同式（6）。

由式（6）结果式可以看出，计算上三角或下三角矩阵求逆过程基本可归纳为以下几步：

（1）首先由分解出的L矩阵和U矩阵首先可以求得逆矩阵R和逆矩阵S的对角元素；

（2）由对角元素，在j=1的情况下，以式（6）中递推式计算出R矩阵第1列或S矩阵第1行的结果元素；

（3）后再在j取2到n的顺序条件下，用通过（1）～（2）及此步骤计算得到的结果元素，按i=j−1至i=1的先后顺序倒序递推，便可得出R矩阵或S矩阵的其他元素，完成该过程。

由上面的计算过程得知，LU分解法的迭代过程多，且迭代次数会随着矩阵阶数的提升而成倍增加，如果想要利用该方法实现高阶矩阵的求逆过程，其迭代计算的复杂程度直接影响硬件设计过程的难易度。因此准确控制求逆中的迭代过程，降低其在设计上的复杂性，是设计的关键。此外，本设计在该过程中将采用并行计算设计，因此保证中间过程计算的准确性是设计的难点。

2.3 矩阵相乘

该过程是将上三角逆矩阵S与下三角逆矩阵R进行矩阵相乘，依靠公式A−1=U−1L−1=SR，来求得原始方阵A最终逆矩阵，进而完成求逆的计算过程。

那么由矩阵乘法原公式容易得到逆矩阵A的结果式为：

因此由式（7）得知，该过程所需的输入数据数量较多，因此控制好输入数的准确性和输入及储存的时效性将是本过程的关键。

2.4 并行运算

因本处理器中矩阵LU分解的输入输出过程及L矩阵与U矩阵的求逆过程是相对独立且互不干涉，若将每个过程都串行连接，在时钟控制下依次执行，虽能得到计算结果，却运算时间长，计算过程繁琐。

因此为了提高运算速度，考虑将独立互不干涉的过程并行执行，即对上述两个过程采用并行结构（Parallel Structure）[8]，形成并行形的数据输入过程及矩阵求逆过程，该过程将原串行过程改为两次并行过程，即可在原串行基础上提速四倍；
同时并行结构还具有并发性和同时性，是各领域尤其是大数据领域达到高性能及高实时性计算需求的理想选择。并行结构的示意图，如图1所示。

图1 LU分解输入输出并行结构示意图Fig.1 Parallel Structural Diagram of Input and Output Process of LU Decomposition

2.5 数的定标

我们知道，在Quartus II开发平台上，浮点数及小数是无法进行综合的，因此在使用vhdl语言描述硬件程序时，是不能在程序中使用和出现带小数点的数据的，但如将浮点型近似为整型后，其数据的精度将大大降低。因此为了能够在设计中使用浮点数和小数，采用定标法给予解决。

所谓定标法（Calibration Method），就是确定一个数小数点位置的方法[9]，而由定标数确定出来的数，称为定点数。若设定标数为ξ，定点数为Φ，原浮点数设为X，则在二进制中，浮点数与定点数的关系为：

由关系式看出，浮点数变为定点数后，其精度将随ξ的增大而增大，但随着ξ的增加，在位数一定的情况下，其原数据的涵盖范围也将减小，即数据的精度和范围是对立的。由式（8），也能推断出，该方法的适用范围并不局限于输入数为浮点数或小数，整数同样可以，若输入数为整数，则式（8）也可变为：

那么当整型数需要进行除法或开方运算想要保住小数位数不丢失时，可采用定标法，这便是思路，输入数据为整数，但因过程数据会产生小数，为保留住过程数据的精度不流失，保留小数位数，因此采用定标法。

2.6 任意高阶满秩方阵求逆法

虽然在设计上仅完成了8阶满秩矩阵的求逆过程，但最终目的是要依靠本设计去完成解决任意阶满秩矩阵的求逆工作。经研究，将本设计与高阶矩阵分解求逆方法结合起来，可以解决任意阶矩阵的求逆。下面将对此方法进行阐述。

首先需要简要介绍在解决高阶满秩矩阵求逆时所要用到的分解求逆方法。假设矩阵为N阶满秩方阵，且其主对角线上的K（K

式中：A、D—σ 阶和（N−σ）阶矩阵；
B、C—σ×(N−σ)阶和(N−σ)×σ阶矩阵，矩阵H−1=(D−CA−1B)−1—δ阶矩阵[10]，若令σ=8，则矩阵A的逆可由设计计算得到，因此只需求得H−1即可；
再将H继续按式(10）分解，又可得到下一个需求逆的8阶矩阵和一个（N−16）阶矩阵，8阶矩阵求解方式同上，同样又可按照式(10）进行分解；
以此类推，最后必然分解得到m个8阶矩阵和1个阶数小于8且大于0的矩阵R，然后再将矩阵R继续按照σ=4，σ=2及σ=1的次序继续进行分解求逆过程后便可得到最终结果。

其中，4阶矩阵、2阶矩阵及1阶矩阵的逆，可以简单通过修改8阶矩阵程序参数来获取，十分方便，因此综上可知对于N阶方阵的分解公式可为：

若令4α+2β+γ=t，则t可视为N除8后的余值，m视为商，便可快速得到分解次数，再按照上述的分解求逆方法便可得到任意阶满秩矩阵的求逆结果了。

本并行矩阵求逆的实现主要采用Altera公司提供的Quartus II 13.0软件开发平台来编写求逆IP核的设计文件，再通过设计好的IP核搭建出最终的处理器，该过程便是设计思路。

与理论算法对应，硬件实现电路包括4个主要模块：总控制模块（Total FSM，TFSM）、LU 分解分模块（LU Decomposition，LU Decp）、上下三角矩阵求逆分模块（LU Inversion，LU Inv）和矩阵乘法分模块（Matrix Multiplication，Matrix Mult）。其中，总控制模块包含总体状态控制器（Total State Controller，TSC）和总体地址控制及状态转换计数器（Total Address&State Counter，TASC）；
LU分解分模块包含LU分解器（LU Resolver，LU RS）及其RAM存储器，其中RAM存储器包含L矩阵存储器L_RAM和U矩阵存储器U_RAM；
上下三角矩阵求逆分模块包含其自身求逆控制器（In⁃version FSM，IFSM）、求逆模块（Inversion Model，Inv Model）及其RAM存储器，其中求逆模块和存储器模块分别包含L矩阵求逆部分和U 矩阵求逆部分；
矩阵乘法分模块包含（8×8）矩阵乘法器（Matrix Multiplication，Mx Mult）及其RAM存储器；
硬件系统的总体结构，如图2所示。

图2 矩阵求逆系统总体结构示意图Fig.2 Overall Structure Diagram of Matrix Inversion System

其中，三个分模块中的主要模块：LU分解器、矩阵求逆模块及矩阵乘法器，在设计上对其程序都进行了参数化设计，因此容易根据实际情况，通过简单程序更改，实现不同阶数的矩阵求逆过程，提升了程序自身的可重构能力和实用性。

另外，对于输入矩阵A 进入本系统的输入办法，采用的是TXT文本输入方式，即将写着输入数据的TXT文本文件由程序进行识别，通过时钟信号的具体控制，使文件中的数据准确稳定的传输到程序当中；
且由于其独特的输入方式，输入数据可根据实际情况随时进行更改，也可将其他软件得出来需要求逆的数据文件，通过格式转换，转换为TXT文件格式，之后便可通过本程序快速得出求逆结果。因此本设计的实用性、稳定性和便利性都是十分突出的。根据由IP核的总体结构示意图以及要实现的功能，此设计的总体工作流程，如图3所示。

（1）首先，在复位信号的控制下，总系统归于初始状态，且TSC处于S0的稳定复位状态。

（2）在时钟信号的总体控制下，TASC中计数信号加一，并发送转换信号q，TSC接收信号后使总体状态转换为S1，待状态平稳后，控制L_RAM从TXT文件中稳定接收输入矩阵A的各元素数据并加以存储，为接下来的计算做准备。

（3）在TASC计数信号的控制下，TSC继续按S2至S4的顺序发生循环状态转换，LU分解过程便依序稳定进行，其中LU分解的输入输出过程采用并行结构，以增加数据传输的效率性；
且由于需要计算的是8阶矩阵的逆，因此该过程需循环8次，当然也可根据实际情况，通过对本程序进行简要修改来更改矩阵阶数，以实现不同阶数矩阵的求逆，从而增加本设计的可重构性能。

此步骤便是LU分解分模块的工作过程，同时当LU分解结束时，TSC的状态信号转换成为S4；

（4）之后，系统进入上下三角求逆过程，由TSC发出的使能信号，控制启动LU 求逆分模块的内部控制器（IFSM），使其以内部状态转换稳定控制内部RAM接收并存储LU分解分模块求得的各数据结果；
之后，在IFSM状态循环转换的控制下，内部的求逆运算过程循环执行，且其循环次数同步骤（3），同样与矩阵阶数一致，且求逆过程同样采用并行结构，使L矩阵求逆和U矩阵求逆双过程同时进行，提高了计算过程的效率性。当执行循环过程后，求逆控制器将给总控制器发送信号，一方面使总控制器继续执行下面步骤，一方面使矩阵乘法分模块再次复位清空数据，使之后矩阵乘法过程能稳定进行；

（5）当LU求逆分模块工作完成后，TSC的内部状态继续由S4转换为S5，待状态稳定后，矩阵乘法分模块在其控制下开始接收LU求逆的结果数据；
之后，状态变量再次自S6至S12循环转换，控制矩阵乘法分模稳定进行计算过程，由式（7），因此乘法器计算采用边相乘累加边储存的方法，按行方向依次依序得到结果数据，并存进结果矩阵，循环往复；
待计算全部完成后，通过输出端口输出最终结果，进而完成整个求逆过程。其流程步骤，如图3所示。

图3 矩阵求逆系统总体工作流程图Fig.3 Overall Work Flow Chart of Matrix Inversion System

其中，由于LU分解过程及LU矩阵求逆过程都包含了除法过程，因此为保证数据的精度不流失，采用了定标法来进行计算，过程数据利用定标数采用先放大后缩小的方式来尽量减少中间过程因小数而产生的误差损失，以提升结果的有效性。且为了保证快速性，在LU分解过程和LU求逆过程都采用了并行结构，最大程度的保证求逆过程的时效性和本设计的实用性。

设计的求逆处理器，是由8×8并行满秩方阵求逆的IP核搭建而成的。其IP核的RTL级描述原理图，可由四个分模块共同构成，其RTL级描述简图，如图4所示。

图4 矩阵求逆系统RTL级描述简图Fig.4 RTL Level Description Diagram of Matrix Inversion System

①为总体控制单元，是本系统的主控制单元，主要通过状态转换，为其他单元提供使能信号和地址信号。

①的上方TSC模块表示总状态转换控制单元，用于执行状态转换，提供使能信号，下方TASC模块表示总地址控制及状态转换计数单元，用于通过计数信号控制状态转换控制单元执行状态转换，并为其他模块提供地址信号；

②为LU分解总单元，用于执行LU分解求逆中的LU分解过程工作，其共分为三个模块，L_ram和U_ram两模块分别用于接收分解得到的L矩阵与U矩阵的结果，其中模块L_ram同时还包含从文件接收输入矩阵A及输出输入矩阵A两项工作，模块LU_Decp为LU分解模块，是LU分解单元的主要功能模块，同时对其采用了并行结构与定标法，以提升运算速度及运算精确度；

③共同形成LU求逆总单元，其中ISC和IASC为该单元中的内部控制部分，其原理与①相同，③中其他模块为该单元中的主要执行部分，其中Inv_ModelU与ram_inv_u两模块和Inv_ModelL与ram_inv_l是并行执行关系，他们分别作为U 矩阵和L 矩阵的求逆执行部分，以L 矩阵求逆部分为例，ram_inv_l 为RAM 存储模块，用于对运算数据和结果数据进行存储、输入和输出，Inv_ModelL 为L 求逆部分的主要功能模块—求逆模块，用于完成其求逆工作；

④中三模块共同形成为乘法器总单元，其中mult_ram 为RAM存储模块，r_sct为行选择模块，matrix8×8为（8×8）矩阵乘法模块，其中，RAM存储器的功能除了包含数据的传输存储，还增添了为乘法模块提供地址的第二用途，行选择模块是为了保证矩阵乘法过程中在RAM循环提供右乘矩阵列方向元素时，其提供的左乘矩阵的行方向元素能准确按时供应，矩阵乘法模块为该总单元的主要功能模块。

为了验证该8×8满秩矩阵求逆处理器的功能性要求，将利用Modelsim仿真平台对设计出来的IP核进行仿真测试。主要测试过程、测试结果及对结果的分析结论简述如下：

设计选用Stratix高端系列芯片为目标开发板，以方便在设计上降低对编译速度、I/O口数、记忆位数（Memory Bits）及LEs逻辑元使用数等各参数的限制，增加在功能实现方面的有效性。芯片选用Altera Stratix III 系列的EP3SL340F1417C2 为目标芯片，通过在Quartus II 13.0开发平台上编译、例化、综合及在Modelsim上以10MHz为时钟频率进行时序仿真分析后，得到的最终结果，如图5所示。

图5 矩阵求逆系统仿真计算结果Fig.5 Simulation Results of Matrix Inversion System

图中的结果的显示方式为，右斜列方向8个相邻数据为结果矩阵的一行，例如图中右斜第一列数：−131731，138855，−30116，88090，16053，−45763，−3422，−74240即为实际逆矩阵的第一行结果，其他行结果类似可得。图中的所示结果是将原结果扩大了5122倍且截断了小数位数后的显示结果，即相当于结果作为了以18为定标的定点数。这样显示的目的是为了解决在硬件系统无法处理浮点数的情况下，浮点数能以更高的精度在硬件中显示和运算的问题。其中，输入矩阵：

若将输入矩阵扩大5122后放入MATLAB中，以LU分解的求逆方法求解结果，得到的结果，如图6所示。

图6 MATLAB中LU分解求逆计算结果Fig.6 Calculation Results of Inverse Method of LU Decomposition in MATLAB

若将本设计的结果矩阵与MATLAB求解矩阵的各元素按行方向顺序展开成行向量后，在MATLAB中进行误差比较，可得到结果对比图和误差百分比图，如图7所示。计算出的百分误差的平均值为−0.43%。可见用本文设计求得的结果值与MATLAB按此方法求得结果值的误差十分小。若用MATLAB的INV函数去获得结果，得到的结果值，如图8所示。

图7 MATLAB中LU分解求逆结果与设计求逆结果误差对比Fig.7 The Error Comparison Between the Result of LU Decomposition in MATLAB and that of the Design

图8 MATLAB中INV函数求逆计算结果Fig.8 Inverse Calculation Results of INV Function in MATLAB

将此结果按行方向顺序展开成行向量后，与本设计结果进行误差比较，同样可得到结果对比图和误差百分比图，如图9所示。

图9 MATLAB中INV函数求逆结果与设计求逆结果误差对比Fig.9 The Error Comparison Between the Result of INV Function in MATLAB and that of the Design

计算出的百分误差的平均值为5.94%。可见用此设计方法求得的结果值与实际结果值的误差较小。此外，使用本设计求得结果耗费的时间与MATLAB 按两种方法求解花费的时间对比，如表1所示。

表1 本设计仿真计算时间与MATLAB求解时间比较Tab.1 The Comparison Between the Simulation Calcula⁃tion Time of the Design and the Time of MATLAB Solution

由上表可知，本设计的计算效率已远超于MATLAB。通过图5的数据结果可知，本设计实现了处理器能够计算高阶矩阵的逆，证明了其具有高速计算能力；
同时通过图5、图6、图8中结果数据的对比分析证明了使用定标法计算出来的数据在扩大了512倍后与MATLAB以两种方式计算得出的数据十分接近；
且由图7、图9中各结果元素的结果对比图和误差百分比图，可以证明该设计的计算结果误差较小精度较高；
且再通过对表1中仿真计算时间与MATLAB 用两种方法计算所花费的时间进行对比，可以得知，使用FPGA及并行结构设计出的矩阵求逆，其速度已远超MATLAB求解的速率，从0.598s缩短到0.645μs，达到了高速求解的目的。通过以上分析可知，本设计能够达到了预期的高速性和高效性。该设计虽然满足求解的高精度要求，但从图7与图9也看出，所求得数据还存在偏差，图9中有极个别数据的偏差较大。存在偏差的原因在于：

在LU分解和LU矩阵求逆过程中，虽然使用了定标法尽量保证了小数位数的保留，但仍有个别较低位的小数部分丢失，遇到较小浮点数与较大整型数的乘积运算时，其结果误差就偏大。因此在后续的设计中，会对其做进一步优化，减小因定标法截断误差产生的计算误差，使本设计更加具有实用价值。

采用LU分解、并行结构、定标等设计方法实现了一个高精度高效率的基于FPGA 支持的8 阶满秩矩阵求逆处理器的研究设计。该处理器利用了FPGA的并行能力，使矩阵求逆过程得以快速实现；
同时该设计还具有计算稳定性高，设计结构清晰，可重构性能强等特点。该设计不仅充分利用了FPGA芯片的高速并行性能，解决了高阶矩阵难以求逆的问题，同时不再依赖繁琐的通用计算机运算，对实际工程具有一定的意义和价值。

猜你喜欢定标高阶乘法算乘法小学生学习指导(低年级)(2022年10期)2022-11-05我国为世界大豆精准选种“定标”今日农业(2021年19期)2022-01-12我们一起来学习“乘法的初步认识”数学小灵通(1-2年级)(2021年10期)2021-11-05有限图上高阶Yamabe型方程的非平凡解数学物理学报(2021年1期)2021-03-29高阶各向异性Cahn-Hilliard-Navier-Stokes系统的弱解数学物理学报(2020年6期)2021-01-14基于恒星的电离层成像仪在轨几何定标空间科学学报(2020年1期)2021-01-14滚动轴承寿命高阶计算与应用哈尔滨轴承(2020年1期)2020-11-03《整式的乘法与因式分解》巩固练习语数外学习·初中版(2020年11期)2020-09-10把加法变成乘法小学生学习指导(低年级)(2019年9期)2019-09-25一类完整Coriolis力作用下的高阶非线性Schrödinger方程的推导数学物理学报(2018年5期)2018-11-16