基于矩阵三角化分解的Cholesky分解及FPGA并行结构设计
刘书勇 , 林俊宇 , 吴艳霞 , 张博为     
哈尔滨工程大学 计算机科学与技术学院, 哈尔滨 150001
摘要:矩阵运算是高性能计算中核心问题之一,矩阵分解是提高矩阵运算并行性的重要途径,飞速发展的FPGA为并行运算结构提供了有力的环境支持。该文基于子矩阵更新同一化算法实现了Cholesky分解,基于FPGA设计了相应的并行结构。实验结果表明:与通用处理器的软件实现相比,本文实现的Cholesky分解的FPGA并行结果在核心计算性能上可以取得10倍以上的加速比,该算法针对矩阵三角化计算过程具有更高的数据和流水并行性。
关键词矩阵三角化分解     Cholesky分解     并行结构     现场可编程门阵列    
Cholesky decomposition and parallel structure design based on matrix triangularization decomposition
LIU Shuyong , LIN Junyu , WU Yanxia , ZHANG Bowei     
Collage of Computer Science and Technology, Harbin Engineering University, Harbin 150001, China
Abstract:Matrix computing is one of the core problems in high performance computing with matrix decomposition being an important way to improve the parallelism of matrix computations. FPGA gives a powerful environment for parallel computing. This study uses Cholesky decomposition based on a hardware-adaptive parallel sub-matrix identity updating algorithm. Its parallel structure is based on FPGA. Tests show that this structure achieves more than 10 fold speedup compared to general-purpose processors in terms of the kernel computational speed because the algorithm has better data-parallelism and pipeline-parallelism during matrix triangularization.
Key words: matrix triangularization decomposition     Cholesky decomposition     parallel structure     field programmable gate array    

当前,对高性能矩阵三角化分解的研究主要从通用计算[1]、脉动阵列[2]和FPGA实现[3]3个方向开展。在高性能通用计算领域内,从计算特征及应用领域考虑,对矩阵三角化分解的研究主要在基于单指令流多数据流(SIMD)或多指令流多数据流(MIMD)的向量机、共享存储的多处理机等技术[4-10]的并行计算机体系结构上开展。一些研究项目在单处理器以及通用计算图形处理器(GPGPU)上也有良好的表现。该领域的研究通过向程序设计者提供具有标准接口和较好移植性的数值计算模块,以在目标体系结构上实现算法的高性能计算。许多高层编程语义数值计算库如基于C++的cvmlib、英特尔数学核心函数库(Intel MKL)和矩阵计算模板库(matrix template library,MTL)等则通过调用1个或若干基础数值计算子程序库来实现。在业界已有的研究成果中,以BLAS、LAPACK和LINPACK等[11-12]为代表的一系列数值线性代数软件包最具有代表性。使用线性或二维脉动阵列机实现三角化分解的研究集中在20世纪80年代和90年代初,往往针对某些特定科学工程领域的应用。该结构是VLSI设计中一种经典体系结构,通过消除阵列中的广播路径实现三角化分解中流水并行的最大化,往往占用较大的片上面积,产生较大的功耗。近年来,更多的研究者在固定面积的FPGA上设计具有较高性能的矩阵三角化分解并行结构。归纳起来,这些设计分别选择从单个处理单元(processing element,PE)设计、存储层次设计、阵列结构设计等方面开展优化,实现FPGA矩阵三角化分解的研究以及针对软件数值线性代数软件包的FPGA设计[13-15]。从实现结构角度来说,它们大多是针对某种具体矩阵分解的一维线性阵列及二维多边形阵列以不同并行算法加以实现[16]。这些设计在科学及工程应用领域已经得到广泛应用。但是,对FPGA矩阵三角化分解的研究仍有可提升的空间。从矩阵计算类属的角度考虑,专用的矩阵分解并行算法及阵列结构很难直接复用到其他矩阵分解的实现中,而以通用化为设计目标的结构则在计算性能上无法与专用阵列结构相比。此外,已有的一些矩阵分解细粒度并行算法往往需要对算法本身进行较大的改动,通过增加硬件操作原语来实现算法指令级的并行[17-18]。在设计一类具有共同计算特征的矩阵分解应用中,对计算过程和算法本身细粒度并行的分析和优化成为提高计算阵列性能的关键。矩阵三角化分解过程中均含有三角化计算,这一数值计算过程同样具有规律性较强的数据排列与流动特征,非常适合进行较细粒度的FPGA并行优化设计。

本文针对矩阵三角化分解中共有的三角化计算迭代过程,通过分析其中子矩阵更新过程的线性规律,提出一种具有一般性的子矩阵更新同一化并行算法,进而提出基于该算法的计算复杂度较低的矩阵三角化分解并行结构,针对Cholesky矩阵三角化分解在并行结构上的高性能实现及优化方法开展研究。

1 研究背景

本文提到的矩阵和向量均在实数范围内,矩阵皆为方阵。

单个矩阵的三角化分解将矩阵A分解为正交矩阵与三角矩阵之积,或分解为一个上三角矩阵与一个下三角矩阵之积。三角化过程可以总结为: 给定n×p(pn)维的矩阵X,根据式(1)求n×n维的非奇异矩阵M以及p×p维的上三角矩阵R

$MX=\left[ \begin{align} & R \\ & 0 \\ \end{align} \right].$ (1)

关于子矩阵更新同一化算法请参考文[19]

2 Cholesky分解及FPGA并行结构设计

Cholesky分解将n×n维对称正定矩阵A(其中元素为aij,1≤i,jn)分解成A=LLT的形式,L(其中元素为lij,1≤jin)为对角线元素为正数的下三角矩阵。LLT如式(2)所示。

$\begin{align} & L=\left[ \begin{matrix} {{l}_{11}} & 0 & \cdots & 0 & \cdots & 0 \\ {{l}_{21}} & {{l}_{22}} & \cdots & 0 & \cdots & 0 \\ \vdots & \vdots & {} & \vdots & {} & \vdots \\ {{l}_{k1}} & {{l}_{k2}} & \cdots & {{l}_{kk}} & \cdots & {} \\ \vdots & \vdots & {} & \vdots & {} & \vdots \\ {{l}_{n1}} & {{l}_{n2}} & \cdots & {{l}_{nk}} & \cdots & {{l}_{nn}} \\ \end{matrix} \right], \\ & {{L}^{T}}=\left[ \begin{matrix} {{l}_{11}} & {{l}_{21}} & \cdots & {{l}_{k1}} & \cdots & {{l}_{n1}} \\ 0 & {{l}_{22}} & \cdots & {{l}_{k2}} & \cdots & {{l}_{n2}} \\ \vdots & \vdots & {} & \vdots & {} & \vdots \\ 0 & 0 & \cdots & {{l}_{kk}} & \cdots & {{l}_{nk}} \\ \vdots & \vdots & {} & \vdots & {} & \vdots \\ 0 & 0 & \cdots & 0 & \cdots & {{l}_{nn}} \\ \end{matrix} \right]. \\ \end{align}$ (2)

L可以通过式(3)和(4)求得。

${{l}_{ij}}=\left\{ \begin{matrix} {{a}_{ij}}/{{l}_{jj}}, & italic>j=1; \\ \left( {{a}_{ij}}-\sum\limits_{k=1}^{j-1}{{{l}_{ik}}{{l}_{jk}}} \right)/{{l}_{jj}}, & italic>j>1; \\ \end{matrix} \right.~$ (3)
${{l}_{ii}}=\sqrt{{{a}_{ii}}-\sum\limits_{k=1}^{j-1}{l_{ik}^{2}}}.$ (4)

将 Cholesky分解的执行过程由计算L的过程转化为计算LT的等价过程(简称LT-SC)。 从矩阵几何特征的角度,LT-SC将更新矩阵求解下三角矩阵过程变为求解上三角矩阵过程,如图 1所示,其中“×”为无关项。

图 1 LT-SC的计算过程

因此,LT-SC基于列k的1次子矩阵更新亦可以表示为式(5)和(6)。

$\begin{align} & {{L}^{T}}_{(1,1)}=\left[ \begin{matrix} {{u}_{11}} & {{u}_{12}} & {{u}_{13}} & {{u}_{14}} \\ \times & a_{22}^{'} & a_{23}^{'} & a_{24}^{'} \\ \times & \times & a_{33}^{'} & a_{34}^{'} \\ \times & \times & \times & a_{44}^{'} \\ \end{matrix} \right] \\ & ,{{\left[ \begin{matrix} {{u}_{11}} \\ {{u}_{12}} \\ {{u}_{13}} \\ {{u}_{14}} \\ \end{matrix} \right]}^{T}}=\left[ \begin{matrix} \sqrt{{{a}_{11}}} \\ {{a}_{12}}/\sqrt{{{a}_{11}}} \\ {{a}_{13}}/\sqrt{{{a}_{11}}} \\ {{a}_{14}}/\sqrt{{{a}_{11}}} \\ \end{matrix} \right], \\ & {{A}^{'}}=\left[ \begin{matrix} {{a}_{11}} & {{a}_{12}} & {{a}_{13}} & {{a}_{14}} \\ \times & a_{22}^{'} & a_{23}^{'} & a_{24}^{'} \\ \times & \times & a_{33}^{'} & a_{34}^{'} \\ \times & \times & \times & a_{44}^{'} \\ \end{matrix} \right]= \\ & {{T}_{1,1}}\left[ \begin{matrix} {{a}_{11}} & {{a}_{12}} & {{a}_{13}} & {{a}_{14}} \\ {{a}_{21}} & {{a}_{22}} & {{a}_{23}} & {{a}_{24}} \\ {{a}_{31}} & \times & {{a}_{33}} & {{a}_{34}} \\ {{a}_{41}} & \times & \times & {{a}_{44}} \\ \end{matrix} \right] \\ \end{align}$ (5)
${T_{1,1}} = \left[ {\matrix{ 1 & 0 & 0 & 0 \cr { - {a_{21}}/{a_{11}}} & 1 & 0 & 0 \cr { - {a_{31}}/{a_{11}}} & 0 & 1 & 0 \cr { - {a_{41}}/{a_{11}}} & 0 & 0 & 1 \cr } } \right].$ (6)

引入置换矩阵P,令:

$P=\left[ \begin{matrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\ \end{matrix} \right]$ (7)

LT-SC的求解过程亦可以表示为式(8)的计算过程。LT-SC每一步骤的计算对象均为上三角矩阵,因此每完成1次矩阵更新,有效上三角矩阵的维数减1。

$\left\{ \begin{align} & {{A}^{(t)}}=P({{T}_{1,1}}{{A}^{(t-1)}})P,0 <t\le n. \\ & {{A}^{(0)}}=A \\ \end{align} \right.$ (8)

Z 看成更新后的Azij为更新后的aij。设:

$z_{~i-1,j-1}^{(t)}={{g}^{(t)}}\left( i,j,k \right)=z_{i,j}^{(t-1)}-\left( \text{ }z_{i,1}^{(t-1)}/\text{ }z_{1,1}^{(t-1)} \right)z_{1,j}^{(t-1)},$

1 <i,jn,1≤kn且0 <tn,Z(0)=Z

则使用子矩阵同一化并行算法后,子矩阵的更新过程可以表示为式(9)。

${T_d}{A_{4 \times 4}} = {T_d}\left[ {\matrix{ {{a_{11}}} & {{a_{12}}} & {{a_{13}}} & {{a_{14}}} \cr \times & {{a_{22}}} & {{a_{23}}} & {{a_{24}}} \cr \times & \times & {{a_{33}}} & {{a_{34}}} \cr \times & \times & \times & {{a_{44}}} \cr } } \right],$ (9a)
$\eqalign{ & {T_d} = \left[ {\matrix{ { - {a_{12}}/{a_{11}}} & 1 & 0 & 0 \cr { - {a_{13}}/{a_{11}}} & 0 & 1 & 1 \cr { - {a_{14}}/{a_{11}}} & 0 & 1 & 1 \cr 1 & 0 & 1 & 1 \cr } } \right] = \cr & p\left[ {\matrix{ 1 & 0 & 0 & 0 \cr { - {a_{12}}/{a_{11}}} & 1 & 0 & 0 \cr { - {a_{13}}/{a_{11}}} & 0 & 1 & 0 \cr { - {a_{14}}/{a_{11}}} & 0 & 0 & 1 \cr } } \right]. \cr} $ (9b)

LT-SC的计算依赖图如图 2所示。在该依赖图中,以k轴的正方向为时间维度,有: ai,j(t)依赖于ai+1,j+1(t),0 <jin。即只更新除去右端边界列的子三角矩阵,且a1,1(t)始终为下一次子矩阵更新计算的主元。根据LT-SC的计算依赖图构建的阵列结构与时空图如图 3所示。将子矩阵更新趟间依赖映射为连结斜上方邻接PE单元的输出边。为首行PE单元(a1,2,a1,3,…,a1,n)增加旁路边(ai,j=aj,i),以获取子矩阵行更新行变换关键向量计算的数据并行。对应的LT-SC并行结构PE的配置和LT-SC总体并行结构分别如图 45所示。

图 2 LT-SC的计算依赖图

图 3 LT-SC的阵列结构与计算时空图

图 4 LT-SC并行结构PE的配置

图 5 LT-SC总体并行结构

3 实验结果

本文提出的设计采用VHDL语言实现。使用Mentor Graphic Modelsim v6.5仿真验证工具进行仿真后在Xilinx ISE 11.1工具中完成综合,并给出资源占用结果; 对设计的性能进行估计与分析; 最后将设计加载到FPGA上进行实际运行与功能验证。所使用的FPGA分别为xc5vlx330t和xc5vlx110t。

3.1 FPGA资源占用

本文涉及的浮点计算如加法、减法、乘法、除法以及开平方等均使用开源算术浮点库libhdlfltp的单精度流水浮点计算部件实现,遵照IEEE754标准的浮点数格式。

表 1给出LT-SC并行结构中PE单元的综合结果,包括BcPE-1、IcPE、BcPE-2和rPE等4种PE单元的资源使用和最大频率。可以看到,因内部没有配置浮点运算部件,rPE的资源占用非常少。

表 1 LT-SC并行结构中PE单元的综合结果
单元寄存器个数查找表个数查找表触发器个数DSP48E个数最大频率/MHz
BcPE-11 1561 2121 2280188.56
IcPE 3388478682169.21
BcPE-2 2 0341 9422 0540150.76
rPE 3440420550.01

3.2 性能比较

用于比较的台式计算机配置为主频2.33 GHz 英特尔酷睿2处理器、4 GB主存及CentOS 5.5操作系统。C语言编译器为gcc 4.3,优化选项设为O2。图 6给出LT-SC阵列并行结构与通用处理器上软件实现相比所获得的加速比和最大频率随矩阵规模(阵列规模)增加的变化。可以看到,在各个矩阵规模上的计算中LT-SC阵列并行结构均可以取得10倍以上的加速比,即使最大频率在降低,加速比也随着矩阵规模的增加逐步提升。需要指出的是,将矩阵规模从10×10维增加到18×18维过程中,随着IcPE数量的大量增加,FPGA上的DSP48E Slice资源全部消耗,故最大频率迅速下降。而矩阵规模接近32×32维时,随着FPGA上LUT Flip-Flop pairs资源的耗尽,最大频率大幅度下降。

图 6 LT-SC阵列并行结构与通用处理器软件实现的加速比和最大频率随矩阵规模增加的变化

4 结 论

本文基于子矩阵更新同一化算法,针对Cholesky分解的求解特征,提出以求解矩阵L的转置阵LT的LT-SC算法,并在分析LT-SC线性代数特征的基础上,给出了LT-SC的并行结构实现方案。在并行结构的具体实现中,采用FPGA协处理器浮点数学库libhdlfltp支持阵列PE的浮点流水运算。理论分析表明,该算法针对矩阵三角化计算过程具有更高的数据并行性与流水并行性; 实验结果表明,针对Cholesky分解,基于LT-SC算法在FPGA上的实现在核心计算性能上优于在通用处理器上的软件实现,可达到10倍以上的加速比。但作为实际运行硬件的可重构计算验证平台中的通信开销所占据整体运行时间比例过大,已经成为整体计算加速的瓶颈,下一步应加以改进。

参考文献
[1] Satchidanand G H, Sotirios G Z. FPGA implementation of a Cholesky algorithm for a shared-memory multiprocessor architecture[J]. Parallel Algorithms and Applications , 2004, 19 (4) : 211–226. DOI:10.1080/10637190412331279957
[2] Kung H T, Leiserson C E. Systolic arrays technical report [R]. Pitsburg, PA, USA: Carneige Mellon University, 1978.
[3] YANG Depeng, Gregory D P, LI Husheng. Compressed sensing and Cholesky decomposition on FPGAs and GPUs[J]. Parallel Computing , 2012, 38 (8) : 421–437. DOI:10.1016/j.parco.2012.03.001
[4] Alfredo B, Julien L, Jakub K, et al. A class of parallel tiled linear algebra algorithms for multicore architectures[J]. Parallel Computing , 2009, 35 (1) : 38–53. DOI:10.1016/j.parco.2008.10.002
[5] DOU Yong, Vassiliadis S, Kuzmanov G, et al. 64-bit floating-point FPGA matrix multiplication [C]//Proceedings of the 13th ACM/SIGDA International Symposium on Field Programmable Gate Arrays(FPGA’05), New York, NY, USA: ACM, 2005: 86-95. http://cn.bing.com/academic/profile?id=2061624656&encoded=0&v=paper_preview&mkt=zh-cn
[6] Vasily V, James W D. LU, QR and Cholesky factorizations using vector capabilities of GPUs [R]. Berkeley, CA, USA: EECS Department, University of California, Berkeley, 2008. http://cn.bing.com/academic/profile?id=2187691968&encoded=0&v=paper_preview&mkt=zh-cn
[7] Oleg M, Volodymyr L, Anatoli S, et al. Parallel implementation of Cholesky LLT-algorithm in FPGA-based processor[J]. Parallel Processing and Applied Mathematics , 2008, 49 (67) : 137–147.
[8] Eunice E S, CHU Peiyun. Efficient and optimal parallel algorithms for Cholesky decomposition[J]. Journal of Mathematical Modeling and Algorithms , 2003, 2 (3) : 217–234. DOI:10.1023/B:JMMA.0000015832.41014.ed
[9] Aatonio R, George A. A high throughput FPGA-based floating point conjugate gradient implementation for dense matrices[J]. ACM Transactions on Reconfigurable Technology and Systems , 2010, 3 (1) : 1–19.
[10] Badia J M, Vidal A M. Parallel algorithms to computer the eigenvalues and eigenvectors of symmetric Toeplitz matrices[J]. Parallel Algorithms and Applications , 1998, 13 (1) : 75–93. DOI:10.1080/01495739808947361
[11] Anderson E, BAI Z, Bischof C, et al. LAPACK Users Guide[M]. Philadelphia, PA, USA: The society of industrial and applied mathematics, 1999 .
[12] Blackford L S, Choi J, Cleary A, et al. ScaLAPACK Users Guide[M]. Philadelphia, PA, USA: The Society of Industrial and Applied Mathematics, 1997 .
[13] WU Guiming, DOU Yong, SUN Junqing, et al. A high performance and memory efficient LU decomposer on FPGAs[J]. Computers, IEEE Transactions , 2012, 61 (3) : 366–378. DOI:10.1109/TC.2010.278
[14] WU Guiming, DOU Yong, Gregory D P. Blocking LU decomposition for FPGAs [C]//Proceedings of the 18th IEEE Symposium on Field-Programmable Custom Computing Machines(FCCM’10). Charlotte, NC,USA: IEEE, 2010, 109-112.
[15] Haridas S G, Sotirios G Z. FPGA implementation of a Cholesky algorithm for a shared memory multiprocessor architecture[J]. International Journal of Parallel, Emergent and Distributed Systems , 2004, 19 (4) : 211–226.
[16] 郭磊, 唐玉华, 周杰, 等. 基于FPGA的Cholesky分解细粒度并行结构与实现[J]. 计算机研究与发展, 48 (Suppl) , 2011, 48 (supp) : 258–265. GUO Lei, TANG Yuhua, ZHOU Jie, et al. Fine-grained architecture and implementation for Cholesky decomposition on FPGA[J]. Journal of Computer Research and Development , 2011, 48 (supp) : 258–265. (in Chinese)
[17] WANG Xiaojun, Leeser M. A truly two-dimensional systolic array FPGA implementation of QR decomposition[J]. ACM Transactions on Embedded Computing Systems (TECS) , 2009, 9 (1) : 17–18.
[18] 邬贵明, 窦勇, 王淼. Cholesky分解细粒度并行算法[J]. 计算机工程与科学 , 2010, 32 (9) : 102–106. WU Guiming, DOU Yong, WANG Miao. A fine-grained parallel algorithm for the Cholesky decomposition[J]. Computer Engineering & Science , 2010, 32 (9) : 102–106. (in Chinese)
[19] 刘书勇, 吴艳霞, 张博文, 等. 基于可重构计算系统的矩阵三角化分解硬件并行结构研究[J]. 电子学报, 43(8) , 2015, 43 (8) : 1642–1650. LIU Shuyong, WU Yanxia, ZHANG Bowen, et al. Research of parallel hardware architecture for matrix triangularization decomposition based on reconfigurable computing system[J]. ACTA Electronica Sinica , 2015, 43 (8) : 1642–1650. (in Chinese)