《矩阵计算》课程教学资源（讲义）第四讲 非对称特征值问题

第四讲非对称特征值问题设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又[;/1<1,i=2,3,...,n,所以limi=2.3.....n故当k趋向于无穷大时，向量（1.%(2)k = 0,1,2,..-收敛到e1=[1,0,,0]T.所以向量(k)=Ak(0) /Ak(0)]l2收敛到±u1,即A的对应于(模)最大的特征值入1的特征向量.而μk=（Ar(k),r())则收敛到Au1=>1.显然，幕送代的收敛快慢取决于[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)