Python实战：用NumPy手写高斯消元法解线性方程组（附完整代码）

从零构建：用NumPy实现工业级高斯消元法，攻克线性方程组实战难题

在数据科学和工程计算的日常工作中，线性方程组就像空气一样无处不在。从物理仿真、金融建模到机器学习算法的优化，我们总在和各种形式的 Ax = b 打交道。很多朋友一遇到这类问题，第一反应就是调用 numpy.linalg.solve 或者 scipy.linalg.solve，这当然没问题。但你是否想过，当库函数返回一个“奇异矩阵”错误时，背后到底发生了什么？或者当矩阵规模巨大、条件数很差时，如何确保解的稳定性和精度？

理解高斯消元法，就是理解这一切的基石。它不仅是线性代数教材里的经典算法，更是众多高级分解方法（如LU、Cholesky分解）的灵魂。今天，我们不满足于理论推导，而是要亲手用NumPy实现一个鲁棒、实用的高斯消元求解器。我会带你穿越从基础实现到处理边界情况的完整路径，让你不仅“知其然”，更能“知其所以然”，在未来的项目中拥有自主调试和优化的能力。

这篇文章面向的是已经熟悉Python和NumPy基础，希望深入算法底层、提升工程实践能力的数据科学家和工程师。我们将采用Jupyter Notebook风格的交互式代码演示，确保每一步都可执行、可验证。

1. 高斯消元法：核心思想与算法骨架

高斯消元法的目标非常直接：将一个普通的系数矩阵 A，通过一系列行变换，转化为上三角矩阵 U。一旦得到上三角矩阵，我们就可以通过“回代”这个简单的过程，从最后一个方程开始，逐个解出所有未知数。

这个过程的本质，是用矩阵的行初等变换，保持方程组的等价性。具体来说，我们只允许三种操作：

1. 交换两行的位置（行交换）。
2. 将某一行乘以一个非零常数（行缩放）。
3. 将一行的倍数加到另一行上（行倍加）。

我们主要利用第三种操作。想象一下，对于一个 n x n 的矩阵，消元过程会进行 n-1 步。在第 k 步，我们的目标是利用第 k 行，将第 k 列中第 k+1 行到第 n 行的元素全部变为零。

算法的核心伪代码可以概括如下：

对于 k 从 0 到 n-2:
    对于 i 从 k+1 到 n-1:
        计算乘子 factor = A[i, k] / A[k, k]
        对于 j 从 k 到 n-1:
            A[i, j] = A[i, j] - factor * A[k, j]
        b[i] = b[i] - factor * b[k]

这个过程结束后，矩阵 A 就变成了上三角矩阵。紧接着是回代：

创建解向量 x，长度为 n
对于 i 从 n-1 到 0 (倒序):
    sum = b[i]
    对于 j 从 i+1 到 n-1:
        sum = sum - A[i, j] * x[j]
    x[i] = sum / A[i, i]

听起来很简单，对吧？但如果我们直接把上面的伪代码翻译成Python，会遇到几个致命的陷阱：

1. 除零错误：当主对角线元素 A[k, k] 为零时，计算 factor 会导致程序崩溃。
2. 数值不稳定：即使 A[k, k] 不为零，但如果它非常接近于零，作为除数会导致 factor 极大，放大舍入误差，最终得到完全错误的解。
3. 效率问题：三重循环的朴素实现，在 n 较大时速度很慢。

我们先给出一个最基础的、不考虑任何问题的“教学版本”实现，感受一下流程。

import numpy as np

def gaussian_elimination_naive(A, b):
    """
    基础版高斯消元法（无主元选择，不稳定，仅用于演示）
    参数:
        A: numpy.ndarray, 形状为 (n, n) 的系数矩阵
        b: numpy.ndarray, 形状为 (n,) 的常数向量
    返回:
        x: numpy.ndarray, 解向量
    """
    # 确保是浮点数类型，避免整数除法问题
    A = A.astype(float)
    b = b.astype(float)
    n = len(A)

    # 消元过程
    for k in range(n-1):
        for i in range(k+1, n):
            factor = A[i, k] / A[k, k]
            A[i, k:] = A[i, k:] - factor * A[k, k:]
            b[i] = b[i] - factor * b[k]

    # 回代过程
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]
    return x

# 测试一个简单的例子
A_test = np.array([[2., 1., -1.],
                   [-3., -1., 2.],
                   [-2., 1., 2.]])
b_test = np.array([8., -11., -3.])

x_naive = gaussian_elimination_naive(A_test.copy(), b_test.copy())
print("基础版求解结果:", x_naive)
print("验证 Ax - b:", np.dot(A_test, x_naive) - b_test)

运行这段代码，对于这个精心挑选的例子，你会得到正确的解 [2., 3., -1.]。但请不要高兴太早，这个版本脆弱得不堪一击。我们马上就会看到它的问题。

2. 核心挑战与实战技巧：主元选择与精度控制

上面“教学版本”失败的核心在于它盲目地将每行的第一个非零元素（即当前列的对角线元素 A[k, k]）作为除数。这引出了数值计算中一个至关重要的概念：主元（Pivot）。主元的选择直接决定了算法的稳定性和成败。

2.1 为什么需要选主元？

考虑下面这个例子：

A_bad = np.array([[1e-15, 1.],
                   [1., 1.]])
b_bad = np.array([1., 2.])

这个方程组的精确解大约是 x1 = 1.0, x2 = 1.0。但如果我们用基础版算法，第一步计算 factor = A[1,0] / A[0,0] = 1.0 / 1e-15 = 1e15。这个巨大的乘子会使得 A[1,1] 在计算 A[1,1] - factor * A[0,1] 时，发生严重的有效数字丢失（Catastrophic Cancellation），因为 1 - 1e15 * 1 在浮点数运算中，1 的信息完全被吞没了。结果会溢出或产生极大的误差。

注意：在浮点数运算中，当两个数量级相差巨大的数相加减时，较小数的精度会完全丢失。这是数值误差的主要来源之一。

解决方案就是选主元：在每一步消元之前，从当前列（第 k 列）的第 k 行到第 n 行中，选取绝对值最大的元素作为主元。然后，通过行交换，将这个最大元素放到 A[k, k] 的位置上。这种方法被称为部分选主元（Partial Pivoting）。

为什么选最大值？因为作为除数的绝对值越大，计算出的乘子 factor 的绝对值就越小（≤1），从而可以最大限度地减少后续计算中舍入误差的放大效应。

2.2 实现带部分选主元的高斯消元

我们需要在消元循环中增加选主元的逻辑。同时，为了后续可能需要的行交换信息（这在LU分解中很有用），我们引入一个记录行交换顺序的数组 pivot。

def gaussian_elimination_partial_pivoting(A, b):
    """
    带部分选主元的高斯消元法
    参数:
        A: numpy.ndarray, 形状为 (n, n) 的系数矩阵
        b: numpy.ndarray, 形状为 (n,) 的常数向量
    返回:
        x: numpy.ndarray, 解向量
        L: numpy.ndarray, 单位下三角矩阵（用于理解，非必需返回）
        U: numpy.ndarray, 上三角矩阵
        pivot: list, 行交换记录
    """
    A = A.astype(float).copy()
    b = b.astype(float).copy()
    n = len(A)
    pivot = list(range(n))  # 初始行顺序

    # 消元过程（含选主元）
    for k in range(n-1):
        # --- 部分选主元 ---
        # 找到第k列中，从第k行开始绝对值最大的行索引
        max_row_index = np.argmax(np.abs(A[k:, k])) + k
        if max_row_index != k:
            # 交换矩阵A的行
            A[[k, max_row_index], :] = A[[max_row_index, k], :]
            # 交换向量b的元素
            b[[k, max_row_index]] = b[[max_row_index], k]
            # 记录行交换
            pivot[k], pivot[max_row_index] = pivot[max_row_index], pivot[k]

        # 检查主元是否为零（选主元后仍可能为零，意味着矩阵奇异）
        if np.abs(A[k, k]) < 1e-15:  # 使用一个很小的容差
            raise ValueError(f"矩阵在消元第{k}步后主元接近零，矩阵可能奇异。")

        # --- 消元 ---
        for i in range(k+1, n):
            factor = A[i, k] / A[k, k]
            # 使用向量化操作提高效率，只计算k列之后的元素
            A[i, k+1:] = A[i, k+1:] - factor * A[k, k+1:]
            # 同时处理b向量
            b[i] = b[i] - factor * b[k]
            # 将乘子存储在已消元的位置，便于后续构建L矩阵（可选）
            A[i, k] = factor

    # 此时A已成为上三角矩阵U（对角线以下存储了L的乘子）
    U = np.triu(A)  # 取上三角部分，包括对角线
    # 构建单位下三角矩阵L（对角线为1，下三角为消元乘子）
    L = np.eye(n) + np.tril(A, -1)

    # 回代求解
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]

    return x, L, U, pivot

# 测试之前的“坏”矩阵
try:
    x_pp, L_pp, U_pp, p_pp = gaussian_elimination_partial_pivoting(A_bad.copy(), b_bad.copy())
    print("部分选主元求解结果:", x_pp)
    print("验证 Ax - b:", np.dot(A_bad, x_pp) - b_bad)
except ValueError as e:
    print(f"求解失败: {e}")

通过选主元，算法会先交换两行，用 1.0 作为主元，从而避免了巨大乘子的产生，得到了稳定的解。

2.3 精度控制与条件数

即使用了选主元，有些问题本质上是“病态”的。这由矩阵的条件数（Condition Number） 决定。条件数 cond(A) 衡量了输出值（解 x）对输入值（A 和 b）微小变化的敏感程度。条件数越大，问题越病态，数值解越不可靠。

我们可以用NumPy快速估算条件数：

cond_number = np.linalg.cond(A_test)
print(f"矩阵A的条件数: {cond_number:.2e}")

对于条件数很大的矩阵（比如大于 1e10），即使使用最稳定的算法，得到的解也可能与真实解相差甚远。在实际工程中，如果遇到条件数过大的情况，我们需要：

1. 重新审视问题：检查物理或数学模型是否合理，量纲是否统一。
2. 使用更高精度：使用 np.float128（如果系统支持）或 decimal 库。
3. 采用正则化方法：例如在机器学习中，通过添加L2正则项（岭回归）来改善矩阵的条件数。

下表对比了不同版本高斯消元法的特性：

特性	基础版 (无主元)	部分选主元版	完全选主元版
稳定性	极差，易因小主元失败	良好，适用于绝大多数情况	最优，但计算开销大
计算复杂度	O(n³)	O(n³)，增加选主元比较	O(n³)，增加行列搜索
实现难度	简单	中等	复杂
适用场景	仅教学演示	通用科学计算、工程应用	极端病态问题、符号计算

提示：完全选主元会在当前右下子矩阵中搜索绝对值最大的元素，并进行行和列的交换。这更稳定，但破坏了矩阵的稀疏结构，且需要跟踪列交换以调整解的顺序，实现复杂。部分选主元在实践中的表现通常足够好。

3. 从高斯消元到LU分解：算法的深度解析

细心的你可能已经注意到，在我们带部分选主元的实现中，消元后的矩阵 A 其实包含了两个信息：上三角部分 U 和下三角部分（对角线以下）的乘子。这正是LU分解的雏形。

高斯消元法本质上在不进行行交换的情况下，隐式地完成了 A = LU 的分解，其中 L 是单位下三角矩阵（对角线为1，下三角是消元乘子），U 是上三角矩阵。当引入行交换后，分解关系变为 PA = LU，其中 P 是一个排列矩阵，代表我们的行交换操作。

理解这种等价性有巨大的实用价值：

1. 高效求解多组方程：当需要反复求解 Ax = b₁, Ax = b₂, ... 时，我们只需做一次 LU 分解（O(n³)），然后对每个不同的 b 执行两次三角矩阵的回代（O(n²)），成本远低于每次都做完整消元。
2. 求逆矩阵：矩阵的逆可以通过求解 AX = I 得到，其中 X 的每一列就是逆矩阵的一列。利用 LU 分解可以高效完成。
3. 计算行列式：det(A) = det(P⁻¹) * det(L) * det(U) = (±1) * 1 * (∏ U[i,i])。符号由行交换次数（排列矩阵的奇偶性）决定。

让我们修改之前的函数，使其显式返回 P, L, U，并演示如何用它们来高效求解。

def lu_decomposition_partial_pivoting(A):
    """
    对矩阵A进行带部分选主元的LU分解 (PA = LU)
    参数:
        A: numpy.ndarray, 形状为 (n, n) 的系数矩阵
    返回:
        P: numpy.ndarray, 排列矩阵（或交换索引）
        L: numpy.ndarray, 单位下三角矩阵
        U: numpy.ndarray, 上三角矩阵
    """
    A = A.astype(float).copy()
    n = len(A)
    L = np.eye(n)
    P = np.eye(n)  # 排列矩阵初始化
    pivot_history = list(range(n))

    for k in range(n-1):
        # 选主元
        max_idx = np.argmax(np.abs(A[k:, k])) + k
        if max_idx != k:
            # 交换A的行
            A[[k, max_idx], :] = A[[max_idx, k], :]
            # 交换P的行（记录累积的排列）
            P[[k, max_idx], :] = P[[max_idx, k], :]
            # 交换L中已计算的部分（关键！）
            if k > 0:
                L[[k, max_idx], :k] = L[[max_idx, k], :k]
            pivot_history[k], pivot_history[max_idx] = pivot_history[max_idx], pivot_history[k]

        if np.abs(A[k, k]) < 1e-15:
            raise ValueError("矩阵奇异或近似奇异，无法进行LU分解。")

        # 计算乘子并消元
        for i in range(k+1, n):
            L[i, k] = factor = A[i, k] / A[k, k]
            A[i, k:] = A[i, k:] - factor * A[k, k:]

    U = np.triu(A)
    return P, L, U

def solve_via_lu(P, L, U, b):
    """
    利用LU分解结果PA=LU求解Ax=b
    步骤: 1. 解 L y = P b   (前向替换)
         2. 解 U x = y     (回代)
    """
    n = len(b)
    # 步骤1: 前向替换解 Ly = Pb
    y = np.zeros(n)
    Pb = np.dot(P, b)  # 应用排列矩阵

    for i in range(n):
        y[i] = Pb[i] - np.dot(L[i, :i], y[:i])
        # L是单位下三角，L[i,i]=1，所以不需要除法

    # 步骤2: 回代解 Ux = y
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
    return x

# 演示LU分解的用法
A_lu = np.array([[3., 17., 10.],
                 [2., 4., -2.],
                 [6., 18., -12.]])
b_lu1 = np.array([1., 2., 3.])
b_lu2 = np.array([-1., 0., 1.])

P, L, U = lu_decomposition_partial_pivoting(A_lu.copy())
print("排列矩阵 P:\n", P)
print("下三角矩阵 L:\n", L)
print("上三角矩阵 U:\n", U)
print("验证 PA - LU (应接近零矩阵):\n", np.dot(P, A_lu) - np.dot(L, U))

# 求解第一个方程组
x1 = solve_via_lu(P, L, U, b_lu1)
print("\n解 x1:", x1)
print("验证 A*x1 - b1:", np.dot(A_lu, x1) - b_lu1)

# 求解第二个方程组，复用LU分解结果！
x2 = solve_via_lu(P, L, U, b_lu2)
print("\n解 x2:", x2)
print("验证 A*x2 - b2:", np.dot(A_lu, x2) - b_lu2)

可以看到，在得到 P, L, U 后，求解新的 b 向量成本极低。这是工业级线性系统求解库（如Intel MKL， OpenBLAS）背后的核心思想之一。

4. 工程实践：处理奇异矩阵、性能优化与完整代码

4.1 检测与处理奇异或病态矩阵

在实际应用中，我们遇到的矩阵不总是可逆的。如何优雅地处理奇异或接近奇异的矩阵？

1. 秩亏检测：在消元过程中，如果选主元后找到的主元绝对值仍然小于一个设定的阈值（例如 1e-12），我们可以认为矩阵是数值奇异的。
2. 条件数预警：计算矩阵的条件数（虽然成本较高），如果超过某个阈值（如 1e12），则给出警告，提示用户结果可能不可信。
3. 提供替代方案：对于欠定方程组（无穷多解）或超定方程组（最小二乘解），可以提示用户使用 np.linalg.lstsq 或 np.linalg.pinv 来求伪逆解。

下面是一个增强版的求解器，包含了更健壮的异常处理：

def robust_linear_solver(A, b, tol=1e-12, warn_cond=1e10):
    """
    鲁棒的线性方程组求解器，集成了高斯消元/LU分解与异常处理。
    参数:
        A: 系数矩阵
        b: 常数向量
        tol: 判断主元为零的容差
        warn_cond: 条件数警告阈值
    返回:
        x: 解向量或最小二乘解
        info: 字典，包含状态、条件数等信息
    """
    A = np.array(A, dtype=float, copy=True)
    b = np.array(b, dtype=float, copy=True).flatten()
    n = A.shape[0]

    # 快速形状检查
    if A.shape[0] != A.shape[1]:
        raise ValueError(f"系数矩阵A必须是方阵，当前形状为 {A.shape}")
    if A.shape[0] != b.shape[0]:
        raise ValueError(f"A的行数({A.shape[0]})必须等于b的长度({b.shape[0]})")

    # 可选：检查条件数（对于大矩阵可能较慢）
    try:
        cond_A = np.linalg.cond(A)
        if cond_A > warn_cond:
            print(f"警告：矩阵条件数较大 ({cond_A:.2e})，解可能不准确。")
    except np.linalg.LinAlgError:
        cond_A = np.inf
        print("警告：无法计算矩阵条件数。")

    pivot = list(range(n))
    # 消元过程
    for k in range(n-1):
        # 部分选主元
        max_idx = k + np.argmax(np.abs(A[k:, k]))
        if max_idx != k:
            A[[k, max_idx], :] = A[[max_idx, k], :]
            b[[k, max_idx]] = b[[max_idx], k]
            pivot[k], pivot[max_idx] = pivot[max_idx], pivot[k]

        # 检查主元
        if np.abs(A[k, k]) < tol:
            # 主元过小，尝试寻找非零列？或者判定为秩亏
            # 这里简单抛出异常，更复杂的实现可以尝试继续消元或返回最小二乘解
            raise np.linalg.LinAlgError(f"在步骤 {k} 检测到数值奇异性（主元 < {tol}）。矩阵可能秩亏。")

        # 消元
        for i in range(k+1, n):
            factor = A[i, k] / A[k, k]
            A[i, k+1:] = A[i, k+1:] - factor * A[k, k+1:]
            b[i] = b[i] - factor * b[k]
            A[i, k] = factor  # 存储乘子

    # 检查最后一个主元
    if np.abs(A[-1, -1]) < tol:
        raise np.linalg.LinAlgError("矩阵数值奇异，无法求解唯一解。")

    # 回代
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]

    info = {
        'success': True,
        'condition_number': cond_A,
        'pivot_history': pivot
    }
    return x, info

# 测试奇异矩阵
A_singular = np.array([[1., 2., 3.],
                       [4., 5., 6.],
                       [7., 8., 9.]])  # 第三行是第一行和第二行的线性组合，行列式为0
b_singular = np.array([1., 2., 3.])

try:
    x_sol, info_sol = robust_linear_solver(A_singular, b_singular)
    print("奇异矩阵求解结果:", x_sol)
except np.linalg.LinAlgError as e:
    print(f"求解失败，原因: {e}")
    # 在实际应用中，这里可以fallback到求最小二乘解
    x_pseudo, residuals, rank, s = np.linalg.lstsq(A_singular, b_singular, rcond=None)
    print(f"退而求其次，最小二乘解为: {x_pseudo}")

4.2 性能优化与向量化

我们之前的循环实现清晰易懂，但在处理大规模矩阵时效率不高。NumPy的威力在于向量化操作。我们可以将内层循环替换为向量化操作，显著提升性能。

def gaussian_elimination_vectorized(A, b):
    """使用更多向量化操作的高斯消元法"""
    n = len(A)
    A = A.astype(float).copy()
    b = b.astype(float).copy().reshape(-1, 1)  # 确保b是列向量

    for k in range(n-1):
        # 选主元
        max_idx = k + np.argmax(np.abs(A[k:, k]))
        if max_idx != k:
            A[[k, max_idx], :] = A[[max_idx, k], :]
            b[[k, max_idx], :] = b[[max_idx, k], :]

        if np.abs(A[k, k]) < 1e-15:
            raise np.linalg.LinAlgError("矩阵奇异")

        # 向量化消元：一次操作多行
        i = k+1
        factors = A[i:, k:k+1] / A[k, k]  # 保持二维形状，便于广播
        A[i:, k+1:] = A[i:, k+1:] - factors * A[k, k+1:]
        b[i:, :] = b[i:, :] - factors * b[k, :]
        # 可选：存储乘子
        A[i:, k:k+1] = factors

    # 向量化回代 (从下往上)
    x = np.zeros((n, 1))
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]
    return x.flatten()

对于追求极致性能的场景，我们可以使用SciPy的scipy.linalg.lu_factor和scipy.linalg.lu_solve，它们底层由高度优化的Fortran或C库（如LAPACK）实现。但自己实现一遍的意义在于，当出现问题时，你知道该从哪里入手调试。

4.3 完整、可复用的工业级代码示例

最后，我将提供一个整合了所有最佳实践的完整函数：部分选主元、LU分解格式输出、基本的异常处理以及清晰的文档。你可以直接将它用于你的项目或作为进一步开发的起点。

import numpy as np

def solve_linear_system(A, b, return_lu=False, pivot_tol=1e-12):
    """
    使用带部分选主元的高斯消元法/LU分解求解线性方程组 Ax = b。

    参数
    ----------
    A : array_like, shape (n, n)
        系数矩阵。
    b : array_like, shape (n,) or (n, 1)
        常数向量。
    return_lu : bool, optional
        如果为True，则额外返回P, L, U矩阵。默认为False。
    pivot_tol : float, optional
        判断主元是否为0的容差。默认为1e-12。

    返回
    -------
    x : ndarray, shape (n,)
        方程组的解向量。
    (P, L, U) : tuple of ndarrays (仅在 return_lu=True 时返回)
        P: 排列矩阵，满足 P A = L U
        L: 单位下三角矩阵
        U: 上三角矩阵

    异常
    ------
    LinAlgError
        如果矩阵A是奇异的，或者形状不匹配。

    示例
    --------
    >>> A = [[2, 1, -1], [-3, -1, 2], [-2, 1, 2]]
    >>> b = [8, -11, -3]
    >>> x = solve_linear_system(A, b)
    >>> print(x)
    [ 2.  3. -1.]
    """
    A = np.asarray(A, dtype=np.float64)
    b = np.asarray(b, dtype=np.float64).flatten()

    n, m = A.shape
    if n != m:
        raise np.linalg.LinAlgError(f"矩阵A必须是方阵，当前形状为 ({n}, {m})")
    if n != b.size:
        raise np.linalg.LinAlgError(f"A的行数 ({n}) 必须等于b的长度 ({b.size})")

    # 初始化
    A_work = A.copy()
    b_work = b.copy()
    P = np.eye(n, dtype=np.float64)  # 排列矩阵
    pivot_idx = np.arange(n)         # 行交换记录

    # 高斯消元（LU分解）主循环
    for k in range(n - 1):
        # 1. 部分选主元
        max_row = k + np.argmax(np.abs(A_work[k:, k]))
        if max_row != k:
            # 交换行
            A_work[[k, max_row], :] = A_work[[max_row, k], :]
            b_work[[k, max_row]] = b_work[[max_row], k]
            P[[k, max_row], :] = P[[max_row, k], :]
            pivot_idx[[k, max_row]] = pivot_idx[[max_row, k]]

        # 2. 检查主元
        pivot = A_work[k, k]
        if np.abs(pivot) < pivot_tol:
            raise np.linalg.LinAlgError(
                f"在消元步 {k} 检测到数值奇异性。主元值为 {pivot:.2e}，小于容差 {pivot_tol}。"
            )

        # 3. 计算乘子并消元（向量化）
        i = k + 1
        multipliers = A_work[i:, k] / pivot
        A_work[i:, k+1:] = A_work[i:, k+1:] - multipliers[:, None] * A_work[k, k+1:]
        b_work[i:] = b_work[i:] - multipliers * b_work[k]
        # 将乘子存储在下三角部分
        A_work[i:, k] = multipliers

    # 最终主元检查
    if np.abs(A_work[-1, -1]) < pivot_tol:
        raise np.linalg.LinAlgError("矩阵数值奇异，无法求解。")

    # 从处理后的A_work中提取L和U
    L = np.eye(n)
    L[np.tril_indices(n, -1)] = A_work[np.tril_indices(n, -1)]
    U = np.triu(A_work)

    # 回代求解
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b_work[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]

    if return_lu:
        return x, (P, L, U)
    else:
        return x

# =============== 演示与验证 ===============
if __name__ == "__main__":
    print("=== 测试1: 标准可逆矩阵 ===")
    A1 = np.array([[4, 3, 2],
                   [3, 2, 1],
                   [2, 1, 1]], dtype=float)
    b1 = np.array([10, 7, 5], dtype=float)
    x1 = solve_linear_system(A1, b1)
    print(f"解 x = {x1}")
    print(f"残差范数 ||Ax - b|| = {np.linalg.norm(A1 @ x1 - b1):.2e}")
    print()

    print("=== 测试2: 需要行交换的矩阵 ===")
    A2 = np.array([[0, 2, 3],
                   [4, 5, 6],
                   [7, 8, 9]], dtype=float)
    b2 = np.array([1, 2, 3], dtype=float)
    x2, (P2, L2, U2) = solve_linear_system(A2, b2, return_lu=True)
    print(f"解 x = {x2}")
    print(f"验证 PA - LU 范数: {np.linalg.norm(P2 @ A2 - L2 @ U2):.2e}")
    print()

    print("=== 测试3: 接近奇异的矩阵（应给出警告或失败）===")
    A3 = np.array([[1, 2],
                   [1, 2.00000001]], dtype=float)  # 两行几乎线性相关
    b3 = np.array([3, 3.00000001], dtype=float)
    try:
        x3 = solve_linear_system(A3, b3, pivot_tol=1e-8)
        print(f"解 x = {x3}")
        print(f"残差范数: {np.linalg.norm(A3 @ x3 - b3):.2e}")
        print("注意：由于矩阵接近奇异，解对误差极其敏感。")
    except np.linalg.LinAlgError as e:
        print(f"求解器正确捕获异常: {e}")
        # 使用最小二乘作为后备
        x3_lstsq, *_ = np.linalg.lstsq(A3, b3, rcond=None)
        print(f"最小二乘解: {x3_lstsq}")

这个 solve_linear_system 函数是一个实用的起点。它处理了形状检查、选主元、奇异性检测，并提供了返回LU分解结果的选项。在实际项目中，你可能还需要根据具体情况添加日志、更复杂的误差估计、或者与迭代法结合来处理稀疏矩阵。

自己动手实现核心算法，能让你在调用 np.linalg.solve 时更加自信，因为你知道黑盒里大概发生了什么。当遇到数值问题时，你也有了深入排查的工具和思路。高斯消元法这座桥梁，连接着线性代数的理论世界和工程计算的现实需求，掌握它，是迈向更高级数值计算领域坚实的一步。