CN104298649B

CN104298649B - 一种低复杂度的快速并行矩阵求逆方法

Info

Publication number: CN104298649B
Application number: CN201410495246.0A
Authority: CN
Inventors: 何世文; 余登高; 黄永明; 王海明; 杨绿溪; 张军
Original assignee: In Jiangsu Emerging Micro-Communication Ceases Science And Technology Ltd
Current assignee: Nanjing etaco Communication Technology Co.,Ltd.
Priority date: 2014-09-24
Filing date: 2014-09-24
Publication date: 2017-05-17
Anticipated expiration: 2034-09-24
Also published as: CN104298649A

Abstract

本发明公开了一种低复杂度的快速并行矩阵求逆方法，该发明包括以下步骤：第一步，对于给定矩阵A，令矩阵E为与矩阵A同阶的单位矩阵，将矩阵A与矩阵E组成扩展矩阵B，对矩阵B使用改进的Givens旋转(MSGR，Modified Squared Givens Rotations)，求得上三角矩阵U和其中矩阵U和的定义如下：使用SGR(Squared Givens Rotations)方法对矩阵A的QR分解进行变形与原始QR分解的关系为Q_A＝QD_R，D_R＝diag(R)，U＝D_RR，D_U＝diag(D_RR)。MSGR方法可以去掉Givens旋转过程中的平方根运算同时减少除法运算，显著降低算法复杂度；第二步，回代法求上三角矩阵U的逆矩阵U^‑1；第三步，矩阵求逆本发明去掉了大量的除法以及平方根运算，降低了算法复杂度，可以应用于无线通信、信号处理以及数值计算等领域的矩阵求逆问题。

Description

一种低复杂度的快速并行矩阵求逆方法

技术领域

本发明涉及一种去平方根运算、减少除法运算的低复杂度的快速并行矩阵求逆方法，主要涉及无线通信、信号处理以及数值计算等领域。

背景技术

随着无线通信技术的不断发展，用户对通信系统的可靠性和有效性要求越来越高。为了提高系统的传输速率、增加系统的频谱效率，多输入多输出(Multiple InputMultiple Output，MIMO)技术作为一种关键技术得到了广泛研究。IEEE802.11n、IEEE802.11ac、IEEE 802.11ad以及3GPP-LTE等无线通信标准，都采用了MIMO技术。

MIMO通信系统，在接收端的设计相当复杂。接收端信道估计会涉及到大量的矩阵求逆运算，尤其随着发射天线以及接收天线数目的增加，矩阵的维数也随着增加，矩阵求逆的运算复杂度也会大量增加，导致矩阵求逆的算法设计以及硬件实现难度大大提高，这也成为当前MIMO系统亟待解决的一个技术难题。因此，提出一个有效的矩阵求逆算法，可以提高通信系统的效率、降低硬件实现的成本，对于通信系统的发展具有很重要的现实意义。

Givens旋转在数值线性代数中的主要用途是在向量或矩阵中介入零，常用于矩阵的QR分解，Givens旋转具有稳定性高、易于映射到脉动阵列、便于并行处理等优点，改进的Givens旋转，去掉Givens旋转过程中的平方根以及除法运算有利于降低算法复杂度。

发明内容

发明目的：本发明目的在于提供一种基于改进Givens旋转的矩阵求逆方法，去掉Givens旋转过程中的平方根、除法运算，大大降低算法复杂度。

为更好的理解本发明的技术方案，首先对技术理论基础说明如下：对于矩阵A进行QR分解，可以分解为一个酉矩阵和一个上三角矩阵的乘积，即A＝QR，其中Q是一个酉矩阵，R是一个上三角矩阵，使用SGR(Squared Givens Rotations)方法，矩阵A的QR分解可以变形为：

(公式1)

其中Q_A＝QD_R，D_R＝diag(R)，U＝D_RR为上三角矩阵，则有：

(公式2)

本发明的主要目标是提供一种省略平方根运算的同时减少除法运算的低复杂度并行矩阵求逆方法，主要是获得上三角矩阵U及相应的逆矩阵和

技术方案：本发明提供一种低复杂度的快速并行矩阵求逆方法，应用于MIMO通信系统接收端的信道估计和接收端的信号均衡处理，对于一个发送天线数为M，接收天线数为N的MIMO通信系统，其接收机信号可以表示为r＝Hs+n，其中，r表示接收信号，是维数为M的列向量；s表示发送信号，是维数为N的列向量；H表示信道矩阵，是维数为M×N的矩阵；n表示加性高斯白噪声，是维数为M的列向量。接收端的均衡器从接收到的信号r估计出发送信号s，常见的均衡算法包括迫零算法和最小均方误差算法。基于迫零算法的均衡表达式为其中表示基于迫零算法的均衡器对发送信号的估计量；基于最小均方误差算法的均衡表达式为其中表示基于最小均方误差算法的均衡器对发送信号的估计量，上标H表示矩阵的共轭转置，I_M表示M维单位矩阵，表示噪声的平均功率。将待求逆矩阵(H^HH)或矩阵记为A，E为与A同阶的单位矩阵，分别表示如下：

(公式3)

本发明所提供的方法主要包括以下步骤：

步骤1：采用MSGR(Modified Squared Givens Rotations)方法求解矩阵A按(公式1)分解时的上三角矩阵U以及

步骤2：采用回代法求上三角矩阵U的逆矩阵U^-1；

步骤3：将步骤2中获得的U^-1和步骤1中获得的相乘，得到矩阵A的逆矩阵。

所述步骤1的具体步骤如下：

步骤1.1：定义矩阵B为矩阵A与E组成的扩展矩阵，即B表示为：

(公式4)

步骤1.2：记矩阵B的维数为M₁×M₂，如果M₁＜2，则执行步骤1.6；如果M₁≥2，对矩阵B调用矩阵消零模块，矩阵消零模块输出矩阵

(公式5)

其中符号×表示非零元素值，符号表示用向量标量对表示矩阵的各行元素，即用<m,n>表示矩阵的第1行元素，用<p_i,q_i>表示矩阵的第i行元素，i＝2,3,…,M₁，各向量标量对的取值由矩阵消零模块计算得到。

步骤1.3：定义向量用来保存向量m，定义标量用来保存n，即令

步骤1.4：将矩阵B重新赋值且维数减1，即M₁＝M₁-1，M₂＝M₂-1且令表示矩阵B由矩阵的第2行到第M₁行以及第2列到第M₂列元素组成。

步骤1.5：对矩阵B重复步骤1.2到步骤1.4。

步骤1.6：执行完上述步骤，矩阵B的维数为1×(N+1)，定义向量h_N，令h_N＝B^*(1,1)·B(1,1:N+1)，其中B(1,1)表示矩阵B的第一行第一列位置的元素值，B(1,1:N)表示矩阵B的第一行元素，定义标量g_N＝1。

步骤1.7：令h_x(j)表示向量h_x的第j个元素，x＝1,2,…,N，j＝1,2,…,2N-x+1，则上三角矩阵U和可以分别表示成：

(公式6)

所述步骤2采用回代法求上三角矩阵U的逆矩阵U^-1的具体方法为：

上三角矩阵U的逆矩阵也是上三角矩阵，U^-1的形式可以表示为：

(公式7)

采用下面的回代公式，按照w₁₁,w₁₂,…w_1n,w₂₂,…,w_2n,w₃₃,…,w_nn的顺序求出W中各个元素的值：

(公式8)

上述步骤1.2中的矩阵消零模块对复数矩阵B进行消零运算的具体方法为：

对任意M₁×M₂维复数矩阵B，记为：

其中b_k表示矩阵B的第k行元素，p_k表示维数与b_k相同的行向量，q_k表示标量，符号表示分别用向量标量对<p_k,q_k>表示矩阵B的第k行；p_k初始化为矩阵B的第k行元素组成的行向量，即p_k＝b_k，q_k初始化为1，k＝1,2,…,M₁。

定义向量s，t为M₁维列向量且元素初始值为s(l)＝1和t(l)＝1，其中s(l)和t(l)分别表示向量s和t的第l个元素，l＝1,2,…,M₁。

对复数矩阵B进行消零运算的具体步骤为：

步骤1.2.1：旋转行向量b₁和b₂将b₂₁消为0。

1)将b₁用向量m和标量n表示为对m，n，进行初始化：p₁(r)表示向量p₁的第r个元素，如果p₁(1)＝0，则m＝s(1)p₁，n＝t(1)q₁，如果p₁(1)≠0，则m＝s(1)p₁(1)^*p₁，

2)使用MSGR旋转，对各变量进行更新：约定表示x更新后的值，则向量b₁和b₂更新为：下面分情况给出以及和的表达式：

当m(1)≠0时

当m(1)＝0，且p₂(1)＝0时

当m(1)＝0，且p₂(1)≠0时 (公式9)

3)溢出处理：exp(c)表示浮点数c的指数值，令x，y为两个复数，约定对于复数num，real(num)表示该复数的实部，imag(num)表示该复数的虚部，则x和y的实部和虚部可以分别用x₁，x₂，y₁，y₂表示如下：x₁＝real(x)，x₂＝imag(x)，y₁＝real(y)，y₂＝imag(y)，令约定语句IFTHEN如果条件A成立，则执行B操作。

对复数对<x,y>进行溢出处理可以用以下语句表示：

IFTHEN

根据上述方法分别对复数对(或复数向量对)进行溢出处理。

4)将更新后的以及和分别赋值给m，n，p₂，q₂，以及s(2)和t(2)，即令

步骤1.2.2：将b₃₁，b₄₁，……，位置元素消为0。

将b₂，p₂，q₂，s(2)，t(2)替换为b_r，p_r，q_r，s(r)，t(r)，r＝3,4,…,M₁，重复步骤1.2.1中2)到4)的操作，可以依次将b₃₁，b₄₁，……，位置元素消为0。

为降低算法复杂度，在不过多损失算法精度的情况下，整个矩阵消零模块可以省略关于向量s和t的运算，此时，公式(9)，用如下公式替换：

当m(1)≠0时

当m(1)＝0，且p₂(1)＝0时 (公式10)

当m(1)＝0，且p₂(1)≠0时

其他步骤涉及到s和t的操作不做即可。

有益效果：本发明提出的低复杂度并行矩阵求逆方法，主要在两个方面提高了算法效率：第一，改进的Givens旋转完全避免了平方根运算，大大降低算法复杂度；第二，改进的Givens旋转，采用存储分子分母的方式，节省了大量的除法运算，降低了算法复杂度。这两点对于矩阵求逆的算法效率都有显著提高，适用于无线通信、信号处理以及数值计算等领域的矩阵求逆问题。

附图说明

图1为本发明的算法流程图；

图2为本发明的MSGR变换过程示意图；

图3为本发明的矩阵消零模块变换示意图。

具体实施方式

下面以3×3矩阵求逆过程为例，结合附图和具体实施方式对本发明做进一步说明。应理解该实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利。

本实施例中，假设给定矩阵为：

其相应的同阶单位矩阵为

如图1所示，一种低复杂度的快速并行矩阵求逆方法包括如下步骤：

步骤1：对于给定矩阵A，将矩阵A与矩阵E组成扩展矩阵B，对矩阵B采用MSGR方法求解矩阵A按(公式1)分解时的上三角矩阵U以及MSGR变换过程示意图如图2所示，具体操作如下：

步骤1.1：定义矩阵B，令矩阵B为矩阵A与E组成的扩展矩阵，即B表示为：

(公式11)

步骤1.2：对矩阵B调用矩阵消零模块进行矩阵消零运算，其变换过程示意图如图3所示，具体操作如下：

将矩阵B记为：

(公式12)

其中b_k表示矩阵B的第k行元素，p_k表示维数与b_k相同的行向量，q_k表示标量，符号表示用向量标量对<p_k,q_k>表示矩阵B的第k行；p_k初始化为矩阵B的第k行元素组成的行向量，即p_k＝b_k，q_k初始化为1，k＝1,2,3。

定义向量s，t为M维列向量且元素初始值为s(l)＝1和t(l)＝1，其中s(l)和t(l)分别表示向量s和t的第l个元素，l＝1,2,3。

步骤1.2.1：将b₂₁位置元素消为0。

1)将b₁用向量m和标量n表示为对m、n进行初始化：p₁(r)表示向量p₁的第r个元素，如果p₁(1)＝0，则m＝s(1)p₁，n＝t(1)q₁，如果p₁(1)≠0，则m＝s(1)p₁(1)^*p₁，

2)使用MSGR旋转，对各变量进行更新：约定表示x更新后的值，则向量b₁，b₂更新为：下面分情况给出以及和的表达式：

当m(1)≠0时

当m(1)＝0，且p₂(1)＝0时

当m(1)＝0，且p₂(1)≠0时 (公式13)

3)溢出处理：exp(v)表示浮点数v的指数值，令x，y为两个复数，约定对于复数num，real(num)表示该复数的实部，imag(num)表示该复数的虚部，则x和y的实部和虚部可以分别用x₁，x₂，y₁，y₂表示如下：x₁＝real(x)，x₂＝imag(x)，y₁＝real(y)，y₂＝imag(y)，令约定语句IFTHEN如果条件A成立，则执行B操作。

对复数对<x,y>进行溢出处理可以用如下语句表示：

IFTHEN

根据上述方法分别对复数对(或复数向量对) 进行溢出处理。

步骤1.2.2：将b₃₁位置元素消为0。

将b₂，p₂，q₂，s(2)，t(2)替换为b₃，p₃，q₃，s(3)，t(3)，重复步骤1.2.1中2)到4)的操作，可以将b₃₁位置元素消为0。

步骤1.2.3：定义向量h₁，用来保存向量m，定义标量g₁用来保存n，

即令h₁＝m，g₁＝n。

步骤1.3：对更新后的矩阵B继续调用矩阵消零模块。

经过步骤1.2，矩阵B更新为

令

步骤1.3.1：对B调用矩阵消零模块，将位置元素消为零。

对矩阵B重复步骤1.2.1操作，可以将b₃₂位置元素消为零。

步骤1.3.2：定义向量h₂，用来保存向量m，定义标量g₂用来

保存n即令h₂＝m，g₂＝n。

步骤1.4：经过步骤1.3，矩阵B更新为

令

定义向量h₃，且同时定义一个标量g₃＝1。

步骤1.5：上三角矩阵U和可以分别表示成：

(公式14)

步骤2：回代法求上三角矩阵U的逆矩阵U^-1。

(公式15)

采用下面的回代公式，按照w₁₁,w₁₂,w₁₃,w₂₂,w₂₃,w₃₃的顺序求出W中各个元素的值：

(公式16)

步骤3：矩阵求逆。

将步骤2获得的U^-1和步骤1获得的相乘，即从而完成整个矩阵求逆过程。

表1，以3×3复数矩阵为例列出了基于MSGR方法的消零模块的复杂度；表2，针对Matlab随机产生的矩阵，对本发明的求逆结果与Matlab自带函数求逆结果进行了对比，可以发现误差数量级都在10^-15以上。

表1 消零模块复杂度说明(3×3矩阵为例)

算法	乘法次数	加法次数	除法次数
				算法一：包含操作因子s,t	627	291	3
算法二：不包含操作因子s,t	474	285	3

表2 矩阵求逆实例(输入矩阵为Matlab随机产生的矩阵)

Claims

1.一种低复杂度的快速并行矩阵求逆方法，应用于MIMO通信系统接收端的信道估计和接收端的信号均衡，对于一个发送天线数为M，接收天线数为N的MIMO通信系统，其接收机信号可以表示为r＝Hs+n，其中，r表示接收信号，是维数为M的列向量；s表示发送信号，是维数为N的列向量；H表示信道矩阵，是维数为M×N的矩阵；n表示加性高斯白噪声，是维数为M的列向量；接收端的均衡器从接收到的信号r估计出发送信号s，均衡算法包括迫零算法和最小均方误差算法，基于迫零算法的均衡表达式为其中表示基于迫零算法的均衡器对发送信号的估计量；基于最小均方误差算法的均衡表达式为其中表示基于最小均方误差算法的均衡器对发送信号的估计量，上标H表示矩阵的共轭转置，I_M表示M维单位矩阵，表示噪声的平均功率，将待求逆矩阵(H^HH)或矩阵记为A，其特征在于，所述矩阵求逆方法将矩阵A的QR分解变形为其中Q_A＝QD_R，D_R＝diag(R)，U＝D_RR，Q是一个酉矩阵，R是一个上三角矩阵，所述矩阵求逆包括如下步骤：

(1)将矩阵A与同阶单位矩阵E组成的扩展矩阵的各行元素分别用向量标量对表示，所述向量标量对表示形式为<a,b>，其中a表示向量，b表示标量，基于向量标量对进行迭代消零运算来求解矩阵A按分解的上三角矩阵U以及

(2)计算上三角矩阵U的逆矩阵U^-1；

(3)将步骤(2)中获得的U^-1和步骤(1)中获得的相乘，得到矩阵A的逆矩阵。

2.根据权利要求1所述的低复杂度的快速并行矩阵求逆方法，其特征在于：所述步骤(1)求解矩阵A按分解的上三角矩阵U以及的具体步骤包括：

(1.1)定义矩阵B为矩阵A与E组成的扩展矩阵，E为与A同阶的单位矩阵，矩阵B表示为：

(1.2)记矩阵B的维数为M₁×M₂，如果M₁＜2，则执行步骤(1.6)；如果M₁≥2，对矩阵B进行消零运算，输出矩阵

其中符号×表示非零元素值，符号表示用向量标量对表示矩阵的各行元素，即用<m,n>表示矩阵的第1行元素，用<p_i,q_i>表示矩阵的第i行元素，i＝2,3,…,M₁；

(1.3)定义向量令定义标量令

(1.4)将矩阵B重新赋值，令并将维数减1，即M₁＝M₁-1，M₂＝M₂-1；

(1.5)对矩阵B重复步骤(1.2)到步骤(1.4)；

(1.6)执行完上述步骤，矩阵B的维数为1×(N+1)，定义向量h_N，令h_N＝B^*(1,1)·B(1,1:N+1)，其中B(1,1)表示矩阵B的第一行第一列位置的元素值，B(1,1:N)表示矩阵B的第一行元素，定义标量g_N＝1；

(1.7)令h_x(j)表示向量h_x的第j个元素，x＝1,2,…,N，j＝1,2,…,2N-x+1，则得到上三角矩阵

和

3.根据权利要求2所述的低复杂度的快速并行矩阵求逆方法，其特征在于，所述步骤(1.2)中对矩阵B进行消零运算的具体步骤包括：

(1.2.1)将b₁用向量m和标量n表示为对m，n，进行初始化：p₁(r)表示向量p₁的第r个元素，如果p₁(1)＝0，则m＝s(1)p₁，n＝t(1)q₁，如果p₁(1)≠0，则m＝s(1)p₁(1)^*p₁，其中b_k表示矩阵B的第k行元素，p_k表示维数与b_k相同的行向量，向量s，t为M₁维列向量，元素初始值为1；

(1.2.2)使用改进的Givens旋转，按照公式9对各变量进行更新；

当m(1)≠0时

当m(1)＝0，且p₂(1)＝0时

当m(1)＝0，且p₂(1)≠0时

向量b₁和b₂更新为：

(1.2.3)分别对复数对进行溢出处理；

(1.2.4)将更新后的以及和分别赋值给m，n，p₂，q₂，以及s(2)和t(2)；

(1.2.5)如果M₁＞2，将b₂，p₂，q₂，s(2)，t(2)依次替换为b_r，p_r，q_r，s(r)，t(r)，r＝3,4,…,M₁，重复步骤(1.2.2)到(1.2.4)的操作，依次将b₃₁，b₄₁，……，元素消为0。

4.根据权利要求2所述的低复杂度的快速并行矩阵求逆方法，其特征在于，所述步骤(1.2)中对矩阵B进行消零运算的具体步骤包括：

(1.2.1)将b₁用向量m和标量n表示为对m，n，进行初始化：p₁(r)表示向量p₁的第r个元素，如果p₁(1)＝0，则m＝p₁，n＝q₁，如果p₁(1)≠0，则m＝p₁(1)^*p₁，其中b_k表示矩阵B的第k行元素，p_k表示维数与b_k相同的行向量；

(1.2.2)使用改进的Givens旋转，按照公式10对各变量进行更新；

当m(1)≠0时

当m(1)＝0，且p₂(1)＝0时

当m(1)＝0，且p₂(1)≠0时

向量b₁和b₂更新为：

(1.2.3)分别对复数对进行溢出处理；

(1.2.4)将更新后的分别赋值给m，n，p₂，q₂；

(1.2.5)将b₂，p₂，q₂依次替换为b_r，p_r，q_r，r＝3,4,…,M₁，重复步骤(1.2.2)到(1.2.4)的操作，依次将b₃₁，b₄₁，……，元素消为0。

5.根据权利要求3或4所述的低复杂度的快速并行矩阵求逆方法，其特征在于，所述步骤(1.2.3)中对复数对进行溢出处理的方法为：设x，y为两个复数，x和y的实部和虚部分别用x₁，x₂，y₁，y₂表示，令对复数对<x,y>进行溢出处理的方法为：

如果z＞4，则令tem＝exp(z)mod2，

exp(x₁)＝exp(x₁)-(tem-1)，exp(x₂)＝exp(x₂)-(tem-1)，

exp(y₁)＝exp(y₁)-(tem-1)，exp(y₂)＝exp(y₂)-(tem-1)；

如果z≤0.25，则令tem＝-|exp(z)|mod2，

exp(x₁)＝exp(x₁)+(tem+1)，exp(x₂)＝exp(x₂)+(tem+1)，

exp(y₁)＝exp(y₁)+(tem+1)，exp(y₂)＝exp(y₂)+(tem+1)，

其中exp()表示浮点数的指数值。

6.根据权利要求1所述的低复杂度的快速并行矩阵求逆方法，其特征在于：所述步骤(2)计算上三角矩阵U的逆矩阵U^-1的具体方法为：对于

采用回代公式

按照w₁₁,w₁₂,…w_1n,w₂₂,…,w_2n,w₃₃,…,w_nn的顺序求出W中各个元素的值。