《矩阵计算》课程教学资源（讲义）第六讲 线性方程组基本迭代法

第六讲线性方程组基本迭代法考虑线性方程组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 multi￾grid); 快速多极子法 (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