基于FPGA 的复数矩阵求逆设计

周金强，凤继锋

（四川大学电子信息学院，成都610065）

0 引言

矩阵运算[1]在现代科学工程领域的基本问题，如何更快，效率更高地进行矩阵运算是科学计算研究的重点。在信号处理和图像识别领域中，矩阵乘法、求逆是最基本的数学工具，大量算法都涉及到这些基本运算，因此如何更快、更节省运算资源来进行矩阵的求逆运算是现代工程计算领域一直在研究的重点。但是，当今信号处理算法大部分是在PC 端的软件上进行实施的，在一些讲究实时性的场景，如无人机电子侦查，电磁频谱感知[2]等领域，这些场景对实时性要求高，并且受到环境影响，无法时刻装配大型计算机来进行信号算法运行，如调制宽带转换器（Modulated Wideband Converter，MWC）[3]，其混频、滤波、采样都在硬件上执行，但恢复算法则是在MATLAB 端进行，无法实时恢复出输入信号的频谱，这对该系统的实际应用形成了限制，使其无法满足特定场景下实时恢复出信号的频谱图，无法满足实时性要求，而FPGA 平台则可以满足这些特点，其功耗，运行速度都有优势，并且能够适应更为复杂的环境因素，通过将其使用的OMP 算法进行FPGA 平台移植可使该系统满足实时性场景，而OMP算法的硬件化的关键问题在于最小二乘部分，该部分在数学上表示为矩阵求逆计算，若用传统定义法或者SVD 分解法之间进行求伪逆需要消耗大量资源，使得硬件运算速度下降，因此如何选择合适的分解求逆法，并在此基础上设计高速，逻辑资源使用合理的求逆模块是本文的重点。

1 矩阵分解求逆数学表示

1.1 改进型Cholesky分解

OMP 算法中所遇到的最小二乘问题矩阵为非方阵，可以通过Moore-Penrose 求解方法将伪逆问题转化为求解共轭对称方阵的逆，而对于该类型的矩阵，改进型Cholesky 方法[4]有着较好的分解效果，相比较传统Cholesky 分解法，改进型Cholesky 分解规避了开方计算，而开方计算在硬件执行上需要消耗大量资源，复杂运算只涉及到倒数计算，并且由于矩阵对角线都为正实数，因此所消耗的除法计算与矩阵维度一致，相比较其他分解法如QR 分解[5]，SVD 分解，复杂计算较少，因此本文选择改进型Cholesky 分解法，下面介绍其数学原理：

假设可逆方阵为A，其元素以aij来表示，将矩阵进行分解为A=LDLH的形式：

其中L 为对角线元素全为1 的下三角矩阵，D 为对角矩阵，LH为L 矩阵的共轭转置，通过矩阵计算，可得计算结果为：

将ljj=1代入，之后进行化解，移项，可得：

为方便计算，引入辅助矩阵G 和C，代入公式可得：

此为改进型Cholesky 分解的数学表达。

1.2 下三角矩阵求逆

根据分解公式，对A=LDLH进行两边求逆，可得A-1=(L-1)HD-1L-1，其中D 的逆即为辅助矩阵C，因此求逆关键在于下三角矩阵L 的求逆[6]，考虑矩阵B 为L 的逆，根据单位矩阵定义，可得：

令P=LB，根据乘法运算：

又因pi,j=0(i ≠j),pii=1，最后将对角线1 代入，解得递推关系：

由此可得求逆的递推公式，根据公式，该计算方法规避了除法计算，在硬件上可节省资源，使运算速度得以提升。

2 模块设计

根据先前的计算公式，本文设计了2 个基本运算模块[7]，一是改进型Cholesky 分解模块，二是下三角矩阵求逆递推模块，最后的三角矩阵相乘则通过矩阵之间的数值关联信息，考虑到最终矩阵也为共轭对称阵，计算时只需计算上三角及对角线即可，本文所使用的数据类型为（1,7,19）的浮点数，因此基本运算模块为ISE 14.7 的浮点数计算IP 核，并且由于是复数，考虑到资源消耗，复数乘法计算使用两个浮点数乘法器，一个加法器一个减法器，消耗两个计算时钟来进行乘法计算，下面介绍子模块的硬件设计图。

图1 为分解模块框图，该模块的主要运算单元为3个复数浮点数乘法器，一组减法器（分别计算实部虚部）和一个浮点数除法器（用于计算倒数），三组RAM用于存储迭代运算过程中的中间变量，RAM_L 存储计算所得的L 矩阵的下三角元素，而RAM_C 则是存储对角矩阵D 的逆矩阵元素。三个Sel 为数据选择器，分别用于选择乘法器、减法器、除法器的输入数据以及将输出数据写入到对应RAM 的地址控制操作，根据所得的迭代顺序，三个乘法器用于并行计算同列标的中间变量值以及L 矩阵元素，减法器则是将计算所得的中间变量进行累减操作，求得下一列的G 矩阵元素，新一列的G 矩阵元素根据行标大小进行有序计算，当求得第一个元素即gj,j后，倒数模块启动，除法器一端输入固定参数(1f80000)16，即十进制1，作为被除数，另一端则输入gj,j进行倒数计算。为了方便地址控制，便于后续更高维矩阵的分解计算，将每次计算的时钟周期设置为35 个，保证浮点数乘除减法模块的运算正确性，以4 维复数矩阵为例，列数为1 的中间变量元素有6 个，3 个复数乘法器用2 个计算时钟即可求得，行标小的元素先进行计算，之后启动累减模块，根据减法次数消耗相应的计算时钟，对角倒数操作则是在减法器求第二个gi,j时并行进行，新的对角倒数求得即可进行新一列的L 矩阵元素计算，此过程只涉及到乘法，此时启动乘法器，根据减法器求得的元素顺序进行下三角及中间变量计算，以此类推，由于只使用一组减法器，时钟消耗主要在于减法计算次数。

图1 分解模块图

分解后L 矩阵输入到下三角求逆模块，根据递推公式，该模块只涉及到乘加计算，因此结构较为简单，只需单个复数乘法器和一组加法器同时计算实部虚部即可，图2 为模块图。

图2 下三角求逆模块图

最后矩阵相乘为基本的点积计算，设计时考虑到数据的关联性[8]，只需计算上三角部分及对角线部分，下三角部分元素只需对对应的上三角部分虚部数据进行取反操作即可，对于本文所使用的浮点数来说，将二进制数据前3 位于100 进行异或计算即可得出相反数，硬件上极易实现：

3 仿真结果分析

本文所用的编译平台为Xilinx 的ISE 14.7，使用的仿真软件为ModelSim 10.1c，分解模块的仿真时序图如图3。

下三角矩阵求逆的时序图如图4。

图3 分解时序图

图4 求逆时序图

最终相乘求逆结果的时序图如图5。

图5 求逆时序图

根据时序图的计算结果，可见其都是按一定的递推顺序求解所得，满足先前所推导的计算公式。将最终计算结果进行数值转换，与MATLAB 端计算结果相比，进行相对误差分析，所设计的模块计算精度在保留5 位小数的情况下，与MATLAB 端计算结果一致，综合相对误差为百分之6.82×10-5，表1 为计算结果对比，考虑输入矩阵维度为4 维方阵，取第一列数据，W 为分解过后求逆过程中最后矩阵相乘的结果，即最终的求逆结果。

表1 误差分析（取第一列数据）

分解模块的逻辑资源消耗情况如如图6。

图6 分解模块逻辑资源

根据时序分析报告，分解模块的运行频率最高可以达到326.535MHz，而下三角矩阵求逆模块的运行频率最高为454.277MHz，图7-图8 为分析时序报告图。

图8 求逆时序报告

根据模块设计的计算时钟数，可计算出最终的硬件模块矩阵求逆运行时间为4.322 μs，本文所使用的PC 平台为CPU Intel i5-7300H，内存为16GB，使用MATLAB 对相同四维复数矩阵进行求逆，即使用inv（）函数，通过tic，toc 测得最短运行时间为37 μs，经计算，硬件模块所消耗的计算时间比MATLAB 端快了8.5 倍，节省了88.32%的计算时间，同时计算精度高，满足一般的工程需求，该模块可以运用到涉及最小二乘法的信号处理硬件算法中，提高效率。

4 结语

本文根据最小二乘类问题的硬件实施上的难点，通过公式简化将伪逆计算难点转化为求解特殊方阵的逆，并选择针对该类共轭对称矩阵分解效果较好，硬件实施上较为合理的改进型Cholesky 法作为分解方法，通过分解完的求逆计算，将下三角矩阵的求逆问题通过数学递推推导转化为只涉及乘加计算的求逆计算，最后的矩阵相乘由于数据关联性及共轭对称的关系，节省了多余的计算步骤，节约了逻辑资源消耗，并根据所推导的数学关系设计了对应的硬件模块，给出了其模块原理图，并通过ModelSim 仿真得出时序图及计算结果，根据时序分析图计算出模块的理论运行时间，两者与MATLAB 计算结果及时间消耗进行比对，可见所设计的硬件模块满足一般工程精度要求，且计算时间也远低于MATLAB 端，能够很好的运用在实时性强的场景中。