3. 高斯消元与LU分解

高斯消元是求解线性方程组最基础也最核心的直接法。它的基本思想并不复杂：通过初等行变换把系数矩阵逐步化为上三角矩阵，再用回代求出未知量。但真正把它写成可靠程序时，问题往往不在“会不会消元”，而在“主元怎么选、分解是否可行、舍入误差会不会被放大”。本节围绕消元过程、主元选择策略以及LU分解的构造与实现展开，并重点分析分解失败与数值不稳定的典型场景。

一、消元过程：从方程组到上三角

设方程组为 Ax = b，其中 A 为 n×n 矩阵。高斯消元在第 k 步以第 k 行第 k 列元素 a_kk 为主元，对第 i 行（i > k）计算乘数 m_ik = a_ik / a_kk，再执行 行i ← 行i − m_ik × 行k，使第 k 列下方元素归零。全部步骤完成后得到上三角矩阵 U，右端项同步变换为 c，最后回代求解。

以一个 3 阶示例说明。取

A = [ 2   1  -1 ]      b = [  8 ]
    [ -3  -1   2 ]          [ -11 ]
    [ -2   1   2 ]          [ -3 ]

第 1 步主元 a_11 = 2，乘数 m_21 = -1.5、m_31 = -1。消元后第 2、3 行变为 [0, 0.5, 0.5 | 1] 与 [0, 2, 1 | 5]。第 2 步主元 a_22 = 0.5，乘数 m_32 = 4，消元后第 3 行变为 [0, 0, -1 | 1]。回代得 x_3 = -1、x_2 = 3、x_1 = 2。这一过程清晰展示了消元与回代的分工，也提示了一个关键点：每一步都依赖主元非零且数值不能过小。

二、主元选择策略：为什么不能“照单全收”

如果某一步主元恰好为 0，消元会因除零而直接失败。例如矩阵首元素为 0 时，按自然顺序消元无法进行。更隐蔽的问题是主元虽非零但绝对值极小，此时乘数会非常大，把原本很小的舍入误差急剧放大，导致结果严重失真。这就是数值不稳定的典型来源。

常用策略有三种。其一是列主元（部分主元）：在第 k 步从第 k 列的第 k 行到第 n 行中选取绝对值最大的元素作为主元，并交换相应两行。它实现简单、额外开销小，是工程实现中的默认选择。其二是全主元：在整个右下子矩阵中选绝对值最大元素，同时交换行和列，稳定性更好但列交换会打乱未知量顺序，需要记录置换，开销也更大。其三是按比例选主元：用该行元素的尺度对候选主元做归一化比较，适合行尺度差异悬殊的问题。

需要强调，列主元并非万能。存在一类矩阵，即使采用列主元，增长因子仍可能很大，例如经典的 Wilkinson 增长矩阵。因此“用了主元就一定稳定”是常见误解，实际还需结合矩阵性质与残差检验判断。

三、LU分解的构造与实现

高斯消元的每一步消元本质上是用初等下三角矩阵左乘 A。把这些乘数保留下来，就得到 A = LU：L 为单位下三角矩阵，其第 i 行第 k 列恰为乘数 m_ik；U 为消元所得上三角矩阵。带行交换时写作 PA = LU，P 为置换矩阵。

实现上通常采用原地存储：把 L 的乘数直接写入 A 的严格下三角部分，U 占用上三角及对角线，从而节省内存。核心循环可概括为：对 k = 1..n-1，先选主元并交换行（同时交换 P 与已存的乘数），再对 i > k 计算并存储乘数、更新子矩阵。分解完成后，求解 Ax = b 转化为三步：解 Ly = Pb（前代）、解 Ux = y（回代）。这样同一矩阵面对多个右端项时只需分解一次，是LU分解相对反复消元的主要优势。

四、分解失败与数值不稳定的典型场景

场景一：主元为零导致分解中断。若矩阵某阶顺序主子式为零，自然顺序消元会在对应步遇到零主元。例如 A = [[0,1],[1,0]]，第一步主元为 0，直接失败；但交换两行后即可正常分解。这说明失败往往不是矩阵本身不可解，而是选主元策略缺失。

场景二：主元极小引发误差放大。考虑 A = [[1e-8, 1],[1, 1]]。若不做行交换，第一步乘数为 1e8，消元后第二行元素量级骤变，在有限精度下有效数字大量丢失，解可能完全错误。采用列主元交换两行后，主元变为 1，乘数约 1e-8，结果稳定。这是“主元绝对值过小必须换行”的最直观例证。

场景三：增长因子过大。即使每步都选列主元，消元过程中元素量级仍可能相对原始矩阵急剧增长。增长因子定义为消元过程中出现的最大元素与原始最大元素之比。增长因子很大时，舍入误差会被显著放大。Wilkinson 矩阵就是构造出来使列主元增长因子达到 2^(n-1) 的经典反例，提醒我们不能把稳定性完全寄托于列主元。

场景四：矩阵接近奇异。当矩阵条件数很大时，即便分解顺利完成，解对右端项扰动也极其敏感。此时应结合条件数估计、残差 ‖Ax − b‖ 以及解的合理性综合判断，而不是只看分解是否“跑通”。

五、实现要点小结

工程实现中，建议默认采用列主元并在分解时同步记录置换；乘数原地写入下三角以节省内存；对主元绝对值设置阈值告警，而不是静默继续；对同一系数矩阵的多次求解复用LU分解结果。理解“消元—主元—分解—求解”这条链路，以及每个环节可能失效的条件，才能真正把高斯消元与LU分解用稳、用对。