第九讲特征值问题的迭代解法当矩阵规模很大时，计算其所有的特征值和特征向量是非常困难的。而在实际应用中，我们通常也只对其中的某些特征值和（或）特征向量感兴趣，因此也没有必要计算所有的特征值和特征向量本讲主要介绍计算部分特征值和（或）特征值向量的迭代解法.这些算法的存储量要远小于O（n2)运算量也远小于O(n3)关于大规模矩阵特征值问题的参考资料有：· B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice Hall, Englewood Cliffs, NJ, 1980. (Republishedby SIAM, Philadelphia, 1998.)Z. Bai, J.Demmel, J.Dongarra, A. Ruhe, and H. van der Vorst, Templates for the Solution of AlgebraicEigenvalue Problems: APractical Guide, SIAM, Philadelphia, 2000.Y. Saad, Numerical Methods for Large Eigenvalue Problems, 2nd revised edition.SIAM, Philadelphia, 2011.. G. H. Golub and Ch. van Loan, Matrix Computations, 4th edition, Johns Hopkins University Press, Balti-more,2013..Z. Bai, R. C. Li, and Y. F. Su, Lecture notes on Matrix Eigenvalue Computations, 2009.课程：Prof.Bai,LargeScaleScientificComputing,2013.http://www.cs.ucdavis.edu/~bai/EcS231/·课程:Prof.Dr.Peter Arbenz,Numerical Methods for Solving Large Scale Eigenvalue Problems, 2014.http://people.inf.ethz.ch/arbenz/ewp/·软件:ARPACK-collection of Fortran subroutines designed to compute a few eigenvalues and eigenvectorsof largescalesparsematricesandpencilshttp://www.caam.rice.edu/software/ARPACK/9.1投影算法最简单的特征值问题就是仅仅计算一个特征值，如计算模最大的特征值.这时我们可以使用幂选代方法，这在前面章节已经介绍.为了后面讨论方便，我们仍将该算法描述出来算法9.1.暴选代：计算最大特征值1: Given 2(0)2: fori =1,2,., until converge do3:y(i) = Ar(i1)() = y() / ly()24:5: end for249
第九讲 特征值问题的迭代解法 当矩阵规模很大时, 计算其所有的特征值和特征向量是非常困难的. 而在实际应用中, 我们通常也 只对其中的某些特征值和 (或) 特征向量感兴趣, 因此也没有必要计算所有的特征值和特征向量. 本讲主要介绍计算部分特征值和 (或) 特征值向量的迭代解法. 这些算法的存储量要远小于 O(n 2 ), 运算量也远小于 O(n 3 ). 关于大规模矩阵特征值问题的参考资料有: • B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice Hall, Englewood Cliffs, NJ, 1980. (Republished by SIAM, Philadelphia, 1998.) • Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide, SIAM, Philadelphia, 2000. • Y. Saad, Numerical Methods for Large Eigenvalue Problems, 2nd revised edition. SIAM, Philadelphia, 2011. • G. H. Golub and Ch. van Loan, Matrix Computations, 4th edition, Johns Hopkins University Press, Baltimore, 2013. • Z. Bai, R. C. Li, and Y. F. Su, Lecture notes on Matrix Eigenvalue Computations, 2009. • 课程: Prof. Bai, Large Scale Scientific Computing, 2013. http://www.cs.ucdavis.edu/~bai/ECS231/ • 课程: Prof. Dr. Peter Arbenz, Numerical Methods for Solving Large Scale Eigenvalue Problems, 2014. http://people.inf.ethz.ch/arbenz/ewp/ • 软件: ARPACK – collection of Fortran subroutines designed to compute a few eigenvalues and eigenvectors of large scale sparse matrices and pencils http://www.caam.rice.edu/software/ARPACK/ 9.1 投影算法 最简单的特征值问题就是仅仅计算一个特征值, 如计算模最大的特征值. 这时我们可以使用幂迭代 方法, 这在前面章节已经介绍. 为了后面讨论方便, 我们仍将该算法描述出来. 算法 9.1. 幂迭代: 计算最大特征值 1: Given x (0) 2: for i = 1, 2, . . ., until converge do 3: y (i) = Ax(i−1) 4: x (i) = y (i)/∥y (i)∥2 5: end for 249
第九讲特征值问题的选代解法.250.幂选代所产生的选代向量r(0),a(1)....,(m-1)生成一个Krylov子空间[α(), Ar()..Am-1()/ = span [(), (1),.., (m-1)Km(A, r(0)) = span 在幂选代中，我们取z(m-1）为近似特征向量.显然，如果我们在Km(A,(0)）中找出“最佳"的近似特征向量，则收敛速度就可能会大大加快下面我们讨论如何在Km=Km（A,a(0)）中寻找“最佳”的近似特征向量.设AERnxn，并设Km和Cm是R"的两个m维子空间.投影算法就是寻找A的近似特征对（入，?),满足下面的Petrov-Galerkin条件findieCandiEKmsuchthatAi-iilCm.(9.1)这样的算法我们称为斜投影算法.如果我们取Cm=Km，则上面的算法就是一个正交投影算法，此时条件(9.1)称为Galerkin条件设vo,1,...,Um-1和wo,w1,...，Wm-1分别是Km和Cm的一组标准正交基令Vm=[o,U1.Um-il,Wm=[wo,W1,...,wm-i].则对任意EKm，存在向量yERn使得=Vmy，即可以由Vo，U1....,Um-1线性表出.根据条件(9.1)，我们可得(AVmy-iVmy,w)=0,i=1,2,...,m,即WIAVmy=AWIVmy.(9.2)这是一个广义特征值问题.如果我们取Lm=Km，并令Wm=Vm，则（9.2）就转化为Tmy=iy(9.3)其中Tm=VTAVmERmxm，这意味着（>，y）是矩阵Tm的一个特征对.由于m通常比较小，因此我们可以使用先前讨论的方法（如QR选代）来计算（入，y）.这样我们就可以计算出A的一个近似特征对(a,a),其中=Vmy.9.2Rayleigh-Ritz方法事实上，我们可以在Km（A，r(0)）中找出m个最佳近似特征向量及相应的最佳近似特征值.这些近似特征值和近似特征向量就称为Ritz值和Ritz向量定义9.1设Km是Rnxn的一个m维子空间，它的一组标准正交基为vo,u1,...，Um-1，并令Vm=[uo,U1....，Um-1]记Tm=VmAVm，设（,y)是Tm的一组特征对，即Tmy=入y且lyll2=1.则称是A的一个Ritz值，=Vmy是A的一个Ritz向量Rayleigh-Ritz方法就是用Ritz值和Ritz向量来近似A的特征值与特征向量
· 250 · 第九讲 特征值问题的迭代解法 幂迭代所产生的迭代向量 x (0) , x (1), . . . , x(m−1) 生成一个 Krylov 子空间 Km(A, x(0)) = span { x (0), Ax(0), . . . , Am−1x (0)} = span { x (0), x(1), . . . , x(m−1)} . 在幂迭代中, 我们取 x (m−1) 为近似特征向量. 显然, 如果我们在 Km(A, x(0)) 中找出 “最佳” 的近似特征 向量, 则收敛速度就可能会大大加快. 下面我们讨论如何在 Km = Km(A, x(0)) 中寻找 “最佳” 的近似特征向量. 设 A ∈ R n×n, 并设 Km 和 Lm 是 R n 的两个 m 维子空间. 投影算法就是寻找 A 的近似特征对 (λ, ˜ x˜), 满足下面的 Petrov-Galerkin 条件 find λ˜ ∈ C and x˜ ∈ Km such that Ax˜ − λ˜x˜ ⊥ Lm. (9.1) 这样的算法我们称为斜投影算法. 如果我们取 Lm = Km, 则上面的算法就是一个正交投影算法, 此时条 件 (9.1) 称为 Galerkin 条件. 设 v0, v1, . . . , vm−1 和 w0, w1, . . . , wm−1 分别是 Km 和 Lm 的一组标准正交基. 令 Vm = [v0, v1, . . . , vm−1], Wm = [w0, w1, . . . , wm−1]. 则对任意 x˜ ∈ Km, 存在向量 y ∈ R n 使得 x˜ = Vmy, 即 x˜ 可以由 v0, v1, . . . , vm−1 线性表出. 根据条件 (9.1), 我们可得 (AVmy − λV˜ my, wi) = 0, i = 1, 2, . . . , m, 即 W⊺ mAVmy = λW˜ ⊺ mVmy. (9.2) 这是一个广义特征值问题. 如果我们取 Lm = Km, 并令 Wm = Vm, 则 (9.2) 就转化为 Tmy = λy, ˜ (9.3) 其中 Tm = V ⊺ mAVm ∈ R m×m. 这意味着 (λ, y ˜ ) 是矩阵 Tm 的一个特征对. 由于 m 通常比较小, 因此我 们可以使用先前讨论的方法 (如 QR 迭代) 来计算 (λ, y ˜ ). 这样我们就可以计算出 A 的一个近似特征对 (λ, ˜ x˜), 其中 x˜ = Vmy. 9.2 Rayleigh-Ritz 方法 事实上, 我们可以在 Km(A, x(0)) 中找出 m 个最佳近似特征向量及相应的最佳近似特征值. 这些近 似特征值和近似特征向量就称为 Ritz 值 和 Ritz 向量. 定义 9.1 设 Km 是 R n×n 的一个 m 维子空间, 它的一组标准正交基为 v0, v1, . . . , vm−1, 并令 Vm = [v0, v1, . . . , vm−1]. 记 Tm = V ⊺ mAVm, 设 (λ, y ˜ ) 是 Tm 的一组特征对, 即 Tmy = λy˜ 且 ∥y∥2 = 1. 则称 λ˜ 是 A 的一个 Ritz 值, x˜ = Vmy 是 A 的一个 Ritz 向量. Rayleigh-Ritz 方法 就是用 Ritz 值和 Ritz 向量来近似 A 的特征值与特征向量
9.2Rayleigh-Ritz方法.251 算法9.2.RayleighRitz算法1：计算Km的一组单位正交基：Vm=[uo,U1....,Um-1]2:计算矩阵Tm=VTAVm的特征值（即A的Ritz值）3:选取其中的k个想要的特征值：入1，入2....,入，其中k≤m4:计算对应的特征向量：y1,y2.·,yk，即Tmyi=入iy5:计算A的Ritz向量：,=Vmi,i=1,2,...,k.下面我们讨论A是对称矩阵时，RayleighRitz方法的性质设VuERnxn-m是一个列正交矩阵，且使得V=[Vm，Vu]ERnxn是一个正交矩阵.于是我们有[VIAVmVIAVuVTAV = [Vm, Vu]-A[Vm, Vu] -VIAVmVIAV由于A对称的，故Tm=VIAVmEIRm×m也是对称的.此时，Ritz值和Ritz向量具有下面的最优性质定理9.1设AERnxn对称，则对任意的对称矩阵REIRmxm，有AVm-VmR2≥[|AVm-VmTml12,即lAVm-VmR2在R=Tm处取最小值，此时llAVm-VmRll2=1VIAVml12证明.令R=Tm+Z,其中ZERmxm是任意一个对称矩阵.由于A和Tm都是对称矩阵，且Tm=VmAVm，因此由IAVm-VmR=(AVm-VmR)(AVm-VmR)= (AVm - Vm(Tm + Z))(AVm - Vm(Tm + Z))=lAVm-VmTml-VmAVmZ+TmVVmZ-ZVmAVm+ZVIVmTm+ZiVIVmZ=AVm-VmTm+Z所以当Z=0时AVm-VmR达到最小，且AVm- VmTm2=1VVTAVm-VmTmll2[VAVmVmT.[(Vm, Vu)VTAVm= IIVu(VIAVm)l2=IIVIAVmll2 口
9.2 Rayleigh-Ritz 方法 · 251 · 算法 9.2. Rayleigh Ritz 算法 1: 计算 Km 的一组单位正交基: Vm = [v0, v1, . . . , vm−1]. 2: 计算矩阵 Tm = V ⊺ mAVm 的特征值 (即 A 的 Ritz 值) 3: 选取其中的 k 个想要的特征值: λ˜ 1, λ˜ 2, . . . , λ˜ k, 其中 k ≤ m. 4: 计算对应的特征向量: y1, y2, . . . , yk, 即 Tmyi = λ˜ iyi 5: 计算 A 的 Ritz 向量: x˜i = Vmyi , i = 1, 2, . . . , k. 下面我们讨论 A 是对称矩阵时, Rayleigh Ritz 方法的性质. 设 Vu ∈ R n×n−m 是一个列正交矩阵, 且使得 V = [Vm, Vu] ∈ R n×n 是一个正交矩阵. 于是我们有 V ⊺AV = [Vm, Vu] ⊺A[Vm, Vu] = [ V ⊺ mAVm V ⊺ mAVu V ⊺ u AVm V ⊺ u AVu ] . 由于 A 对称的, 故 Tm = V ⊺ mAVm ∈ R m×m 也是对称的. 此时, Ritz 值和 Ritz 向量具有下面的最优性质. 定理 9.1 设 A ∈ R n×n 对称, 则对任意的对称矩阵 R ∈ R m×m, 有 ∥AVm − VmR∥2 ≥ ∥AVm − VmTm∥2, 即 ∥AVm − VmR∥2 在 R = Tm 处取最小值, 此时 ∥AVm − VmR∥2 = ∥V ⊺ u AVm∥2. 证明. 令 R = Tm + Z, 其中 Z ∈ R m×m 是任意一个对称矩阵. 由于 A 和 Tm 都是对称矩阵, 且 Tm = V ⊺ mAVm, 因此由 ∥AVm − VmR∥ 2 2 = (AVm − VmR) ⊺ (AVm − VmR) = ( AVm − Vm(Tm + Z) )⊺( AVm − Vm(Tm + Z) ) = ∥AVm − VmTm∥ 2 2 − V ⊺ mAVmZ + TmV ⊺ mVmZ − ZV ⊺ mAVm + ZV ⊺ mVmTm + Z ⊺V ⊺ mVmZ = ∥AVm − VmTm∥ 2 2 + ∥Z∥ 2 2 . 所以当 Z = 0 时 ∥AVm − VmR∥ 2 2 达到最小, 且 ∥AVm − VmTm∥2 = ∥V V ⊺AVm − VmTm∥2 = (Vm, Vu) [ V ⊺ mAVm V ⊺ u AVm ] − VmTm 2 = ∥Vu(V ⊺ u AVm)∥2 = ∥V ⊺ u AVm∥2 . □
第九讲特征值问题的选代解法.252.+定理9.1中的2-范数可以改成任意的酉不变范数，如F-范数.定理9.2设AERnxn对称，并设Tm=UAUT是Tm=VIAVm的特征值分解.设QERnxm是满足span(Q)=K的任意单位列正交矩阵，DERmxm是任意对角矩阵，我们有AQ-QD/2≥IAVm-VmTml12,且当Q=VmU,D=A时等式成立证明。由于span(Q)=K=span(Vm)，所以存在矩阵WERmxm使得Q=VmW.又Q是单位列正交的,因此I = Q'Q = (VmW)'VmW = WIVIVmW =- WIW.这表明W是正交矩阵.设WDWT=Tm+Z,则AQ-QD=AVmW-VmWD2=AVm-VmWDWTJ2=AVm-VmTm2+Z≥AVm-VmTml2口如果取W=U和D=A则Z=WDWT-Tm=UAUT-Tm=0.此时上式中的等号成立定理9.2表明，在所有满足span(Q)=K单位列正交矩阵QERnxm和任意的对角矩阵DERm中，当Q=VmU和D=A时，AQ-QD2取到极小值定理9.1和定理9.2表明，在IAQ－QD2极小的意义下，Ritz值是特征值的“最佳”近似.所以我们用Ritz值作为特征值的近似是有道理的9.3Lanczos方法设AERnxn是对称矩阵.Lanczos方法就是利用Lanczos过程来计算Km的基和Tm=VIAVm，然后计算A的Ritz值和Ritz向量算法9.3.Lanczos算法1: Choose a vector vo such that voll = 1, and set βBo = 02:for j=0,1,...doCompute w= Aug -Bjuj-13:4:Qj+1 = (w, u,)5:w=w-j+1uj6:βj+1 = [|27:if Bj+1=0then8:stop
· 252 · 第九讲 特征值问题的迭代解法 † 定理 9.1 中的 2-范数可以改成任意的酉不变范数, 如 F-范数. 定理 9.2 设 A ∈ R n×n 对称, 并设 Tm = UΛU ⊺ 是 Tm = V ⊺ mAVm 的特征值分解. 设 Q ∈ R n×m 是满 足 span(Q) = K 的任意单位列正交矩阵, D ∈ R m×m 是任意对角矩阵. 我们有 ∥AQ − QD∥2 ≥ ∥AVm − VmTm∥2, 且当 Q = VmU, D = Λ 时等式成立. 证明. 由于 span(Q) = K = span(Vm), 所以存在矩阵 W ∈ R m×m 使得 Q = VmW. 又 Q 是单位列正交 的, 因此 I = Q ⊺Q = (VmW) ⊺VmW = W⊺V ⊺ mVmW = W⊺W. 这表明 W 是正交矩阵. 设 W DW⊺ = Tm + Z, 则 ∥AQ − QD∥ 2 2 = ∥AVmW − VmW D∥ 2 2 = ∥AVm − VmW DW⊺ ∥ 2 2 = ∥AVm − VmTm∥ 2 2 + ∥Z∥ 2 2 ≥ ∥AVm − VmTm∥ 2 2 . 如果取 W = U 和 D = Λ, 则 Z = W DW⊺ − Tm = UΛU ⊺ − Tm = 0. 此时上式中的等号成立. □ 定理 9.2 表明, 在所有满足 span(Q) = K 单位列正交矩阵 Q ∈ R n×m 和任意的对角矩阵 D ∈ R m 中, 当 Q = VmU 和 D = Λ 时, ∥AQ − QD∥2 取到极小值. 定理 9.1 和定理 9.2 表明, 在 ∥AQ − QD∥2 极小的意义下, Ritz 值是特征值的 “最佳” 近似. 所以我 们用 Ritz 值作为特征值的近似是有道理的. 9.3 Lanczos 方法 设 A ∈ R n×n 是对称矩阵. Lanczos 方法就是利用 Lanczos 过程来计算 Km 的基和 Tm = V ⊺ mAVm, 然 后计算 A 的 Ritz 值和 Ritz 向量. 算法 9.3. Lanczos 算法 1: Choose a vector v0 such that ∥v0∥ = 1, and set β0 = 0 2: for j = 0, 1, . . . do 3: Compute w = Avj − βjvj−1 4: αj+1 = (w, vj ) 5: w = w − αj+1vj 6: βj+1 = ∥w∥2 7: if βj+1 = 0 then 8: stop
9.3Lanczos方法·253-endif9:10:Ui+1 = w/Bi+111:Compute the eigenvalues and eigenvectors of T12:Checktheconvergence13: end for在Lanczos算法9.3中，选代m步后，向量uo.ut.....Um-1构成子空间Km(A, vo) = span fvo, Auo,..., Am-1 o) ,的一组基,并且有AVm=VmTm+BmUmem,其中 em = [0,0,...,0,1]T E Rmβ1[a1B1Tm = VIAVm -βm-βm-1am设(1,9)是 Tm的一个特征对,则有A(Vmy)=VmTmy+BmUmemy=iVmy+Bm(emy)um于是Ar-Ar=Bm(emy)um即Ar-Xrll2=Bm(em)l其中a=Vmy.如果Iβmemy)|很小，则我们就认为是A的某个特征值的很好的近似.事实上，关于Ritz值入，我们有下面的性质引理9.1设AERnxn是对称矩阵，设r=Ar-入，其中≠0.则eia, - lI1:12E(A)其中α（A）表示A的谱，即所有特征值组成的集合证明设A=UAUT是矩阵A的特征值分解.则T=(A-)T=U(A-)UT.记A=diag(入1,2,..,入n).对任意向量z=[21,z2,..,zn]TER",有I(-) =(, -)22-1
9.3 Lanczos 方法 · 253 · 9: end if 10: vj+1 = w/βj+1 11: Compute the eigenvalues and eigenvectors of Tj 12: Check the convergence 13: end for 在 Lanczos 算法 9.3 中, 迭代 m 步后, 向量 v0, v1, . . . , vm−1 构成子空间 Km(A, v0) = span { v0, Av0, . . . , Am−1 v0 } , 的一组基, 并且有 AVm = VmTm + βmvme ⊺ m, 其中 em = [0, 0, . . . , 0, 1]⊺ ∈ R m, Tm = V ⊺ mAVm = α1 β1 β1 . . . . . . . . . . . . βm−1 βm−1 αm . 设 (λ, y ˜ ) 是 Tm 的一个特征对, 则有 A(Vmy) = VmTmy + βmvme ⊺ my = λV˜ my + βm(e ⊺ my)vm. 于是 Ax − λx˜ = βm(e ⊺ my)vm, 即 ∥Ax − λx˜ ∥2 = |βm(e ⊺ my)|, 其中 x = Vmy. 如果 |βm(e ⊺ my)| 很小, 则我们就认为 λ˜ 是 A 的某个特征值的很好的近似. 事实上, 关于 Ritz 值 λ˜, 我们有下面的性质. 引理 9.1 设 A ∈ R n×n 是对称矩阵, 设 r = Ax − λx˜ , 其中 x ̸= 0. 则 min λ∈σ(A) |λ − λ˜| ≤ ∥r∥2 ∥x∥2 , 其中 σ(A) 表示 A 的谱, 即所有特征值组成的集合. 证明. 设 A = UΛU ⊺ 是矩阵 A 的特征值分解. 则 r = (A − λI˜ )x = U(Λ − λI˜ )U ⊺x. 记 Λ = diag(λ1, λ2, . . . , λn). 对任意向量 z = [z1, z2, . . . , zn] ⊺ ∈ R n, 有 ∥(Λ − λI˜ )z∥ 2 2 = ∑n i=1 (λi − λ˜) 2 z 2 i
.254.第九讲特征值问题的选代解法min.-i12zAEO(A)i=1= /2min、->2XEGA所以l12=[UTrl2=1(-)UTl12≥UT2,min-=[l2min-26口由引理9.1可知，存在A的某个特征值入，使得I- 1 Bm(erw)mll _ Bml enl=[Bml-[emyl-(9.4)1Vmyl2lyll2十在前面的讨论中，我们并没有考虑实际计算时可能的舍人误差，在实际计算中，由于浮点运算的舍人误差，即使m很小（如m=10或m=20)，也可能会导致向量（u失去正交性这时我们必须采取一些补救措施，最简单的方法就是对它们重新来一次正交化，即在算法9.3的第5步后加上一条语句(w,wi)viW=这个过程就称为带全正交过程的Lanczos算法.显然，这个过程是非常费时的另外一个可行的方法就是选择性正交9.4Arnoldi方法这里考虑非对称情形，即计算非对称矩阵A的特征值.与Lanczos算法相类似，我们可以通过Arnoldi过程来计算Km的标准正交基vo，U1...,Um-1和上Hessenberg矩阵Hm=VIAVm，使得AVm=VmHm+hm+1.mUmemandVmAVm=Hm.但此时Hm只是上Hessenberg，而不是对称三对角.但我们同样可以通过计算Hm的特征值和特征向量来得到A的Ritz值的Ritz向量，并用它们来近似A的特征值和特征向量设（,y)是Hm的一个特征对，其中yll2=1,则A(Vmy)=VmHmy+hm+1mUmemy=iVmy+hm+1,m(emy)um所以Ar-rl2=[/hm+1.m(emy)uml2=[hm+1,ml·lemyl其中=Vmy.若|hm+1,ml·lemyl足够小,我就认为(,a)是A的某个特征对的近似
· 254 · 第九讲 特征值问题的迭代解法 ≥ ∑n i=1 min λ∈σ(A) |λ − λ˜| 2 z 2 i = ∥z∥ 2 2 min λ∈σ(A) |λ − λ˜| 2 . 所以 ∥r∥2 = ∥U ⊺ r∥2 = ∥(Λ − λI˜ )U ⊺x∥2 ≥ ∥U ⊺x∥2 min λ∈σ(A) |λ − λ˜| = ∥x∥2 · min λ∈σ(A) |λ − λ˜|. □ 由引理 9.1 可知, 存在 A 的某个特征值 λ, 使得 |λ − λ˜| ≤ ∥βm(e ⊺ my)vm∥2 ∥Vmy∥2 = |βm| · |e ⊺ my| ∥y∥2 = |βm| · |e ⊺ my|. (9.4) † 在前面的讨论中, 我们并没有考虑实际计算时可能的舍入误差. 在实际计算中, 由于浮点运算的 舍入误差, 即使 m 很小 (如 m = 10 或 m = 20), 也可能会导致向量 {vi} 失去正交性. 这时我们必 须采取一些补救措施, 最简单的方法就是对它们重新来一次正交化, 即在算法 9.3 的第 5 步后加 上一条语句 w = w − ∑ j i=1 (w, vi)vi . 这个过程就称为带全正交过程的 Lanczos 算法. 显然, 这个过程是非常费时的. 另外一个可行的方法就是选择性正交. 9.4 Arnoldi 方法 这里考虑非对称情形, 即计算非对称矩阵 A 的特征值. 与 Lanczos 算法相类似, 我们可以通过 Arnoldi 过程来计算 Km 的标准正交基 v0, v1, . . . , vm−1 和上 Hessenberg 矩阵 Hm = V ⊺ mAVm, 使得 AVm = VmHm + hm+1,mvme ⊺ m and V ⊺ mAVm = Hm. 但此时 Hm 只是上 Hessenberg, 而不是对称三对角. 但我们同样可以通过计算 Hm 的特征值和特征向量 来得到 A 的 Ritz 值的 Ritz 向量, 并用它们来近似 A 的特征值和特征向量. 设 (λ, y ˜ ) 是 Hm 的一个特征对, 其中 ∥y∥2 = 1, 则 A(Vmy) = VmHmy + hm+1,mvme ⊺ my = λV˜ my + hm+1,m(e ⊺ my)vm, 所以 ∥Ax − λx˜ ∥2 = ∥hm+1,m(e ⊺ my)vm∥2 = |hm+1,m| · |e ⊺ my|, 其中 x = Vmy. 若 |hm+1,m| · |e ⊺ my| 足够小, 我就认为 (λ, x ˜ ) 是 A 的某个特征对的近似
9.5非对称Lanczos方法.255.算法9.4.Arnoldi算法1: Choose a vector Uo such that lvoll2 = 12: for j = 0, 1....do3:Compute wj+1= Auj4:fori=o,l,...jdo5:hij = (wj+1,i)6:Wi+1=Wj+1-hijVi7:end for8:hj+1,3 = wj+1l29:if hj+1j = 0 then10:stopend if11:12:Uj+1=Wj+1/hj+1,j13:Compute the eigenvalues and eigenvectors of T14:Checkthe convergence15:end for由于A是非对称的,其特征值可能是复的，或者是坏条件的,此时Lanczos方法的一些最优性质就不再成立.尽快如此，目前还是存在一些有效Arnoldi方法的实现方式，可参见[50,63,65]9.5非对称Lanczos方法非对称Lanczos方法就是Lanczos方法在非对称矩阵上的推广，它是基于Lanczos双正交化过程设o和wo是任意的非零向量，并设Km(A, vo) = span(vo, Avo,..,Am-1vo)和Km(AT,wo)= spanfwo,ATwo.....(AT)m-1wo)Lanczos双正交化过程就是计算Km（A,o）和Km（A,wo）的基【u}和（w;]，满足【u}和[w}相互正交即ij(vi,wi)l, i=j所以WIVm = I,其中Vm=[vo,Ui....,Um--],Wm=[wo,Wi,...,Wm-]注意，通常【ui=。或w]-。本身并不一定正交，故Vm和Wm通常并不列正交
9.5 非对称 Lanczos 方法 · 255 · 算法 9.4. Arnoldi 算法 1: Choose a vector v0 such that ∥v0∥2 = 1 2: for j = 0, 1, . . . do 3: Compute wj+1 = Avj 4: for i = 0, 1, . . . , j do 5: hij = (wj+1, vi) 6: wj+1 = wj+1 − hijvi 7: end for 8: hj+1,j = ∥wj+1∥2 9: if hj+1,j = 0 then 10: stop 11: end if 12: vj+1 = wj+1/hj+1,j 13: Compute the eigenvalues and eigenvectors of Tj 14: Check the convergence 15: end for 由于 A 是非对称的, 其特征值可能是复的, 或者是坏条件的, 此时 Lanczos 方法的一些最优性质就不 再成立. 尽快如此, 目前还是存在一些有效 Arnoldi 方法的实现方式, 可参见 [50, 63, 65]. 9.5 非对称 Lanczos 方法 非对称 Lanczos 方法就是 Lanczos 方法在非对称矩阵上的推广, 它是基于 Lanczos 双正交化过程. 设 v0 和 w0 是任意的非零向量, 并设 Km(A, v0) = span{v0, Av0, . . . , Am−1 v0} 和 Km(A ⊺ , w0) = span{w0, A⊺w0, . . . ,(A ⊺ ) m−1w0}. Lanczos 双正交化过程就是计算 Km(A, v0) 和 Km(A, w0) 的基 {vi} 和 {wi}, 满足 {vi} 和 {wi} 相互正 交, 即 (vi , wj ) = 0, i ̸= j 1, i = j . 所以 W⊺ mVm = I, 其中 Vm = [v0, v1, . . . , vm−1], Wm = [w0, w1, . . . , wm−1]. 注意, 通常 {vi} m i=0 或 {wj} m j=0 本身并不一定正交, 故 Vm 和 Wm 通常并不列正交
.256.第九讲特征值问题的选代解法根据Lanczos双正交化过程，我们可得AVm=VmTm+mUmemATWm=WmTI+BmWmem其中em =[0,...,0,1]T e Rm,β1[a]Y1Tm =Bm-m-1S所以WIAVm=WIVmTm+m(WIUm)em=Tm设入是Tm的特征值,其对应的右特征向量和左特征向量分别为y和z,且 Il3ll2=zll2=1,即Tmy=Xy and 21Tm=AzT.于是有AVmy=VmTmy+mUmemy=iVmy+melyum,(Wmz)TA=(ATWmz)T=(WmTz)T+Bmemzwm=i(Wmz)T+Bmemzwm若hm(emg)|和[Bm(emz)|足够小，我们就认为入是A的某个特征值的近似，而Vmy和Wmz就是相应的右特征向量和左特征向量的近似算法9.5.非对称Lanczos算法1:Choose two vectors voand wosuchthat (vo,wo)=12: Set Bo = 0 and % = 03: for j=0,1,...do4:ComputeQj+i=(Auj,w)5:Uj+1=Avj-aj+1wj-Biuj-16:ij+1= ATwj -αj+1wj -wj-1+1=[(j+1,j+1)/27:8:if j+1=0 then9:stopendif10:11:Bj+1=(oj+1,ij+1)/+112:Uj+1 = 0j+1/%i+113:Wj+1 = wj+1/Bi+114:Computetheeigenvalues and eigenvectorsof T
· 256 · 第九讲 特征值问题的迭代解法 根据 Lanczos 双正交化过程, 我们可得 AVm = VmTm + γmvme ⊺ m, A ⊺Wm = WmT ⊺ m + βmwme ⊺ m, 其中 em = [0, . . . , 0, 1]⊺ ∈ R m, Tm = α1 β1 γ1 . . . . . . . . . . . . βm−1 γm−1 αm . 所以 W⊺ mAVm = W⊺ mVmTm + γm(W⊺ mvm)e ⊺ m = Tm. 设 λ˜ 是 Tm 的特征值, 其对应的右特征向量和左特征向量分别为 y 和 z, 且 ∥y∥2 = ∥z∥2 = 1, 即 Tmy = λy˜ and z ⊺Tm = λz˜ ⊺ . 于是有 AVmy = VmTmy + γmvme ⊺ my = λV˜ my + γme ⊺ myvm, (Wmz) ⊺A = (A ⊺Wmz) ⊺ = (WmT ⊺ mz) ⊺ + βme ⊺ mzw⊺ m = λ˜(Wmz) ⊺ + βme ⊺ mzw⊺ m. 若 |γm(e ⊺ my)| 和 |βm(e ⊺ mz)| 足够小, 我们就认为 λ˜ 是 A 的某个特征值的近似, 而 Vmy 和 Wmz 就是相应 的右特征向量和左特征向量的近似. 算法 9.5. 非对称 Lanczos 算法 1: Choose two vectors v0 and w0 such that (v0, w0) = 1 2: Set β0 = 0 and γ0 = 0 3: for j = 0, 1, . . . do 4: Compute αj+1 = (Avj , wj ) 5: v˜j+1 = Avj − αj+1vj − βjvj−1 6: w˜j+1 = A⊺wj − αj+1wj − γjwj−1 7: γj+1 = |(˜vj+1, w˜j+1)| 1/2 8: if γj+1 = 0 then 9: stop 10: end if 11: βj+1 = (˜vj+1, w˜j+1)/γj+1 12: vj+1 = ˜vj+1/γj+1 13: wj+1 = ˜wj+1/βj+1 14: Compute the eigenvalues and eigenvectors of Tj
9.5非对称Lanczos方法·257.15:Checktheconvergence16:endfor非对称Lanczos算法的显著优点就是节省运算量，缺点是更容易被中断
9.5 非对称 Lanczos 方法 · 257 · 15: Check the convergence 16: end for 非对称 Lanczos 算法的显著优点就是节省运算量, 缺点是更容易被中断