第四讲非对称特征值问题设A是一个非对称的稠密矩阵，本讲主要讨论如何计算A的全部特征值和特征向量.为了讨论方便，本讲同样只考虑实矩阵情形记AERnxn的特征值为入1，入2...,入n.本讲中我们总是假定[] ≥[2] ≥ ≥[n/ ≥0,即A的特征值按模降序排列关于稠密矩阵特征值计算的相关参考资料J.H.Wilkinson，TheAlgebraicEigenvalueProblem,1965.[135](有中文翻译)B.N.Parlett, The Symmetric Eigenvalue Problem, 2nd Eds., 1998.[100]G. W. Stewart, Matrix Algorithms, Vol II:Eigensystems, 2001. [120]G. H. Golub and C. F. Van Loan, Matrix Computations, 2013. [57] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, Templates for the Solution of AlgebraicEigenvalue Problems: A Practical Guide, 2000. [7] G. H. Golub and H. A. van der Vorst, Eigenvalue computation in the 20th century, 2000. [56]4.1幂送代法幂选代法是计算特征值和特征向量的一种简单易用的算法，幂迭代法虽然简单，但它却建立了计算特征值和特征向量的算法的一个基本框架4.1.1算法介绍任意给定一个非零向量z(0),不断用A左乘该向量，得到向量序列2(0), A2(0), A2(0), A3(0),...在一定条件下，可以证明该向量序列收敛到模最大特征值入，所对应的特征向量，然后通过计算Rayleigh商（(Rayleigh商的定义见（5.10))就可以得到入1.这就是幕幂迭代法，算法描述如下，为了防止溢出，每次计算后对向量进行单位化处理算法4.1.幂送代法(PowerIteration)1: Choose an initial guess r(0) with Ir(0)|l2 = 12: set k = 03:while not convergencedoy(k+1) = Ar(k)4:129
第四讲 非对称特征值问题 设 A 是一个非对称的稠密矩阵, 本讲主要讨论如何计算 A 的全部特征值和特征向量. 为了讨 论方便, 本讲同样只考虑实矩阵情形. 记 A ∈ R n×n 的特征值为 λ1, λ2, . . . , λn. 本讲中我们总是假定 |λ1| ≥ |λ2| ≥ · · · ≥ |λn| ≥ 0, 即 A 的特征值按模降序排列. 关于稠密矩阵特征值计算的相关参考资料 ▶ J. H. Wilkinson, The Algebraic Eigenvalue Problem, 1965. [135] (有中文翻译) ▶ B. N. Parlett, The Symmetric Eigenvalue Problem, 2nd Eds., 1998. [100] ▶ G. W. Stewart, Matrix Algorithms, Vol II: Eigensystems, 2001. [120] ▶ G. H. Golub and C. F. Van Loan, Matrix Computations, 2013. [57] ▶ Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide, 2000. [7] ▷ G. H. Golub and H. A. van der Vorst, Eigenvalue computation in the 20th century, 2000. [56] 4.1 幂迭代法 幂迭代法是计算特征值和特征向量的一种简单易用的算法. 幂迭代法虽然简单, 但它却建立 了计算特征值和特征向量的算法的一个基本框架. 4.1.1 算法介绍 任意给定一个非零向量 x (0) , 不断用 A 左乘该向量, 得到向量序列: x (0), Ax(0), A2x (0), A3x (0) , . . . . 在一定条件下, 可以证明该向量序列收敛到模最大特征值 λ1 所对应的特征向量, 然后通过计算 Rayleigh 商 (Rayleigh 商的定义见 (5.10)) 就可以得到 λ1. 这就是幂迭代法, 算法描述如下, 为了防 止溢出, 每次计算后对向量进行单位化处理. 算法 4.1. 幂迭代法 (Power Iteration) 1: Choose an initial guess x (0) with ∥x (0)∥2 = 1 2: set k = 0 3: while not convergence do 4: y (k+1) = Ax(k) 129
-130.第四讲非对称特征值问题z(k+1)=y(k+1)/ly(+1)l2%单位化，防止溢出5:μk+1 = (Ar(+1), r(k+1)6:%计算Rayleigh商7:k=k+18: end while4.1.2收敛性分析下面讨论幂送代法的收敛性，我们假定(1)ARnxn是可对角化的，即A=VAV-1,其中A=diag(入1,入2,...,入n）ECnxn,V=[u1,v2,...,n] ECnxn,且 uil2 =1,i=1,2,..,n.(2) [1/ >[2/ ≥ [3/ ≥ - ≥ [n].由于VeCnxn非奇异,所以它的列向量组构成Cn的一组基.因此迭代初始向量r()可表示为2(0) =Qi1 + Q222 +...+QnUn = V[α1,2,.QnJT我们假定α1≠0,即(0)不属于spanu2,u3..,Un)(由于r(0)是随机选取的，从概率意义上讲,这个假设通常是成立的).于是我们可得[a1ia1%()02Q.2A*(0) =(VAV-1)*V=VA.+..:LanAhan (An)anan又[;/11.显然，幕送代的收敛快慢取决于[2/>1的大小，[2/入1|越小，收敛越快。结论(1)(K)收敛到±i；(2)收敛到入1；(3)>2/入i|越小,收敛越快通过上面的分析可知，幂迭代只能用于计算矩阵的模最大的特征值和其相应的特征向量如果A的模最大的特征值是唯一的，则称其为主特征值.当入2/入1接近于1时，收敛速度会非常慢同时，如果A的模最大特征值不唯一，比如一对共轭复特征值，则幂送代就可能就会失效
· 130 · 第四讲 非对称特征值问题 5: x (k+1) = y (k+1)/∥y (k+1)∥2 % 单位化, 防止溢出 6: µk+1 = Ax(k+1), x(k+1) % 计算 Rayleigh 商 7: k = k + 1 8: end while 4.1.2 收敛性分析 下面讨论幂迭代法的收敛性. 我们假定 (1) A ∈ R n×n 是可对角化的, 即 A = V ΛV −1 , 其中 Λ = diag(λ1, λ2, . . . , λn) ∈ C n×n , V = [v1, v2, . . . , vn] ∈ C n×n , 且 ∥vi∥2 = 1, i = 1, 2, . . . , n. (2) |λ1| > |λ2| ≥ |λ3| ≥ · · · ≥ |λn|. 由于 V ∈ C n×n 非奇异, 所以它的列向量组构成 C n 的一组基. 因此迭代初始向量 x (0) 可表示为 x (0) = α1v1 + α2v2 + · · · + αnvn = V [α1, α2, . . . , αn] T . 我们假定 α1 ̸= 0, 即 x (0) 不属于 span{v2, v3, . . . , vn} (由于 x (0) 是随机选取的, 从概率意义上讲, 这个假设通常是成立的). 于是我们可得 A kx (0) = (V ΛV −1 ) kV α1 α2 . . . αn = V Λ k α1 α2 . . . αn = V α1λ k 1 α2λ k 2 . . . αnλ k n = α1λ k 1V 1 α2 α1 λ2 λ1 k . . . αn α1 λn λ1 k . 又 |λi/λ1| < 1, i = 2, 3, . . . , n, 所以 lim k→∞ λi λ1 k = 0, i = 2, 3, . . . , n. 故当 k 趋向于无穷大时, 向量 " 1, α2 α1 λ2 λ1 k , . . . , αn α1 λn λ1 k #T , k = 0, 1, 2, . . . 收敛到 e1 = [1, 0, . . . , 0]T. 所以向量 x (k) = Akx (0)/∥Akx (0)∥2 收敛到 ±v1, 即 A 的对应于 (模) 最 大的特征值 λ1 的特征向量. 而 µk = Ax(k) , x(k) 则收敛到 v ∗ 1Av1 = λ1. 显然, 幂迭代的收敛快慢取决于 |λ2/λ1| 的大小, |λ2/λ1| 越小, 收敛越快. 结论 (1) x (k) 收敛到 ±v1; (2) µk 收敛到 λ1; (3) |λ2/λ1| 越小, 收敛越快. 通过上面的分析可知, 幂迭代只能用于计算矩阵的模最大的特征值和其相应的特征向量. 如 果 A 的模最大的特征值是唯一的, 则称其为主特征值. 当 |λ2/λ1| 接近于 1 时, 收敛速度会非常慢. 同时, 如果 A 的模最大特征值不唯一, 比如一对共轭复特征值, 则幂迭代就可能就会失效
4.1幂选代法+131 .思考：如果A不可以对角化，则能否得到类似的结论?如果需要计算其他特征值，比如模第二大特征值入2，则可以在模最大特征值入1计算出来后，采用收缩(Deflation）技术：构造西矩阵U，使得A1 A12U*AU0A22然后将幂选代作用到A22上，就可以求出入2.以此类推，可以依次求出所有特征值（这里假定特征值互不相同）思考：上面的收缩技术中的U怎么选取？4.1.3位移策略前面已经提到，幕选代法的收敛速度取决于42/入1|的大小，当它的值接近于1时，收敛速度会非常缓慢.因此，为了加快幂送代法的收敛速度，我们希望{入2/入1|的值越小越好一个简单易用的方法就是使用位移策略，即将计算A的特征值转化为计算A一αI的特征值，即对A做一个移位.这里α是一个给定的数，称α为位移(shift).为了使得幂选代作用到A-αI时具有更快的收敛速度，我们要求。满足下面的两个条件：(1）-是A-I的模最大的特征值[入, -0尽可能地小(2)max2<i<n1其中第一个条件保证最后所求得的特征值是我们所要的，第二个条件用于加快幂代的收敛速度显然，在实际应用中，的取值并不是一件容易的事，例4.1设A=XAX-1,其中A为对角矩阵，分别用幂迭代法和带位移的幕迭代法计算A的主特征值(见Eig_Power_shift.m)位移策略在特征值计算中非常重要，特别是在反选代法和QR选代法中
4.1 幂迭代法 · 131 · 思考：如果 A 不可以对角化, 则能否得到类似的结论? b 如果需要计算其他特征值, 比如模第二大特征值 λ2, 则可以在模最大特征值 λ1 计算出来 后, 采用 收缩 (Deflation) 技术: 构造酉矩阵 U, 使得 U ∗AU = " λ1 A12 0 A22# . 然后将幂迭代作用到 A22 上, 就可以求出 λ2. 以此类推, 可以依次求出所有特征值 (这里假 定特征值互不相同). 思考：上面的收缩技术中的 U 怎么选取? 4.1.3 位移策略 前面已经提到, 幂迭代法的收敛速度取决于 |λ2/λ1| 的大小. 当它的值接近于 1 时, 收敛速度 会非常缓慢. 因此, 为了加快幂迭代法的收敛速度, 我们希望 |λ2/λ1| 的值越小越好. 一个简单易用的方法就是使用位移策略, 即将计算 A 的特征值转化为计算 A − σI 的特征值, 即对 A 做一个移位. 这里 σ 是一个给定的数, 称 σ 为位移 (shift). 为了使得幂迭代作用到 A − σI 时具有更快的收敛速度, 我们要求 σ 满足下面的两个条件: (1) λ1 − σ 是 A − σI 的模最大的特征值; (2) max 2≤i≤n λi − σ λ1 − σ 尽可能地小. 其中第一个条件保证最后所求得的特征值是我们所要的, 第二个条件用于加快幂迭代的收敛速 度. 显然, 在实际应用中, σ 的取值并不是一件容易的事. 例 4.1 设 A = XΛX−1 , 其中 Λ 为对角矩阵, 分别用幂迭代法和带位移的幂迭代法计算 A 的主 特征值. (见 Eig_Power_shift.m) b 位移策略在特征值计算中非常重要, 特别是在反迭代法和 QR 迭代法中
-132.第四讲非对称特征值问题4.2反迭代法如果我们将幂选代法作用在A-1上，则可求出A的模最小的特征值，此时的幂选代法称为反代法.4.2.1算法介绍事实上，结合位移策略，我们就可以计算矩阵的任意一个特征值.下面给出带位移反迭代法的算法描述算法4.2.带位移的反选代法(InverseIteration)1: Choose a scalar α and an initial vector r(0) with [(0) l2 = 12: set k = 03: while not convergence dog(++1) =(A- 0I)-1,()4:2(k+1) = g(k+1)/ lg(k+1) /25:μk+1 = (A2(k+1),2(k+1)6:7:k=k+18: end while显然，在反迭代法中，收敛到距离α最近的特征值，而(）则收敛到其对应的特征向量设距离最近的特征值为入，则算法的收敛速度取决于入-0max1<i<n,i+kAi-g的大小，显然，越接近于入，其值越小，即算法收敛越快若。～入，则选代几步就可以了反选代法的另一个优点是，只要选取合适的位移，就可以计算A的任意一个特征值但反选代法的缺点也很明显：每步选代需要解一个线性方程组（A一gI)y(k+1)=2(k),这就需要对A－αI做一次LU分解.另外，与幕送代一样，反送代法一次只能求一个特征值4.2.2Rayleigh商送代在反送代法中，有一个很重要的问题，就是位移的选取显然，选取的基本原则是使得其与所求的特征值越靠近越好.这里需要指出的是，在每步迭代中，位移。可以不一样，即在选代过程中可以不断更新位移，由于在反迭代法4.2中，是收敛到所求的特征值的，所以我们可以选取作为第k步的位移这时所得到的反送代法就称为Rayleigh商迭代（RayleighQuotientIteration)，简记为RQI算法4.3.Rayleigh商选代法(RayleighQuotient Iteration (RQI)1: Choose an initial vector r(0) with /r(0)/2 = 1
· 132 · 第四讲 非对称特征值问题 4.2 反迭代法 如果我们将幂迭代法作用在 A−1 上, 则可求出 A 的模最小的特征值, 此时的幂迭代法称为反 迭代法. 4.2.1 算法介绍 事实上, 结合位移策略, 我们就可以计算矩阵的任意一个特征值. 下面给出带位移反迭代法的 算法描述. 算法 4.2. 带位移的反迭代法 (Inverse Iteration) 1: Choose a scalar σ and an initial vector x (0) with ∥x (0)∥2 = 1 2: set k = 0 3: while not convergence do 4: y (k+1) = (A − σI) −1x (k) 5: x (k+1) = y (k+1)/∥y (k+1)∥2 6: µk+1 = (Ax(k+1), x(k+1)) 7: k = k + 1 8: end while 显然, 在反迭代法中, µk 收敛到距离 σ 最近的特征值, 而 x (k) 则收敛到其对应的特征向量. 设距离 σ 最近的特征值为 λk, 则算法的收敛速度取决于 max 1≤i≤n,i̸=k λk − σ λi − σ 的大小. 显然, σ 越接近于 λk, 其值越小, 即算法收敛越快. 若 σ ≈ λk, 则迭代几步就可以了. 反迭代法的另一个优点是, 只要选取合适的位移 σ, 就可以计算 A 的任意一个特征值. 但反迭代法的缺点也很明显: 每步迭代需要解一个线性方程组 (A−σI)y (k+1) = x (k) , 这就需 要对 A − σI 做一次 LU 分解. 另外, 与幂迭代一样, 反迭代法一次只能求一个特征值. 4.2.2 Rayleigh 商迭代 在反迭代法中, 有一个很重要的问题, 就是位移 σ 的选取. 显然, 选取 σ 的基本原则是使得其与所求的特征值越靠近越好. 这里需要指出的是, 在每步迭 代中, 位移 σ 可以不一样, 即在迭代过程中可以不断更新位移 σ. 由于在反迭代法 4.2 中, µk 是收敛到所求的特征值的, 所以我们可以选取 µk 作为第 k 步的位 移. 这时所得到的反迭代法就称为 Rayleigh 商迭代 (Rayleigh Quotient Iteration), 简记为 RQI. 算法 4.3. Rayleigh 商迭代法 (Rayleigh Quotient Iteration (RQI)) 1: Choose an initial vector x (0) with ∥x (0)∥2 = 1
4.2反选代法+133.2: set k = 03: compute g = (r(0), Ar(0)4: while not converge doy(k+1) = (A -αI)-1(k)5:2(k+1) = g(+1)/lg(k+1) /26:μk+1 = (Ar(++1), 2(++1)7:8:g=μk+19:k=k+110: end while-般来说，如果Rayleigh商迭代收敛到A的一个单特征值，则至少是二次收敛的，即具有局部二次收敛性.如果A是对称的，则能达到局部三次收敛在Rayleigh商选代中，由于每次选代的位移是不同的，因此每次选代需要求解一个不同的线性方程组，这使得运算量大大增加.例4.2设A=XAX-1,其中A为对角矩阵，用Rayleigh商选代计算A的特征值(见Eig_Rayleigh.m)
4.2 反迭代法 · 133 · 2: set k = 0 3: compute σ = (x (0), Ax(0)) 4: while not converge do 5: y (k+1) = (A − σI) −1x (k) 6: x (k+1) = y (k+1)/∥y (k+1)∥2 7: µk+1 = (Ax(k+1), x(k+1)) 8: σ = µk+1 9: k = k + 1 10: end while 一般来说, 如果 Rayleigh 商迭代收敛到 A 的一个单特征值, 则至少是二次收敛的, 即具有局部 二次收敛性. 如果 A 是对称的, 则能达到局部三次收敛. 在 Rayleigh 商迭代中, 由于每次迭代的位移是不同的, 因此每次迭代需要求解一个不同的线 性方程组, 这使得运算量大大增加. 例 4.2 设 A = XΛX−1 , 其中 Λ 为对角矩阵, 用 Rayleigh 商迭代计算 A 的特征值. (见 Eig_Rayleigh.m)
-134.第四讲非对称特征值问题4.3正交送代法幂选代和反迭代都只能同时计算一个特征对.如果想同时计算多个特征对，我们可以采用多个初始向量进行选代.而正交选代算法就是基于这种思想，它能够计算A的一个不变子空间，从而可以同时计算出多个特征值算法4.4.正交选代法(Orthogonal Iteration)1: Choose an n × p column orthogonal matrix Zo2: set k = 03: while not convergence do4:computeYk+1=AZk5:Yk+1=Zk+1Rk+1%QR分解6:k=k+17: end while正交选代方法有时也称为子空间选代方法（subspaceiterationmethod）和同步选代方法（si-multaneous iteration method).在算法中使用QR分解是为了保持Zk的列正交性，使得其列向量组构成子空间spanAZo的一组正交基。一方面提高算法的数值稳定性，另一方面避免所有列都收敛到最大特征值所对应的特征向量.下面我们分析该算法的收敛性质.假设A是可对角化的，即A=VAV-1,其中A=diag(入1,入2,..,入n),且[^1|≥≥[入p|>[p+1/≥·.≥入nl.则可得span[Zk] = span[Yi] = span[AZk-1],k = 1,2, ...由此可知span[Zk} = span[A* Zo] = span[VA*-1 Zo].我们注意到[(A1 / ^,)k[w()]kV-1Z0=V-1 Zo>[w(b,](An/入p)k由于当i>p时有α/>pl[Vp, Vn-p]W(k)
· 134 · 第四讲 非对称特征值问题 4.3 正交迭代法 幂迭代和反迭代都只能同时计算一个特征对. 如果想同时计算多个特征对, 我们可以采用多 个初始向量进行迭代. 而正交迭代算法就是基于这种思想, 它能够计算 A 的一个不变子空间, 从 而可以同时计算出多个特征值. 算法 4.4. 正交迭代法 (Orthogonal Iteration) 1: Choose an n × p column orthogonal matrix Z0 2: set k = 0 3: while not convergence do 4: compute Yk+1 = AZk 5: Yk+1 = Zk+1Rˆ k+1 % QR 分解 6: k = k + 1 7: end while b 正交迭代方法有时也称为子空间迭代方法 (subspace iteration method) 和同步迭代方法 (simultaneous iteration method). 在算法中使用 QR 分解是为了保持 Zk 的列正交性, 使得其列向量组构成子空间 span{AkZ0} 的一组正交基. 一方面提高算法的数值稳定性, 另一方面避免所有列都收敛到最大特征值所对应 的特征向量. 下面我们分析该算法的收敛性质. 假设 A 是可对角化的, 即 A = V ΛV −1 , 其中 Λ = diag(λ1, λ2, . . . , λn), 且 |λ1| ≥ · · · ≥ |λp| > |λp+1| ≥ · · · ≥ |λn|. 则可得 span{Zk} = span{Yk} = span{AZk−1}, k = 1, 2, . . . , 由此可知 span{Zk} = span{A kZ0} = span{V Λ kV −1Z0}. 我们注意到 Λ kV −1Z0 = λ k p (λ1/λp) k . . . 1 . . . (λn/λp) k V −1Z0 ≜ λ k p " W (k) p W (k) n−p # . 由于当 i > p 时有 |λi/λp| < 1, 所以当 k 趋于无穷大时, W (k) n−p 趋向于 0. 令 V = [Vp, Vn−p], 则 V Λ kV −1Z0 = λ k p [Vp, Vn−p] " W (k) p W (k) n−p # = λ k p VpW(k) p + Vn−pW (k) n−p
4.3正交送代法?135.所以当k→8时,有span[Zk) = span(VA*V-1 Z0) = span[Vw() + Vn-pw(-)]→ span[V,W(*)] = span[Vp],即span[Zk}趋向于A的一个p维不变子空间span(Vp].定理41给定正整数p(1≤p≤n),考虑算法4.4.假设A是可对角化的,且[^i/ ≥.…≥/^pl>IΛp+1l≥··≥[^nl.则span[Zk}收敛到A的一个p维不变子空间当A不可对角化时，利用Jordan标准型，我们可以得到同样的结论，见[132，133].在正交选代中，如果我们取Zo=I，则可得到一类特殊的正交选代法.此时，在一定条件下正交迭代会收敛到A的Schur分解
4.3 正交迭代法 · 135 · 所以当 k → ∞ 时, 有 span{Zk} = span{V Λ kV −1Z0} = span{VpW(k) p + Vn−pW (k) n−p } → span{VpW(k) p } = span{Vp}, 即 span{Zk} 趋向于 A 的一个 p 维不变子空间 span{Vp}. 定理 4.1 给定正整数 p (1 ≤ p ≤ n), 考虑算法 4.4 . 假设 A 是可对角化的, 且 |λ1| ≥ · · · ≥ |λp| > |λp+1| ≥ · · · ≥ |λn|. 则 span{Zk} 收敛到 A 的一个 p 维不变子空间. 当 A 不可对角化时, 利用 Jordan 标准型, 我们可以得到同样的结论, 见 [132, 133]. 在正交迭代中, 如果我们取 Z0 = I, 则可得到一类特殊的正交迭代法. 此时, 在一定条件下, 正交迭代会收敛到 A 的 Schur 分解
-136第四讲非对称特征值问题4.4QR选代法OR选代法的基本思想是通过不断的正交相似变换，使得A趋向于一个上三角形式（或拟上三角形式）.算法形式非常简洁，描述如下：算法4.5.QR选代法(QRIteration)l:Set Ai=Aand k-12: while not convergence do%QR分解3:[Qk, Rk] = QR(Ak)4:computeAk+1=R:Qk5:k=k+16: end while在该算法中，我们有Ak+1=RQk=(QTQ)RQk=QT(QR)Qk=QTAQk由这个递推关系可得Ak+1=QTAkQk= QTQT-1Ak-1Qk-1Qk=-.-- =QTQT-1-.QIAQ1-·Qk-1Qk记Qk=Q1.Q-1Qk=[a..,(]则Ak+1=QTAQk,(4.1)即Ak+1与A正交相似.4.4.1OR选代与幂选代的关系记R=RRk-1R1则有QRk=Qk-1(QR)Rk-1=Qk-1(A)Rk-1=Qk-1(QT-,AQk-1)Rk-1=AQk-1Rk-1由此递推下去，即可得QRK= Ak-1QiR1= Ak-1QiR1= Ak(4.2)故QkRke1 = Ahe1,这说明QR选代与幕选代有关假设il>[2|≥≥[nl,则当k充分大时，Aei收敛到A的模最大特征值>所对应的特征向量，故Q的第一列(）也收敛到入1所对应的特征向量。因此，当k充分大时,Ag(k)→^1g(),此时由(4.1)可知,Ak+1的第一列Ak+1(,1) =QTAg()→A1QTq()=A1e1,即Ak+1的第一列的第一个元素收敛到入1，而其它元素都趋向于0,收敛速度取决于[入2/入1|的大小
· 136 · 第四讲 非对称特征值问题 4.4 QR 迭代法 QR 迭代法的基本思想是通过不断的正交相似变换, 使得 A 趋向于一个上三角形式 (或拟上 三角形式). 算法形式非常简洁, 描述如下: 算法 4.5. QR 迭代法 (QR Iteration) 1: Set A1 = A and k = 1 2: while not convergence do 3: [Qk, Rk] = QR(Ak) % QR 分解 4: compute Ak+1 = RkQk 5: k = k + 1 6: end while 在该算法中, 我们有 Ak+1 = RkQk = (Q T k Qk)RkQk = Q T k (QkRk)Qk = Q T k AkQk. 由这个递推关系可得 Ak+1 = Q T k AkQk = Q T k Q T k−1Ak−1Qk−1Qk = · · · = Q T k Q T k−1 · · · Q T 1 AQ1 · · · Qk−1Qk. 记 Q˜ k = Q1 · · · Qk−1Qk = h q˜ (k) 1 , q˜ (k) 2 , . . . , q˜ (k) n i , 则 Ak+1 = Q˜T k AQ˜ k, (4.1) 即 Ak+1 与 A 正交相似. 4.4.1 QR 迭代与幂迭代的关系 记 R˜ k = RkRk−1 · · · R1, 则有 Q˜ kR˜ k = Q˜ k−1(QkRk)R˜ k−1 = Q˜ k−1(Ak)R˜ k−1 = Q˜ k−1(Q˜T k−1AQ˜ k−1)R˜ k−1 = AQ˜ k−1R˜ k−1, 由此递推下去, 即可得 Q˜ kR˜ k = A k−1Q˜ 1R˜ 1 = A k−1Q1R1 = A k . (4.2) 故 Q˜ kR˜ ke1 = A k e1, 这说明 QR 迭代与幂迭代有关. 假设 |λ1| > |λ2| ≥ · · · ≥ |λn|, 则当 k 充分大时, Ak e1 收敛到 A 的模最大特征值 λ1 所 对应的特征向量, 故 Q˜ k 的第一列 q˜ (k) 1 也收敛到 λ1 所对应的特征向量. 因此, 当 k 充分大时, Aq˜ (k) 1 → λ1q˜ (k) 1 , 此时由 (4.1) 可知, Ak+1 的第一列 Ak+1(:, 1) = Q˜T k Aq˜ (k) 1 → λ1Q˜T k q˜ (k) 1 = λ1e1, 即 Ak+1 的第一列的第一个元素收敛到 λ1, 而其它元素都趋向于 0, 收敛速度取决于 |λ2/λ1| 的大 小
4.4QR送代法·137.4.4.2QR选代与反选代的关系下面观察Q的最后一列.由(4.1)可知AQk=QkAk+1=QQk+1Rk+1=Q+1Rk+1,所以有Qk+1=AQR+1由于Qk+1和Q都是正交矩阵，上式两边转置后求逆，可得Qk+1=(QT+1)=((R+)QTAT)=(AT)QRI+1观察等式两边矩阵的最后一列，可得(+1) =c (AT)- (),其中c1为某个常数.依此类推,可知d(+1) =c(AT)-"dl),其中c为某个常数.假设A的特征值满足/>i|≥[>2|≥：≥/入n-1]>[n>0,则入-1是（AT)-的模最大的特征值。由幂迭代的收敛性可知，q(+1）收敛到（AT)-1的模最大特征值入-1所对应的特征向量，即当k充分大时，有(AT)(k+1) →入"g(+1),所以ATa(k+1) → ng(k+1),由（4.1)可知，AT+,的最后一列为AT+1(c,n) = QTATg(k)→ nQTg(k) = Anen即Ak+1的最后一行的最后一个元素收敛到入n，而其它元素都趋向于0,收敛速度取决于/n/入n-1l的大小4.4.3QR选代与正交选代的关系下面的定理给出了QR送代法与正交迭代法（取Zo=I)之间的关系定理4.2设正交选代法4.4和OR算法4.5中所涉及的OR分解都是唯一的.设Ak是由OR迭代法4.5生成的矩阵，Zk是由正交选代法4.4（取Zo=I）生成的矩阵，则有Ak+1=ZTAZk(板书)证明.我们用归纳法证明该结论当 k=0时,A1=A,Zo=I.结论显然成立.设A=ZT-,AZk-1.由于Zk-1是正交矩阵，我们有Z,Rk= Y= AZk-1 = (Zk-1ZT-1)AZk-1=Zk-1Ak =(Zk-1Qk)Rk
4.4 QR 迭代法 · 137 · 4.4.2 QR 迭代与反迭代的关系 下面观察 Q˜ k 的最后一列. 由 (4.1) 可知 AQ˜ k = Q˜ kAk+1 = Q˜ kQk+1Rk+1 = Q˜ k+1Rk+1, 所以有 Q˜ k+1 = AQ˜ kR −1 k+1. 由于 Q˜ k+1 和 Q˜ k 都是正交矩阵, 上式两边转置后求逆, 可得 Q˜ k+1 = Q˜T k+1−1 = R −1 k+1T Q˜T k A T −1 = A T −1 Q˜ kR T k+1. 观察等式两边矩阵的最后一列, 可得 q˜ (k+1) n = c1 A T −1 q˜ (k) n , 其中 c1 为某个常数. 依此类推, 可知 q˜ (k+1) n = c A T −k q˜ (1) n , 其中 c 为某个常数. 假设 A 的特征值满足 |λ1| ≥ |λ2| ≥ · · · ≥ |λn−1| > |λn| > 0, 则 λ −1 n 是 AT −1 的模最大的特征值. 由幂迭代的收敛性可知, q˜ (k+1) n 收敛到 AT −1 的模最大特征值 λ −1 n 所对应 的特征向量, 即当 k 充分大时, 有 A T −1 q˜ (k+1) n → λ −1 n q˜ (k+1) n . 所以 A T q˜ (k+1) n → λnq˜ (k+1) n . 由 (4.1) 可知, AT k+1 的最后一列为 A T k+1(:, n) = Q˜T k A T q˜ (k) n → λnQ˜T k q˜ (k) n = λnen, 即Ak+1 的最后一行的最后一个元素收敛到λn, 而其它元素都趋向于0, 收敛速度取决于|λn/λn−1| 的大小. 4.4.3 QR 迭代与正交迭代的关系 下面的定理给出了 QR 迭代法与正交迭代法 (取 Z0 = I) 之间的关系. 定理 4.2 设正交迭代法 4.4 和 QR 算法 4.5 中所涉及的 QR 分解都是唯一的. 设 Ak 是由 QR 迭 代法 4.5 生成的矩阵, Zk 是由正交迭代法 4.4 (取 Z0 = I) 生成的矩阵, 则有 Ak+1 = Z T k AZk. (板书) 证明. 我们用归纳法证明该结论. 当 k = 0 时, A1 = A, Z0 = I. 结论显然成立. 设 Ak = Z T k−1AZk−1. 由于 Zk−1 是正交矩阵, 我们有 ZkRˆ k = Yk = AZk−1 = (Zk−1Z T k−1 )AZk−1 = Zk−1Ak = (Zk−1Qk)Rk
-138第四讲非对称特征值问题即ZR和（Zk-1Q)R都是Y的QR分解.由QR分解的唯一性可知，Z=Zk-1QkR=Rk所以ZTAZK=(Zk-1Qh)TA(Zk-1Qk)=QTAkQk=QT(QkR)Qk=RKQk= Ak+1)口即Ak+1=ZTAZk.由归纳法可知，定理结论成立.4.4.4QR选代的收敛性定理43 设A=VAV-1 e Rnxn,其中 A=diag(>1,>2,...,入n),且[>i|>[2|>..>[n|>0若V-1的所有顺序主子矩阵都非奇异（即V-1存在LU分解），则A的对角线以下的元素均收敛到0.(板书)证明.设V=Q,R,是V的QR分解，V-1=LU，是V-1的LU分解，其中L是单位下三角矩阵.则Ak=VA*V-1=Q,RA*L,U=QR(A*LA-)A*U注意到矩阵ALA-k是一个下三角矩阵，且其（i）位置上的元素为[0,ij.由于当>时有/入<1,故当充分大时，/趋向于0.所以我们可以把^L-写成A*LuA-k = I + Ek,其中E满足limE=0.于是A=QuR(I +E)A*U=Q(I+R,ER-1)RA*U(4.3)对矩阵I+R,ER1做QR分解：I+R,ER=1=QERE·由于Ek→0,所以QE→I,REx→I.将其代人(22)可得A=QQERERA*Uu=QQEDk(D-'RERuA*U.),(4.4)其中Dk是对角矩阵，其对角线元素的模均为1，它使得上三角矩阵D-1RE,R,A*U，的对角线元素均为正.这样，(22)就构成A的QR分解.又由（4.2)可知A=QRk，根据QR分解的唯一性，我们可得Qu=QQEDk，RK=D-'RERAUu所以由(4.1)可知Ak+1=QTAQk=(QQED)TVAV-1QQE,Dk=DTQEQTQ,RAR-'Q-QQEDk=DIQERuAR-'QEDk
· 138 · 第四讲 非对称特征值问题 即 ZkRˆ k 和 (Zk−1Qk)Rk 都是 Yk 的 QR 分解. 由 QR 分解的唯一性可知, Zk = Zk−1Qk, Rˆ k = Rk. 所以 Z T k AZk = (Zk−1Qk) TA(Zk−1Qk) = Q T k AkQk = Q T k (QkRk)Qk = RkQk = Ak+1, 即 Ak+1 = Z T k AZk. 由归纳法可知, 定理结论成立. □ 4.4.4 QR 迭代的收敛性 定理 4.3 设 A = V ΛV −1 ∈ R n×n , 其中 Λ = diag(λ1, λ2, . . . , λn), 且 |λ1| > |λ2| > · · · > |λn| > 0. 若 V −1 的所有顺序主子矩阵都非奇异（即 V −1 存在 LU 分解）, 则 Ak 的对角线以下的元素均 收敛到 0. (板书) 证明. 设 V = QvRv 是 V 的 QR 分解, V −1 = LvUv 是 V −1 的 LU 分解, 其中 Lv 是单位下三角矩 阵. 则 A k = V Λ kV −1 = QvRvΛ kLvUv = QvRv(ΛkLvΛ −k )ΛkUv. 注意到矩阵 Λ kLvΛ −k 是一个下三角矩阵, 且其 (i, j) 位置上的元素为 Λ kLvΛ −k (i, j) = 0, i j. 由于当 i > j 时有 |λi/λj | < 1, 故当 k 充分大时, λ k i /λ k j 趋向于 0. 所以我们可以把 Λ kLvΛ −k 写成 Λ kLvΛ −k = I + Ek, 其中 Ek 满足 lim k→∞ Ek = 0. 于是 A k = QvRv(I + Ek)ΛkUv = Qv(I + RvEkR −1 v )RvΛ kUv. (4.3) 对矩阵 I + RvEkR−1 v 做 QR 分解: I + RvEkR−1 v = QEkREk . 由于 Ek → 0, 所以 QEk → I, REk → I. 将其代入 (??) 可得 A k = QvQEkREkRvΛ kUv = QvQEkDk(D −1 k REkRvΛ kUv), (4.4) 其中 Dk 是对角矩阵, 其对角线元素的模均为 1, 它使得上三角矩阵 D −1 k REkRvΛ kUv 的对角线元 素均为正. 这样, (??) 就构成 Ak 的 QR 分解. 又由 (4.2) 可知 Ak = Q˜ kR˜ k, 根据 QR 分解的唯一性, 我们可得 Q˜ k = QvQEkDk, R˜ k = D −1 k REkRvΛ kUv. 所以由 (4.1) 可知 Ak+1 = Q˜T k AQ˜ k = (QvQEkDk) TV ΛV −1QvQEkDk = DT k Q T EkQ T v QvRvΛR −1 v Q −1 v QvQEkDk = DT k Q T EkRvΛR −1 v QEkDk