《数值分析特征值Q.ppt》由会员分享，可在线阅读，更多相关《数值分析特征值Q.ppt（59页珍藏版）》请在三一办公上搜索。
1、3.3 QR方法,在特征值计算问题上，QR方法具有里程碑意义。在1955年的时候，人们还觉得特征值的计算是十分困扰的问题，到1965年它的计算基于QR方法的程序已经完全成熟。直到今天QR方法仍然是特征值计算的有效方法之一。,QR 方法是 1961 年由(弗兰西斯)和 设计的。QR 分解是 QR 算法的基础。,QR 方法是一种变换方法，是计算一般矩阵（中小型矩阵）全部特征值问题的最有效方法之一。QR方法具有收敛快，算法稳定等特点.,计算特征值问题的 QR 方法，实际上总是分成 2 个阶段：,Householder方法（正交变换）,QR方法,QR方法,Householder方法（正交变换）,提问：
2、对于一般实矩阵ARnn,可以通过正交相似变换约化到什么程度？,补定理：设ARnn，则存在一个正交矩阵R，使,其中对角块为一阶或二阶方阵，每一个一阶对角块即为A 的实特征值，每一个二阶对角块的两个特征值是A 的一对共轭复特征值。,定义:拟上三角(上Hessenberg)阵指一个 n 阶方阵 B，B 满足:当 i j+1 时，bij=0。拟上三角矩阵有如下的形式:,本节讨论两个问题:1、用正交相似变换约化一般实矩阵为上 Hessenberg 阵。2、用正交相似变换约化对称矩阵为三对角阵。这样，求原矩阵特征值问题，就转化为求上 Hessenberg 阵或对称三对角阵矩阵的特征值问题。,3.3.1 矩
3、阵的QR分解理论,定义1 把矩阵 A 分解为一个正交阵 Q 与一个上三角阵 R 的乘积,称 A 的正交三角分解,简称QR分解.定义2 设 vRn是单位向量,即 vTv=1,令 H=I-2vvT(3.18)其中,I 是 nn 单位矩阵,则矩阵 H 满足:HT=H;HTH=(I-2vvT)(I-2vvT)=I.即,H 是对称的正交阵.式(3.18)定义的矩阵 H 称为Householder(豪斯荷尔德)矩阵,又称镜面映射矩阵.,注：镜面映射矩阵又称初等反射矩阵.简称H-矩阵，它被单位向量 v 唯一确定。,设向量u0，则显然H 是一个Householder(豪斯荷尔德)矩阵。,(补)定理 设x，y为
4、两个不相等的n维向量，则存在Householder矩阵H，使Hx=y,任给非零向量x，令y=Hx，则。正交变换是保模变换。,证明：令,引理3.1 设有非零向量 sRn 和单位向量 eRn,必存在 Householder 矩阵 H，使得Hs=e其中 是实数，并且。,证 取单位向量,其中,用单位向量 v 确定 一个 H-矩阵，则有,又因,故,代入 v 的表达式(3.19),例 给定向量x=(2,2,1)T，确定H，使Hx的后两个分量为零。,解：,防止a与x1抵消，以利于算法稳定。,书P66.9 设有向量s=(-2,1,2)T，r=(4,3,0)T，确定H，使u=Hs成为长度s与相等、方向与r同向的
5、向量。,解：,定理3.2 任何 n 阶方阵 A 总可以分解为一个 正交阵 Q 与一个上三角阵 R 的乘积。,证:正交阵 Q 与上三角阵 R 的构造过程如下：设 n 阶方阵 A=aij;er=(0,0,1,0,0)T 是 n 维单位坐标向量。（1）设 ai1(i=2,3,n)不全为零,令,A的第一列的长度,防止c1与a11抵消，以利于算法稳定。,构造 H-矩阵,由引理3.1,有,故,si 是 A 的列向量,若 ai1(i=2,3,n)全是零,令 H1=I,并且 A2=H1A=A.,（2）设 ai2(2)(i=3,4,n)不全为零,令,构造 H-矩阵,n-1阶H 矩阵,再用引理3.1,有,故,若
6、ai2(2)(i=3,4,n)全是零,令 H2=I,并且 A3=H2 A2=A2.,一般地,若依上法已经得到矩阵,又设不全为零,则令,构造 H-矩阵,有,及,若 air(r)(i=r+1,r+2,n)全是零,令 Hr=I,有 Ar+1=Hr Ar=Ar.,于是,当 r=n-1 时,得到 H-矩阵 H1,H2,Hn-1,使得An=Hn-1Hn-2H1 A.(3.20)是一个上三角阵.由上式可以得到A=(Hn-1Hn-2H1)-1 An=(H1Hn-2Hn-1)An 令 Q=H1Hn-2 Hn-1,R=An,Q 是正交阵,R 是上三角阵,且有等式 A=QR.,补例 用Householder方法将下
7、述矩阵作QR分解。,A=-4-3 7;2 3 2;4 2 7A=-4-3 7 2 3 2 4 2 7 Q,R=QRself(A)num=1Hr=-0.6667 0.3333 0.6667 0.3333 0.9333-0.1333 0.6667-0.1333 0.7333,Ar=6.0000 4.3333 0.6667 0.0000 1.5333 3.2667 0.0000-0.9333 9.5333Sr=0 1.5333-0.9333num=2Hr=1.0000 0 0 0-0.8542 0.5199 0 0.5199 0.8542,Ar=6.0000 4.3333 0.6667 0.0000
8、-1.7951 2.1664 0.0000 0 9.8419Q=-0.6667 0.0619 0.7428 0.3333-0.8666 0.3714 0.6667 0.4952 0.5571R=6.0000 4.3333 0.6667 0.0000-1.7951 2.1664 0.0000 0 9.8419,初始化:记 A1=A,并记 n 阶方阵 Ar=aij(r),令 Q1=I.对于 r=1,2,n-1 执行(1)若 air(r)(i=r+1,r+2,n)全是零,令 Qr+1=Qr,Ar+1=Ar.转(5),否则,转(2).(2)计算,QR 分解具体算法:,(3)令,(4)计算,(5)继续.
9、,计算停止时,正交阵 Q=Qn,上三角阵 R=An.,给定实矩阵 A,QR 法可求出A的全部特征值，过程如下：(1)令 A1=A,(2)做 A1的 QR 分解,A1=Q1R1,交换乘法次序,做 A2=R1Q1,则 A2=Q1TA1Q1,A2与A1相似.(3)对 A2再进行分解,A2=Q2R2,交换乘法次序,做 A3=R2Q2,则 A3与A2相似.,QR 方法的基本思想,Q1是正交阵,R1是上三角阵.,一般地,有(4)对 Ak 进行分解,Ak=QkRk,交换乘法次序,做 Ak+1=Rk Qk,则 Ak+1与Ak相似.Ak+1=QkTAkQk,所以,Ak+1与 A 相似,它们有相同的特征值.最后,
10、矩阵序列 Ak+1 基本上收敛于一个上三角阵 RA,RA 主对角线上的元素是其全部特征值,从而求出 A 的全部特征值.,对于给定的 n 阶实矩阵 A,(1)令 A1=A,(2)Ak=QkRk,(对 Ak作QR分解)(3)Ak+1=Rk Qk,(交换乘法次序)(k=1,2,n)Ak+1=QkTAkQk,故 Ak+1 与 A 相似.它们有相同的特征值.,基本 QR 法的迭代公式:,定理3.3 若 A 的特征值的重数均是 1，则由式(3.22)产生的矩阵序列 Ak+1 基本上收敛于一个分块上三角阵 RA,RA 的对角块均是 1 阶或 2 阶子块。特别,若A为实对称阵,则Ak+1收敛于对角阵 D=di
11、ag(1,n),i 是 A 的全部特征值.基本上收敛是指对角块上方的元素可能不收敛！,关于QR法的收敛性，有下面的定理,计算特征值问题的 QR 方法，实际上总是分成 2 个阶段：,为了减少计算量，一般先把矩阵 A 化为拟上三角阵 A（n-1），再用QR法计算 A（n-1）的全部特征值，而A（n-1）的特征值就是A的特征值。,3.3.2 矩阵的拟上三角化,对一般 n 阶矩阵，QR算法的每一个迭代步需要 O(n3)次乘法运算.如果矩阵阶数稍大，这个算法几乎没有实际的应用价值.通常采用的方法是先对矩阵作相似变换化为拟上三角矩阵(又称上Hessenberg 矩阵)在此基础上应用QR迭代.这时，QR 迭
12、代步的乘法运算次数只需 O(n2)次.,定义 拟上三角(Hessenberg)阵指一个 n 阶方阵 B，B 满足:当 i j+1 时，bij=0。拟上三角矩阵有如下的形式:,(1)设 ai1(i=3,4,n)不全为零,令,构造 H-矩阵,对于实矩阵 A 可做相似变换,把它化为拟上三角阵.变换矩阵用Householder 矩阵,过程如下:,由引理3.1,有,故,若 ai1(i=3,4,n)全是零,令 H1=I,并且 A(2)=H1AH1=A.,(2)设 ai2(2)(i=4,5,n)不全为零,令,构造 H-矩阵,n-2阶H-矩阵,再用引理3.1,有,故,若 ai2(2)(i=4,5,n)全是零,
13、令 H2=I,并且 A(3)=H2 A(2)H2=A(2).,(3)一般地,若依上法已经得到矩阵,又设 不全为零,则令,构造 H-矩阵,n-r阶H-矩阵,有,及,若 air(r)(i=r+2,r+3,n)全是零,令 Hr=I,有 A(r+1)=Hr A(r)Hr=A(r).,(4)当 r=n-2 时,得到拟上三角矩阵A(n-1)=Hn-2Hn-3H1 A H1Hn-3Hn-2.(3.21)注意到，Hr 是对称正交阵，故 Hr=HrT=Hr-1,令 P=H1Hn-3Hn-2,则有 PT=Hn-2Hn-3H1,A(n-1)=PTAP.故 A(n-1)与 A 相似。A(n-1)的形式如下:,特例：当
14、 A 是实对称阵时，A(n-1)也是实对称阵,因此它是对称三角阵：,初始化:记 A(1)=A,并记 A(r)的第r列到第n列元素为 aij(r),对于 r=1,2,n-2 执行(1)若 air(r)(i=r+2,r+3,n)全是零,令 A(r+1)=A(r),转(5),否则,转(2).(2)计算,拟上三角算法:,(3)令,(4)计算,(5)继续.,停止时,得到拟上三角阵 A(n-1).,例 用Householder方法将下述矩阵化为上Hessenberg阵。,上机计算的结果：A=-4-3-7;2 3 2;4 2 7A=-4-3-7 2 3 2 4 2 7 nssj(A)num=1,Hr=1.0
15、000 0 0 0-0.4472-0.8944 0-0.8944 0.4472Ar=-4.0000 7.6026-0.4472-4.4721 7.8000-0.4000 0.0000-0.4000 2.2000ans=1.0000 0 0 0-0.4472-0.8944 0-0.8944 0.4472,带双步位移的QR 法,3.3.3 双步位移的QR方法,其中 s1(k)，s2(k)是一对共轭复数或一对实数，称为位移量。s1(k)，s2(k)是 Ak 的右下角二阶子阵 Dk 的两个特征值。,又因为，,及,知，由式(3.23)产生的矩阵 Ak 均与A相似。此外，当 A 是拟上三角阵时，每个 Ak
16、 都是拟上三角阵，此事实，由式(3.23)中的迭代公式可知。,书P66第10题,A=1 2 1 2;2 1 2-1;1 2 0 3;2-1 3 1A=1 2 1 2 2 1 2-1 1 2 0 3 2-1 3 1 nssj(A)num=1Hr=1.0000 0 0 0 0-0.6667-0.3333-0.6667 0-0.3333 0.9333-0.1333 0-0.6667-0.1333 0.7333,Ar=1.0000-3.0000 0.0000 0.0000-3.0000 2.2222-2.7556 0.1556 0.0000-2.7556-1.9511 1.2311 0.0000 0.
17、1556 1.2311 1.7289num=2Hr=1.0000 0 0 0 0 1.0000 0 0 0 0-0.9984 0.0564 0 0 0.0564 0.9984,Ar=1.0000-3.0000-0.0000 0.0000-3.0000 2.2222 2.7599 0.0000-0.0000 2.7599-2.0780-1.0162 0.0000 0.0000-1.0162 1.8558ans=1.0000 0 0 0 0-0.6667 0.2952-0.6844 0-0.3333-0.9394-0.0805 0-0.6667 0.1745 0.7247,A=2 1 3 4;1-
18、1 2 1;-1 2 1 2;1 0-1 3A=2 1 3 4 1-1 2 1-1 2 1 2 1 0-1 3 nssj(A)num=1Hr=1.0000 0 0 0 0-0.5774 0.5774-0.5774 0 0.5774 0.7887 0.2113 0-0.5774 0.2113 0.7887,书P66第11题,Ar=2.0000-1.1547 3.7887 3.2113-1.7321-0.3333 0.7560-1.9107 0.0000-1.1874 2.5327 1.9880-0.0000-2.8541-0.3214 0.8006num=2Hr=1.0000 0 0 0 0 1
19、.0000 0 0 0 0-0.3841-0.9233 0 0-0.9233 0.3841,Ar=2.0000-1.1547-4.4203-2.2645-1.7321-0.3333 1.4737-1.4319 0.0000 3.0912 1.6473 0.0470-0.0000 0 2.3564 1.6860ans=1.0000 0 0 0 0-0.5774 0.3113-0.7548 0 0.5774-0.4981-0.6470 0-0.5774-0.8093 0.1078,书P66第12题,A=3 1 0;1 2 1;0 1 1A=3 1 0 1 2 1 0 1 1 QRself(A)nu
20、m=1Hr=-0.9487-0.3162 0-0.3162 0.9487 0 0 0 1.0000,Ar=-3.1623-1.5811-0.3162 0.0000 1.5811 0.9487 0 1.0000 1.0000Sr=0 1.5811 1.0000num=2,Hr=1.0000 0 0 0-0.8452-0.5345 0-0.5345 0.8452Ar=-3.1623-1.5811-0.3162-0.0000-1.8708-1.3363-0.0000 0 0.3381ans=-0.9487 0.2673 0.1690-0.3162-0.8018-0.5071 0-0.5345 0.8452,
