基于单边雅克比奇异值分解的FPGA加速实现方法

摘要

本发明公开一种基于单边雅克比奇异值分解的FPGA加速实现方法，该方法首先将输入矩阵平均分为n/2对列向量，并计算每对列向量的范数和内积，然后计算每对列向量的旋转矩阵，并执行正交变换，然后利用round‑robin调度机制将执行正交变换后得到的列向量写入到对应的相邻列向量进行替换，实现每一轮单边Jacobi计算均是在同一份电路上进行循环迭代的效果，简化了数据通道和控制通道的复杂设计，避免了FPGA实现过程中海量级信号布线资源的使用，降低了FPGA资源消耗，提升了电路工作时钟频率，从而显著提升整体性能。

说明书

技术领域

本发明涉及信号处理领域，具体涉及一种基于单边雅克比奇异值分解的FPGA加速实现方法。

背景技术

奇异值分解是线性代数中的一大亮点，在矩阵分解中扮演重要角色，广泛应用于毫米波雷达、无线通信、图像压缩和深度学习等领域。目前已有的研究中，主要以软件程序的方式采用CPU或GPU的实现为主。基于FPGA的Jacobi算法实现矩阵的奇异值分解，虽然结合了FPGA和Jacobi算法两者共同具有的高并行属性优点，但是由于矩阵的奇异值分解涉及行列元素的读写访问、频繁的数据调度替换、大量的循环迭代，造成FPGA实现时仍存在较大的困难。随着输入矩阵尺寸的增大，计算复杂度以O(n

在已有的研究中，基于FPGA的并行双边Jacobi算法主要采用CORDIC模块进行实现，虽然具有占用资源少、计算速度快等优点，但一般只适合矩阵尺寸较小或适中的场景，且要求输入为n行×n列的方阵，数据类型以定点数为主，这对于动态范围和精度要求都有较高需求的场合，则往往无法满足要求。

在已公开的发明中，申请号为CN2019102853514的专利提及了基于FPGA的并行Jacobi计算加速实现方法，主要靠通过优化内部CORDIC电路的计算周期来提高并行流水能力，实现加速，输入的数据类型为定点数。该方法要求输入的矩阵是实对称的协方差矩阵，且对大型尺寸矩阵奇异值分解面临的难题并没有提及。

申请号为CN 2017101340362的专利提及了使用高效的串行控制方案，只使用一个CORDIC模块来实现Jacobi变换，但该方法同样要求输入矩阵是n行×n列的方阵，适用于较小尺寸的矩阵场景需求。

申请号为CN2017107659507的专利侧重于降低运算过程的功从而实现降低成本，将奇异值分解过程中计算密集型的运算搬移到FPGA实现，更多地还需要依赖CPU的调度处理，FPGA在奇异值分解中只是作为加速引擎而已。

并行单边Jacobi算法不需要输入矩阵必须满足方阵的要求，具有更广泛的应用场合。因此，如何设计出一种基于单边Jacobi奇异值分解的FPGA加速实现方法，用以实时、高效地解决矩阵的奇异值分解，尤其是针对大型尺寸矩阵的奇异值分解将很有现实意义。

发明内容

针对现有技术的不足，本发明提出一种基于单边雅克比奇异值分解的FPGA加速实现方法，具体技术方案如下：

一种基于单边Jacobi奇异值分解的FPGA加速实现方法，所述的FPGA有n块BRAM和n/2份计算处理单元；所述方法包括如下步骤：

S1：将输入矩阵的每个列向量一一对应写入到FPGA相应的n块BRAM中；当原有矩阵的列数为奇数时，在原有矩阵的最后补充1列全0元素的列向量，凑成偶数列后作为输入矩阵，输入矩阵为m行×n列，输入数据类型为单精度浮点数；

S2：根据单边Jacobi算法的round-robin调度机制，从输入矩阵的第一个列向量开始，依次将每两个列向量组成一对列向量，从而将输入矩阵平均分为n/2对列向量，并一一对应地存储在FPGA的n/2份计算处理单元的BRAM中；

S3：执行单边Jacobi旋转变换，计算出每对列向量的范数和内积以及Jacobi旋转矩阵，然后按照round-robin调度机制，将正交变换后的列向量写入相应的相邻列向量，覆盖该位置原来的元素；

S4：按照BRAM的固定序号，重复S3，顺序执行单边Jacobi旋转变换，即循环迭代；当n/2对列向量中最大的内积γ的绝对值小于预设的门限时，判定奇异值分解收敛，此时得到各组

为了简化数据通道和控制通道，避免FPGA实现过程中海量级信号布线资源的使用，所述S3具体包括如下子步骤：

S3.1：FPGA的n/2份计算处理单元并行计算每对列向量中的每个列向量的范数α、β，以及它们的内积γ；

S3.2：根据S3.1得到的α、β和γ，FPGA的n/2份计算处理单元并行计算每对列向量的Jacobi旋转矩阵

S3.3：对于第1行元素，FPGA的n/2份计算处理单元并行执行如下操作：

（1）分别从奇数的第i块BRAM和偶数的第i+1块BRAM同时读取第1行元素即u(1,i)和u(1,i+1)，i取值范围是奇数的1,3,5，…，n-1；

（2）执行Jacobi旋转，得到正交变换后的u(1,i)’和u(1,i+1)’；

（3）当i=3,5，…，n-3时，将u(1,i)’写入第i+2块BRAM第1行的位置，即覆盖原来的元素u(1,i+2)；将u(1,i+1)’写入第i-1块BRAM第1行的位置，即覆盖原来的元素u(1,i-1)；

当i=1时，将u(1,i)’写入第i+2块BRAM第1行的位置，即覆盖原来的元素u(1,i+2)；将u(1,i+1)’写入第i块BRAM第1行的位置，即覆盖原来的元素u(1,i)；

当i=n-1时，将u(1,i)’写入第i-1块BRAM第1行的位置，即覆盖原来的元素u(1,i-1)；将u(1,i+1)’写入第i+1块BRAM第1行的位置，即覆盖自身u(1,i+1)；

S3.4：重复S3.3的操作，顺序执行各块BRAM第2行元素，第3行元素，…，直至第m行元素。

为了节省FPGA资源开销，采用分时复用乘法器和加法器执行所述S3中的α、β、γ计算、Jacobi旋转矩阵的计算和正交变换操作。

进一步地，根据如下公式计算所述Jacobi旋转矩阵中的cosθ和sinθ：

其中，sign为符号位，其根据γ与β-α决定，若γ≥0且β-α≥0，或者γ＜0且β-α＜0，则sinθ取正号，反之则取负号。

本发明的有益效果如下：

本发明通过将m行×n列的矩阵平均分为n/2对列向量，在FPGA上实现n/2份PE电路且并行工作，采用保持名义存储列向量序号不变的处理方式，使得每一轮单边Jacobi计算简化为同一份硬件电路针对不同输入的循环迭代，从而简化了数据通道和控制通道的复杂设计，避免了FPGA实现过程中海量级信号布线资源的使用，提高了FPGA布通率和时钟工作频率，从而显著提升整体性能；充分利用范数、内积的计算，Jacobi旋转矩阵的求解以及正交变换操作均有使用浮点乘、加运算且三者间存在串行执行的特点，利用分时复用技术实现乘法器和加法器硬件电路满负荷运行，进一步降低了FPGA计算资源的消耗。

附图说明

图1为512行×128列的矩阵在FPGA内BRAM存储示意图；

图2为列向量对范数和内积计算电路实现图；

图3为列向量对执行正交变换及更新的电路实现图；

图4为512行×128列的矩阵round-robin调度示意图；

图5为利用分时复用技术实现节省计算资源电路图。

具体实施方式

下面根据附图和优选实施例详细描述本发明，本发明的目的和效果将变得更加明白，应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明中使用到一些技术术语，在此给出解释：

（1）FPGA：Field Programmable Gate Array 现场可编程门阵列

（2）BRAM：Block RAM，FPGA内部块状RAM

（3）Jacobi：本发明中特指单边雅克比旋转，常用于基于FPGA的矩阵奇异值分解

（4）round-robin：轮询调度，单边Jacobi旋转奇异值分解常用的一种调度机制；

（5）CORDIC ：Coordinate Rotation Digital Computer，坐标旋转数字计算方法；

（6）PE：Process Element，计算处理单元。

本发明的基于单边雅克比奇异值分解的FPGA加速实现方法，FPGA有n块BRAM和n/2份计算处理单元（以下简称PE电路）；该方法包括如下步骤：

S1：将输入矩阵的每个列向量一一对应写入到FPGA相应的n块BRAM中；其中，当原有矩阵的列数为奇数时，在原有矩阵的最后补充1列全0元素的列向量，凑成偶数列后作为输入矩阵，输入矩阵为m行×n列，输入数据类型为单精度浮点数；写入BRAM中时，第1列的各行写入到第1块BRAM，第2列的各行写入到第2块BRAM，…，第n列的各行写入到第n块BRAM。

S2：根据单边Jacobi算法的round-robin调度机制，从输入矩阵的第一个列向量开始，依次将每两个列向量组成一对列向量，即第1列与第2列组成第1对列向量，第3列与第4列组成第2对列向量，…，第n-1列与第n列组成第n/2对列向量，从而将输入矩阵平均分为n/2对列向量，并一一对应地存储在FPGA的n/2份计算处理单元的BRAM中；

S3：执行单边Jacobi旋转变换，具体操作如下：

S3.1：FPGA的n/2份计算处理单元并行计算每对列向量中的每个列向量的范数α、β，以及它们的内积γ；以第1对列向量为例，分别从第1块BRAM和第2块BRAM同时提取每行元素，分别执行平方或乘法以及累加计算，m行元素经乘累加计算完毕后分别得到α

S3.2：根据S3.1得到的α、β和γ，FPGA的n/2份计算处理单元并行计算每对列向量的Jacobi旋转矩阵

其中，sign为符号位，其根据γ与β-α决定，若γ≥0且β-α≥0，或者γ＜0且β-α＜0，则sinθ取正号，反之则取负号。从而得出n/2个正交变换旋转矩阵，旋转矩阵的角度分别为θ

S3.3：对于第1行元素，FPGA的n/2份计算处理单元并行执行如下操作：

（1）分别从奇数的第i块BRAM和偶数的第i+1块BRAM同时读取第1行元素即u(1,i)和u(1,i+1)，i取值范围是奇数的1,3,5，…，n-1；

（2）执行Jacobi旋转，得到正交变换后的u(1,i)’和u(1,i+1)’；

以第1对列向量的第i行元素为例，根据公式（3）和公式（4）执行正交变换操作：

其中u(j,1)表示原来的第1列第j行元素，u(j,2)表示原来的第2列第j行元素，u(j,1)’表示执行正交变换后得到的新的第1列第j行元素，u(j,2)’表示执行正交变换后得到的新的第2列第j行元素（其中j=1，2，…，m）。

（3）将正交变换后的u(1,i)’写入到相邻奇数列的第i+2列变量的所在的第i+2块BRAM的第1行元素进行覆盖，将正交变换后的u(1,i+1)’写入到相邻偶数列的第i-1列变量的所在的第i-1块BRAM的第1行元素进行覆盖；

特殊情况为：

当上述的奇数的第i列与偶数的第i+1列为第一对列向量即第1列和第2列时，需要特殊处理：执行正交变换后得到的第2列的第1行元素u(1,2)’写入到第1列的第1行元素进行覆盖，即原先覆盖u(1,1)所在BRAM的存储位置；

当上述的奇数的第i列与偶数的第i+1列为最后一对列向量即第n-1列和第n列时，需要特殊处理：执行正交变换后得到的第n-1列的第1行元素u(1,n-1)’写入到偶数的第n-2列的第1行元素进行覆盖，即覆盖原先u(1,n-2)所在BRAM的存储位置；执行正交变换后得到的第n列的第1行元素u(1,n)’覆盖自己原先本身u(1,n)所在BRAM的存储位置。

上述存储过程总结为：

当i=3,5，…，n-3时，将u(1,i)’写入第i+2块BRAM第1行的位置，即覆盖原来的元素u(1,i+2)；将u(1,i+1)’写入第i-1块BRAM第1行的位置，即覆盖原来的元素u(1,i-1)；

当i=1时，将u(1,i)’写入第i+2块BRAM第1行的位置，即覆盖原来的元素u(1,i+2)；将u(1,i+1)’写入第i块BRAM第1行的位置，即覆盖原来的元素u(1,i)；

当i=n-1时，将u(1,i)’写入第i-1块BRAM第1行的位置，即覆盖原来的元素u(1,i-1)；将u(1,i+1)’写入第i+1块BRAM第1行的位置，即覆盖自身u(1,i+1)；

S3.4：重复S3.3，顺序执行各块BRAM第2行元素，第3行元素，…，直至第m行元素；

S4：重复S3，根据BRAM的固定序号顺序执行多轮单边Jacobi旋转变换，即循环迭代；

当n/2对列向量中最大的内积γ的绝对值小于预设的门限时，判定奇异值分解收敛，此时得到各组

实际上，经过S3的相邻奇偶列向量更新操作后，此时各块BRAM实际存储的列向量序号是变化的，即遵循round-robin调度机制的更新规律，为简化FPGA电路实现时数据通道和控制通道设计，采用保持名义存储列向量序号不变（即按照BRAM的固定序号）的处理方式，即第1块BRAM始终存储的是名义上的第1列列向量，第2块BRAM始终存储的是名义上的第2列列向量，…，第n块BRAM始终存储的是名义上的第n列列向量，进行循环迭代，而不区分各块BRAM实际存储的列向量序号。即：仍旧保持第1块BRAM存储的列向量与第2块BRAM存储的列向量进行单边Jacobi旋转变换，第3块BRAM存储的列向量与第4块BRAM存储的列向量进行单边Jacobi旋转变换，…，第n-1块BRAM存储的列向量与第n块BRAM存储的列向量进行单边Jacobi旋转变换。实现每一轮单边Jacobi计算均是在同一份电路上进行循环迭代的效果，达到简化数据通道和控制通道复杂设计的目的，降低FPGA资源消耗同时提升电路工作时钟频率。

因S3中的α、β、γ计算、Jacobi旋转矩阵的计算和正交变换操作，均会使用浮点数乘法和加法操作，且这些计算存在串行执行特点，采用分时复用乘法器和加法器，以节省FPGA资源开销。

下面以512行×128列的矩阵作为输入矩阵，对本发明的方法进一步举例说明。

该矩阵中的每个元素的数据类型为符合IEEE754标准的单精度浮点数，位宽32比特，FPGA型号选择的是Xilinx公司的XC7VX690T-2FFG1761，FPGA内部每一块BRAM基本存储容量为18Kb（其中2Kb用于奇偶校验），可以配置为深度512、位宽32，刚好实现紧凑利用，512行×128列单精度浮点数矩阵使用128块BRAM即可满足要求。根据当前实施例具体应用，设置内积收敛门限为γ

如图1所示，为本发明512行×128列矩阵在FPGA内BRAM存储示意图，具体的执行过程如下：

步骤一：每个列向量分别写入到FPGA相应的BRAM中，即第1列的每一行按顺序写入到第1块BRAM，第2列的每一行按顺序写入到第2块BRAM，…，第128列的每一行按顺序写入到第128块BRAM。

步骤二：根据单边Jacobi算法的round-robin调度机制将其平均分为64对列向量，具体为第1列与第2列组成第1对列向量，第3列与第4列组成第2对列向量，…，奇数第i列与偶数第i+1列组成第（i+1）/2对列向量，第127列与第128列组成第64对列向量。

步骤三：FPGA上的64份PE电路并行计算每对列向量各自的范数及内积即α、β、γ：以第1对列向量为例，在时钟节拍驱动下，同时从第1块BRAM和第2块BRAM分别提取每行32比特的单精度浮点数据，分别执行平方或乘法以及累加计算，经过512拍的时钟之后，列向量的所有元素输出，并且经乘累加计算完毕后分别得到α

根据S1得到的64组的α、β、γ，FPGA上的64份PE电路并行计算每对列向量的Jacobi旋转矩阵

FPGA上的64份PE电路并行执行正交变换，具体操作如下：

（1）分别从奇数的第i块BRAM和偶数的第i+1块BRAM同时读取第1行元素即u(1,i)和u(1,i+1)，i取值范围是奇数的1,3,5，…，127；

（2）执行Jacobi旋转，得到正交变换后的u(1,i)’和u(1,i+1)’；

以第2对列向量即第3列列向量和第4列列向量的第i行元素为例执行正交变换操作如下：

其中u(j,3)表示原来的第3列第j行元素，u(j,4)表示原来的第4列第j行元素，u(j,3)’表示执行正交变换后得到的新的第3列第j行元素，u(j,4)’表示执行正交变换后得到的新的第4列第j行元素（其中j=1，2，…，512），其FPGA实现电路如图3所示；

（3）每一轮的单边Jacobi旋转变换执行如图4所示的调度机制：

第1列列向量执行一次单边Jacobi旋转变换后，向右移动替代原先第3列列向量所在位置，同理，第3列列向量执行一次单边Jacobi旋转变换后，向右移动替代原先第5列列向量所在位置，依次类推处理，直到第127列列向量在执行一次单边Jacobi旋转变换后替换原先的第126列列向量所在位置；第126列列向量在执行一次单边Jacobi旋转变换后替换原先的第124列列向量所在位置，第124列列向量在执行一次单边Jacobi旋转变换后替换原先的第122列列向量所在位置，依次类推处理，直到第2列列向量在执行一次单边Jacobi旋转变换后替换原先的第1列列向量所在位置；第128列列向量执行一次单边Jacobi旋转变换后替换本身第128列列向量所在位置，即保持原先位置不变。

64份PE电路并行执行上述的（1）~（3）操作，对列向量正交变换以及更新操作在硬件上均是同时并行工作，不会造成数据破坏，同时节省不必要的中间结果缓存。

（4）同理，接下来依次顺序执行各块BRAM第2行元素，第3行元素，…，直至第512行元素。

经过步骤（1）~（4）相邻奇偶列向量更新操作后，此时各块BRAM实际存储的列向量序号是变化的，即遵循round-robin调度机制的更新规律，为简化FPGA电路实现时数据通道和控制通道设计，采用保持名义存储列向量序号不变的处理方式，即第1块BRAM始终存储的是名义上的第1列列向量，第2块BRAM始终存储的是名义上的第2列列向量，…，第128块BRAM始终存储的是名义上的第128列列向量。步骤四：重复步骤三，根据BRAM的固定序号顺序执行单边Jacobi旋转变换，即循环迭代；根据预先设置的门限γ

由于64对列向量正交变换以及更新操作在硬件上同时并行进行，同时并行提取原有元素数据和同时并行替换更新相邻列向量数据，不会造成数据破坏，且不必额外增加数据缓存开销，从而实现FPGA硬件资源节省；图3以第3列和第4列列向量为例描述了正交变换及其更新到相邻列向量的电路实现图。

为了兼顾宽动态范围和高精度要求，在实现时采用了单精度浮点数，相对于定点数实现而言，增加了很多的FPGA逻辑资源的消耗，尤其是浮点数乘法器和加法器，会在设计中大量被使用。在本发明的实施例中，范数、内积的计算，Jacobi旋转矩阵

FPGA综合结果显示，512行×128列的单精度浮点矩阵在XC7VX690T-2FFG1761上使用了128块18Kb Block RAM、1600 DSP48、212K LUT，在250MHz时钟运行下实现在3.65ms内快速完成奇异值分解。

本领域普通技术人员可以理解，以上所述仅为发明的优选实例而已，并不用于限制发明，尽管参照前述实例对发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在发明的精神和原则之内，所做的修改、等同替换等均应包含在发明的保护范围之内。