一种矩阵奇异值分解的FPGA加速实现方法与流程

一种矩阵奇异值分解的fpga加速实现方法
技术领域
1.本发明涉及信号处理领域，具体涉及一种矩阵奇异值分解的fpga加速实现方法。


背景技术：

2.奇异值分解是线性代数中一种重要的矩阵分解，广泛应用于信号处理、图像压缩和深度学习等领域。目前已有的研究中，主要以软件程序的方式采用cpu或gpu进行实现，近年来随着fpga技术迅猛发展，采用fpga实现矩阵奇异值分解逐渐成为一项热门技术，尤其是一些已经部署了fpga的应用场景，基于fpga实现矩阵奇异值分解来替代gpu方案可以达到减少成本、降低功耗，相比cpu方案则可以获得更加实时、低延迟的性能。
3.但是，奇异值分解涉及大量的数学运算和循环迭代，特别是对于大型尺寸矩阵，具有计算密集型和存储密集型的鲜明特征，对fpga内部的计算资源、逻辑资源、bram容量及片内外传输带宽都提出了苛刻需求，同时也导致了fpga开发的异常困难和巨大的工作量。在已公开的研究和发明中，受制于fpga资源有限以及奇异值分解本身高复杂性等原因，目前基于fpga实现的奇异值分解的矩阵尺寸普遍偏小、实时性差、只能支持固定尺寸大小的矩阵输入等不足。


技术实现要素：

4.针对现有技术的不足，本发明提出一种矩阵奇异值分解的fpga加速实现方法，基于单边jacobi算法，采用分而治之的策略，对输入的大尺寸矩阵以列向量形式进行分块，将原本异常复杂的大尺寸矩阵的奇异值分解，转化成各种相对易于实现的小尺寸矩阵分解，使得fpga开发难度得到显著下降，对fpga资源的苛刻需求明显降低。
5.本发明的目的通过如下的技术方案来实现：一种矩阵奇异值分解的fpga加速实现方法，所述fpga有3k个bram，所述矩阵为m行
×
n列，所述方法包括如下步骤：s1：将所述矩阵按每k列列向量为一组，平均分成p=n/k个子块，若p=n/k不能整除时，则预先对所述矩阵末尾进行列向量补齐以达到整除，且新添加的列向量所有元素值取0；s2：将第1子块与第2子块进行组合得到m行
×
2k列的新矩阵，并将每列列向量一一对应写入到fpga相应的2k个bram中，每个bram对应其中一列列向量；s3：对m行
×
2k列的新矩阵执行单边jacobi旋转变换，以round
‑
robin调度机制共执行2k
‑
1轮，同时将第3子块的k列列向量预先写入到剩余的k个bram；s4：将s3中对应第1子块的k列列向量中间计算结果与预先写入的第3子块进行组合，得到新的m行
×
2k列矩阵，执行2k
‑
1轮的单边jacobi旋转变换；同时将s3中对应第2子块的k列列向量的中间结果a'2写回片外dram，等到这部分的k个bram的空间释放后，再将第4子块的k列列向量写入这部分的k个bram中；以此类推，按照如下子块间的组合规则和顺序(1,2)
→
(1,3)
→
(1,4)
→
...
ꢀ→
(1,
p)；(2,p)
→
(3,p)
→
(4,p)
→
...
ꢀ→
(p
‑
1,p)；(p
‑
1,2)
→
(p
‑
1,3)
→
(p
‑
1,4)
→
...
ꢀ→
(p
‑
1,p
‑
2)；...；(3,2)；p个子块两两组合，共计种情况，完成整个矩阵一整轮的单边jacobi旋转变换；s5：将每个子块经过单边jacobi旋转变换后的中间结果作为新一整轮的迭代输入，重复s2~s4相同的组合和计算操作，直至满足收敛条件，即m行
×
n列矩阵奇异值分解完毕。
6.进一步地，所述s4中，完成整个矩阵最后一个子块组合(3,2)的2k
‑
1轮的单边jacobi旋转变换后，将第3个子块的中间结果写回片外dram，将第2个子块的中间结果保留在片内k个bram中，即复用第2个子块的中间计算结果。
7.本发明的有益效果如下：本发明的方法特别适用于大尺寸矩阵，通过将m行
×
n列的大尺寸矩阵以每k列列向量为一组平均分成p=n/k个子块，各子块间通过两两组合形成m行
×
2k列的小尺寸矩阵，通过该分而治之策略使得原本异常复杂的奇异值分解得到显著简化，并可支持任意尺寸矩阵奇异值分解，子块间组合规则简单清晰，易于fpga开发实现；fpga内部的3组bram，当任意两组bram用于jacobi旋转变换操作时，剩余的一组bram用于与片外dram的读写缓存，形成乒乓结构工作机制，可以提高整体电路的并行效率，提高流水线工作能力；子块组合内以round
‑
robin方式执行单边jacobi旋转变换，提高了数据复用率，子块组合间以交替组合的方式则更进一步地提高了数据复用，减少了中间计算结果的搬运量，降低了片内外高带宽传输需求。
附图说明
8.图1为1024行
×
1024列的矩阵平均分为16个子块及组合示意图；图2为第1、2子块组成的矩阵2k
‑
1轮单边jacobi变换示意图；图3为第1、2、3子块与片外dram关系图；图4为第1、2子块中间计算结果、第3子块及片外dram关系图；图5为jacobi变换、bram及片外dram乒乓工作机制示意图。
具体实施方式
9.下面根据附图和优选实施例详细描述本发明，本发明的目的和效果将变得更加明白，应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
10.首先给出技术术语解释：（1）fpga：field programmable gate array 现场可编程门阵列（2）bram：block ram，fpga内部块状ram（3）jacobi：本发明中特指单边雅克比旋转，常用于基于fpga的矩阵奇异值分解（4）round
‑
robin：轮询调度，单边jacobi旋转奇异值分解常用的一种调度机制（5）dram：dynamic random access memory，这里特指片外dram， 如ddr3（或ddr4）sdram。
11.本发明的矩阵奇异值分解的fpga加速实现方法，要求fpga有3k个bram，定义待分解的矩阵为m行
×
n列，则该方法包括如下步骤：
s1：将所述矩阵按每k列列向量为一组，平均分成p=n/k个子块，若p=n/k不能整除时，则预先对所述矩阵末尾进行列向量补齐以达到整除，且新添加的列向量所有元素值取0。
12.s2：将第1子块与第2子块进行组合得到m行
×
2k列的新矩阵，并将每列列向量一一对应写入到fpga相应的2k个bram中，每个bram对应其中一列列向量。
13.s3：对m行
×
2k列的新矩阵执行单边jacobi旋转变换，以round
‑
robin调度机制共执行2k
‑
1轮，同时将第3子块的k列列向量预先写入到剩余的k个bram。
14.s4：将s3中对应第1子块的k列列向量中间计算结果与预先写入的第3子块进行组合，得到新的m行
×
2k列矩阵，执行2k
‑
1轮的单边jacobi旋转变换；同时将s3中对应第2子块的k列列向量的中间结果a'2写回片外dram，等到这部分的k个bram的空间释放后，再将第4子块的k列列向量写入这部分的k个bram中；以此类推，按照如下的组合规则和顺序(1,2)
→
(1,3)
→
(1,4)
→
...
ꢀ→
(1,p)；(2,p)
→
(3,p)
→
(4,p)
→
...
ꢀ→
(p
‑
1,p)；(p
‑
1,2)
→
(p
‑
1,3)
→
(p
‑
1,4)
→
...
ꢀ→
(p
‑
1,p
‑
2)；...；(3,2)；p个子块两两组合，共计种情况，完成整个矩阵一整轮的单边jacobi旋转变换。
15.定义第i个子块为a
i
，a
i
与其他子块执行完jacobi旋转变换后的k列列向量的中间计算结果为a'
i
，该组合和写入、写出规则用一种更为通用的表述方式可以表述为：a'1与预先写入的a
i
（i≥3）组成新的m行
×
2k列矩阵，执行2k
‑
1轮的单边jacobi旋转变换；同时将前一子块中间计算结果a'
i
‑1写回片外dram，等到这部分的k个bram的空间释放后，再将后一子块a
i 1
写入这部分的k个bram中；直到完成a'1和a
p
的2k
‑
1轮的单边jacobi旋转变换，至此子块a1与其他子块的组合和计算全部遍历一次，保持a'
p
在fpga内部bram不变；并将中间计算结果a'1写回片外dram，等到这部分的k个bram的空间释放后，将a'2从片外dram写入这部分的k个bram中，与前述a'
p
组合，重复上述过程，完成a'
p
与其他子块的组合和计算；依次类推，每次保留一个子块的中间结果在片内，依次写入未与该子块组合和计算过的子块，执行2k
‑
1轮的单边jacobi旋转变换；直至实现p个子块间两两组合，共计种情况，完成整个矩阵一整轮的单边jacobi旋转变换计算。
16.s5：将每个子块经过单边jacobi旋转变换后的中间结果作为新一整轮的迭代输入，重复s2~s4相同的组合和计算操作，直至满足收敛条件，即m行
×
n列矩阵奇异值分解完毕。
17.另外，在完成s4后，可以将最后得到的第3个子块和第2个子块的k列列向量的中间计算结果均写回片外dram，然后s5的下一轮重新开始写入、组合和计算操作。但为了进一步提高数据复用，减少中间计算结果的搬运量，此时也可以仅将第3个子块的k列列向量的中间计算结果写回片外dram，保留第2个子块的k列列向量的中间计算结果仍在片内的k个bram中。这样开始新的一整轮的单边jacobi旋转变换时，就不用重复把第2个子块的k列列向量的中间计算结果导出导入，即第2个子块的中间计算结果得到复用，省去该中间结果在fpga内部的bram与片外dram间的来回搬运操作，减少数据搬移次数。新的一整轮的单边jacobi旋转变换时，刚开始的组合依次为(2,1)、(1,3)、(1,4) ...
ꢀ→
(1,p)。
18.下面以一个具体的实施例对本发明的方法予以解释和说明。
19.该具体实施例以1024行
×
1024列矩阵奇异值分解进行阐述，矩阵元素数据类型为单精度浮点数即位宽32比特，选择xilinx公司提供的vc707作为开发板，fpga型号是xc7vx485t
‑
2ffg1761c，内含1015块36kb的bram，每块bram的4kb用于奇偶校验实际可用为32kb，如果不考虑分块处理方法，仅将整个矩阵数据存储到内部bram就需1024块36kb的bram，超出了该fpga所能提供的内部bram资源。为解决该问题，本实施例中，以64列列向量作为一个子块，将原始输入矩阵分为p=1024/64即16个子块，任意两个子块进行组合形成1024行
×
128列的新矩阵。每一列列向量由1024个单精度浮点数组成刚好占用1块bram，64列列向量作为一组共占用64块bram，本实施例中一共使用3组bram，共192块bram。如图1所示，为本发明1024行
×
1024列矩阵平均分成16个子块及其两两交替组合示意图，图1中的m=1024，第1列置第64列为第1子块，第65列至第128列为第2子块，依次类推，第961列至第1024列为第16子块；其组合规则和顺序如下：(1,2)
→
(1,3)
→
(1,4)
→
...
ꢀ→
(1,16)；(2,16)
→
(3,16)
→
(4,16)
→
...
ꢀ→
(15,16)；(15,2)
→
(15,3)
→
(15,4)
→
...
ꢀ→
依次类推直至(3,2)，共计=120种组合情况。
20.本实施例具体的执行过程如下：步骤1：将第1个子块和第2个子块进行组合得到1024行
×
128列新矩阵，并将第1个子块即第1列至第64列列向量从片外dram写入到第1组bram，第2子块即第65列至第128列列向量从片外dram写入到第2组bram，每列列向量对应一块bram。
21.步骤2：对步骤1的1024行
×
128列新矩阵，执行单边jacobi旋转变换，如图2所示，图中的k=64，以round
‑
robin调度机制共执行127轮，实现子块组合内数据的高度复用，且在执行完127轮后128列列向量的顺序与执行前完全一致，即第1列至第64列列向量仍旧更新写回到第1组bram的第1块至第64块bram中，第65列至第128列列向量仍旧更新写回到第2组bram的第1块至第64块bram中；在执行单边jacobi旋转变换的同时，预先将下一组合即第3子块的第129列至第192列列向量写入到第3组的第1块至第64块bram中，从而达到了乒乓结构的工作机制，可以提高整体电路的并行流水工作性能，如图3所示，图中的k=64。
22.步骤3：将步骤2中第1子块的64列列向量的中间计算结果即第1组bram与第3组bram进行组合形成新的1024行
×
128列新矩阵，同理执行步骤2中的127轮单边jacobi旋转变换；与此同时，将步骤2中第2组bram的中间计算结果写回片外dram进行保存，达到乒乓结构工作机制的效果；具体电路如图4所示，图中的k=64，矩阵元素a'
1,1
表示原来第1子块第1行第1列元素a
1,1
执行了jacobi变换后的中间结果，a'
2,k 1
表示第2子块原来第2行第1列元素a
2,k 1
执行了jacobi变换后的中间结果，其余元素同理不再赘述，图5给出了整体电路中jacobi变换、bram及片外dram乒乓工作机制。
23.步骤4：依次类推，直至第1子块的中间结果与第16个子块组合并执行单边jacobi旋转变换完毕，此时需要将第1子块列向量中间结果即第1组bram数据写回片外dram保存，第16子块的列向量中间结果与下一组合即第2子块进行组合，即从子块组合(1,16)变为子块组合(2,16)继续执行类似步骤2和步骤3中的操作，子块16的中间计算结果得到复用，省去该中间结果传片外dram操作，减少数据搬移次数。
24.步骤5：依次类推，直至16个子块两两交替组合共计120中情况，每个组合分别执行上述同类操作，直至完成整个大尺寸矩阵一整轮的单边jacobi旋转变换，所有的子块顺序
与最初始的大尺寸矩阵顺序完全一致，并且所有的中间计算结果均写回片外dram进行暂存，上述步骤1至步骤5操作，称作一次sweep；步骤6：重复上述所有步骤，对于1024行
×
1024列的矩阵，执行8次sweep可以认定达到收敛条件，从而实现1024行
×
1024列大尺寸矩阵奇异值分解。
25.fpga综合结果显示，1024行
×
1024列的单精度浮点矩阵在xc7vx485t
‑
2ffg1761c上使用了384块bram、253k
‑
lut，在200mhz时钟运行下实现在0.690秒内快速完成奇异值分解。相比于xilinx公司公布的矩阵奇异值分解solver库中的结果，实现512行
×
512列实对称单精度浮点矩阵奇异值分解，在alveo u250加速卡上使用了128块uram 307块bram，65k
‑
lut，用时1.687秒，根据矩阵尺寸规模折算，本发明实施例可以实现接近20倍的实时性能提升。
26.通过本发明实施例可以发现，对于m行
×
n列大尺寸矩阵，列数量的增加对逻辑、存储资源消耗几乎没有明显增加，只是增加了两两子块组合情况，而对于不同的行数而言，只需对bram的深度进行加深，就可适应不同行数的矩阵分解。因此，本发明可以实现任意尺寸矩阵奇异值分解的效果。
27.本领域普通技术人员可以理解，以上所述仅为发明的优选实例而已，并不用于限制发明，尽管参照前述实例对发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在发明的精神和原则之内，所做的修改、等同替换等均应包含在发明的保护范围之内。