Chapter9-2-QR算法

QR方法在特征值计算问题的发展上具有里程碑意义。在1955年的时候人们还觉得特征值的计算是十分困扰的问题，到1965年它的计算——基于QR方法的程序已经完全成熟。直到今天QR方法仍然是特征值计算的有效方法之一。§7QR方法一Householder变换定义1设,且,则初等矩阵nRw12w称为Householder变换矩阵，也称初等镜面反射矩阵。Householder矩阵基本性质性质1H是对称正交矩阵，即性质2设HxyRxn,则性质2的意义是任意向量,在Householder变换作用下,Euclid长度不变.此外，由nRx知wxwyx)2(T由于xwT2是实数,上式表示向量x-y与向量w平行2)(TwwIwH1HHHT22xyxwwxHxyT2性质3设nRyx,yx22yx，则由向量2yxyxw确定的Household矩阵H(w),使得Hx=y。yyxxyxyxxxwwxHxTT22))((22xyxxyxxyyxyyxxxyxyxyxTTTTTTTT22)(222)()(例2设Tx)12,4,3(试求H矩阵,使TyHx)0,0,13(Tyxu)13,4,16(43-12-3-124-12-4-3-13112,4,16124164162100010001222uuuIHT直接验证THx)0,0,13(计算THxy)0,...,0,(的算法如下：21121))((signniixxTnTnuuuxxxu),...,,(),...,,(212111)(uxuxy二、矩阵的QR分解定理1设矩阵nnRA矩阵Q，上三角矩阵R，使A=QR且当要求R的主对角元素均为正数时，则分解式是唯一的。且非奇异，则一定存在正交A的非奇异性及Householder变换矩阵的性质知，一定可构造n-1个H矩阵)1,,2,1(1nkAHAkkkAA1)()(11)(22)(1)(121nnnnnnnnnnnnnaaσaσaaσAR例3设矩阵542112111---A试作矩阵A=QR分解。(1)求1H，作AHA12。;3))((sign1213121111iiaa;)2,2,4(,2,4221111Tuuau;12433111u21312222131111-------uuIHT433030033312----AHA(2)求2H，作RAHA223);1)0(sign(3))((sign122)2(2)2(2222约定iiaa;)3,3,0(,3,3,0223232)2(2221Tuauauu;93222u01010000131122uuIHT4R----AHA3003303332231222122213121HHQ由矩阵乘法可直接验证QRA。通常采用的方法是先将矩阵相似约化为上Hessenberg形式的矩阵,在此基础上应用QR迭代.这时,一个QR迭代步的乘法运算次数只需O(n２)次.三、相似约化为上Hessenberg矩阵对一般n阶矩阵,QR算法的每一个迭代步需要O(n３)次乘法运算.如果矩阵阶数稍大，这个算法几乎没有实际的应用价值.例如：一个5阶的上Hessenberg矩阵具有如下的形式：下面介绍QR方法时，都假设矩阵A是一个上Hessenberg矩阵.所谓上Hessenberg矩阵是指一个n阶矩阵A,如果当ij+1时，aij=0，称A为上Hessenberg矩阵.**000***00****0**********A定义1设A是n阶矩阵且有QR分解A＝QR，这里，Q是酉矩阵，R是上三角矩阵.如果A是满秩并规定R有正对角元，这个分解是惟一的.五、QR算法QR方法是1961年由作者J.G.F.Francis和V.N.Kublanovskaya设计的QR分解是QR算法的基础(1)QR算法的基本思想记A＝A1且有A1＝Q1R1.将等号右边两个矩阵因子的次序交换，得A2＝R1Q1且11112QAQA即12~AA不难证明:kkkkkkQQAQQQAQA1111111即11~~~AAAkk矩阵序列{Ak}有相同的特征值.因为上Hessenberg矩阵次对角线以下的元素全为0,因此,只要证明,当k→∞时,由迭代格式产生的矩阵Ak的次对角元趋向于零就可以了.记kkQQQ1~kkRRR1~容易得到是Ak的一个QR分解kkkRQA~~如果A是一个满秩的上Hessenberg矩阵,可以证明,经过一个QR迭代步得到的A2＝Q-11A1Q1仍然是上Hessenberg矩阵.这里,基本收敛的含义指{Ak}的元素中除对角线以下的元素趋于零外,可以不收敛于R的元素.(2)QR算法的收敛性设n阶矩阵A的n个特征值满足|λ1||λ2|…|λn|0，其相应的n个线性无关特征向量为x1,x2,…,xn.记X=(x1,x2,…,xn),Y=X-1.如果Y存在LU分解,那么,矩阵Ak基本收敛于上三角矩阵R.定理3(3)QR算法的迭代过程1.一个QR迭代步的计算①对l=1,2,…,n-1,构造n-1个平面旋转矩阵Pl,l+1,使A1的次对角元全部零化实现A1的QR分解的计算lljlllllllllllllaaaacssc,1,11,1,1,1,nllj,,2,12,121,llllllllaaac这里2,12,11,llllllllaaas②用Pl,l+1右乘，所得结果也放回矩阵A相应的元素中.1,1,1,1,1,1,lilillllllllliliaacsscaa1,,2,1li2.QR算法的迭代控制当迭代步数k充分大时,由迭代格式产生的Ak的次对角元趋于0.在实际计算中,控制迭代次数常用的一种办法是,预先给定一个小的正数ε,在一个迭代步的计算结束后,对l=n-1,n-2,…,1,依次判别次对角元的绝对值是否满足||||||,1Aall或更严格的准则|}||,min{|||1,1,1llllllaaa或不太严格的准则|}||{|||1,1,1llllllaaa如果上面三个不等式中有一个成立,把看做实际上为零.lla,1例4设矩阵410121012A试用QR算法求它的特征值。111286335.300449490.2449490.20447214.0236068.25912871.0408248.00365148.0816497.0447214.0182574.0408248.0894427.0RQA于是.000033416.103416.1.00003.095410.09541.00003112QRA同理作222QRA，又有.77271.973800.97380.52143.9558009558.0.70593223QRA如此下去，可得.26801.004800.00480.00873.1299001299.0.72334999QRA.26801.004800.00200.00353.0781000781.0.72824101010QRA从10A可以看出，已近似接近对角矩阵，即有特征值,2680.1,0035.3,7282.4321与矩阵A的三个精确解2679.133,3,7321.433321相比，已有良好精确度。随着迭代次数增加，nA将收敛到矩阵A的三个精确特征值。(4)带原点位移的ＱＲ算法由ＱＲ算法收敛性证明可以看出，ＱＲ算法的收敛速度依赖于矩阵相邻特征值的比值.为了加快算法的收敛速度,在迭代过程中,对矩阵Ak确定一个原点位移量sk,称Ak-skI为带原点位移量的矩阵.再对Ak-skI应用QR算法.这时,迭代格式改为kkkkRQIsAIsQRAkkkk1称为带原点位移的QR算法计算特征值问题的QR方法，实际上总是分成2个阶段：对称矩阵三对角矩阵对角矩阵一般矩阵上Hessenberg矩阵上三角矩阵