CN102129420B

CN102129420B - 基于Cholesky分解解决最小二乘问题的FPGA实现装置

Info

Publication number: CN102129420B
Application number: CN 201110053248
Authority: CN
Inventors: 彭宇; 刘大同; 乔立岩; 王少军; 刘琦; 仲雪洁; 王建民
Original assignee: Harbin Institute of Technology
Current assignee: Harbin Institute of Technology
Priority date: 2011-03-07
Filing date: 2011-03-07
Publication date: 2013-03-20
Anticipated expiration: 2031-03-07
Also published as: CN102129420A

Abstract

基于Cholesky分解解决最小二乘问题的FPGA实现装置，涉及基于Cholesky分解解决最小二乘问题的FPGA实现装置，适用于最小二乘问题的求解，解决了PC机的计算效率不能满足实时和嵌入式应用的问题，它包括待求矩阵输入接口模块、分解模块和求解模块，待求矩阵输入接口模块的输出端连接在分解模块的输入端，分解模块的输出端连接在求解模块的输入端，用于满足实时、低功耗和嵌入式应用。

Description

基于Cholesky分解解决最小二乘问题的FPGA实现装置

技术领域

本发明涉及基于Cholesky分解解决最小二乘问题的FPGA实现装置。

背景技术

求解线性方程组的解的问题可以看作最小二乘问题的求解，目前，主要在冯·诺依曼结构的PC机上实现，PC机的计算效率不能满足实时和嵌入式应用的需求；采用ASIC(专用集成电路)的方法可以提高运算效率，但是适用性差，且成本较高。

发明内容

本发明的目的是为了解决现有PC机的计算效率不能满足实时和嵌入式应用的问题，提供一种基于Cholesky分解解决最小二乘问题的FPGA实现装置。

基于Cholesky分解解决最小二乘问题的FPGA实现装置，它包括待求矩阵输入接口模块6、分解模块1和求解模块2，待求矩阵输入接口模块6的输出端连接在分解模块1的输入端，分解模块1的输出端连接在求解模块2的输入端。

利用FPGA实现256维矩阵的改进Cholesky分解的运算时间与PC机平台实现同一矩阵的改进Cholesky分解的对比实验情况如下表所示：

计算平台	计算时间
		PC机	101.563ms
FPGA	12.263ms

其中FPGA的工作频率是100MHZ，实验用的PC机配置为：Pentium Dual core CPU，2.60GHz，2G DDR2存储器。由上表可以看出，利用FPGA开发实现解256维的最小二乘问题，可以实现比PC机平台的运算效率提高8倍以上，计算精度能精确到10-9。

附图说明

图1为本发明的系统结构示意图，图2为本发明的分解模块的结构示意图，图3为本发明的PE_D模块的组成和结构示意图，图4为本发明的PE_L运算模块的组成和结构示意图，图5为本发明的求解模块的结构示意图，图6为本发明的PE模块的组成和结构示意图。

具体实施方式

具体实施方式一：

结合图1说明本实施方式，本实施方式包括待求矩阵输入接口模块6、分解模块1和求解模块2，待求矩阵输入接口模块6的输出端连接在分解模块1的输入端，分解模块1的输出端连接在求解模块2的输入端。

对于一个n维的线性方程组，设为：Ax＝b其中A为n×n维矩阵，x为n维待求的解向量，b为n维列向量。若想求解线性方程组的解向量x，则需要求解n×n维矩阵A的逆，求逆矩阵的方法有很多，如线性代数中介绍的伴随矩阵法、初等变换法、分块矩阵法等，还有一些工程中常用的求逆方法，如Gauss-Jordan消去法求逆阵、矩阵分解求逆等，上述各种矩阵求逆的方法中，计算量大，对存储空间的需求也较大，不利于硬件实现。

矩阵分解实现求逆克服了上述方法的缺点，三角矩阵求逆硬件实现简单，可采用硬件平台的并行结构实现，且运算速度快，因此可以采用三角分解的方法实现矩阵求逆的运算，可以采用LU分解、QR分解、Cholesky分解(为该领域常用算法)等方法来求解，最终将矩阵求逆问题转化为最小二乘问题的求解。矩阵分解算法中，QR分解的计算复杂度比较高，硬件实现的代价比较高，因此，主要对比LU分解和Cholesky分解，主要从矩阵适用范围、计算复杂度、硬件实现资源占用情况等三方面来进行比较分析：

适用范围：LU分解适合所有非奇异矩阵，即行列式不为零的矩阵，适用范围比较大；Cholesky分解适合对称正定矩阵，适用条件比较苛刻。但是对于线性方程组：Ax＝b，如果矩阵A不是对称正定矩阵，则可以通过变换转化为对称正定阵，对方程组做如下变换：A^TAx＝A^Tb，设B＝A^TA，则B为对称正定阵，A^Tb为n维列向量。经过变换后的线性方程组就可以用Cholesky分解来实现了。

计算复杂度：Cholesky分解和LU分解的计算格式较简单，但是Cholesky算法的运算量是LU分解法的一半，逻辑操作和数据移动都比LU分解少，程序化实现容易。

硬件实现资源占用情况：对一个分块后4×4的小矩阵，需要一个LU分解模块，L矩阵求逆模块，U矩阵求逆模块，4×4矩阵乘法模块，两个除法器单元，需要的PE(processelement运算单元)比较多，占用资源较大。Cholesky分解的计算比较简单，只需要分解模块和求解模块。每个模块采用十个PE单元并行的流水线结构和一个除法器。需要的DSP处理单元数据较LU分解少了一半，占用FPGA内部运算单元少。

综上所述，可以得出以下结论：由于本发明的应用背景，处理的矩阵为对称正定阵，因此基于FPGA的矩阵分解采用改进的Cholesky分解算法。

改进Cholesky分解的基本原理

设A＝(a_ij)∈R^n×n是对称正定矩阵，则可以对矩阵A进行Cholesky分解，直接的Cholesky分解需要进行开方运算难于硬件实现，因此引入Cholesky分解的改进算法。

令A＝LDL^T，其中L为单位下三角矩阵，D为对角阵，L^T为L的转置矩阵。

容易得到D和L中的元素：

\{\begin{matrix} d_{r} = (a_{rr} - Σ_{k = 1}^{r - 1} l_{rk}^{2} d_{k}), \\ l_{ir} = (a_{ir} - Σ_{k = 1}^{r - 1} l_{ik} d_{k} l_{rk}) / d_{r} . \end{matrix}

其中，r＝1，2，…，n；i＝r+1，r+2，…，n。

具体实施方式二：

结合图2、图3和图4说明本实施方式，本实施方式的分解模块1包括PE_D模块3、多个PE_L运算模块4、多个分解结果产生模块L_ij 15、控制模块7和开关模块8，待求矩阵输入接口模块6的输出端连接在PE_D模块3的一个输入端，PE_D模块3的输出端分别连接在多个PE_L运算模块4的输入端，多个PE_L运算模块4的数据传送端分别与对应的分解结果产生模块L_ij 15的数据传送端连通，控制模块7的输出端连接在开关模块8的控制信号输入端，开关模块8的数据输入端一次与一个分解结果产生模块L_ij 15的输出端连通，实现每开关一次把一个分解结果产生模块L_ij 15的数据传送端连通，开关模块8的输出端连接在PE_D模块3的另一个输入端。其它组成和连接关系与实施方式一相同。

由于分解模块需要计算对角阵D和下三角阵L，因此需要设计两个PE单元：计算对角阵D的PE_D模块3以及计算下三角阵L的PE_L运算模块4，两个模块交替、并行运行可计算出结果。

由于下三角阵L需要计算的值比较多，为提升计算效率可以采用多个PE_L运算模块4并行运算的方式，PE_L运算模块4并行的个数越多计算效率的提升就越明显，综合片内资源的占用情况及计算效率等各方面因素，选用8个PE_L运算模块4并行的方式计算下三角矩阵L中的元素。

矩阵L中的每一个元素都需要与对角阵D中的元素做除法运算，但是由于除法运算的延迟时间较长，每一个元素的除法运算时间会使整个模块的计算效率急速下降，因此考虑在计算对角阵元素后先做除法运算得到对角阵元素d_r的倒数，计算矩阵L元素时的除法运算就可以转变为乘法运算，从而提高计算速度。

PE_D模块3的组成和结构：

PE_D模块3主要用来计算对角矩阵D中的元素，计算公式为：

其中r＝1，2，…，n；a_rr为待分解矩阵的对角线上的元素。

计算d_r涉及到乘法运算、加法运算及减法运算，为计算矩阵L还需要计算1/d_r，因此需要用到乘法器、加法器、减法器及除法器，两个乘法器用来计算

由于乘法器运算时有延迟时间，因此采用FIFO1作为第二个乘法器输入的缓存，FIFO1的深度只要大于乘法器的延迟时钟的个数即可。

加法器用来实现累加算式

的计算，由于加法器延迟时间的存在，使得累加计算的速度取决于加法器的延迟时间，每隔加法器延迟时钟的时间才能进行下一步累加的计算。得到的累加结果与矩阵A的对角线元素做减法可得到d_r的值，通过除法器的除法运算可以得到1/d_r，将1/d_r储存在分解结果产生模块1/dr14中。

PE_L运算模块4的组成和结构：

PE_L运算模块主要用来计算下三角矩阵L中的元素，计算公式为：

l_{ir} = (a_{ir} - Σ_{k = 1}^{r - 1} l_{ik} d_{k} l_{rk}) / d_{r .},

其中，r＝1，2，…，n；i＝r+1，r+2，…，n

涉及到乘法运算、累加运算、减法运算及除法运算，1/d_r已由PE_D模块算出，因此将除法运算转化为乘法运算，算式中需要计算l_ikd_kl_rk，d_kl_rk已由PE_D模块3在计算d_r时由第一个乘法器算出，因此在PE_L运算模块4中放置一个FIFO存放由PE_D模块3算出的d_kl_rk，这样做的好处是可以省掉一个乘法器且节省计算时间。

FIFO_pe的深度取决于待分解矩阵的维数以及使用的PE_L运算模块4的个数，例如，如果矩阵的维数是256维，使用8个PE_L运算模块4，则FIFO_pe的深度为32(256/8)。每计算一轮l_ir，即每完成矩阵L的一列的计算，需要对FIFO_pe复位一次，以便进行下一轮l_ir的计算。

减法运算：与乘、加模块并行运算，实时检测加法器的rdy信号，判断何时进行减法运算。减法器输出的计算结果需要与PE_D模块3计算出的1/d_r做乘法运算，即可得到l_ir，将结果存入到分解结果产生模块L_ij 15中。

对角阵D的计算需要l_rk，因此分解结果产生模块L_ij 15中的数据需要切换输出给PE_D和各自对应的PE_L运算模块4输入端，因此需要设计一个32位宽度，8端口的总线开关模块8，实现每个分解结果产生模块L_ij 15与两个PE单元的切换。

模块调用IP核的设定

PE单元的设计需要调用IP核，IP核有多种生成方式，本方案采用使用DSP48E来搭建，每个IP核计算时的延迟时间可以设定，不同的延迟时间所对应的时钟的上限不同，综合占用资源及计算效率等方面的考虑，设定乘法器的延迟时间为6个时钟，加法器的延迟时间和减法器的延迟时间设为5个时钟，除法器的延迟时间为10个时钟。

IP核	DSP48E的个数	延迟时间
			乘法器	3	6
加法器	2	5
			减法器	2	5
除法器	0(logic only)	10

具体实施方式三：

结合图5和图6说明本实施方式，本实施方式的求解模块2包括控制单元16、多个PE单元5、减法器9、RAM_z10、选通开关11、RAM_b12、乘法器13和分解结果产生模块1/dr14，多个PE单元5的一个输入端分别对应连接在分解结果产生模块L_ij 15的结果输出端，控制单元16的输出端连接在选通开关11的控制输入端，选通开关11每次连通一个PE单元5的数据传送端，选通开关11的输出端连接在减法器9的一个数据输入端，减法器9的数据输出端分别连接在RAM_z10和乘法器13的一个数据输入端，分解结果产生模块1/dr14的输出端连接在乘法器13的另一个数据输入端，分解结果产生模块1/dr14的输入端连接在PE_D模块3的1/dr输出端，乘法器13的数据输出端连接在RAM_b12的输入端，RAM_b12的输出端连接在减法器9的另一个数据输入端，RAM_z10的输出端分别连接在每个PE单元5的另一个输入端。其它组成和连接关系与实施方式一相同。

求解模块2的组成和结构：

设计原理和说明

在进行完矩阵的Cholesky分解后，可将原矩阵方程化为如下模式：

故求解部分根据以上等式分为三个部分Lz＝b，Dr＝z，L^Tη＝r。

第一部分：Lz＝b

原方程可展开为：

z₁＝b₁

z₂＝(b₂-l₂₁z₁)

z₃＝(b₃-l₃₁z₁-l₃₂z₂)(1)

z_{n} = (b_{n} - Σ_{i = 1}^{n - 1} l_{ni} z_{i})

具体求解过程下述。

第二部分：Dr＝z

由于d矩阵为对角阵，其求解过程等价于求解r_n＝z_n/d_n。利用矩阵分解部分求取并存储的1/d_n，可直接对结果进行乘法运算，提高计算效率。

第三部分：L^Tη＝r

在第二部分得到向量r后可以通过下列三角线性方程组求得变量η

η_n＝r_n

η_n-1＝r_n-1-u_n-1nη_n (2)

η_{i} = r_{i} - Σ_{p = i + 1}^{n} u_{ip} η_{p}

(1)和(2)从本质上来说是一样的三角计算式，以(1)式为例，硬件实现的步骤按下表进行。

表1

采用流水线结构设计，共需八个并行的PE单元。

图6的PE单元5中的l_ij为分解结果产生模块L_ij 15中存储的8块L矩阵，八路并行通过加法器输出累加结果，控制单元16控制选通开关11切换，保证有效时序内正确结果的采集。

分解结果产生模块1/dr14中的

为分解模块1中的FIFO存储结果，将其与减法器输出结果(即第一部分中的结果z，此结果亦同步写入RAM_z10中作为下一循环输入端口的输入)相乘，所得即为第二部分结果r，再存入RAM_b12的相应地址中(RAM_b12为双端口RAM，控制时序不产生读写冲突)，即可实现第二部分的运算。这样设计的优势在于可以不必预留出第二部分的计算时间，使第一部分计算和第二部分计算同时进行，提升效率。

第三部分和第一部分只在控制单元16有所改变，将开关的正序切换改为倒序切换，最终存储的RAM_z10为求解模块2的最终结果，并输出结束标志位。

程序的时序由控制单元16控制。控制单元16产生起始标志位，计数器定时反馈控制单元16数据信息，通过控制单元16实现开关的切换、地址的变化和产生最终结束标志位。以计数器控制数据的选取，即通过记录RAM_1的数据数目来标志运算轮数，并且给乘法器相应标志信号，经加法器传递，到达减法器来实现正确结果的输出。

存储设备按照256维矩阵设计。FIFO深度为8192，RAM_b12为外置RAM，最终结果存在RAM_z10中，两个RAM深度均为257。

资源占用指标如下所示：

IP核	DSP48E的个数	延迟时间
			乘法器	3	3
加法器	2	3
			减法器	2	3

Claims

1.基于Cholesky分解解决最小二乘问题的FPGA实现装置，其特征是它包括待求矩阵输入接口模块(6)、分解模块(1)和求解模块(2)，待求矩阵输入接口模块(6)的输出端连接在分解模块(1)的输入端，分解模块(1)的输出端连接在求解模块(2)的输入端；

分解模块(1)包括PE_D模块(3)、多个PE_L运算模块(4)、多个分解结果产生模块L_ij(15)、控制模块(7)和开关模块(8)，待求矩阵输入接口模块(6)的输出端连接在PE_D模块(3)的一个输入端，PE_D模块(3)的输出端分别连接在多个PE_L运算模块(4)的输入端，多个PE_L运算模块(4)的数据传送端分别与对应的分解结果产生模块L_ij(15)的数据传送端连通，控制模块(7)的输出端连接在开关模块(8)的控制信号输入端，开关模块(8)的数据输入端一次与一个分解结果产生模块L_ij(15)的输出端连通，实现每开关一次把一个分解结果产生模块L_ij(15)的数据传送端连通，开关模块(8)的输出端连接在PE_D模块(3)的另一个输入端；

求解模块(2)包括控制单元(16)、多个PE单元(5)、减法器(9)、RAM_z(10)、选通开关(11)、RAM_b(12)、乘法器(13)和分解结果产生模块1/dr(14)，多个PE单元(5)的一个输入端分别对应连接在分解结果产生模块L_ij(15)的结果输出端，控制单元(16)的输出端连接在选通开关(11)的控制输入端，选通开关(11)每次连通一个PE单元(5)的数据传送端，选通开关(11)的输出端连接在减法器(9)的一个数据输入端，减法器(9)的数据输出端分别连接在RAM_z(10)和乘法器(13)的一个数据输入端，分解结果产生模块1/dr(14)的输出端连接在乘法器(13)的另一个数据输入端，分解结果产生模块1/dr(14)的输入端连接在PE_D模块(3)的1/dr输出端，乘法器(13)的数据输出端连接在RAM_b(12)的输入端，RAM_b(12)的输出端连接在减法器(9)的另一个数据输入端，RAM_z(10)的输出端分别连接在每个PE单元(5)的另一个输入端。