CN102882579A - 一种用于多天线系统的并行矩阵求逆方法 - Google Patents

Abstract

本发明公开了一种用于多天线系统的并行矩阵求逆方法，用于多天线系统中接收端的检测与自适应传输，该发明包括以下步骤：第一步，求取旋转矩阵：以快速Givens旋转矩阵代替常规的Givens旋转矩阵，可使乘法运算量减半；采用极限、泰勒展开式的方法进一步降低开方运算的复杂度。第二步，进行旋转变换：在计算列向量的模值时采用迭代更新的方法，相较于直接计算向量模值，显著提高了算法效率。第三步，反复迭代。第四步，矩阵求逆。本发明具有很高的效率和很好的并行性，也可以应用于无线通信、信号处理以及数值计算等领域中的矩阵求逆问题。

Description

一种用于多天线系统的并行矩阵求逆方法

技术领域

本发明涉及一种用于多天线(MIMO)系统的并行矩阵求逆方法，尤其涉及MIMO系统检测及自适应传输时的快速矩阵求逆运算。

背景技术

移动用户数的增长及移动互联网业务的增加，对通信系统的容量和质量提了更高的要求。MIMO技术可以提供空分复用增益、分集增益及天线阵列增益，从而显著提高通信系统的频谱效率并改善通信质量。事实上，MIMO技术已经成为许多无线通信标准，如3GPP-LTE、IEEE 802.11n以及IEEE 802.16e等标准中最为关键的技术之一。因此，MIMO技术一直是当前无线通信技术领域研究的热点。

MIMO即多输入多输出，它的定义很简单：对于无线通信系统，如果其发射天线和接收天线都是多根的话，则就是一个MIMO系统。MIMO系统的数学模型涉及矩阵运算。尤其是MIMO系统接收端的检测和MIMO系统的自适应传输，均涉及大量矩阵求逆运算。因此用于MIMO系统中的矩阵求逆算法的好坏，将极大的影响无线通信系统的成本和效率。而随着天线数目的增加，矩阵维数随之上升，矩阵求逆的算法复杂度更是迅速增加，这已成为MIMO技术领域中的难点与瓶颈。因此设计出更为高效的矩阵求逆算法，适应新一代宽带无线通信系统的应用需求，具有重要的现实意义。

在MIMO检测算法中，经典的矩阵求逆算法有QR分解算法和Jacobi算法。尽管QR分解算法收敛速度更快，但是由于Jacobi算法具有很高的并行性和更好的数值精度，因而Jacobi算法在矩阵求逆中应用更为广泛。Jacobi方法有两种类型，即单边Jacobi和双边Jacobi。传统的双边Jacobi算法被广泛应用于特征值求解问题中。然而，作为一种双边Jacobi方法，每次迭代都同时涉及行和列的更新，行和列之间存在很大的数据相关性，这种数据相关性在算法的并行实现时将带来很大的核间通信代价，降低算法效率。

发明内容

本发明提出了一种用于多天线系统的并行矩阵求逆方法，基于单边Jacobi，每次迭代只涉及列的数据更新，数据相关性大大降低，因而在算法并行实现时将不需要如此巨大的核间通信代价，因此各处理器之间可以更加独立的工作，算法的并行实现将更加容易，算法效率也将更高。

本发明采用的技术方案为：一种用于多天线系统的并行矩阵求逆方法，包括以下步骤：

1)一个发送天线数为M、接收天线数为N的多天线系统，其信号的输入输出关系可表示为r＝HWs+n，其中：r表示接收信号，是维度为M的列向量；s表示发送信号，是维度为L的列向量；W表示发送预编码矩阵，是维度为N×L的矩阵；H表示信道矩阵，是维度为M×N的矩阵；n表示加性高斯白噪声，是维度为M的列向量；接收端进行检测，从接收到信号r中估计出发送信号s；多天线系统主流的检测算法包括迫零(ZF)算法和最小均方误差(MMSE)算法，其检测器表达式可分别表示为 ${\hat{s}}_{\mathrm{ZF}}={\left(H_e^H{H}_e\right)}^{-1}{H}_e^{H}r,{\hat{s}}_{\mathrm{MMSE}}={(H_e^H{H}_e+{\mathrm{\σ}}_n^2{I}_M)}^{-1}{H}_e^{H}r,$ 其中，

分别表示检测器采用迫零算法、最小均方误差算法对发送信号的估计，上标H表示矩阵的共轭转置，I_M表示M维单位矩阵，

表示噪声的平均功率，H_e＝HW表示等效信道矩阵；

2)为了使多天线系统的传输信号能自适应于信道环境的变化，需要在发送端根据信道条件，在给定准则下自适应地选取发送预编码矩阵W，最大化MMSE和速率是一种常用的选取预编码矩阵的准则，其计算公式为

$R_{\mathrm{MMSE}}\left(W\right)=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{\log }_2\frac{1}{{\left[{({\mathrm{\σ}}_n^{-2}{H}_e^H{H}_e+I_M)}^{-1}\right]}_{j,j}}$

其中R_MMSE(W)表示预编码矩阵为W时接收端接收信号产生的均方误差，下标j，j表示矩阵的对角线元素。

3)计算旋转矩阵：

将上述步骤1)和步骤2)中多天线系统中待求逆的矩阵

或矩阵记为A，并记初始矩阵A₀＝A，逐次选择修正矩阵R_k，使A_k的i，j两列正交化。修正矩阵R_k采用快速Givens矩阵的形式：

R_k，A_k分别表示第k次迭代时的旋转矩阵和待修正的矩阵。

参数α计算如下：设I＝<a_i ^(k)，a_i ^(k)>，J＝<a_j ^(k)，a_j ^(k)>，e＝<a_i ^(k)，a_j ^(k)>，其中a_j ^(k)表示矩阵A_k的第j列的列向量。这里<x，y>表示向量x和y的内积，若e＝0，则说明两列已经正交，取α＝0，否则记

其中sgn(t)表示符号函数。α这样取值可以保证相应的旋转角满足

以提高算法的数值稳定性。注意到参数

这里涉及到开方运算，复杂度较高，考虑进行一定的近似，从而降低复杂度。注意到α关于t是奇对称的，所以仅需讨论t＞0的情形，即利用极限的知识以及Taylor(泰勒)展开式，易得：

$\underset{t\&RightArrow;\&infin;}{\lim }\sqrt{t^2+1}-t=\underset{t\&RightArrow;\&infin;}{\lim }\frac{1}{\sqrt{t^2+1}+t}=\underset{t\&RightArrow;\&infin;}{\lim }\frac{1}{2t},$

$\underset{t\&RightArrow;0}{\lim }\sqrt{t^2+1}-t\&ap;\underset{t\&RightArrow;0}{\lim }(1+\frac{t^2}{2})-t$

所以，在实际计算时，可作如下近似：

$\widetilde{\mathrm{\&alpha;}}=\left\{\begin{array}{cc}1-t+\frac{t^2}{2}& 0<t\&le;0.25\\ -t+\sqrt{t^2+1}& 0.25<t\&le;8\\ \frac{1}{2t}& 8<t\&le;1024\\ 0& 1024<t\end{array}\right.$

其中表示对α的近似。分段点的选取是随用户要求的精度变化而变化的，上述分段点的选取是假设控制误差在10^-4数量级内得到的。

4)进行旋转变换：A_k+1＝A_kR_k(k＝0，1，…)：

注意到在每次旋转变换时都需要计算第i，j列的模值，即I，J的值，常规的单侧Jacobi算法采用直接求取列向量模值的方法，共需要4n个flop，n为矩阵维数。本发明提出一种效率更高的列向量模值计算的方法，如下：

从表达式中易知 $\left\{\begin{array}{c}{a_i}^{(k+1)}={a_i}^{\left(k\right)}-{{\mathrm{\&alpha;a}}_j}^{\left(k\right)}\\ {a_j}^{(k+1)}={{\mathrm{\&alpha;a}}_i}^{\left(k\right)}+{a_j}^{\left(k\right)}\end{array}\right.,$ 所以

$\left\{\begin{array}{c}{I}^{(k+1)}=<{a_i}^{(k+1)},{a_i}^{(k+1)}>={({a_i}^{\left(k\right)}-{{\mathrm{\&alpha;a}}_j}^{\left(k\right)})}^T({a_i}^{\left(k\right)}-{{\mathrm{\&alpha;a}}_j}^{\left(k\right)})=I^{\left(k\right)}-{2\mathrm{\&alpha;e}}^{\left(k\right)}+{\mathrm{\&alpha;}}^2{J}^{\left(k\right)}\\ J^{(k+1)}=<{a_j}^{(k+1)},{a_j}^{(k+1)}>={({{\mathrm{\&alpha;a}}_i}^{\left(k\right)}+{a_j}^{\left(k\right)})}^T({{\mathrm{\&alpha;a}}_i}^{\left(k\right)}+{a_j}^{\left(k\right)})={\mathrm{\&alpha;}}^2{I}^{\left(k\right)}+2\mathrm{\&alpha;}{e}^{\left(k\right)}+J^{\left(k\right)}\end{array}\right.$

其中，上标(k)表示迭代k次后的结果。由此，可以在每次旋转变换结束前对列的模值按照上述表达式进行更新，这样只需要12个flop，仅为原来的

n为矩阵维数。在n较大的情况下，效率可以得到比较大的提升。

5)反复迭代：重复进行上述两步，选择Robin Ring Ordering作为列的正交化顺序，则A_k的各列趋向于两两正交，即A_k趋向于Q，而V＝R₀R₁R₂…R_k，即有AV＝Q，即A＝QV^-1，其中Q表示各列两两正交的矩阵，V表示各次旋转矩阵的乘积。

6)矩阵求逆：将Q矩阵各列化为单位向量，得Q＝Q₀Λ，Q₀为单位正交阵，Λ为对角阵，所以A＝Q₀ΛV^-1，则A^-1＝VΛ^-1Q₀ ^-1＝VΛ^-1Q₀ ^T，上标T表示矩阵的转置，此即完成了矩阵求逆。

本发明在单边Jacobi方法的基础上提出用快速Givens旋转矩阵代替常规的Givens旋转矩阵。采用快速Givens旋转矩阵，相较于常规的Givens旋转矩阵将节省一半的乘法运算量。采用快速Givens旋转矩阵虽然使矩阵失去了正交性，但是并不影响利用迭代的结果求取原矩阵的逆。针对单边Jacobi方法中每次迭代需要计算列向量的模，本发明提出利用迭代更新的方法，在每次矩阵修正后，及时更新改变了的列向量的模值，这样比直接计算列向量的模值效率更高。由于在求取每次的修正矩阵，即快速Givens矩阵时涉及复杂的开方运算，为了提高算法效率，减少开方运算，本发明提出采用极限逼近、泰勒展开式近似的方法来减少算法的开方运算。

有益效果：本发明提供了一种用于多天线系统的并行矩阵求逆方法，该方法基于单边Jacobi。相较于常规的单边Jacobi算法，本发明主要在三个方面提高了算法效率：第一，采用快速Givens旋转矩阵作为修正矩阵，节省了一半的乘法运算量；第二，采用极限和泰勒展开式近似，减少复杂的开方运算；第三，在计算列向量的模值时采用迭代的方法，对列向量模值进行更新。从这三个方面很大的提高了算法效率，以用于无线通信MIMO系统中的检测以及自适应传输技术中。

附图说明

图1是该发明的具体算法流程图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步说明。

如图1所示，一种用于多天线系统的并行矩阵求逆方法，包括以下步骤：

1)一个发送天线数为M、接收天线数为N的多天线系统，其信号的输入输出关系可表示为r＝HWs+n，其中：r表示接收信号，是维度为M的列向量；s表示发送信号，是维度为L的列向量；W表示发送预编码矩阵，是维度为N×L的矩阵；H表示信道矩阵，是维度为M×N的矩阵；n表示加性高斯白噪声，是维度为M的列向量；接收端进行检测，从接收到信号r中估计出发送信号s；多天线系统主流的检测算法包括迫零(ZF)算法和最小均方误差(MMSE)算法，其检测器表达式可分别表示为 ${\hat{s}}_{\mathrm{ZF}}={\left(H_e^H{H}_e\right)}^{-1}{H}_e^{H}r,$ ${\hat{s}}_{\mathrm{MMSF}}={(H_e^H{H}_e+{\mathrm{\&sigma;}}_n^2{I}_M)}^{-1}{H}_e^{H}r,$ 其中，

分别表示检测器采用迫零算法、最小均方误差算法对发送信号的估计，上标H表示矩阵的共轭转置，I_M表示M维单位矩阵，

表示噪声的平均功率，H_e＝HW表示等效信道矩阵；

2)为了使多天线系统的传输信号能自适应于信道环境的变化，需要在发送端根据信道条件，在给定准则下自适应地选取发送预编码矩阵W，最大化MMSE和速率是一种常用的选取预编码矩阵的准则，其计算公式为

$R_{\mathrm{MMSE}}\left(W\right)=\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{\log }_2\frac{1}{{\left[{({\mathrm{\&sigma;}}_n^{-2}{H}_e^H{H}_e+I_M)}^{-1}\right]}_{j,j}}$

其中R_MMSE(W)表示预编码矩阵为W时接收端接收信号产生的均方误差，下标j，j表示矩阵的对角线元素。

3)计算旋转矩阵：

将上述步骤1)和步骤2)中多天线系统中待求逆的矩阵

或矩阵

记为A，设初始矩阵A₀＝A，逐次选择修正矩阵R_k，使A_k的i，j两列正交化。修正矩阵R_k采用快速Givens矩阵的形式：

R_k，A_k分别表示第k次迭代时的旋转矩阵和待修正的矩阵。参数α计算如下：设I＝<a_i ^(k)，a_i ^(k)>，J＝<a_j ^(k)，a_j ^(k)>，e＝<a_i ^(k)，a_j ^(k)>，其中a_j ^(k)表示矩阵A_k的第j列的列向量。这里<x，y>表示向量x和y的内积，若e＝0，则说明两列已经正交，取α＝0，否则记

其中sgn(t)表示符号函数。α这样取值可以保证相应的旋转角满足

以提高算法的数值稳定性。注意到参数

这里涉及到开方运算，复杂度较高，考虑进行一定的近似，从而降低复杂度。注意到α关于t是奇对称的，所以仅需讨论t＞0的情形，即利用极限的知识以及Taylor(泰勒)展开式，易得：

$\underset{t\&RightArrow;\&infin;}{\lim }\sqrt{t^2+1}-t=\underset{t\&RightArrow;\&infin;}{\lim }\frac{1}{\sqrt{t^2+1}+t}=\underset{t\&RightArrow;\&infin;}{\lim }\frac{1}{2t},$

$\underset{t\&RightArrow;0}{\lim }\sqrt{t^2+1}-t\&ap;\underset{t\&RightArrow;0}{\lim }(1+\frac{t^2}{2})-t$

所以，在实际计算时，可作如下近似：

$\widetilde{\mathrm{\&alpha;}}=\left\{\begin{array}{cc}1-t+\frac{t^2}{2}& 0<t\&le;0.25\\ -t+\sqrt{t^2+1}& 0.25<t\&le;8\\ \frac{1}{2t}& 8<t\&le;1024\\ 0& 1024<t\end{array}\right.$

其中表示对α的近似。分段点的选取是随用户要求的精度变化而变化的，上述分段点的选取是假设控制误差在10^-4数量级内得到的。

4)进行旋转变换：A_k+1＝A_kR_k(k＝0，1，…)：

注意到在每次旋转变换时都需要计算第i，j列的模值，即I，J的值，常规的单侧Jacobi算法采用直接求取列向量模值的方法，共需要4n个flop，n为矩阵维数。本发明提出一种效率更高的列向量模值计算的方法，如下：

从表达式中易知 $\left\{\begin{array}{c}{a_i}^{(k+1)}={a_i}^{\left(k\right)}-{{\mathrm{\&alpha;a}}_j}^{\left(k\right)}\\ {a_j}^{(k+1)}={{\mathrm{\&alpha;a}}_i}^{\left(k\right)}+{a_j}^{\left(k\right)}\end{array}\right.,$ 所以

$\left\{\begin{array}{c}{I}^{(k+1)}=<{a_i}^{(k+1)},{a_i}^{(k+1)}>={({a_i}^{\left(k\right)}-{{\mathrm{\&alpha;a}}_j}^{\left(k\right)})}^T({a_i}^{\left(k\right)}-{{\mathrm{\&alpha;a}}_j}^{\left(k\right)})=I^{\left(k\right)}-{2\mathrm{\&alpha;e}}^{\left(k\right)}+{\mathrm{\&alpha;}}^2{J}^{\left(k\right)}\\ J^{(k+1)}=<{a_j}^{(k+1)},{a_j}^{(k+1)}>={({{\mathrm{\&alpha;a}}_i}^{\left(k\right)}+{a_j}^{\left(k\right)})}^T({{\mathrm{\&alpha;a}}_i}^{\left(k\right)}+{a_j}^{\left(k\right)})={\mathrm{\&alpha;}}^2{I}^{\left(k\right)}+2\mathrm{\&alpha;}{e}^{\left(k\right)}+J^{\left(k\right)}\end{array}\right.$

其中，上标(k)表示迭代k次后的结果。由此，可以在每次旋转变换结束前对列的模值按照上述表达式进行更新，这样只需要12个flop，仅为原来的

在n较大的情况下，效率可以得到比较大的提升。

5)反复迭代：重复进行上述两步，选择Robin Ring Ordering作为列的正交化顺序，则A_k的各列趋向于两两正交，即A_k趋向于Q，而V＝R₀R₁R₂…R_k，即有AV＝Q，即A＝QV^-1，其中Q表示各列两两正交的矩阵，V表示各次旋转矩阵的乘积。

6)矩阵求逆：将Q矩阵各列化为单位向量，得Q＝Q₀Λ，Q₀为单位正交阵，Λ为对角阵，所以A＝Q₀ΛV^-1，则A^-1＝VΛ^-1Q₀ ^-1＝VΛ^-1Q₀ ^T，上标T表示矩阵的转置，此即完成了矩阵求逆。

应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (3)

1.一种用于多天线系统的并行矩阵求逆方法，其特征在于：包括以下步骤：

1)一个发送天线数为M、接收天线数为N的多天线系统，其信号的输入输出关系表示为r＝HWs+n，其中：r表示接收信号，是维度为M的列向量；s表示发送信号，是维度为L的列向量；W表示发送预编码矩阵，是维度为N×L的矩阵；H表示信道矩阵，是维度为M×N的矩阵；n表示加性高斯白噪声，是维度为M的列向量；接收端进行检测，从接收到信号r中估计出发送信号s；多天线系统主流的检测算法包括迫零算法和最小均方误差算法，迫零算法检测器表达式表示为

其中，

表示检测器采用迫零算法对发送信号的估计，最小均方误差算法检测器表达式表示为

其中，

表示检测器采用最小均方误差算法对发送信号的估计，上标H表示矩阵的共轭转置，I_M表示M维单位矩阵，

表示噪声的平均功率，H_e＝HW表示等效信道矩阵；

2)为了使多天线系统的传输信号能自适应于信道环境的变化，需要在发送端根据信道条件，在给定准则下自适应地选取发送预编码矩阵W，最大化MMSE和速率是一种常用的选取预编码矩阵的准则，其计算公式为

$R_{\mathrm{MMSE}}\left(W\right)=\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{\log }_2\frac{1}{{\left[{({\mathrm{\&sigma;}}_n^{-2}{H}_e^H{H}_e+I_M)}^{-1}\right]}_{j,j}}$

其中R_MMSE(W)表示预编码矩阵为W时接收端接收信号产生的均方误差，下标j，j表示矩阵的对角线元素；

3)计算旋转矩阵：

将上述步骤1)和步骤2)中多天线系统中待求逆的矩阵

或矩阵

记为A，同时初始矩阵记为A₀＝A，选择旋转矩阵R_k，使A_k的i，j两列正交化，旋转矩阵R_k采用快速Givens矩阵的形式，其中R_k，A_k分别表示第k次迭代时的旋转矩阵和待修正的矩阵；

4)进行旋转变换：将求取的旋转矩阵作用于矩阵进行旋转变换：A_k+1＝A_kR_k(k＝0，1，…)；

5)反复迭代：选择Robin Ring Ordering作为列的正交化顺序，则A_k的各列趋向于两两正交，即A_k趋向于Q，而V＝R₀R₁R₂…R_k，即有AV＝Q，即A＝QV^-1，其中Q表示各列两两正交的矩阵，V表示每次旋转矩阵的乘积；

6)矩阵求逆：将Q矩阵各列化为单位向量，得Q＝Q₀Λ，Q₀为单位正交阵，Λ为对角阵，所以A＝Q₀ΛV^-1，则A^-1＝VΛ^-1Q₀ ^-1＝VΛ^-1Q₀ ^T，上标T表示矩阵的转置，此即完成了矩阵求逆。

2.根据权利要求1所述的一种用于多天线系统的并行矩阵求逆方法，其特征在于：所述步骤3)中旋转矩阵采用快速Givens矩阵的形式、参数的求取采用极限和泰勒展开式近似，来避免过多的开方运算。

3.根据权利要求1所述的一种用于多天线系统的并行矩阵求逆方法，其特征在于：所述步骤4)中旋转变换中在计算列向量的模值时采用迭代更新的方法来减少运算量。

Images