2/78 非对称矩阵特征值 / 特征向量的计算基本约定 1: A R n n 非对称稠密基本约定 2: λ 1 λ 2 λ n 0 本讲主要讨论如何计算 A 的全部特征值和 / 或特征向量. 主要介绍以下方法 : 幂迭代方法反迭代方法 ( 位移策略,Rayleigh 商迭代 ) 正交迭代方法

Size: px

Start display at page:

Download "2/78 非对称矩阵特征值 / 特征向量的计算基本约定 1: A R n n 非对称稠密基本约定 2: λ 1 λ 2 λ n 0 本讲主要讨论如何计算 A 的全部特征值和 / 或特征向量. 主要介绍以下方法 : 幂迭代方法反迭代方法 ( 位移策略,Rayleigh 商迭代 ) 正交迭代方法"

亡刘
5 years ago
Views:

1 第四讲非对称特征值问题 1 幂迭代 2 反迭代 3 正交迭代 4 QR 迭代 5 带位移的隐式 QR 迭代 6 特征向量的计算 7 广义特征值问题 8 应用 : 多项式求根

2 2/78 非对称矩阵特征值 / 特征向量的计算基本约定 1: A R n n 非对称稠密基本约定 2: λ 1 λ 2 λ n 0 本讲主要讨论如何计算 A 的全部特征值和 / 或特征向量. 主要介绍以下方法 : 幂迭代方法反迭代方法 ( 位移策略,Rayleigh 商迭代 ) 正交迭代方法 QR 方法

3 关于稠密矩阵特征值计算的参考资料有 : J. H. Wilkinson, The Algebraic Eigenvalue Problem, 1965 B. N. Parlett, The Symmetric Eigenvalue Problem, 2nd Eds., 1998 G. W. Stewart, Matrix Algorithms, Vol II: Eigensystems, 2001 G. H. Golub and C. F. Van Loan, Matrix Computations, 2013 P. Arbenz, The course G, Numerical Methods for Solving Large Scale Eigenvalue Problems, ( 该课程的主页 )

4 4/78 1 幂迭代幂迭代是计算特征值和特征向量的一种简单易用的算法. 虽然简单, 但它却建立了计算特征值和特征向量的算法的一个基本框架. 算法 1.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) 5: x (k+1) = y (k+1) / y (k+1) 2 6: µ k+1 = (x (k+1), Ax (k+1) ) % 内积 7: k = k + 1 8: end while

5 5/78 幂迭代的收敛性假设 1: A R n n 可对角化, 即 A = V ΛV 1, 其中 Λ = diag(λ 1,..., λ n ), V = [v 1,..., v n ] C n n, v i 2 = 1 假设 2: λ 1 > λ 2 λ 3 λ n. 由于 V 的列向量组构成 C n 的一组基, 因此 x (0) 可表示为 x (0) = α 1 v 1 + α 2 v α n v n = V [α 1, α 2,..., α n ]. 我们假定 α 1 0, 即 x (0) 不属于 span{v 2, v 3,..., v n } ( 由于 x (0) 是随机选取的, 从概率意义上讲, 这个假设通常是成立的 ).

6 于是我们可得 α 1 α 1 α 1 λ k 1 α 2 A k x (0) = (V ΛV 1 ) k V. = V α 2 Λk. = V α 2 λ k 2. α 2 = α 1 λ k α 1V 1 α n α 1 1 α n ( ) k λ2 λ 1. ( ) k λn λ 1 又 λ i /λ 1 < 1, i = 2, 3,..., n, 所以. α n α n λ k n ( λi lim k λ 1 ) k = 0, i = 2, 3,..., n. 6/78

7 7/78 故当 k 趋向于无穷大时, 向量 [ 1, α ( ) k 2 λ2,..., α n α 1 λ 1 α 1 收敛到 e 1 = [1, 0,..., 0]. ( λn λ 1 ) k ], k = 0, 1, 2,... 所以向量 x (k) = A k x (0) / A k x (0) 2 收敛到 ±v 1, 即 λ 1 的特征向量. 而 µ k = (x (k) ) Ax (k) 则收敛到 v 1Av 1 = λ 1. 幂迭代的收敛快慢取决于 λ 2 /λ 1 的大小, λ 2 /λ 1 越小, 收敛越快. 幂迭代只能用于计算 ( 模 ) 最大的特征值和其相应的特征向量. 当 λ 2 /λ 1 接近于 1 时, 收敛速度会非常慢. 如果模最大的特征值是一对共轭复数, 则幂迭代可能会失效.

8 1.1 位移策略出发点 : 加快幂迭代算法的收敛速度尽可能地减小 λ 2 /λ 1 策略 : 位移策略 (shift), 即计算 A σi 的特征值我们称 σ 为位移 (shift), 满足 (1) λ 1 σ 是 A σi 的模最大的特征值 ; (2) max λ i σ 2 i n λ 1 σ 尽可能地小. 其中第一个条件保证最后所求得的特征值是我们所要的, 第二个条件用于加快幂迭代的收敛速度. 缺点 : (1) σ 很难选取 ; (2) 加速效果有限改进方法 : 与反迭代相结合, 能起到很好的加速效果 8/78

9 9/78 2 反迭代用幂迭代求 A 1 的模最小特征值, 这就是反迭代算法 2.1 反迭代算法 (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) 1 x (k) 5: x (k+1) = y (k+1) / y (k+1) 2 6: µ k+1 = (x (k+1), Ax (k+1) ) 7: σ = µ k+1, k = k + 1 8: end while 显然 : µ k 收敛到距离 σ 最近的特征值, x (k) 收敛到对应的特征向量

10 10/78 事实上, 反迭代 + 位移策略, 可以计算矩阵的任意一个特征值优点 : 若 σ 与某个特征值 λ k 非常接近, 则反迭代算法的收敛速度非常快. 只要选取合适的位移 σ, 就可以计算 A 的任意一个特征值. 缺点 : 每步迭代需要解一个线性方程组 (A σi)y (k+1) = x (k), 这需要对 A σi 做 LU 或 PLU 分解. 与幂迭代一样, 反迭代算法一次只能求一个特征值. 怎样选取位移 σ? Rayleigh 商 : 动态选取, 自动调整

11 11/ Rayleigh 商迭代出发点 : 使得 σ 与所求的特征值越靠近越好. 期望能直接给出一个理想位移是不太现实的. 比较现实的方法就是动态调整, 使得位移逐渐靠近某个特征值. Rayleigh 商迭代 : 以 Rayleigh 商 µ k 为第 k 步的位移理由 : µ k 会逐渐收敛到某个特征值.

12 12/78 算法 2.2 Rayleigh 商迭代 (Rayleigh Quotient Iteration, RQI) 1: Choose an initial vector x (0) with x (0) 2 = 1 2: set k = 0 3: compute σ = (x (0) ) Ax (0) 4: while not converge do 5: y (k+1) = (A σi) 1 x (k) 6: x (k+1) = y (k+1) / y (k+1) 2 7: µ k+1 = (x (k+1), Ax (k+1) ) 8: σ = µ k+1 9: k = k : end while

13 13/78 RQI 算法的收敛性一般来说, 如果 Rayleigh 商迭代收敛到 A 的一个单特征值, 则至少是二次收敛的, 即具有局部二次收敛性. 如果 A 是对称的, 则能达到局部三次收敛, 详情见后面的对称特征值问题. 缺点 : 由于每次迭代的位移是不同的, 因此每次迭代需要求解一个不同的线性方程组, 这使得运算量大大增加. 因此通常应用于三对角矩阵的特征值计算.

14 14/78 3 正交迭代出发点 : 同时计算多个特征值 / 特征向量策略 : 同时采用多个初始向量, 希望收敛到 A 的一个不变子空间算法 3.1 正交迭代算法 (Orthogonal Iteration) 1: Choose an n p column orthogonal matrix Z 0 2: set k = 0 3: while not convergence do 4: compute Y k+1 = AZ k 5: Y k+1 = Z k+1 ˆRk+1 % QR 分解 6: k = k + 1 7: end while

15 说明 : 在算法中使用 QR 分解是为了保持 Z k 的列正交性, 使得其列向量组构成子空间 span{a k Z 0 } 的一组正交基. 一方面提高算法的数值稳定性, 另一方面避免所有列都收敛到最大特征值所对应的特征向量. 15/78

16 收敛性分析假设 A 是可对角化的, 即 A = V ΛV 1, 其中 Λ = diag(λ 1, λ 2,..., λ n ), 且 λ 1 λ p > λ p+1 λ n. 则可得由此可知 span{z k } = span{y k } = span{az k 1 }, k = 1, 2,..., span{z k } = span{a k Z 0 } = span{v Λ k V 1 Z 0 }. 我们注意到 (λ 1 /λ p ) k... Λ k V 1 Z 0 = λ k p 1... (λ n /λ p ) k V 1 Z 0 = λ k p [ W (k) p W (k) n p ]. 16/78

17 17/78 由于当 i > p 时有 λ i /λ p < 1, 所以当 k 趋于无穷大时, W (k) n p 趋向于 0. 令 V = [V p, V n p ], 则 [ ] V Λ k V 1 Z 0 = λ k W p (k) p[v p, V n p ] W (k) = λ k p n p 所以当 k 时, 有 span{z k } = span{v Λ k V 1 Z 0 } = span{v p W (k) p ( ) V p W p (k) + V n p W (k) n p. + V n p W (k) n p} span{v p W (k) p } = span{v p }, 即 span{z k } 趋向于 A 的一个 p 维不变子空间 span{v p }. 定理给定正整数 p (1 p n), 考虑算法 3.1. 假设 A 是可对角化的, 且 λ 1 λ p > λ p+1 λ n. 则 span{z k } 收敛到 A 的一个 p 维不变子空间.

18 18/78 说明 : 如果 A 不可对角化, 利用 Jordan 标准型, 可以到同样的结论, 见 [Watkins 2007, Watkins-Elsner 1991]. 在正交迭代中, 如果我们取 Z 0 = I, 则可得到一类特殊的正交迭代算法. 此时, 在一定条件下, 正交迭代会收敛到 A 的 Schur 标准型.

19 19/78 4 QR 迭代 4.1 算法介绍 4.2 QR 迭代与幂迭代的关系 4.3 QR 迭代与反迭代的关系 4.4 QR 迭代与正交迭代的关系 4.5 QR 迭代的收敛性 4.6 带位移的 QR 迭代

20 20/ 算法介绍基本思想 : 通过不断的正交相似变换, 将 A 转化为 ( 拟 ) 上三角形式算法 4.1 QR 迭代算法 (QR Iteration) 1: Set A 1 = A and k = 1 2: while not convergence do 3: A k = Q k R k % QR 分解 4: compute A k+1 = R k Q k 5: k = k + 1 6: end while

21 21/78 正交相似性在 QR 迭代算法中, 我们有 A k+1 = R k Q k = (Q k Q k)r k Q k = Q k (Q kr k )Q k = Q k A kq k. 由这个递推关系可得 A k+1 = Q k A kq k = = Q k Q k 1 Q 1 AQ 1 Q k 1 Q k. 记 Qk = Q 1 Q k 1 Q k = [ ] q (k) 1, q(k) 2,..., q(k) n, 则 A k+1 = Q k A Q k (4.1) 即 A k+1 与 A 正交相似.

22 22/ QR 迭代与幂迭代的关系记 Rk = R k R k 1 R 1, 则有 Q k Rk = Q k 1 (Q k R k ) R k 1 = Q k 1 (A k ) R k 1 = Q k 1 ( Q k 1 A Q k 1 ) R k 1 = A Q k 1 Rk 1, 由此递推下去, 即可得 Q k Rk = A k 1 Q1 R1 = A k 1 Q 1 R 1 = A k 故 Q k Rk e 1 = A k e 1

23 23/78 假设 λ 1 > λ 2 λ n, 则当 k 充分大时, A k e 1 收敛到 A 的模最大特征值 λ 1 所对应的特征向量. 故 Qk 的第一列 q (k) 1 也收敛到 λ 1 所对应的特征向量因此, 当 k 充分大时, A q (k) 1 λ 1 q (k) 1 由 A k+1 = Q k A Q k 可知, A k+1 的第一列 A k+1 (:, 1) = Q k A q(k) 1 λ 1 Q k q(k) 1 = λ 1 e 1 结论 A k+1 的第一列的第一个元素收敛到 λ 1, 而其它元素都趋向于 0. 收敛速度取决于 λ 2 /λ 1 的大小.

24 24/ QR 迭代与反迭代的关系观察 Qk 的最后一列. 由 A k+1 = Q k A Q k 可知所以有 A Q k = Q k A k+1 = Q k Q k+1 R k+1 = Q k+1 R k+1, Q k+1 = A Q k R 1 k+1. 由于 Qk+1 和 Qk 都是正交矩阵, 上式两边转置后求逆, 可得 ( ) 1 ( (R ) Q k+1 = Q k+1 = 1 k+1 Q k A ) 1 = (A ) 1 Qk R k+1. 观察等式两边矩阵的最后一列, 可得依此类推, 可知 q (k+1) n = c 1 (A ) 1 q (k) n, (c 1 为某个常数 ) q (k+1) n = c (A ) k q (1) n (c 为某个常数 )

25 25/78 假定 λ 1 λ n 1 > λ n > 0, 则 λ 1 n 是 (A ) 1 的模最大特征值. 由幂迭代可知, q (k+1) n 所以收敛到 λ 1 n (A ) 1 q (k+1) n 所对应的特征向量, 即 λ 1 n q n (k+1) (k ) (k) A q n λ n q n (k) (k ) 由 A k+1 = Q k A Q k 可知, A k+1 的最后一列 A k+1 (:, n) = Q (k) ka q n λ n Q k q(k) n = λ n e n. 结论 A k+1 的最后一行的最后一个元素收敛到 λ n, 而其它元素都趋向于 0. 收敛速度取决于 λ n /λ n 1 的大小.

26 26/ QR 迭代与正交迭代的关系下面的定理给出了 QR 迭代算法与正交迭代算法 (Z 0 = I) 之间的关系. 定理假定正交迭代算法 3.1 和 QR 算法 4.1 中所涉及的 QR 分解都是唯一的. A k 是由 QR 迭代算法 4.1 生成的矩阵, Z k 是由正交迭代算法 3.1 ( 取 Z 0 = I) 生成的矩阵, 则有 A k+1 = Z k AZ k. ( 板书 )

27 27/ QR 迭代的收敛性定理设 A = V ΛV 1 R n n, 其中 Λ = diag(λ 1, λ 2,..., λ n ), 且 λ 1 > λ 2 > > λ n. 若 V 1 的所有主子矩阵都非奇异 ( 即 V 1 存在 LU 分解 ), 则 A k 的对角线以下的元素均收敛到 0. ( 板书 ) 说明 : 需要指出的是, 由于 D k 的元素不一定收敛, 故 A k+1 对角线以上 ( 不含对角线 ) 的元素不一定收敛, 但这不妨碍 A k+1 的对角线元素收敛到 A 的特征值 ( 即 A k+1 的对角线元素是收敛的 ).

28 28/78 例 QR 迭代算法演示 ( 见 Eig_QR_demo.m). 设 9 5 A = X 3 X 1, 1 其中 X 是由 MATLAB 随机生成的非奇异矩阵. 在迭代过程中, 对于 A k 的下三角部分中元素, 如果其绝对值小于某个阈值 tol, 则直接将其设为 0, 即这里我们取 tol = 10 6 a (k) ij = 0 if i > j and a (k) ij < tol. max 1 i,j n { a(k) ij }, 迭代过程如下 :

29 A = e e e e e e e e e e e e e e e e-01 A_7 = e e e e e e e e e e e e e e e-01 A_8 = e e e e e e e e e e e e e e-01 29/78

30 A_12 = e e e e e e e e e e e e e+00 A_13 = e e e e e e e e e e e e+00 A_22 = e e e e e e e e e e e+00 30/78

31 A_28 = e e e e e e e e e e+00 31/78

32 32/ 带位移的 QR 迭代为了加快 QR 迭代的收敛速度, 可以采用位移策略和反迭代思想. 算法 4.2 带位移的 QR 迭代算法 (QR Iteration with shift) 1: Set A 1 = A and k = 1 2: while not convergence do 3: Choose a shift σ k 4: A k σ k I = Q k R k % QR 分解 5: Compute A k+1 = R k Q k + σ k I 6: k = k + 1 7: end while

33 33/78 正交相似性 : A k+1 = R k Q k + σ k I = (Q k Q k)r k Q k + σ k I = Q k (A k σ k I)Q k + σ k I = Q k A kq k

34 34/78 位移 σ k 的选取在前面的分析可知, A k+1 (n, n) 收敛到 A 的模最小特征值. 若 σ k 就是 A 的一个特征值, 则 A k σ k I 的模最小特征值为 0, 故 QR 算法迭代一步就收敛. 此时 [ ] A (n 1) (n 1) k+1 A k+1 = R k Q k + σ k I =. 0 σ k A 的其它特征值可通过对 A (n 1) (n 1) k+1 使用带位移 QR 迭代算法得到. 通常, 如果 σ k 与 A 的某个特征值非常接近, 则收敛速度通常会很快. 由于 A k (n, n) 收敛到 A 的一个特征值, 所以在实际使用中, 一个比较直观的位移选择策略是 σ k = A k (n, n). 事实上, 这样的位移选取方法通常会使得 QR 迭代算法有二次收敛速度.

35 35/78 例带位移的 QR 迭代算法演示 ( 见 Eig_QR_shift.m). 所有数据和设置与例 4.1 相同, 在迭代过程中, 取 σ k = A k (n, n). 如果 A k (n, n) 已经收敛, 则取 σ k = A k (n 1, n 1). A = e e e e e e e e e e e e e e e e-01 A_5 = e e e e e e e e e e e e e+00

36 A_7 = e e e e e e e e e e e+00 A_10 = e e e e e e e e e e+00 36/78

37 37/78 5 带位移的隐式 QR 迭代直接实施 QR 方法的困难 : 运算量每一步迭代需要做一次 QR 分解和矩阵乘积, 运算量为 O(n 3 ). 即使每计算一个特征值只需迭代一步, 则总运算量为 O(n 4 ). 我们的目标 : 从 O(n 4 ) 减小到 O(n 3 ) 实现方法 : 两个步骤 (1) 首先通过相似变化将 A 转化成一个上 Hessenberg 矩阵 (2) 对这个 Hessenberg 矩阵实施隐式 QR 迭代

38 38/78 隐式 QR 迭代 : 在 QR 迭代算法中, 并不进行显式的 QR 分解和矩阵乘积, 而是通过特殊手段来实现从 A k 到 A k+1 的迭代, 并且将运算量控制在 O(n 2 ) 量级, 从而将总运算量降到 O(n 3 ).

39 39/ 上 Hessenberg 矩阵上 Hessenberg 矩阵 : H = [h ij ] R n n, 当 i > j + 1 时, 有 h ij = 0 定理设 A R n n, 则存在正交矩阵 Q R n n, 使得 QAQ 是上 Hessenberg 矩阵. 下面我们以一个 5 5 的矩阵 A 为例, 给出具体的转化过程, 采用的工具为 Householder 变换.

40 40/78 第一步 : 令 Q 1 = diag(i 1 1, H 1 ), 其中 H 1 是对应于向量 A(2 : 5, 1) 的 Householder 矩阵. 于是可得 Q 1 A = 由于用 Q 1 右乘 Q 1 A 时, 不会改变 Q 1 A 第一列元素的值, 故 A 1 = Q 1 AQ 1 =

41 41/78 第二步 : 令 Q 2 = diag(i 2 2, H 2 ), 其中 H 2 是对应于向量 A 1 (3 : 5, 2) 的 Householder 矩阵, 则用 Q 2 左乘 A 1 时, 不会改变 A 1 的第一列元素的值. 用 Q 2 右乘 Q 2 A 1 时, 不会改变 Q 2 A 1 前两列元素的值. 因此, Q 2 A 1 = 0 和 A 2 = Q 2 A 1 Q 2 =

42 42/78 第三步 : 令 Q 3 = diag(i 3 3, H 3 ), 其中 H 3 是对应于向量 A 2 (4 : 5, 3) 的 Householder 矩阵, 则有 Q 3 A 2 = 0 和 A 3 = Q 3 A 2 Q 3 = 这时, 我们就将 A 转化成一个上 Hessenberg 矩阵, 即 QAQ = A 3 其中 Q = Q 3 Q 2 Q 1 是正交矩阵, A 3 是上 Hessenberg 矩阵.

43 43/78 上 Hessenberg 化算法算法 5.1 上 Hessenberg 化算法 (Upper Hessenberg Reduction) 1: Set Q = I 2: for k = 1 to n 2 do 3: compute Householder matrix H k with respect to A(k + 1 : n, k) 4: A(k + 1 : n, k : n) = H k A(k + 1 : n, k : n) = A(k + 1 : n, k : n) β k v k (v ka(k + 1 : n, k : n)) 5: A(1 : n, k + 1 : n) = A(1 : n, k + 1 : n) H k = A(1 : n, k + 1 : n) β k A(1 : n, k + 1 : n)v k v k 6: Q(k + 1 : n, k : n) = H k Q(k + 1 : n, k : n) = Q(k + 1 : n, k : n) β k v k (v kq(k + 1 : n, k : n)) 7: end for

44 44/78 说明 : 在实际计算时, 我们不需要显式地形成 Householder 矩阵 H k. 上述算法的运算量大约为 14 3 n3 + O(n 2 ). 如果不需要计算特征向量, 则正交矩阵 Q 也不用计算, 此时运算量大约为 10 3 n3 + O(n 2 ). 上 Hessenberg 矩阵的一个很重要的性质就是在 QR 迭代中保持形状不变定理设 A R n n 是非奇异上 Hessenberg 矩阵, 其 QR 分解为 A = QR, 则 Ã RQ 也是上 Hessenberg 矩阵. ( 板书 ) 若 A 是奇异的, 也可以通过选取适当的 Q, 使得上述结论成立.

45 45/78 由这个性质可知, 如果 A R n n 是上 Hessenberg 矩阵, 则 QR 迭代中的每一个 A k 都是上 Hessenberg 矩阵. 这样, 在进行 QR 分解时, 运算量可大大降低. Hessenberg 矩阵另一个重要性质 : 在 QR 迭代过程中保持不可约性. 定理设 A R n n 是不可约上 Hessenberg 矩阵, 其 QR 分解为 A = QR, 则 Ã RQ 也是不可约上 Hessenberg 矩阵. ( 证明留作练习 ) 推论设 A R n n 是不可约上 Hessenberg 矩阵, 则在带位移的 QR 迭代中, 所有的 A k 均是不可约上 Hessenberg 矩阵.

46 46/ 隐式 QR 迭代在 QR 迭代中, 我们要先做 QR 分解 A k = Q k R k, 然后计算 A k+1 = R k Q k. 但事实上, 我们可以将这个过程进行简化, 即在不对 A k 进行 QR 分解的前提下, 直接计算出 A k+1. 这就是隐式 QR 迭代. 不失一般性, 我们假定 A 是不可约的上 Hessenberg 矩阵. 隐式 QR 迭代的理论基础就是下面的隐式 Q 定理. 定理 (Implicit Q Theorem) 设 H = Q AQ R n n 是一个不可约上 Hessenberg 矩阵, 其中 Q R n n 是正交矩阵, 则 Q 的第 2 至第 n 列均由 Q 的第一列所唯一确定 ( 可相差一个符号 ). ( 板书 )

47 47/78 由于 Q k 的其它列都由 Q k 的第一列唯一确定 ( 至多相差一个符号 ), 所以我们只要找到一个正交矩阵 Qk 使得其第一列与 Q k 的第一列相等, 且 Q k A k Q k 为上 Hessenberg 矩阵, 则由隐式 Q 定理可知 Qk = W Q k, 其中 W = diag(1, ±1,..., ±1). 于是 Q k A k Q k = W Q k A kq k W = W A k+1 W. 又 W A k+1 W 与 A k+1 相似, 且对角线元素相等, 而其它元素也至多相差一个符号, 所以不会影响 A k+1 的收敛性, 即下三角元素收敛到 0, 对角线元素收敛到 A 的特征值. 在 QR 迭代算法中, 如果我们直接令 A k+1 = Q k A k Q k, 则其收敛性与原 QR 迭代算法没有任何区别! 这就是隐式 QR 迭代的基本思想. 由于 A 是上 Hessenberg 矩阵, 因此在实际计算中, 我们只需 Givens 变换.

48 48/78 下面我们举一个例子, 具体说明如何利用隐式 Q 定理, 由 A 1 得到 A 2. 设 A R 5 5 是一个不可约上 Hessenberg 矩阵, 即 A 1 = A =

49 49/78 第一步 : 构造一个 Givens 变换 c 1 s 1 G 1 G(1, 2, θ 1) = s 1 c 1 (c 1, s 1 待定 ) I 3 于是有 G 1 A = 和 A (1) G 1 AG 1 = 与 A 1 相比较, A (1) 在 (3, 1) 位置上多出一个非零元, 我们把它记为 +, 并称之为 bulge. 在下面的计算过程中, 我们的目标就是将其赶出矩阵, 从而得到一个新的上 Hessenberg 矩阵, 即 A 2.

50 50/78 第二步 : 为了消去这个 bulge, 我们可以构造 Givens 变换 1 G 2 G(2, 3, θ c 2 s 2 2) = 使得 G 2 s 2 c 2 A(1) = I 为了保持与原矩阵的相似性, 需要再右乘 G 2, 所以 A (2) G 2 A(1) G 2 = 此时, bugle 从 (3, 1) 位置被赶到 (4, 2) 位置.

51 51/78 第三步 : 与第二步类似, 构造 Givens 变换 I 2 G 3 G(3, 4, θ c 3 s 3 3) = 使得 G 3 s 3 c 3 A(2) = 这时 A (3) G 3 A(2) G 3 = 于是, bugle 又从 (4, 2) 位置又被赶到 (5, 3) 位置.

52 52/78 第四步 : 再次构造 Givens 变换 1 1 G 4 G(4, 5, θ 4) = 1 使得 G 4 A(3) = 0 c 4 s s 4 c 于是 A (4) G 4 A(3) G 4 = 现在, bulge 被赶出矩阵, A (4) 就是我们所要的矩阵!

53 53/78 算法分析, 以及 c 1, s 1 的取值常规 QR 迭代 : A 1 = Q 1 R 1, A 2 = R 1 Q 1 = A 2 = Q 1 A 1Q 1 根据前面的计算过程, 有 A (4) = G 4 G 3 G 2 G 1 A 1G 1 G 2 G 3 G 4 = Q 1 A 1 Q 1, 其中 Q1 = G 1 G 2 G 3 G 4 = A (4) = Q 1 A 1 Q 1 通过直接计算可知, Q 1 的第一列为 [c 1, s 1, 0, 0, 0]. 如果将其取为 A 1 的第一列 [a 11, a 21, 0,..., 0] 单位化后的向量, 则 Q1 的第一列与 Q 1 的第一列相同! = A (4) = W A 2 W

54 54/78 针对带位移的 QR 方法, 我们取 A 1 σ 1 I 的第一列 [a 11 σ 1, a 21, 0,..., 0] 单位化后的向量作为 G 1 的第一列即可. 运算量 : 如果 A R n n 是上 Hessenberg 矩阵, 则使用上面的算法, 带位移 QR 迭代中每一步的运算量为 6n 2 + O(n).

55 55/ 位移的选取通常, 位移越离某个特征值越近, 则收敛速度就越快. 由习题 4.10 可知, 如果位移 σ 与某个特征值非常接近, 则 A k (n, n) σ 就非常接近于 0. 这说明 A k (n, n) 通常会首先收敛到 A 的一个特征值. 所以 σ = A k (n, n) 是一个不错的选择. 但是, 如果这个特征值是复数, 这种位移选取方法就可能失效.

56 56/78 双位移策略设 σ C 是 A 的某个复特征值 λ 的一个很好的近似, 则其共轭 σ 也应该是 λ 的一个很好的近似. 因此我们可以考虑双位移策略, 即先以 σ 为位移迭代一次, 然后再以 σ 为位移迭代一次, 如此不断交替进行迭代. 这样就有 A 1 σi =Q 1 R 1, A 2 =R 1 Q 1 + σi, A 2 σi =Q 2 R 2, A 3 =R 2 Q 2 + σi. 容易验证 A 3 = Q 2 A 2Q 2 = Q 2Q 1A 1 Q 1 Q 2 = Q A 1 Q, 其中 Q = Q 1 Q 2.

57 57/78 我们注意到 σ 可能是复的, 所以 Q 1 和 Q 2 都可能是复矩阵. 但我们却可以选取适当的 Q 1 和 Q 2, 使的 Q = Q 1 Q 2 是实正交矩阵. 引理在双位移 QR 迭代中, 我们可以选取酉矩阵 Q 1 和 Q 2 使得 Q = Q 1 Q 2 是实矩阵. ( 板书 )

58 58/78 双位移策略的实现由前面的结论可知, 存在 Q 1 和 Q 2, 使得 Q = Q 1 Q 2 是实矩阵, 从而 A 3 = Q A 1 Q 也是实矩阵. 因此我们希望不计算 A 2, 而是直接从 A 1 得到 A 3. 实现方法 : 根据隐式 Q 定理 : 只要找到一个实正交矩阵 Q, 使得其第一列与 A 2 1 2Re(σ)A 1 + σ 2 I 的第一列平行, 并且 A 3 = Q A 1 Q 是上 Hessenberg 矩阵即可.

59 59/78 易知, A 2 1 2Re(σ)A 1 + σ 2 I 的第一列为 a a 12 a 21 2Re(σ)a 11 + σ 2 a 21 (a 11 + a 22 2Re(σ)) a 21 a 32. (4.5) 0 所以 Q 的第一列是上述向量的单位化.. 其它过程可以通过隐式 QR 迭代来实现. 但此时的 bulge 是一个 2 2 的小矩阵. 因此, 在双位移隐式 QR 迭代过程中, 需要使用 Householder 变换. 需要指出的是, 双位移 QR 迭代算法中的运算都是实数运算.

60 60/78 下面通过一个例子来说明如何在实数运算下实现双位移隐式 QR 迭代. 设 A R 6 6 是一个不可约上 Hessenberg 矩阵, 即 A 1 = A =

61 61/78 [ ] H 1 第一步 : 构造一个正交矩阵 H 1 = 0, 其中 H1 R 3 3, 使得其第 0 I 3 一列与 A 2 1 2Re(σ)A 1 + σ 2 I 的第一列平行. 于是有 H 1 A = 和 A (1) H 1 AH 1 = 与 A 1 相比较, A (1) 在 (3, 1), (4, 1) 和 (4, 2) 位置上出现 bulge. 在下面的计算过程中, 我们的目标就是把它们赶出矩阵, 从而得到一个新的上 Hessenberg 矩阵.

62 62/ 第二步 : 令 H 2 = 0 H 2 0, 其中 H2 R 3 3 是对应于 A(2 : 4, 1) 的 0 0 I 2 Householder 变换, 使得 H 2 A(1) = 和 A (2) H 2 A(1) H 2 = 这时, 我们将 bugle 向右下角方向赶了一个位置.

63 63/78 I 第三步 : 与第二步类似, 令 H 3 = 0 H 3 0, 其中 H3 R 3 3 是对应于 A(3 : 5, 2) 的 Householder 变换, 使得 H 3 A(2) = 和 A (3) H 3 A(2) H 3 = 此时, bugle 又被向右下角方向赶了一个位置.

64 64/78 [ ] I 3 0 第四步 : 令 H 4 = 0 H, 其中 H4 R 3 3 是对应于 A(4 : 6, 3) 的 4 Householder 变换, 使得 H 4 A(3) = 和 A (4) H 4 A(3) H 4 =

65 65/78 第五步 : 只需构造一个 Givens 变换 G 5 = [ ] I G(4, 5, θ), 使得 G 5 A(4) = 和 A (5) G 5 A(4) G 5 = 现在, bulge 已经被全部消除, 且 A (5) = Q AQ, 其中 Q = H 1 H 2 H 3 H 4 G 5. 通过直接计算可知, Q 的第一列即为 H 1 的第一列. 根据隐式 Q 定理, 可以直接令 A 3 A (5) = Q AQ

66 66/78 位移的具体选取在单位移 QR 迭代算法中, 若 A 的特征值都是实的, 则取 σ k = A k (n, n). 推广到复共轭特征值上, 我们可以取 A k 的右下角矩阵 [ ] A k (n 1, n 1) A k (n 1, n) A k (n, n 1) A k (n, n) 的复共轭特征值作为双位移. 这样选取的位移就是 Francis 位移. 如果上述矩阵的两个特征值都是实的, 则选取其中模较小的特征值做单位移.

67 采用 Francis 位移的 QR 迭代会使得 A k 的右下角收敛到一个上三角矩阵 ( 两个实特征值 ) 或一个 2 阶的矩阵 ( 一对复共轭特征值 ), 而且通常会有二次收敛性. 在实际计算中, 一个特征值一般平均只需迭代两步. 收敛性判断 : 判断收敛主要是看 A k (n 1, n 2) ( 或 A k (n, n 1)) 是否趋向于 0. 需要指出的是, QR 迭代并不是对所有矩阵都收敛. 例如 : A = 对于上面的矩阵, 采用 Francis 位移的 QR 迭代算法无效. 另外, 也可以考虑多重位移策略, 参见 [Watkins 2007]. 67/78

68 68/ 收缩 Deflation 收缩 (deflation) 技术是实用 QR 迭代中的一个非常重要概念. 在隐式 QR 迭代过程中, 当矩阵 A k+1 的某个下次对角线元素 a i+1,i 很小时, 我们可以将其设为 0. 由于 A k+1 是上 Hessenberg 矩阵, 这时 A k+1 就可以写成分块上三角形式, 其中两个对角块都是上 Hessenberg 矩阵. 因此我们可以将隐式 QR 迭代作用在这两个规模相对较小的矩阵上, 从而可以大大节约运算量.

69 69/78 6 特征向量的计算设 A 的特征值都是实的, R = Q AQ 是其 Schur 标准型. 若 Ax = λx, 则 Ry = λy, 其中 y = Q x 或 x = Qy. 故只需计算 R 的特征向量 y 即可. 因为 R 的对角线元素即为 A 的特征值, 不妨设 λ = R(i, i). 假定 λ 是单重特征值, 则方程 (R λi)y = 0 即为 R 11 λi R 12 R 13 y R 23 y 2 = 0, 0 0 R 33 λi y 3

70 即 (R 11 λi)y 1 + R 12 y 2 + R 13 y 3 = 0, (4.6) R 23 y 3 = 0, (4.7) (R 33 λi)y 3 = 0, (4.8) 其中 R 11 R (i 1) (i 1), R 33 R (n i) (n i). 由于 λ 是单重特征值, 故 R 33 λi 非奇异, 因此 y 3 = 0. 令 y 2 = 1, 则可得 y 1 = (R 11 λi) 1 R 12. 因此计算特征向量 y 只需求解一个上三角线性方程组. 若 λ 是多重特征值, 则计算方法类似. 但如果 A 有复特征值, 则需要利用实 Schur 标准型, 计算较复杂.

71 71/78 7 广义特征值问题设 A, B R n n, 若存在 λ C 和非零向量 x C n 使得 Ax = λbx, 则称 λ 为矩阵对 (A, B) 的特征值, x 为相应的特征向量. 计算矩阵对 (A, B) 的特征值和特征向量就是广义特征值问题当 B 非奇异时, 广义特征值问题就等价于标准特征值问题其中 y = Bx. B 1 Ax = λx 或 AB 1 y = λy,

72 容易看出, λ 是 (A, B) 的一个特征值当且仅当 det(a λb) = 0. (4.9) 若 (4.9) 对所有 λ C 都成立, 则称矩阵对 (A, B) 是奇异矩阵对, 否则称为正则矩阵对. 当 B 非奇异时, 特征方程 (4.9) 是一个 n 次多项式, 因此恰好有 n 个特征值. 当 B 奇异时, 特征方程 (4.9) 的次数低于 n, 因此方程的解的个数小于 n. 但是, 注意到 λ 0 是 (A, B) 的特征值当且仅当 µ = 1 λ 是 (B, A) 的特征值. 因此, 当 B 奇异时, µ = 0 是 (B, A) 的特征值, 于是我们自然地把 λ = 1 µ = 当作是 (A, B) 的特征值. 所以, 广义特征值不是分布在 C 上, 而是分布在 C { } 上. 容易验证, 若 U, V 非奇异, 则矩阵对 (U AV, U BV ) 的特征值与 (A, B) 是一样的. 因此我们称这种变换为矩阵对的等价变换. 如果 U, V 是酉矩阵, 则称为酉等价变换.

73 73/ 广义 Schur 分解广义 Schur 分解是矩阵对在酉等价变化下的最简形式. 定理 ( 广义 Schur 分解 ) 设 A, B C n n, 则存在酉矩阵 Q, Z C n n, 使得 Q AZ = R A, Q BZ = R B, (4.10) 其中 R A, R B C n n 都是上三角矩阵. 此时矩阵对 (A, B) 的特征值为 R A 和 R B 的对角线元素的比值, 即 λ i = R A(i, i), i = 1, 2,..., n. R B (i, i) 当 R B (i, i) = 0 时, 对应的特征值 λ i =. 证明. 参见 [Xu-Qian 2011].

74 74/78 与实 Schur 分解类似, 当 A, B 都是实矩阵时, 我们有相应的广义实 Schur 分解. 定理 ( 广义实 Schur 分解 ) 设 A, B R n n, 则存在正交矩阵 Q, Z R n n, 使得 Q AZ = T A, Q BZ = T B, (4.11) 其中 T A, T B R n n 都是拟上三角矩阵. 证明. 参见 [Xu-Qian 2011].

75 75/ QZ 迭代 QZ 迭代是用于计算 (A, B) 的广义 Schur 分解的算法, 是 QR 算法的自然推广, 实质上可以看作是将 QR 算法作用到矩阵 AB 1 上. 详细算法可参见 [Kressner 2005, Xu-Qian 2011].

76 76/78 8 应用 : 多项式求根考虑 n 次多项式 q n (x) = x n + c n 1 x n c 1 x + c 0, c i R. 由代数学基本定理可知, p n (x) 在复数域中有且仅有 n 的零点 n 5 时, 不存在求根公式非线性迭代方法求解 MATLAB 中的 roots 命令 : 通过特征值计算方法求出所有零点

77 77/78 友矩阵 0 c c 1 A = c n 1. 多项式 q n (x) 的零点 A 的特征值无需上 Hessenberg 化 A 非常稀疏, 但经过一步 QR 迭代后, 上三角部分的零元素会消失, 总运算量仍是 O(n 3 ) 快速 QR 方法 : 利用 A 的特殊结构, 运算量降为 O(n 2 ) - 将 A 写成一个酉矩阵与秩一矩阵之差, 具体实现参见相关文献

78 课后习题见讲义

슬라이드 1

슬라이드 1 08-09 年度第一学期 0050 00503 计算方法 (B) 童伟华管理科研楼 05 室 E-mail: togwh@ustc.edu.c 中国科学技术大学数学科学学院 http://math.ustc.edu.c/ 第七章计算矩阵的特征值与特征向量特征值与特征向量在实际工程计算中, 经常会遇到特征值和特征向量的计算, 如 : 机械结构或电磁振动中的固有值问题 ; 物理学中的各种临界值等