Jacobi与Gauss-Seidel迭代的比较及r算法的MATLAB实现

Jacobi与Gauss-Seidel迭代的比较及r算法的MATLAB
实现
胡志成
【摘要】探讨比较了Jacobi迭代和Gauss-Seidel迭代,尤其是它们的可并行性.给出了它们在MATLAB中的一个实现.通过算例说明,算法的并行化实现可以大幅度地提升Jacobi迭代的效率.
【期刊名称】《高师理科学刊》
【年(卷),期】2018(038)003
【总页数】4页(P59-61,84)
【关键词】Jacobi迭代;Gauss-Seidel迭代;并行性
【作者】胡志成
【作者单位】南京航空航天大学理学院,江苏南京 210016
【正文语种】中文
【中图分类】O241.6;G642.0
Jacobi迭代和Gauss-Seidel迭代是解线性方程组的2个经典的迭代法，简单又易于实现，是大学计算方法相关课程必讲的知识点．然而教材[1-3]往往只注重理论而轻视算法的程序实现．就这2个迭代法而言，表现为着重讨论迭代法的收敛性和收敛速度，而甚少给出或者仅给出一个粗糙低效的算法实现．至于Jacobi迭代的并行性，探讨的更是少之又少．由于课时的限制，在工科计算方法教学中更是如
此．这一方面是由于一般程序语言关于算法的并行化实现并不那么容易，另一方面则是因为在课堂上短时间内举例说明算法并行化实现带来的效率提升较为困难．然而在课堂上适当地给出算法的实现和具体算例，显然有助于学生对算法本身的理解．本文首先回顾了Jacobi迭代和Gauss-Seidel迭代的基本格式和传统比较，然后探讨了它们的可并行性，并给出了在MATLAB中的一个实现．文中给出的算例不仅可以在课堂的短时间内展示结果，也能说明算法的并行化实现可以大幅度地提升Jacobi迭代的效率．
设是阶非奇异方阵，为维列向量．若的主对角元，，则线性方程组可以用Jacobi 迭代和Gauss-Seidel迭代求解．其迭代格式分别为
若再将表示成，其中：；；，则Jacobi迭代可用矩阵表示为
相应地，Gauss-Seidel迭代的矩阵形式为
关于Jacobi迭代和Gauss-Seidel迭代的分析和比较，在众多的文献中都有提及[1，4-6]，主要归结为：
（1）收敛性．当且仅当各自的迭代矩阵（分别是和）的谱半径小于1时，该迭代法收敛．若是严格对角占优阵，则2个迭代法都收敛．但一般地， Jacobi迭代收敛不能保证同一问题的Gauss-Seidel迭代也收敛，反之亦然．
（2）收敛速度．当Gauss-Seidel迭代和Jacobi迭代两者都收敛时，多数情况下Gauss-Seidel迭代比Jacobi迭代收敛得更快．对于一些模型问题，Gauss-Seidel迭代的收敛速度甚至是Jacobi迭代的2倍以上[5]58．也就是说，从同一初值出发，为得到满足精度要求的解，Gauss-Seidel迭代所需的步数可能比Jacobi迭代所需步数的一半还要少．
（3）每步计算量．由于矩阵的求逆并不那么容易，具体的算法实现多采用分量形式，即式（1）和式（2）．比较式（1）和式（2）可知，2个迭代法每步的计算量相等．
（4）存储量．Gauss-Seidel迭代只需1个向量来存储或的分量，而Jacobi迭代需要2个向量来分别存储和．
比较可知，在传统的算法实现中，Gauss-Seidel迭代比Jacobi迭代往往更具优势．随着多核处理器的迅猛发展，并行计算渐次成为提升计算效率的重要技术手段，因而一个算法的可并行性也成为了衡量其优劣的重要指标[7-9]．
进一步比较迭代法（1）和（2），易知Jacobi迭代中，各分量的计算是相互独立的，因而具有高度的并行性；而Gauss-Seidel迭代中，各分量是依次顺序计算的，因而本质上是一个串行算法．这就意味着，利用算法的可并行性，如果实现得当，即使收敛速度较慢，Jacobi迭代仍然有可能在更短的时间内得到满足精度要求的
数值解．
为了加深学生对此的理解，并让学生体会到算法实现，特别是并行化实现的重要性，有必要给出适当的实例和程序实现．
MATLAB是一种面向科学与工程计算的数学软件，因其强大的数值计算和图形绘
制能力、广泛的工具包支持以及使用的便利性，深受教学和科研工作者的喜爱，是工科大学生需要掌握的重要软件之一[10]．借助其灵活的矩阵向量运算，给出Jacobi迭代（1）在MATLAB环境中的一个实现（见算法1）和Gauss-Seidel迭代（2）的一个MATLAB实现（见算法2）．可以看到，算法1完全避免了Jacobi迭代各分量之间的for循环，而算法2则无法避免Gauss-Seidel迭代各分量之间的for循环，这是导致2个算法执行效率巨大差异的直接原因．
算法1（Jacobi 迭代）：
function [x，k，res] = Jacobi（A，b，x0，tol，kmax ）
k = 0； res = 1；
while（k <= kmax）& res > tol）
x =（b - A * x0 + diag（A）．* x0 ）．/ diag（A）；
k = k+1；
res = max（abs（b - A * x ））；
x0 = x；
end
算法2（ Gauss-Seidel 迭代）：
function [x，k，res] = GaussSeidel（A，b，x0，tol，kmax ）
k = 0；res = 1；
x = x0；
while （k <= kmax）&（res > tol）
for i = 1：length（x）
x（i）=（b（i）- A（i，1：end）* x + A（i，i）* x（i））/ A（i，i）；
end
k = k+1；
res = max（abs（b - A * x））；
end
以零向量为初值，以残量的无穷范数为收敛准则，考察比较这2个算法求解线性方程组的数值结果，包括迭代步数、每步耗时和总耗时．为了使测试算例具有一定的普遍意义，且能在课堂上较短的时间内展示比较结果．用算法3来随机生成阶线性方程组的系数矩阵．
算法3（随机生成一个严格对角占优矩阵）：
function A = randSDDMatrix（n）
nzmax = randi（[n*5，round（n*10）]，1）；
i = [randi（[1，n]，nzmax，1）]；
j = [randi（[1，n]，nzmax，1）]；
s = rand（nzmax，1）；
A = sparse（i，j，s，n，n）；
x = randi（[1，1e1]，n，1）+ sum（ abs（A），2）；
ii = 1：n；
A = sparse（ii，ii，x，n，n）+ A；
该矩阵是稀疏的，并且是严格对角占优的，可以确保阶数较大时，2个迭代法仍能快速收敛．当然，右端项也是随机生成的．对每一个阶数，都测试100组和，并取其平均数据（见表1）．
由表１可以看出，虽然Jacobi迭代（1）和Gauss-Seidel迭代（2）的每步计算量相等，但是算法1的每步计算效率显著优于算法2．这最终导致：虽然Jacobi 迭代的总迭代步数是Gauss-Seidel迭代的２倍以上，但是总耗时却远少于Gauss-Seidel迭代．而且随着矩阵阶数增加，在本文的算法实现下，这个优势将进一步扩大．
本文比较分析了Jacobi迭代和Gauss-Seidel迭代，强调说明Jacobi迭代是一个高度并行的算法，而Gauss-Seidel迭代则是一个串行算法．因而如果算法实现得当，即使收敛速度更慢的时候，Jacobi迭代仍然可能比Gauss-Seidel迭代更快地给出数值解．在教学实践中，课堂上展示这些MATLAB实现和数值算例，明显地加深了学生对这2个迭代法的理解，也使学生意识到了算法实现的重要性，获得了较好的教学效果．
【相关文献】
[1] 倪勤，王正盛，刘皞．数值计算方法[M]．北京：高等教育出版社，2012
[2] 李庆杨，王能超，易大义．数值分析[M]．5版．北京：清华大学出版社，2008
[3] 李有法，李晓勤．数值计算方法[M]．2版．北京：高等教育出版社，2005
[4] 杜衡吉，徐昆良．Jacobi和Gauss-Seidel迭代法求解线性方程组的分析及应用[J]．曲靖师范学院学报，2011，30（3）：46-50
[5] 白红梅．Jacobi迭代法与Gauss-Seidel迭代法的收敛性比较分析[J]．呼伦贝尔学院学报，2009，17（6）：55-58
[6] 包丽君．“数值分析”中实践教学的探讨[J]．宁波广播电视大学学报，2008，6（2）：91-93
[7] 徐树方，高立，张平文．数值线性代数[M]．北京：北京大学出版社，2000
[8] 杨学军．并行计算六十年[J]．计算机工程与科学，2012，34（8）：1-10
[9] 韩建勋．并行计算在矩阵运算中的应用[J]．信息与电脑：理论版，2017（5）：93-96
[10] 韦慧，蒋利华．基于MATLAB的大学工科《计算方法》课程教学方法探讨[J]．内江科技，2017（4）：72-73。