贝斯特全球最豪华最新版3个德国人创造的线性迭代法超越了一个时代

　　本文为澎湃号作者或机构在澎湃新闻上传并发布■◆★◆■，仅代表该作者或机构观点，不代表澎湃新闻的观点或立场，澎湃新闻仅提供信息发布平台。申请澎湃号请用电脑访问。

　　用D表示由A中的对角元素a11■◆★★, a22■■★■★★, …■★★, ann保持原先位置和次序而得到的相应对角矩阵★◆，L表示A中严格位于对角线之下的所有元素aij (i j)按照原先位置而构成的严格下三角矩阵◆★◆★，U表示A中严格位于对角线之上的所有元素aij (i j)按照原先位置而构成的严格上三角矩阵。这里习惯采用的三个字母分别是三个英文单词diagonal（对角）、lower（下）和 upper（上）的大写首字母■◆■★★，一目了然。于是★■■★，A = D + L + U = N - P。请注意L和U的所有对角元素为0，此结构的一个推论是Ln = Un = 0，这将在后面用到◆◆★■★■。

　　根据如上的谱分析结论◆◆，迭代法对所有初始点收敛当且仅当ρ(N-1P) 1。

　　然而，有一个收敛的充分必要条件等待着我们去挖掘，这个条件不再用到矩阵的范数★■■■◆，而是矩阵的谱。谱给出了矩阵不依赖于矩阵范数的一个内蕴性质。方阵的谱是它所有特征值全体组成的一个有限的复数集合。给定n阶方阵M■◆◆★，我们称复数λ为M的一个特征值，如果复矩阵M - λI不是可逆的■★◆■，即存在一个复的非零列向量x使得Mx = λx。这个x被称为M对应于特征值λ的一个特征向量。我们曾经提到■★★■■，一个矩阵是奇异的当且仅当它的行列式等于零，故λ是M的特征值当且仅当det(M - λI) = 0★■■，其中符号det表示行列式。如果把这个等式左边中的λ看成是变元，根据行列式的定义★■◆★，det(M - λI)的展开结果是关于λ的一个n阶多项式，所以一个n阶方阵M顶多有n个相异的特征值◆■◆★。我们把M的所有特征值绝对值中的最大值称作 M的谱半径★◆，记为ρ(M)。

　　为了证明这个断言◆◆◆◆★★，先设迭代法对所有的初始列向量x0都收敛，则极限列向量满足不动点方程x = Mx + c◆■■。又因为如上所述I - M是非奇异矩阵，这个极限向量必定是F唯一的不动点p◆★★■■★。令ek = xk - p为迭代到第k步时的误差向量，那么对所有的列向量e0，极限 ek = Mek-1 = Mke0 → 0都成立。如果ρ(M) ≥ 1■■■，则存在一个特征值λ满足λ ≥ 1★◆◆■■，也存在它所对应的一个特征向量e0◆★◆◆★。因为e0是非零向量◆★★■，然后我们就会发现■◆■★◆★，ek = Mke0 = λke0当k趋向于无穷大时不可能趋向于0，这与假设矛盾。所以ρ(M) 1■★，必要性得证。

　　我们把高斯-赛德尔方法另一个独特的收敛性命题端上餐桌■★■，作为本文的最后一道数学菜肴：

　　在证明它前，我们指出，对实对称正定矩阵A和所有复向量x ≠ 0，也有xHAx 0★★■◆，其中xH是x的共轭转置。在此强调之★★◆★◆，是因为后面的特征向量一般是复的。在下列的推理中，我们常用到简单事实■◆■◆◆■：对两个可乘的复矩阵B和C◆■★◆，有(BC)H = CHBH；对任意复向量x★■★◆，二次型xHBx是个复数；对实矩阵B★◆★，有BH = BT；对一阶复矩阵[b]，有

　　自然，无论是雅可比迭代法还是高斯-赛德尔迭代法■★★■，对矩阵A缺乏某些额外的假设条件，都不能保证其收敛性■★。那么，当“A的对角线元素个个不能为零”这个先决条件被满足后◆★，什么是保证它们收敛的充分条件呢？

　　我们将雅可比方法和高斯-赛德尔方法的迭代矩阵分别记为M1 = I - D-1A = -D-1(L + U) = -D-1L - D-1U和M2 = -(D + L)-1U = -(I + D-1L)-1D-1U。因为M1的每行所有元素之和都不大于它们中的最大值M1∞，而一个矩阵乘以特殊列向量e的结果就是将该矩阵每行加起来的和所构成的列向量，故有向量不等式

　　这里，保证雅可比迭代法收敛性的严格对角占优的要求可以被改成不可约弱对角占优，即上面的所有严格不等式换为一般不等式 “≥”，且至少有一个不等式严格成立★★★■，但增加了矩阵不可约的附加条件。如果一个方阵的行和列能以同样的方式重新排列（不再排列也算“重新排列”◆◆★，因为单位矩阵也是排列矩阵）★■，结果变成一个2×2块上三角矩阵，即矩阵的左下角有至少为1×1的0矩阵，那么我们说这个矩阵是可约的，否则就说它是不可约的。

　　之一满足，则雅可比迭代法对所有的初始列向量都收敛到Ax = b的唯一解p = A-1b。

　　我在前一篇《返朴》文章《》中，为了能与适用于一般完备距离空间的压缩映射定理挂上钩★■■，只给出了迭代法求解线性方程组的一个简单的收敛性充分条件■★★，即若要迭代格式xk = Mxk-1 + c, k = 1, 2, 3★★◆■◆■, …对所有的初始列向量x0都收敛◆■◆■◆★，一个对迭代矩阵简单易懂的要求是：Rn上的向量2-范数所诱导出的矩阵2-范数M2小于1。这时，迭代向量序列{xk}当k趋向于无穷大时将收敛到线性方程组x = Mx + c的唯一解■■★◆★。在上述假设下，一个直接的推论是I - M为一可逆矩阵。

　　这样一来，n维向量空间Rn在随便哪个范数•下都变成一个完备的距离空间。给定一个n阶的方阵M，这个被选用的向量范数诱导出对应的矩阵范数M，它等于连续函数Mx在(n – 1)维单位球面上的最大值。对同一个矩阵，不同的向量范数诱导出不同的矩阵范数，因而如果某一个向量范数诱导出的矩阵范数小于1，那么根据巴拿赫不动点定理，线性迭代法xk = Mxk-1 + c★■◆, k = 1, 2★■◆◆, 3, …对任一初始点x0都收敛到同一个极限。

　　有两个对于计算矩阵范数而言是最方便的向量范数，其中之一是上面的q-范数取q = 1时的1-范数，即n维向量x = (a1★★★◆, a2, …★■, an)T的1-范数为

　　至于高斯-赛德尔迭代法★◆◆■■■，我们下面花点笔墨，证明在A为严格对角占优的假设下◆★■■◆，它的迭代矩阵-(D + L)-1U的∞-范数不大于雅可比迭代矩阵I - D-1A的∞-范数。这个证明主要基于“序 (order) ◆★★★★”的概念。法国布尔巴基数学学派认为纯粹数学主要研究三种结构■★：关于代数运算的代数结构、关于连续性质的拓扑结构、关于大小关系的有序结构。两个同一尺寸的实矩阵B和C若被写成B ≤ C◆◆■★■，指的是B的每一个元素bij小于或等于C的对应元素cij◆◆■◆★■。对于B■★★，我们记B为它的绝对值矩阵，即B的每一个元素为B的对应元素bij的绝对值bij。注意，B是非负矩阵，而B是非负数。另外，我们用字母e表示所有分量均为1的任意维数的列向量★★◆■■■。

　　上述条件(ii)可以用来定义一类矩阵。一个n阶方阵A被称为严格对角占优的条件是■◆，对所有的i = 1, 2, …◆◆★■★■, n■◆◆■★★，严格不等式

　　德国数学家在两百年前点燃的线性迭代法火炬，百年之后远渡重洋，传递到了科技腾飞、电脑出世的美洲大陆。为了优化迭代效率，从经典的高斯-赛德尔方法出发◆★■，在1950年，美国数学家贝斯特全球最豪华最新版、计算机科学家杨大卫 (David M. Young Jr.★★■◆，1923-2008，“返朴”将另行介绍他) 和美国物理学家、计算机科学家弗兰克尔 (Stanley Phillips Frankel，1919-1978) 几乎同时引入了一个松弛因子ω进行某种仿射组合◆◆◆，引出了“逐次超松弛迭代法 (method of successive over-relaxation◆◆◆，SOR) ★◆■★”■★。前者在哈佛大学数学系的博士学位论文Iterative Methods for Solving Partial Difference Equations of Elliptic Type（《求解椭圆型偏差分方程的迭代方法》）中提出了SOR方法；后者在美国数学会期刊（十年后改名为Mathematics of Computation[《计算数学》]）上发表论文Convergence rates of iterative treatments of partial differential equations（《偏微分方程迭代处理的收敛率》），其中对他设计并称之为“Extrapolated Liebmann method（外推利伯曼法）的SOR方法进行了全面的分析。这两项提出SOR方法的先驱性研究，都与在科学和工程中大量出现的偏微分方程有关，用电子计算机求解这些连续方程离散化后的大型稀疏线性方程组，迭代法是首要之选，这就是SOR方法应运而生的时代背景■★◆。

　　当松弛因子ω取值为1时，SOR方法就还原成高斯-赛德尔方法。我们不细表此法，就此结束本文，因为我们已经通过了解两个古典的迭代法而知悉线性迭代法收敛性准则的基本思想■■◆。总而言之◆★◆■■，十九世纪三名德国数学家的发明创造◆◆■■■★，是现代大型、稀疏、结构化线性方程组迭代解研究浪潮的源头■■■★◆。

　　范数的等价性有什么妙用呢■★◆◆★★？答案是它们各自定义出的序列收敛性彼此等价，也就是说Rn中的向量序列{xk}在由范数•所定义的距离下收敛于一个向量，当且仅当它在由范数•所定义的距离下收敛于同一个向量。我们给出这个事实的证明◆◆■★★★，因为它不难：设当k趋向于无穷大时xk - x → 0■★■◆■，则xk - x ≤ β xk - x → 0；反之■◆■◆■★，若xk - x → 0，则xk - x ≤ (1/α) xk - x → 0。顺便提及，该性质是有穷维空间的特征，即一旦线性空间的维数变成无穷大，那么赋予它的两个范数不一定导出同等的收敛性◆◆。所以，从★★“有穷”到★◆“无穷”■◆★，它们之间常有“万里之遥”★★。

　　我们想要提醒读者注意的是，上面的收敛性谱半径等价条件为真的一个前提是逆矩阵(I - M)-1存在，否则的话，反例多得很，比如令M为一的对角矩阵，其对角元素为1和1/2，则它的谱半径不小于1。易知I - M是一个对角线的奇异矩阵■◆■★■。另一方面，如果让c = (0◆■◆★, 1)T，则对任意的初始向量x0 = (a★★◆★■■, b)T★◆◆■，迭代xk = Mxk-1 + c■■■★◆★, k = 1■◆★◆, 2◆★◆★, …产生的向量序列xk = (a, b/2k + 2 - 1/2k-1)T → (a, 2)T，说明这个迭代对所有的初始向量都收敛■★★，但极限向量不唯一， 它们构成xy-平面内的一维仿射集{(a, 2)T: a是任意实数}★★★。

　　如果雅可比迭代法的收敛条件(ii)成立，则高斯-赛德尔迭代法对任意的初始列向量都收敛到Ax = b的唯一解p = A-1b。

　　反过来，如果ρ(M) 1，则存在一个足够小的正数ε 使得ρ(M) + ε 1。这一看似简单的事实是实数的基本性质。另一方面，对于M的每一个特征值λ和Rn上的任一个范数•，设x为λ所对应的一个特征向量，由不等式M x ≥ Mx = λx = λ x以及x ≠ 0可知，刚刚获得的不等式两边可以除以正数x◆■◆■◆■，而得到不等式M ≥ λ，从而可知：矩阵M的谱半径ρ(M)是所有向量范数诱导出的矩阵范数在M的值M的一个下界。事实上，这个下界是对应于固定矩阵M的全部范数值M构成的实数集合的“下确界”，即这个集合的所有下界中的“最大下界”。“上、下确界”的概念是微积分学的最基础★◆★、最重要的概念，真正学通了它，极限理论以及随之而来的连续性◆■■■、导数、积分、级数等等概念学起来甚至可以“无师自通”了。在稍微高等一点的矩阵理论中★◆■◆◆，可以利用M的★■“若尔当标准型”◆★◆，构造Rn上的一个范数•，使得它诱导出的矩阵范数满足条件M ≤ ρ(M) + ε 1（这个对任意正数ε的矩阵范数构造，恰恰证明了上述的“谱半径是所有矩阵范数值的下确界”断言）。故对任何初始列向量e0有

　　这些范数都满足前一篇文章《》里列出的欧几里得范数符合的那三项■★◆◆“范数公理■■”(i)-(iii)。这样◆★★★■，Rn就有无穷多个q-范数，此外还有无穷多个不属于这个范畴的其他范数。如此数目众多的范数，会不会令人看得眼花缭乱？不同的范数会影响线性迭代向量序列的收敛性吗？幸好有一个关于“有穷维线性赋范空间”的基本事实让我们放心：Rn上的任意两个范数•和•都是所谓★■★★“彼此等价◆◆■”的范数，即存在两个正数α和β，使得对Rn中的所有向量x都有

　　好在我们还有其他的范数可用。作为通常直线段长度概念的直接推广，n维欧几里得空间Rn中向量的2-范数只是一种范数，在多元微积分里用得最多，但它是更一般的q-范数取q = 2时的特例★★■◆。选取一个大于或等于1的实数q，对于分量为a1★■, a2, …, an的任意n维向量x，它的q-范数被定义为

　　因为上述等式的右端是正数◆◆★，所以左端的1 - λ2 0，即-λ = λ 1★■◆■。由于-λ是M2的任一特征值，这就证明了ρ(M2) 1。迭代矩阵谱半径小于1等价于线性迭代法收敛，命题得证★◆■■◆。

　　可以验证，但细节比较繁琐，不在这里表出，Rn上的向量1-范数所诱导出的n阶方阵M = [mij] 的1-范数为

　　如果线性方程组Ax = b的系数矩阵A是对称正定的，那么高斯-赛德尔迭代法收敛◆■■■★。

　　研究线性迭代的主要目的是数值求解线性方程组Ax = b，其中的系数矩阵A为非奇异的，这样保证对所有的右端常向量b★■■，该方程组有并且仅有一个解，它就是p = A-1b■◆★◆◆。为了设计一个迭代法■★■★■◆，首先将A分裂成N - P的形式，其中N也是非奇异的◆■。然后原方程组等价于不动点线性方程组x = Mx + c，其中M = N-1P和c = N-1b■◆◆◆■。基本的A非奇异性假设推出I - M也是非奇异的，故由F(x) = Mx + c定义的仿射向量函数F◆■★◆◆★: Rn → Rn有并且仅有一个不动点向量p。现在我们可以证明：线性迭代

　　当A为严格对角占优矩阵时◆◆★■★，高斯-赛德尔迭代矩阵的∞-范数不大于雅可比迭代矩阵的∞-范数，而这个范数越小■◆◆■，收敛速度通常就越快，因而就收敛快慢而言，前一种方法不亚于后一种方法。

　　我们试图写出上述两个迭代法的分量迭代公式。对i◆■◆◆★, j = 1, 2◆◆, …■◆★◆★, n■★◆■, 记A的第i行、第j列元素为aij★★■。为了避免矩阵求逆运算，将雅可比迭代的循环公式写成如下形式

　　今年12月26日是德国数学家高斯发明历史上第一个线周年■★◆■◆★。在前不久发表的《》的一文中，我们从基础概念出发，厘清了迭代法求解线性方程组的准备知识和具体过程◆◆■■。这次我们将探讨历史上最有名气、也最简单实用的两种一阶定常迭代法：雅可比迭代法和高斯-赛德尔迭代法的基本思想和收敛性★■。

　　所以★◆■◆◆■，只要迭代矩阵M满足M1 1或M∞ 1，则迭代法对任意的初始向量都收敛■★。

　　都成立。由此可见★◆■★，如果线性方程组Ax = b中的矩阵A是严格对角占优的，那么雅可比迭代法收敛★◆■★◆。这时，A的非奇异性自动由严格对角占优的假设保证◆◆★◆★，因为A = D[I - (I - D-1A)]是两个非奇异矩阵的乘积，后一因子矩阵I - (I - D-1A)的非奇异性是由于，从不等式ρ(I - D-1A) ≤ I - D-1A∞ 1可以推出：1非I - D-1A的特征值◆■◆，故而I - (I - D-1A)没有特征值0。此外还可证明，若严格对角占优矩阵A是实对称（即AT = A）的并且对角元素均为正数，那么A必定是正定矩阵，即对所有的n维非零列向量x■★■，严格不等式xTAx 0都成立★■■。

　　我们已经知道★★◆★，它们收敛的充分必要条件分别是ρ(I - D-1A) 1和ρ[-(D + L)-1U] 1。然而，矩阵的谱半径不是那么好得到的，因为它用到所有特征值★■。但是矩阵的1-范数和∞-范数的计算可以说是举手之劳。由于矩阵I - D-1A的所有对角元素均是0，我们马上发现它的1-范数等于

　　现在我们可以集中讨论求解线性方程组Ax = b的雅可比迭代法和高斯-赛德尔迭代法了，其中A为n阶可逆方阵■◆★★。这两种产自德意志的经典迭代法都是基于在矩阵分解A = N - P中对N和P的特殊选取，其迭代格式都可写为

　　。谱半径是矩阵固有的内蕴性质，而范数却是“强加于★★■★”矩阵身上的一个外在标尺，通过极限之桥梁，本质透过外观而表达。更神的是，这个关系式与维数无关，对无穷维巴拿赫空间上的有界线性算子同样正确。

　　另一方面★■★，由于-D-1L和-D-1U分别为严格下三角矩阵和严格上三角矩阵，显然有等式

　　另外可以证明，迭代法对所有初始向量都收敛的另一个充分必要条件是：矩阵幂次序列{Mk}当k趋向于无穷大时收敛到零矩阵★■◆◆■★，但这个等价条件远不及迭代矩阵谱半径小于1的条件简单实用。

　　因为A 为对称正定矩阵★★◆，U = LT。令-λ为高斯-赛德尔迭代矩阵M2 = -(L + D)-1LT的一个复特征值，且x为对应的一个复特征向量，则LTx = λ(L + D)x，因而

　　我们比较一下雅可比方法与高斯-赛德尔方法在迭代点分量计算过程中的显著差别★■■◆。由雅可比迭代法中n个分量的迭代公式可见■◆★■◆，在第k - 1个迭代向量所有的分量都计算完毕后，第k个迭代向量的各个分量才一个接着一个地计算出来；换句线步迭代获得的全部分量都被用于对第k步迭代所有分量的计算。而对于高斯-赛德尔迭代法，每当需要算出第k个迭代向量的第i个分量时■◆，不仅需要已完成计算的第k - 1个迭代向量从第i + 1到第n个分量帮忙，而且还需要正在进行中的本次迭代已得到的第k个迭代向量的第1到第i - 1个分量相助；或言之，高斯-赛德尔法比雅可比法“更急于求成”■◆■，命令当前迭代步骤中刚刚收入囊中的迭代点分量取代上一迭代步骤被■★★★■“擒获”的同下标分量，◆◆★■◆★“披挂上阵★★◆”◆■◆。这说明，高斯和赛德尔比雅可比更会“用兵◆★★★◆◆”，不肯浪费半点士气。

　　版权说明■★◆◆：欢迎个人转发■★■◆◆◆，任何形式的媒体或机构未经授权，不得转载和摘编。转载授权请在「返朴」微信公众号内联系后台■■★。

　　然而，这里又产生一个技术问题★★■◆■★：M的欧几里得2-范数极难计算，至少它不能从M各行各列的具体元素中立即通过简单代数运算算出。事实上★★◆■，线性代数中的二次型理论告诉我们，实数方阵M的2-范数等于M的转置矩阵MT与M的乘积MTM这个所谓“半正定矩阵”（意指二次型xTMTMx = （Mx2）2 对所有的实列向量x都是非负实数）的最大特征值之平方根（因为MTM的所有特征值均为非负数◆■，故平方根存在）。上句话里包含了好几个数学概念，可想而知计算出M2绝对不是一件轻而易举的事情◆★★◆★。