第五讲对称特征值问题设AERnxn是对称矩阵.在计算A的特征值和特征值向量时，我们可以充分利用A的对称结构，一方面尽可能地减少运算量，另一方面也能构造出更加快速高效的算法，本讲主要介绍以下方法·Jacobi选代法：较古老的方法，收敛速度较慢，但能达到很高的计算精度，且非常适合并行·Rayleigh商选代法：利用Rayleigh商作为动态位移的反选代算法，一般具有全局线性收敛和局部三次收敛性·对称QR选代法：针对对称矩阵的带位移隐式QR送代算法如果只需计算一个对称三对角矩阵的所有特征值，则该算法是目前最快的方法，运算量为（n2).如果需要计算所有的特征值和特征向量，则运算量约为6n3分而治之法：同时计算对称矩阵的特征值和特征向量的一种快速算法.基本思想是将大矩阵分解成小矩阵，然后利用递归思想求解，是目前求解对称矩阵的所有特征值和特征向量的最快方法之一·对分法和反选代法：对分法主要用于求解对称矩阵在某个区间中的特征值，反送代法用于计算特征向量除了Jacobi选代和Rayleigh商代外,其余算法都需要先将对称矩阵三对角化.这个过程大约需花费n3的工作量，如果需要计算特征向量的话，则总运算量约为5.1Jacobi送代法该算法的基本思想是通过一系列的Jacobi旋转Js，将A正交相似于一个对角矩阵，即A(O) = A,A(h+1)= JkA(k)JT,k=0,1,...)且A()收敛到一个对角矩阵，其中Js为Jacobi旋转，通常选取J为Givens变换，即1cos OksinOkJh = G(ik, Jk, Ou) =cosOk sin Ok1157
第五讲 对称特征值问题 设 A ∈ R n×n 是对称矩阵. 在计算 A 的特征值和特征值向量时, 我们可以充分利用 A 的对称 结构, 一方面尽可能地减少运算量, 另一方面也能构造出更加快速高效的算法. 本讲主要介绍以下方法. • Jacobi 迭代法: 较古老的方法, 收敛速度较慢, 但能达到很高的计算精度, 且非常适合并行. • Rayleigh 商迭代法: 利用 Rayleigh 商作为动态位移的反迭代算法, 一般具有全局线性收敛和 局部三次收敛性. • 对称 QR 迭代法: 针对对称矩阵的带位移隐式 QR 迭代算法. 如果只需计算一个对称三对角 矩阵的所有特征值, 则该算法是目前最快的方法, 运算量为 O(n 2 ). 如果需要计算所有的特征 值和特征向量, 则运算量约为 6n 3 . • 分而治之法: 同时计算对称矩阵的特征值和特征向量的一种快速算法. 基本思想是将大矩阵 分解成小矩阵, 然后利用递归思想求解, 是目前求解对称矩阵的所有特征值和特征向量的最 快方法之一. • 对分法和反迭代法: 对分法主要用于求解对称矩阵在某个区间中的特征值, 反迭代法用于计 算特征向量. b 除了 Jacobi 迭代和 Rayleigh 商迭代外, 其余算法都需要先将对称矩阵三对角化. 这个过程大 约需花费 4 3 n 3 的工作量, 如果需要计算特征向量的话, 则总运算量约为 8 3 n 3 . 5.1 Jacobi 迭代法 该算法的基本思想是通过一系列的 Jacobi 旋转 Jk, 将 A 正交相似于一个对角矩阵, 即 A (0) = A, A(k+1) = JkA (k)J T k , k = 0, 1, . . . , 且 A(k) 收敛到一个对角矩阵, 其中 Jk 为 Jacobi 旋转, 通常选取 Jk 为 Givens 变换, 即 Jk = G(ik, jk, θk) = 1 . . . 1 cos θk sin θk . . . − sin θk cos θk 1 . . . 1 157
-158第五讲对称特征值问题易知,在A()两边分别左乘J和右乘T时，只会修改A(K)的第i和第行,以及第i和第j列.由于A()是对称矩阵，由下面的引理可知,通过选取适当的k,可以将A()(i,jk)和A()(jki）同时化为0.引理5.1设AER2x2是对称矩阵，则存在Givens变换GER2×2,使得GAGT为对角矩阵(板书)证明.设cosOsinesindhFOS则COsasineCos 6sinGAGTsin0sinocosOcosO/acos?+csin?+bsin20(c-a)sin 20 + b cos 20asin20+ccos9-bsin20(c-a)sin 20+bcos20令（c-a)sin20+bcos20=0,可得1 - tan? 0a-c=cot20262tan解得sign(←)tang:26I/ + V1++2口故引理结论成立为了使得A()收敛到一个对角矩阵，其非对角线元素必须趋向于0.记of(A)为所有非对角线元素的平方和,即Aaoff(A) =α, = IIAl=1讲我们的目标就是使得off(A）尽快趋于0
· 158 · 第五讲 对称特征值问题 易知, 在 A(k) 两边分别左乘 Jk 和右乘 J T k 时, 只会修改 A(k) 的第 ik 和第 jk 行, 以及第 ik 和第 jk 列. 由于 A(k) 是对称矩阵, 由下面的引理可知, 通过选取适当的 θk, 可以将 A(k) (ik, jk) 和 A(k) (jk, ik) 同时化为 0. 引理 5.1 设 A ∈ R 2×2 是对称矩阵, 则存在 Givens 变换 G ∈ R 2×2 , 使得 GAGT 为对角矩阵. (板书) 证明. 设 A = " a b b c# , G = " cos θ sin θ − sin θ cos θ # , 则 GAGT = " cos θ sin θ − sin θ cos θ # " a b b c# " cos θ sin θ − sin θ cos θ #T = " a cos2 θ + c sin2 θ + b sin 2θ 1 2 (c − a)sin 2θ + b cos 2θ 1 2 (c − a)sin 2θ + b cos 2θ a sin2 θ + c cos2 θ − b sin 2θ # 令 1 2 (c − a)sin 2θ + b cos 2θ = 0, 可得 a − c 2b = cot 2θ = 1 − tan2 θ 2 tan θ . 解得 tan θ = sign(τ ) |τ | + √ 1 + τ 2 , τ = a − c 2b . 故引理结论成立. □ 为了使得 A(k) 收敛到一个对角矩阵, 其非对角线元素必须趋向于 0. 记 off(A) 为所有非对角 线元素的平方和, 即 off(A) = X i̸=j a 2 ij = ∥A∥ 2 F − Xn i=1 a 2 ii , 我们的目标就是使得 off(A) 尽快趋于 0
5.1Jacobi送代法·159.引理5.2设A=[ai] eRnxn是对称矩阵,A=[ai]=JAJT,J=G(i,j,0),其中 a的选取使得ai=aj=0,则off(A)=off(A)-2az，i≠j.(板书)证明.记A=[a1,a2,...,an].令A=JA=[ai]nxn.由于J是正交阵，故IJakll2 = akl|2, k = 1,2,..., n.又J左乘ak时，只影响其第i和第个元素的值，故由Jai2=aill2和Jai2=aill2可得最+=喵+哦，喝+哆=喝+喝(5.1)同理，由A=AJT可得+=+，+=+(5.2)又a=a=0故+=++喝+=喵++2由于JAJT只影响A的第ii行和第i,列,故对角线元素中只有au和aii受影响.所以Zaix =ai + 2dig,=1k=1故= AI-off(A) = IIA Zak-2a,= off(A)- 2aj-k=1k=1口即引理结论成立.由此可知，off(A()）总是不断减小的.下面给出Jacobi选代算法算法5.1.Jacobi选代算法1: Given a symmetric matrix Ae IRnxn2: if eigenvectors are desired then3:set J = I and flag = 14:endif5: while not converge dochoose an index pair (i, i) such that aij + 06:7:T = (a-ajs)/(2ai)t = sign()/(I| + V1+T2)% 计算 tan g8:c=1/V1+t,s=ct%计算cosの和sine9:A = G(i, j, 0)AG(i,j,0)T10:11:if flag = 1 then12:J=G(i,j,0)J
5.1 Jacobi 迭代法 · 159 · 引理 5.2 设 A = [aij ] ∈ R n×n 是对称矩阵, Aˆ = [ˆaij ] = JAJT, J = G(i, j, θ), 其中 θ 的选取使 得 aˆij = ˆaji = 0, 则 off(Aˆ) = off(A) − 2a 2 ij , i ̸= j. (板书) 证明. 记 A = [a1, a2, . . . , an]. 令 A˜ = JA = [˜aij ]n×n. 由于 J 是正交阵, 故 ∥Jak∥2 = ∥ak∥2, k = 1, 2, . . . , n. 又 J 左乘 ak 时, 只影响其第 i 和第 j 个元素的值, 故由 ∥Jai∥2 = ∥ai∥2 和 ∥Jaj∥2 = ∥aj∥2 可得 a˜ 2 ii + ˜a 2 ji = a 2 ii + a 2 ji, a˜ 2 ij + ˜a 2 jj = a 2 ij + a 2 jj . (5.1) 同理, 由 Aˆ = AJ˜ T 可得 aˆ 2 ii + ˆa 2 ij = ˜a 2 ii + ˜a 2 ij , aˆ 2 ji + ˆa 2 jj = ˜a 2 ji + ˜a 2 jj . (5.2) 又 aˆij = ˆaji = 0, 故 aˆ 2 ii + ˆa 2 jj = a 2 ii + a 2 jj + a 2 ij + a 2 ji = a 2 ii + a 2 jj + 2a 2 ij . 由于 JAJT 只影响 A 的第 i, j 行和第 i, j 列, 故对角线元素中只有 aii 和 ajj 受影响. 所以 Xn k=1 aˆ 2 kk = Xn k=1 a 2 kk + 2a 2 ij , 故 off(Aˆ) = ∥Aˆ∥ 2 2 − Xn k=1 aˆ 2 kk = ∥A∥ 2 2 − Xn k=1 a 2 kk − 2a 2 ij = off(A) − 2a 2 ij , 即引理结论成立. □ 由此可知, off(A(k) ) 总是不断减小的. 下面给出 Jacobi 迭代算法. 算法 5.1. Jacobi 迭代算法 1: Given a symmetric matrix A ∈ R n×n 2: if eigenvectors are desired then 3: set J = I and flag = 1 4: end if 5: while not converge do 6: choose an index pair (i, j) such that aij ̸= 0 7: τ = (aii − ajj )/(2aij ) 8: t = sign(τ )/(|τ | + √ 1 + τ 2) % 计算 tan θ 9: c = 1/√ 1 + t 2, s = ct % 计算 cos θ 和 sin θ 10: A = G(i, j, θ)AG(i, j, θ) T 11: if flag = 1 then 12: J = G(i, j, θ)J
-160.第五讲对称特征值问题end if13:14:end while该算法涉及到ai的选取问题，一种直观的选取方法就是使得ai为所有非对角线元素中绝对值最大的一个，这就是经典Jacobi选代算法算法5.2.经典Jacobi选代算法1: Given a symmetric matrix Ae IRnxn2: if eigenvectors are desired thenset J = I and flag = 13:4: end if5: while off(A)> tol do6:choose (i, j) such that [ail = maxk# [al%选取绝对值最大的元素7:T= (an -aji)/(2ai)8:t = sign(↑)/(I| + V1+ ↑2)9:c=1/V1+t2,s=ct%计算cos0和sine10:A = G(i,j, 0)AG(i, j, 0)T11:if flag =1 then12:J=G(i,j,O)J13:end if14: end while可以证明，经典Jacobi算法至少是线性收敛的定理5.3对于经典Jacobi算法5.2,有)of(A(), N= n(n -1)of(A(t+1)≤ (1- 1)故k步选代后，有of(4()≤ (1-) of(A()=(1-) (4),(板书)证明.由于在经典Jacobi算法 5.2中,[ail=maxk+ilaxil,故 off(A()≤n(n+1)(a)),即2()≥f(4(), N=n(n_1),2所以由引理5.2可知of(A(+1) =of(A() - ()~≤(1-) of(A()口事实上，经典Jacobi算法最终是（渐进）二次收敛的[30，100]
· 160 · 第五讲 对称特征值问题 13: end if 14: end while 该算法涉及到 aij 的选取问题, 一种直观的选取方法就是使得 aij 为所有非对角线元素中绝 对值最大的一个, 这就是经典 Jacobi 迭代算法. 算法 5.2. 经典 Jacobi 迭代算法 1: Given a symmetric matrix A ∈ R n×n 2: if eigenvectors are desired then 3: set J = I and flag = 1 4: end if 5: while off(A)> tol do 6: choose (i, j) such that |aij | = maxk̸=l |akl| % 选取绝对值最大的元素 7: τ = (aii − ajj )/(2aij ) 8: t = sign(τ )/(|τ | + √ 1 + τ 2) 9: c = 1/√ 1 + t 2, s = ct % 计算 cos θ 和 sin θ 10: A = G(i, j, θ)AG(i, j, θ) T 11: if flag = 1 then 12: J = G(i, j, θ)J 13: end if 14: end while 可以证明, 经典 Jacobi 算法至少是线性收敛的. 定理 5.3 对于经典 Jacobi 算法 5.2, 有 off(A (k+1)) ≤ 1 − 1 N off(A (k) ), N = n(n − 1) 2 . 故 k 步迭代后, 有 off(A (k) ) ≤ 1 − 1 N k off(A (0)) = 1 − 1 N k off(A). (板书) 证明. 由于在经典 Jacobi 算法 5.2 中, |aij | = maxk̸=l |akl|, 故 off(A(k) ) ≤ n(n + 1) a (k) ij 2 , 即 2 a (k) ij 2 ≥ 1 N off(A (k) ), N = n(n − 1) 2 . 所以由引理 5.2 可知 off(A (k+1)) = off(A (k) ) − a (k) ij 2 ≤ 1 − 1 N off(A (k) ). □ 事实上, 经典 Jacobi 算法最终是 (渐进) 二次收敛的 [30, 100]
5.1Jacobi选代法+ 161 .定理5.4经典Jacobi算法5.2是N步（渐进)二次收敛的，即对足够大的k，有off(A(k+N) = O(ff(A())由于在经典Jacobi算法中，每一步都要寻找绝对值最大的非对角元，比较费时，因此实用性较差.我们可以通过逐行扫描来选取（i，i），这就是循环Jacobi选代算法算法5.3.循环Jacobi选代算法（逐行扫描）1:Given a symmetric matrixAeIRnxn2: if eigenvectors are desired then3:set J = I and flag = 14:end if5: while off(A)> tol do6:for i = 1 to n -1 do7:for j - i+ 1 to n doif aij≠0then9:T= (ai -aj)/(2ag)10:t = sign()/(l/ + V1 +T2)c= 1/V1+ t211:12:s=c-t13:A = G(i,j, 0)TAG(i,j, 0)14:if flag = 1 then15:J=J.G(i,j,0)16:end ifendif17:18:endfor19:end for20:endwhile循环Jacobi也具有（渐进)二次收敛性[135,page270].Jacobi选代法的优缺点·优点：能够达到很高的计算精度（特别是小特征值）；同时非常适合并行计算·缺点：计算速度较慢；矩阵稀疏性得不到充分的利用
5.1 Jacobi 迭代法 · 161 · 定理 5.4 经典 Jacobi 算法 5.2 是 N 步 (渐进) 二次收敛的, 即对足够大的 k, 有 off A (k+N) = O off2 A (k) . 由于在经典 Jacobi 算法中, 每一步都要寻找绝对值最大的非对角元, 比较费时, 因此实用性较 差. 我们可以通过逐行扫描来选取 (i, j), 这就是循环 Jacobi 迭代算法. 算法 5.3. 循环 Jacobi 迭代算法 (逐行扫描) 1: Given a symmetric matrix A ∈ R n×n 2: if eigenvectors are desired then 3: set J = I and flag = 1 4: end if 5: while off(A)> tol do 6: for i = 1 to n − 1 do 7: for j = i + 1 to n do 8: if aij ̸= 0 then 9: τ = (aii − ajj )/(2aij ) 10: t = sign(τ )/(|τ | + √ 1 + τ 2) 11: c = 1/√ 1 + t 2 12: s = c · t 13: A = G(i, j, θ) TAG(i, j, θ) 14: if flag = 1 then 15: J = J · G(i, j, θ) 16: end if 17: end if 18: end for 19: end for 20: end while 循环 Jacobi 也具有 (渐进) 二次收敛性 [135, page 270]. Jacobi 迭代法的优缺点 • 优点: 能够达到很高的计算精度 (特别是小特征值); 同时非常适合并行计算. • 缺点: 计算速度较慢; 矩阵稀疏性得不到充分的利用
-162第五讲对称特征值问题5.2Rayleigh商送代法在反送代算法中，a()的Rayleigh商是特征值入的近似，因此我们可以把它作为位移，于是就得到下面的Rayleigh商选代法算法5.4.Rayleigh商选代算法(RQI,RayleighQuotientIterations)1: Given an initial guess r(0) with [r(0)Il2 = 12: compute the Rayleigh quotient Po = (r(0), Ar(0)3: set k = 14: while not converge do5:g = Pk-1g(k) = (A -αI)-12(k-1)6:() = y(k) / lg() 27:8:Pk = (r(k), Ar(k)9:k=k+110:end while关于Rayleigh商送代的收敛性，我们有下面的结论定理5.5如果特征值是单重的，则当误差足够小时,Rayleigh商选代法中每步选代所得的正确数(板书)字的位数增至三倍，即Rayleigh商选代是局部三次收敛的。证明.设A=QAQT,令全()=QTa(k),则在Rayleigh商送代算法中Pk = (2(k)TAr(k) = (全(k)TQTAQ(k) = ((k)TA(k)令g(k)= QTy(K),则g(k) = QT(A - Pk-11)-1(k) = (QTAQ - Pk-1I)-1(k-1) = (A - Pk-1I)-1(k-1),即，“以初始向量r(0)对A做Rayleigh商选代"等价于“以初始向量(0)对A做Rayleigh商选代即它们有相同的收敛性.因此,不失一般性，我们可以假定A=A为对角阵，此时A的特征向量为ei,i= 1,2,...,n.我们假定r(k)收敛到e1.令dk=r(k)-e1,则dkll2→0.为了证明算法具有局部三次收敛，我们需要证明：当ek=[ld:/2充分小时,有 ek+1= dk+1ll2=I2(k+1)-e1ll2=O(ex)。我们注意到1 = (r(k))Ta(k) = (e1 + de)T(e1 + dk) = 1 + 2de(1) + dfdk = 1 +2dk(1)+e其中dk(1)表示ds的第一个元素.故di(1)=-/2.所以Pk = (r()TA(k) = (e1 + de)TA(ei + dk)=eTAe1 + 2eTAd +dAd, = 1 -n,其中=-(2eTAdk+dk)=-2d(1)-Adk=-dk于是≤ + A12 d2A2
· 162 · 第五讲 对称特征值问题 5.2 Rayleigh 商迭代法 在反迭代算法中, x (k) 的 Rayleigh 商是特征值 λk 的近似, 因此我们可以把它作为位移, 于是就 得到下面的 Rayleigh 商迭代法. 算法 5.4. Rayleigh 商迭代算法 (RQI, Rayleigh Quotient Iterations) 1: Given an initial guess x (0) with ∥x (0)∥2 = 1 2: compute the Rayleigh quotient ρ0 = (x (0), Ax(0)) 3: set k = 1 4: while not converge do 5: σ = ρk−1 6: y (k) = (A − σI) −1x (k−1) 7: x (k) = y (k)/∥y (k)∥2 8: ρk = (x (k) , Ax(k) ) 9: k = k + 1 10: end while 关于 Rayleigh 商迭代的收敛性, 我们有下面的结论. 定理 5.5 如果特征值是单重的, 则当误差足够小时, Rayleigh 商迭代法中每步迭代所得的正确数 字的位数增至三倍, 即 Rayleigh 商迭代是局部三次收敛的. (板书) 证明. 设 A = QΛQT, 令 xˆ (k) = QTx (k) , 则在 Rayleigh 商迭代算法中 ρk = (x (k) ) TAx(k) = (ˆx (k) ) TQ TAQxˆ (k) = (ˆx (k) ) TΛˆx (k) . 令 yˆ (k) = QTy (k) , 则 yˆ (k) = Q T (A − ρk−1I) −1x (k) = (Q TAQ − ρk−1I) −1xˆ (k−1) = (Λ − ρk−1I) −1xˆ (k−1) , 即, “以初始向量 x (0) 对 A 做 Rayleigh 商迭代” 等价于 “以初始向量 xˆ (0) 对 Λ 做 Rayleigh 商迭代”, 即它们有相同的收敛性. 因此, 不失一般性, 我们可以假定 A = Λ 为对角阵, 此时 A 的特征向量为 ei , i = 1, 2, . . . , n. 我们假定 x (k) 收敛到 e1. 令 dk = x (k) − e1, 则 ∥dk∥2 → 0. 为了证明算法具有局部三次收敛, 我们需要证明: 当 εk = ∥dk∥2 充分小时, 有 εk+1 = ∥dk+1∥2 = ∥x (k+1) − e1∥2 = O(ε 3 k ). 我们注意到 1 = (x (k) ) Tx (k) = (e1 + dk) T (e1 + dk) = 1 + 2dk(1) + d T k dk = 1 + 2dk(1) + ε 2 k , 其中 dk(1) 表示 dk 的第一个元素. 故 dk(1) = −ε 2 k /2. 所以 ρk = (x (k) ) TΛx (k) = (e1 + dk) TΛ(e1 + dk) = e T kΛe1 + 2e T 1 Λdk + d T kΛdk ≜ λ1 − η, 其中 η = −(2e T 1 Λdk + d T kΛdk) = −2λ1dk(1) − d T kΛdk = λ1ε 2 k − d T kΛdk. 于是 |η| ≤ |λ1|ε 2 k + ∥Λ∥2 · ∥dk∥ 2 2 ≤ 2∥Λ∥2 ε 2 k
5.2Rayleigh商选代法·163.由Rayleigh商算法5.4可知g(k+1) = (A - pμI)-1()[()(1) 2(k)(2)(k)(n)入1 - pk2 - pkAn-pkd;(n) Td(2)[1 + dk(1)入1 - pkX2-PkAn-pk[1 -/2dk(2)dk(n)入2-1 + n入n-+nn1 -%/2d;(2)nd:(n)n/2)(2—入1 +n)n (1 -/2)(入n -1 +n)α 1-e/2 (e1 + dx+1),n其中dk(2)ndk(n)ndk+1(1-% /2)(入2 -入1 +n)(1-e/2)(n-1+n)因为入1是单重特征值，所以gap(>1, A) min [ ->1l >0,故对于i=2.3...,n，当ek足够小时有[ - + l ≥ [; -] -[l ≥ gap(1, ) -[l ≥ gap(1, ) - 2]l2 > 0.于是我们有2/A2 Ild;l/2 nldk+1l1(1/2)(gap(A1,))I(1 -e /2)(gap(A1,A) -nl)即dx+1,=O(e). 又1 x+1l, ≤le1 + dx+1, ≤1 + dk+1l即[1-e1 + dx+1=d+1由于g(k+1)ei + dk+2(k+1)Il(k+1)12e1 + dk+ll所以(1 le1 + dk+1/) e + d+IIdk+1l2 = I/2(k+1) - ell/ei+dk+1[1Ile1 + d+/ + d+ll2dk+1ei+dk+1e1 + dk+1口又dk+1=0(ex),故ek+1=1dk+12=0(e元)
5.2 Rayleigh 商迭代法 · 163 · 由 Rayleigh 商算法 5.4 可知 y (k+1) = (Λ − ρkI) −1x (k) = " x (k) (1) λ1 − ρk , x (k) (2) λ2 − ρk , . . . , x (k) (n) λn − ρk #T = 1 + dk(1) λ1 − ρk , dk(2) λ2 − ρk , . . . , dk(n) λn − ρk T = 1 − ε 2 k /2 η , dk(2) λ2 − λ1 + η , . . . , dk(n) λn − λ1 + η T = 1 − ε 2 k /2 η 1, dk(2)η (1 − ε 2 k /2)(λ2 − λ1 + η) , . . . , dk(n)η (1 − ε 2 k /2)(λn − λ1 + η) T ≜ 1 − ε 2 k /2 η · (e1 + ˆdk+1). 其中 ˆdk+1 = 0, dk(2)η (1 − ε 2 k /2)(λ2 − λ1 + η) , . . . , dk(n)η (1 − ε 2 k /2)(λn − λ1 + η) T . 因为 λ1 是单重特征值, 所以 gap(λ1,Λ) ≜ min i̸=1 |λi − λ1| > 0, 故对于 i = 2, 3, . . . , n, 当 εk 足够小时有 |λi − λ1 + η| ≥ |λi − λ1| − |η| ≥ gap(λ1,Λ) − |η| ≥ gap(λ1,Λ) − 2∥Λ∥2 ε 2 k > 0. 于是我们有 ˆdk+1 2 ≤ ∥dk∥2 |η| (1 − ε 2 k /2)(gap(λ1,Λ) − |η|) ≤ 2∥Λ∥2 ε 3 k (1 − ε 2 k /2)(gap(λ1,Λ) − |η|) , 即 ˆdk+1 2 = O(ε 3 k ). 又 1 − ˆdk+1 2 ≤ e1 + ˆdk+1 2 ≤ 1 + ˆdk+1 2 , 即 1 − e1 + ˆdk+1 2 ≤ ˆdk+1 2 . 由于 x (k+1) = y (k+1) ∥y (k+1)∥2 = e1 + ˆdk+1 e1 + ˆdk+1 2 , 所以 ∥dk+1∥2 = ∥x (k+1) − e1∥2 = 1 − e1 + ˆdk+1 2 e1 + ˆdk+1 2 e1 + ˆdk+1 2 ≤ 1 − e1 + ˆdk+1 2 + ˆdk+1 2 e1 + ˆdk+1 2 ≤ 2 ˆdk+1 2 e1 + ˆdk+1 2 . 又 ˆdk+1 2 = O(ε 3 k ), 故 εk+1 = ∥dk+1∥2 = O(ε 3 k ). □
164第五讲对称特征值问题否RQI算法具有局部三次收敛性，但无法确定收敛到哪个特征向量(特征值),因此可以作为其他算法的加速手段，即先使用其他算法（比如幂送代）计算出所需特征值的近似值，然后再使用RQI算法加速.白实际计算时，判断（pk，r(k)）是否收敛可以观察残量rk=（A-PkI)r(k）是否趋于零下面是关于RQI算法的全局收敛性，可参见文献[100]定理5.6在RQI算法中,设rk=(A-PkI)(K),则有rk+1l≤lrk:l其中等号成立当且仅当pk+1=Pk且2(K)是(A一pkI)的特征向量
· 164 · 第五讲 对称特征值问题 b RQI 算法具有局部三次收敛性, 但无法确定收敛到哪个特征向量 (特征值), 因此可以作为其 他算法的加速手段, 即先使用其他算法 (比如幂迭代) 计算出所需特征值的近似值, 然后再 使用 RQI 算法加速. b 实际计算时, 判断 (ρk, x(k) ) 是否收敛可以观察残量 rk = (A − ρkI)x (k) 是否趋于零. 下面是关于 RQI 算法的全局收敛性, 可参见文献 [100]. 定理 5.6 在 RQI 算法中, 设 rk = (A − ρkI)x (k) , 则有 ∥rk+1∥ ≤ ∥rk∥, 其中等号成立当且仅当 ρk+1 = ρk 且 x (k) 是 (A − ρkI) 2 的特征向量
·165.5.3对称QR选代法5.3对称QR选代法将带位移的隐式OR方法运用到对称矩阵，就得到对称OR迭代法.由于此时A是对称的，所以上Hessenberg化后就转化为一个对称三对角矩阵，相应的过程就称为对称三对角化任何一个对称矩阵AERnxn都可以通过正交变换转化成一个对称三对角矩阵T.这个过程可以通过Householder变换或Givens变换来实现对称QR选代算法基本步骤1.对称三对角化：利用Householder变换，将A化为对称三对角矩阵，即寻找正交矩阵Q使得T=QAQT为对称三对角矩阵；2.使用带（单）位移的隐式QR迭代算法计算T的特征值与特征值向量；3.计算A的特征向量.对称QR选代算法的运算量(1）三对角化需要4n3/3+O(n2)），如果需要计算特征向量，即形成矩阵Q，则需要额外4n3/3+O(n2)，因此总运算量为8n3/3+O(n2)；也可以暂时不计算Q而是将所有Householder向量存放在A的严格上三角部分（A对称，因此A只需保存下三角部分）(2)对T做带位移的隐式QR选代，每次选代的运算量为6n;(3）计算T的特征值时，假定每个平均送代2步，则计算所有特征值的运算量为12n2；（每求出一个特征值后，通过defation，矩阵规模也会相应减小，因此运算量可能会更少）(4)若要计算T的所有特征值和特征向量,则运算量为6n3+O(n2);(5）若只要计算A的所有特征值，运算量为4n3/3+O（n2);(6）若需要计算A的所有的特征值和特征向量，则运算量为26n3/3+O(n2);位移的选取一Wilkinson位移位移的选取直接影响到算法的收敛速度.我们可以通过下面的方式来选取位移.设6(k)[a(k)6(k).Ak=6(k)b(-1ak)一种简单的位移选取策略就是令k=an).这种位移选取方法几乎对所有的矩阵都有三次渐进收敛速度，但也存在不收敛的例子，故我们需要对其做一些改进图事实上，a(k）就是收敛到特征向量的选代向量的Rayleigh商：如果在QR选代算法中使用位移gk=am,而在Rayleigh商选代算法5.4中取o=en=[0..,0,1]T,则利用QR送代与
5.3 对称 QR 迭代法 · 165 · 5.3 对称 QR 迭代法 将带位移的隐式 QR 方法运用到对称矩阵, 就得到对称 QR 迭代法. 由于此时 A 是对称的, 所 以上 Hessenberg 化后就转化为一个对称三对角矩阵, 相应的过程就称为 对称三对角化. b 任何一个对称矩阵 A ∈ R n×n 都可以通过正交变换转化成一个对称三对角矩阵 T. 这个过 程可以通过 Householder 变换或 Givens 变换来实现. 对称 QR 迭代算法基本步骤 1. 对称三对角化: 利用 Householder 变换, 将 A 化为对称三对角矩阵, 即寻找正交矩阵 Q 使得 T = QAQT 为对称三对角矩阵; 2. 使用带 (单) 位移的隐式 QR 迭代算法计算 T 的特征值与特征值向量; 3. 计算 A 的特征向量. 对称 QR 迭代算法的运算量 (1) 三对角化需要 4n 3/3 + O(n 2 ), 如果需要计算特征向量, 即形成矩阵 Q, 则需要额外 4n 3/3 + O(n 2 ), 因此总运算量为 8n 3/3 + O(n 2 ); 也可以暂时不计算 Q, 而是将所有 Householder 向 量存放在 A 的严格上三角部分 (A 对称, 因此 A 只需保存下三角部分). (2) 对 T 做带位移的隐式 QR 迭代, 每次迭代的运算量为 6n; (3) 计算 T 的特征值时, 假定每个平均迭代 2 步, 则计算所有特征值的运算量为 12n 2 ; (每求出 一个特征值后, 通过 deflation, 矩阵规模也会相应减小, 因此运算量可能会更少) (4) 若要计算 T 的所有特征值和特征向量, 则运算量为 6n 3 + O(n 2 ); (5) 若只要计算 A 的所有特征值, 运算量为 4n 3/3 + O(n 2 ); (6) 若需要计算 A 的所有的特征值和特征向量, 则运算量为 26n 3/3 + O(n 2 ); 位移的选取 — Wilkinson 位移 位移的选取直接影响到算法的收敛速度. 我们可以通过下面的方式来选取位移. 设 Ak = a (k) 1 b (k) 1 b (k) 1 . . . . . . . . . . . . b (k) n−1 b (k) n−1 a (k) n , 一种简单的位移选取策略就是令 σk = a (k) n . 这种位移选取方法几乎对所有的矩阵都有三次渐进 收敛速度, 但也存在不收敛的例子, 故我们需要对其做一些改进. b 事实上, a (k) n 就是收敛到特征向量的迭代向量的 Rayleigh 商: 如果在 QR 迭代算法中使用位 移 σk = a (k) n , 而在 Rayleigh 商迭代算法 5.4 中取 x0 = en = [0, . . . , 0, 1]T, 则利用 QR 迭代与
-166.第五讲对称特征值问题反送代之间的关系可以证明：QR选代算法中的k与Rayleigh商选代算法中的pk相等[a(-]]6(k)的最接近a的特征值作为一种有效的位移是Wilkinson位移，即取子矩阵[62, ]位移.通过计算可得Wilkinson位移为=()+8-sign()82+(62),其中 8=(a(21-a%),出于稳定性方面的考虑，我们通常用下面的计算公式(6(k) )g = a(k)8 + sign(8) /82 + (6(2.)定理5.7[57,100]采用Wilkinson位移的QR选代是整体收敛的，且至少是线性收敛.事实上几乎对所有的对称矩阵都是渐进三次收敛的例5.1带Wilkinson位移的隐式QR选代算法收敛性演示(Eig_TriQR.m)T1.9345e-01000-1.1495e+00001.9345e-01-5.7144e-01-3.5163e+0000-3.5163e+001.4138e+00-1.2639e+0000-1.2639e+00-2.0125e-014.3216e+000004.3216e+001.9285e+00iter = 1A人00-1.1606e+002.0488e-0100-1.0240e+002.0488e-01-3.1370e+0000-3.5560e+00-1.0240e+001.2439e+0000-3.5560e+00-1.1053e+003.1938e-010003.1938e-015.5790e+00iter = 2T =000-1.1748e+002.6621e-01002.6621e-01-3.3005e+00-6.4187e-0100-6.4187e-01-3.1162e+00-1.8413e+000e-1.8413e+003.4052e+001.3658e-030001.3658e-035.6064e+00iter = 3
· 166 · 第五讲 对称特征值问题 反迭代之间的关系可以证明: QR 迭代算法中的 σk 与 Rayleigh 商迭代算法中的 ρk 相等. 一种有效的位移是 Wilkinson 位移, 即取子矩阵 " a (k) n−1 b (k) n−1 b (k) n−1 a (k) n # 的最接近 a (k) n 的特征值作为 位移. 通过计算可得 Wilkinson 位移为 σ = a (k) n + δ − sign(δ) r δ 2 + b (k) n−1 2 , 其中 δ = 1 2 (a (k) n−1 − a (k) n ). 出于稳定性方面的考虑, 我们通常用下面的计算公式 σ = a (k) n − b (k) n−1 2 δ + sign(δ) r δ 2 + b (k) n−1 2 . 定理 5.7 [57, 100] 采用 Wilkinson 位移的 QR 迭代是整体收敛的, 且至少是线性收敛. 事实上, 几乎对所有的对称矩阵都是渐进三次收敛的. 例 5.1 带 Wilkinson 位移的隐式 QR 迭代算法收敛性演示. (Eig_TriQR.m) T = ‐1.1495e+00 1.9345e‐01 0 0 0 1.9345e‐01 ‐5.7144e‐01 ‐3.5163e+00 0 0 0 ‐3.5163e+00 1.4138e+00 ‐1.2639e+00 0 0 0 ‐1.2639e+00 ‐2.0125e‐01 4.3216e+00 0 0 0 4.3216e+00 1.9285e+00 iter = 1 T = ‐1.1606e+00 2.0488e‐01 0 0 0 2.0488e‐01 ‐3.1370e+00 ‐1.0240e+00 0 0 0 ‐1.0240e+00 1.2439e+00 ‐3.5560e+00 0 0 0 ‐3.5560e+00 ‐1.1053e+00 3.1938e‐01 0 0 0 3.1938e‐01 5.5790e+00 iter = 2 T = ‐1.1748e+00 2.6621e‐01 0 0 0 2.6621e‐01 ‐3.3005e+00 ‐6.4187e‐01 0 0 0 ‐6.4187e‐01 ‐3.1162e+00 ‐1.8413e+00 0 0 0 ‐1.8413e+00 3.4052e+00 1.3658e‐03 0 0 0 1.3658e‐03 5.6064e+00 iter = 3