第四讲非对称特征值问题设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又[/入1<1,i=2,3,...,n,所以lin0i=2.3.....n.故当k趋向于无穷大时，向量[.%()（)k=0,1,2...收敛到e1=[1,0...,0].所以向量(K)=A(0) /A*(0)2收敛到±1,即A的对应于（模)最大的特征值入的特征向量.而=((k))*Az(）则收敛到Au1=入1.显然，幂选代的收敛快慢取决于[>2/入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<1,所以当k趋于无穷大时，W(K),趋向于0.令V=[Vp,Vn-pl,则W(k)=(vpw(k) + Vn-pw())VA*V-1Zo =>,[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| < 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 ) . 所以当 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 标准型, 我们可以得到同样的结论, 见 [73, 74]. 在正交迭代中, 如果我们取 Z0 = I, 则可得到一类特殊的正交迭代算法. 此时, 在一定条件下, 正交 迭代会收敛到 A 的 Schur 标准型. 4.4 QR 迭代方法 QR 迭代算法的基本思想是通过不断的正交变换, 将 A 转化为上三角形式 (或拟上三角形式). 算法 形式非常简单, 描述如下: