第五讲对称特征值问题设AERnxn是对称矩阵.在计算A的特征值和特征值向量时，我们可以充分利用A的对称结构，一方面尽可能地减少运算量，另一方面也能构造出更加快速高效的算法关于对称矩阵的特征值和特征值向量，目前常用算法有：·Jacobi送代：最古老的方法，收敛速度较慢，但精度较高，且很适合并行计算·Rayleigh商选代：利用Rayleigh商作为位移的反送代算法，一般具有三次收敛性·对称QR迭代：计算对称矩阵的特征值和特征向量的QR算法.如果只需计算对称三对角矩阵的所有特征值，则该算法是目前最快的方法（运算量为O（n)).如果需要计算所有的特征值和特征向量，则运算量约为6n3.·分而治之法（Divide-and-Conquer）：计算对称三对角矩阵的特征值和特征值向量的一种快速算法基本思想是将大矩阵分解成小矩阵，然后利用递推思想求特征值和特征向量.在最坏的情形下，运算量为O(n3），但在实际应用中，平均为O(n2.3).如果使用快速多极子算法(FMM)后，理论上的运算量可降低到O(nlog’n)，其中p是一个较小的整数，这使得分而治之算法成为目前求解对称三对角矩阵的所有特征值和特征向量的最快方法之一。·对分法和反选代：对分法主要用于求解对称三对角矩阵在某个区间中的特征值，运算量约为O(kn)，其中k为所需计算的特征值的个数.反选代用于计算特征向量，在最佳情况下，即特征值“适当分离”时，运算量约为O(kn），但在最差情况下，即特征值成串地紧靠在一起时，运算量约为O(2n)）而且不能保证特征向量的精度（虽然实际上它几乎是精确的）十除了Jacobi送代和Rayleigh商选代外，其余算法都需要先将对称矩阵三对角化，这个过程大约需花费n的工作量，如果需要计算特征向量的话，则运算量约为n3福5.1Jacobi迭代该算法的基本思想是通过一系列的Jacobi旋转J将A正交相似于二个对角矩阵，即A(O) = A, A(+1) =JIA(K)J,k= 0,1,...且A(）收敛到一个对角矩阵，其中Js为正交矩阵135
第五讲 对称特征值问题 设 A ∈ R n×n 是对称矩阵. 在计算 A 的特征值和特征值向量时, 我们可以充分利用 A 的对称结构, 一方面尽可能地减少运算量, 另一方面也能构造出更加快速高效的算法. 关于对称矩阵的特征值和特征值向量, 目前常用算法有: • Jacobi 迭代: 最古老的方法, 收敛速度较慢, 但精度较高, 且很适合并行计算. • Rayleigh 商迭代: 利用 Rayleigh 商作为位移的反迭代算法, 一般具有三次收敛性. • 对称 QR 迭代: 计算对称矩阵的特征值和特征向量的 QR 算法. 如果只需计算对称三对角矩阵的所 有特征值, 则该算法是目前最快的方法 (运算量为 O(n 2 )). 如果需要计算所有的特征值和特征向量, 则运算量约为 6n 3 . • 分而治之法 (Divide-and-Conquer) : 计算对称三对角矩阵的特征值和特征值向量的一种快速算法. 基本思想是将大矩阵分解成小矩阵, 然后利用递推思想求特征值和特征向量. 在最坏的情形下, 运 算量为 O(n 3 ), 但在实际应用中, 平均为 O(n 2.3 ). 如果使用快速多极子算法 (FMM) 后, 理论上的运 算量可降低到 O(n logp n), 其中 p 是一个较小的整数, 这使得分而治之算法成为目前求解对称三对 角矩阵的所有特征值和特征向量的最快方法之一. • 对分法和反迭代: 对分法主要用于求解对称三对角矩阵在某个区间中的特征值, 运算量约为O(kn), 其中 k 为所需计算的特征值的个数. 反迭代用于计算特征向量, 在最佳情况下, 即特征值 “适当分 离” 时, 运算量约为 O(kn), 但在最差情况下, 即特征值成串地紧靠在一起时, 运算量约为 O(k 2n), 而且不能保证特征向量的精度 (虽然实际上它几乎是精确的). † 除了 Jacobi 迭代和 Rayleigh 商迭代外, 其余算法都需要先将对称矩阵三对角化. 这个过程大约需 花费 4 3 n 3 的工作量, 如果需要计算特征向量的话, 则运算量约为 8 3 n 3 . 5.1 Jacobi 迭代 该算法的基本思想是通过一系列的 Jacobi 旋转 Jk 将 A 正交相似于一个对角矩阵, 即 A (0) = A, A(k+1) = J ⊺ k A (k)Jk, k = 0, 1, . . . , 且 A(k) 收敛到一个对角矩阵, 其中 Jk 为正交矩阵. 135
.136.第五讲对称特征值问题Jacobi旋转J的构造我们通常选取Jk为Givens变换，即iJk1cosOk-sinkihJ&= G(ik,Jk,Ou)=jksinOkcOsOk1由于A()是对称矩阵，所以可以选取适当的，将A()(i,)和A(K)(j,)化为0.引理5.1设AeIR2x2是对称矩阵，则存在Givens变换GER2x2使得GTAG为对角阵证明.设-sinecoseGsinecose则cos6-singcos6-singGTAG=sinecososinecosaacos+csin?0+bsin20(c-a)sin20+bcos20(c-a)sin20+bcos20asin?0+ccos?-bsin20令号（c-a)sin20+bcos20=0即得1- tan?0a-cCot20262tang解得sign(↑)a-ctang:/+ V1++226口故引理结论成立为了使得A()收敛到一个对角矩阵，其非对角元素必须趋向于0.记of(A)为所有非对角元素的平方和，即off(A)=, =IIA -Ea,详i=1我们的目标就是使得of(A）尽快趋向于0
· 136 · 第五讲 对称特征值问题 Jacobi 旋转 Jk 的构造 我们通常选取 Jk 为 Givens 变换, 即 Jk = G(ik, jk, θk) = ik jk 1 . . . 1 ik cos θk − sin θk . . . jk sin θk cos θk 1 . . . 1 . 由于 A(k) 是对称矩阵, 所以可以选取适当的 θk, 将 A(k) (i, j) 和 A(k) (j, i) 化为 0. 引理 5.1 设 A ∈ R 2×2 是对称矩阵, 则存在 Givens 变换 G ∈ R 2×2 使得 G⊺AG 为对角阵. 证明. 设 A = [ a b b c] , G = [ cos θ − sin θ sin θ cos θ ] , 则 G ⊺AG = [ cos θ − sin θ sin θ cos θ ]⊺ [ a b b c] [cos θ − sin θ sin θ cos θ ] = [ 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) = ∑ i̸=j a 2 ij = ∥A∥ 2 F − ∑n i=1 a 2 ii , 我们的目标就是使得 off(A) 尽快趋向于 0
5.1Jacobi选代.137.引理5.2 设 A =[ai]nxn E Rnxn 是对称矩阵, A =[ai]nxn = JTAJ, J = G(i,j,0),其中 的选取使得a=a=0,则ff(A) = of(A) - 2azj证明.设A=[a1,a2,...,an].令A=JTA=[aij]nxn·由于J是正交阵,故JTal2 = al|2, k = 1, 2,..-, n.又JT左乘a时，只影响其第i个元素的值，故由IJTaill2=al2和JTaill2=lall2可得+=++=+(5.1)同理，由 A=AJ可得a+a=a+a+ai=a+ai(5.2)又=a=0,故+aj=a+++dji=+j+2由于JTAJ只影响A的第ij行和第i,j列，故对角线元素中只有aii和aji受影响.所以nLa=)Zakk+2ag,k=1k=1故of(A)=IIAl2-ax=IlAll2 -a%k-2a, = of(A)-2a)k=1k=1口即引理结论成立由此可知,of(A(k)）总是不断减小的.下面给出[acobi迭代算法算法5.1.Jacobi选代算法l:Given a symmetric matrix Ae IRnxn2:if eigenvectors are desired thenset J= I and shift = 13:4: end if5:while not converge do6:choose an index pair (i, j) such that aij 07:T= (au- a)/(2ai)8:t = sign()/(I| + V1+r2)O%计算tan6C=1/V1+t2%计算cosg9:s=c·t%计算sing10:11:A=G(i,j,0)TAG(i,j,0)%实际计算时不需要做矩阵乘积12:if shift =1 then13:J= J·G(i,J,0)
5.1 Jacobi 迭代 · 137 · 引理 5.2 设 A = [aij ]n×n ∈ R n×n 是对称矩阵, Aˆ = [ˆaij ]n×n = J ⊺AJ, J = G(i, j, θ), 其中 θ 的选取 使得 aˆij = ˆaji = 0, 则 off(Aˆ) = off(A) − 2a 2 ij . 证明. 设 A = [a1, a2, . . . , an]. 令 A˜ = J ⊺A = [˜aij ]n×n. 由于 J 是正交阵, 故 ∥J ⊺ ak∥2 = ∥ak∥2, k = 1, 2, . . . , n. 又 J ⊺ 左乘 ak 时, 只影响其第 i, j 个元素的值, 故由 ∥J ⊺ai∥2 = ∥ai∥2 和 ∥J ⊺aj∥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˜ 可得 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 . 由于 J ⊺AJ 只影响 A 的第 i, j 行和第 i, j 列, 故对角线元素中只有 aii 和 ajj 受影响. 所以 ∑n k=1 aˆ 2 kk = ∑n k=1 a 2 kk + 2a 2 ij , 故 off(Aˆ) = ∥Aˆ∥ 2 2 − ∑n k=1 aˆ 2 kk = ∥A∥ 2 2 − ∑n 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 shif t = 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 % 计算 cos θ 10: s = c · t % 计算 sin θ 11: A = G(i, j, θ) ⊺AG(i, j, θ) % 实际计算时不需要做矩阵乘积 12: if shif t = 1 then 13: J = J · G(i, j, θ)
第五讲对称特征值问题.138.endif14:15: end while该算法涉及到a的选取问题，一种直观的选取方法就是使得ai为所有非对角元素中绝对值最大的一个，于是我们就得到下面的经典Jacobi算法算法5.2.经典Jacobi选代算法1:Givena symmetricmatrixAERnxn2: if eigenvectors are desired thensetJ-Iand shift=13:4: endif5:while off(A)> tol do6:choose (i, j) such that aii=maxk±ak7:T= (ai-aj)/(2ai)t = sign(T)/(I/ + V1+ T2)8:c=1/V1+t29:s=c-t10:11:A=G(i,j,0)TAG(i,j,0)12:if shift =1 then13:J= J.G(i,j,0)endif14:15: end while可以证明，经典Jacobi算法至少是线性收敛的定理5.1对于经典Jacobi算法5.2,有om(4(+)≤(1-)of(4(),N=(n=1),-故k步选代后，有of(4()≤(1-)om(4)-(1-) om(4),证明.由于在经典Jacobi算法5.2中,[ail=maxk=iail,故of(A(k)≤n(n+1)(a()）即2(8)≥(4),2所以由引理5.2可知of(4(+1) = f(A() - (a)≤ (1- ) f(A()口事实上，经典Jacobi算法最终是二次局部收敛的
· 138 · 第五讲 对称特征值问题 14: end if 15: 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 shif t = 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 10: s = c · t 11: A = G(i, j, θ) ⊺AG(i, j, θ) 12: if shif t = 1 then 13: J = J · G(i, j, θ) 14: end if 15: end while 可以证明, 经典 Jacobi 算法至少是线性收敛的. 定理 5.1 对于经典 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 算法最终是二次局部收敛的
5.2Rayleigh商迭代.139 .定理5.2经典Jacobi算法5.2是N步局部二次收敛的，即对足够大的k，有off(A(k+N) = 0(off(A())由于在经典Jacobi算法中，每一步都要寻找绝对值最大的非对角元，比较费时，因此它并不实用.我们可以通过逐行扫描来选取（i，），这就是循环Jacobi迭代算法算法5.3.循环Jacobi选代算法（逐行扫描）I:Given a symmetricmatrix ARnxn2: if eigenvectors are desired thenset J=I and shift=13:4: end if5: while off(A)> tol do6:for i=l ton-ldo7:forj=i+ltondo8:ifaij≠0then9:T=(ai-a)/(2ai)t = sign()/(I + V1+ 2)10:c=1/V1+t211:s=ct12:13:A= G(i,j,0)TAG(i,j, 0)14:if shift =1 then15:J=J.G(i,j,0)endif16:end if17:end for18:end for19:20: end while循环Jacobi也具有局部二次收敛性[75,page270]5.2Rayleigh商送代在反迭代算法中，工；的Rayleigh商是特征值入；的近似，因此我们可以把它作为位移，于是就得到下面的Rayleigh商迭代算法算法5.4.Rayleigh商选代I:Givenan initialguessowithloll2=1
5.2 Rayleigh 商迭代 · 139 · 定理 5.2 经典 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 shif t = 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, θ) ⊺AG(i, j, θ) 14: if shif t = 1 then 15: J = J · G(i, j, θ) 16: end if 17: end if 18: end for 19: end for 20: end while 循环 Jacobi 也具有局部二次收敛性 [75, page 270]. 5.2 Rayleigh 商迭代 在反迭代算法中, xi 的 Rayleigh 商是特征值 λi 的近似, 因此我们可以把它作为位移, 于是就得到下 面的 Rayleigh 商迭代算法. 算法 5.4. Rayleigh 商迭代 1: Given an initial guess x0 with ∥x0∥2 = 1
.140.第五讲对称特征值问题2: compute the Rayleigh quotient po =rAro3: set i=14:while not converge do5:o = pi-16:yi =(A-αI)-1i-17:=yi/llyill28:Pi=rTAri9:i=i+l10: end while关于Rayleigh商选代的收敛性，我们有下面的结论定理5.3如果特征值是单重的，则当误差足够小时，Rayleigh商选代中每步选代所得的正确数字的位数增至三倍，即Rayleigh商选代是局部三次收敛的.证明设A=QAQT，令至，=QT工，则在Rayleigh商选代算法中P=A=QTAQ=A令=QTyi，则9i=QT(A-Pi-1I)-1;= (QTAQ-pi-1I)-1;-1= (A-pi-1)-1i-1即，“以初始向量To对A做Rayleigh商迭代”等价于“以初始向量o对A做Rayleigh商选代”，即它们有相同的收敛性.因此，不失一般性，我们可以假定A=A为对角阵，此时A的特征向量为ei，i=1,2,..，n.我们假定；收敛到e1.令d;=i-e1,则dill2→0.为了证明算法具有局部三次收敛，我们需要证明：当=d:2充分小时，有+1=di+12=+1e12=O(e)我们注意到1 = r;ri = (e1 + d;)T(ei + d,) =1 +2d,(1) + d,d, = 1 + 2d,;(1) +e?其中d;(1)表示d;的第一个元素.故d(1)=-ε/2.所以P==(e+d,)(ei+)=ee+2ed, +d,-n其中=-(2eAd,+dAd)=-2d;(1)-dAd=>i=?-dAdg.于是22由Rayleigh商算法5.4可知i+1 = (A - p;I)-1r;;(1);(2);(n) TAn- pi1-i2pi7d;(2)d;(n)1 + d;(1)An - pi]-pi入2 - Pi
· 140 · 第五讲 对称特征值问题 2: compute the Rayleigh quotient ρ0 = x ⊺ 0Ax0 3: set i = 1 4: while not converge do 5: σ = ρi−1 6: yi = (A − σI) −1xi−1 7: xi = yi/∥yi∥2 8: ρi = x ⊺ i Axi 9: i = i + 1 10: end while 关于 Rayleigh 商迭代的收敛性, 我们有下面的结论. 定理 5.3 如果特征值是单重的, 则当误差足够小时, Rayleigh 商迭代中每步迭代所得的正确数字的位 数增至三倍, 即 Rayleigh 商迭代是局部三次收敛的. 证明. 设 A = QΛQ⊺ , 令 xˆi = Q⊺xi , 则在 Rayleigh 商迭代算法中 ρi = x ⊺ i Axi = ˆx ⊺ i Q ⊺AQxˆi = ˆx ⊺ i Λˆxi . 令 yˆi = Q⊺yi , 则 yˆi = Q ⊺ (A − ρi−1I) −1xi = (Q ⊺AQ − ρi−1I) −1xˆi−1 = (Λ − ρi−1I) −1xˆi−1, 即, “以初始向量 x0 对 A 做 Rayleigh 商迭代” 等价于 “以初始向量 xˆ0 对 Λ 做 Rayleigh 商迭代”, 即它们有相 同的收敛性. 因此, 不失一般性, 我们可以假定 A = Λ 为对角阵, 此时 A 的特征向量为 ei , i = 1, 2, . . . , n. 我们假定 xi 收敛到 e1. 令 di = xi − e1, 则 ∥di∥2 → 0. 为了证明算法具有局部三次收敛, 我们需要 证明: 当 εi = ∥di∥2 充分小时, 有 εi+1 = ∥di+1∥2 = ∥xi+1 − e1∥2 = O(ε 3 i ). 我们注意到 1 = x ⊺ i xi = (e1 + di) ⊺ (e1 + di) = 1 + 2di(1) + d ⊺ i di = 1 + 2di(1) + ε 2 i , 其中 di(1) 表示 di 的第一个元素. 故 di(1) = −ε 2 i /2. 所以 ρi = x ⊺ i Λxi = (e1 + di) ⊺Λ(e1 + di) = e ⊺ i Λe1 + 2e ⊺ 1Λdi + d ⊺ i Λdi ≜ λ1 − η, 其中 η = −(2e ⊺ 1Λdi + d ⊺ i Λdi) = −2λ1di(1) − d ⊺ i Λdi = λ1ε 2 i − d ⊺ i Λdi . 于是 |η| ≤ |λ1|ε 2 i + ∥Λ∥2 · ∥di∥ 2 2 ≤ 2∥Λ∥2 ε 2 i . 由 Rayleigh 商算法 5.4 可知 yi+1 = (Λ − ρiI) −1xi = [ xi(1) λ1 − ρi , xi(2) λ2 − ρi , . . . , xi(n) λn − ρi ]⊺ = [ 1 + di(1) λ1 − ρi , di(2) λ2 − ρi , . . . , di(n) λn − ρi ]⊺
5.2Rayleigh商迭代.141 .d;(2)[1-/2d;(n)入21+nn-i+nSE2/2d;(n)nd;(2)n(1 -2/2)(2->1+n)*"(1-=2/2)(入n-入1 +nn1-/22. (e1 + di+1).n其中d;(2)nd;(n)ndi+1(1 -e /2)(入2- 入1 + m)*" (1 -/2)(入n- 入1 +n)因为入1是单重特征值，所以gap(入1, A) min /; - 1/ > 0,故当e，足够小时，[ + ≥/ - ≥gap(1,) -[l ≥gap(1,) -2[2 > 0于是我们有Ild;2 (l2A2[a++l, (1- /2((, ) - /) (1-/2)(ep(a, ) -1 l = 0(e). 又N1+,e1 +d+1++即1e +d+ll/+由于1+d.yi+1Ti+1Ili+1ll2en + di+1l,所以(1-Ile + d+ll)e + d+llId+1l2 =[i+1-e1ll2ei + di+ll1-e + d+ ++ei+di+1口= 0(e),故ei+1=IIdi+1ll2= 0(e)d:关于RQI算法的全局收敛性，可见参考文献[57]
5.2 Rayleigh 商迭代 · 141 · = [ 1 − ε 2 i /2 η , di(2) λ2 − λ1 + η , . . . , di(n) λn − λ1 + η ]⊺ = 1 − ε 2 i /2 η [ 1, di(2)η (1 − ε 2 i /2)(λ2 − λ1 + η) , . . . , di(n)η (1 − ε 2 i /2)(λn − λ1 + η) ]⊺ ≜ 1 − ε 2 i /2 η · (e1 + ˆdi+1). 其中 ˆdi+1 = [ 0, di(2)η (1 − ε 2 i /2)(λ2 − λ1 + η) , . . . , di(n)η (1 − ε 2 i /2)(λn − λ1 + η) ]⊺ . 因为 λ1 是单重特征值, 所以 gap(λ1,Λ) ≜ min i̸=1 |λi − λ1| > 0, 故当 εi 足够小时, |λj − λ1 + η| ≥ |λj − λ1| − |η| ≥ gap(λ1,Λ) − |η| ≥ gap(λ1,Λ) − 2∥Λ∥2 ε 2 i > 0. 于是我们有 ˆdi+1 2 ≤ ∥di∥2 |η| (1 − ε 2 i /2)(gap(λ1,Λ) − |η|) ≤ 2∥Λ∥2 ε 3 i (1 − ε 2 i /2)(gap(λ1,Λ) − |η|) , 即 ˆdi+1 2 = O(ε 3 i ). 又 1 − ˆdi+1 2 ≤ e1 + ˆdi+1 2 ≤ 1 + ˆdi+1 2 , 即 1 − e1 + ˆdi+1 2 ≤ ˆdi+1 2 . 由于 xi+1 = yi+1 ∥yi+1∥2 = e1 + ˆdi+1 e1 + ˆdi+1 2 , 所以 ∥di+1∥2 = ∥xi+1 − e1∥2 = ( 1 − e1 + ˆdi+1 2 ) e1 + ˆdi+1 2 e1 + ˆdi+1 2 ≤ 1 − e1 + ˆdi+1 2 + ˆdi+1 2 e1 + ˆdi+1 2 ≤ 2 ˆdi+1 2 e1 + ˆdi+1 2 . 又 ˆdi+1 2 = O(ε 3 i ), 故 εi+1 = ∥di+1∥2 = O(ε 3 i ). □ 关于 RQI 算法的全局收敛性, 可见参考文献 [57]
.142.第五讲对称特征值问题5.3对称QR迭代将带位移的隐式QR算法运用到对称矩阵，就得到对称QR选代算法，基本步骤1.对称三对角化：利用Householder变换，将A化为对称三对角矩阵，即寻找正交矩阵Q使得T=QAQT为对称三对角矩阵；2.使用带（单）位移的隐式QR选代算法计算T的特征值与特征值向量3.计算A的特征向量对称三对角化任何一个对称矩阵AERnxn都可以通过正交变换转化成一个对称三对角矩阵T.这个过程可以通过Householder变换来实现，也可以通过Givens变换来实现如果A是一个稠密矩阵，则Givens变换的总运算量大约是Householder变换的两倍，因此建议使用Householder变换对称QR送代算法的运算量·三对角化需要4n3/3+O（n2)，如果需要计算特征向量，则运算量为8n3/3+O(n2)·对T做带位移的隐式QR选代，每次迭代的运算量为6n;·计算T的特征值时，假定每个平均选代2步，则计算所有特征值的运算量为12n2；·若要计算T的所有特征值和特征向量，则运算量为6n3十O(n2）·若只要计算A的所有特征值，运算量为4n3/3+O(n2)；·若需要计算A的所有的特征值和特征向量，则运算量为26n3/3+O（n2);位移的选取一Wilkinson位移位移的好坏直接影响到算法的收敛速度.我们可以通过下面的方式来选取位移.设6(k)[a()5'.A(K) =6(4)6(1a(k)一种简单的位移选取策略就是令α=a*)。事实上，ak）就是收敛到特征向量的选代向量的Rayleigh商，这种位移选取方法几乎对所有的矩阵都有三次渐进收敛速度，但也存在不收敛的例子，故我们需要对其做一些改进+如果在QR选代算法中使用位移ak=a，而在Rayleigh商选代算法5.4中取ro=en[0,,0,1]T，则利用QR选代与反选代之间的关系可以证明：QR选代算法中的a与Rayleigh
· 142 · 第五讲 对称特征值问题 5.3 对称 QR 迭代 将带位移的隐式 QR 算法运用到对称矩阵, 就得到对称 QR 迭代算法. 基本步骤 1. 对称三对角化: 利用 Householder 变换, 将 A 化为对称三对角矩阵, 即寻找正交矩阵 Q 使得 T = QAQ⊺ 为对称三对角矩阵; 2. 使用带 (单) 位移的隐式 QR 迭代算法计算 T 的特征值与特征值向量; 3. 计算 A 的特征向量. 对称三对角化 任何一个对称矩阵 A ∈ R n×n 都可以通过正交变换转化成一个对称三对角矩阵 T. 这个过程可以 通过 Householder 变换来实现, 也可以通过 Givens 变换来实现. 如果 A 是一个稠密矩阵, 则 Givens 变换的总运算量大约是 Householder 变换的两倍, 因此建议使用 Householder 变换. 对称 QR 迭代算法的运算量 • 三对角化需要 4n 3/3 + O(n 2 ), 如果需要计算特征向量, 则运算量为 8n 3/3 + O(n 2 ); • 对 T 做带位移的隐式 QR 迭代, 每次迭代的运算量为 6n; • 计算 T 的特征值时, 假定每个平均迭代 2 步, 则计算所有特征值的运算量为 12n 2 ; • 若要计算 T 的所有特征值和特征向量, 则运算量为 6n 3 + O(n 2 ); • 若只要计算 A 的所有特征值, 运算量为 4n 3/3 + O(n 2 ); • 若需要计算 A 的所有的特征值和特征向量, 则运算量为 26n 3/3 + O(n 2 ); 位移的选取 — Wilkinson 位移 位移的好坏直接影响到算法的收敛速度. 我们可以通过下面的方式来选取位移. 设 A (k) = a (k) 1 b (k) 1 b (k) 1 . . . . . . . . . . . . b (k) n−1 b (k) n−1 a (k) n , 一种简单的位移选取策略就是令 σk = a (k) n . 事实上, a (k) n 就是收敛到特征向量的迭代向量的 Rayleigh 商. 这种位移选取方法几乎对所有的矩阵都有三次渐进收敛速度, 但也存在不收敛的例子, 故我们需要 对其做一些改进. † 如果在 QR 迭代算法中使用位移 σk = a (k) n , 而在 Rayleigh 商迭代算法 5.4 中取 x0 = en = [0, . . . , 0, 1]⊺ , 则利用 QR 迭代与反迭代之间的关系可以证明: QR 迭代算法中的 σk 与 Rayleigh
5.4分而治之法.143 -商送代算法中的p相等a(k)6(k)1的最接近a）的特征值作为位移一种有效的位移是Wilkinson位移，即取子矩阵通过计算可得Wilkinson位移为0=a()+8-sgn(0)82+(6(2)其中8=(a(21-()出于稳定性方面的考虑，我们通常用下面的计算公式(6(2.)g=a(k)8 + sign(0) /82 + (6(2.)定理5.4采用Wilkinson位移的QR选代是整体收敛的，且至少是线性收敛.事实上，几乎对所有的对称矩阵都是渐进三次收敛的口证明.见参考文献[33,57].例5.1带Wilkinson位移的隐式QR选代算法收敛性演示(见[22,page214),相应的MATLAB代码为TriQR.m5.4分而治之法分而治之(Divide-and-Conquer)算法是由Cuppen[21]于1981年首次提出，但直到1995年才出现稳定的实现方式[39].该算法是目前计算维数大于25的矩阵的所有特征值和特征向量的最快算法.下面我们介绍该算法考虑不可约对称三对角矩阵b1ab1am-1bm-1bmbm-1amT=bmbm+1am+1bm+1bn-1
5.4 分而治之法 · 143 · 商迭代算法中的 ρ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(δ) √ δ 2 + ( b (k) n−1 )2 , 其中 δ = 1 2 (a (k) n−1 − a (k) n ). 出于稳定性方面的考虑, 我们通常用下面的计算公式 σ = a (k) n − ( b (k) n−1 )2 δ + sign(δ) √ δ 2 + ( b (k) n−1 )2 . 定理 5.4 采用 Wilkinson 位移的 QR 迭代是整体收敛的, 且至少是线性收敛. 事实上, 几乎对所有的对 称矩阵都是渐进三次收敛的. 证明. 见参考文献 [33, 57]. □ 例 5.1 带 Wilkinson 位移的隐式 QR 迭代算法收敛性演示 (见 [22, page 214]), 相应的 MATLAB 代码为 TriQR.m 5.4 分而治之法 分而治之 (Divide-and-Conquer) 算法是由 Cuppen [21] 于 1981 年首次提出, 但直到 1995 年才出现稳 定的实现方式 [39]. 该算法是目前计算维数大于 25 的矩阵的所有特征值和特征向量的最快算法. 下面 我们介绍该算法. 考虑不可约对称三对角矩阵 T = a1 b1 b1 . . . . . . . . . am−1 bm−1 bm−1 am bm bm am+1 bm+1 bm+1 . . . . . . . . . . . . bn−1 bn−1 an
.144.第五讲对称特征值问题bta1b1bm-am-1bmbm-1-bmbmam-bmbm-bmbm+1am+1bm+1bnh.bmuwt,12其中=[0..,0,1,1,0,..,0Jf假定T和T2的特征值分解已经计算出来了，即T=Q1A1QI,T2=Q2A2Q2，下面考虑T的特征值分解首先介绍一个引理引理5.3设,yERn,则det(I+ry)=1+yTr口证明.见习题5.2.我们首先考虑T的特征值与T和T2的特征值之间的关系T0+bmUuTT20QiAiQI0+bmuT0Q2A2Q1bmuu其中[Q 的最后一列Q11Q2Q的第一列令α=bm,D=diag(A1,A2)=diag(di,d2,,dn),并假设di≥d2≥.≥dn.则T的特征值与D+QuuT的特征值相同下面计算D十guuT的特征值，设入是D十ouuT的一个特征值，若D一>I非奇异，则det(D + αuuT - ^) = det(D - ^) - det(I + α(D - ^I)-1uuT)故det(I+α(D-入I)-1uu)=0.又由引理5.3可知de(I+a(D-AI)-lu)=1+au(D-)-lu=1+aZ (A),d--
· 144 · 第五讲 对称特征值问题 = a1 b1 b1 . . . . . . . . . am−1 bm−1 bm−1 am − bm am+1 − bm bm+1 bm+1 . . . . . . . . . . . . bn−1 bn−1 an + bm bm bm bm = [ T1 0 0 T2 ] + bmvv⊺ , 其中 v = [0, . . . , 0, 1, 1, 0, . . . , 0]⊺ . 假定 T1 和 T2 的特征值分解已经计算出来了, 即 T1 = Q1Λ1Q ⊺ 1 , T2 = Q2Λ2Q ⊺ 2 , 下面考虑 T 的特征值分解. 首先介绍一个引理. 引理 5.3 设 x, y ∈ R n, 则 det(I + xy⊺ ) = 1 + y ⊺x . 证明. 见习题 5.2. □ 我们首先考虑 T 的特征值与 T1 和 T2 的特征值之间的关系. T = [ T1 0 0 T2 ] + bmvv⊺ = [ Q1Λ1Q ⊺ 1 0 0 Q2Λ2Q ⊺ 2 ] + bmvv⊺ = [ Q1 0 0 Q2 ] ([Λ1 0 0 Λ2 ] + bmuu⊺ ) [Q1 0 0 Q2 ]⊺ , 其中 u = [ Q1 0 0 Q2 ]⊺ v = [ Q ⊺ 1 的最后一列 Q ⊺ 2 的第一列 ] . 令 α = bm, D = diag(Λ1,Λ2) = diag(d1, d2, . . . , dn), 并假设 d1 ≥ d2 ≥ · · · ≥ dn. 则 T 的特征值与 D + αuu⊺ 的特征值相同. 下面计算 D + αuu⊺ 的特征值. 设 λ 是 D + αuu⊺ 的一个特征值, 若 D − λI 非奇异, 则 det(D + αuu⊺ − λI) = det(D − λI) · det(I + α(D − λI) −1uu⊺ ). 故 det(I + α(D − λI) −1uu⊺ ) = 0. 又由引理 5.3 可知 det(I + α(D − λI) −1uu⊺ ) = 1 + αu⊺ (D − λI) −1u = 1 + α ∑n i=1 u 2 i di − λ ≜ f(λ)