做了这么多年数值计算我越来越觉得线性方程组求解这一关是绕不过去的。它看起来只是个小小的 Axb但你往深了看最小二乘拟合、有限元分析、电路仿真、流体计算、最优化内点法哪个底层不是一遍又一遍在解方程组。可以说解好线性方程组基本上就等于拿到了科学计算的半张入场券。这次想聊聊两类最典型的解法三角分解法和高斯迭代法。前者属于直接法咔咔一顿分解几步回代就能拿到精确解后者属于迭代法先给个初始猜测再一圈圈转着逼近真值。很多初学者会纠结哪个更厉害其实这俩根本不是一个赛道就跟搬家找搬家公司、短途运货找同城快递一样各有各的适用场景用对了才顺手。这篇文章适合正在学《计算方法》《数值分析》这门课的学生也适合工作中偶尔要解方程但不想每次都调现成库的工程师。我会把两种方法背后的原理、手算步骤、编程时的注意点以及我踩过的那些坑都掰开揉碎讲一遍文末还会用同一个例题把两种方法各跑一遍让你直观看到它们在速度、精度和使用条件上的差异。1. 先把问题摆到桌面上我们要解的到底是什么1.1 线性方程组从哪里冒出来的在讨论算法之前先明确一下对象。我们面对的是这样一类问题给定一个 n×n 的系数矩阵 A 和右端向量 b求满足 Axb 的未知向量 x。这里的 A 不一定是玩具矩阵。我在实际项目里见过从几千阶到几十万阶的矩阵来源五花八门有限元/有限差分把连续问题离散化刚度矩阵就是 A节点位移就是 x电路节点分析对每个节点列基尔霍夫方程导纳矩阵就是 A最小二乘拟合求法方程 (AᵀA)x Aᵀb最优化内点法每次迭代都要解一个牛顿方向的线性系统。你可以把 A 想象成一张地图b 是目的地坐标x 就是我们要规划的路线。地图的结构决定了该用什么交通工具也就是用什么算法。1.2 两条路线两种哲学直接法理论上经过有限步算术运算就能得到精确解误差主要来自舍入。最经典的就是高斯消元法以及它的变体LU 分解/三角分解。迭代法从一个近似解出发构造一个不断修正的序列让它在极限上收敛到精确解。最基础的是雅可比迭代Jacobi和高斯-赛德尔迭代Gauss-Seidel也就是标题里说的“高斯迭代法”的常用实现。选直接法还是迭代法不能拍脑袋。我的经验是中等规模、矩阵较稠密用直接法大规模稀疏矩阵尤其是来自偏微分方程离散化的那类问题几乎只能用迭代法——因为直接法在分解过程中会产生“填充”把稀疏矩阵弄成一个稠密的因子内存直接爆炸。迭代法只要矩阵能拿来算乘法就行存储压力小得多。2. 三角分解法把复杂问题拆成两个简单问题2.1 为什么非要把 A 拆成 LU三角分解法俗称LU 分解核心思想非常朴素把一个复杂的矩阵 A 拆成下三角矩阵 L和上三角矩阵 U的乘积即 A LU。拆完之后解 Axb 就变成解两步先解 Lyb这叫前代再解 Uxy这叫回代。三角形方程组的求解是极度机械的一个方程只含一个未知数顺着来就行。三角分解的精髓在于把一次复杂的消元过程变成一个可以重复利用的“拆解结果”。2.2 拆解原理消元法的另一种解读你可能在高斯消元里见过这种操作用第一行消去下面各行的第一列元素再用第二行消去第三行及以下的第二列元素……消元过程中每次做的行变换其实等价于在矩阵左侧乘一个初等矩阵。把所有初等操作矩阵乘在一起得到的就是 L而消元结束后剩下的阶梯矩阵就是 U。整个高斯消元的流程天然就是一个 LU 分解。反过来理解这句话LU 分解就是高斯消元的“存档”版本把过程记录下来之后不管右端项 b 怎么变都不用重新消元。2.3 手算三步走Doolittle 分解最常用的 LU 分解叫Doolittle 分解它规定 L 的对角线元素全是 1。手算的流程我推荐按这个顺序来先写 U 的第一行 A 的第一行再算 L 的第一列公式是 lᵢ₁ aᵢ₁ / u₁₁i 2, 3, …, n按“交替算一行 U、一列 L”的顺序推进直到分解完成。具体到第 k 步U 的第 k 行第 j 列计算公式是u(k,j) a(k,j) - Σ(t1..k-1) l(k,t)·u(t,j)L 的第 i 行第 k 列计算公式则是l(i,k) [a(i,k) - Σ(t1..k-1) l(i,t)·u(t,k)] / u(k,k)光看公式容易晕我拿一个 3 阶例子走一遍。设A [4, 1, -1; 1, 4, 2; -1, 2, 4]对应的 b 我取 [3, 15, 15]这个方程组精确解是 x [1, 2, 3]方便验证。第一步U 的第一行就是 [4, 1, -1]。第二步算 L 第一列l21 1 / 4 0.25l31 -1 / 4 -0.25第三步算 U 第二行u22 4 - l21·u12 4 - 0.25×1 3.75u23 2 - l21·u13 2 - 0.25×(-1) 2.25第四步算 L 第二列l32 (2 - l31·u12) / u22 (2 0.25) / 3.75 0.6第五步算 U 第三行u33 4 - l31·u13 - l32·u23 4 - 0.25 - 0.6×2.25 4 - 0.25 - 1.35 2.4于是分解结果为L [[1, 0, 0], [0.25, 1, 0], [-0.25, 0.6, 1]] U [[4, 1, -1], [0, 3.75, 2.25], [0, 0, 2.4]]你随手验一下 L×U会发现等于 A分毫不差。接着解 Lyby1 3y2 15 - 0.25×3 14.25y3 15 - (-0.25×3) - 0.6×14.25 7.2再解 Uxyx3 7.2 / 2.4 3x2 (14.25 - 2.25×3) / 3.75 2x1 (3 - 1×2 1×3) / 4 1完美得到 [1, 2, 3]。2.4 程序实现要点手算只是为了理解实际编程一般按“原地覆盖”的思路实现L 和 U 直接存在 A 的位置里对角线以下存 L对角线及其以上存 U三重循环外层对 k 从 1 到 n-1中间行 i 从 k1 到 n内层算 L 和 U每算一个元素就用 A 原位置覆盖省内存也方便后面做多次右端项求解。注意这个基础版本有一个很大的问题——它不选主元。一旦遇到对角线元素接近 0u(k,k) 会非常小除出来的 L 元素会爆炸后面误差会成倍放大。所以实际工程里几乎都会用列主元 LU 分解每步先在第 k 列及以下找绝对值最大的元素把那一行调到第 k 行再做分解。选主元这个操作看着不起眼但它是直接法数值稳定性的保命符。我做实验时经常看到同样的方程组不选主元跑出来解完全对不上选了主元立刻正常。原因不是算法错了而是浮点运算下小主元会放大舍入误差。2.5 LU 分解的隐藏优势一次分解多次使用工程里很多场景是 A 不变、b 变。比如你需要对 100 个不同的载荷工况计算结构响应左端刚度矩阵始终相同。这种情况 LU 分解的优势就非常明显第一次做分解的代价是 O(n³)之后每个 b 只需要做一次前代和一次回代代价是 O(n²)。相比之下如果你每次都重新做高斯消元那就是 100 次 O(n³)。所以很多软件内部在解线性方程组时用户真正调用的其实是“先分解后求解”两步接口。3. 高斯迭代法让答案自己慢慢“转”出来3.1 雅可比迭代和高斯-赛德尔迭代迭代法的思路完全不同。它不给有限步精确解而是构造一个序列 x⁽ᵏ⁾使它收敛到精确解。以 3 阶方程组为例雅可比迭代的格式是x1⁽ᵏ⁺¹⁾ (b1 - a12·x2⁽ᵏ⁾ - a13·x3⁽ᵏ⁾) / a11 x2⁽ᵏ⁺¹⁾ (b2 - a21·x1⁽ᵏ⁾ - a23·x3⁽ᵏ⁾) / a22 x3⁽ᵏ⁺¹⁾ (b3 - a31·x1⁽ᵏ⁾ - a32·x2⁽ᵏ⁾) / a33注意等号右边的全部来自上一轮的结果这一轮的每个新值只用旧值算这叫“同时更新”。高斯-赛德尔迭代则更“勤劳”一点它算完 x1 的新值立刻就用新值去算 x2x2 的新值又立刻传给 x3。也就是x1⁽ᵏ⁺¹⁾ (b1 - a12·x2⁽ᵏ⁾ - a13·x3⁽ᵏ⁾) / a11 x2⁽ᵏ⁺¹⁾ (b2 - a21·x1⁽ᵏ⁺¹⁾ - a23·x3⁽ᵏ⁾) / a22 x3⁽ᵏ⁺¹⁾ (b3 - a31·x1⁽ᵏ⁺¹⁾ - a32·x2⁽ᵏ⁺¹⁾) / a33Gauss-Seidel 的直观优势是信息利用更充分通常收敛更快但实现时要注意更新顺序且它不适合直接并行化因为每一步都依赖最新的分量。3.2 什么时候迭代法会失灵迭代法的“命门”是收敛性。你有没有遇到过这样的情形明明方程组很简单迭代起来却发散了值越算越大最后直接溢出。我当年第一次跑 Jacobi 就遇到过原因是矩阵不满足严格对角占优。严格对角占优的定义是每一行对角线元素的绝对值必须大于同行所有非对角线元素绝对值之和即 |a(i,i)| Σ(j≠i) |a(i,j)|。这只是一个充分条件不是必要条件。很多时候矩阵不满足这个条件也能收敛但初学者先记住这一条最安全。判定收敛性的黄金标准是看迭代矩阵的谱半径ρ(B)也就是 B 的模最大的特征值的模。谱半径小于 1迭代必收敛大于等于 1迭代发散或震荡。提示直接用特征值判断收敛性需要额外付出 O(n³) 的求模运算这在大规模问题里本身就不便宜。所以工程上更常用的做法是先判断对角占优再用迭代过程的残差变化来观测收敛趋势。3.3 怎么判断迭代到什么时候该停迭代是无限逼近的过程总得有个停下来的标准。常用的有三类绝对误差准则||x⁽ᵏ⁺¹⁾ - x⁽ᵏ⁾|| ε比如小于 1e-6 就停相对误差准则||x⁽ᵏ⁺¹⁾ - x⁽ᵏ⁾|| / ||x⁽ᵏ⁾|| ε更适合量级差异大的未知量残差准则||Ax⁽ᵏ⁾ - b|| ε直接评估“没满足方程的程度”。实际代码里我一般同时用相对误差和残差两个中有一个达标就停。另外要设置最大迭代次数比如 1000 次防止某些病态矩阵在错误条件下无限循环。3.4 伪代码参考以 Gauss-Seidel 为例伪代码相当简洁输入: A, b, 初始猜测 x, 容差 eps, 最大迭代次数 maxIt for k 1 to maxIt: for i 1 to n: 新值 b[i] for j 1 to n: if j ! i: 新值 新值 - A[i][j] * x[j] x[i] 新值 / A[i][i] end for 计算残差 r ||A*x - b|| if r eps: break end for有一点很容易踩坑在 Gauss-Seidel 中ji 时要用本轮已经更新的 x[j]ji 时才用上一轮的旧值。有些同学实现的时候没注意这个细节写成了全程用旧值那就变成 Jacobi 了。4. 实操对比同一个方程组两种方法谁更快4.1 例题设定还是用第二节那个矩阵A [4, 1, -1; 1, 4, 2; -1, 2, 4] b [3, 15, 15]这个矩阵对角线是 4每行非对角和的绝对值分别为 2、3、3严格对角占优成立迭代法必然收敛。4.2 三角分解法实操记录直接用前面算好的 L 和 U分解过程 5 步搞定回代 3 步出结果总共 8 步拿到精确解。全程没有任何近似得到的就是 [1, 2, 3]。如果用代码跑核心逻辑大致是import numpy as np A np.array([[4.0, 1.0, -1.0], [1.0, 4.0, 2.0], [-1.0, 2.0, 4.0]]) b np.array([3.0, 15.0, 15.0]) # Doolittle 分解 n A.shape[0] L np.eye(n) U A.copy() for k in range(n): for i in range(k1, n): factor U[i, k] / U[k, k] L[i, k] factor U[i, k:] - factor * U[k, k:] # 前代解 Ly b y np.zeros(n) for i in range(n): y[i] (b[i] - L[i, :i] y[:i]) / L[i, i] # 回代解 Ux y x_lu np.zeros(n) for i in range(n-1, -1, -1): x_lu[i] (y[i] - U[i, i1:] x_lu[i1:]) / U[i, i]这段代码虽然简单但分解、前代、回代分得很清楚适合用来验证手算结果。实际开发中我一般直接调用成熟库里的带选主元接口至少保证数值稳定性。4.3 高斯迭代法实操记录用 Gauss-Seidel 迭代从 x⁰ [0, 0, 0] 出发第 1 轮x1 (3 - 0 0) / 4 0.75x2 (15 - 1×0.75 - 2×0) / 4 3.5625x3 (15 1×0.75 - 2×3.5625) / 4 2.15625第 2 轮x1 (3 - 3.5625 2.15625) / 4 0.3984375x2 (15 - 0.3984375 - 2×2.15625) / 4 2.572265625x3 (15 0.3984375 - 2×2.572265625) / 4 2.5634765625此时已经能看到解在朝 [1, 2, 3] 靠拢。大概迭代 8 轮左右各分量误差就小于 1e-6。如果用 Jacobi 迭代从同样的初始值出发速度略慢大概要多转几轮但也能到同样的精度。4.4 两个方法各花了多少代价我把两者算账对比一下见下表对比项三角分解法高斯迭代法计算代价分解 O(n³/3)每次回代 O(n²)每轮 O(n²)总代价 O(m·n²)m 为迭代轮数内存占用需要存 L 和 U或原地覆盖 A只需存 A 本身稀疏时极省存储适解规模n 在几千以内比较舒服几十万阶稀疏矩阵的主力方法精度理论上有限步精确解只能逼近可控精度对稀疏结构分解过程会产生填充破坏稀疏性天然保留稀疏性依赖稀疏矩阵乘法并行友好度分解过程有顺序依赖Jacobi 天然可并行Gauss-Seidel 较难看这个表就明白了如果你只解一次中等规模的稠密方程LU 分解又快又准如果矩阵来自三维有限元几十万阶还想在普通电脑上跑迭代法几乎是唯一选择。5. 避坑清单那些年我在线性方程组上踩过的雷5.1 矩阵奇异LU 解不出来迭代发“飘”有一次调试一个模型系数矩阵对角线出了个很大的数LU 分解时某一步主元变成 0程序直接报错。当时第一反应是 “哎呀分解公式写错了”查了半天才发现是矩阵本身奇异问题出在模型的约束条件不足。判断矩阵是否奇异可以靠 LU 分解时是否出现零主元来发现。迭代法在奇异问题上表现更隐蔽可能看起来收敛但解随初值漂移或者残差永远降不到零。所以拿到一个问题我习惯先求一下条件数。条件数巨大比如 1e12时即使矩阵非奇异解也会对右端项 b 的微小扰动极其敏感。5.2 不满足对角占优时怎么办遇到不占优的矩阵迭代法大概率不收敛。我的建议按优先级试重排行方程数值上对角占优可能只是行序不好交换行之后可能就占优了行缩放把每一行乘以适当常数改善对角线的相对大小换用直接法小规模问题时直接用 LU 分解别硬刚迭代法改用 SOR逐次超松弛给 Gauss-Seidel 加一个松弛因子 ω调参得当可以显著加速甚至让本不收敛的问题收敛。SOR 里的 ω 一般是 (0, 2) 区间最优取值比较微妙。我的经验是从 ω1 开始此时就是 Gauss-Seidel逐步微调看哪一组收敛最快。5.3 迭代法发散程序崩溃发散常常表现为解的分量指数级增长。遇到这种事最该做的不是调循环次数而是先检查三个地方对角线上有没有零元素如果 a(i,i)0公式直接除以 0程序必炸这种情况需要先交换行或列是否满足对角占优不满足收敛条件发散是“正常”的初始值是否离谱有些问题对初值敏感虽然理论上任意初值都收敛但一个差的初值可能让前几百轮看起来很“炸”容易误判。我一个比较笨但有效的习惯是在迭代循环里加上对 x 的模的监控发现增长超过一定阈值立刻中断并打印诊断信息。这样即便发散也能快速定位是在哪一轮、哪个分量开始崩的。5.4 收敛慢得像蜗牛爬收敛性满足但就是慢也是常见问题。如果你的矩阵是非对称、非正定的Gauss-Seidel 可能比 Jacobi 还慢这在理论上都是可能的并不稀奇。提速思路主要是预处理。实际工程里常用的有对角缩放把 A 的每行按对角元素归一化让对角线元素都变成 1改善谱分布不完全 LU 预处理对稀疏矩阵做一个“阉割版”LU 分解作为迭代前的变换能大幅压低谱半径共轭梯度法CG如果矩阵对称正定别用 Jacobi/Gauss-Seidel直接上 CG收敛速度快一个等级。我在某公司的仿真项目里就遇到过原始 Gauss-Seidel 要迭代上万次才收敛的稀疏矩阵做了简单的不完全 LU 预处理后几十次就收敛了。这个差距在工程上是决定性的。5.5 误差判断两件套无论用哪种方法我都建议最终输出同时附上残差范数和解的相对误差估算而不是只看迭代停止条件。迭代停止条件只能说明“迭代步之间变化小”不能说明“解真的正确”。如果矩阵病态相邻迭代步变化很小但离真实解还很远。我的检查“两件套”是计算 r Ax - b看它是零附近还是偏大换一个不同的初始值重新迭代看会不会得到不同的极限值如果不同基本可以说矩阵病态需要警惕。结束前再说几句体己话这两类方法我用了很长时间最大的一个体会是别崇拜哪一个也别贬低哪一个。LU 分解稳定高效是中小规模问题的默认答案迭代法存储友好、擅长稀疏是大规模工程问题的生存之道。真正的高手不是死记算法步骤而是能在一开始就判断出问题属于哪一类然后选择合适的工具。关于细节我还有一个私心推荐不管题目给的是不是对角占优矩阵你要真想快速验证一个算法的实现有没有问题最好的办法就是拿一个已知精确解的方程组来试错。就像我上面选的例子解是整数 [1, 2, 3]每一步手算都能对得上。等算法在简单例子上跑通了再换真实工程矩阵也心里有底。如果你正准备课程设计或者面试算法岗位建议把 LU 分解的程序背熟把 Jacobi 和 Gauss-Seidel 的公式推导自己走一遍再对着本文第五节把几个坑都踩一遍这个知识点就算真的吃透了。数值计算的功夫说到底是靠动手算、动手写、动手调试堆出来的光看书是看不出感觉的。
