其他
怎样迭代求解线性方程组？
加星标，才能不错过每日推送！方法见文末动图
在科学计算和工程应用领域，经常涉及求解线性方程组问题。迭代法是解决这类问题的有力武器，通过多次迭代逐步逼近精确解，在处理大型稀疏线性方程组时具有重要优势。《返朴》之前已刊发数篇函数迭代的相关文章，本文从基础概念讲起，由浅入深介绍了迭代法求解线性方程组的准备知识和具体过程。
记得有个大科学家曾经说过，这个世界是非线性的。然而非线性方程一般不能直接求解，即解析解虽可证存在却无具体表达式，因而迭代法几乎是唯一可行的办法。比如说，从介值定理可知，方程x = cos x在区间(0, 1)内定有一解，但没有一步到位的法子找到它，人们只能用基于介值定理的二分法或基于切线逼近的牛顿法，来求得此方程的迭代近似解。这样，从最古老的巴比伦平方根迭代法，到今日非线性方程组数值解的最重要方法——牛顿迭代法，人们一直热衷于迭代法的理论探索和创新实践。
既然有非线性方程，就必然有线性方程，而且后者在科学技术和日常生活中同样随处可见。电话公司的客户安排、航空公司的航班分配、运输公司的运货路线……都可以用线性规划来优化成本。四十年前，当我在南京大学数学系时，读到一篇关于线性规划的综述性文章，其中写道：全世界的计算机大约有百分之六十的时间，都在求解线性规划问题！
另一方面，微分概念告诉我们，非线性函数可以局部地用线性函数来近似，构造牛顿迭代法本质上就基于这个简单的观察。只要非线性函数f在一个点a存在导数，那么在a点的附近，函数值f(x)就约等于线性函数值f(a) + f'(a)(x-a)，其几何直观则是函数图像的切线在切点附近与曲线相差无几。
此外，在物理科学中，许多描述自然规律的常微分方程或偏微分方程对于未知函数以及它们的导数或偏导数都呈现出线性关系，如热传导方程或波动方程。当用计算数学中的差商逼近导数技术或有限元思想来离散化这些连续方程时，获得的差分方程也是线性的。它们可以写成一般形式Ax = b，其中A是一个n行、n列的系数矩阵，b是一个n维的已知列向量，x是一个n维的未知列向量。这个线性方程看上去像一元一次方程ax = b一样简单，但如果按照矩阵乘法的法则将方程左边每个分量的代数表达式全部写出来，结果就是一组含有n个未知数x1, x2, …, xn的n个n元一次方程。如果将方矩阵A中第i行、第j列的元素记为aij，将列向量b的第i个分量记为bi，那么线性方程组Ax = b展开后的第i个方程为
ai1 x1 + ai2 x2 + … + ain xn = bi, i = 1, 2, …, n。
计算数学中的一个子学科“数值线性代数”有个基本的职责：怎样有效求解上面的线性方程组？当n = 2或3，最多不超过4时，我们早已在中学课堂上学会了怎样求解二元一次方程组、三元一次方程组，等等。采用的方法无非是两种，一种叫做代入法，另一种称为加减法。它们都是基于同一个思想：通过消去方程中的一个或几个未知数，设法获得只含有一个未知数的方程，比方说-2y = 3，从而解出这个未知数。
中学生所学的上述方法，实际上是大名鼎鼎的高斯消去法的简单实现，发明者正是德国天才高斯 (Carl Friedrich Gauss，1777-1855) 。高斯消去法的基本思想还是如上所说：将含有多个未知数的方程化简到只剩下一个未知数，因为只含一个变元的线性方程如-2y = 3总是可以一步到位解出y = -3/2的。基于这个非常简单的思路，高斯想出了一个点子，将线性方程组的系数方阵程式化系统性地化约成主对角线下方的元素全为零的上三角矩阵，这样新的等价的线性方程组的最后一个方程只含一个未知数，倒数第二个方程含有两个未知数，如此等等。因而，通过所谓的“回代法”，从最后一个方程解出唯一的未知数，然后将它的值代入到倒数第二个方程，结果该方程也只含一个未知数，马上可以解出，再将这两个未知数的值代入到倒数第三个方程，同法解出下一个的未知数。如此这般，就可以依次解出所有的未知数。
高斯消去法在经过有限次代数运算后就可以完全找出给定线性方程组的精确解（如果每步运算的误差为零的话）。所以它属于求解线性方程组的第一类数值方法——直接法，这是解线性方程组优越于解非线性方程组的一个地方。然而，对有巨量（如成千上万个或更大数量级）变元个数或者特别多的系数为零（比如系数矩阵为三对角或其他稀疏类型）的那种方程组，直接法常常没有第二类方法——迭代法有效。对线性偏微分方程直接用差商代替偏导数的差分方法，或将微分方程问题变为数学上等价的变分问题而进行数值优化的有限元方法，所得到的大型线性方程组一般都带有特定结构的稀疏系数矩阵，这时，迭代法几乎是最有效的算法了。
我们下面只讲解线性方程组的迭代法。让我们回忆求解单变量非线性方程的迭代法，一般形式是xn = f(xn-1), n = 1, 2, …，其中f是将定义域区间映到自身的一个函数，x0是迭代所取的初始点。然而，对于线性迭代法，迭代函数不再是一个自变量的线性函数，而是有n个自变量的线性向量函数。由于字母n现在另有他用，我们将用字母k代表迭代次数的下标，而将多变量线性向量函数用y = L(x)表示，其中L(x)的表达式是Mx + c，M是一个有n行和n列的给定矩阵（也称为n阶正方矩阵或n阶方阵），c是一个给定的n维列向量，其分量是c1, c2, …, cn，x = (x1, x2, …, xn)T是n维变元列向量，其中的上标字母T表示转置运算，即一般的有m行和n列的矩阵A的转置矩阵AT有n行和m列，其第j行的元素是原矩阵第j列的对应元素。
本章将要面对的线性迭代程序可以写成
xk = L(xk-1) = Mxk-1 + c, k = 1, 2, …; x0是给定的。
我在以前的科普文章里写过，单变量非线性方程的迭代法收敛的一个充分条件，是迭代函数在不动点的导数绝对值小于1。那么，如果被迭代的是有多个变量的线性向量函数，什么样的条件可以充分到保证迭代法收敛呢？所谓序列的收敛，本质上是用距离这个概念来定义的，这个概念我们在以前提到巴拿赫不动点定理或压缩映射定理时已经见到过。我们再来回忆有关内容。在一个抽象的距离空间(X, d)里，一个由X中的元素所组成的无穷序列{xk}被说成是收敛到X中的一个元素x，指的是由距离函数d所确定的非负数列{d(xk, x)}当k趋向于无穷大时趋向于0。对于单变量函数迭代法的收敛性问题，我们实际上选取了一个简单而具体的距离空间，它就是所有实数全体构成的一维欧几里得空间R，其任意两个实数之间的距离就是实数轴上两点x和y之间的通常距离|x – y|。
上面两数之间的距离可用线段长度的说法来等价地定义。任意一个实数a的绝对值|a|在几何上的意义是以数轴上代表a的那个点和原点0为两端点的线段的长度。不难看出，两个数x和y之间的距离就是数x – y所对应的那条线段的长度。
对于现在要考察的高维欧几里得空间Rn，两个n维向量之间的距离可如法炮制地定义。对于二维平面R2，设第一个行向量为x = (a, b)，第二个行向量是y = (c, d)，则它们之间的距离就是平面上这两点之间的通常距离
它也就是空间向量x – y的实际长度。
三维以上的欧几里得空间，我们的眼力无法感知，但我们的想象力可以穿越界限，所以将上面直观的距离公式推广开来，就得到n维空间Rn的欧几里得距离定义：设向量x的分量为a1, a2, …, an，向量y的分量为b1, b2, …, bn，则它们之间的距离是
d(x, y) = {(a1 - b1)2 + … + (an - bn)2}1/2。
我们将三维以上空间中的向量“长度”改称为“欧几里得范数”，简称范数或2-范数，则n维向量x的范数||x||2被定义为x的所有分量的平方之和再开平方根。这样，x和y之间的距离就等于向量x – y的欧几里得范数||x – y||2。
作为长度概念的推广，范数继承了长度的三个基本性质：
1. 向量x的范数||x||2总是非负实数，且||x||2 = 0当且仅当x是零向量，即它的每个分量都为0。
2. 一个数α乘以向量x得到的向量αx，其范数等于α的绝对值乘以x的范数，即||αx||2 = |α| ||x||2。这里αx的各个分量等于α乘以x的相应分量。
3. 向量x和向量y的和的范数小于或等于向量x的范数加上向量y的范数，即||x + y||2 ≤ ||x||2 + ||y||2。这里x + y的每个分量为x和y的对应分量之和。
性质(iii)中的不等式一般被称为三角不等式，因为在二维和三维的情形，它就是平面几何中的一条定理：三角形的任意一边的长度不大于其余两边的长度之和。
类似于关于绝对值的三角不等式，范数的三角不等式可以导出不等式| ||x||2 - ||y||2 | ≤ ||x - y||2。由此可知将x映到 ||x||2的范数函数|| ||2是定义域Rn上的一个连续函数。
有了n维欧几里得空间Rn的向量范数概念，定义在Rn上的多变量线性向量函数L可以按照线性代数或泛函分析的习惯叫法被称为仿射算子；特别地，当常向量c = 0时，L叫做线性算子。这里，“线性”意味着对所有的数α和β以及向量x和y，恒有等式L(αx + βy) = αL(x) + βL(y)，而“仿射”则意味着对所有的满足α + β = 1的数α和β以及任意向量x和y，恒有等式L(αx + βy) = αL(x) + βL(y)。这样，我们可以用上由向量范数所诱导出的矩阵范数，结合著名的“压缩映射定理”来检验线性迭代法的收敛性。为此目的，我们将Rn中所有其欧几里得范数不大于1的那些向量组成的集合称为Rn的n维单位闭球，它的“表面”，即所有范数等于1的向量全体称为(n - 1)维单位球面，用Sn-1表示。
给定一个n阶方阵M，让x取遍Rn的(n – 1)维单位球面中的所有向量，则定义了一个非负实函数，它把x映到||Mx||2。由于线性运算x → Mx和范数运算u → ||u||2都是连续的，它们的复合函数||Mx||2在定义域Sn-1上是处处连续的。又因为有限维的单位球面不仅是有界集合（因为集合中所有元素的范数一致有界；1就是一个上界），而且是闭集（因为只要Sn-1中的任意一个向量序列{xk}当k趋向于无穷大时收敛到向量x，则在恒等式||xk||2 ≡ 1中两边取k趋向于无穷大时的极限，考虑到范数的连续性，就有||x||2 = 1，即x也属于Sn-1），因此根据微积分中的一条定理：定义在欧几里得空间的有界闭集上的连续函数一定有最大值，我们得出结论：非负函数||Mx||2在Sn-1上的某一点取到最大值。没有运气学过连续函数最大值定理的读者肯定在中学学过代数，可以通过画出站在有界闭区间[-1, 1]上的一条包含两个端点的抛物线，马上感知：此连续曲线上有一点最高（当然也有一点最低）。
这个最大值就被定义为由Rn上的欧几里得范数所诱导出的矩阵范数||M||2。从矩阵运算的性质马上得知如下关于矩阵、向量范数的常用不等式
||Mx||2 ≤ ||M||2 ||x||2 对所有的n维列向量x都成立。
另外，被向量范数诱导出的矩阵范数也和向量范数一样满足上述三项“范数”的基本性质，即||M||2是非负实数，它的值为零当且仅当M是零矩阵；一个数α和M的数乘αM的范数等于α的绝对值乘以M的范数；两个同阶矩阵之和的范数不大于它们各自的范数之和。由于两个同阶方阵可以相乘，矩阵范数比向量范数多了一个基本性质：两个同阶方阵之积的范数不大于它们各自的范数之积。特别地，||Mk||2 ≤ (||M||2)k。
是时候寻找巴拿赫不动点定理的新应用了。对于线性迭代向量函数L(x) = Mx + c，我们有
||L(x) - L(y)||2 = ||Mx + c - (My + c)||2 = ||M(x - y)||2 ≤ ||M||2 ||x - y||2。
如果矩阵M满足严格不等式
||M||2 < 1，
那么，L在完备的距离空间(Rn, d)上是压缩的；这里的距离d(x, y) = ||x – y||2。空间Rn的“完备性”是由实数集R的完备性继承而来。根据巴拿赫不动点定理，方程x = L(x)有唯一的不动点列向量p，并且从任一个初始列向量x0出发，迭代点列xk = L(xk-1), k = 1, 2, …当k趋向于无穷大时将收敛到p。这就是线性迭代法的一个收敛性条件。自然，这仅仅是一个充分条件，不是必要条件。在后续的一篇文章里，我们将想方设法推导出线性迭代收敛的一个既充分又必要的条件。
我们将如上的充分条件用到迭代求解线性方程组Ax = b，其中的方阵A假定为非奇异的，这保证了对任意右端列向量b，该方程组有唯一解。先解释什么叫非奇异矩阵。一个有n行和n列的矩阵M可以被看成是定义在n维欧几里得空间Rn上的一个线性算子，它将每一个n维列向量x映到n维列向量Mx。自然它将n维零向量0映到它自己（故向量0是M的一个不动点）。如果M不能将一个非零列向量（即该向量至少有一个分量为非零数）映成向量0，或等价地说，如果M将所有的非零向量映成非零的向量，则我们称M为一个非奇异矩阵。非奇异矩阵又被称为可逆矩阵，因为一个方阵是非奇异的当且仅当它一一对应地将Rn中的每一个元素映到Rn中的每一个元素。换言之，可逆矩阵不仅是单射的（即不同的向量被映到不同的向量），而且是满射的（即Rn中的每一个列向量都是该矩阵乘以某一个列向量的结果）。
如果n阶方阵M是非奇异的，那么它是单射和满射的。单射和满射这两个性质同时成立时就叫做是双射的，因而M作为定义域为Rn的线性算子，它具有所谓的逆算子，其定义域是M的值域Rn。故存在一个与M同阶的方阵，称为M的逆矩阵，记为M-1，满足条件：对任意的n维列向量x和y，等式Mx = y成立当且仅当等式M-1y = x成立。如果把算子M视为一个函数，那么逆矩阵实质上无异于中学代数课本里讲到的反函数概念。如此一来，可逆矩阵和它的逆矩阵同时满足两个等式MM-1 = I及M-1M = I，这和可逆函数f及其反函数f-1满足的两个恒等式f-1(f(x)) ≡ x，f(f-1(y)) ≡ y一模一样，其中x取f的定义域中的所有元素，y取f的值域中的所有元素。这些都是最简单的算术等式aa-1 = a-1a = 1(a ≠ 0)的自然推广。我们也顺便指出，从矩阵理论可知，上述关于矩阵和逆矩阵的两个等式并非彼此独立，其中的一个成立就能推出另一个也成立。如果读者熟悉行列式的概念，可能也会知道方阵为非奇异的一个简单的等价条件：矩阵的行列式不等于0。
有了上面的准备知识，我们可以回到设计迭代法求解具有标准形式的线性代数方程组Ax = b这个文章主题。第一个需要考虑的问题是怎样将此向量方程写成迭代法的标准形式x = L(x)。方法很简单，只要将矩阵A分解为两个矩阵之差就行：A = N - P，但需要一个额外条件，那就是矩阵N必须是非奇异的。然后原方程Ax = b就可写成Nx = Px + b。由于N是可逆矩阵，前述方程等价于下列不动点方程
x = N-1Px + N-1b = Mx + c，
其中M = N-1P和c = N-1b。因此，对应的求解线性方程组Ax = b的迭代公式是
xk = N-1Pxk-1 + N-1b, k = 1, 2, 3, …; 给定x0。
读者自然会问，怎么选取关键的矩阵N？这个矩阵要有逆矩阵，要保证迭代法收敛，并且是收敛得越快越好。这是计算数学家们一直追求的目标。我们现在潜入历史，看看有哪几位数学家在线性迭代法方面贡献卓越，成为先驱。
历史上第一个提出迭代法的数学家又是高斯，所以他不仅是以高斯消去法为代表的直接法的鼻祖，也是迭代法的鼻祖。1823年12月26日，高斯在给他昔日博士生戈林 (Christian Ludwig Gerling，1788-1864) 的一封信中，提出了这个他称之为“非直接消去法”的迭代法，用于求解关于四个城市的一个应用问题。他首先利用也是他发明的最小二乘法得到一组线性方程，然后用自己设计出的迭代法求解对应的法方程。这走出了迭代法求解线性方程组历史进程的第一步。然而，由于高斯的许多数学发现和发明并没有被他正式发表，所以他的发明权不得不常与后人分享，他也不在乎，没有当今“不发表便灭亡”的“科学家烦恼”。半个世纪后，他的同胞数学家冯·赛德尔 (Phillip Ludwig von Seidel，1821-1896) 于1874年发表了包含同样思想的迭代法。
1845年，另一名德国数学家雅可比 (Carl Gustov Jacob Jacobi，1804-1851) 也提出了一个以他的名字命名的迭代法。这个方法比较简单，在数值代数的教科书中一般是作为迭代法求解线性方程组的第一个例子被介绍给学生。雅可比的人生轨迹虽然只持续了不到四十七年，他却是一位多产的数学大师，其最伟大的贡献在椭圆函数方面。他被上世纪的美国数学史家贝尔 (Eric Temple Bell，1883-1960) 写进影响了几代读者的《大数学家》 (Men of Mathematics)，该书总共只描述了到上世纪初为止留名数学史的三十多个“大数学家”。雅可比的职业生涯中有十三年是在柯尼斯堡大学担任正教授，这所大学的杰出校友有数学家希尔伯特 (David Hilbert，1862-1943) 和闵可夫斯基 (Hermann Minkowski，1864-1909)。
求解线性方程组历史上最有名气也最简单实用的两个迭代法都是德国人的杰作，今日关于它们的基本理论也很漂亮。但是本文已够长了，所以我将再写一文，专门介绍雅可比迭代法和高斯-赛德尔迭代法的基本思想和收敛性。
写于2023年11月3日星期五美国哈蒂斯堡夏日山庄
本文受科普中国·星空计划项目扶持
出品：中国科协科普部
监制：中国科学技术出版社有限公司、北京中科星河文化传媒有限公司
相关阅读
近期推荐
4 有机化学界的一场“华山论剑”，结下了丰硕的科学和精神果实
特 别 提 示
1. 进入『返朴』微信公众号底部菜单“精品专栏“，可查阅不同主题系列科普文章。
2. 『返朴』提供按月检索文章功能。关注公众号，回复四位数组成的年份+月份，如“1903”，可获取2019年3月的文章索引，以此类推。
版权说明：欢迎个人转发，任何形式的媒体或机构未经授权，不得转载和摘编。转载授权请在「返朴」微信公众号内联系后台。
长按下方图片关注「返朴」，查看更多历史文章