LU分解

LU分解

小象实战讲义 · 人工智能数学基础

在线性代数的计算中，求解线性方程组是最核心的任务之一。高斯消元法是大家熟知的方法，但其计算过程缺乏“可复用性”。本节介绍的LU分解，本质上是将高斯消元法的过程“固化”下来，将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积。这种分解不仅揭示了矩阵的内在结构，更在求解多个具有相同系数矩阵的方程组、计算行列式等方面展现出巨大的效率优势。学完本节，你将掌握LU分解的算法原理、实现细节及其核心应用场景。

💡 核心导读

本节将围绕以下几个核心要点展开：

1. 定义与唯一性：理解LU分解（$A = LU$）的数学定义，并认识到在未加约束时，这种分解并不唯一。
2. 算法本质：揭示LU分解与高斯消元法的等价关系，理解通过一系列下三角初等矩阵左乘，将矩阵$A$化为上三角矩阵$U$的过程，而$L$正是这些初等矩阵乘积的逆。
3. 选主元（Pivoting）：解决基本LU分解中主元可能为零的“尴尬”问题，引入行置换操作，得到更稳定、通用的PLU分解（$PA = LU$）。
4. 复杂度与应用：分析LU分解的时间与空间复杂度，并掌握其在线性方程组求解、行列式计算等任务中的高效应用方法。

LU分解的定义与基本思想

对于一个 $n$ 阶方阵 $A$，如果存在一个下三角矩阵 $L$ 和一个上三角矩阵 $U$，使得： $$A = LU$$ 则称此等式为矩阵 $A$ 的 LU分解。其中 $L$ 代表 Lower triangular matrix，$U$ 代表 Upper triangular matrix。

例如，一个 $3 \times 3$ 的矩阵可以分解为： $$A = \begin{pmatrix} a_{11} & a_{12} & a_{13} \ a_{21} & a_{22} & a_{23} \ a_{31} & a_{32} & a_{33} \end{pmatrix} = \begin{pmatrix} l_{11} & 0 & 0 \ l_{21} & l_{22} & 0 \ l_{31} & l_{32} & l_{33} \end{pmatrix} \begin{pmatrix} u_{11} & u_{12} & u_{13} \ 0 & u_{22} & u_{23} \ 0 & 0 & u_{33} \end{pmatrix}$$

如果仅凭这个定义，通过对应元素相等来列方程求解，我们会发现未知数个数（$L$和$U$的非零元素之和）多于方程个数（$A$的元素个数）。这意味着，如果不加任何额外约束，LU分解不是唯一的。例如，你可以将$L$的对角线元素全部缩放，同时将$U$的对角线元素进行相反的缩放，仍然满足 $A=LU$。

LU分解的算法本质：高斯消元法

LU分解的构造性算法直接源于高斯消元法。回顾高斯消元过程：我们通过一系列初等行变换（用某行的倍数加到另一行），将系数矩阵 $A$ 化为上三角矩阵 $U$。

每一步消元操作都对应左乘一个下三角初等矩阵。例如，要消去第一列中 $a_{21}$ 这个元素，我们执行的操作是：第二行加上第一行的 $(-\frac{a_{21}}{a_{11}})$ 倍。这个行变换对应的初等矩阵 $E_{21}$ 为： $$E_{21} = \begin{pmatrix} 1 & 0 & 0 & \dots & 0 \ -\frac{a_{21}}{a_{11}} & 1 & 0 & \dots & 0 \ 0 & 0 & 1 & \dots & 0 \ \vdots & \vdots & \vdots & \ddots & \vdots \ 0 & 0 & 0 & \dots & 1 \end{pmatrix}$$ 显然，$E_{21}$ 是一个下三角矩阵。

依次用主元消去其下方的所有元素，相当于左乘一系列这样的下三角初等矩阵 $E_{21}, E_{31}, \dots, E_{n1}, E_{32}, \dots$。设它们的乘积为 $E$，则有： $$E A = U$$ 由于所有 $E_{ij}$ 都是下三角矩阵，它们的乘积 $E$ 仍然是下三角矩阵。因此，我们可以得到： $$A = E^{-1} U$$ 令 $L = E^{-1}$。可以证明，下三角矩阵的逆仍然是下三角矩阵，且下三角矩阵的乘积也是下三角矩阵。因此，$L$ 是一个下三角矩阵。这样，我们就从高斯消元的过程中自然地得到了 $A$ 的一个 LU 分解：$A = LU$。

在实际算法实现中，我们通常采用一种“原地”（in-place）的紧凑格式。初始化 $U = A$，$L$ 为单位矩阵 $I$。然后按列进行高斯消元：

1. 对于第 $k$ 列（$k = 1, 2, \dots, n-1$）。
2. 对于第 $i$ 行（$i = k+1, k+2, \dots, n$），计算乘数 $l_{ik} = u_{ik} / u_{kk}$，并将 $l_{ik}$ 存储在原 $A$ 矩阵的 $(i, k)$ 位置（这个位置在消元后 $U$ 中应为0）。
3. 用这个乘数更新第 $i$ 行：$U[i, k:] = U[i, k:] - l_{ik} \cdot U[k, k:]$。

算法结束后，$U$ 的上三角部分（包括对角线）存储了 $U$ 矩阵，而 $L$ 矩阵的严格下三角部分（即对角线以下）存储在了 $U$ 原来这些被消为零的位置。由于我们约定 $L$ 的对角线元素为1（单位下三角矩阵），所以不需要额外存储。这样，空间复杂度为 $\Theta(n^2)$，恰好存储原矩阵 $A$。

import numpy as np

def lu_decomposition_basic(A):
    """
    基本的LU分解（不选主元）。
    假设A是numpy数组表示的方阵。
    返回 L (单位下三角), U (上三角) 使得 A = L @ U。
    """
    n = A.shape[0]
    U = A.copy().astype(float) # 工作矩阵，最终存放U
    L = np.eye(n)              # 初始化L为单位矩阵

    for k in range(n-1):       # 第k列为主元列
        if U[k, k] == 0:
            # 主元为零，基本算法失效
            raise ValueError(f"Zero pivot encountered at position ({k},{k}). Basic LU decomposition fails.")
        for i in range(k+1, n):
            # 计算乘数
            L[i, k] = U[i, k] / U[k, k]
            # 更新U的第i行
            U[i, k:] = U[i, k:] - L[i, k] * U[k, k:]
            # 显式地将U中已消元位置置零（非必须，但更清晰）
            # U[i, k] = 0.0
    return L, U

# 示例
A = np.array([[2, 4, 1],
              [4, 9, 3],
              [1, 3, 2]], dtype=float)

L, U = lu_decomposition_basic(A)
print("矩阵 A:")
print(A)
print("\n分解得到的 L:")
print(L)
print("\n分解得到的 U:")
print(U)
print("\n验证 A - L@U (应接近零矩阵):")
print(A - L @ U)

选主元LU分解（PLU分解）

基本LU分解算法存在一个致命问题：当主对角线元素 $u_{kk}$（主元）为0时，计算乘数 $l_{ik}$ 会导致除以零的错误。即使主元不为零但非常接近零，也会在数值计算中引入巨大的舍入误差，导致结果不稳定。

为了解决这个问题，我们引入选主元（Pivoting）技术。最常用的是部分选主元：在消去第 $k$ 列之前，从该列第 $k$ 行到第 $n$ 行中找出绝对值最大的元素，将其所在行与第 $k$ 行交换。这个行交换操作对应左乘一个置换矩阵 $P$。

因此，整个消元过程变为：我们寻找一个置换矩阵 $P$，以及单位下三角矩阵 $L$ 和上三角矩阵 $U$，使得： $$PA = LU$$ 这就是 PLU分解。它对于任何非奇异矩阵都是存在的，并且通过选主元保证了数值稳定性。

改进后的算法流程如下：

1. 初始化 $U = A$, $L = I$, $P = I$。
2. 对于每一列 $k$： a. 选主元：在第 $k$ 列中，从第 $k$ 行到最后一行，找到绝对值最大的元素所在的行 $p$。 b. 交换：交换 $U$ 的第 $k$ 行和第 $p$ 行；交换 $L$ 的第 $k$ 行和第 $p$ 行（但只交换 $k$ 列之前的已计算部分）；记录置换，更新 $P$。 c. 消元：进行与基本LU分解相同的消元步骤。
3. 算法结束，得到 $P, L, U$。

选主元操作（找最大值和交换行）的复杂度是 $O(n^2)$，而消元的主体部分是 $O(n^3)$。因此，PLU分解的总体时间复杂度仍然是 $\frac{2}{3}n^3$ 量级（精确常数约为 $\frac{2}{3}n^3$）。

def lu_decomposition_with_pivoting(A):
    """
    带部分选主元的LU分解 (PLU分解)。
    返回 P (置换矩阵), L (单位下三角), U (上三角) 使得 P @ A = L @ U。
    """
    n = A.shape[0]
    U = A.copy().astype(float)
    L = np.eye(n)
    P = np.eye(n) # 置换矩阵

    for k in range(n-1):
        # 1. 选主元：在第k列，从第k行开始找最大绝对值
        pivot_row = k + np.argmax(np.abs(U[k:, k]))
        if pivot_row != k:
            # 2. 交换行
            # 交换U
            U[[k, pivot_row], k:] = U[[pivot_row, k], k:]
            # 交换L中已计算的部分（前k-1列）
            if k > 0:
                L[[k, pivot_row], :k] = L[[pivot_row, k], :k]
            # 记录置换
            P[[k, pivot_row], :] = P[[pivot_row, k], :]

        # 检查主元是否仍为0（理论上选主元后不会，但数值上可能）
        if np.abs(U[k, k]) < 1e-15:
            print(f"警告：第{k}列选主元后主元仍接近零，矩阵可能奇异。")

        # 3. 消元
        for i in range(k+1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] = U[i, k:] - L[i, k] * U[k, k:]

    return P, L, U

# 演示选主元的必要性
A_bad = np.array([[0, 2, 3],
                  [1, 1, 1],
                  [4, 5, 6]], dtype=float)

print("尝试对‘坏’矩阵进行基本LU分解:")
try:
    L_basic, U_basic = lu_decomposition_basic(A_bad)
except ValueError as e:
    print(f"基本LU分解失败: {e}")

print("\n使用选主元PLU分解:")
P, L_pivot, U_pivot = lu_decomposition_with_pivoting(A_bad)
print("置换矩阵 P:")
print(P)
print("单位下三角矩阵 L:")
print(L_pivot)
print("上三角矩阵 U:")
print(U_pivot)
print("\n验证 P@A - L@U (应接近零矩阵):")
print(P @ A_bad - L_pivot @ U_pivot)

LU分解的应用

1. 求解线性方程组

对于方程组 $A\mathbf{x} = \mathbf{b}$，如果已经得到 $A$ 的 PLU 分解 $PA = LU$，则原方程等价于： $$LU\mathbf{x} = P\mathbf{b}$$ 令 $\mathbf{y} = U\mathbf{x}$，则可以分两步求解：

1. 前向代入：解下三角方程组 $L\mathbf{y} = P\mathbf{b}$。
   • 因为 $L$ 是单位下三角矩阵，求解非常高效，复杂度为 $O(n^2)$。 $$\begin{aligned} y_1 &= (P\mathbf{b})$$1 \ y_2 &= (P\mathbf{b})2 - l{21}y_1 \ &\vdots \ y_n &= (P\mathbf{b})n - \sum{j=1}^{n-1} l{nj} y_j \end{aligned}
2. 后向代入：解上三角方程组 $U\mathbf{x} = \mathbf{y}$。
   • 同样高效，复杂度为 $O(n^2)$。 $$\begin{aligned} x_n &= y_n / u_{nn} \ x_{n-1} &= (y_{n-1} - u_{n-1,n}x_n) / u_{n-1,n-1} \ &\vdots \ x_1 &= (y_1 - \sum_{j=2}^{n} u_{1j}x_j) / u_{11} \end{aligned}$$

当需要求解多个具有相同系数矩阵 $A$ 但不同右端项 $\mathbf{b}$ 的方程组时，LU分解的优势尤为明显。我们只需进行一次 $O(n^3)$ 的分解，之后对每个新的 $\mathbf{b}$，只需进行两次 $O(n^2)$ 的代入求解即可。

2. 计算行列式

由于 $\det(A) = \det(L)\det(U)$（对于PLU分解，$\det(P) = \pm 1$，需考虑符号），且三角矩阵的行列式等于其对角线元素的乘积。对于单位下三角矩阵 $L$，$\det(L)=1$。因此：

• 对于基本LU分解：$\det(A) = \prod_{i=1}^{n} u_{ii}$。
• 对于PLU分解：$\det(A) = \det(P^{-1}LU) = \det(P^{-1})\det(L)\det(U) = \operatorname{sign}(P) \prod_{i=1}^{n} u_{ii}$，其中 $\operatorname{sign}(P)$ 由置换的奇偶性决定（行交换次数为奇数次则为 -1，偶数次则为 1）。

这提供了另一种计算行列式的方法，其复杂度与LU分解相同，为 $O(n^3)$，比直接按定义计算高效得多。

def solve_via_plu(A, b):
    """使用PLU分解求解线性方程组 A x = b."""
    P, L, U = lu_decomposition_with_pivoting(A)
    n = A.shape[0]
    # 步骤1: 前向代入解 L y = P b
    y = np.zeros(n)
    Pb = P @ b
    for i in range(n):
        # 注意L是单位下三角，对角线为1
        y[i] = Pb[i] - np.dot(L[i, :i], y[:i])
    # 步骤2: 后向代入解 U x = 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

# 应用示例：求解方程组并计算行列式
A_sys = np.array([[1, 2, 4],
                  [3, 8, 14],
                  [2, 6, 13]], dtype=float)
b_sys = np.array([3, 13, 4])

x_solution = solve_via_plu(A_sys, b_sys)
print("线性方程组 A x = b 的解 x:")
print(x_solution)
print("\n验证 A @ x - b (应接近零向量):")
print(A_sys @ x_solution - b_sys)

# 计算行列式
P_det, L_det, U_det = lu_decomposition_with_pivoting(A_sys)
det_A = np.linalg.det(P_det) * np.prod(np.diag(U_det)) # det(P) = +/-1, det(L)=1
print(f"\n通过PLU分解计算的行列式 det(A) = {det_A}")
print(f"使用numpy.linalg.det验证: {np.linalg.det(A_sys)}")

📝 动手练一练

1. 手动推导：对以下矩阵 $A$ 进行基本的LU分解（即 $A=LU$，$L$ 为单位下三角矩阵）。请写出每一步消元过程对应的乘数，并给出最终的 $L$ 和 $U$ 矩阵。 $$A = \begin{pmatrix} 2 & 1 & 1 \ 4 & 3 & 3 \ 8 & 7 & 9 \end{pmatrix}$$

   参考答案： 消元步骤：

   • 第1列：用第1行消去第2行和第3行。乘数 $l_{21}=4/2=2$, $l_{31}=8/2=4$。 $U$ 更新为 $\begin{pmatrix}2 & 1 & 1 \ 0 & 1 & 1 \ 0 & 3 & 5\end{pmatrix}$。
   • 第2列：用第2行消去第3行。乘数 $l_{32}=3/1=3$。 $U$ 最终为 $\begin{pmatrix}2 & 1 & 1 \ 0 & 1 & 1 \ 0 & 0 & 2\end{pmatrix}$。 因此， $$L = \begin{pmatrix}1 & 0 & 0 \ 2 & 1 & 0 \ 4 & 3 & 1\end{pmatrix}, \quad U = \begin{pmatrix}2 & 1 & 1 \ 0 & 1 & 1 \ 0 & 0 & 2\end{pmatrix}$$ 验证 $L \times U = A$。

2. 编程验证：使用上面提供的 lu_decomposition_with_pivoting 函数，对一个 $5 \times 5$ 的随机矩阵（np.random.randn(5,5)）进行PLU分解。然后：

   • 验证 $PA - LU$ 是否接近零矩阵。
   • 随机生成5个不同的右端项向量 $\mathbf{b}_1, \dots, \mathbf{b}_5$，利用已分解好的 $P, L, U$，分别求解 $A\mathbf{x}_i = \mathbf{b}_i$。比较直接调用 np.linalg.solve 和你的分步求解函数 solve_via_plu 的结果差异（计算范数误差）。

   参考答案：（代码框架）

   import numpy as np
   np.random.seed(42) # 固定随机种子以便复现
   A = np.random.randn(5, 5)
   P, L, U = lu_decomposition_with_pivoting(A)
   # 验证分解
   print("PA - LU 的Frobenius范数:", np.linalg.norm(P@A - L@U))
   # 求解多个方程组
   solutions_plu = []
   solutions_np = []
   for i in range(5):
       b = np.random.randn(5)
       x_plu = solve_via_plu(A, b) # 需要先定义这个函数
       x_np = np.linalg.solve(A, b)
       solutions_plu.append(x_plu)
       solutions_np.append(x_np)
       err = np.linalg.norm(x_plu - x_np)
       print(f"第{i+1}个方程解误差: {err}")

本章小结

本节深入探讨了线性代数中一个强大而实用的工具——LU分解。

要点回顾：

• 定义：LU分解将方阵 $A$ 表示为下三角矩阵 $L$ 和上三角矩阵 $U$ 的乘积，即 $A = LU$。基本分解不唯一，常约定 $L$ 为单位下三角矩阵。
• 算法本质：LU分解是高斯消元法的“矩阵化”表述。消元过程等价于左乘一系列下三角初等矩阵，其逆的乘积即为 $L$，而最终得到的上三角矩阵就是 $U$。
• 选主元（PLU）：为避免零主元或小主元导致的数值问题，引入了行交换（置换矩阵 $P$），得到稳定的分解 $PA = LU$。这是实际计算中的标准方法。
• 复杂度：分解的算术复杂度为 $O(n^3)$（约 $\frac{2}{3}n^3$ 次浮点运算），空间复杂度为 $O(n^2)$。
• 核心应用：
  1. 高效求解线性方程组：将 $O(n^3)$ 的消元过程固化为分解，后续对每个新的右端项只需 $O(n^2)$ 的前向/后向代入。
  2. 便捷计算行列式：$\det(A) = \pm \prod u_{ii}$（符号由置换决定）。

行动清单：

1. 理解与推导：任选一个 $3\times3$ 矩阵，手动模拟一遍带选主元的PLU分解过程，写出每一步的置换和消元乘数，加深对算法流程的理解。
2. 代码实现：不依赖示例代码，尝试自己独立编写 lu_decomposition_with_pivoting 函数。挑战自己处理行交换时 $L$ 矩阵已计算部分的正确交换逻辑。
3. 拓展思考：LU分解要求矩阵是方阵。思考一下，对于非方阵的 $m \times n$ 矩阵，是否存在类似的三角分解？这引出了下一节可能讨论的QR分解。

— 小象教研组