正定Hermite矩阵Cholesky分解高速脉动阵列的设计方法

摘要

本发明公开了正定Hermite矩阵的Cholesky分解高速脉动阵列的设计方法。本发明将处理单元划分为圆形处理单元和方形处理单元，并对不同类别的处理单元进行了功能划分，保证硬件设计具有规律性，实现数据的并行化流水处理。所述设计方法适合所有维数为M×M的正定Hermite矩阵的Cholesky分解，具有通用性，适用范围广。

说明书

技术领域

本发明涉及正定Hermite矩阵的分解方法，具体是指正定Hermite矩阵 Cholesky分解脉动阵列的方法。

背景技术

在科学计算和工程领域中，Cholesky分解是一个非常重要的线性代数法。 尤其是在第三代移动通信TD-SCDMA系统中，多种联合检测算法都要涉及到系 统矩阵A求逆问题，而在对矩阵A求逆时，首先要对矩阵A进行分解，当用户 比较多且接收天线比较多时，将矩阵A分解的计算量将会比较庞大，复杂度也 比较高。Cholesky分解因为利用了Hermite矩阵的性质，在对矩阵A进行分解时， 有效地降低了运算复杂度。

而基于Cholesky分解正定Hermite矩阵的方法大多是根据Cholesky分解公 式直接编写的，这种方法因为把分解当作一个整体，没有体现分块设计思想， 特别是矩阵维数较大时，数据之间的关系复杂而混乱，难以保证设计的正确性， 并且因为数据计算较复杂，因此硬件的并行处理速度并不高。

发明内容

本发明的目的在于克服通过Cholesky分解公式直接编写的硬件设计不适用 于维数较大的矩阵、并且硬件并行处理速度不高的缺陷，提供一种涉及M×M维 的正定Hermite矩阵的Cholesky分解高速脉动阵列的设计方法。

本发明的目的通过下述技术方案实现：

正定Hermite矩阵Cholesky分解高速脉动阵列的设计方法，包括以下步骤：

(a)根据维数为M×M正定Hermite矩阵确定脉动阵列的维数以及构架；

(b)根据脉动阵列的架构对脉动阵列的处理单元进行分类；

(c)根据处理单元类别设计处理单元结构；

(d)根据脉动阵列的架构连接处理单元；

(e)输出脉动阵列。同时，所述步骤(a)中脉动阵列为M×M维的包含 斜边的上三角处理单元阵列，且所述脉动阵列共包含(1+M)M/2个处理单元。

并且，所述步骤(a)中脉动阵列为M×M维的上三角单元处理阵列，且所 述脉动阵列共有(1+M)M/2个处理单元。

进一步的，所述步骤(c)中对圆形处理单元结构的设计具体包括以下步 骤：

(c1)将脉动阵列上位于第M列、第M行的处理单元划分为第一类圆形处 理单元，所述第一类圆形处理单元内置有一个开方器和一个寄存器，且包含一 个输入端和输出端；

(c2)将脉动阵列斜边上的其他处理单元划分为第二类圆形处理单元，所 述第二类圆形处理单元内置有一个开方器、一个除法器和两个寄存器，且包含 一个输入端和两个输出端。

更进一步的，所述步骤(c)中对方形处理单元结构的设计包括以下步骤：

(I)将平行于圆形处理单元且与圆形处理单元相邻的处理单元划分为第 一类方形处理单元，所述第一类方形处理单元由两个乘法器、一个加法器、一 个选择器、两个移位寄存器和两个寄存器组成，且包含两个输入端和输出端；

(II)将平行于第一类方形处理单元且与第一类方形处理单元相邻的处理 单元划分为第二类方形处理单元，所述第二类方形处理单元由两个乘法器、一 个加法器、三个选择器两个移位寄存器和两个寄存器组成，且包含三个输入端 和两个输出端；

(III)以此类推，将平行于第K-1类方形处理单元且与第K-1类方形处理 单元相邻的处理单元划分为第K类方形处理单元，所述第K类方形处理单元由 两个乘法器、一个加法器、三个选择器、两个移位寄存器和两个寄存器组成， 其中第K类方形单元有K+1个输入端，两个输出端。

并且所述K≤M。

为了更好的实现本发明，所述步骤(d)中在连接处理单元时，脉动阵列上 同时连接有延迟单元。

本发明较现有技术相比，具有以下优点及有益效果：

(1)本发明采用脉动阵列结构，数据按预先确定的方式在脉动阵列的处理 单元间有规律地传输，所有处理单元能同时并行地对流经它的数据进行处理， 解决了常规设计中硬件并行处理速度不高的缺陷，提高了处理单元数据处理效 率；

(2)本发明对整个矩阵进行了分解处理，将脉动阵列结构划分为不同类型 的处理单元，实现处理单元的分工合作，保证硬件设计具有规律性，不仅适合 维数较小的矩阵，同时适合维数较大的矩阵计算，从而在提高计算正确率的同 时，扩大了适用范围，使其具有了更为广阔的应用领域；

(3)本发明在方形处理单元中尽量复用乘法器，避免了方形处理单元使用 乘法器数量随分解矩阵维数的增加而快速增加的缺陷，因而有效的减少了硬件 资源的开销；

(4)本发明中，将数据输入处理单元阵列后，处理单元阵列能自动完成数 据处理，脉动阵列的运行速度大大提高，同时简化了逻辑控制；

(5)本发明的处理单元结构和阵列简单，并且具有一致的规则，可以达到 很高的模块化程度，非常适合超大规模集成电路的设计和制造；

(6)本发明处理单元的规则性强，适合所有维数为M×M的正定Hermite 矩阵的Cholesky分解，具有通用性，适用范围广。

附图说明

图1为本发明-实施例4×4维的脉动阵列结构图。

图2为本发明-实施例的输入输出时序图。

图3为本发明的第一类圆形处理单元的结构图。

图4为本发明的第二类圆形处理单元的结构图。

图5为本发明的第一类方形处理单元的结构图。

图6为本发明的第二类方形处理单元的结构图。

图7为本发明的第三类方形处理单元的结构图。

图8为本发明-实施例的Verilog HDL定点仿真与Matlab浮点仿真的误差 曲线图。

图9为本发明-实施例的流程图。

具体实施方式

下面结合实施例对本发明作进一步的详细说明，但本发明的实施方式不限 于此。

实施例

本实施例中，使用Matlab仿真Cholesky分解脉动阵列实现正定Hermite矩 阵的分解。

正定Hermite矩阵A为4×4维矩阵，由Hermite矩阵的性质：A＝A^H，即 矩阵的元素将A表示为：

$A=\left(\begin{array}{cccc}{a}_{11}& a_{21}^{*}& a_{31}^{*}& a_{41}^{*}\\ a_{21}& a_{22}& a_{32}^{*}& a_{42}^{*}\\ a_{31}& a_{32}& a_{33}& a_{43}^{*}\\ a_{41}& a_{42}& a_{43}& a_{44}\end{array}\right)$

显然，矩阵A满足Cholesky分解的要求，因此待分解矩阵处于斜边上半部分 的元素完全可以用斜边左下方的元素表示，输入矩阵A包含斜边元素在内的下 三角的元素就足以表示整个矩阵A。

下三角矩阵L中每个元素l_ij(i＝1，2，3，4；j＝1，2，3，4)的计算公式如下：

$l_{\mathrm{ij}}=\left(\begin{array}{ccc}\sqrt{a_{\mathrm{ii}}-\underset{k=1}{\overset{i-1}{\Sigma }}{l}_{\mathrm{ik}}{l}_{\mathrm{ik}}^{*}},& \mathrm{for}& i=j\\ \frac{1}{l_{\mathrm{ij}}}(a_{\mathrm{ij}}-\underset{k=1}{\overset{j-1}{\Sigma }}{l}_{\mathrm{ik}}{l}_{\mathrm{jk}}^{*}),& \mathrm{for}& i>j\\ 0,& \mathrm{for}& i<j\end{array}\right)$

其中，下三角矩阵L的主对角元素l_ij均为正实数。

根据上述公式的展开可得矩阵L第一列的所有元素：

$\left(\begin{array}{c}{l}_{11}=\sqrt{a_{11}}\\ l_{21}=\frac{a_{21}}{l_{11}}\\ l_{31}=\frac{a_{31}}{l_{11}}\\ l_{41}=\frac{a_{41}}{l_{11}}\end{array}\right)$

分析上式可知，作为一个模块，包含开方和除法两个部分，而该模 块的输出l₁₁作为l₂₁、l₃₁、l₄₁的输入，因此便可以得到矩阵L第一列的所有元素。

矩阵L第二列的计算公式如下：

$\left(\begin{array}{c}\\ l_{22}=\sqrt{a_{22}-l_{21}{l}_{21}^{*}}\\ l_{32}=\frac{1}{l_{22}}(a_{32}-l_{31}{l}_{21}^{*})\\ l_{42}=\frac{1}{l_{22}}(a_{42}-l_{41}{l}_{21}^{*})\end{array}\right)$

矩阵L第三列的计算公式如下：

$\left(\begin{array}{c}{l}_{33}=\sqrt{a_{33}-l_{31}{l}_{31}^{*}-l_{32}{l}_{32}^{*}}\\ l_{43}=\frac{1}{l_{43}}(a_{43}-l_{41}{l}_{31}^{*}-l_{42}{l}_{32}^{*})\end{array}\right)$

根据第二列和第三列的计算公式可得，也同样包含开方和除法两部 分。矩阵L第一列的计算结果l₂₁、l₃₁、l₄₁都作为第二列、第三列计算的输入，具体 说来，l₂₁作为求l₂₂的输入，l₃₁作为求l₃₂的输入，l₄₁作为求l₄₂的输入。

同时，第四行第四列的元素计算公式如下：

$l_{44}=\sqrt{a_{44}-l_{41}{l}_{41}^{*}-l_{42}{l}_{42}^{*}-l_{43}{l}_{43}^{*}}$

基于以上数据的流动规则，可得上述元素计算公式的脉动阵列结构图，如 图1所示，本发明-实施例的脉动阵列是一个4×4维的包含斜边在内的上三角 处理单元阵列，所述脉动阵列共包含10个处理单元，处理单元简称PE。

根据图1所示，脉动阵列上的PE分为圆形PE和方形PE，根据上述公式规 律可得，脉动阵列斜边上的PE为圆形PE，其余为方形PE。

本发明-实施例旨在实现4×4维正定Hermite矩阵通过Cholesky分解高速脉 动阵列的方法实现矩阵的分解，输入矩阵A经过一定处理后按每一行并行输出 其分解矩阵L。

对应的输入输出时序如图2所示，在输出时序图2后，如图9所示，由下 列步骤实现本发明-实施例：

步骤一：确定脉动阵列的维数以及构架

其中矩阵H是均值为0、方差为1的4阶随机复矩阵，由H生成4×4维正定 Hermite矩阵A，通过Cholesky分解设计出其脉动阵列，并用FPGA实现，写成 MATLAB代码为：

n＝4；

H＝(randn(n，n)+1j*randn(n，n))/sqrt(2)；

A＝H’*H；

用(y1，y2，y3)表示一个定点数，其中y1为符号位，取值为0时表示无 符号，取值为1时表示有符号，y2表示整数部分的位宽，y3表示小数部分的位 宽。在本发明-实施例中，将输入矩阵A定点为(1，5，19)，即表示矩阵为有 符号位，5位整数位和19位小数位。

根据矩阵L每个元素的计算公式确定脉动阵列的PE构架，并根据矩阵元素 的表示方式对各个PE进行编号，如图1所示。

步骤二：对圆形PE结构及功能进行设计。

斜边上PE为圆形PE。其中PE₄₄只包含开方运算，因此被划分为第一类圆 形PE；斜边上其余处理单元包含开方以及除法运算，因此被划分为第二类圆形 处理单元，包括PE₁₁、PE₂₂、PE₃₃。

第一类圆形处理单元结构如图3所示，所述第一类圆形PE有一个输入、一 个输出。该类圆形PE输入数据rx_yin，然后进行开方，开方的结果作为该圆形 PE的第一个输出tx_xout，计算公式如下：

$\mathrm{tx}\_\mathrm{xout}=\sqrt{\mathrm{rx}\_\mathrm{yin}}$

第二类圆形PE结构如图4所示，所述第二类圆形PE有一个输入和两个输 出，输入数据rx_yin进行开方后输出结果tx_xout，同时tx_xout进行求倒数， 计算公式如下：

$\mathrm{tx}\_\mathrm{xout}\_d=\frac{1}{\mathrm{tx}\_\mathrm{xout}}$

步骤三：对方形PE的结构和功能进行设计。

方形PE被划分为：第一类方形PE，包括PE₁₂、PE₂₃、PE₃₄；第二类方形 PE，包括PE₁₃、PE₂₄；第三类方形PE，包括PE₁₄。

第一类方形PE的结构框图如图5所示，该类方形PE的输入rx_xin连接的 是与它处于同一行的第一类圆形PE的输出tx_xout_d；所述第一类圆形PE一定 位于第一类方形PE的左侧，同时所述第一类方形PE的输入rx_yin与它处于同 一列的第二类方形PE的输出tx_xout相连。输入端rx_xin每M个时钟输入一个 数据；输入端rx_yin每M个时钟输入两个数据，第一个数据为rx_yin(1)，第二 个数据为rx_yin(2)。rx_xin与rx_yin(1)作乘法，然后将乘积作为tx_xout输出， 计算公式如下：

tx_xout＝rx_yin(1)×rx_xin

输出的同时，将tx_xout反馈回另一个乘法器的输入端，同时另一个输入端 输入tx_xout的共轭tx_xout^*，然后将该乘法器输出的结果取反后与rx_yin(2)作 加法运算，计算公式如下：

tx_yout＝-tx_xout^*×tx_xout+rx_yin(2)

因此第一类方形PE由两个乘法器、一个加法器和一个选择器组成，计算过 程中复用了一次乘法器，减少硬件资源的开销。

第二类方形PE的结构框图如图6所示，该类方形PE的输入rx_xin1连接 的是与它处于同一行的第二类圆形PE的输出tx_xout_d，输入rx_xin2连接的是 与它处于同一行的第一个方形PE的输出tx_xout，另一个输入rx_yin与它处于 同一列上方的第三类方形PE的输出tx_yout相连接，其中输入端rx_xin1与输入 端rx_xin2每M个时钟输入一个数据，输入端rx_yin每M个时钟输入三个数据， 分别为rx_yin(1)，rx_yin(2)，rx_yin(3)。

rx_xin1的第一个数据与rx_yin的第一个数据rx_yin(1)作乘法，然后将乘积 作为tx_xout输出，计算公式如下：

tx_xout＝rx_yin(1)×rx_xin1

输出的同时，将tx_xout反馈回另一个选择器的输入端，一个选择器的输入 端输入全部tx_xout，另一个选择器的一个输入端输入tx_xout^*，选择器将rx_xin2 的共轭rx_xin2^*以及tx_xout输入到另一个乘法器输入端。然后将该乘法器输出 的结果取反后与rx_yin(2)作加法运算，计算公式如下所示：

tx_yout(1)＝-rx_xin2^*×tx_xout+rx_yin(2)

选择器再将tx_xout^*以及tx_xout输入到乘法器输入端，乘法器将输出的结 果取反后与rx_yin(3)作加法运算，计算公式如下：

tx_yout(2)＝-tx_xout^*×tx_xout+rx_yin(3)

因此第二类方形PE包括两个乘法器、一个加法器和三个选择器。

第三类方形PE的结构框图如图7所示，该类方形PE共有4个输入，包括 rx_xin1、rx_xin2，rx_xin3和rx_yin，第二类圆形PE的输出tx_xout_d与rx_xin1 相连，第一类方形PE的输出tx_xout与rx_xin2相连，同时第二类方形处理单元 的输出tx_xout与rx_xin3相连，并且rx_yin依次输入a₄₁，a₄₂，a₄₃，a₄₄。

其中输入端rx_xin1、rx_xin2与rx_xin3每M个时钟输入一个数据，rx_yin 每M个时钟输入4个数据，分别为rx_yin(1)，rx_yin(2)，rx_yin(3)和rx_yin (4)。rx_xin1的第一个数据与rx_yin的第一个数据rx_yin(1)作乘法，然后将乘 积作为tx_xout输出，计算公式如下：

tx_xout＝rx_yin(1)×rx_xin1

输出的同时，将tx_xout反馈回另一个到选择器的输入端，一个选择器的输 入端输入全部tx_xout，另一个选择器的一个输入端输入tx_xout^*，选择器将 rx_xin2的共轭rx_xin2^*以及tx_xout输入到另一个乘法器输入端。然后将该乘法 器输出的结果取反后与输入rx_yin的第二个数据rx_yin(2)作加法运算，计算公 式如下：

tx_yout(1)＝-rx_xin2^*×tx_xout+rx_yin(2)

依次类推，根据上述公式可以得到tx_yout(2)：

tx_yout(2)＝-rx_xin3^*×tx_xout+rx_yin(3)

在计算tx_yout(3)时，选择器再将tx_xout^*以及tx_xout输入到乘法器输入端， 乘法器输出的结果取反后与输入rx_yin的第三个数据rx_yin(3)作加法运算，计 算公式如下：

tx_yout(3)＝-tx_xout^*×tx_xout+rx_yin(4)

所述第三类方形PE包括两个乘法器、一个加法器和三个选择器。

这样就完成了各个PE的结构设计。

步骤四：按照脉动阵列架构连接各个PE。

按照图1所示，将各个PE互联，形成满足要求的脉动阵列。在连接各个 PE的时候，根据各个PE实际的输入和输出时序，添加必要的延迟单元，使得 各个PE满足输入与输出的时序要求，实现数据的并行流水化处理。

步骤五：输出脉动阵列。

添加延迟单元后保证每一行PE的输出tx_xout作为脉动阵列在该行的有次 序的输出，输出的顺序与PE所处的位置一致，从左到右按顺序依次进行，便于 后续模块的接入或者数据观察。

按照上述步骤完成数据处理，便设计出完整的脉动阵列。

最后输出矩阵仍为(1，5，19)的定点形式，说明脉动阵列中的数据计算 并没有出现错误。

在Verilog HDL的仿真中，选取仿真工具为Xilinx ISE 13.3及Modelsim SE 6.6f，测试得到矩阵从完全输入到它的逆矩阵完全输出共需经历355个时钟。 输入500组符合本发明实例的4阶正定Hermit随机复矩阵，将Verilog HDL仿 真得到的Cholesky分解的结果与Matlab浮点仿真结果进行比较，计算每个矩 阵中16个元素的平均相对误差，绘制出本发明方案下Verilog HDL定点仿真与 Matlab浮点仿真的误差曲线，如图8所示。对这500个统计结果进一步取均值， 得到平均误差为5.0545e-004。可以看到，采用本发明的设计方法，并且选择 合适的定点方案，计算得到的结果具有比较高的精确度。

在Xilinx ISE中选取硬件平台XC7VX485T对本发明实例的工程进行综合， 得到布局布线后的综合报告并将其占用的主要硬件资源列于表1中，其中，设 计中的各乘法器以及除法器使用ISE中相应的IP核，且采用Mults方式实现。 综合得到的最大时钟频率在280.513MHz。

表1本发明实例综合报告中主要硬件资源占有情况

 资源类型   使用数量   可用数量   占有率  Slice Registers   29,790   607,200   4％  Slice LUTs   14,572   303,600   4％  bonded IOBs   385   600   64％  Block RAM/FIFO   39   1,030   3％  DSP48E1s   72   2,800   2％

如上所述，便可以很好的实现本发明。