把不收敛的迭代拉回安全区：M矩阵与校正矩阵实战

做迭代求解的人大概都经历过这种场面残差曲线好不容易下来一段突然卡在某个平台期开始震荡左调右调松弛因子都没用。后来我养成一个习惯遇到不收敛先不碰参数直接判断系数矩阵的结构。见过太多次问题根源不是算法而是矩阵里藏着几个正的非对角元——它们让整个迭代过程失去了收敛保证。这时候真正该做的是从收敛条件出发把矩阵修成M矩阵也就是俗称的校正矩阵策略。这篇文章是写给两类人的一类是和当年的我一样被不收敛问题折磨得睡不着觉的求解器开发者和数值计算工程师另一类是学过数值代数但一直没想明白M矩阵到底有什么用的人。我尽量不用教科书式的写法而是把我反复用过的一套逻辑讲清楚先看收敛条件的本质是什么再把它翻译成矩阵元素层面的不等式最后根据这个不等式去构造校正矩阵。中间会用一个3x3的反例手推一遍完整流程再用对流占优方程的离散矩阵做一次工程实战最后聊聊我踩过的几个坑。1. 先搞懂迭代收敛到底卡在矩阵的哪个部位1.1 谱半径条件一句话版本的收敛判据几乎所有线性方程组的迭代解法最后都能写成同一个形式x^{k1} T x^k c这里的T叫迭代矩阵。比如Jacobi迭代里假设把系数矩阵A拆成对角部分D和非对角部分LU那T就是I - D^{-1}A。Gauss-Seidel和SOR本质上也是换了一个T只是换的方式不同它们背后的判据完全一样。误差怎么变设每次迭代的误差是e^k x^k - x^*代入迭代式很容易得到e^{k1} T e^k也就是说每一步迭代误差向量都会被同一个矩阵T乘一次。那收敛不收敛就看这个T的谱半径——也就是所有特征值模的最大值——是不是小于1。如果ρ(T) 1误差就会按几何级数衰减只要有一个特征值模超过了1迭代必然发散。用个不太严谨但很好记的类比把谱半径看成每次迭代误差的最大放大倍数。谱半径小于1相当于存款利率是负的钱越存越少谱半径大于1相当于银行在按复利抽你的钱越滚越多最后必然爆掉。这里有个工程上容易被忽略的点很多人习惯用残差的绝对值有没有下降来判断迭代是否可行这在接近收敛边界时是很不可靠的。某些误差分量可能恰好落在特征值模接近1的子空间里残差曲线会先骗你下降一段时间等那个慢分量主导了才原形毕露。所以判断收敛性最靠谱的还是把谱半径的量级估出来而不是凭残差曲线的前半段感觉。1.2 M矩阵为什么拥有免检资格M矩阵这个概念在数值代数里被神化了但真正搞懂它为什么有用的人不多。一个非奇异的M矩阵最简单的定义有两部分一是非对角元全部小于等于0这样的矩阵叫做Z矩阵二是它的逆矩阵所有元素都大于等于0。等价地M矩阵可以写成A sI - B的形式其中B的每个元素都非负且s大于B的谱半径。这个写法看起来抽象但它的真正威力藏在下面这条定理里如果A是非奇异M矩阵那么对它做任意正则分裂A M - N其中M可逆、M^{-1}非负、N非负迭代矩阵T M^{-1}N的谱半径必然小于1迭代必然收敛。反过来也成立。这意味着什么意味着只要矩阵是M矩阵你真的不需要费心去挑分裂方式Jacobi、Gauss-Seidel、各种块迭代基本都能收敛。这在工程上是一个极强的好消息。你会发现椭圆型偏微分方程用标准差分格式离散后比如五点格式的Laplace矩阵行和是零、非对角元为负、对角元为正天然就是M矩阵。所以传统扩散问题的求解器很少出现结构型发散。问题往往出在那些偏离标准结构的场景强对流项、非物理数值格式、大变形网格这些场景会把矩阵从M矩阵的安全区推出去迭代发不发散就成了碰运气的事。这时候自然会产生一个想法既然M矩阵有免检资格那能不能把给定的坏矩阵通过一个修正项重新拉回这个安全区这就是校正矩阵策略的出发点。2. 从谱半径反推校正量按需加药的构造路线2.1 把收敛条件翻译成矩阵元素不等式前面说的谱半径判据虽然精准但没法直接用来施工——你不能对着一个谱半径大于1的结论反推该改哪个元素。所以工程上需要一个中间步骤把ρ(T) 1这个条件松弛成一个可直接验证、可逐行修补的充分条件。以Jacobi迭代为例T I - D^{-1}A。把T写开它的第i行第j列元素是-a_{ij}/a_{ii}其中j不等于i。根据Gershgorin圆盘定理T的所有特征值都会落在以0为中心、以第i行绝对值之和为半径的一族圆盘里。每个圆盘半径是r_i Σ_{j≠i} |a_{ij}| / |a_{ii}|如果对每一行来说r_i都小于1那么所有Gershgorin圆盘都落在单位圆内谱半径自然小于1。而这个条件翻译回来就是严格对角占优条件对于任意一行i非对角元绝对值之和小于对角元的绝对值。这一步翻译极其关键。它把一个需要知道所有特征值才能验证的判据变成了一个逐行检查的数值不等式。而且它的逻辑是充分不必要——满足一定收敛不满足未必发散。但从构造校正矩阵的角度我们要的就是这个充分条件只要把矩阵修成严格对角占优的Z矩阵收敛性就有了保证。2.2 正的非对角元一切不稳定的源头现在看这个不等式左边是非对角元绝对值之和右边是对角元绝对值。如果矩阵的非对角元都是负的那么它们以绝对值的形式贡献到左边和对角元做对抗物理上可以理解成相邻节点之间的耗散项。一旦某个非对角元变成了正数它的绝对值同样贡献到左边但本质上它不再消耗能量反而在相邻节点之间引入了一股反向输送。在偏微分方程离散的语境里正的非对角元常常对应着数值格式允许信息从下游传到上游这本身就是违反物理常识的。在纯代数语境里它让矩阵丢失了Z矩阵的性质也就失去了成为M矩阵的资格。更麻烦的是它还会让某些迭代矩阵的特征值跑到单位圆外。处置一个正的非对角元a_{ij} 0最自然的办法不是直接把它抹掉而是把它吸收进第i行和第j行的对角元。具体做法是对每对i、j满足a_{ij} 0选定一个权重α ∈ [0, 1]然后执行a_{ii} α * a_{ij} a_{jj} (1 - α) * a_{ij} a_{ij} 0这样正的非对角元消失了它的能量被分配到两个相关的对角元上。这步操作不会改变方程组的整体平衡但会让矩阵恢复Z矩阵结构。α这个权重是留给用户的手术自由度如果第i行更危险非对角元已经快吞掉对角元了就把α调大一点让这一行分得更多的对角增量。如果矩阵原本是对称的、希望保持对称那就固定取α 1/2。2.3 校正矩阵的两种模式符号吸收与对角增补如果把整个构造过程拆成两个动作那就是下面两种模式模式修改对象解决的病症本质符号吸收正的非对角元Z矩阵条件不满足把违规的交叉项转移到对角元对角增补对角元本身对角占优不等式不满足直接增加对角元的安全余量这两步是顺序执行的。第一步先把矩阵里所有正的非对角元处理干净得到Z矩阵第二步再逐行检查严格对角占优条件哪一行不满足就在那一行补一个对角增量γ_i直到不等式成立。整个校正矩阵C就是这两步所有修改量的总和它记录了我们对原始矩阵动过的手术。后面会看到这个C不能丢它在不动点迭代里要显式参与计算否则改完矩阵却解错了方程。这套按需加药的思路我特意强调一下它和我们平时拍脑袋加大对角元做对角占优化有本质区别每一步增量都是从一个明确的不等式出发算出来的不是靠经验试出来的。所以整个流程的每一步都可以被审计、被复现出了问题也容易定位是哪个环节没算对。3. 一个3x3反例走完全流程发散矩阵如何变成收敛M矩阵3.1 病根诊断Jacobi迭代矩阵的谱半径为什么是2理论说再多不如一个手算的例子。看下面这个矩阵A [[1, 2, 0], [0, 1, 2], [2, 0, 1]]这个矩阵的符号模式一眼就能看出问题位置(1,2)、(2,3)、(3,1)上的元素都是正数指定不是Z矩阵更不可能是M矩阵。而且每行的非对角元绝对值之和都是2对角元是1严格的严格对角占优完全不满足。把它代进Jacobi迭代因为对角元全部是1D是单位矩阵迭代矩阵T I - A直接就是T [[0, -2, 0], [0, 0, -2], [-2, 0, 0]]用Gershgorin圆盘定理看每一行的半径都是2也就是说特征值必然落在半径2的大圆盘里。更直接一点这个T是三个循环置换的组合它的特征值是-2乘以三个立方根模长全是2。ρ(T) 2远大于1所以Jacobi迭代必然发散。在开始校正之前先把病根说透这不是松弛因子选得不好也不是初值给得差而是矩阵的结构本身不允许迭代收敛。后续所有操作都应该是针对这个结构的确定性修改而不是碰运气调参数。3.2 逐元吸收正元看矩阵如何被捋顺现在按第二套策略走一遍α统一取1/2。第一步处理(1,2)位置上的正元a_{12} 2。执行a_{11} 1 a_{22} 1 a_{12} 0这一步结束后矩阵变成[[2, 0, 0], [0, 2, 2], [2, 0, 1]]注意原来(2,3)上的2和(3,1)上的2还没有处理矩阵仍然是坏的。第二步处理(2,3)位置上的正元a_{23} 2a_{22} 1 a_{33} 1 a_{23} 0矩阵变成[[2, 0, 0], [0, 3, 0], [2, 0, 2]]第三步处理(3,1)位置上的正元a_{31} 2a_{33} 1 a_{11} 1 a_{31} 0最终得到A_tilde [[3, 0, 0], [0, 3, 0], [0, 0, 3]]一个原始发散矩阵经过三个等价的符号吸收步骤变成了3倍单位矩阵。它当然严格对角占优也当然是M矩阵收敛性没有任何悬念。那校正矩阵C是多少用A_tilde减去原始AC A_tilde - A [[2, -2, 0], [0, 2, -2], [-2, 0, 2]]可以看到C的构造完全符合预期对角上是正的增量非对角上是用来清零原正元的负量。3.3 验证与代码不动点迭代确实保持原解校正这一步做完最容易犯的错误是直接拿着A_tilde去解方程。记住我们想要的是Ax b的解而不是A_tilde x b的解。为了保持原方程的解不变正确的做法是用不动点迭代A x b等价改写为(A C) x C x b移项得到x^{k1} (A C)^{-1} (C x^k b)在这个3x3例子里A C 3I所以迭代化成极其简单的标量形式x^{k1} (C x^k b) / 3迭代矩阵是C/3。前面已经算过C的特征值模是2所以迭代矩阵的谱半径是2/3。也就是说这个原来发散的问题现在变成了每一步误差缩小到0.667倍的绝对收敛问题。下面是完整的验证代码可以复制到本地跑一下import numpy as np A np.array([[1., 2., 0.], [0., 1., 2.], [2., 0., 1.]]) b np.array([1., 0., 1.]) # 符号吸收正的非对角元按一半一半并入对角元 C np.zeros_like(A) for i in range(3): for j in range(3): if i ! j and A[i, j] 0: C[i, i] 0.5 * A[i, j] C[j, j] 0.5 * A[i, j] C[i, j] - A[i, j] A_tilde A C print(修正后的矩阵 AC:) print(A_tilde) # 不动点迭代求解 A x b x np.zeros(3) for k in range(30): x_next C x / 3.0 b / 3.0 res np.linalg.norm(A x_next - b) x x_next if k % 5 0: print(fk{k:2d}, 原方程残差{res:.3e}) print(迭代解:, np.round(x, 6)) print(参考解:, np.linalg.solve(A, b))这里的核心是把C保留下来放进迭代格式里而不是把C只当作修矩阵的一次性工具。迭代收敛后极限点满足(A C) x C x b两边消去C x项正好还原成Ax b所以解是原问题的精确解没有引入偏差。这一点是整个策略里最容易出bug的地方我在第5节还会单独展开讲。4. 工程实战对流占优离散矩阵的人工扩散校正4.1 中心差分矩阵里冒出来的正元局部Peclet数代数例子讲完了这个例子太小可能有人会觉得是雕虫小技。真正让这套策略发挥价值的地方在偏微分方程离散。拿一维稳态对流扩散方程来说-ε u w u f, x ∈ (0, 1)边界u(0) u(1) 0。如果对流速度w 0用中心差分离散对流项w u对网格步长h在内部节点i处得到下面这组系数u_{i-1}的系数-ε/h^2 - w/(2h) u_i的系数2ε/h^2 u_{i1}的系数-ε/h^2 w/(2h)前两项都是负的没有问题。问题出在第三项当w/(2h)超过ε/h^2时u_{i1}的系数就会变成正数。定义局部Peclet数Pe w h / (2ε)当Pe大于1时中心差分格式离散出的矩阵非对角元会出现正值偏离Z矩阵结构。这和前面3x3反例里看到的正元本质上是同一种病。在物理上它对应着一种违反直观的现象下游节点的信息被数值格式反向输送给了上游从而引发出解的振荡和迭代的不稳定。这个现象在CFD里被讨论了很多年但很少有人提它的代数本质。用本文的框架来看这就是一个典型的离散格式引入结构性病态的案例修法也不一定是换成更复杂的格式而是给矩阵做一次校正。4.2 人工扩散的本质一种结构化的校正矩阵处理对流占优问题的经典手段是迎风差分也就是把原来中心差分的对流项改成单侧差分。w 0时迎风格式得到u_{i-1}的系数-ε/h^2 - w/h u_i的系数2ε/h^2 w/h u_{i1}的系数-ε/h^2可以看到这一组系数里非对角元全是负的对角元明显变大而且每行都严格对角占优因为2ε/h^2 w/h ε/h^2 w/h ε/h^2右边的差正好是ε/h^2这个正数。所以迎风矩阵不仅是Z矩阵是严格对角占优的M矩阵Jacobi和Gauss-Seidel迭代都有收敛保证。那迎风格式和中心差分之间差了什么做一下等价推导在中心差分的基础上额外加入一个人工扩散项κ(u_{i-1} - 2u_i u_{i1})/h^2取κ w h / 2得到的新系数正好就是迎风格式的系数。换句话说迎风格式 中心差分 人工扩散而人工扩散项就是我们要讨论的校正矩阵C。C的作用不是把物理问题改成另一个问题而是把离散矩阵重新拉回M矩阵的安全区。三种系数放一张表格里看一目了然离散格式u_{i-1}系数u_i系数u_{i1}系数中心差分-ε/h^2 - w/(2h)2ε/h^2-ε/h^2 w/(2h)人工扩散项-κ/h^22κ/h^2-κ/h^2迎风格式(中心差分扩散)-ε/h^2 - w/h2ε/h^2 w/h-ε/h^24.3 校正量的下界与参数选择经验既然人工扩散是校正矩阵那扩散系数κ到底取多大合适从恢复Z矩阵的最低要求出发u_{i1}的系数要小于等于0得到条件-ε/h^2 w/(2h) - κ/h^2 ≤ 0解出来就是κ ≥ w h / 2 - ε所以最小校正量是κ_min max(0, w h / 2 - ε)。当Pe ≤ 1时κ_min 0矩阵本来就没问题不需要任何校正当Pe 1时必须至少加上这个量才能把正元压回非正区域。迎风格式取κ w h / 2比κ_min多了ε所以总是安全。但这里有一个实战中的权衡人工扩散的本质是一个数值粘性加得越小越贴近原始方程加得越大解越被抹平尤其是边界层这类剧烈变化的区域扩散大了几乎看不出梯度。我自己的做法是先把κ_min算出来然后从κ_min开始往上扫描观察两个指标一个是迭代步数或残差下降率另一个是解的最大值变化。通常在κ_min的1.2到1.5倍这个区间就能找到一个很稳的工作点没必要一上来就把迎风系数拉满。还有一个工程上容易混淆的点当网格加密时h减小Pe会随h线性下降原来发散的对流占优矩阵可能自动回到收敛区。所以有时候加密网格之后迭代突然正常了并不完全是因为精度变好而是矩阵结构自己恢复了良性。这时候再回头看校正策略它其实给了你一个预判手段不用等到加密网格试算光看Pe数就知道这个离散矩阵会不会在迭代法上制造麻烦。5. 做完这套策略之后我踩过的五个坑5.1 别把Gershgorin圆盘当精确判据这套策略的第一步需要判断矩阵到底坏没坏很多人图省事只做一次Gershgorin圆盘检查如果发现不严格对角占优立刻就开始校正。问题是严格对角占优是充分的但不是必要的真实世界里有大量矩阵不满足占优条件谱半径却依然小于1迭代照样收敛。比如对称正定矩阵配合CG类方法完全不依赖占优结构。所以我现在的流程是先用幂法跑大概500次矩阵向量乘粗估一下Jacobi迭代矩阵的谱半径再决定要不要启动校正流程。幂法只要一个矩阵向量乘函数很多求解器里本来就有现成的成本很低。这一步能拦住大量不必要的手术毕竟矩阵能不动刀就尽量不动刀。5.2 校正之后忘了补偿解得的是错误的方程这个坑我最少见过三次包括我自己早期也犯过。做完A C之后有些人会把修正后的矩阵直接扔给求解器解A_tilde x b结果解出来的东西自然不对。因为你要的是Ax b的解A_tilde x b的解是另一回事。正确的做法是第3节里的不动点格式求解每一步时右端要带上C x^k项。很多迭代求解器接口只支持固定右端向量这时候就得多写一层包装把迭代格式拆解成矩阵向量乘和线性系统求解两步。实现上多花十分钟但忘记补偿的后果是整个结果不可信而且很难排查因为残余可能会随着迭代缓慢下降给人的直觉是大概快收敛了实际解早偏了。5.3 行缩放治不了收敛病矩阵行缩放也就是对每一行乘以一个正数在预处理里很常见。不少人遇到收敛问题第一反应是我这矩阵量纲不对缩放一下就好了。但数学上有个很好玩的事实行缩放不改变Jacobi迭代矩阵。设缩放后的矩阵是A R A其中R是对角正矩阵。那么A的对角部分D R DJacobi迭代矩阵T I - D^{-1} A I - (R D)^{-1} R A I - D^{-1} A和原来一模一样。所以行缩放对Jacobi迭代的谱半径没有任何影响指望它修复收敛性纯粹是白费力气。Gauss-Seidel也类似行缩放不改变迭代矩阵的本质结构。这个结论经常让人惊讶但它恰恰说明收敛性是矩阵的内在结构决定的不是量纲问题。列缩放则会改变迭代矩阵但那样做需要配合等价变换通常更适合放在预条件设计里而不是当作收敛修补工具。5.4 校正量越大越好的错觉校正矩阵的目的是把矩阵拉回收敛区但过量的校正会引入另一个问题它会改变原方程解的物理特性。在对流扩散方程里人工扩散系数加得过大边界层整个被抹平在一般代数矩阵里过大的对角增补会让系统变得过于刚性条件数劣化迭代收敛速度反而可能下降。我自己调试时习惯写一个很小的扫描脚本把校正量从κ_min开始递增每隔几步记录一次残差下降率。结果几乎总是存在一个最优工作点过了这个点残差下降率开始变差或者解开始漂移。所以校正策略的格言应该是够用就好而不是越多越好。5.5 浮点误差与小对角元留出安全带的余量最后一条是关于数值实现的细节。在判断对角占优条件时不要用严格的小于号而是留一个安全带余量。比如每行非对角绝对值之和小于0.95倍的对角元绝对值再判定为安全。这样做主要是因为浮点计算和矩阵组装过程中的舍入误差可能让本应刚好成立的占优条件变得不稳定尤其是近奇异矩阵一行小小的误差扰动就可能让收敛性质反转。另外在步进迭代过程中残差评估本身也有精度上限。如果你发现残差下降到10^{-8}左右不再动不要急着怀疑校正策略先检查一下是不是浮点累加已经到极限精度了。改用Kahan求和或者更高精度来评估残差很多时候只是虚惊一场。把这几条经验放在一起其实是在强调同一件事校正矩阵策略不是机械套公式它是一套需要理解、诊断、验证的完整流程。每一步都要知道自己在修什么、为什么修、修完之后怎么验证才能真正把不收敛问题变成收敛问题。这也是我从最初只会调松弛因子到现在遇到迭代问题先看矩阵结构最大的一个转变。