颜力刚Ligang Yan

· NSCC

手算 LU 分解：把一个 3×3 矩阵拆成 L·U

LU 分解就是把高斯消元的过程存下来：U 是消元的结果，L 记录每一步减了多少倍。拿一个 3×3 的矩阵从头算一遍，再用它解方程、算行列式，讲清楚为什么要换行（PA = LU），最后说说它到底有什么用。

线性代数LU分解机器学习NumPy

English version: LU Decomposition by Hand: Splitting a 3×3 Matrix into L·U

LU 分解说的是：一个方阵 A 可以拆成两个三角矩阵相乘。

A   =   L   ·   U
3×3    3×3     3×3
部分 形状 含义
L 下三角，对角线全是 1 记录消元时每一步减了多少倍
U 上三角 高斯消元做完后剩下的那个阶梯形矩阵

一句话：LU 分解 = 把高斯消元的过程存下来。 U 是消元的结果，L 是消元的操作记录。

题目

A = [ 2  1  1 ]
    [ 4  3  3 ]
    [ 8  7  9 ]

求 L 和 U，使得 A = LU。

思路：高斯消元，把倍数存下来

先用第 1 行把第 1 列下面的数消成 0，再用第 2 行把第 2 列下面的数消成 0。每次「第 i 行减去第 j 行的 m 倍」里的 m，就是 L 的 (i, j) 位置上的数。

m_ij = 要消掉的那个数 / 主元      →  填进 L 的 (i, j)

第 1 步：消第 1 列

主元是 a₁₁ = 2。

m₂₁ = 4 / 2 = 2     行2 − 2·行1 = [4 3 3] − [4 2 2] = [0 1 1]
m₃₁ = 8 / 2 = 4     行3 − 4·行1 = [8 7 9] − [8 4 4] = [0 3 5]

消完之后：

[ 2  1  1 ]
[ 0  1  1 ]
[ 0  3  5 ]

第 2 步：消第 2 列

主元是第 2 行第 2 列的 1。

m₃₂ = 3 / 1 = 3     行3 − 3·行2 = [0 3 5] − [0 3 3] = [0 0 2]

消完就得到 U：

U = [ 2  1  1 ]
    [ 0  1  1 ]
    [ 0  0  2 ]

第 3 步：把倍数填进 L

对角线是 1，下三角位置放刚才的三个倍数：

L = [ 1  0  0 ]     ← m₂₁ = 2 在 (2,1)
    [ 2  1  0 ]     ← m₃₁ = 4 在 (3,1)，m₃₂ = 3 在 (3,2)
    [ 4  3  1 ]

检查：乘回去。

LU 第 1 行 = 1·[2 1 1]                          = [2 1 1]   ✓
LU 第 2 行 = 2·[2 1 1] + 1·[0 1 1]              = [4 3 3]   ✓
LU 第 3 行 = 4·[2 1 1] + 3·[0 1 1] + 1·[0 0 2]  = [8 7 9]   ✓

为什么 L 能直接填进去？ LU 的第 3 行 = 4·(U 第 1 行) + 3·(U 第 2 行) + 1·(U 第 3 行)。这正好是把消元的操作倒着做一遍：消元时从第 3 行里减掉了 4 倍的第 1 行和 3 倍的第 2 行，现在再加回去，就还原成原来的第 3 行。L 的每一行，就是「原来这一行是哪几个 U 行的组合」。

用 LU 解方程 Ax = b

A = LU，所以 Ax = b 变成 L(Ux) = b。令 y = Ux，分两步：

  1. 解 Ly = b。L 是下三角，从上往下代入，叫前代。
  2. 解 Ux = y。U 是上三角，从下往上代入，叫回代。

取 b = (4, 10, 24)ᵀ：

前代 Ly = b：
  y₁ = 4
  2·y₁ + y₂ = 10          →  y₂ = 10 − 8 = 2
  4·y₁ + 3·y₂ + y₃ = 24   →  y₃ = 24 − 16 − 6 = 2

回代 Ux = y：
  2·x₃ = 2                →  x₃ = 1
  x₂ + x₃ = 2             →  x₂ = 1
  2·x₁ + x₂ + x₃ = 4      →  x₁ = 1

x = (1, 1, 1)ᵀ。检查：A·(1, 1, 1) = (2+1+1, 4+3+3, 8+7+9) = (4, 10, 24)，对。

顺带的好处：行列式

det(A) = det(L) · det(U) = 1 · (2·1·2) = 4。L 的对角线全是 1，所以 det(A) 就是 U 对角线的乘积。

遇到主元是 0：要换行

A = [ 0  2 ]
    [ 1  3 ]

第一个主元是 0，m₂₁ = 1/0 算不出来。办法是先把两行对调，再消元：

PA = [ 0  1 ] [ 0  2 ]  =  [ 1  3 ]  =  [ 1  0 ] [ 1  3 ]  =  L · U
     [ 1  0 ] [ 1  3 ]     [ 0  2 ]     [ 0  1 ] [ 0  2 ]

P 是置换矩阵，记录换了哪些行。这时分解写成 PA = LU。

主元不是 0 但很小时，除以它会把误差放大，所以软件默认每一步都在当前列里挑绝对值最大的数当主元（部分主元法）。也就是说，实际用的分解几乎总是 PA = LU。

容易出错的地方

  • m 填错位置。 m₃₂ 是「第 3 行减第 2 行的倍数」，放在 (3,2)，不是 (2,3)。
  • 用原矩阵的数算 m。 第 2 步的分子 3 是第 1 步消完之后的 a₃₂，不是原矩阵里的 7。
  • 忘了 L 的对角线是 1。 对角线上的数已经被「提」到 U 里了。
  • 前代和回代的方向搞反。 L 是从上往下，U 是从下往上。

用 NumPy 核对

import numpy as np
from scipy.linalg import lu

A = np.array([[2, 1, 1], [4, 3, 3], [8, 7, 9]], dtype=float)
P, L, U = lu(A)          # scipy 的约定：A = P @ L @ U

P    # [[0. 1. 0.]
     #  [0. 0. 1.]
     #  [1. 0. 0.]]
L    # [[1.     0.     0.]
     #  [0.25   1.     0.]
     #  [0.5    0.6667 1.]]
U    # [[ 8.    7.    9.   ]
     #  [ 0.   -0.75 -1.25 ]
     #  [ 0.    0.   -0.6667]]

结果和手算不一样，这是正常的：scipy 用了部分主元，第 1 列最大的是 8，它先把第 3 行换到最上面。分解不唯一，但乘回去都是 A，行列式也对得上：8·(−0.75)·(−2/3) = 4，再乘上 P 的符号 +1（三行循环移动，是偶置换），还是 4。

两个小坑：

  • scipy 的约定是 A = P @ L @ U，不是 PA = LU。两个写法差一个转置：P 要换成 Pᵀ。
  • 想要不换行的 L 和 U，只能自己写。三行代码的事：对每一列 j，对它下面的每一行 i，算 m = U[i,j] / U[j,j]，存进 L[i,j]，再做 U[i] -= m * U[j]。

算这个有什么用

A = LU 本身不是目的。它的价值是：把「消元」这件最贵的事做一次，之后反复使用。

1. 同一个 A，很多个 b

解 Ax = b 最直接的办法是每次都从头消元，每次约 n³/3 次乘法。有了 LU，消元只做一次，之后每个 b 只要一次前代加一次回代，约 n² 次。

n = 1000，要解 100 个不同的 b：

做法 运算量（乘法次数的量级）
每次都重新消元 100 × 3.3×10⁸ ≈ 3.3×10¹⁰
LU 一次 + 100 次代入 3.3×10⁸ + 100 × 10⁶ ≈ 4.3×10⁸

差了大约 80 倍。实际场景很常见：时间序列里每一步的矩阵相同、只有右边在变；迭代算法里每轮都要解同一个系统。

2. 求逆其实就是解 n 个方程

A⁻¹ 的第 j 列，就是解 Ax = eⱼ（单位矩阵的第 j 列）。所以求逆 = 做一次 LU + n 次代入。

更要紧的是反过来的结论：如果只是想算 A⁻¹b，不要真的去求逆，直接 LU 解方程。 又快，数值上也更稳。np.linalg.solve(A, b) 内部做的就是 LU，不是先求逆再乘。

3. 行列式几乎白送

det(A) = U 对角线的乘积，再乘一个 +1 或 −1（取决于换了几次行）。用余子式展开是 O(n!)，用 LU 是 O(n³)。

4. 它是更多算法的零件

  • 线性回归的正规方程 (XᵀX)β = Xᵀy：解一个线性方程组。XᵀX 是对称正定的，可以用 LU 的特殊版本 Cholesky 分解，再快一倍。
  • 牛顿法：每一步要解 J·Δ = −f，J 是雅可比矩阵。如果同一个 J 重复用好几步（弦法这类变体），就只分解一次。
  • 高斯过程、卡尔曼滤波：里面都有「用协方差矩阵解方程」这一步，底下是 Cholesky / LU。
  • 微分方程的隐式求解：每个时间步都要解一个形状相同的系统。

写成一句话：只要你看到「解 Ax = b」，而且 A 会重复出现，就该想到 LU。

小作业

  1. 手算。 对 A = [[1,2,1],[3,8,7],[2,7,9]] 求 L、U，乘回去验证，并算出 det(A)。
  2. 解方程。 用第 1 题的 L、U 解 Ax = (1, 2, 3)ᵀ，先前代再回代，再代回 A 检查。
  3. 写代码。 用纯 NumPy 实现 my_lu(A)，返回不换行的 (L, U)，用本文的例子测试 np.allclose(L @ U, A)。再试「遇到主元是 0」那一节的矩阵，看看会发生什么。
参考答案（先自己做，再看）
  1. m₂₁ = 3，m₃₁ = 2；消完第 1 列后第 2、3 行是 [0,2,4]、[0,3,7]；m₃₂ = 3/2。 L = [[1,0,0],[3,1,0],[2,1.5,1]]，U = [[1,2,1],[0,2,4],[0,0,1]]，det = 1·2·1 = 2。
  2. 前代 y = (1, −1, 2.5)，回代 x = (9.5, −5.5, 2.5)。检查：A·x = (1, 2, 3)。
  3. 主元是 0 的矩阵会在 U[j,j] == 0 处除零，得到 inf 或 nan。