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迭代法,·当w<1时,称为低松弛(underrelaxation)迭代法,当w>1时,称为超松弛(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 时, 称为 低松弛 (under relaxation) 迭代法, • 当 ω > 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