QR分解及其应用

1、矩阵分析与应用专题报告QR分解及应用学生姓名：卢楠、胡河群、朱浩2015年11月25日目录1 引言 32 QR 分解 42.1QR分解的性质 42.2 QR分解算法 52.2.1 采用修正Gram-Schmidt法的QR分解 52.2.2 Householder QR 分解 62.2.3 采用Give ns旋转的QR分解 83 QR分解在参数估计中的应用 931 基于 QR 分解的参数估计问题 93. 2 基于 Householder 变换的快速时变参数估计 123. 3基于Give ns旋转的时变参数估计 144 QR分解在通信系统中的应用 164.1基于QR分解的稳健干扰对齐算法 164.

2、2基于QR分解的MIMO置信传播检测器 19总结 21参考文献 221 引言矩阵分解是指将一个矩阵表示为结构简单或具有特殊性质的若干矩阵之积或 之和，大体上可以分为满秩分解、 QR 分解和奇异值分解。矩阵分解在矩阵分析 中占有很重要的地位，常用来解决各种复杂的问题。而QR分解是工程中应用最为广泛的一类矩阵分解。QF分解是目前求一般矩阵全部特征值的最有效并广泛 应用的方法，一般矩阵先经过正交相似变换成为 Hessenberg 矩阵，然后再应用 QF分解求特征值和特征向量。它是将矩阵分解成一个正交矩阵Q与上三角矩阵R, 所以称为QR分解。参数估计是在已知系统模型结构时， 用系统的输入与输出数据计算

3、系统模型参 数的过程。它在系统辨识和无线通信领域有着广泛的应用。 18 世纪末德国数学 家 C.F. 高斯首先提出参数估计的方法，他用最小二乘法计算天体运行的轨道。 20 世纪 60 年代，随着电子计算机的普及，参数估计有了迅猛的发展。参数估计 有很多方法，如矩估计、极大似然法、一致最小方差无偏估计、最小风险估计、 同变估计、最小二乘法、贝叶斯估计、极小极大熵法等。其中最基本的是最小二 乘法和极大似然法。本文将重点介绍QR分解及其在参数估计和通信系统中的应用2 QR 分解2.1QR分解的性质定理2.1.1 (QR分解)若a Rmn，且m n ,则存在列正交矩阵Q Rm n和上三角矩阵R Rm

4、n使得A = QR。当m n时，Q是正交矩阵。如果A是非 奇异的 n n 矩阵，则 R 的所有对角线元素均为正， 并且在这种情况下 Q 和 R 二者是唯一的。若 A 是复矩阵，则 Q 和 R 取复值。注意到 ATA =(QR)T(QR) =RTR , 因此可以得出结论： G = RT 是 ATA 的 下三角 Cholelskey 因子。由于这个原因，在关于估计的文献中，矩阵 R 常称为 平方根滤波器(算子) 。下面的引理称为矩阵分解引理，它在矩阵的QR分解的应用中是一个很有结果。 引理 2.2.1 若 A 和 B 是任意两个 m n 矩阵，则AHA = BHB(2.1.1)当且仅当存在一个 m

5、 m 酉矩阵 Q ，使得QA =B(2.1.2)证 明 充 分 性 证 明 ： 若 QA = B ， 并 且 Q 是 酉 矩 阵 ， 则 BHB = AHQHQA = AHA 。必要性证明：令A和B的奇异值分解分别为A = UA A VAHB = UB B VBH式中， UA和UB均为m m酉矩阵；Va和Vb都是n n酉矩阵；而m n矩阵a和b分别包含了矩阵A和B的非负奇异值。由于H H HA A = VA A A V AHBHB = VBHBQ = UbUaHhHh易知 QA = UbUa A = UbUa Ua a Va = Ub bVb = B这就证明了引理的必要条件10。2.2 QR分

6、解算法2.2.1 采用修正Gram-Schmidt法的QR分解 矩阵A 的QR分解可以利用Gram-Schmidt正交化方法实现。Gram-Schmidt正交1的向量q1 ,q 2 ,.,q n的方法。将向量a1标准正交化的结果取作q1，即R1q1a1q仁Rn(2.2.1)然后，从a2中除去与a1平行的向量,再进行标准正交化，并将结果取作q2，化方法原本是一种由 n个向量a1,a2,.,an构造互相正交且范数为则有进而，又从a 3除去与a1Hq1 a2R22I a2q1 R12(a2 q1R2)/R22R12q2和a?(2.2.2)平行的两个分量，再进行标准正交化，并使用该结果作q3，即有R|

7、3Hq1 a3R23Hq2 a3R33a3q 1 R13q2R23q3(a3 q1 R3q 2 尺3 ) R33(2.2.3)如此继续，则对于qk(2 k n)有(224)RjkqjHak,1 j k 1Rkkakqk (akAm n a1a2,., an容易验证，qi是标准正交基，即满足Hqi qj u(225)其中，ij为Kronecker 函数。如果令m n矩阵A的列向量ai，a2,.，an ,则 以qi,q2,.,qn为列向量的矩阵Q与A之间有下列关系：A =QR(2.2.6)又由于qi组成标准正交基，所以QHQ = In将A与Q重写在同一矩阵，应用以上Gram-Schmidt正交化的

8、方法叫做经典Gram-Schmidt正交化法。2.2.2 Householder QR 分解Householder变换可以实现任意m n矩阵A的QR分解，其原理是使用变维 向量的Householder变换，使得该向量除第一个元素外，其他元素皆为0。根据Householder变换的相关知识，欲使一个p维向量xx1,x2,.,x 的第1个元素后面的所有元素变为 0，则p维的Householder向量应取xe1X1)(2.2.7)X1X1式中(2.2.8)假定m n矩阵A的列分块形式为可以计算得到U1Wm。此时，H1 Iu1u1TA1 H1A a1(1),a2(1)a (1) n(2.2.9)变换后

9、，矩阵A1的第1列a1的第一个元素等于a1212 2a21. am11 2，而该列的其他元素全部为0。第二步针对矩阵A1的第2列aj，令P m 1和xa2；),a3；),.a(1) T,am2又可按照式(227)和式(2.2.8)求出(m-1)维向量Wm 1。此时，取U20Wm 1又可得到h2 IT典u?u 2A 2H 2A1 H2H1A a1(1), (2)a 2, ., a n(2.2.10)首先令x aiai1,a21,.,am1T，并取P m，则按照式(227)和式(228),变换后，矩阵A2的第1列与Ai的第1列相同，而第2列ai的第一个元素等于aj,第二个元素等于a22)am21

10、2，而该列的其他元素全部为0。类似地，又可针对矩阵A2的第3列设计Householder变换矩阵出，使得A 2的第一、二个元素保持不变，其他元素组成的m-2 维向量 xa3? , a43，,换为除第一个元素外的全部元素变为0假定矩阵A经过k-1次Householder变换后，已变成A(k 1)，即A(k 1) Hk1A(k2)Hk 1.H1A(k 1) (k 1)(k 1)a 1, a2 ,，a n并且其前k-1列具有以下变换结果：ajk1)a；：1,ajo,.,。因此，第k次Householder变换的目的就是保持前k列的下述变换：k=2,3，j=1,2,k-1k-1列不变，实现A(k1)列

11、第akkk0M0(k 1) ak, k(k 1)H% % 1,kk Ma(k 1)m,k这相当于对矩阵 A(k 1) 进行 Householder 变换 H kA(k 1) 时取I k 1 00H%kn次Householder变换后，即可实现 QR分解。223采用Give ns旋转的QR分解Givens 旋转也可以用来计算分解的思想：QR分解。这里以43 矩阵为例，说明 Givens QR(3,4)(2,3)(1,2) 0(3,4)000000(2,3)0(3,4)0000000000000其中，表示用 Givens旋转进行变化你的元素。从上述说明中易得出结论： 如果令 G j 代表约化过程中

12、的第 j 次 Givens 旋转，则QtA = R是上三角矩阵，其中Q = GtGtLG i，而t是总的旋转次数。3 QR分解在参数估计中的应用3. 1基于QR分解的参数估计问题现在以系统辨识为例，说明如何利用矩阵的 QR分解进行系统参数的递推估计令系统在k时刻的输入为x k，系统输出的观测值由卷积方程piX ki0i e kxT 6 e k(3.1.1)给出，其中，表示离散卷积,代表k时刻的观测误差，且Xk xk ,x1丄，k 1 ,L ,x kTP(3.1.2)若将k1,2,L ,n的所有观测数据组成一向量，则ynAn 6en(3.1.3)An x， yn y 1 , y 2 丄，y nT

13、1 ,x 2 ,L ,x n 0系统辨识问题的提法是：已知系统输入k 1,2,L ,n，估计系统参数向量670ene 1 ,e 2 丄，ek和输出观测值在时变系统的辨识中,n时刻的系统参数向量6的情况下，使用增加的x n 1 ,y，其中,则要求在已估计n 1值，通过简单的运算，递推出n 1时刻的系统参数向量n 1 0 n时刻的系统辨识问题可以 简化为最小二乘问题m6n An 6 仆:(3.1.4)求解，并且其解由“法方程”AnTAn 6 AnTyn 或 R xx 6n 心(3.1.5)确定。式中，Rxx AnTAn代表系统输入X k的协方差矩阵，rn人.丁丫直接求解式（3.1.5）的方法叫做协

14、方差方法。例如，先计算协方差矩阵Rxx的Cholesky分解Rxx GG,然后利用回带法解三角矩阵 G也G Rn直接 得到6n。然而，由于Rxx AnTAn的条件数是An的条件数的平方，因此，直接计算式（3.1.5）的得到的解有可能是严重病态的（即条件数很大），即使An本 身的条件数并不大，不是严重病态的。在系统参数向量9的自适应递推辨识中，标准的递推最小二乘RLS法和分解（其中,UtDU分解法都是针对协方差矩阵Rxx进行更新的。虽然utdu但是这些递推辨识方U为上三角矩阵，D为对角矩阵）在数值上比较稳定, 法也同样存在条件数变大的毛病。相比之下,An的QR分解可以保持原问题的条件数不变不妨令

15、式中，QnQnTAnRnO(3.1.6)是n n正交矩阵，R上三角矩阵，而0为p 1维零矩阵。由于正交变换可以保持被变换向量的Euclidean长度或范数不变，所以式（3.1.4）的最小二乘冋题可等价与作min Qn An 9QnTyn(3.1.7)(3.1.8)QnTyn且yn为p 11向量，为n P11向量，它们可以从QnTyn直接分块min Rn9 y式中(3.1.9)得到。一旦获得了入，即可由RnOnyn得到O。解此方程需要n n 1 /2次计算，并且最小残差值等于|%|假定增加了两个已知数值x n 1和y n 1，我们来讨论如何更新系统参数的估计，即使用已估计的参数向量3和简单的运算

16、，得到1时刻的新估计On 1。为了减少过去数据数据对参数估计的影响， 对数据X和y k 米取指数加权，即n 1时刻的数据矩阵和观测数据向量分别取作AnTXn 1,yn 1yny n 1(3.1.10)式中，01称为遗忘因子，且Xn 1 X n 1 ,x n丄，X n。于是,可以写出式(3.1.4)在n 1时刻的形式为Bn 1argmOn An 10 yn 1arg minAn yny n 1TXn 1乍一看，上式似乎没有什么特别吸引人之处，其实不然。这是因为,引理所述，式(3.1.11)的极小化变量等价为下述式的极小化变量,合于递推更新。引理 3.1.2r 若An Qn RO，QnTyn(3.

17、1.11)如同下面的而后者非常适y:，其中，Qn是正交矩阵，Rn是%n上三角矩阵，则式(3.1.11)的极小化变量等同于下式的极小化变量:Bn 1argmOn An 100yn 1arg minAnTXn 10 yn1y n 1(3.1.12)证明见。如果将式(3.1.12)的极小化变量记作0n 1，则以上讨论可总结为On 1的自适应递推估计算法如下。算法1 (系统参数的自适应估计算法)步骤1对矩阵RRnTXn 1进行QR分解，得qTF q：1 RnXn 1Rn 1O(3.1.13)式中，Qn 1是n1正交矩阵,Rn上三角矩阵,零矩阵。步骤二进行分块运算yn 1%1(3.1.14)其中，yn

18、1为p 11向量， 1为1向量。步骤三 切结三角矩阵方程Rn 1 9n 1 yn 1得到Bn 1。3.2基于Householder变换的快速时变参数估计考察n p 1矩阵a11a12La1,p 1Aa21a22La2,p 1A nMMMan,1an,2Lan, p 1的 Householder QR 分解，即*ana12La1,p 10*a22L*a2, p 1MMMHnAn00L*a p 1, p 1(3.2.1)00L0MMM00L0显然，只需要进行P次Householder变换即可。换言之，为了得到上述 QR分解，应该选择Hn为p个Householder变换矩阵之积，即Hn HnPHnP

19、lLHnl(322)式中Hn j I UjU； / j, j 1,2,L ,p(3.2.3):.T是对矩阵AnjHn j 1 Hn 2 Hn 1 An第j列向量印j禺,L曲 进行的Householder变换矩阵，其参数选择方法为2jaijjjajj,j 1,2L ,p(3.2.4)0,i jjUj i ajjsgn a”j,iaijj,i j其中并且Anj 1AnTUjqj(3.2.5)Tqj(3.2.6)递推的Householder QR分解算法如下:（1）递推更新基于QR分解的自适应参数估计算法一般由两个分开的过程组成:QR分解QtA R中的上三角矩阵R ；（2）用回代法求解三角矩阵方程。

20、由 于直接的回代需要O m2次运算（m为数据长度）。因此，即便Householder变 换再快速，整个自适应算法也至少需要O m2次运算。文献将上述快速Householder QR分解算法和求解三角矩阵方程的回代法综合起来考虑，提出了只具有0 m复杂度的快速自适应算法3.3基于Give ns旋转的时变参数估计现在考虑另外一种递推方法，递推求解0n的变化量,而不是直接递推求0n也1本身。换句话说，令(3.3.1)冋题是如何更新4 。假定正交矩阵Q为已知,它满足Q RnXn 1Rn 1O(332)由式(3.1.11)，式(3.3.1)和式(3.3.2)易知,4是下式的极小化变量:argmin QR

21、nTXn 1YnRn 0nXn 1(3.3.3)此式又可化简为式中，u n 1 y解出，其中，满足argmi nXn 1 0n。Rn 10因此,Rn 1 4Yn 1(3.3.4)可以从三角矩阵方程(3.3.5)(3.3.6)为了求出满足式(3.3.2)的Q ，可以使用Give ns平面旋转进行清零，将式(3.3.2)中的行向量X； 1的全部元素变成零。由于Q必须左乘式(3.3.6)，所以Yn 1r n 1RnTXn 15对增广的矩阵(3.3.7)执行所需要的清零。综合以上分析，在每一步递推更新中需要的步骤如下1）计算预测误差 yk 1（ 2）形成式 （3.3.7） 中的 n 1 n 1 矩阵；

22、（ 3）利用一系列 Givens 旋转将上述矩阵最底一行的左边 n 个元素扫除为零； （ 4）解上三角矩阵方程 （3.3.5） 得到 k 。利用Give ns旋转求解方程A 0 y的递推最小二乘算法的程序见文献。该算法中，同时对矩阵A和向量y应用Give ns旋转，因此无需存储正交矩阵 q。4 QR分解在通信系统中的应用4.1基于QF分解的稳健干扰对齐算法考虑K用户MIMO干扰信道，每个发送端的天线数为 Nt，每个接收端的天线 数为M ,每个用户对应的自由度为 dd2,L ,dk ，此处的自由度代表每个用 户能使用的独立数据流个数。为了让系统自由度达到最大值，即 Kmin(Nr,Nt)：2，那

23、么每个发送端所提供的信号空间的维数应该相等，故此处 不妨设di d2 L dK d ，并假设在同一时刻同一频率上的各个发送接收对之 间的信道是平坦衰落的，且信道系数独立同分布。在一个特定的时频资源上，接 收端i的接收信号可以表示为Kyi H ii Ws iHjiWjSj ni(4.1.1)j 1,j i其中维数为Nr Nt的Hii和Hji分别是发送端i和j到接收端i的信道矩阵。Wi和Wj分别是发送端i和j对应接收端i和j的预编码矩阵，且满足WiHWi Id，WjHWj Idj。维数为di 1的S是接收端i的下行数据矢量信 号，且满足功率约束E sHS P(i)。维数为Nr 1的ni是均值为0，

24、方差为1 的加性高斯白噪声噪声，且E ni 1比。干扰对齐往往要求完美的CSI，但在实际通信系统中，发送端得到CSI常常 是有误差的。为了构建稳健的干扰对齐算法，此处引入信道误差变量 EjiH ji H ji ，H ji表示真实的信道矩阵，Hji表示具有误差的信道矩阵，并且假设Er的元素服从均值为0,方差为2的循环对称复高斯分布(GSC)即H2满足 E EjiEjiel Nr。故式(4.1.1)变为_Kyi (Hii Eii)WSii(H ji Eji )WjSj ni(4.1.2)此时整个系统的联合接收信号可以表示为y1H11H21LHk1W1S1YyH12H 22LH k2V?S2IMMM

25、OMMykH1kH 2kLHkkWk Sk(4.1.3)E11E21LE k1W1s1E12E 22LEk2W2s2r2MMOMMME1kE2kLEkkWKSKnK对得到的误差联合信道矩阵H进行QR分解有HiiH21LHk1R11R21LRk1-Hi2H22LHk2R22LR k2H12QQR(4.1.4)MMOMOMH1kH2kLHkkRkk其中，Q是维数为KNrKNr的酉矩阵，R是维数为KNrKNt的上三角矩阵。因为Q是酉矩阵，根据矩阵理论可知 R和H有相同的统计特性，所以定义 R为 系统的误差等效联合信道矩阵。根据式(4.1.4)，式(4.1.3)可以改写为R11 R21LR k1W1s

26、1Y QR22LR k2W2s2OMMR kkWKSKE11E21LEk1W1s1n1E12E 22LEk2W2s2n?MMOMMME1kE2kLEkk1W KSKnK(4.1.5)这是考虑联合接收，对式(4.1.5)的联合信号接收信号进行左乘 Q 1的预处理, 得到如下的联合接收信号：R11R21LR k1W1s1E11E21LE k1W1s1YR22LR k2W2s2Q 1 E12E22LEk2W2s21OMMMMOMMRkkWKsKE1kE2kLEkkWKsKR11 R21 LR k1W S1rhR22 LR k 2W2s2Q(4.1.6)MOMMnKRkkWKsKE11E21LEk1W

27、1s1n1E12E22LEk2W2s2n2MMOMMME1kE2kLEkkWkSknK因为Q 1是酉矩阵,于是Ej和Ej有相同的统计特性，同理n j和n j有相同的统计特性。通过式（4.1.6），在接收端i经过干扰抑制矩阵Ui处理后，接收端i的接收信号为W1s1HUi 0,0, L ,Rii,R(i1)i ,L R KiW2s2MWKSKIIUiE1i,E2i,L EKiW1s1W2s2M(4.1.7)WKSK通过稳健干扰对齐算法得到最优干扰对齐矩阵uOpt和Wpt，具体算法流程总结为：（1）初始化Wi ，这里可以随机选择均值为0，方差为 1的矩阵Wi, i 1,2, L ,K ,WiHWi

28、Idi o（2）计算出Ui, i 1,2,L ,K，并且单位化Ui（3）计算出Wi, i 1,2,L ,K，并且单位化Wi（4）重复步骤（2）和（3），直到收敛。4.2基于QR分解的MIMO置信传播检测器在一个N个发射天线和M 个接收天线的MIMO系统中，s和K,snT是 N 1传输信号向量。系统的输入输出关系可以写为y Hs n(4.2.1)其中，y %,K ,ym 丁是M 1接收信号向量，n是M 1噪声向量，n的元素是零均值、方差2的独立同分布(i.i.d )复高斯随机变量。H是M N的MIMO信道矩阵，rank H N 。T由QR分解，MIMO信道矩阵H可以写为H Q RtrT ,其中Q

29、是M M酉矩阵，R。是MN N非零矩阵，R是NN上三角矩阵AG斤,2Lr1,Nr20RMMr2,2OLOr2,NM(4.2.2)50L0rN ,N由于Q Q I M M ,R。是非零矩阵，I NRtN 0QHnR。式(4.2.1)可以写为%I N N R 0QHyRs%(4.2.3)这样，R相比于由H得出的全连接的偶图就是一个包含更少的边数和环数的 的偶图。如下图所示。图中，N M 4 , Q 1。在第一次迭代前，全部0m初始化为0,其中1 n NQ且1 m g n 。在第l次迭代时，每一个 鳥可以由最大对数近似得到NQlmnmaxa AN m 1:bn 1l 1 kmf f m : bk1,

30、k nNQ(424)maxa An m 1:bn 0f f m:bk 1 ,k nl 1 km其中，bn和bk中的第n f m比特和第k f m 比特对应的s中的第n比特和第k比特,l 1km表示第l1次迭代时第k节点发出、第m节点检测的消息。计算出丨 后mn 丿口,lmng nlknk 1,k m(4.2.5)其中1n NQm g n 。第 l次迭代后bn的软输出L：是g n(4.2.6)lknk 1其中,1 n NQ。这样，给出的QR BP检测器就由上三角信道矩阵 R得lmn出的偶图来进行操作。上述QR BP的计算复杂度主要由式(4.2.3)中的线性变换和式(4.2.4)中的计算决定。线性

31、变换的复杂度开销主要来自于QR分解，它的复杂度是O MN 2 3。另外，mn的计算的大部分复杂度开销来自%ma的计算。M文献1给出了次迭代后该算法的计算复杂度是 O MN2m2mQ .所以m 1给出的QR BP检测器的计算复杂度小于文献1中标准BP检测器的复杂度OMN2nq 。总结通过此次文献调研，我们获益匪浅。首先，我们了解了矩阵分解和矩阵变换， 学习了 QR分解、Householder变换、Give ns旋转。其次，我们通过文献调研， 了解了 QR分解算法在实际中的应用，如参数估计，实际的MIMOS统等。最后,我们也从中学到了团队分工和协作的重要性。参考文献1 Hu, J., and Duma n, T.M.: Grap-based detect ion algorithms for layeredspacetime arc