-134.第四讲非对称特征值问题4.3正交送代法幂选代和反迭代都只能同时计算一个特征对.如果想同时计算多个特征对，我们可以采用多个初始向量进行选代.而正交选代算法就是基于这种思想，它能够计算A的一个不变子空间，从而可以同时计算出多个特征值算法4.4.正交选代法(Orthogonal Iteration)1: Choose an n × p column orthogonal matrix Zo2: set k = 03: while not convergence do4:computeYk+1=AZk5:Yk+1=Zk+1Rk+1%QR分解6:k=k+17: end while正交选代方法有时也称为子空间选代方法（subspaceiterationmethod）和同步选代方法（si-multaneous iteration method).在算法中使用QR分解是为了保持Zk的列正交性，使得其列向量组构成子空间spanAZo的一组正交基。一方面提高算法的数值稳定性，另一方面避免所有列都收敛到最大特征值所对应的特征向量.下面我们分析该算法的收敛性质.假设A是可对角化的，即A=VAV-1,其中A=diag(入1,入2,..,入n),且[^1|≥≥[入p|>[p+1/≥·.≥入nl.则可得span[Zk] = span[Yi] = span[AZk-1],k = 1,2, ...由此可知span[Zk} = span[A* Zo] = span[VA*-1 Zo].我们注意到[(A1 / ^,)k[w()]kV-1Z0=V-1 Zo>[w(b,](An/入p)k由于当i>p时有α/>pl<1,所以当k趋于无穷大时,w(-,趋向于0.令V=[Vp,Vn-p],则W(k)= X (vpw(t) + Vn-pW())VA*V-1 Zo = >[Vp, Vn-p]W(k)
· 134 · 第四讲 非对称特征值问题 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 分解 6: k = k + 1 7: end while b 正交迭代方法有时也称为子空间迭代方法 (subspace iteration method) 和同步迭代方法 (simultaneous iteration method). 在算法中使用 QR 分解是为了保持 Zk 的列正交性, 使得其列向量组构成子空间 span{AkZ0} 的一组正交基. 一方面提高算法的数值稳定性, 另一方面避免所有列都收敛到最大特征值所对应 的特征向量. 下面我们分析该算法的收敛性质. 假设 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
4.3正交送代法?135.所以当k→8时,有span[Zk) = span(VA*V-1 Z0) = span[Vw() + Vn-pw(-)]→ span[V,W(*)] = span[Vp],即span[Zk}趋向于A的一个p维不变子空间span(Vp].定理41给定正整数p(1≤p≤n),考虑算法4.4.假设A是可对角化的,且[^i/ ≥.…≥/^pl>IΛp+1l≥··≥[^nl.则span[Zk}收敛到A的一个p维不变子空间当A不可对角化时，利用Jordan标准型，我们可以得到同样的结论，见[132，133].在正交选代中，如果我们取Zo=I，则可得到一类特殊的正交选代法.此时，在一定条件下正交迭代会收敛到A的Schur分解
4.3 正交迭代法 · 135 · 所以当 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 标准型, 我们可以得到同样的结论, 见 [132, 133]. 在正交迭代中, 如果我们取 Z0 = I, 则可得到一类特殊的正交迭代法. 此时, 在一定条件下, 正交迭代会收敛到 A 的 Schur 分解
-136第四讲非对称特征值问题4.4QR选代法OR选代法的基本思想是通过不断的正交相似变换，使得A趋向于一个上三角形式（或拟上三角形式）.算法形式非常简洁，描述如下：算法4.5.QR选代法(QRIteration)l:Set Ai=Aand k-12: while not convergence do%QR分解3:[Qk, Rk] = QR(Ak)4:computeAk+1=R:Qk5:k=k+16: end while在该算法中，我们有Ak+1=RQk=(QTQ)RQk=QT(QR)Qk=QTAQk由这个递推关系可得Ak+1=QTAkQk= QTQT-1Ak-1Qk-1Qk=-.-- =QTQT-1-.QIAQ1-·Qk-1Qk记Qk=Q1.Q-1Qk=[a..,(]则Ak+1=QTAQk,(4.1)即Ak+1与A正交相似.4.4.1OR选代与幂选代的关系记R=RRk-1R1则有QRk=Qk-1(QR)Rk-1=Qk-1(A)Rk-1=Qk-1(QT-,AQk-1)Rk-1=AQk-1Rk-1由此递推下去，即可得QRK= Ak-1QiR1= Ak-1QiR1= Ak(4.2)故QkRke1 = Ahe1,这说明QR选代与幕选代有关假设il>[2|≥≥[nl,则当k充分大时，Aei收敛到A的模最大特征值>所对应的特征向量，故Q的第一列(）也收敛到入1所对应的特征向量。因此，当k充分大时,Ag(k)→^1g(),此时由(4.1)可知,Ak+1的第一列Ak+1(,1) =QTAg()→A1QTq()=A1e1,即Ak+1的第一列的第一个元素收敛到入1，而其它元素都趋向于0,收敛速度取决于[入2/入1|的大小
· 136 · 第四讲 非对称特征值问题 4.4 QR 迭代法 QR 迭代法的基本思想是通过不断的正交相似变换, 使得 A 趋向于一个上三角形式 (或拟上 三角形式). 算法形式非常简洁, 描述如下: 算法 4.5. QR 迭代法 (QR Iteration) 1: Set A1 = A and k = 1 2: while not convergence do 3: [Qk, Rk] = QR(Ak) % QR 分解 4: compute Ak+1 = RkQk 5: k = k + 1 6: end while 在该算法中, 我们有 Ak+1 = RkQk = (Q T k Qk)RkQk = Q T k (QkRk)Qk = Q T k AkQk. 由这个递推关系可得 Ak+1 = Q T k AkQk = Q T k Q T k−1Ak−1Qk−1Qk = · · · = Q T k Q T k−1 · · · Q T 1 AQ1 · · · Qk−1Qk. 记 Q˜ k = Q1 · · · Qk−1Qk = h q˜ (k) 1 , q˜ (k) 2 , . . . , q˜ (k) n i , 则 Ak+1 = Q˜T k AQ˜ k, (4.1) 即 Ak+1 与 A 正交相似. 4.4.1 QR 迭代与幂迭代的关系 记 R˜ k = RkRk−1 · · · R1, 则有 Q˜ kR˜ k = Q˜ k−1(QkRk)R˜ k−1 = Q˜ k−1(Ak)R˜ k−1 = Q˜ k−1(Q˜T k−1AQ˜ k−1)R˜ k−1 = AQ˜ k−1R˜ k−1, 由此递推下去, 即可得 Q˜ kR˜ k = A k−1Q˜ 1R˜ 1 = A k−1Q1R1 = A k . (4.2) 故 Q˜ kR˜ ke1 = A k e1, 这说明 QR 迭代与幂迭代有关. 假设 |λ1| > |λ2| ≥ · · · ≥ |λn|, 则当 k 充分大时, Ak e1 收敛到 A 的模最大特征值 λ1 所 对应的特征向量, 故 Q˜ k 的第一列 q˜ (k) 1 也收敛到 λ1 所对应的特征向量. 因此, 当 k 充分大时, Aq˜ (k) 1 → λ1q˜ (k) 1 , 此时由 (4.1) 可知, Ak+1 的第一列 Ak+1(:, 1) = Q˜T k Aq˜ (k) 1 → λ1Q˜T k q˜ (k) 1 = λ1e1, 即 Ak+1 的第一列的第一个元素收敛到 λ1, 而其它元素都趋向于 0, 收敛速度取决于 |λ2/λ1| 的大 小
4.4QR送代法·137.4.4.2QR选代与反选代的关系下面观察Q的最后一列.由(4.1)可知AQk=QkAk+1=QQk+1Rk+1=Q+1Rk+1,所以有Qk+1=AQR+1由于Qk+1和Q都是正交矩阵，上式两边转置后求逆，可得Qk+1=(QT+1)=((R+)QTAT)=(AT)QRI+1观察等式两边矩阵的最后一列，可得(+1) =c (AT)- (),其中c1为某个常数.依此类推,可知d(+1) =c(AT)-"dl),其中c为某个常数.假设A的特征值满足/>i|≥[>2|≥：≥/入n-1]>[n>0,则入-1是（AT)-的模最大的特征值。由幂迭代的收敛性可知，q(+1）收敛到（AT)-1的模最大特征值入-1所对应的特征向量，即当k充分大时，有(AT)(k+1) →入"g(+1),所以ATa(k+1) → ng(k+1),由（4.1)可知，AT+,的最后一列为AT+1(c,n) = QTATg(k)→ nQTg(k) = Anen即Ak+1的最后一行的最后一个元素收敛到入n，而其它元素都趋向于0,收敛速度取决于/n/入n-1l的大小4.4.3QR选代与正交选代的关系下面的定理给出了QR送代法与正交迭代法（取Zo=I)之间的关系定理4.2设正交选代法4.4和OR算法4.5中所涉及的OR分解都是唯一的.设Ak是由OR迭代法4.5生成的矩阵，Zk是由正交选代法4.4（取Zo=I）生成的矩阵，则有Ak+1=ZTAZk(板书)证明.我们用归纳法证明该结论当 k=0时,A1=A,Zo=I.结论显然成立.设A=ZT-,AZk-1.由于Zk-1是正交矩阵，我们有Z,Rk= Y= AZk-1 = (Zk-1ZT-1)AZk-1=Zk-1Ak =(Zk-1Qk)Rk
4.4 QR 迭代法 · 137 · 4.4.2 QR 迭代与反迭代的关系 下面观察 Q˜ k 的最后一列. 由 (4.1) 可知 AQ˜ k = Q˜ kAk+1 = Q˜ kQk+1Rk+1 = Q˜ k+1Rk+1, 所以有 Q˜ k+1 = AQ˜ kR −1 k+1. 由于 Q˜ k+1 和 Q˜ k 都是正交矩阵, 上式两边转置后求逆, 可得 Q˜ k+1 = Q˜T k+1−1 = R −1 k+1T Q˜T k A T −1 = A T −1 Q˜ kR T k+1. 观察等式两边矩阵的最后一列, 可得 q˜ (k+1) n = c1 A T −1 q˜ (k) n , 其中 c1 为某个常数. 依此类推, 可知 q˜ (k+1) n = c A T −k q˜ (1) n , 其中 c 为某个常数. 假设 A 的特征值满足 |λ1| ≥ |λ2| ≥ · · · ≥ |λn−1| > |λn| > 0, 则 λ −1 n 是 AT −1 的模最大的特征值. 由幂迭代的收敛性可知, q˜ (k+1) n 收敛到 AT −1 的模最大特征值 λ −1 n 所对应 的特征向量, 即当 k 充分大时, 有 A T −1 q˜ (k+1) n → λ −1 n q˜ (k+1) n . 所以 A T q˜ (k+1) n → λnq˜ (k+1) n . 由 (4.1) 可知, AT k+1 的最后一列为 A T k+1(:, n) = Q˜T k A T q˜ (k) n → λnQ˜T k q˜ (k) n = λnen, 即Ak+1 的最后一行的最后一个元素收敛到λn, 而其它元素都趋向于0, 收敛速度取决于|λn/λn−1| 的大小. 4.4.3 QR 迭代与正交迭代的关系 下面的定理给出了 QR 迭代法与正交迭代法 (取 Z0 = I) 之间的关系. 定理 4.2 设正交迭代法 4.4 和 QR 算法 4.5 中所涉及的 QR 分解都是唯一的. 设 Ak 是由 QR 迭 代法 4.5 生成的矩阵, Zk 是由正交迭代法 4.4 (取 Z0 = I) 生成的矩阵, 则有 Ak+1 = Z T k AZk. (板书) 证明. 我们用归纳法证明该结论. 当 k = 0 时, A1 = A, Z0 = I. 结论显然成立. 设 Ak = Z T k−1AZk−1. 由于 Zk−1 是正交矩阵, 我们有 ZkRˆ k = Yk = AZk−1 = (Zk−1Z T k−1 )AZk−1 = Zk−1Ak = (Zk−1Qk)Rk
-138第四讲非对称特征值问题即ZR和（Zk-1Q)R都是Y的QR分解.由QR分解的唯一性可知，Z=Zk-1QkR=Rk所以ZTAZK=(Zk-1Qh)TA(Zk-1Qk)=QTAkQk=QT(QkR)Qk=RKQk= Ak+1)口即Ak+1=ZTAZk.由归纳法可知，定理结论成立.4.4.4QR选代的收敛性定理43 设A=VAV-1 e Rnxn,其中 A=diag(>1,>2,...,入n),且[>i|>[2|>..>[n|>0若V-1的所有顺序主子矩阵都非奇异（即V-1存在LU分解），则A的对角线以下的元素均收敛到0.(板书)证明.设V=Q,R,是V的QR分解，V-1=LU，是V-1的LU分解，其中L是单位下三角矩阵.则Ak=VA*V-1=Q,RA*L,U=QR(A*LA-)A*U注意到矩阵ALA-k是一个下三角矩阵，且其（i）位置上的元素为[0,i<j,(A*LA-) (i, ) =3 1,i=j,L(i,)/,i>j.由于当>时有/入<1,故当充分大时，/趋向于0.所以我们可以把^L-写成A*LuA-k = I + Ek,其中E满足limE=0.于是A=QuR(I +E)A*U=Q(I+R,ER-1)RA*U(4.3)对矩阵I+R,ER1做QR分解：I+R,ER=1=QERE·由于Ek→0,所以QE→I,REx→I.将其代人(22)可得A=QQERERA*Uu=QQEDk(D-'RERuA*U.),(4.4)其中Dk是对角矩阵，其对角线元素的模均为1，它使得上三角矩阵D-1RE,R,A*U，的对角线元素均为正.这样，(22)就构成A的QR分解.又由（4.2)可知A=QRk，根据QR分解的唯一性，我们可得Qu=QQEDk，RK=D-'RERAUu所以由(4.1)可知Ak+1=QTAQk=(QQED)TVAV-1QQE,Dk=DTQEQTQ,RAR-'Q-QQEDk=DIQERuAR-'QEDk
· 138 · 第四讲 非对称特征值问题 即 ZkRˆ k 和 (Zk−1Qk)Rk 都是 Yk 的 QR 分解. 由 QR 分解的唯一性可知, Zk = Zk−1Qk, Rˆ k = Rk. 所以 Z T k AZk = (Zk−1Qk) TA(Zk−1Qk) = Q T k AkQk = Q T k (QkRk)Qk = RkQk = Ak+1, 即 Ak+1 = Z T k AZk. 由归纳法可知, 定理结论成立. □ 4.4.4 QR 迭代的收敛性 定理 4.3 设 A = V ΛV −1 ∈ R n×n , 其中 Λ = diag(λ1, λ2, . . . , λn), 且 |λ1| > |λ2| > · · · > |λn| > 0. 若 V −1 的所有顺序主子矩阵都非奇异（即 V −1 存在 LU 分解）, 则 Ak 的对角线以下的元素均 收敛到 0. (板书) 证明. 设 V = QvRv 是 V 的 QR 分解, V −1 = LvUv 是 V −1 的 LU 分解, 其中 Lv 是单位下三角矩 阵. 则 A k = V Λ kV −1 = QvRvΛ kLvUv = QvRv(ΛkLvΛ −k )ΛkUv. 注意到矩阵 Λ kLvΛ −k 是一个下三角矩阵, 且其 (i, j) 位置上的元素为 Λ kLvΛ −k (i, j) = 0, i < j, 1, i = j, Lv(i, j)λ k i /λ k j , i > j. 由于当 i > j 时有 |λi/λj | < 1, 故当 k 充分大时, λ k i /λ k j 趋向于 0. 所以我们可以把 Λ kLvΛ −k 写成 Λ kLvΛ −k = I + Ek, 其中 Ek 满足 lim k→∞ Ek = 0. 于是 A k = QvRv(I + Ek)ΛkUv = Qv(I + RvEkR −1 v )RvΛ kUv. (4.3) 对矩阵 I + RvEkR−1 v 做 QR 分解: I + RvEkR−1 v = QEkREk . 由于 Ek → 0, 所以 QEk → I, REk → I. 将其代入 (??) 可得 A k = QvQEkREkRvΛ kUv = QvQEkDk(D −1 k REkRvΛ kUv), (4.4) 其中 Dk 是对角矩阵, 其对角线元素的模均为 1, 它使得上三角矩阵 D −1 k REkRvΛ kUv 的对角线元 素均为正. 这样, (??) 就构成 Ak 的 QR 分解. 又由 (4.2) 可知 Ak = Q˜ kR˜ k, 根据 QR 分解的唯一性, 我们可得 Q˜ k = QvQEkDk, R˜ k = D −1 k REkRvΛ kUv. 所以由 (4.1) 可知 Ak+1 = Q˜T k AQ˜ k = (QvQEkDk) TV ΛV −1QvQEkDk = DT k Q T EkQ T v QvRvΛR −1 v Q −1 v QvQEkDk = DT k Q T EkRvΛR −1 v QEkDk