第四讲非对称特征值问题设AERnxn是一个非对称的稠密矩阵，本章主要讨论如何计算A的全部特征值和特征向量.记A的特征值为入1，入2，.….，入n.本讲中我们总是假定[1/ ≥ [α2] ≥... ≥[n/ ≥0,即A的特征值按绝对值（模）降序排列本章主要介绍以下方法：·幂选代算法·位移策略与反选代技巧·正交送代法·QR算法与实用的QR算法关于稠密矩阵特征值计算的参考资料有：·J.H.Wilkinson,TheAlgebraicEigenvalueProblem,1965[75](有中文翻译)B.N.Parlett,TheSymmetricEigenvalueProblem,2ndEds., 1998[57].G. W. Stewart, Matrix Algorithms, Vol II: Eigensystems, 2001 [67]. G.H.Golub and C.F. Van Loan, Matrix Computations, 2013 [33]P. Arbenz, The course 252-0504-00 G, Numerical Methods for Solving Large Scale Eigenvalue Problems, 2014.4.1幂选代方法幂送代方法是计算特征值和特征向量的一种简单易用的算法，幂选代虽然简单，但它却建立了计算特征值和特征向量的算法的一个基本框架.算法4.1.幂选代算法(PowerIteration)I: Choose an initial guess r(0) with (0) l2 = 12: set k = 03:while not convergencedog(k+1) = Az(k)4:2(+1) = y(#+1) /ly(k+1)I25:μk+1=(2(+1), Ar(k+1))%内积6:111
第四讲 非对称特征值问题 设 A ∈ R n×n 是一个非对称的稠密矩阵, 本章主要讨论如何计算 A 的全部特征值和特征向量. 记 A 的特征值为 λ1, λ2, . . . , λn. 本讲中我们总是假定 |λ1| ≥ |λ2| ≥ · · · ≥ |λn| ≥ 0, 即 A 的特征值按绝对值（模）降序排列. 本章主要介绍以下方法： • 幂迭代算法 • 位移策略与反迭代技巧 • 正交迭代法 • QR 算法与实用的 QR 算法 关于稠密矩阵特征值计算的参考资料有： • J. H. Wilkinson, The Algebraic Eigenvalue Problem, 1965 [75] (有中文翻译) • B. N. Parlett, The Symmetric Eigenvalue Problem, 2nd Eds., 1998 [57] • G. W. Stewart, Matrix Algorithms, Vol II: Eigensystems, 2001 [67] • G. H. Golub and C. F. Van Loan, Matrix Computations, 2013 [33] • P. Arbenz, The course 252-0504-00 G, Numerical Methods for Solving Large Scale Eigenvalue Problems, 2014. 4.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) 5: x (k+1) = y (k+1)/∥y (k+1)∥2 6: µk+1 = (x (k+1), Ax(k+1)) % 内积 111
第四讲非对称特征值问题.112.k=k+17:8: end while下面讨论幂选代的收敛性假设(1)ARnxn是可对角化的，即A=VAV-1,其中A=diag(A1,2,...,An）ECnxn,=[u1,u2,..*, n] e Cnxn, 且 lluill2 = 1, i = 1,2,...,n.(2)同时,我们还假设[|>[2|≥[入3|≥.≥[入]由于VECnxn非奇异，所以它的列向量组构成Cn的一组基.因此迭代初始向量r(0)可表示为2(0) = Q1U1 +Q22 + ... + QnUn = V[a1,Q2,..,an]T.我们假定α1≠0,即z(0)不属于span[u2,U3,...,Un)（由于z(0)是随机选取的，从概率意义上讲,这个假设通常是成立的).于是我们可得1[Q1[ai[()a21(a=VAAA*(0) = (VAV-1)KV..:.anahanan又[/入12/入1|的大小，[>2/入1|越小，收敛越快通过上面的分析可知，幂迭代只能用于计算矩阵的模最大的特征值和其相应的特征向量，如果A的模最大的特征值是唯一的，则称其为主特征值.当^2/入1接近于1时，收敛速度会非常慢.同时，如果A的模最大特征值是一对共轭复数，则幂选代就可能就会失效4.2反选代方法前面已经提到，幂选代算法的收敛速度取决于/入2/入|的大小，当它的值接近于1时，收敛速度会非常缓慢.因此，为了加快幂送代算法的收敛速度，我们希望入2/入1的值越小越好个简单易用的方法就是使用位移策略，即将计算A的特征值转化为计算A－αI的特征值，即对A做一个移位.这里α是一个给定的数，称α为位移(shift).为了使得幂选代作用到A一αI时具有更快的收敛速度，我们要求。满足下面的两个条件：
· 112 · 第四讲 非对称特征值问题 7: k = k + 1 8: end while 下面讨论幂迭代的收敛性. 假设 (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] ⊺ . 我们假定 α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 ]⊺ , k = 0, 1, 2, . . . 收敛到 e1 = [1, 0, . . . , 0]⊺ . 所以向量 x (k) = Akx (0)/∥Akx (0)∥2 收敛到 ±v1, 即 A 的对应于 (模) 最大的 特征值 λ1 的特征向量. 而 µk = (x (k) ) ∗Ax(k) 则收敛到 v ∗ 1Av1 = λ1. 显然, 幂迭代的收敛快慢取决于 |λ2/λ1| 的大小, |λ2/λ1| 越小, 收敛越快. 通过上面的分析可知, 幂迭代只能用于计算矩阵的模最大的特征值和其相应的特征向量. 如果 A 的模最大的特征值是唯一的, 则称其为 主特征值. 当 |λ2/λ1| 接近于 1 时, 收敛速度会非常慢. 同时, 如 果 A 的模最大特征值是一对共轭复数, 则幂迭代就可能就会失效. 4.2 反迭代方法 前面已经提到, 幂迭代算法的收敛速度取决于 |λ2/λ1| 的大小. 当它的值接近于 1 时, 收敛速度会非 常缓慢. 因此, 为了加快幂迭代算法的收敛速度, 我们希望 |λ2/λ1| 的值越小越好. 一个简单易用的方法就是使用位移策略, 即将计算 A 的特征值转化为计算 A − σI 的特征值, 即对 A 做一个移位. 这里 σ 是一个给定的数, 称 σ 为位移 (shift). 为了使得幂迭代作用到 A − σI 时具有更快 的收敛速度, 我们要求 σ 满足下面的两个条件:
4.2反选代方法.113.(1)入1-α是A-αI的模最大的特征值；i-0尽可能地小(2)22icnl其中第一个条件保证最后所求得的特征值是我们所要的，第二个条件用于加快幂送代的收敛速度显然，在实际应用中，的取值并不是一件容易的事，例4.1设A=XAX-1,其中A为对角矩阵，分别用幂迭代方法和带位移的幂迭代方法计算A的主特征值(见Eig_Power_shift.m)十位移策略在特征值计算中非常重要，主要是用在反选代算法和QR选代算法中4.2.1反选代算法如果我们将幂代算法作用在A-1上，则可求出A的模最小的特征值.事实上，结合这种思想和位移策略，我们就可以计算矩阵的任意一个特征值算法4.2.反选代算法(lInverseIteration)1: Choose a scalar g and an initial vector (0) with (0)ll2 = 12: set k = 03:while not convergence dog(k+1) = (A - gI)-1(k)4:2(+1) = gy(++1) / lg(k+1 125:μk+1 = (r(k+1), Ar(k+1)6:7:k=k+18:end while该算法称为反送代算法.显然，在反选代算法中，收敛到距离最近的特征值，而（）则收敛到其对应的特征向量设距离最近的特征值为入，则算法的收敛速度取决于入01-0的大小显然，α约接近于入，其值越小，即算法收敛越快若α～入，则选代几步就可以了反迭代算法的另一个优点是，只要选取合适的位移，就可以计算A的任意一个特征值但反选代算法的缺点也很明显：每步选代需要解一个线性方程组（A－αI)y(+1)=r(k)，这就需要对A一αI做一次LU分解，另外，与幂选代一样，反选代算法一次只能求一个特征值，4.2.2Rayleigh商选代在反迭代算法中，有一个很重要的问题，就是位移。的选取
4.2 反迭代方法 · 113 · (1) λ1 − σ 是 A − σI 的模最大的特征值; (2) max 2≤i≤n λi − σ λ1 − σ 尽可能地小. 其中第一个条件保证最后所求得的特征值是我们所要的, 第二个条件用于加快幂迭代的收敛速度. 显然, 在实际应用中, σ 的取值并不是一件容易的事. 例 4.1 设 A = XΛX−1 , 其中 Λ 为对角矩阵, 分别用幂迭代方法和带位移的幂迭代方法计算 A 的主特 征值. (见 Eig_Power_shift.m) † 位移策略在特征值计算中非常重要, 主要是用在反迭代算法和 QR 迭代算法中. 4.2.1 反迭代算法 如果我们将幂迭代算法作用在 A−1 上, 则可求出 A 的模最小的特征值. 事实上, 结合这种思想和位 移策略, 我们就可以计算矩阵的任意一个特征值. 算法 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 = (x (k+1), Ax(k+1)) 7: k = k + 1 8: end while 该算法称为反迭代算法. 显然, 在反迭代算法中, µk 收敛到距离 σ 最近的特征值, 而 x (k) 则收敛到 其对应的特征向量. 设距离 σ 最近的特征值为 λk, 则算法的收敛速度取决于 max 1≤i≤n λk − σ λi − σ 的大小. 显然, σ 约接近于 λk, 其值越小, 即算法收敛越快. 若 σ ≈ λk, 则迭代几步就可以了. 反迭代算法的另一个优点是, 只要选取合适的位移 σ, 就可以计算 A 的任意一个特征值. 但反迭代算法的缺点也很明显: 每步迭代需要解一个线性方程组 (A − σI)y (k+1) = x (k) , 这就需要 对 A − σI 做一次 LU 分解. 另外, 与幂迭代一样, 反迭代算法一次只能求一个特征值. 4.2.2 Rayleigh 商迭代 在反迭代算法中, 有一个很重要的问题, 就是位移 σ 的选取
.114.第四讲非对称特征值问题显然，选取的基本原则是使得其与所求的特征值越靠近越好这里需要指出的是，在每步迭代中位移可以不一样，即在迭代过程中可以选取不同的位移由于在反迭代算法4.2中，是收敛到所求的特征值的，所以我们可以选取作为第步的位移这时所得到的反选代算法就称为Rayleigh商选代（RayleighQuotientIteration），简记为RQI算法4.3.Rayleigh商送代算法(RayleighQuotientIteration(RQD)1: Choose an initial vector (0)with [(0)12 = 12:setk=03: compute g = ((0))*Ar(0)4: while not converge doy(h+1) = (A - αI)-1r(h)5:(k+1) = 3y(k+1) /lg(k+1) 26:μk+1 = ((k+1), Ar(+1)7:8:=μk+19:k=k+110:end while一般来说，如果Rayleigh商迭代收敛到A的一个单特征值，则至少是二次收敛的，即具有局部二次收敛性.如果A是对称的，则能达到局部三次收敛在Rayleigh商迭代中，由于每次迭代的位移是不同的，因此每次迭代需要求解一个不同的线性方程组，这使得运算量大大增加例4.2设A=XAX-1,其中A为对角矩阵，用Rayleigh商选代计算A的特征值(见Eig_Rayleigh.m)4.3正交选代方法幂选代和反选代都只能同时计算一个特征对，如果想同时计算多个特征对，我们可以采用多个初始向量进行选代，而正交选代算法就是基于这种思想，它能够计算A的一个不变子空间，从而可以同时计算出多个特征值算法4.4.正交选代算法(OrthogonalIteration):Choosean n ×pcolumn orthogonal matrixZo2:setk=03: while not convergence do4:compute Yk+1=AZkYk+1=Zk+1Rk+1%QR分解5:
· 114 · 第四讲 非对称特征值问题 显然, 选取 σ 的基本原则是使得其与所求的特征值越靠近越好. 这里需要指出的是, 在每步迭代中, 位移 σ 可以不一样, 即在迭代过程中可以选取不同的位移 σ. 由于在反迭代算法 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 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 = (x (k+1), Ax(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) 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 分解
4.4QR选代方法.115.k=k+16:7: end while在算法中使用OR分解是为了保持Zk的列正交性，使得其列向量组构成子空间spanA'Zo!的一组正交基，以确保算法的数值稳定性，下面我们分析该算法的收敛性质。假设A是可对角化的，即A=VAV-1,其中Λ=diag(入1，入2,..,入),且[≥..≥[>p/>[p+1/≥.≥[>则可得span[Z] = span[Y] = span[AZk-1],k =1, 2,...由此可知span[Z] = span[A*Zo] = span[VA^V-1 Zo]我们注意到[(Ai/p)k[w(k)V-1Zo=入V-1Zo≤入[w(k)n(An/入p)由于当i>时有[入/入pl,[Vp, Vn-p]W(k)所以当k→8时，有span[Z,] = span[VA*V-1 Zo] = span[V,W(k) + Vn-pW(-),]→ span[V,W('} = span[Vp]即span(Z]趋向于A的一个p维不变子空间span(V,]定理4.1给定正整数p(1≤P≤n)，考虑算法4.4.假设A是可对角化的，且/>il≥..·≥Ipl>[p+1/≥..≥[^nl.则span[Z}收敛到A的一个p维不变子空间当A不可对角化时，利用Jordan标准型，我们可以得到同样的结论，见[73，74]在正交迭代中，如果我们取Zo=I，则可得到一类特殊的正交迭代算法.此时，在一定条件下，正交选代会收敛到A的Schur标准型4.4QR送代方法QR迭代算法的基本思想是通过不断的正交变换，将A转化为上三角形式(或拟上三角形式).算法形式非常简单，描述如下：
4.4 QR 迭代方法 · 115 · 6: k = k + 1 7: end while 在算法中使用 QR 分解是为了保持 Zk 的列正交性, 使得其列向量组构成子空间 span{AiZ0} 的一 组正交基, 以确保算法的数值稳定性. 下面我们分析该算法的收敛性质. 假设 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| |λp+1| ≥ · · · ≥ |λn|. 则 span{Zk} 收敛到 A 的一个 p 维不变子空间. 当 A 不可对角化时, 利用 Jordan 标准型, 我们可以得到同样的结论, 见 [73, 74]. 在正交迭代中, 如果我们取 Z0 = I, 则可得到一类特殊的正交迭代算法. 此时, 在一定条件下, 正交 迭代会收敛到 A 的 Schur 标准型. 4.4 QR 迭代方法 QR 迭代算法的基本思想是通过不断的正交变换, 将 A 转化为上三角形式 (或拟上三角形式). 算法 形式非常简单, 描述如下:
.116.第四讲非对称特征值问题算法4.5.QR选代算法（QRIteration)1: Set Ar = A and k=12:while not convergence do%QR分解3:A=QRk4:compute Ak+1 = R,Qk5:k=k+16:endwhile在该算法中，我们有Ak+1=RQk=(QQ)RQ=Q(QR)Qk=QAQk由这个递推关系可得Ak+1=QLA+Qk=QIQ-,Ak-1Qk-1Q=...=QIQ1-1*--QIAQ1--.Qk-1Qk记Q=Q1.Q-1Q=a*,.,]则Ak+1=QTAQk(4.1)即Ak+1与A正交相似4.4.1QR送代与幂选代的关系记R=RRk-1..R1,则有Q,Rk=Qk-1(QRk)Rk-1=Qk-1(A)Rk-1=Qk-1(Q1-,AQk-1)Rk-1=AQk-1Rk-1:由此递推下去，即可得Q,R=Ak-1Q,Ri=Ak-1QiR1=Ak(4.2)故Q,Rkei = Akei,这说明QR送代与幂选代有关假设il>[2|≥..≥/nl,则当充分大时，Ae1收敛到A的模最大特征值入对应的特征向量，故Q的第一列g()也收敛到入1对应的特征向量.因此，当充分大时，Ag(）→入1g(),此时由(4.1)可知，Ak+1的第一列为A+1(,1)=QIAa()→AQa()=1e1,即Ak+1的第一列的第一个元素收敛到入1，而其它元素都趋向于0,收敛速度取决于/2/入1|的大小
· 116 · 第四讲 非对称特征值问题 算法 4.5. QR 迭代算法 (QR Iteration) 1: Set A1 = A and k = 1 2: while not convergence do 3: Ak = QkRk % QR 分解 4: compute Ak+1 = RkQk 5: k = k + 1 6: end while 在该算法中, 我们有 Ak+1 = RkQk = (Q ⊺ kQk)RkQk = Q ⊺ k (QkRk)Qk = Q ⊺ kAkQk. 由这个递推关系可得 Ak+1 = Q ⊺ kAkQk = Q ⊺ kQ ⊺ k−1Ak−1Qk−1Qk = · · · = Q ⊺ kQ ⊺ k−1 · · · Q ⊺ 1AQ1 · · · Qk−1Qk. 记 Q˜ k = Q1 · · · Qk−1Qk = [ q˜ (k) 1 , q˜ (k) 2 , . . . , q˜ (k) n ] , 则 Ak+1 = Q˜⊺ kAQ˜ 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˜⊺ 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˜⊺ kAq˜ (k) 1 → λ1Q˜⊺ k q˜ (k) 1 = λ1e1, 即 Ak+1 的第一列的第一个元素收敛到 λ1, 而其它元素都趋向于 0, 收敛速度取决于 |λ2/λ1| 的大小
4.4QR选代方法.117.4.4.2QR选代与反选代的关系下面观察Qk的最后一列.由（4.1）可知AQk=QkAk+1=QQk+1Rk+1=Qk+1Rk+1所以有Qk+1=AQRt+1由于Qk+1和Q都是正交矩阵，上式两边转置后求逆，可得Qk+1 = (Q1+)- = (R+)QIAt)- = (AT)-" QxRI+1观察等式两边矩阵的最后一列，可得G(k+1) = C1 (AT)- a(k),其中c1为某个常数，依此类推、可知g(k+1) = c(AT)-k g(1),其中c为某个常数.假设A的特征值满足[i|≥2|≥.≥[n-1>[n|>0,则>=1是(AT)-1的模最大的特征值.由幂选代的收敛性可知,G(+1）收敛到（AT)-的模最大特征值入"所对应的特征向量，即当充分大时，有(AT)- g(k+1) → -"g(k+1)所以ATa(k+1) →入ng(k+1),由(4.1)可知，AI+1的最后一列为AI+1(c,n)= QIATa(k)→ AnQIa(k) = Anen即Ak+1的最后一行的最后一个元素收敛到入n，而其它元素都趋向于0,收敛速度取决于[/入n-1/的大小，4.4.3QR选代与正交选代的关系下面的定理给出了QR选代算法与正交代算法（取Zo=D)之间的关系定理4.2设正交选代算法4.4和QR算法4.5中所涉及的QR分解都是唯一的.Ak是由QR选代算法4.5生成的矩阵，Zk是由正交选代算法4.4（取Zo=I)生成的矩阵，则有Ak+1 = ZIAZk.证明.我们用归纳法证明该结论当k=0时，A1=A,Zo=I.结论显然成立
4.4 QR 迭代方法 · 117 · 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˜⊺ k+1)−1 = (( R −1 k+1)⊺ Q˜⊺ kA ⊺ )−1 = (A ⊺ ) −1 Q˜ kR ⊺ k+1. 观察等式两边矩阵的最后一列, 可得 q˜ (k+1) n = c1 (A ⊺ ) −1 q˜ (k) n , 其中 c1 为某个常数. 依此类推, 可知 q˜ (k+1) n = c (A ⊺ ) −k q˜ (1) n , 其中 c 为某个常数. 假设 A 的特征值满足 |λ1| ≥ |λ2| ≥ · · · ≥ |λn−1| > |λn| > 0, 则 λ −1 n 是 (A⊺ ) −1 的模 最大的特征值. 由幂迭代的收敛性可知, q˜ (k+1) n 收敛到 (A⊺ ) −1 的模最大特征值 λ −1 n 所对应的特征向量, 即当 k 充分大时, 有 (A ⊺ ) −1 q˜ (k+1) n → λ −1 n q˜ (k+1) n . 所以 A ⊺ q˜ (k+1) n → λnq˜ (k+1) n . 由 (4.1) 可知, A ⊺ k+1 的最后一列为 A ⊺ k+1(:, n) = Q˜⊺ kA ⊺ q˜ (k) n → λnQ˜⊺ 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 ⊺ kAZk. 证明. 我们用归纳法证明该结论. 当 k = 0 时, A1 = A, Z0 = I. 结论显然成立
第四讲非对称特征值问题.118.设Ak=ZT-AZk-1.由于Zk-1是正交矩阵，我们有Z,R=Y = AZk-1=(Zk-1ZT-1)AZk-1= Zk-1A, =(Zk-1Qk)Rk即ZR和（Zk-1Qk)R都是Y的QR分解.由QR分解的唯一性可知Zk=Zk-1Qk，R=Rk.所以ZIAZk=(Zk-1Q)TA(Zk-1Q)=QLAQ=QI(QR)Qh=RQk=Ak+1)口即Ak+1=ZLAZk.由归纳法可知，定理结论成立.4.4.4QR选代的收敛性定理43设A=VAV-1 Rnxn,其中A=diag(入1,>2...,入n),且[i/>[2/ >...>[nl.若V-1的所有主子矩阵都非奇异（即V-1存在LU分解），则A的对角线以下的元素均收敛到0证明.设V=QR是V的QR分解，V-1=LU，是V-1的LU分解，其中L是单位下三角矩阵.则A*=VA*V-1=QuRwA*LuUu=QuR(A*LuA-k)A*Uu.注意到矩阵A^LA-k是一个下三角矩阵，且其（i，）位置上的元素为(o,ij.由于当i>时有[/l<1，故当充分大时，/入趋向于0.所以我们可以把-写成A*LuA-* = I + Ek,其中E满足limE=0.于是Ah=QR(I+Ex)A*U,=Qu(I+RER-")RuA*Uu(4.3)对矩阵I+RER-1做QR分解：I+R,ER-1=QRE.由于E→0,所以Q→I,RE→I将其代人(4.3)可得Ah=QuQERER,A*Uu=Q,QED(D-'RERuA*U,),(4.4)其中Dk是对角矩阵，其对角线元素的模均为1，它使得上三角矩阵D-1ReR,AU，的对角线元素均为正.这样，（4.4)就构成A*的QR分解。又由（4.2)可知A*=QRk，根据QR分解的唯一性，我们可得Q=QQEDk,R=D-'RER,A*Uu.所以由(4.1)可知Ak+1= QIAQk
· 118 · 第四讲 非对称特征值问题 设 Ak = Z ⊺ k−1AZk−1. 由于 Zk−1 是正交矩阵, 我们有 ZkRˆ k = Yk = AZk−1 = (Zk−1Z ⊺ k−1 )AZk−1 = Zk−1Ak = (Zk−1Qk)Rk, 即 ZkRˆ k 和 (Zk−1Qk)Rk 都是 Yk 的 QR 分解. 由 QR 分解的唯一性可知, Zk = Zk−1Qk, Rˆ k = Rk. 所以 Z ⊺ kAZk = (Zk−1Qk) ⊺A(Zk−1Qk) = Q ⊺ kAkQk = Q ⊺ k (QkRk)Qk = RkQk = Ak+1, 即 Ak+1 = Z ⊺ kAZk. 由归纳法可知, 定理结论成立. □ 4.4.4 QR 迭代的收敛性 定理 4.3 设 A = V ΛV −1 ∈ R n×n, 其中 Λ = diag(λ1, λ2, . . . , λn), 且 |λ1| > |λ2| > · · · > |λn|. 若 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. 将其代入 (4.3) 可得 A k = QvQEkREkRvΛ kUv = QvQEkDk(D −1 k REkRvΛ kUv), (4.4) 其中 Dk 是对角矩阵, 其对角线元素的模均为 1, 它使得上三角矩阵 D −1 k REkRvΛ kUv 的对角线元素均 为正. 这样, (4.4) 就构成 Ak 的 QR 分解. 又由 (4.2) 可知 Ak = Q˜ kR˜ k, 根据 QR 分解的唯一性, 我们可得 Q˜ k = QvQEkDk, R˜ k = D −1 k REkRvΛ kUv. 所以由 (4.1) 可知 Ak+1 = Q˜⊺ kAQ˜ k
4.4QR选代方法.119.=(Q.QEDa)TVAV-1QuQE,D=DIQEQTQR,AR-'Q-Q.QE,Dk=DIQE,RAR-'QEDk由于QE→I，所以当k→o时，Ak+1收敛到一个上三角矩阵.收敛速度取决于入i+1max口的大小，需要指出的是，由于D的元素不一定收敛，故A+1对角线以上（不含对角线）的元素不一定收敛，但这不妨碍Ak+1的对角线元素收敛到A的特征值（即Ak+1的对角线元素是收敛的）例4.3QR迭代算法演示.设95r-1A=X31其中X是由MATLAB随机生成的非奇异矩阵.在选代过程中，对于A的下三角部分中元素，如果其绝对值小于某个阈值tol，则直接将其设为0,即a=0 if i>jand [a|<tol.这里我们取tol=10-6max, (lag,'n.(见Eig_QR.m)1≤i.j≤r4.4.5带位移的QR送代为了加快QR送代的收敛速度，我们可以采用位移策略和反选代思想.算法4.6.带位移的QR选代算法（QRIterationwithshift）l;Set Ai=Aand k =12:while not convergencedo3:Choose a shiftakA-oI=QRk%QR分解4:5:ComputeAk+1=RQk+0klk=k+16:7: end while与不带位移的QR代一样，我们有A+1=RQk+0I=(QQ)RQk+I=Q(Ak-0)Q+0I=QLAQ
4.4 QR 迭代方法 · 119 · = (QvQEkDk) ⊺V ΛV −1QvQEkDk = D ⊺ kQ ⊺ Ek Q ⊺ vQvRvΛR −1 v Q −1 v QvQEkDk = D ⊺ kQ ⊺ Ek RvΛR −1 v QEkDk. 由于 QEk → I, 所以当 k → ∞ 时, Ak+1 收敛到一个上三角矩阵. 收敛速度取决于 max 1≤i j and |a (k) ij | < tol. 这里我们取 tol = 10−6 max 1≤i,j≤n {|a (k) ij |}. (见 Eig_QR.m) 4.4.5 带位移的 QR 迭代 为了加快 QR 迭代的收敛速度, 我们可以采用位移策略和反迭代思想. 算法 4.6. 带位移的 QR 迭代算法 (QR Iteration with shift) 1: Set A1 = A and k = 1 2: while not convergence do 3: Choose a shift σk 4: Ak − σkI = QkRk % QR 分解 5: Compute Ak+1 = RkQk + σkI 6: k = k + 1 7: end while 与不带位移的 QR 迭代一样, 我们有 Ak+1 = RkQk + σkI = (Q ⊺ kQk)RkQk + σkI = Q ⊺ k (Ak − σkI)Qk + σkI = Q ⊺ kAkQk
第四讲非对称特征值问题.120所以，带位移的QR算法中所得到的矩阵A仍然与A1=A正交相似在带位移的QR迭代算法中，一个很重要的问题就是位移k的选取，在前面的分析中我们已经知道，As+1(n，n）将收敛到A的模最小的特征值，且收敛速度取决于模最小特征值与模第二小特征值之间的比值.显然，若k就是A的一个特征值，则A一kI的模最小特征值为0,故QR算法选代一步就收敛.此时[A(n-1)x(n-1)A+1=RQ+I:0(a-1)x(n-1)使用带位移的QR选代算法如果需要计算A的其它特征值，则可对子矩阵A通常，如果k与A的某个特征值非常接近，则收敛速度通常会很快.由于A(n，n）收敛到A的一个特征值，所以在实际使用中，一个比较直观的位移选择策略是k=Ak(（n，n).事实上，这样的位移选取方法通常会使得OR选代有二次收敛速度例4.4带位移的QR选代算法演示.所有数据与设置与例4.3相同，在选代过程中，取ak=A（n,n).如果A(n,n)已经收敛，则取k=A(n-1,n-1)(见Eig_QR_shift.m)4.5带位移的隐式QR选代QR选代算法中需要考虑的另一个重要问题就是运算量：每一步送代都需要做一次QR分解和矩阵乘积，运算量为O(n3).即使每计算一个特征值只需选代一步，则计算所有特征值也需要(n4）的运算量.这是令人无法忍受的.下面我们就想办法将总运算量从On*）减小到O(n3)为了实现这个目标，我们需要利用Hessenberg矩阵.具体步骤如下：首先通过相似变化将A转化成一个上Hessenberg矩阵，然后再对这个Hessenberg矩阵实施隐式QR选代.所谓隐式QR迭代，就是在QR选代中，我们不需要进行显式的QR分解.这样就可以将QR选代的每一步运算量从O(n3）降低到O（n2).从而将总的运算量降低到O(n3)4.5.1上Hessenberg矩阵设H=（hi)ERnxn,若当i>j+1时,有hii=0,则称H为上Hessenberg矩阵定理4.4设AeRnxn，则存在正交矩阵QeIRnxn，使得QAQT是上Hessenberg矩阵下面我们给出具体的转化过程.这里我们主要使用Householder变换：对任意向量rERn，总存在个Householder矩阵H=I-BuvT使得Hr=rl2e1,其中e1=[1,0,...,0]T.我们以一个5×5的矩阵A为例第一步：令Q1=diag(Iix1,Hi),其中Hi是对应于向量A(2：5,1)的Householder矩阵.于是可得Q1A=0
· 120 · 第四讲 非对称特征值问题 所以, 带位移的 QR 算法中所得到的矩阵 Ak 仍然与 A1 = A 正交相似. 在带位移的 QR 迭代算法中, 一个很重要的问题就是位移 σk 的选取. 在前面的分析中我们已经知 道, Ak+1(n, n) 将收敛到 A 的模最小的特征值, 且收敛速度取决于模最小特征值与模第二小特征值之间 的比值. 显然, 若 σk 就是 A 的一个特征值, 则 Ak − σkI 的模最小特征值为 0, 故 QR 算法迭代一步就收 敛. 此时 Ak+1 = RkQk + σkI = [ A (n−1)×(n−1) k+1 ∗ 0 σk ] . 如果需要计算 A 的其它特征值, 则可对子矩阵 A (n−1)×(n−1) k+1 使用带位移的 QR 迭代算法. 通常, 如果 σk 与 A 的某个特征值非常接近, 则收敛速度通常会很快. 由于 Ak(n, n) 收敛到 A 的一 个特征值, 所以在实际使用中, 一个比较直观的位移选择策略是 σk = Ak(n, n). 事实上, 这样的位移选 取方法通常会使得 QR 迭代有二次收敛速度. 例 4.4 带位移的 QR 迭代算法演示. 所有数据与设置与例 4.3 相同, 在迭代过程中, 取 σk = Ak(n, n). 如果 Ak(n, n) 已经收敛, 则取 σk = Ak(n − 1, n − 1). (见 Eig_QR_shift.m) 4.5 带位移的隐式 QR 迭代 QR 迭代算法中需要考虑的另一个重要问题就是运算量: 每一步迭代都需要做一次 QR 分解和矩阵 乘积, 运算量为 O(n 3 ). 即使每计算一个特征值只需迭代一步, 则计算所有特征值也需要 O(n 4 ) 的运算 量. 这是令人无法忍受的. 下面我们就想办法将总运算量从 O(n 4 ) 减小到 O(n 3 ). 为了实现这个目标, 我们需要利用 Hessenberg 矩阵. 具体步骤如下: 首先通过相似变化将 A 转化成 一个上 Hessenberg 矩阵, 然后再对这个 Hessenberg 矩阵实施隐式 QR 迭代. 所谓隐式 QR 迭代, 就是在 QR 迭代中, 我们不需要进行显式的 QR 分解. 这样就可以将 QR 迭代的每一步运算量从 O(n 3 ) 降低到 O(n 2 ). 从而将总的运算量降低到 O(n 3 ). 4.5.1 上 Hessenberg 矩阵 设 H = (hij ) ∈ R n×n, 若当 i > j + 1 时, 有 hij = 0, 则称 H 为 上 Hessenberg 矩阵. 定理 4.4 设 A ∈ R n×n, 则存在正交矩阵 Q ∈ R n×n, 使得 QAQ⊺ 是上 Hessenberg 矩阵. 下面我们给出具体的转化过程. 这里我们主要使用 Householder 变换: 对任意向量 x ∈ R n, 总存在 一个 Householder 矩阵 H = I − βvv⊺ 使得 Hx = ∥x∥2e1, 其中 e1 = [1, 0, . . . , 0]⊺ . 我们以一个 5 × 5 的矩阵 A 为例. 第一步: 令 Q1 = diag(I1×1, H1), 其中 H1 是对应于向量 A(2 : 5, 1) 的 Householder 矩阵. 于是可得 Q1A = ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ ∗