《数值分析（原书第2版）》—— 2.1 高斯消去法

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

2.1　高斯消去法

考虑方程组x+y=3
3x-4y=2(2.1)
71具有两个方程两个变量的系统可以从代数方面考虑，也可以从几何方面考虑.从几何观点，每一个线性方程表示xy平面中的一条直线，如图2.1所示.在点x=2，y=1处两条直线相交，交点满足两个方程，这正是我们要寻找的解.

几何观点有助于我们可视化求解系统，但是为了计算具有很多精确数位的解，我们回到代数方法.这种称为高斯消去法的方法，可以有效地求解具有n个未知数的n个方程.在下面的几节中，我们将研究用于典型问题最优的高斯消去法.

2.1.1　朴素的高斯消去法

我们首先描述最简单的高斯消去法.事实上，该方法过于简单，以至于有时难以完成求解，更不要说找到精确解.“朴素”方法的改进将会在下一节中进行介绍.
对于线性方程组系统可以施加三个有用的操作，使用该操作可以生成等价的系统.这意味着二者具有等价解.这些操作如下:
（1） 两个方程彼此交换位置.
（2） 在一个方程上加上或者减去另外一个方程的倍数.
（3） 对于一个方程乘上一个非零的常数.
对于方程（2.1），我们可以从第二个方程中减去第一个方程的3倍，从而从第二个方程中消去x变量.从第二个方程中减去3·［x+y=3］得到方程组x+y=3
-7y=-7(2.2)从底部的方程开始，我们可以“从后向前”求解以得到所有的解，过程如下-7y=-7→y=1以及x+y=3→x+(1)=3→x=2因而（2.1）的解是（x，y）=（2，1）.72
相同的消去过程可以在没有变量参与的情况下进行，把系统写成如下的表格形式:11|3
3-4|2→从第2行中减去第1行的3倍→11|3
0-7|-7(2.3)表格形式的优势在于在消去过程中隐藏了变量.当表格左边的方块数组变成“三角”时，我们可以从底部开始，由后向前求解对应方程的解.
例2.1　对表格形式应用高斯消去法求解三个未知数、三个方程的系统：x+2y-z=3
2x+y-2z=3
-3x+y+z=-6(2.4)表格形式如下：12-1|3
21-2|3
-311|-6(2.5)需要两步消去第1列:12-1|3
21-2|3
-311|-6→从第2行中减第1行的2倍→12-1|3
0-30|-3
-311|-6
→从第3行中减去第1行的-3倍→12-1|3
0-30|-3
07-2|3还需要一步消去第2列：12-1|3
0-30|-3
07-2|3→第3行中减去第2行的-73倍→12-1|3
0-30|-3
00-2|-4返回方程组为x+2y-z=3
-3y=-3
-2z=-4(2.6)我们可以按z，y，x的顺序求解变量x=3-2y+z
-3y=-3
-2z=-4(2.7)后面的一个过程称为回代或者后向求解.这是由于完成消去后，方程组可以很容易地从最下面的一个方程开始求解.对应的解为x=3,y=1,z=2.73

2.1.2　操作次数

在本节中，我们对高斯消去法中的两部分：消去步骤和回代步骤的计算中涉及的计算次数进行估计.为了进行估计，把前面两个例子中讲述的高斯消去法写做一般的情况会有所帮助.首先，回忆两个整数求和对应的事实.
引理2.1　对于任何正整数n，（a）1+2+3+4+…+n=n（n+1）/2，（b）12+22+32+42+…+n2=n(n+1)(2n+1)/6.
关于n个未知数的n个方程的表格形式如下：a11a12…a1n|b1
a21a22…a2n|b2
|
an1an2…ann|bn为了进行消去，我们需要使用允许的行操作把矩阵的下三角部分都转换为0.
我们可以使用循环进行消去

其中，“消去列j”的意思是“使用行变换把0放在主对角线下的每个位置，包括aj+1,j,aj+2,j,…,anj”.例如，为了消去第1列，需要把0放在a21，…，an1的位置.这可以写出如下的循环，该循环在前面的循环内部:

在这个双重循环的内部，应用行操作把aij元素设置为0.例如，第一个要消去的元素是a21.为了做到这一点，假设a11≠0，第1行乘上a21/a11和第2行相减，前面的两行从a11a12…a1n|b1
a21a22…a2n|b2变为a11a12…a1n|b1
0a22-a21a11a12…a2n-a21a11a1n|b2-a21a11b1在这些计算中，需要一次除法（找到乘子a21/a11），以及n步乘法和n步加法.用于消去第1列中ai1的行操作，即a11a12…a1n|b1
|
0ai2-ai1a11a12…ain-ai1a11a1n|bi-ai1a11b1需要相似的操作.74
上面描述的过程中，只要a11不是0就可以进行.a11和其他的aii在高斯消去中作为除数，这些数字称为主元.正如我们前面所解释的，0主元会使算法终止.现在我们忽略这个问题，然后在2.4节中再认真考虑这个问题.
回到计算次数的统计，注意到每次在消去第1列中ai1时使用1次除法、n次乘法和n次加/减法，放在一起一共有2n+1次运算.把0放在第1列需要重复这2n+1次运算n-1次.
当第1列被消去后，主元a22以相同的方式消去第2列以及余下的列.例如，用于消去aij的行操作如下：00ajjaj,j+1…ajn|bj
|
000ai,j+1-aijajjaj,j+1…ain-aijajjajn|bi-aijajjbj在我们的标记中，a22指的是经过第1列消去后在该位置上得到的数值，这和原始的a22并不相同.消去aij的操作需要1次除法、n-j+1次乘法以及n-j+1次加/减法.
把这一步插入相同的双重循环中得到

在这个代码段中有两点需要讨论.第一，令k从j遍历到n会在aij位置生成0;但是从j+1遍历到n是最有效的代码.后者不会在aij位置上放上0，这些非0位置这恰好是我们要消去数字的位置!尽管这看起来是一个错误，但是注意到在余下的高斯消去过程与回代过程中我们不会用到这些数值，因而实际上，从效率的观点上来看，在这些位置上放上0是一种浪费.第二，我们使用MATLAB的error命令在遇到0主元的时候终止程序.正如前面提到的，对于这种可能性，应该更加认真地对待，在2.4节讨论行置换时会再讨论这个问题.
我们可以数出来在高斯消去中所有计算的次数.消去每个aij需要如下所示的操作次数，包含除法、乘法和加/减法:0
2n+10
2n+12(n-1)+10
2n+12(n-1)+12(n-2)+10

2n+12(n-1)+12(n-2)+1…2(3)+10
2n+12(n-1)+12(n-2)+1…2(3)+12(2)+10

75以相反的顺序把所有的操作次数加起来更简单.从右边开始，操作次数是n-1j=1ji=12(j+1)+1=n-1j=12j(j+1)+j
=2n-1j=1j2+3n-1j=1j=2(n-1)n(2n-1)6+3(n-1)n2
=(n-1)n2n-13+32=n(n-1)(4n+7)6
=23n3+12n2-76n其中使用了引理2.1.
1.高斯消去法中消去步骤的操作次数.
n个方程n个未知数的消去计算，可以在23n3+12n2-76n次操作后完成.
一般来说，由于在不同的计算机处理器上实施的细节不同，计算次数对应的数量级比精确的操作次数更重要.重要的是，执行运算的次数和算法所需的时间成比例.一般地，我们在消去步骤中将进行23n3次操作.当n很大时，这个估计相对准确.
在所有的消去完成后，表格变为上三角形式:a11a12…a1n|b1
0a22…a2n|b2
|
00…ann|bn方程形式如下a11x1+a12x2+…+a1nxn=b1
a22x2+…+a2nxn=b2

annxn=bn(2.8)其中，aij指的是修正后的，而不是原始的矩阵元素.为了完成求解x的运算，我们必须进行回代过程，这可以通过简单重写（2.8）实现:x1=b1-a12x2-…-a1nxna11
x2=b2-a23x3-…-a2nxna22

xn=bnann(2.9)

76复杂度　运算次数表明直接使用高斯消去法求解关于n个未知数的n个方程所需要的操作次数是O（n3），这对于求解大规模系统是一个非常有用的事实.例如为了计算n=500的方程组在一个特定的计算机上所花费的时间，我们可以求解一个n=50的方程组然后将时间乘上103=1000.
由于方程组是非0系数的三角形状，我们首先从底部开始，然后逐渐向上求解对应的方程.在这种方式下，再计算下一个位置变量时所需要的xi的值已经知道了.操作计数如下：1+3+5+…+(2n-1)=ni=12i-1=2ni=1i-ni=11=2n(n+1)2-n=n2MATLAB句法中，回代步骤如下：

2.高斯消去法中回代过程的操作次数.
n个方程ｎ个未知数的三角形系统的回代过程可以使用n2次操作完成.
两种操作计数放在一起，表明高斯消去由两个不等同的部分组成：相对计算代价庞大的消去过程，和相对计算代价小的回代过程.如果忽略具有更低阶乘除法计算次数的部分，我们发现消去的操作次数是2n3/3，回代的操作次数则是n2.
我们将使用缩写的术语“大写O”即相似数量级表示“算法的复杂度，”消去是一个O（n3）的算法，而回代的复杂度是O（n2）.
这种用法表明在n很大时，其中n的低阶在比较中可以忽略.例如，如果n=100，只有大约1%的高斯消去的计算花在回代步骤.总体来说，高斯消去用去2n3/3+n2≈2n3/3次的计算.换句话说，对于n，在复杂度计算中的低阶项对于算法运行时间的估计没有大的影响，并可以忽略.
例2.2　估计在一个规模为500个方程500个未知变量的系统中，进行回代的时间代价，已知该计算机系统中消去过程花了1秒时间.
由于我们知道消去计算比回代计算更加花费时间，所以结果会是1秒的一小部分.使用近似数字2(500)3/3作为消去步骤中乘法/除法计算的次数，(500)2是回代过程的计算次数，我们估计回代所花费的时间是(500)22(500)3/3=32(500)=0.003秒该例子表明两点:1）n的更低阶的幂在算法次数统计中可以安全地忽略.2）高斯消去77法的两部分在运行时间上可能差得很远，全部的运行时间为1.003秒，其中几乎所有的时间都被消去步骤所消耗.下个例子则显示了第三个要点：尽管回代时间有时可以忽略，但它也可能是计算的一个重要的组成部分.
例2.3　在特定的计算机上，一个5000×5000三角矩阵的回代时间是0.1秒.估计一下使用高斯消去法求解3000个方程3000个未知变量的一般系统所花费的时间.
计算机可以在0.1秒中进行(5000)2次操作，或者做(5000)2(10)=2.5×108次操作/秒.求解一个一般的（非三角）系统需要大约2（3000）3/3次操作，这些操作可以在大约2(3000)3/3(5000)2(10)≈72秒完成.
2.1节习题
1.使用高斯消去法求解方程组:
(a) 2x-3y=25x-6y=8　　　　　　　　(b) x+2y=-12x+3y=1　　　　　　　　　　(c) -x+y=23x+4y=15
2.使用高斯消去法求解方程组:
(a) 2x-2y-z=-2
4x+y-2z=1　
-2x+y-z=-3(b) x+2y-z=2
3y+z=4
2x-y+z=2(c) 2x+y-4z=-7
x-y+z=-2
-x+3y-2z=6　
3.使用回代求解:
(a) 3x-4y+5z=2　
3y-4z=-1
5z=5　(b) x-2y+z=2
4y-3z=1
-3z=3
4.求解表格形式
(a) 3-4-2|3
6-61|2
-382|-1　(b)21-1|2
62-2|8
46-3|5
5.使用高斯消去法的近似操作次数为2n3/3，如果n乘上3倍，需要花多长时间求解n个方程n个未知变量的系统？
6.假设你的计算机做5000次方程回代需要0.005秒.如果回代的近似操作次数为n2以及消去的近似操作次数为2n3/3，估计使用多长时间才能在该规模的问题上做一次完整的高斯消去.把你的结果舍入到最接近的秒数.
7.假设一个给定的计算机需要0.002秒计算4000×4000的上三角矩阵方程的回代，估计求解9000个方程9000个未知变量的一般方程组需要多长时间.把你的结果舍入到最接近的秒数.
8.如果一个3000个方程3000个未知变量的系统在一个给定的计算机上，使用高斯消去的时间是5秒，那么每秒可以做多少次相同规模的回代？78
2.1节编程问题
1.把本节中所有的MATLAB程序代码段放在一起生成一个“朴素”的高斯消去法（不允许进行行交换）的完整M程序.使用该程序求解习题2中的问题.
2.令H表示n×n的希尔伯特矩阵，其中（i，j）元素是1/（i+j-1）.使用编程问题1中的MATLAB程序求解Hx=b，其中b是一个元素全为1的向量，（a） n=2，（b） n=5，（c） n=10.