数值分析特征值Q.pptVIP

3.3 QR方法,在特征值计算问题上，QR方法具有里程碑意义。在1955年的时候，人们还觉得特征值的计算是十分困扰的问题，到1965年它的计算基于QR方法的程序已经完全成熟。直到今天QR方法仍然是特征值计算的有效方法之一。,QR 方法是 1961 年由 J.G.F.Francis (弗兰西斯)和V.N.Kublanovskaya 设计的。QR 分解是 QR 算法的基础。,QR 方法是一种变换方法，是计算一般矩阵（中小型矩阵）全部特征值问题的最有效方法之一。QR方法具有收敛快，算法稳定等特点.,计算特征值问题的 QR 方法，实际上总是分成 2 个阶段：,Householder方法（正交变换）,QR方法,QR方法,Householder方法（正交变换）,提问：对于一般实矩阵ARnn ,可以通过正交相似变换约化到什么程度？,补定理：设ARnn ，则存在一个正交矩阵R，使,其中对角块为一阶或二阶方阵，每一个一阶对角块即为A 的实特征值，每一个二阶对角块的两个特征值是A 的一对共轭复特征值。,定义: 拟上三角 (上Hessenberg)阵指一个 n 阶方阵 B， B 满足: 当 i j +1 时，bij= 0。 拟上三角矩阵有如下的形式:,本节讨论两个问题: 1、用正交相似变换约化一般实矩阵为上 Hessenberg 阵。 2、用正交相似变换约化对称矩阵为三对角阵。 这样，求原矩阵特征值问题，就转化为求上 Hessenberg 阵或对称三对角阵矩阵的特征值问题。,3.3.1 矩阵的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为两个不相等的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同向的向量。,解：,定理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, 有,故,若 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方法将下述矩阵作QR分解。, A=-4 -3 7;2 3 2;4 2 7 A = -4 -3 7 2 3 2 4 2 7 Q,R=QRself(A) num = 1 Hr = -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.5333 Sr = 0 1.5333 -0.9333 num = 2 Hr = 1.0000 0 0 0 -0.8542 0.5199 0 0.5199 0.8542,Ar = 6.0000 4.3333 0.6667 0.0000 -1.7951 2.1664 0.0000 0 9.8419 Q = -0.6667 0.0619 0.7428 0.3333 -0.8666 0.3714 0.6667 0.4952 0.5571 R = 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) 继续.,计算停止时, 正交阵 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 相似, 它们有相同的特征值. 最后, 矩阵序列 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= diag(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 迭代步的乘法运算次数只需 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) 全是零, 令 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)的形式如下:,特例：当 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 7 A = -4 -3 -7 2 3 2 4 2 7 nssj(A) num = 1 ,Hr = 1.0000 0 0 0 -0.4472 -0.8944 0 -0.8944 0.4472 Ar = -4.0000 7.6026 -0.4472 -4.4721 7.8000 -0.4000 0.0000 -0.4000 2.2000 ans = 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 都是拟上三角阵，此事实，由式(3.23)中的迭代公式可知。,书P66 第10题, A=1 2 1 2;2 1 2 -1;1 2 0 3 ;2 -1 3 1 A = 1 2 1 2 2 1 2 -1 1 2 0 3 2 -1 3 1 nssj(A) num = 1 Hr = 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.1556 1.2311 1.7289 num = 2 Hr = 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.8558 ans = 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 -1 2 1;-1 2 1 2;1 0 -1 3 A = 2 1 3 4 1 -1 2 1 -1 2 1 2 1 0 -1 3 nssj(A) num = 1 Hr = 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.8006 num = 2 Hr = 1.0000 0 0 0 0 1