第七章-特征值问题的迭代解法

第七章特征值问题的迭代解法当矩阵规模很大时,计算其所有的特征值和特征向量是非常困难的.而在实际应用中,我们通常也只对其中的某些特征值和特征向量感兴趣,因此也没必要计算所有的特征值和特征向量.本章讨论计算部分特征值和特征值向量的迭代解法.这些算法的存储量要远小于O(n2),运算量也远小于O(n3).
7-2
7.1投影算法最简单的特征值问题就是仅仅计算一个特征值,如计算模最大的特征值.这时我们可以使用幂迭代算法.算法7.1幂迭代:计算最大特征值1:Givenx(0)2:fori=1;2;:::,untilconvergedo3:y(i)=Ax(i1)4:x(i)=y(i)=∥y(i)∥25:endfor幂迭代所产生的迭代向量x(0),x(1);:::;x(m1)生成一个Krylov子空间Km(A;x(0))=span{x(0);Ax(0);:::;Am1x(0)}:在幂迭代中,我们取x(m1)为近似特征向量.显然,如果我们在Km(A;x(0))中找出“最佳”的近似特征向量,则收敛速度就可能会大大加快.下面我们讨论如何在Km=Km(A;x(0))中寻找“最佳”的近似特征向量.设A2Rnn,并设Km和Lm是Rn的两个m维子空间.投影算法就是在寻找A的近似特征对(~;~x),满足下面的Petrov-Galerkin条件find~2Cand~x2KmsuchthatA~x~~x?Lm:(7.1)
7-3
这样的算法我们称为斜投影算法.如果我们取Lm=Km,则上面的算法就是一个正交投影算法,此时条件(7.1)称为Galerkin条件.设v0;v1;:::;vm1和w0;w1;:::;wm1分别是Km和Lm的一组标准正交基.令Vm=[v0;v1;:::;vm1],Wm=[w0;w1;:::;wm1].则对任意~x2Km,存在向量y2Rn使得~x=Vmy,即~x可以由v0;v1;:::;vm1线性表出.根据条件(7.1),我们可得(AVmy~Vmy;wi)=0;i=1;2;:::;m;即WTmAVmy=~WTmVmy:(7.2)这是一个广义特征值问题.如果我们取Lm=Km,并令Wm=Vm,则(7.2)就化为Tmy=~y;(7.3)其中Tm=VTmAVm2Rmm.这意味着(~;y)是矩阵Tm的一个特征对.由于m通常比较小,因此我们可以使用先前讨论的方法(如QR迭代)来计算(~;y).这样我们就可以计算出A的一个近似特征对(~;~x),其中~x=Vmy.
7-4
7.2Rayleigh-Ritz算法事实上,我们可以在Km(A;x(0))中找出m个最佳近似特征向量及相应的最佳近似特征值.这些近似特征值和近似特征向量就是Ritz值和Ritz向量.定义7.1设Km是Rnn的一个m维子空间,它的一组标准正交基为v0;v1;:::;vm1,并令Vm=[v0;v1;:::;vm1].记Tm=VTmAVm,设(~;y)是Tm的一组特征对,即Tmy=~y且∥y∥2=1.则我们成~是A的一个Ritz值,~x=Vmy是A的一个Ritz向量.Rayleigh-Ritz算法就是用Ritz值和Ritz向量来近似A的特征值与特征向量.算法7.2RayleighRitzprocedure1:ComputeanorthonormalbasisofK:Vm=[v0;v1;:::;vm1].2:ComputetheeigenvaluesofthematrixTm=VTmAVm,i.e.,theRitzvaluesofA.3:Selectkdesiredones:~1;~2;:::;~kwherekm.4:ComputetheeigenvectorsyiofTmassociatedwith~i,i=1;2;:::;k.5:ApproximatetheeigenvectorsofAbytheRitzvectors~xi=Vmyi,i=1;2;:::;k.7.2.1对称矩阵这里我们讨论A对称时,RayleighRitz算法的性质.
7-5
设Vu2Rnnm是一个列正交矩阵,且使得V=[Vm;Vu]2Rnn是一个正交矩阵.于是我们有VTAV=[Vm;Vu]TA[Vm;Vu]=[VTmAVmVTmAVuVTuAVmVTuAVu]:由于A对称的,故Tm=VTmAVm2Rmm也是对称的.此时,Ritz值和Ritz向量具有下面的最优性质.定理7.1设A2Rnn对称,则对任意的对称矩阵R2Rmm,有∥AVmVmR∥2∥AVmVmTm∥2;即∥AVmVmR∥2在R=Tm处取最小值,此时∥AVmVmR∥2=∥VTuAVm∥2.证明.LetR=Tm+ZwhereZ2Rmmissymmetric.NotethatbothAandTmaresymmetricandTm=VTmAVm,wehave∥AVmVmR∥22=(AVmVmR)T(AVmVmR)=(AVmVm(Tm+Z))T(AVmVm(Tm+Z))=∥AVmVmTm∥22VTmAVmZ+TmVTmVmZZVTmAVm+ZVTmVmTm+ZTVTmVmZ=∥AVmVmTm∥22+∥Z∥22:Itfollowsthat∥AVmVmR∥22isminimizedwhenZ=0and∥AVmVmTm∥2=VVTAVmVmTm2
7-6
=(Vm;Vu)[VTmAVmVTuAVm]VmTm2=Vu(VTuAVm)2=VTuAVm2:2注:定理7.1中的2-范数可以改成任意的酉不变范数,如F-范数.定理7.2设A2Rnn对称,并设Tm=UUT是Tm=VTmAVm的特征值分解.设Q2Rnm是满足span(Q)=K的任意列正交矩阵,D2Rmm是任意对角矩阵.我们有∥AQQD∥2∥AVmVmTm∥2;且当Q=VmU,D=时等式成立.证明.Becausespan(Q)=K=span(Vm),itfollowsthatthereexistamatrixW2RmmsuchthatQ=VmW.AsQiscolumnorthogonal,wehaveI=QTQ=(VmW)TVmW=WTVTmVmW=WTW;whichimpliesthatWisorthogonal.LetWDWT=Tm+Z.enwehave∥AQQD∥22=∥AVmWVmWD∥22
7-7
=∥AVmVmWDWT∥22=∥AVmVmTm∥22+∥Z∥22∥AVmVmTm∥22:IfwechooseW=UandD=,thenZ=WDWTTm=UUTTm=0and,hence,theaboveequalityholds.2isresultshowsthattheminimumof∥AQQD∥2overalln-by-kcolumnorthogonalmatricesQwithspan(Q)=KandoveralldiagonalmatricesDisattainedbyQ=VmUandD=.定理7.1和定理7.2表明,在∥AQQD∥2极小的意义下,Ritz值是特征值的“最佳”近似.所以我们用Ritz值作为特征值的近似是有道理的.
7-8
7.3Lanczos算法设A2Rnn是对称矩阵.Lanczos就是利用Lanczos算法来计算Km的基和Tm=VTmAVm,然后计算A的Ritz值和Ritz向量.算法7.3LanczosAlgorithm1:Chooseavectorv0suchthat∥v0∥=1,andset0=02:forj=0;1;:::do3:Computew=Avjjvj14:j+1=(w;vj)5:w=wj+1vj6:j+1=∥w∥27:ifj+1=0then8:stop9:endif10:vj+1=w=j+111:ComputetheeigenvaluesandeigenvectorsofTj12:Checktheconvergence13:endfor
7-9
在Lanczos算法??中,迭代m步后,向量v0;v1;:::;vm1构成子空间Km(A;v0)=span{v0;Av0;:::;Am1v0};的一组基,并且有AVm=VmTm+mvmeTm;其中em=[0;0;:::;0;1]T2Rm,Tm=VTmAVm=2666664111............m1m1m3777775:设(~;y)是Tm的一个特征对,则有A(Vmy)=VmTmy+mvmeTmy=~Vmy+m(eTmy)vm:于是Ax~x=m(eTmy)vm;即∥Ax~x∥2=jm(eTmy)j;其中x=Vmy.如果jm(eTmy)j很小,则我们就认为~是A的某个特征值的很好的近似.事实上,关于Ritz值~,我们有下面的性质.
7-10
引理7.1设A2Rnn是对称矩阵,设r=Ax~x,其中x̸=0.则min2(A)j~j∥r∥2∥x∥2;其中(A)表示A的谱,即所有特征值组成的集合.证明.LetA=UUTbetheeigendecompositionofA.enwehaver=(A~I)x=U(~I)UTx:Denote=diag(1;2;:::;n).Foranyvectorz=[z1;z2;:::;zn]T2Rn,itholdsthat∥(~I)z∥22=n∑i=1(i~)2z2in∑i=1min2(A)j~j2z2i=∥z∥22min2(A)j~j2:erefore,wehave∥r∥2=∥UTr∥2=∥(~I)UTx∥2∥UTx∥2min2(A)j~j=∥x∥2
min2(A)j~j;
7-11
whichcompletetheproof.2由引理7.1可知,存在A的某个特征值,使得j~j∥m(eTmy)vm∥2∥Vmy∥2=jmj
jeTmyj∥y∥2=jmj
jeTmyj:(7.4)注:在前面的讨论中,我们都没有考虑实际计算时可能的误差.在实际计算中,由于浮点运算的误差,即使m很小(如m=10或m=20),也可能会导致向量fvig失去正交性.这时我们必须采取一些补救措施,最简单的方法就是对它们重新来一次正交话,即在算法??的第5步后加上一条语句w=wj∑i=1(w;vi)vi:这个过程就称为带全正交过程的Lanczos算法.显然,这个过程是非常费时的.另外一个可行的方法就是选择性正交.
7-12
7.4Arnoldi算法这里考虑非对称情形,即计算非对称矩阵A的特征值.与Lanczos算法的思想相类似,我们可以使用Arnoldi算法,即通过Arnoldi算法计算Km的标准正交基v0;v1;:::;vm1和上Hessenberg矩阵Hm=VTmAVm,使得AVm=VmHm+hm+1;mvmeTmandVTmAVm=Hm:但此时Hm只是上Hessenberg,而不是对称三对角.但我们同样可以通过计算Hm的特征值和特征向量来得到A的Ritz值的Ritz向量,并用它们来近似A的特征值和特征向量.设(~;y)是Hm的一个特征对,其中∥y∥2=1,则A(Vmy)=VmHmy+hm+1;mvmeTmy=~Vmy+hm+1;m(eTmy)vm;所以∥Ax~x∥2=∥hm+1;m(eTmy)vm∥2=jhm+1;mj
jeTmyj;其中x=Vmy.若jhm+1;mj
jeTmyj足够小,我就认为(~;x)是A的某个特征对的近似.算法7.4ArnoldiAlgorithm1:Chooseavectorv0suchthat∥v0∥2=12:forj=0;1;:::do3:Computewj+1=Avj
7-13
4:fori=0;1;:::;jdo5:hij=(wj+1;vi)6:wj+1=wj+1hijvi7:endfor8:hj+1;j=∥wj+1∥29:if