CN105426345A - 一种矩阵求逆运算方法 - Google Patents

Abstract

本发明涉及一种矩阵求逆运算方法，包括如下步骤：1）首先进行列选主元LU分解：根据公式PA=LU，将源矩阵A分解为单位下三角矩阵L、上三角矩阵U和置换矩阵P；2）接着进行三角矩阵求逆：对L矩阵求逆得到其逆矩阵L^-1，将U矩阵的转置矩阵求逆后再转置得到U^-1；3）最后进行矩阵相乘：将矩阵U^-1和矩阵？L^-1相乘，并根据置换矩阵P将矩阵乘法结果进行列变换得到源矩阵A^-1。有益效果为：通过使用列选主元LU分解算法，有效地降低矩阵求逆算法的时间复杂度，增加矩阵求逆运算的可并行性，减少矩阵求逆运算的时间，并且支持任意阶数的矩阵求逆运算，可以根据运算点数需求增加或减少硬件资源，更好地满足实际应用的需求。

Description

一种矩阵求逆运算方法

技术领域

本发明涉及矩阵求逆方法，尤其涉及基于时分复用技术的矩阵求逆运算的VLSI设计方法。

背景技术

矩阵求逆的方法有很多，包括伴随矩阵法、初等变换法、分块矩阵法、Gauss.Jordan消去法、矩阵分解法等等。上述各种矩阵求逆的方法中，伴随矩阵法中需要求大量的行列式，每个行列式都几乎要计算到所有的矩阵元素，计算量大，对存储空间的需求也较大，不利于实现。

矩阵分解是将矩阵分解为一些较简单矩阵的乘积，如三角矩阵或酉矩阵，他们都具有某种特性，较容易得到逆矩阵。对分解得到的矩阵求逆后相乘，得到原矩阵的逆矩阵。矩阵分解常见的有三种：1)LU分解法，2)QR分解法，3)奇异值分解法。

LU分解是一种科学计算中的矩阵分解方法，常用于线性方程的求解、矩阵求逆和行列式计算，在图形图像处理、科学计算等诸多领域也有广泛应用。LU分解法是将原正方矩阵分解成一个上三角形矩阵和一个下三角形矩阵。

例如若将矩阵A作LU分解成两个三角矩阵，即A＝LU，则A^-1＝U^-1L^-1。其中三角矩阵L和U几乎一半数据都为0，所以求逆简单。将分解得到的特殊矩阵分别求逆后再相乘，即最终得到原矩阵的逆。比较三种分解方法，发现LU分解法的可并行度最高。

当矩阵A为n阶非奇异矩并且所有顺序主子式不为0，则矩阵A可以分解为一个主对角元素全为1的下三角矩阵L和一个上三角矩阵U的乘积，即A＝LU，且分解是唯一的。

L和U的元素可由下面的递推式求出：

$\left\{\begin{array}{c}{u}_{1j}=a_{1j},(j=1,2,...n);\\ l_{i1}=\frac{a_{i1}}{u_{11}},(i=1,2,...,n);\\ u_{rj}=a_{rj}-\underset{k=1}{\overset{r-1}{\Σ}}{l}_{rk}{u}_{kj},(r=1,2,...,n;j=r,...,n)\\ l_{ir}=\frac{a_{ir}-\underset{k=1}{\overset{r-1}{\Σ}}{l}_{ik}{u}_{kr}}{u_{rr}},(r=1,2,...,n-1;i=r+1,...,n)\end{array}\right.$

由于主元出现在分母中，因此只有当u_rr≠0(k＝1,,n)时，分解才能进行到底。此外，在实际的计算过程中，即使不为零，但当它很小的时候，由于计算机的精度限制，会导致下溢，也会使分解运算不稳定。因此需要借助选主元的方法，即适当交换矩阵的行(或列)，以得到较大的非零主元。

如果分解不进行主元选取，会影响到计算精度和正确性，在碰到一些奇异矩阵时，甚至会导致计算失败。

发明内容

本发明的目的在于提供一种矩阵求逆方法，以克服现有技术领域的不足，本发明的方法包括如下步骤：

1)首先进行列选主元LU分解：根据公式PA＝LU，将源矩阵A分解为单位下三角矩阵L、上三角矩阵U和置换矩阵P；

2)接着进行三角矩阵求逆：对L矩阵求逆得到其逆矩阵L^-1，将U矩阵的转置矩阵求逆后再转置得到U^-1；

3)最后进行矩阵相乘：将矩阵U^-1和矩阵L^-1相乘，并根据置换矩阵P将矩阵乘法结果进行列变换得到源矩阵A^-1。

所述的矩阵求逆运算方法的进一步设计在于，步骤1)包括如下步骤：

I)进行选主元：从第一列中选取最大主元，与该列的第一个元素的值进行交换，并存储最大主元的行下标；

II)进行归一化：将第一行与最大主元所在行互换，得到U分解因子的一行元素u_1j，用交换后的第一列第一个元素依次除以该列的其它元素，得到L分解因子的一列元素l_i1，其中1≤i≤n，1≤j≤n，n表示矩阵阶数；

III)进行数据更新：对矩阵进行更新操作：a¹ _ij＝a_ij-l_i1*u_1j，得到新子矩阵a¹ _ij，每次更新操作后的子矩阵的阶数比原矩阵的阶数少一阶，其中，a_ij为原矩阵中的值，2≤i≤n，2≤j≤n，n表示矩阵阶数。对新的子矩阵循环重复以上三个步骤，最终得到L矩阵和U矩阵。

所述的矩阵求逆运算方法的进一步设计在于，所述步骤I)中根据当前循环次数k，将矩阵的主列元素读出，先将所述主列元素通过两个实数乘法器和一个实数加法器，流水运算得出主列元素的实部虚部平方和，再通过一个实数减法器，找出主元。

所述的矩阵求逆运算方法的进一步设计在于，所述步骤II)中根据当前循环中找到的最大主元行下标，将源矩阵的第k行与最大主元所在行互换，接着将数据交换后的主列对角线数的实部和虚部分别送入两个实数除法器中用于算出主元的倒数，再将主列对角线以下元素读出，通过一个复数乘法器与主元的倒数相乘，完成归一操作，其中k表示当前循环次数。

所述的矩阵求逆运算方法的进一步设计在于，所述步骤III)中采用四路并行运算，每路运算使用一个复数乘法器和一个复数加法器，将矩阵第k行第k+1列、第k行第k+2列、第k行第k+3列、第k行第k+4列的数据，分别送入四路运算的对应复数乘法器中，与步骤2)中归一化的结果进行复数乘法运算；复数乘法器的结果直接送入复数加法器，与矩阵第k+1列、第k+2列、第k+3列、第k+4列中第k+1行至第n行数据分别相减，当所述四列数据更新完后，再对所述四列数据后续四列数据进行数据更新操作，依此类推直至矩阵的最后一列完成数据更新操作，即完成了一次数据更新，其中k表示当前循环次数，n为矩阵阶数。

所述的矩阵求逆运算方法的进一步设计在于，所述步骤2)通过取倒单元、乘累加单元以及取反乘单元实现矩阵的求逆、转置，其中，

A.取倒单元：根据下三角矩阵求逆公式，首先将对角线元素进行取倒，对于L矩阵，结果为1；对于转置后的U矩阵，将对角线元素读出，通过两个实数乘法器、一个实数加法器、一个实数除法器和一个复数乘法器得到取倒运算的结果；

B.乘累加单元：根据下三角矩阵求逆公式分别分析L矩阵和转置后的U矩阵的数据依赖关系，通过四组由一个复数乘法器和一个复数加法器所组成的乘累加器进行乘累加运算，所述乘累加运算从第2行开始按行顺序计算非对角线元素，L矩阵和转置后的U矩阵同时进行乘累加运算；

C.取反乘单元：通过复数乘法器对步骤B)中乘累加单元的结果，与步骤A)中取倒单元的结果进行取反乘运算，即得到相应行的取反结果。

所述的矩阵求逆运算方法的进一步设计在于，所述步骤A)中将复数的实部虚部分别输出到所述两个实数乘法器，两个实数乘法器的输出端分别与所述实数加法器的两个输入端连接，实数加法器的输出端与所述实数除法器的输入端连接，实数除法器的另一输入端固定地设置为1，实数除法器的输出端与所述复数乘法器的输入端连接。

所述的矩阵求逆运算方法的进一步设计在于，所述步骤B)中L矩阵和转置后的U矩阵进行乘累加运算时，分别使用两组乘累加器，每组乘累加器包括

复数乘法器完成L、U递推公式 $s_{ij}=-\left(\underset{k=j}{\overset{i-1}{\Σ}}{l}_{ik}{s}_{kj}\right)/l_{ii},(j=1,2,...,n-1;i=j+1,...,n)$ 中l_ik与S_kj的复数相乘操作，l_ik表示L矩阵中的数，S_kj表示L矩阵的逆矩阵中的数，i表示行号，j表示列号，k表示累加循环次数；

复数加法器和延迟控制逻辑，根据复数乘法器的结果完成累加操作。

所述的矩阵求逆运算方法的进一步设计在于，所述步骤C)中对于L矩阵，将乘累加单元结果的符号位取反作为取反结果；对于U矩阵，需将乘累加单元的结果取反后再与对应行的对角线元素相乘作为取反结果。

所述的矩阵求逆运算方法的进一步设计在于，所述步骤3)中将U^-1矩阵按行划分为四等份，同时使用四路乘累加运算单元分别与L^-1矩阵进行乘累加运算，并根据步骤1)中的P矩阵，对乘法结果矩阵进行列变换得到最终的矩阵求逆结果。

本发明的优点如下：

(1)本发明可以实现任意阶可逆矩阵的求逆运算，可以通过增减存储资源和增减运算资源，满足不同阶数的性能需求。

(2)本发明将矩阵求逆运算分为列选主元LU分解、下三角矩阵求逆、矩阵乘法模块三步骤执行，通过一个状态机控制每步的执行，每个步骤都可以有效地并行计算，且具有可扩展性。

(3)本发明有效的重复利用运算和存储资源，相对于其他求逆方法有效地减少运算的资源需求量，通过分时复用的方式最大限度地提高硬件资源的利用率以及运算效率。

(4)本发现在LU分解过程中通过列选主元方式，有效地提高了矩阵求逆的计算精度和正确性，在基于单精度浮点运算单元条件下，128阶矩阵求逆的误差矩阵均方根统计误差(RMS)小于10^-4。

附图说明

图1为矩阵求逆整体架构示意图。

图2为4×4的矩阵进行列选主元LU分解示意图。

图3为列选主元LU分解流程图。

具体实施方式

下面结合本发明实施例中的附图，对本发明的矩阵求逆运算方法进行详细的描述。

矩阵求逆分三步骤进行，分别是列选主元LU分解、三角矩阵求逆、矩阵相乘，使用一个状态机控制先后顺序，并且通过控制复选器，分时复用运算资源和存储资源，整体架构如图1。

(1)列选主元LU分解

计算矩阵A的部分选主元LU分解：PA＝LU，其中P为置换矩阵，L为单位下三角矩阵(对角线元素为1)，U为上三角矩阵。置换矩阵P为单位矩阵的行重新排列后的矩阵，每行每列有且仅有一个为1的元素，所以只需使用一个向量表示，向量中的第i个元素为置换矩阵中第i行中唯一一个1的列数。

下面举例介绍4×4的矩阵A上进行部分选主元LU分解的过程。如图2所示，首先从第一列中选取最大主元u₁₁，与该列的第一个元素的值进行交换，并存储最大主元的行下标i；然后进行归一化，用u₁₁去除该列其它元素，得到L分解因子的一列元素l₂₁、l₃₁和l₄₁。后续每一列第一个元素的值与行下标i指向的元素的值进行交换，分别得到U分解因子的一行元素u₁₂、u₁₃和u₁₄，同时得到新的子矩阵a_ij(2≤i,j≤4)，对子矩阵进行更新操作：a¹ _ij＝a_ij-l_i1*u_1j。然后再对新的子矩阵重复以上步骤，最终得到L矩阵和U矩阵。另外存储下来的最大主元行下标要对三角矩阵L和U求逆后的乘法结果进行列变换。

列选主元LU分解循环重复三个步骤：选主元、归一化、数据更新。对于一个n阶矩阵，需要循环n-1次，k表示当前循环次数，如图3所示。

I.选主元：根据当前循环次数k，将矩阵的主列元素(第k列对角线及以下)读出，通过两个实数乘法器和一个实数加法器，流水运算得出主列元素的实部虚部平方和，再通过一个实数减法器，找出主元(即第k列对角线及以下元素中模最大的点)。具体实现过程是：首先读出两个平方和结果送到实数减法器，通过判断减法结果浮点数的符号位，确定并记录较大数及其行下标后，读入下一个平方和结果与当前最大数进行比较，直至找出当前循环中n-k+1个主列元素中主元(即模最大的数)。每次循环找出的主元行下标即组成置换矩阵P。

II.归一化：根据当前循环中找到的最大主元行下标i，将源矩阵的第k行与第i行互换。随后将数据交换后的主列(第k列)对角线数的实部和虚部，分别送入两个实数除法器算出主元的倒数，再将主列(第k列)对角线以下元素读出，通过一个复数乘法器与主元的倒数相乘(即归一化操作)。

III.数据更新：数据更新操作是列选主元LU中计算量最大的一步，针对此特点本设计采用四路并行运算，每路运算单元使用一个复数乘法器和一个复数加法器，将矩阵第k行第k+1列、第k行第k+2列、第k行第k+3列、第k行第k+4列的数据，分别送入四路运算单元中复数乘法器，同时与归一化的结果进行复数乘法运算；复数乘法器的结果直接送入复数加法器，与矩阵第k+1列、第k+2列、第k+3列、第k+4列中第k+1行至第n行数据分别相减。当这四列数据更新完后，再对矩阵第k行第k+5列、第k行第k+6列、第k行第k+7列、第k行第k+8列并行的进行数据更新操作，依此类推直至矩阵的最后一列完成数据更新操作，即完成了一次循环运算。

每次循环完成之后，对k进行加1操作，进行下一次循环，重复执行上述三大步骤，直至k＝n-1，即得到矩阵分解的结果，上三角矩阵U和单位下三角矩阵L。

(2)三角矩阵求逆

对L矩阵求逆得到其逆矩阵L^-1，将U矩阵的转置矩阵求逆后再转置得到U^-1。对于分解后的下三角矩阵L的求逆，有以下的计算公式：

l_ij≠0(i＝1,2,…n-1；i＝j+1,j+2,…,n)

l_ij＝0(j＞i)的逆矩阵为其中：

$s_{ij}=\left\{\begin{array}{cc}\frac{1}{l_{ii}},& (i=j)\\ -\frac{1}{l_{ii}}\underset{k=j}{\overset{i-1}{\Σ}}{l}_{ik}{s}_{kj},& (j=1,2,...n-1;i=j+1,...n)\\ 0,& (j>i)\end{array}\right.---\left(2\right)$

上三角矩阵U的求逆过程可以利用U^-1＝((U^-1)^H)^H＝((U^H)^-1)^H，先转置成下三角矩阵求逆，然后再转置回U的逆矩阵。

本实施例中，根据上述步骤2)，提供一种下三角矩阵求逆运算模块，该模块划分为取倒单元、乘累加单元和取反乘单元。对于一个n阶下三角矩阵求逆，取倒单元只需1次，乘累加单元和取反乘单元需要循环n-1次。

A.取倒单元：根据下三角矩阵求逆公式，首先将对角线元素进行取倒。对于L矩阵，其对角线元素全为1，所以结果也都为1。对于转置后的U矩阵，将对角线元素读出，通过两个实数乘法器、一个实数加法器、一个实数除法器和一个复数乘法器得到其取倒的结果。具体实施如下：将复数的实部虚部分别输出到两个实数乘法器，两个实数乘法器的输出端分别与实数加法器的两个输入端连接，实数加法器的输出端与实数除法器的输入端连接，实数除法器的另一输入端固定地设置为1，最后实数除法器的输出端与复数乘法器的输入端连接。

B.乘累加单元：根据下三角矩阵求逆公式分析其数据依赖关系，通过四组由一个复数乘法器和一个复数加法器所组成的乘累加器，按行顺序并行计算非对角线元素，从第2行开始执行。L矩阵和转置后的U矩阵同时计算，分别使用两组乘累加器。乘累加器的具体实施如下：复数乘法器完成公式 $s_{ij}=-\left(\underset{k=j}{\overset{i-1}{\Σ}}{l}_{ik}{s}_{kj}\right)/l_{ii},(j=1,2,...,n-1;i=j+1,...,n)$ 中l_ik与S_kj的复数相乘操作，其结果通过复数加法器和延迟控制逻辑完成累加操作。l_ik表示L矩阵中的数，S_kj表示L矩阵的逆矩阵中的数，i表示行号，j表示列号，k表示累加循环次数。

C.取反乘单元：通过复数乘法器对前一步乘累加单元的结果，与取倒单元的结果(对应行的对角线元素倒数)进行取反乘运算，即得到相应行的求逆结果。对于L矩阵，因其对角线元素为1，只需将乘累加单元结果的符号位取反即为求逆结果。对于U矩阵，需将乘累加单元的结果取反后再与对应行的对角线元素相乘。

经过上述二、三步骤循环执行n-1次后，即可求得L、U矩阵的逆矩阵。

(3)矩阵相乘

最后A矩阵的求逆要通过计算公式A^-1＝U^-1*L^-1*P来完成。矩阵乘法模块采用四路并行运算，每路运算单元使用一个复数乘法器和两个复数加法器，其中两个复数加法器负责累加操作。整个乘法的过程当中，乘法器和累加器保持流水，累加器每隔n个周期出一个结果数据，其中n是矩阵的阶数。乘法的结果矩阵根据列选主元LU分解步骤中的置换矩阵P，对结果矩阵进行列变换得到最终的矩阵求逆结果。

本实施例根据上述技术方案(参见图1)提供了一种硬件模块。该硬件模块具体采用仿真/综合工具和40nmCMOS工艺，主频达到1GHz。基于1GHz时钟频率下，复数乘法器、复数加法器、实数乘法器、实数加法器均为4拍延迟的流水运算单元，实数除法器为19拍延迟的非流水运算单元，所有运算单元均为单精度浮点运算单元，存储器的取数延迟为6拍，128阶矩阵求逆的运行时间为1.42ms，结果误差矩阵均方根统计误差(RMS)小于10^-4，增加了整体的可靠性。另一方面，该技术方案是基于单端口Memory和单精度浮点运算单元的并行流水处理，并支持任意阶数矩阵求逆的硬件模块有效加快了运算速度，提高了硬件利用率，且具有良好的可扩展性，适合于VLSI的设计实现。

Claims (10)

1.一种矩阵求逆运算方法，其特征在于包括如下步骤：

1)首先进行列选主元LU分解：根据公式PA＝LU，将源矩阵A分解为单位下三角矩阵L、上三角矩阵U和置换矩阵P；

2)接着进行三角矩阵求逆：对L矩阵求逆得到其逆矩阵L^-1，将U矩阵的转置矩阵求逆后再转置得到U^-1；

3)最后进行矩阵相乘：将矩阵U^-1和矩阵L^-1相乘，并根据置换矩阵P将矩阵乘法结果进行列变换得到源矩阵A^-1。

2.根据权利1所述的矩阵求逆运算方法，其特征在于，步骤1)包括如下步骤：

I)进行选主元：从第一列中选取最大主元，与该列的第一个元素的值进行交换，并存储最大主元的行下标；

II)进行归一化：将第一行与最大主元所在行互换，得到U分解因子的一行元素u_1j，用交换后的第一列第一个元素依次除以该列的其它元素，得到L分解因子的一列元素l_i1，其中1≤i≤n，1≤j≤n，n表示矩阵阶数；

III)进行数据更新：对矩阵进行更新操作：a¹ _ij＝a_ij-l_i1*u_1j，得到新子矩阵a¹ _ij，每次更新操作后的子矩阵的阶数比原矩阵的阶数少一阶，其中，a_ij为原矩阵中的值，2≤i≤n，2≤j≤n，n表示矩阵阶数。对新的子矩阵循环重复以上三个步骤，最终得到L矩阵和U矩阵。

3.根据权利要求2所述的矩阵求逆运算方法，其特征在于，所述步骤I)中根据当前循环次数k，将矩阵的主列元素读出，先将所述主列元素通过两个实数乘法器和一个实数加法器，流水运算得出主列元素的实部虚部平方和，再通过一个实数减法器，找出主元。

4.根据权利要求2所述的矩阵求逆运算方法，其特征在于，所述步骤II)中根据当前循环中找到的最大主元行下标，将源矩阵的第k行与最大主元所在行互换，接着将数据交换后的主列对角线数的实部和虚部分别送入两个实数除法器中用于算出主元的倒数，再将主列对角线以下元素读出，通过一个复数乘法器与主元的倒数相乘，完成归一操作，其中k表示当前循环次数。

5.根据权利要求2所述的矩阵求逆运算方法，其特征在于，所述步骤III)中采用四路并行运算，每路运算使用一个复数乘法器和一个复数加法器，将矩阵第k行第k+1列、第k行第k+2列、第k行第k+3列、第k行第k+4列的数据，分别送入四路运算的对应复数乘法器中，与步骤2)中归一化的结果进行复数乘法运算；复数乘法器的结果直接送入复数加法器，与矩阵第k+1列、第k+2列、第k+3列、第k+4列中第k+1行至第n行数据分别相减，当所述四列数据更新完后，再对所述四列数据后续四列数据进行数据更新操作，依此类推直至矩阵的最后一列完成数据更新操作，即完成了一次数据更新，其中k表示当前循环次数，n为矩阵阶数。

6.根据权利要求2所述的矩阵求逆运算方法，其特征在于，所述步骤2)通过取倒单元、乘累加单元以及取反乘单元实现矩阵的求逆、转置，其中，

A.取倒单元：根据下三角矩阵求逆公式，首先将对角线元素进行取倒，对于L矩阵，结果为1；对于转置后的U矩阵，将对角线元素读出，通过两个实数乘法器、一个实数加法器、一个实数除法器和一个复数乘法器得到取倒运算的结果；

B.乘累加单元：根据下三角矩阵求逆公式分别分析L矩阵和转置后的U矩阵的数据依赖关系，通过四组由一个复数乘法器和一个复数加法器所组成的乘累加器进行乘累加运算，所述乘累加运算从第2行开始按行顺序计算非对角线元素，L矩阵和转置后的U矩阵同时进行乘累加运算；

C.取反乘单元：通过复数乘法器对步骤B)中乘累加单元的结果，与步骤A)中取倒单元的结果进行取反乘运算，即得到相应行的取反结果。

7.根据权利要求6所述的矩阵求逆运算方法，其特征在于，所述步骤A)中将复数的实部虚部分别输出到所述两个实数乘法器，两个实数乘法器的输出端分别与所述实数加法器的两个输入端连接，实数加法器的输出端与所述实数除法器的输入端连接，实数除法器的另一输入端固定地设置为1，实数除法器的输出端与所述复数乘法器的输入端连接。

8.根据权利要求6所述的矩阵求逆运算方法，其特征在于，所述步骤B)中L矩阵和转置后的U矩阵进行乘累加运算时，分别使用两组乘累加器，每组乘累加器包括

复数乘法器完成L、U递推公式 $s_{ij}=-\left(\underset{k=j}{\overset{i-1}{\Σ}}{l}_{ik}{s}_{kj}\right)/l_{ii},(j=1,2,...,n-1;i=j+1,...,n)$ 中l_ik与S_kj的复数相乘操作，l_ik表示L矩阵中的数，S_kj表示L矩阵的逆矩阵中的数，i表示行号，j表示列号，k表示累加循环次数；

复数加法器和延迟控制逻辑，根据复数乘法器的结果完成累加操作。

9.根据权利要求6所述的矩阵求逆运算方法，其特征在于，所述步骤C)中对于L矩阵，将乘累加单元结果的符号位取反作为取反结果；对于U矩阵，需将乘累加单元的结果取反后再与对应行的对角线元素相乘作为取反结果。

10.根据权利要求1所述的矩阵求逆运算方法，其特征在于，所述步骤3)中将U^-1矩阵按行划分为四等份，同时使用四路乘累加运算单元分别与L^-1矩阵进行乘累加运算，并根据步骤1)中的P矩阵，对乘法结果矩阵进行列变换得到最终的矩阵求逆结果。