第六讲线性方程组基本迭代法考虑线性方程组Ar=b,AeRnxn,beRn目前，求解线性方程组的方法有：●直接法：PLU分解，LDLT分解，Cholesky分解等·送代法：基本迭代法：Jacobi,Gauss-Seidel,SOR,SSOR,Richardson,ADI等DKrylov子空间迭代法：CG,MINRES,GMRES,BiCGStab等·快速方法：D基于各种快速变换的求解方法，如FFT,DCT,DST等；代数多重网格法（Algebraicmultigrid)；快速多极子法(Fastmultipole),等等快速方法通常只适用于某些具有特殊结构的方程组实际应用中，这些方法常常结合使用，如混合(hybrid）方法，预处理（preconditioning）方法等直接法的优点是稳定可靠，能在有限步内得到近似解（如果不考虑误差，则得到精确解），而且所需存储量和运算量都是可知的.缺点是所需运算量约为O（n3），这对于大规模线性方程组来说是非常巨大的而且在实际应用中，很多问题中需要求解的大规模线性方程组都是稀疏的，如偏微分方程的有限差分/有限元离散，但直接法很难有效地利用问题的稀疏性来降低总运算量，而送代法则可以很好地利用问题的稀疏性，大大降低运算量从历史上看，最早的迭代法可以追溯到十九世纪Gauss，Jacobi，Seidel和Nekrasov等的工作[13，76].但是针对迭代法的系统研究主要还是在计算机出现以后，大约是从二十世纪五十年代开始在开始阶段主要研究的是基本选代法[6,30，139]（也称经典迭代法[57])，典型代表有Jacobi,GS,SOR,SSOR,ADI,Chebyshev迭代法等.在这期间，有两本非常有名的经典著作，一本是Varga的“MatrixIterativeAnalysis”（1962)[139]，另一本是Yong的“IterativeSolutionofLargeLinearSystems”(1971)[141].基本迭代法的收敛性有着非常完美的理论分析，但在实际使用中却存在着许多不足，比如收敛速度较慢，最优参数估计困难，等等从二十世纪七十年代中期开始，研究重点慢慢转向Krylov子空间选代法.事实上，早在1952年，Lanczos[85]】和Hestenes&Stiefel[68]就同时独立地提出了求解对称正定线性方程组的共轭梯度法（CG).对于一个n阶的线性方程组，如果不考虑舍人误差的影响，则共轭梯度法在n步后就一定会得到精确解。因此共轭梯度法一开始被认作是直接法。但在实际使用中发现，由于舍人误差的影响，迭代步数可能会超过n，特别是对于坏条件问题.而对于条件数较小的线性方程组，在给定精度下，所需送代步数则可能远远小于n，这使得共轭梯度法具有一定的吸引力.但由于种种原因[55]，共轭梯度法提出后并没有受到重视，在其出现后的近二十年里，主流方法仍然是Gauss198
第六讲 线性方程组基本迭代法 考虑线性方程组 Ax = b, A ∈ R n×n , b ∈ R n . 目前, 求解线性方程组的方法有: • 直接法: PLU 分解, LDLT 分解, Cholesky 分解等 • 迭代法: ▷ 基本迭代法: Jacobi, Gauss-Seidel, SOR, SSOR, Richardson, ADI 等 ▷ Krylov 子空间迭代法: CG, MINRES, GMRES, BiCGStab 等 • 快速方法: ▷ 基于各种快速变换的求解方法, 如 FFT, DCT, DST 等; 代数多重网格法 (Algebraic multigrid); 快速多极子法 (Fast multipole), 等等. b 快速方法通常只适用于某些具有特殊结构的方程组. b 实际应用中, 这些方法常常结合使用, 如混合 (hybrid) 方法, 预处理 (preconditioning) 方法等. 直接法的优点是稳定可靠, 能在有限步内得到近似解 (如果不考虑误差, 则得到精确解), 而且 所需存储量和运算量都是可知的. 缺点是所需运算量约为 O(n 3 ), 这对于大规模线性方程组来说 是非常巨大的. 而且在实际应用中, 很多问题中需要求解的大规模线性方程组都是稀疏的, 如偏微 分方程的有限差分/有限元离散, 但直接法很难有效地利用问题的稀疏性来降低总运算量, 而迭代 法则可以很好地利用问题的稀疏性, 大大降低运算量. 从历史上看, 最早的迭代法可以追溯到十九世纪 Gauss, Jacobi, Seidel 和 Nekrasov 等的工作 [13, 76]. 但是针对迭代法的系统研究主要还是在计算机出现以后, 大约是从二十世纪五十年代开 始. 在开始阶段主要研究的是基本迭代法 [6, 30, 139] (也称经典迭代法 [57]), 典型代表有 Jacobi, GS, SOR, SSOR, ADI, Chebyshev 迭代法等. 在这期间, 有两本非常有名的经典著作, 一本是 Varga 的 “Matrix Iterative Analysis” (1962) [139], 另一本是 Yong 的 “Iterative Solution of Large Linear Systems” (1971) [141]. 基本迭代法的收敛性有着非常完美的理论分析, 但在实际使用中却存在着许多不足, 比如收敛速度较慢, 最优参数估计困难, 等等. 从二十世纪七十年代中期开始, 研究重点慢慢转向 Krylov 子空间迭代法. 事实上, 早在 1952 年, Lanczos [85] 和 Hestenes & Stiefel [68] 就同时独立地提出了求解对称正定线性方程组的共轭梯 度法 (CG). 对于一个 n 阶的线性方程组, 如果不考虑舍入误差的影响, 则共轭梯度法在 n 步后就 一定会得到精确解. 因此共轭梯度法一开始被认作是直接法. 但在实际使用中发现, 由于舍入误 差的影响, 迭代步数可能会超过 n, 特别是对于坏条件问题. 而对于条件数较小的线性方程组, 在 给定精度下, 所需迭代步数则可能远远小于 n, 这使得共轭梯度法具有一定的吸引力. 但由于种种 原因 [55], 共轭梯度法提出后并没有受到重视, 在其出现后的近二十年里, 主流方法仍然是 Gauss 198
+199.消去法，SOR选代法和Chebyshev迭代法1971年，Reid[106]指出，对于好条件的大规模稀疏线性方程组，共轭梯度法能在很少的迭代步数内得到一个很好的近似解（事实上，Engeli等[41]在1959年就发现了该现象，但并没有引起关注).特别是预处理方法的引入[94]，大幅提升了共轭梯度法的收敛速度，这极大地促发了大家对共轭梯度法的研究兴趣，包括各种改进和推广，如求解对称不定线性方程组的MINRES迭代法和SYMMLQ送代法[99],求解一般线性方程组的GMRES送代法[111],QMR送代法[48],BiCGSTAB送代法[129]，等等.目前，带预处理的Krylov子空间迭代法已成为求解大规模稀疏线性方程组的主流方法本讲介绍常用的基本迭代法，关于Krylov子空间选代法，我们将在下一讲介绍.关于线性方程组基本选代法的相关参考资料G.H.Goluband C.F.VanLoan,MatrixComputations,2013.[57]R.S.Varga, MatrixIterativeAnalysis, 2nd edition,2000.[139]D.M.Young,IterativeSolutionofLargeLinearSystems,1971.[14]]O. Axelsson, Iterative Solution Methods, 1994. [] R. Barrett, et.al, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods,1994. [11]徐树方，矩阵计算的理论与方法，1995.[150] Y. Saad and H. A. van der Vorst, Iterative solution oflinear systems in the 20th century, 2000. [114]随着矩阵规模的增大，直接法的运算量也随之快速增长.对于大规模的线性方程组，由于运算量太大，往往会采用送代法.当直接求解Ar=b比较困难时，我们可以求解一个比较容易求解的近似等价方程组Ma=b,其中M可以看作是A在某种意义下的近似.设M=b的解为(1).易知它与原方程组的解*=A-1b之间的差距满足A(r+- 2(1)) =b- Ar(1),如果r(1）已经满足精度要求，即非常接近真解*，则可以停止计算，否则需要修正。记△会*-(1),则△r满足方程AAa=b-Ar(1).但由于直接求解该方程组比较困难(与求解原方程组一样困难），因此我们还是通过求解近似方程组MAr(1) = b- Ar(1),得到一个修正量△z(1).于是修正后的近似解为r(2) = 2(1) + △r(1) = r(1) + M-1(b - Ar(1)如果z(2)已经满足精度要求，则停止计算，否则继续按以上的方式进行修正，即求解M△r(2）=b-Ar(2)得到修正量△r(2),然后加到r(2)上得到z(3)：2(3) = 2(2) + r(2) = r(2) + M-1(b - Ar(2)不断重复以上步骤，于是，我们就得到一个序列2(1), r(2),.*,2()
· 199 · 消去法, SOR 迭代法和 Chebyshev 迭代法. 1971 年, Reid [106] 指出, 对于好条件的大规模稀疏线性方程组, 共轭梯度法能在很少的迭代 步数内得到一个很好的近似解 (事实上, Engeli 等 [41] 在 1959 年就发现了该现象, 但并没有引起关 注). 特别是预处理方法的引入 [94], 大幅提升了共轭梯度法的收敛速度, 这极大地促发了大家对 共轭梯度法的研究兴趣, 包括各种改进和推广, 如求解对称不定线性方程组的 MINRES 迭代法和 SYMMLQ 迭代法 [99], 求解一般线性方程组的 GMRES 迭代法 [111], QMR 迭代法 [48], BiCGSTAB 迭代法 [129], 等等. 目前, 带预处理的 Krylov 子空间迭代法已成为求解大规模稀疏线性方程组的 主流方法. 本讲介绍常用的基本迭代法, 关于 Krylov 子空间迭代法, 我们将在下一讲介绍. 关于线性方程组基本迭代法的相关参考资料 ▶ G. H. Golub and C. F. Van Loan, Matrix Computations, 2013. [57] ▶ R. S. Varga, Matrix Iterative Analysis, 2nd edition, 2000. [139] ▶ D. M. Young, Iterative Solution of Large Linear Systems, 1971. [141] ▶ O. Axelsson, Iterative Solution Methods, 1994. [6] ▶ R. Barrett, et.al, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 1994. [11] ▶ 徐树方, 矩阵计算的理论与方法, 1995. [150] ▷ Y. Saad and H. A. van der Vorst, Iterative solution of linear systems in the 20th century, 2000. [114] 随着矩阵规模的增大, 直接法的运算量也随之快速增长. 对于大规模的线性方程组, 由于运算 量太大, 往往会采用迭代法. 当直接求解 Ax = b 比较困难时, 我们可以求解一个比较容易求解的近似等价方程组 Mx = b, 其中 M 可以看作是 A 在某种意义下的近似. 设 Mx = b 的解为 x (1) . 易知它与原方程组的解 x∗ = A−1 b 之间的差距满足 A x∗ − x (1) = b − Ax(1) . 如果 x (1) 已经满足精度要求, 即非常接近真解 x∗, 则可以停止计算, 否则需要修正. 记 ∆x ≜ x∗ − x (1) , 则 ∆x 满足方程 A∆x = b − Ax(1) . 但由于直接求解该方程组比较困难 (与求解原方程 组一样困难), 因此我们还是通过求解近似方程组 M∆x (1) = b − Ax(1) , 得到一个修正量 ∆x (1) . 于是修正后的近似解为 x (2) = x (1) + ∆x (1) = x (1) + M−1 (b − Ax(1)). 如果 x (2) 已经满足精度要求, 则停止计算, 否则继续按以上的方式进行修正, 即求解 M∆x (2) = b − Ax(2) 得到修正量 ∆x (2) , 然后加到 x (2) 上得到 x (3): x (3) = x (2) + ∆x (2) = x (2) + M−1 (b − Ax(2)). 不断重复以上步骤, 于是, 我们就得到一个序列 x (1), x(2), . . . , , x(k) , . . .
200.第六讲线性方程组基本选代法满足以下递推关系r(k+1) = 2(k) + M-1(b - Ar(k), k = 1, 2, ...这就构成了一个选代法注记通常，构造一个好的迭代法，需要考虑以下两点：(1）以M为系数矩阵的线性方程组要比原线性方程组更容易求解；(2）M应该是A的一个很好的近似，而且选代序列(())要收敛。常用的基本选代法主要包括：·Richardson迭代法·Jacobi送代法·Gauss-Seidel (G-S)迭代法·超松弛（SOR，SuccessiveOver-Relaxation）送代法·对称超松弛（SSOR,SymmetricSOR）送代法·加速超松弛（AOR，AcceleratedOver-Relaxation）送代法·Chebyshev（加速）送代法·交替方向（ADI,AlternatingDirectionImplicit)送代法
· 200 · 第六讲 线性方程组基本迭代法 满足以下递推关系 x (k+1) = x (k) + M−1 (b − Ax(k) ), k = 1, 2, . . . . 这就构成了一个迭代法. 注记 通常, 构造一个好的迭代法, 需要考虑以下两点: (1) 以 M 为系数矩阵的线性方程组要比原线性方程组更容易求解; (2) M 应该是 A 的一个很好的近似, 而且迭代序列 {x (k)} 要收敛. 常用的基本迭代法主要包括: • Richardson 迭代法 • Jacobi 迭代法 • Gauss-Seidel (G-S) 迭代法 • 超松弛 (SOR, Successive Over-Relaxation) 迭代法 • 对称超松弛 (SSOR, Symmetric SOR) 迭代法 • 加速超松弛 (AOR, Accelerated Over-Relaxation) 迭代法 • Chebyshev (加速) 迭代法 • 交替方向 (ADI, Alternating Direction Implicit) 迭代法
6.1矩阵分裂与迭代法-201.6.1矩阵分裂与迭代法首先给出矩阵分裂的定义定义6.1（矩阵分裂MatrixSplitting）设AEIRnxn非奇异，称(6.1)A=M-N为A的一个矩阵分裂，其中M 非奇异考虑线性方程组Ar=b,(6.2)其中AEIRnxn非奇异.选代法的基本思想：给定一个选代初始值r(0),通过一定的迭代格式生成一个选代序列（(k)=0,使得lim 2()=会 A-1b.给定一个矩阵分裂（6.1),则原方程组（6.2)就等价于MT=Na+6.于是我们就可以构造出以下的迭代格式Mr(k+1) = Nr(k) +b, k = 0, 1,2,..,或2(k+1) = M-1Nα(k) + M-1b ≤ Gr(k) + g, k = 0, 1,2,..(6.3)其中G会M-1N称为选代矩阵.这就是基于矩阵分裂（6.1)的选代法.易知，选取不同的M，就可以构造出不同的迭代法6.1.1 Jacobi送代法将矩阵A写成A-D-L-U其中D为A的对角线部分，-L和-U分别为A的严格下三角和严格上三角部分在矩阵分裂A=M-N中取M=D,N=L+U,则可得Jacobi选代法2(k+1) = D-1(L+U)(k) + D-1b, k = 0, 1,2,..(6.4)对应的迭代矩阵为G) = D-1(L + U).写成分量形式即为n1(k)(k+1)N.aijr=1.2.....nCaiij=1.jti
6.1 矩阵分裂与迭代法 · 201 · 6.1 矩阵分裂与迭代法 首先给出 矩阵分裂 的定义. 定义 6.1 (矩阵分裂 Matrix Splitting) 设 A ∈ R n×n 非奇异, 称 A = M − N (6.1) 为 A 的一个矩阵分裂, 其中 M 非奇异. 考虑线性方程组 Ax = b, (6.2) 其中 A ∈ R n×n 非奇异. 迭代法的基本思想: 给定一个迭代初始值 x (0) , 通过一定的迭代格式生成 一个迭代序列 {x (k)}∞ k=0, 使得 lim k→∞ x (k) = x∗ ≜ A −1 b. 给定一个矩阵分裂 (6.1), 则原方程组 (6.2) 就等价于 Mx = Nx + b. 于是我们就可以构造出以下 的迭代格式 Mx(k+1) = Nx(k) + b, k = 0, 1, 2, . . . , 或 x (k+1) = M−1Nx(k) + M−1 b ≜ Gx(k) + g, k = 0, 1, 2, . . . , (6.3) 其中 G ≜ M−1N 称为 迭代矩阵. 这就是基于矩阵分裂 (6.1) 的迭代法. 易知, 选取不同的 M, 就 可以构造出不同的迭代法. 6.1.1 Jacobi 迭代法 将矩阵 A 写成 A = D − L − U, 其中 D 为 A 的对角线部分, −L 和 −U 分别为 A 的严格下三角和严格上三角部分. 在矩阵分裂 A = M − N 中取 M = D, N = L + U, 则可得 Jacobi 迭代法: x (k+1) = D−1 (L + U)x (k) + D−1 b, k = 0, 1, 2, . . . . (6.4) 对应的迭代矩阵为 GJ = D−1 (L + U). 写成分量形式即为 x (k+1) i = 1 aii bi − Xn j=1,j̸=i aijx (k) j , i = 1, 2, . . . , n
-202.第六讲线性方程组基本选代法白由于Jacobi选代法中a(k+1)的更新顺序与i无关，即可以按顺序i=1,2..,n计算，也可以按顺序i=n,n－1,,2,1计算，或者乱序计算.因此Jacobi送代法非常适合并行计算.算法6.1.Jacobi选代法1: Given an initial guess 2(0)2: while not converge dofori=ltondo3:1n2(+1) (k)4:b;aiiaij=ljtiendfor5:6: end while我们有时也将Jacobi送代格式写为2(k+1) = 2() + D-1(b - Ar(k) = () + D-1rk, k = 0, 1,2, ..其中rkb-Az()是k次选代后的残量.这表明，2(k+1)是通过对r(k)做一个修正得到的6.1.2Gauss-Seidel迭代法在分裂A=M-N中取M=D-L,N=U,即可得Gauss-Seidel(G-S)迭代法：r(k+1) = (D - L)-1Ua(k) + (D - L)-1b.(6.5)对应的迭代矩阵为GGs = (D - L)-1U.将G-S选代法改写为Dr(k+1) = Lr(h+1) + Ua() + b,即可得分量形式i-1n12(+),(k+1)(k)T>b1,2,...,n.aijraijtaij=1j=i+1算法6.2.Gauss-Seidel选代法1: Given an initial guess (0)2: while not converge dofor i= l to n do3:i-11(k+1)(k+1)4:bZaijjCiaijtai5end for6: end while
· 202 · 第六讲 线性方程组基本迭代法 b 由于 Jacobi 迭代法中 x (k+1) i 的更新顺序与 i 无关, 即可以按顺序 i = 1, 2, . . . , n 计算, 也可 以按顺序 i = n, n − 1, . . . , 2, 1 计算, 或者乱序计算. 因此 Jacobi 迭代法非常适合并行计算. 算法 6.1. Jacobi 迭代法 1: Given an initial guess x (0) 2: while not converge do 3: for i = 1 to n do 4: x (k+1) i = 1 aii bi − Pn j=1,j̸=i aijx (k) j ! 5: end for 6: end while 我们有时也将 Jacobi 迭代格式写为 x (k+1) = x (k) + D−1 (b − Ax(k) ) = x (k) + D−1 rk, k = 0, 1, 2, . . . , 其中 rk ≜ b − Ax(k) 是 k 次迭代后的残量. 这表明, x (k+1) 是通过对 x (k) 做一个修正得到的. 6.1.2 Gauss-Seidel 迭代法 在分裂 A = M − N 中取 M = D − L, N = U, 即可得 Gauss-Seidel (G-S) 迭代法: x (k+1) = (D − L) −1Ux(k) + (D − L) −1 b. (6.5) 对应的迭代矩阵为 GGS = (D − L) −1U. 将 G-S 迭代法改写为 Dx(k+1) = Lx(k+1) + Ux(k) + b, 即可得分量形式 x (k+1) i = 1 aii bi − X i−1 j=1 aijx (k+1) j − Xn j=i+1 aijx (k) j , i = 1, 2, . . . , n. 算法 6.2. Gauss-Seidel 迭代法 1: Given an initial guess x (0) 2: while not converge do 3: for i = 1 to n do 4: x (k+1) i = 1 aii bi − i P−1 j=1 aijx (k+1) j − Pn j=i+1 aijx (k) j ! 5: end for 6: end while
6.1矩阵分裂与迭代法-203.G-S选代法的主要优点是充分利用了已经获得的最新数据但在G-S选代法中，未知量的更新必须按自然顺序进行，因此不适合并行计算6.1.3SOR送代法在G-S送代法的基础上，我们可以通过引入一个松弛参数w来加快收敛速度这就是SOR(SuccessiveOverrelaxation)方法[140].该方法的基本思想是将G-S迭代法中的第k十1步近似解与第K步近似解做一个加权平均，从而给出一个新的近似解，即(k+1) = (1 - w)() +wD-1 (L(K+1) + U(k) +b)(6.6)整理后即为(h+1) = (D -wL)-1(1 -w)D + wU)() +w(D -wL)-1b,其中w称为松弛参数(relaxationparameter)当w=1时,SOR即为G-S迭代法，·当w1时，称为超松弛（overrelaxation）送代法在大多数情况下，当w>1时会取得比较好的收敛效果SOR选代法曾经在很长一段时间内是科学计算中求解线性方程组的首选方法SOR的迭代矩阵为GsOR = (D -wL)-1 ((1 -w)D +wU),对应的矩阵分裂为11-WD+U.M=D-L.Nww由（6.6)可知SOR迭代法的分量形式为i~1w(++1) = (1 - w)() +(k+1)2kTaijraijjaii1j=i+1WiS(k+1)= 2() +Fbaiiai算法6.3.求解线性方程组的SOR选代法1: Given an initial guess r(0) and parameter w2:whilenot converge dofor i =l to n do3:
6.1 矩阵分裂与迭代法 · 203 · b G-S 迭代法的主要优点是充分利用了已经获得的最新数据. b 但在 G-S 迭代法中, 未知量的更新必须按自然顺序进行, 因此不适合并行计算. 6.1.3 SOR 迭代法 在 G-S 迭代法的基础上, 我们可以通过引入一个松弛参数 ω 来加快收敛速度. 这就是 SOR (Successive Overrelaxation) 方法 [140]. 该方法的基本思想是将 G-S 迭代法中的第 k + 1 步近似解与 第 k 步近似解做一个加权平均, 从而给出一个新的近似解, 即 x (k+1) = (1 − ω)x (k) + ωD−1 Lx(k+1) + Ux(k) + b . (6.6) 整理后即为 x (k+1) = (D − ωL) −1 (1 − ω)D + ωU x (k) + ω(D − ωL) −1 b, 其中 ω 称为 松弛参数 (relaxation parameter). • 当 ω = 1 时, SOR 即为 G-S 迭代法, • 当 ω 1 时, 称为 超松弛 (over relaxation) 迭代法. 在大多数情况下, 当 ω > 1 时会取得比较好的收敛效果. b SOR 迭代法曾经在很长一段时间内是科学计算中求解线性方程组的首选方法. SOR 的迭代矩阵为 GSOR = (D − ωL) −1 ((1 − ω)D + ωU), 对应的矩阵分裂为 M = 1 ω D − L, N = 1 − ω ω D + U. 由 (6.6) 可知 SOR 迭代法的分量形式为 x (k+1) i = (1 − ω)x (k) i + ω aii bi − X i−1 j=1 aijx (k+1) j − Xn j=i+1 aijx (k) j = x (k) i + ω aii bi − X i−1 j=1 aijx (k+1) j − Xn j=i aijx (k) j 算法 6.3. 求解线性方程组的 SOR 迭代法 1: Given an initial guess x (0) and parameter ω 2: while not converge do 3: for i = 1 to n do
-204.第六讲线性方程组基本选代法fa二n(++1) = (1 - w)(*) +(k+1)(k)ais4aiij=1j=i+1end for5:6: end whileSOR选代法最大的优点是引入了松弛参数w：通过选取适当的w就可以大大提高方法的收敛速度.但是SOR选代法最大的难点就是如何选取最优的参数6.1.4SSOR送代法将SOR送代法中的L和U相互交换位置，则可得送代格式2(+1) = (D -wU)-1 (1 - w)D +wL) 2(k) + w(D -wU)-1b将这个迭代格式与SOR相结合，就可以得到下面的两步选代法[r(h+) = (D -wL)-1[(1 - w)D +wU]r(k) + w(D -wL)-1b,((h+1) = (D - wU)-I[(1 - w)D + wL)r(+) + w(D - wU)-1b.这就是对称超松弛（SSOR）送代法，相当于将L与U同等看待，交替做两次SOR迭代法消去中间选代向量(k+),可得2(+1) = GssORr(k) + g,其中迭代矩阵GssOR = (D -wU)-1[(1 -w)D +wL](D -wL)-1[(1 -w)D +wU)对应的矩阵分裂为1[D-w(L+U) +w?LD-1U]M =w(2 -w) 1(D -wL)D-1(D - wU),w(2 -w)1N:w(2-][1 - w)D +wL] -[1 -)D +wu] 对于某些特殊问题,SOR选代法不收敛，但仍然可能构造出收敛的SSOR选代法一般来说，SOR收敛速度要快于SSOR，但SOR的收敛速度对参数更敏感算法6.4.SSOR送代法1: Given an initial guess (0) and parameter w2: while not converge dofori=l tondo3:
· 204 · 第六讲 线性方程组基本迭代法 4: x (k+1) i = (1 − ω)x (k) i + ω aii bi − i P−1 j=1 aijx (k+1) j − Pn j=i+1 aijx (k) j ! 5: end for 6: end while b SOR 迭代法最大的优点是引入了松弛参数 ω: 通过选取适当的 ω 就可以大大提高方法的收 敛速度. 但是 SOR 迭代法最大的难点就是如何选取最优的参数. 6.1.4 SSOR 迭代法 将 SOR 迭代法中的 L 和 U 相互交换位置, 则可得迭代格式 x (k+1) = (D − ωU) −1 ((1 − ω)D + ωL) x (k) + ω(D − ωU) −1 b. 将这个迭代格式与 SOR 相结合, 就可以得到下面的两步迭代法 x (k+ 1 2 ) = (D − ωL) −1 (1 − ω)D + ωU x (k) + ω(D − ωL) −1 b, x (k+1) = (D − ωU) −1 (1 − ω)D + ωL x (k+ 1 2 ) + ω(D − ωU) −1 b. 这就是 对称超松弛 (SSOR ) 迭代法, 相当于将 L 与 U 同等看待, 交替做两次 SOR 迭代法. 消去中间迭代向量 x (k+ 1 2 ) , 可得 x (k+1) = GSSORx (k) + g, 其中迭代矩阵 GSSOR = (D − ωU) −1 (1 − ω)D + ωL (D − ωL) −1 (1 − ω)D + ωU . 对应的矩阵分裂为 M = 1 ω(2 − ω) D − ω(L + U) + ω 2LD−1U = 1 ω(2 − ω) (D − ωL)D−1 (D − ωU), N = 1 ω(2 − ω) (1 − ω)D + ωL D−1 (1 − ω)D + ωU . b 对于某些特殊问题, SOR 迭代法不收敛, 但仍然可能构造出收敛的 SSOR 迭代法. b 一般来说, SOR 收敛速度要快于 SSOR, 但 SOR 的收敛速度对参数 ω 更敏感. 算法 6.4. SSOR 迭代法 1: Given an initial guess v (0) and parameter ω 2: while not converge do 3: for i = 1 to n do
6.1矩阵分裂与迭代法-205.(k+)2(+) = (1 - w)a() +w4:biCaijrj4aij=i+1end for5:for i = n to 1 do6:INWw2(+1) = (1 - w)a(+),(k+)(k+1)bi7:aijrjaijtjaij=i+18:end for9: end whileAOR选代法*6.1.5Hadjidimos于1978年提出了加速超松弛（AcceleratedOver-Relaxation，AOR）送代法，送代矩阵为GAOR = (D -L)-1[(1 -w)D + (w -)L +wU]其中和w为松弛参数.对应的矩阵分裂为I(D-L), N=M =[(1-w)D + (w-)L+wU]易知：(1）当=w时，AOR即为SOR送代法；(2)当=w=1时,AOR即为G-S送代法：(3）当=0,w=1时,AOR即为Jacobi送代法AOR选代法中含有两个参数.因此在理论上，通过选取合适的参数，AOR选代法会收敛得更快.但也是因为含有两个参数，使得参数的选取变得更加困难，因此较少使用凸与SSOR类似，我们也可以定义SAOR选代法6.1.6Richardson迭代法Richardson选代法是一类形式非常简单的算法，其选代格式为r(k+1) = 2() +w(b - Ar(k),k=0.1.2.它可以看作是基于以下矩阵分裂的选代法1N:M =Cw对应的迭代矩阵为GR = I - wA.在每次选代时可以选取不同的参数，即(#+1) = 2() + wk(b - Ar(k), k= 0, 1,2,
6.1 矩阵分裂与迭代法 · 205 · 4: x (k+ 1 2 ) i = (1 − ω)x (k) i + ω aii bi − i P−1 j=1 aijx (k+ 1 2 ) j − Pn j=i+1 aijx (k) j ! 5: end for 6: for i = n to 1 do 7: x (k+1) i = (1 − ω)x (k+ 1 2 ) i + ω aii bi − i P−1 j=1 aijx (k+ 1 2 ) j − Pn j=i+1 aijx (k+1) j ! 8: end for 9: end while 6.1.5 AOR 迭代法 * Hadjidimos 于 1978 年提出了 加速超松弛 (Accelerated Over-Relaxation, AOR) 迭代法, 迭代矩阵 为 GAOR = (D − γL) −1 (1 − ω)D + (ω − γ)L + ωU , 其中 γ 和 ω 为松弛参数. 对应的矩阵分裂为 M = 1 ω (D − γL), N = 1 ω [(1 − ω)D + (ω − γ)L + ωU]. 易知: (1) 当 γ = ω 时, AOR 即为 SOR 迭代法; (2) 当 γ = ω = 1 时, AOR 即为 G-S 迭代法; (3) 当 γ = 0, ω = 1 时, AOR 即为 Jacobi 迭代法. b AOR 迭代法中含有两个参数. 因此在理论上, 通过选取合适的参数, AOR 迭代法会收敛得 更快. 但也是因为含有两个参数, 使得参数的选取变得更加困难, 因此较少使用. b 与 SSOR 类似, 我们也可以定义 SAOR 迭代法. 6.1.6 Richardson 迭代法 Richardson 迭代法是一类形式非常简单的算法, 其迭代格式为 x (k+1) = x (k) + ω(b − Ax(k) ), k = 0, 1, 2, . . . . 它可以看作是基于以下矩阵分裂的迭代法: M = 1 ω I, N = 1 ω I − A. 对应的迭代矩阵为 GR = I − ωA. b 在每次迭代时可以选取不同的参数, 即 x (k+1) = x (k) + ωk(b − Ax(k) ), k = 0, 1, 2, . . . .
206.第六讲线性方程组基本选代法6.1.7分块送代法如果A的对角线中出现零，则Jacobi,G-S,SOR等方法就不再有定义.此时我们可以采用分块送代格式.另外，分块选代法也能提升算法的计算效率An将A写成如下的分块形式（如右图所示）：A22A11A12A1pA21A22..A2pA=Ap1Ap2.. App]Ap则相应的分块Jacobi，分块G-S和分块SOR选代法分别定义为·分块Jacobi送代：之Aua(k+1)Aijafe,= bii=1,2,...,p.j=1.jti·分块Gauss-seidel送代：i-1Pr(k+1)Agafe,Aja(+1) = b; ZEAgai=1,2,....p-j=1j=i+1·分块SOR送代：-a(h+1) = (1 - w)r(k),(k+1)+wAAijaAijcbej=1j=i+1i=1,2,...,p.分块选代法采用更多的3级BLAS运算，因此更有利于发挥现代计算机的性能
· 206 · 第六讲 线性方程组基本迭代法 6.1.7 分块迭代法 如果 A 的对角线中出现零, 则 Jacobi, G-S, SOR 等方法就不再有定义. 此时我们可以采用分块 迭代格式. 另外, 分块迭代法也能提升算法的计算效率. 将 A 写成如下的分块形式 (如右图所示): A = A11 A12 · · · A1p A21 A22 · · · A2p . . . . . . . . . . . . Ap1 Ap2 · · · App , 则相应的分块 Jacobi, 分块 G-S 和分块 SOR 迭代法分别定义为: • 分块 Jacobi 迭代: Aiix (k+1) i = bi − X p j=1,j̸=i Aijx (k) j , i = 1, 2, . . . , p. • 分块 Gauss-seidel 迭代: Aiix (k+1) i = bi − X i−1 j=1 Aijx (k+1) j − X p j=i+1 Aijx (k) j , i = 1, 2, . . . , p. • 分块 SOR 迭代: x (k+1) i = (1 − ω)x (k) i + ωA−1 ii bi − X i−1 j=1 Aijx (k+1) j − X p j=i+1 Aijx (k) j , i = 1, 2, . . . , p. b 分块迭代法采用更多的 3 级 BLAS 运算, 因此更有利于发挥现代计算机的性能
6.2收敛性分析-207.6.2收敛性分析6.2.1基本概念首先给出向量序列收敛的定义定义 6.2 (向量序列的收敛)设[r(k)) =, 是 R" (或 Cn)中的一个向量序列. 如果存在向量 =[E1,2,...,an]TeRn(或Cn)使得Jlim a(k) = ri, i=1,2,...,n,-其中r()表示2(h)的第i个分量，则称(k)) (按分量)收敛到r,即是2(k)的极限,记为lim 2(k) = 2.0类似地，我们可以给出矩阵序列收敛的定义定义6.3(矩阵序列的收敛)设【A()=[a(]是Rmxn（或Cmxn）中的一个矩阵序列.如果存在矩阵 A=[ai] e Rmxn(或 Cmxn)使得(k)limaisi=1,2,...,m,j=1,2,...,n,=aj，k则称A()收敛到A,即A是A(K)的极限,记为lim A(k) = A.k-→00关于向量序列和矩阵序列的收敛性，我们有下面的基本判别方法[144]（天）C Rmxn (或定理6.1 设向量序列 [r(k)%= C Rn (或 Cn),矩阵序列[A(k)=a%Cmxn),则[lim 2(k)=→,lim ll2(k)-2l= 0, 其中 I 为任一向量范数;(1)lim A(k)=Alim IIA(k)-All=0,其中 II·II为任一矩阵范数;(2) limA()=0limA()=0,VrERn(或Cn)(3)-x42由定理6.1和引理1.39，我们可以立即得到下面的结论，(板书）定理6.2 设矩阵 A Rnxn(或 Cnxn),则 lim Ak=0 当且仅当 p(A)<1.证明.充分性：设p(A)<1,则由引理1.39可知,存在某个矩阵范数I·II使得Al<1.因此由0AA→0(→)可得lim IIA'll= 0.根据定理6.1,我们有limAk=0.必要性：设limAk=0,则由定理6.1可知lim II/Al = 0
6.2 收敛性分析 · 207 · 6.2 收敛性分析 6.2.1 基本概念 首先给出向量序列收敛的定义. 定义 6.2 (向量序列的收敛) 设 x (k) ∞ k=1 是 R n (或 C n ) 中的一个向量序列. 如果存在向量 x = [x1, x2, . . . , xn] T ∈ R n (或 C n ) 使得 lim k→∞ x (k) i = xi , i = 1, 2, . . . , n, 其中 x (k) i 表示 x (k) 的第 i 个分量, 则称 x (k) (按分量) 收敛到 x, 即 x 是 x (k) 的极限, 记为 lim k→∞ x (k) = x. 类似地, 我们可以给出矩阵序列收敛的定义. 定义 6.3 (矩阵序列的收敛) 设 n A(k) = h a (k) ij io∞ k=0 是 R m×n (或 C m×n ) 中的一个矩阵序列. 如 果存在矩阵 A = [aij ] ∈ R m×n (或 C m×n ) 使得 lim k→∞ a (k) ij = aij , i = 1, 2, . . . , m, j = 1, 2, . . . , n, 则称 A(k) 收敛到 A, 即 A 是 A(k) 的极限, 记为 lim k→∞ A (k) = A. 关于向量序列和矩阵序列的收敛性, 我们有下面的基本判别方法 [144]. 定理 6.1 设向量序列 {x (k)}∞ k=0 ⊂ R n (或 C n ), 矩阵序列 n A(k) = h a (k) ij io∞ k=0 ⊂ R m×n (或 C m×n ), 则 (1) lim k→∞ x (k) = x ⇐⇒ lim k→∞ ∥x (k) − x∥ = 0, 其中 ∥ · ∥ 为任一向量范数; (2) lim k→∞ A(k) = A ⇐⇒ lim k→∞ ∥A(k) − A∥ = 0, 其中 ∥ · ∥ 为任一矩阵范数; (3) lim k→∞ A(k) = 0 ⇐⇒ lim k→∞ A(k)x = 0, ∀ x ∈ R n (或 C n ). 由定理 6.1 和引理 1.39, 我们可以立即得到下面的结论. 定理 6.2 设矩阵 A ∈ R n×n (或 C n×n ), 则 lim k→∞ Ak = 0 当且仅当 ρ(A) < 1. (板书) 证明. 充分性: 设 ρ(A) < 1, 则由引理 1.39 可知, 存在某个矩阵范数 ∥ · ∥ 使得 ∥A∥ < 1. 因此由 0 ≤ ∥Ak∥ ≤ ∥A∥ k → 0 (k → ∞) 可得 lim k→∞ ∥A k ∥ = 0. 根据定理 6.1, 我们有 lim k→∞ Ak = 0. 必要性: 设 lim k→∞ Ak = 0, 则由定理 6.1 可知 lim k→∞ ∥A k ∥ = 0