做迭代求解的人，大概都经历过这种场面：残差曲线好不容易下来一段，突然卡在某个平台期开始震荡，左调右调松弛因子都没用。后来我养成一个习惯，遇到不收敛先不碰参数，直接判断系数矩阵的结构。见过太多次问题根源不是算法，而是矩阵里藏着几个正的非对角元——它们让整个迭代过程失去了收敛保证。这时候真正该做的，是从收敛条件出发，把矩阵修成M矩阵，也就是俗称的校正矩阵策略。
这篇文章是写给两类人的：一类是和当年的我一样，被不收敛问题折磨得睡不着觉的求解器开发者和数值计算工程师；另一类是学过数值代数但一直没想明白M矩阵到底有什么用的人。我尽量不用教科书式的写法，而是把我反复用过的一套逻辑讲清楚：先看收敛条件的本质是什么，再把它翻译成矩阵元素层面的不等式，最后根据这个不等式去构造校正矩阵。中间会用一个3x3的反例手推一遍完整流程，再用对流占优方程的离散矩阵做一次工程实战，最后聊聊我踩过的几个坑。
1. 先搞懂：迭代收敛到底卡在矩阵的哪个部位
1.1 谱半径条件：一句话版本的收敛判据
几乎所有线性方程组的迭代解法，最后都能写成同一个形式：
x^{k+1} = T x^k + c
这里的T叫迭代矩阵。比如Jacobi迭代里，假设把系数矩阵A拆成对角部分D和非对角部分L+U，那T就是I - D^{-1}A。Gauss-Seidel和SOR本质上也是换了一个T，只是换的方式不同，它们背后的判据完全一样。
误差怎么变？设每次迭代的误差是e^k = x^k - x^*，代入迭代式很容易得到：
e^{k+1} = 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迭代，因为对角元全部是1，D是单位矩阵，迭代矩阵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} = 2：
a_{22} += 1 a_{33} += 1 a_{23} = 0
矩阵变成：
[[2, 0, 0], [0, 3, 0], [2, 0, 2]]
第三步，处理(3,1)位置上的正元a_{31} = 2：
a_{33} += 1 a_{11} += 1 a_{31} = 0
最终得到：
A_tilde = [[3, 0, 0], [0, 3, 0], [0, 0, 3]]
一个原始发散矩阵，经过三个等价的符号吸收步骤，变成了3倍单位矩阵。它当然严格对角占优，也当然是M矩阵，收敛性没有任何悬念。
那校正矩阵C是多少？用A_tilde减去原始A：
C = 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^{k+1} = (A + C)^{-1} (C x^k + b)
在这个3x3例子里，A + C = 3I，所以迭代化成极其简单的标量形式：
x^{k+1} = (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("修正后的矩阵 A+C:") 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(f"k={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_{i+1}的系数：-ε/h^2 + w/(2h)
前两项都是负的，没有问题。问题出在第三项：当w/(2h)超过ε/h^2时，u_{i+1}的系数就会变成正数。定义局部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_{i+1}的系数：-ε/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_{i+1})/h^2，取κ = w h / 2，得到的新系数正好就是迎风格式的系数。换句话说，迎风格式 = 中心差分 + 人工扩散，而人工扩散项就是我们要讨论的校正矩阵C。C的作用不是把物理问题改成另一个问题，而是把离散矩阵重新拉回M矩阵的安全区。
三种系数放一张表格里看，一目了然：
| 离散格式 | u_{i-1}系数 | u_i系数 | u_{i+1}系数 |
|---|---|---|---|
| 中心差分 | -ε/h^2 - w/(2h) | 2ε/h^2 | -ε/h^2 + w/(2h) |
| 人工扩散项 | -κ/h^2 | +2κ/h^2 | -κ/h^2 |
| 迎风格式(中心差分+扩散) | -ε/h^2 - w/h | 2ε/h^2 + w/h | -ε/h^2 |
4.3 校正量的下界与参数选择经验
既然人工扩散是校正矩阵，那扩散系数κ到底取多大合适？从"恢复Z矩阵"的最低要求出发，u_{i+1}的系数要小于等于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 D，Jacobi迭代矩阵：
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求和或者更高精度来评估残差，很多时候只是虚惊一场。
把这几条经验放在一起，其实是在强调同一件事：校正矩阵策略不是机械套公式，它是一套需要理解、诊断、验证的完整流程。每一步都要知道自己在修什么、为什么修、修完之后怎么验证，才能真正把不收敛问题变成收敛问题。这也是我从最初只会调松弛因子，到现在遇到迭代问题先看矩阵结构，最大的一个转变。
