3. 线性方程组直接解法与选型

工程计算中大量问题最终都归结为求解线性方程组 Ax = b：结构有限元中的刚度方程、电路网络的节点方程、数据拟合中的正规方程，都是这一形式。直接解法的共同思路是把矩阵 A 逐步化为易于回代的形式，理论上在无舍入误差时给出精确解。本节围绕“怎么算”和“怎么选”两条线索展开：先看清高斯消去法与 LU 分解的等价关系，再理解 Cholesky 分解对对称正定矩阵的专门优化，然后讨论主元选择如何决定数值稳定性，最后把稀疏性与矩阵规模纳入选型判断。

高斯消去法与 LU 分解：同一过程的两种视角

高斯消去法分两步。第一步是消元：对 k = 1, 2, …, n-1，用第 k 行消去下方各行第 k 列的元素，乘数为 m_ik = a_ik / a_kk；第二步是回代，从最后一个未知量开始向上求解上三角方程组。以三阶系统为例，若消元后得到上三角矩阵 U，则回代公式为 x_i = (b_i - Σ_{j>i} u_ij x_j) / u_ii。

关键观察是：消元过程中记录下来的乘数 m_ik，恰好构成一个单位下三角矩阵 L，使得 A = LU。这就是 LU 分解。因此高斯消去法可以理解为“边分解边求解”，而 LU 分解把分解与求解拆开：分解一次得到 L、U，之后对每个右端项只需一次前代 Ly = b 和一次回代 Ux = y。当右端项有多个时（例如同一刚度矩阵对应多种荷载工况），这一拆分的价值非常明显：分解的 O(n³) 代价只付一次，每个右端项只付 O(n²)。

需要强调，朴素高斯消去法在遇到 a_kk = 0 时会直接失败，即使矩阵非奇异；更隐蔽的问题是 a_kk 虽不为零但绝对值很小，会导致乘数急剧放大，舍入误差被严重放大。这正是主元选择要解决的问题。

Cholesky 分解：对称正定矩阵的专用通道

当 A 对称正定时，存在唯一的分解 A = LLᵀ，其中 L 是下三角矩阵且对角元为正。计算公式为：

for j = 1..n:
    L[j][j] = sqrt(A[j][j] - Σ_{k<j} L[j][k]^2)
    for i = j+1..n:
        L[i][j] = (A[i][j] - Σ_{k<j} L[i][k]*L[j][k]) / L[j][j]

与 LU 分解相比，Cholesky 分解只需存储和计算一半的矩阵元素，运算量约为 LU 的一半，且无需选主元——对称正定性本身保证了数值稳定性（只要矩阵条件数不过分恶劣）。判断矩阵是否对称正定，实践中常用尝试分解法：若在计算 L[j][j] 时根号内出现非正数，说明矩阵不是正定的，应改用带主元的 LU 分解。有限元中的刚度矩阵在施加足够约束后通常是对称正定的，因此 Cholesky 是这类问题的首选。

主元选择策略：稳定性的守门人

主元选择的核心是在第 k 步消元前，从第 k 列（或第 k 行）的候选元素中挑一个绝对值较大的作为主元，通过行交换（或行列交换）把它换到对角位置。常见策略有三种：

• 部分主元（行交换）：在第 k 列的第 k 行及以下找绝对值最大者，交换行。实现简单，是通用 LU 分解的默认选择，对应分解形式 PA = LU。
• 完全主元（行列交换）：在整个右下子矩阵中找绝对值最大者，稳定性最好，但每次都要搜索子矩阵并交换行列，代价高，实际工程中较少使用。
• 无主元：仅当矩阵具有特殊结构（对称正定、对角占优）时才安全，例如 Cholesky 分解。

一个直观的示例：若矩阵首行为 [1e-18, 1]、次行为 [1, 1]，不做主元交换时乘数高达 1e18，在双精度下有效数字几乎全部丢失；交换两行后乘数约为 1e-18，结果恢复正常。这说明主元选择不是可选项，而是通用直接解法的必要组成部分。

稀疏矩阵直接解法简介

工程中的大规模矩阵往往绝大多数元素为零，例如一维问题离散后每行只有两三个非零元。若按稠密矩阵存储和分解，O(n³) 的运算量和 O(n²) 的存储量会迅速变得不可接受。稀疏直接解法的思路是：只存储非零元，并在消元过程中尽量保持稀疏性。

但消元会引入“填充”（fill-in）：原本为零的位置在消元后变为非零。填充量强烈依赖消元顺序。因此稀疏直接解法的核心工序是重排序，常用方法包括最小度排序（Minimum Degree）、嵌套剖分（Nested Dissection）等，目的是减少填充、降低运算量。工程上通常直接调用成熟的稀疏求解库（如 SuiteSparse 中的 UMFPACK、CHOLMOD，或 Eigen 的 SparseLU、SimplicialLLT），而不是手写分解。选型时要关注：矩阵是否对称、是否正定、重排序策略是否与问题结构匹配、内存是否足够容纳填充。

不同规模与结构下的算法选型

选型不是背结论，而是根据矩阵的四个属性做判断：规模、稀疏性、对称性、正定性。下图把这一判断过程组织为一条决策路径，并标注了各分支的典型代价量级。

结合上图，可以把选型经验归纳为几条可操作的判断。第一，稠密且规模不大（例如 n 在几百到一两千以内）时，带部分主元的 LU 分解是稳妥的通用选择，代码成熟、行为可预期。第二，对称正定时优先 Cholesky，既省一半运算又免去选主元逻辑，但务必先验证正定性，否则根号内出现负数会直接中断。第三，右端项多时一定要把分解与求解分离，避免对每个 b 重复做 O(n³) 的分解。第四，大规模稀疏问题不要自己写稠密分解，应使用稀疏求解库并重视重排序；此时内存往往比运算时间更早成为瓶颈。

最后给出一个可复用的选型检查清单，供实际编码前逐项确认：

需要提醒的是，直接解法给出的是“在浮点运算下尽可能精确”的解，其误差仍受矩阵条件数制约。若矩阵条件数极大，即使算法选型正确，结果也可能不可信，此时应结合残差 ‖Ax - b‖ 与解的合理性做交叉验证，而不是仅凭算法名称判断结果可靠。