发明内容
发明目的:本发明目的在于提供一种基于改进Givens旋转的矩阵求逆方法,去掉Givens旋转过程中的平方根、除法运算,大大降低算法复杂度。
为更好的理解本发明的技术方案,首先对技术理论基础说明如下:对于矩阵A进行QR分解,可以分解为一个酉矩阵和一个上三角矩阵的乘积,即A=QR,其中Q是一个酉矩阵,R是一个上三角矩阵,使用SGR(Squared Givens Rotations)方法,矩阵A的QR分解可以变形为:
(公式1)
其中QA=QDR,DR=diag(R),U=DRR为上三角矩阵,则有:
(公式2)
本发明的主要目标是提供一种省略平方根运算的同时减少除法运算的低复杂度并行矩阵求逆方法,主要是获得上三角矩阵U及相应的逆矩阵和
技术方案:本发明提供一种低复杂度的快速并行矩阵求逆方法,应用于MIMO通信系统接收端的信道估计和接收端的信号均衡处理,对于一个发送天线数为M,接收天线数为N的MIMO通信系统,其接收机信号可以表示为r=Hs+n,其中,r表示接收信号,是维数为M的列向量;s表示发送信号,是维数为N的列向量;H表示信道矩阵,是维数为M×N的矩阵;n表示加性高斯白噪声,是维数为M的列向量。接收端的均衡器从接收到的信号r估计出发送信号s,常见的均衡算法包括迫零算法和最小均方误差算法。基于迫零算法的均衡表达式为其中表示基于迫零算法的均衡器对发送信号的估计量;基于最小均方误差算法的均衡表达式为其中表示基于最小均方误差算法的均衡器对发送信号的估计量,上标H表示矩阵的共轭转置,IM表示M维单位矩阵,表示噪声的平均功率。将待求逆矩阵(HHH)或矩阵记为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的维数为M1×M2,如果M1<2,则执行步骤1.6;如果M1≥2,对矩阵B调用矩阵消零模块,矩阵消零模块输出矩阵
(公式5)
其中符号×表示非零元素值,符号表示用向量标量对表示矩阵的各行元素,即用<m,n>表示矩阵的第1行元素,用<pi,qi>表示矩阵的第i行元素,i=2,3,…,M1,各向量标量对的取值由矩阵消零模块计算得到。
步骤1.3:定义向量用来保存向量m,定义标量用来保存n,即令
步骤1.4:将矩阵B重新赋值且维数减1,即M1=M1-1,M2=M2-1且令表示矩阵B由矩阵的第2行到第M1行以及第2列到第M2列元素组成。
步骤1.5:对矩阵B重复步骤1.2到步骤1.4。
步骤1.6:执行完上述步骤,矩阵B的维数为1×(N+1),定义向量hN,令hN=B*(1,1)·B(1,1:N+1),其中B(1,1)表示矩阵B的第一行第一列位置的元素值,B(1,1:N)表示矩阵B的第一行元素,定义标量gN=1。
步骤1.7:令hx(j)表示向量hx的第j个元素,x=1,2,…,N,j=1,2,…,2N-x+1,则上三角矩阵U和可以分别表示成:
(公式6)
所述步骤2采用回代法求上三角矩阵U的逆矩阵U-1的具体方法为:
上三角矩阵U的逆矩阵也是上三角矩阵,U-1的形式可以表示为:
(公式7)
采用下面的回代公式,按照w11,w12,…w1n,w22,…,w2n,w33,…,wnn的顺序求出W中各个元素的值:
(公式8)
上述步骤1.2中的矩阵消零模块对复数矩阵B进行消零运算的具体方法为:
对任意M1×M2维复数矩阵B,记为:
其中bk表示矩阵B的第k行元素,pk表示维数与bk相同的行向量,qk表示标量,符号表示分别用向量标量对<pk,qk>表示矩阵B的第k行;pk初始化为矩阵B的第k行元素组成的行向量,即pk=bk,qk初始化为1,k=1,2,…,M1。
定义向量s,t为M1维列向量且元素初始值为s(l)=1和t(l)=1,其中s(l)和t(l)分别表示向量s和t的第l个元素,l=1,2,…,M1。
对复数矩阵B进行消零运算的具体步骤为:
步骤1.2.1:旋转行向量b1和b2将b21消为0。
1)将b1用向量m和标量n表示为对m,n,进行初始化:p1(r)表示向量p1的第r个元素,如果p1(1)=0,则m=s(1)p1,n=t(1)q1,如果p1(1)≠0,则m=s(1)p1(1)*p1,
2)使用MSGR旋转,对各变量进行更新:约定表示x更新后的值,则向量b1和b2更新为:下面分情况给出以及和的表达式:
当m(1)≠0时
当m(1)=0,且p2(1)=0时
当m(1)=0,且p2(1)≠0时 (公式9)
3)溢出处理:exp(c)表示浮点数c的指数值,令x,y为两个复数,约定对于复数num,real(num)表示该复数的实部,imag(num)表示该复数的虚部,则x和y的实部和虚部可以分别用x1,x2,y1,y2表示如下:x1=real(x),x2=imag(x),y1=real(y),y2=imag(y),令约定语句IFTHEN如果条件A成立,则执行B操作。
对复数对<x,y>进行溢出处理可以用以下语句表示:
IFTHEN
IFTHEN
根据上述方法分别对复数对(或复数向量对)进行溢出处理。
4)将更新后的以及和分别赋值给m,n,p2,q2,以及s(2)和t(2),即令
步骤1.2.2:将b31,b41,……,位置元素消为0。
将b2,p2,q2,s(2),t(2)替换为br,pr,qr,s(r),t(r),r=3,4,…,M1,重复步骤1.2.1中2)到4)的操作,可以依次将b31,b41,……,位置元素消为0。
为降低算法复杂度,在不过多损失算法精度的情况下,整个矩阵消零模块可以省略关于向量s和t的运算,此时,公式(9),用如下公式替换:
当m(1)≠0时
当m(1)=0,且p2(1)=0时 (公式10)
当m(1)=0,且p2(1)≠0时
其他步骤涉及到s和t的操作不做即可。
有益效果:本发明提出的低复杂度并行矩阵求逆方法,主要在两个方面提高了算法效率:第一,改进的Givens旋转完全避免了平方根运算,大大降低算法复杂度;第二,改进的Givens旋转,采用存储分子分母的方式,节省了大量的除法运算,降低了算法复杂度。这两点对于矩阵求逆的算法效率都有显著提高,适用于无线通信、信号处理以及数值计算等领域的矩阵求逆问题。
具体实施方式
下面以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)
其中bk表示矩阵B的第k行元素,pk表示维数与bk相同的行向量,qk表示标量,符号表示用向量标量对<pk,qk>表示矩阵B的第k行;pk初始化为矩阵B的第k行元素组成的行向量,即pk=bk,qk初始化为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:将b21位置元素消为0。
1)将b1用向量m和标量n表示为对m、n进行初始化:p1(r)表示向量p1的第r个元素,如果p1(1)=0,则m=s(1)p1,n=t(1)q1,如果p1(1)≠0,则m=s(1)p1(1)*p1,
2)使用MSGR旋转,对各变量进行更新:约定表示x更新后的值,则向量b1,b2更新为:下面分情况给出以及和的表达式:
当m(1)≠0时
当m(1)=0,且p2(1)=0时
当m(1)=0,且p2(1)≠0时 (公式13)
3)溢出处理:exp(v)表示浮点数v的指数值,令x,y为两个复数,约定对于复数num,real(num)表示该复数的实部,imag(num)表示该复数的虚部,则x和y的实部和虚部可以分别用x1,x2,y1,y2表示如下:x1=real(x),x2=imag(x),y1=real(y),y2=imag(y),令约定语句IFTHEN如果条件A成立,则执行B操作。
对复数对<x,y>进行溢出处理可以用如下语句表示:
IFTHEN
IFTHEN
根据上述方法分别对复数对(或复数向量对) 进行溢出处理。
4)将更新后的以及和分别赋值给m,n,p2,q2,以及s(2)和t(2),即令
步骤1.2.2:将b31位置元素消为0。
将b2,p2,q2,s(2),t(2)替换为b3,p3,q3,s(3),t(3),重复步骤1.2.1中2)到4)的操作,可以将b31位置元素消为0。
步骤1.2.3:定义向量h1,用来保存向量m,定义标量g1用来保存n,
即令h1=m,g1=n。
步骤1.3:对更新后的矩阵B继续调用矩阵消零模块。
经过步骤1.2,矩阵B更新为
令
步骤1.3.1:对B调用矩阵消零模块,将位置元素消为零。
对矩阵B重复步骤1.2.1操作,可以将b32位置元素消为零。
步骤1.3.2:定义向量h2,用来保存向量m,定义标量g2用来
保存n即令h2=m,g2=n。
步骤1.4:经过步骤1.3,矩阵B更新为
令
定义向量h3,且同时定义一个标量g3=1。
步骤1.5:上三角矩阵U和可以分别表示成:
(公式14)
步骤2:回代法求上三角矩阵U的逆矩阵U-1。
上三角矩阵U的逆矩阵也是上三角矩阵,U-1的形式可以表示为:
(公式15)
采用下面的回代公式,按照w11,w12,w13,w22,w23,w33的顺序求出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随机产生的矩阵)