《数值分析（原书第2版）》—— 2.5　迭代方法

本节书摘来自华章出版社《数值分析（原书第2版）》一 书中的第2章，第2.5节，作者：（美）Timothy Sauer，更多章节内容可以访问云栖社区“华章计算机”公众号查看。

2.5　迭代方法

高斯消去是一个包含O（n3）次有穷序列的浮点操作，并得到一个解.正因为如此，高斯消去被称为求解线性方程组的直接方法.理论上，直接方法在有穷步里可以得到精确解（当然在有穷精度的计算机上进行计算，只能得到近似解.正如我们在前面的章节中所见，精度的损失可以通过条件数进行度量.）直接方法和第1章讲过的根求解方法完全不同，这是由于第1章中的方法是迭代求解.
迭代方法也可以用来求解线性方程组系统.和不动点迭代类似，方法从初始估计开始，然后在每步中不断精化估计，最后收敛到解向量.

2.5.1　雅可比方法

雅可比（Jacobi）方法是方程组系统中的一种形式的不动点迭代.在FPI中，第一步是重写方程，进而求解未知量.雅可比方法以如下标准方式进行该重写步骤:求解第i个方程得到第i个未知变量.然后使用不动点迭代，从初始估计开始，进行迭代.
例2.19　对系统3u+v=5，u+2v=5使用雅可比方法.
首先求解第一个方程得到u，求解第二个方程得到v.我们将使用初始估计（u0，v0）=（0，0）.我们有u=5-v3
v=5-u2(2.35)两个方程进行迭代:u0
v0=0
0
u1
v1=5-v03
5-u02=5-03
5-02=53
52
u2
v2=5-v13
5-u12=5-5/23
5-5/32=56
53
u3
v3=5-5/33
5-5/62=109
2512(2.36)雅可比方法更多的步骤显示该方法将收敛，并得到解［1，2］.
现在假设方程以相反的顺序给出.
例2.20　对系统u+2v=5，3u+v=5使用雅可比方法.
首先求解第一个方程得到u，求解第二个方程得到v.u=5-2v
v=5-3u(2.37)
106两个方程如前进行迭代，但是结果全然不同:u0
v0=0
0
u1
v1=5-2v0
5-3u0=5
5
u2
v2=5-2v1
5-3u1=-5
-10
u3
v3=5-2(-10)
5-3(-5)=25
20(2.38)在这种情况下，迭代发散，雅可比方法失效了.
由于雅可比方法并不总能成功得到解，了解该方法可以工作的条件有助于我们的理解.下面的定义给出一个重要的条件.
定义2.9　n×n的矩阵A=(aij)是严格对角占优矩阵，要求对于每个1≤i≤n，aii>∑j≠iaij.换句话说，对角线元素在对应行占优，其对应值在数量上比该行其他元素的和还要大.
定理2.10　如果n×n矩阵A是严格对角占优矩阵，则（1）A是非奇异矩阵，（2）对于所有向量b和初始估计，对Ax=b应用雅可比方法都会收敛到（唯一）解.
定理2.10指出，如果A是严格对角占优矩阵，则对方程Ax=b应用雅可比方法，将对于每一个初始估计都收敛.该事实的证明在2.5.3节中给出.在例2.19中，开始时的系数矩阵是A=31
12这是一个严格对角占优矩阵，由于3>1，2>1.在这种情况下保证收敛.而在例2.20中，雅可比方法应用在矩阵A=12
31该矩阵不是对角占优矩阵，对于收敛的保证也不存在.注意，严格对角占优仅仅是一个充分条件.不满足对角占优时，雅可比方法依然可能收敛.
例2.21　确定矩阵A=31-1
2-52
168,B=326
181
92-2是否严格对角占优.
矩阵A是对角占优，这是因为3>1+-1，-5>2+2，8>1+6.B不是对角占优，因为，例如3>2+6不成立.但是如果B的第一行和第三行交换，B是严格对角占优，并且雅可比方法保证收敛.107
雅可比方法是不动点迭代的一种形式，令D表示A的主对角线矩阵，L表示矩阵A的下三角矩阵（主对角线以下的元素），U表示上三角矩阵（主对角线以上的元素），则A=L+D+U，求解的方程变为Lx+Dx+Ux=b.注意，这里对于L和U的使用和LU分解中不同，当前L和U的主对角线元素都是零.方程组Ax=b可以写成如下不动点迭代形式:Ax=b
(D+L+U)x=b
Dx=b-(L+U)x
x=D-1(b-(L+U)x)(2.39)由于D是对角线矩阵，它的逆矩阵中主对角线元素是A的对角线元素的倒数.雅可比方法就是（2.39）中的不动点迭代:
雅可比方法x0=初始向量
xk+1=D-1(b-(L+U)xk)，k=0，1，2，…(2.40)对于例2.19，31
12=u
v=5
5（2.40）的不动点迭代如下，其中xk=uk
vkuk+1
vk+1=D-1(b-(L+U)xk)
=1/30
01/25
5-01
10uk
vk
=(5-vk)/3
(5-uk)/2这和我们原来计算的版本一致.

2.5.2　高斯塞德尔方法和SOR

和雅可比方法紧密相关的迭代方法是高斯塞德尔（Gauss-Seidel）方法.高斯塞德尔方法和雅可比方法的唯一差异是在前者中，最近更新的未知变量的值在每一步中都使用，即使更新发生在当前步骤.回到例2.19，我们发现高斯塞德尔方法如下:u0
v0=0
0
u1
v1=5-v03
5-u12=5-03
5-5/32=53
53
u2
v2=5-v13
5-u22=5-5/33
5-10/92=109
3518u3
v3=5-v23
5-u32=5-35/183
5-55/542=5554
215108(2.41)

108注意,高斯塞德尔方法和雅可比方法的差异:对于v1的定义中使用u1，而不是u0.我们看到该方法如同雅可比方法一样趋近于解［1，2］，但是在使用相同步骤的情况下，该近似更加精确.如果收敛，高斯塞德尔方法常常比雅可比方法收敛更快.定理2.11证实高斯塞德尔方法和雅可比方法一样，只要系数矩阵是严格对角占优矩阵，该方法就收敛.
高斯塞德尔方法可以写成矩阵形式，并和不动点迭代一致，其中将方程（L+D+U）x=b拆为(L+D)xk+1=-Uxk+b注意，通过把矩阵A的下三角部分放在方程左侧，迭代使用了最近确定的元素xk+1，将方程进行重新组织，得到高斯塞德尔方法.
高斯塞德尔方法x0=初始向量
xk+1=D-1(b-Uxk-Lxk+1),k=0,1,2,…　　例2.22　对如下系统应用高斯塞德尔方法：31-1
241
-125u
v
w=4
1
1高斯塞德尔迭代如下uk+1=4-vk+wk3
vk+1=1-2uk+1-wk4
wk+1=1+uk+1-2vk+15初始x0=［u0，v0，w0］=［0，0，0］，我们计算得到u1
v1
w1=4-0-03=43
1-8/3-04=-512
1+4/3+5/65=1930≈1.3333
-0.4167
0.6333
u2
v2
w2=10160
-34
251300≈1.6833
-0.7500
0.8367系统严格对角占优，因而迭代会收敛到解［2，-1，1］.
连续过松弛（SOR）方法使用高斯塞德尔的求解方向，并使用过松弛109以加快收敛速度.令ω是一个实数，并将新的估计中的每个元素xk+1定义为ω乘上高斯塞德尔公式和1-ω乘上当前估计xk的平均.数字ω被称为松弛参数，而当ω>1时被称为过松弛.
例2.23　对例2.22的系统使用SOR方法，ω=1.25.
连续过松弛得到uk+1=(1-ω)uk+ω4-vk+wk3
vk+1=(1-ω)vk+w1-2uk+1-wk4
wk+1=(1-ω)wk+ω1+uk+1-2vk+15从［u0，v0，w0］=［0，0，0］开始，我们计算得到u1
v1
w1≈1.6667
-0.7292
1.0312
u2
v2
w2≈1.9835
-1.0672
1.0216在这个例子中，SOR迭代比雅可比和高斯塞德尔方法收敛更快，并收敛到［2，-1，1］.
正如雅可比和高斯塞德尔方法，另一种导出SOR的方法是将该系统看做不动点迭代.问题Ax=b可以写成（L+D+U）x=b，乘上ω并重新组织方程，(ωL+ωD+ωU)x=ωb
(ωL+D)x=ωb-ωUx+(1-ω)Dx
x=(ωL+D)-1［(1-ω)Dx-ωUx］+ω(D+ωL)-1b连续过松弛（SOR）x0=初始向量
xk+1=(ωL+D)-1［(1-ω)Dxk-ωUxk］+ω(D+ωL)-1b,k=0,1,2,…当ω=1时，SOR就是高斯塞德尔方法，在连续欠松弛方法中，参数ω还可以比1小.
例2.24　比较雅可比、高斯塞德尔和SOR求解6个方程6个未知变量的方程组:3-100012
-13-10120
0-13-100
00-13-10
0120-13-1
12000-13u1
u2
u3
u4
u5
u6=52
32
1
1
32
52(2.42)

110　　该问题的解是x=［1，1，1，1，1，1］.使用三种方法，经过6步迭代得到的解向量x6，如下表所示:
雅可比高斯塞德尔SOR雅可比高斯塞德尔SOR雅可比高斯塞德尔SOR0.98790.99500.99890.98460.99460.99930.96740.99691.00040.96740.99961.00090.98461.00161.00090.98791.00131.0004
过松弛方法的参数ω被设为1.1.SOR在这个例子上看起来性能最好.
图2.3比较了例2.24中，对于不同ω，6步迭代后的无穷范数的误差，尽管没有通用的理论描述ω的最优选择，但是在当前例子中很显然有图2.3　例2.24中，SOR 6步迭代后关于过松弛参数ω的无穷范数误差.高斯塞德尔对应ω=1.最小误差发生在ω≈1.13处最好的选择.参看Ortega［1972］中关于在一些常见的特例中最优ω的讨论.

2.5.3　迭代方法的收敛

在本节中我们证明雅可比和高斯塞德尔方法对于严格对角占优矩阵收敛.这是定理2.10和定理2.11的内容.
雅可比方法可以写做xk+1=-D-1(L+U)xk+D-1b(2.43)附录A中的定理A.7决定这类迭代的收敛.根据这个理论，我们需要知道谱半径ρ(D-1(L+U))<1以保证雅可比方法收敛.如下所示，这正是严格对角占优矩阵所包含的意思.
定理2.10的证明　令R=L+U表示矩阵的非对角线部分.为了检验ρ(D-1R)<1，令λ是矩阵D-1R的特征值，对应的特征向量是v.选择v使得‖v‖∞=1，因而对于一些1≤m≤n，元素vm=1，而其他的元素不大于1.（任何特征向量除去最大的元素都可以满足上面的向量要求.任何常数乘上特征向量后依然是特征向量，并具有相同的特征值.）特征值的定义意味着D-1Rv=λv或者Rv=λDv111由于rmm=0，把这个向量方程的第m个元素的绝对值去掉意味着rm1v1+rm2v2+…+rm,m-1vm-1+rm,m+1vm+1+…+rmnvn=λdmmvm=λ‖dmm由于所有vi≤1，左侧最多是∑j≠mrmj,根据严格对角占优假设，其小于dmm.这意味着λdmm把高斯塞德尔方法写成（2.43）的形式得到xk+1=-(L+D)-1Uxk+(L+D)-1b很显然，高斯塞德尔方法是否收敛取决于矩阵(L+D)-1U(2.44)的谱半径是否小于1.下一个定理表明对于严格对角占优矩阵，意味着特征值满足谱半径的要求.
定理2.11　如果n×n矩阵A是严格对角占优矩阵，则（1）A是非奇异矩阵，（2）对每一个向量b和初始估计，对Ax=b使用高斯塞德尔方法收敛到解.
证明　令λ是（2.44）的特征值，对应的特征向量是v.如同前面的证明，选择特征向量使得vm=1且其他元素在数量上更小.注意到L的元素是aij，i>j，U的元素是aij（imami<λamm-∑i≤λamm+∑imamivi≤∑i>mami可以得到λ<1，完成证明.112

2.5.4　稀疏矩阵计算

基于高斯消去的直接方法，通过有限步数的计算就可以得到解.为什么我们要研究并使用迭代方法，特别是迭代方法只能得到近似解，并需要多步计算以得到收敛.
有两个主要的原因让我们使用诸如高斯塞德尔这样的迭代方法.两个原因都是基于如下的事实：迭代方法中的一步计算，仅仅需要LU分解的浮点计算所需要时间的一小部分.正如我们在本章前面所确立的，对于n×n的矩阵做一次高斯消去需要n3次的操作.雅可比方法的一步，仅需要大约n2次的乘法（对于每个矩阵元素做一次乘法操作）以及大约相同数量的加法.问题是需要多少步才能在用户提供的容差里得到收敛.
如果已知解的较好的近似，在这种特定的情况下支持使用迭代技术.例如，假设知道Ax=b的解，随后A与b同时或者单个仅仅发生小的变化.我们可以想象一个动态系统，当A和b改变时，对A和b进行持续的度量，因而需要持续更新精确的解向量x.如果把前面问题的解作为新的相似问题的初始估计，使用雅可比或者高斯塞德尔可以得到更快的收敛.
假设问题（2.42）中的b和最初的b=［2.5，1.5，1，1，1.5，2.5］相比发生微小变化，得到新的b=［2.2，1.6，0.9，1.3，1.4，2.45］.我们检查后发现真实解从［1，1，1，1，1，1］变为［0.9，1，1，1.1，1，1］.假设从前面的表中，在内存中已经有了高斯塞德尔迭代的第6步x6，用于初始估计.继续进行高斯塞德尔迭代，其中使用新的b以及有帮助的初始估计x6，仅仅多做一步就可以得到好的结果.后面的两步如下:
x7x8x7x8x7x80.89800.89940.96590.99270.99711.00050.99800.98891.08921.09660.99931.0003
这种技术通常称为修饰，这是由于该问题从一个近似解开始，该近似解可能对应前面一个相关问题的解，然后仅仅修饰近似解使其更加精确.修饰在实时应用中很常见，随着时间的流逝，数据会进行更新，因而相似的问题需要不断重复求解.如果系统很大，时间很短，可能难以在给定的时间内进行完整的高斯消去或者仅仅进行回代.如果解的变化不大，相对代价低的迭代方法的少数几步，当解随着时间流逝可能保持足够的精度.
第二个使用迭代方法的主要原因是用于求解稀疏的方程组.当矩阵中的很多元素都是0系数矩阵被称为稀疏矩阵.通常，对于稀疏矩阵中的n2个矩阵元素，只有O（n）个非零元素.一个完全的矩阵恰恰相反，其中没有几个元素可以假设为0.对稀疏矩阵使用高斯消去经常会导致填充，这使得其中系数矩阵由于必要的行交换从稀疏变为不稀疏.由于这个原因，高斯消去以及PA=LU实现的效率对于稀疏矩阵变得可疑，在这种情况下迭代方法是一种合理的选择.113
例2.24可以扩展为如下的稀疏矩阵。
例2.25　雅可比方法求解例2.24中包含100000个方程的版本.
令n是一个偶数，考虑n×n矩阵A对角线元素为3，-1出现在上对角线和下对角线，1/2出现在（i，n+1-i）位置，其中i=1，…，n，除了i=n/2以及n/2+1.对于n=12，A=3-100000000012
-13-10000000120
0-13-1000001200
00-13-100012000
000-13-10120000
0000-13-100000
00000-13-10000
0000120-13-1000
00012000-13-100
001200000-13-10
0120000000-13-1
12000000000-13(2.45)定义向量b=（2.5，1.5，…，1.5，1.0，1.0，1.5，…，1.5，2.5），其中有n-4次重复的1.5和2次重复的1.0.注意到，如果n=6，A和b定义了例2.24的系统.对于一般的n，系统的解为［1，…，1］.A中没有任何一行的非零元素的个数大于4.由于在n2的元素中只有小于4n的元素非零，我们可以称该矩阵为稀疏矩阵.
如果我们要求解这样的方程组，其中n=100000或者更大，可以选择什么方式来做？把稀疏矩阵A看做一个完全矩阵意味着保存n2=1010个元素，其中每个都是浮点双精度数，并需要8字节来保存.注意到8×1010字节大约是80G字节.根据你的计算机配置，可能难以把所有n2个元素放进RAM.
不仅空间大小是个问题，时间也是一个问题.高斯消去所需要的操作次数的数量级是n3≈1015.如果你的计算机运行速度是几个GHz（109次时钟每秒），每秒进行浮点操作的次数大约是108.因而1015/108＝107，是对高斯消去所需要的秒数的一个合理估计.一年有3×107秒.尽管这只是一个很粗略的估计，但是很显然对于这个问题不是一天就能干完的.
另一方面，一步迭代大约需要2×4n=800000次操作，对于每个非零矩阵元素进行两次.我们可以做大约100步雅可比迭代，所需的计算次数远小于108，这在一个现代PC上，需要1秒或者更短的时间就可以完成.对于刚才定义的系统，其中n=100000，后面的雅可比代码jacobi.m仅需要50步就可以从初始估计（0，…，0）收敛到解（1，…，1），精确到小数点后6位.50步迭代在一般计算机上需要的时间不到1秒.

114
注意到前面的代码中一些有趣的方面.程序sparsesetup.m使用MATLAB的spdiags命令，把矩阵A定义为一个稀疏矩阵结构.本质上，这意味着把矩阵表示为一组三元数（i，j，d），其中d是矩阵（i，j）位置上的实数.而不需要在内存中保存n2个元素，仅仅保存需要的元素.spdiags命令拿出矩阵的一列，并把它们放在主对角线上，或者是主对角线的上侧对角线与下侧对角线.
MATLAB的矩阵乘法命令可以无缝地对稀疏矩阵进行计算.例如，另一种处理代码的方式是使用MATLAB的命令lu直接求解系统.但是，对于那个例子，即使A是稀疏矩阵，从高斯消去中得到的上三角矩阵U会面临填充的问题.例如，从n=12规模的矩阵A使用高斯消去法得到的上三角矩阵U为3-1.00000000000.500
02.7-1.000000000.5000.165
002.6-1.0000000.5000.1870.062
0002.6-1.0000000.5000.1910.0710.024
00002.618-1.00000.5000.1910.0730.0270.009
000002.618-1.0000.1910.0730.0280.0100.004
0000002.618-0.9270.0280.0110.0040.001
00000002.562-1.032-0.012-0.005-0.001
000000002.473-1.047-0.018-0.006
0000000002.445-1.049-0.016
00000000002.440-1.044
000000000002.458由于U变成完全矩阵，前面提到的内存约束又出现了，并成为一个求解的限制.在求解过程中需要n2的内存来保存U.使用迭代方法，在执行时间和存储方面，效率会得到数量级的提高.
2.5节习题
1.计算雅可比和高斯塞德尔方法的前两步，初始向量是［0，…，0］.
(a) 3-1
-12u
v=5
4　　(b)2-10
-12-1
0-12u
v
w=0
2
0　　(c)311
131
113u
v
w=6
3
5

115
2.对方程组重新组织，得到一个严格对角占优矩阵.使用雅可比和高斯塞德尔方法做两步，初始向量是［0，…，0］.
(a) 　u+3v=-1
5u+4v=6(b) u-8v-2w=1
　u+v+5w=4
　3u-v+w=-2(c) 　u+4v=5
　v+2w=2
4u+3w=0
3.在习题1的系统上使用SOR方法做两步.初始向量［0，…，0］，ω=1.5.
4.在习题2的重新排列的系统上使用SOR方法做两步.初始向量［0，…，0］，ω=1和1.2.
5.令λ是一个n×n的矩阵A的特征值.（a）证明Gershgorin Circle定理:存在一个对角线元素Amm，满足Amm-λ≤∑j≠mAmj（提示:如定理2.10的证明一样，从满足‖v‖∞=1的特征向量v开始.）（b） 证明严格对角占优矩阵没有0特征值.这是定理2.10（1）的另一个证明.
2.5节编程问题
1.使用雅可比方法求解稀疏系统，精确到小数点后6位（用无穷范数表示的前向误差），其中n=100和n=100000.正确解是［1，…，1］.报告所需要的步数和后向误差.系统如下：3-1
-13-1

-13-1
-13x1

xn=2
1

1
22.使用雅可比方法求解稀疏系统，精确到小数点后3位（用无穷范数表示的前向误差），其中n=100.正确解是［1，-1，1，-1，…，1，-1］.报告所需要的步数和后向误差.系统如下：21
121

121
12x1

xn=1
0

0
-13.重写程序2.2进行高斯塞德尔迭代.求解例2.24中的问题验证你的工作.
4.重写程序2.2进行SOR.使用ω=1.1，再次验证例2.24.
5.执行编程问题1中的步骤，n=100.（a） 高斯塞德尔方法，（b） SOR，ω=1.2.
6.执行编程问题2中的步骤，（a） 高斯塞德尔方法，（b） SOR，ω=1.5.116
7.使用编程问题3的程序，确定形如（2.38）的系统，使用高斯塞德尔方法在1秒内所能够精确求解的系统最大规模.报告对于不同的n所需的时间，以及前向误差.