笔记6-《程序员的数学3：线性代数》LU分解

目录

1. LU 分解
2. 用 Ruby 验证
3. 数值计算里的两个细节
4. LU 分解在书中的位置

LU 分解

• 3.3 LU分解

  • 核心思想

    • 把矩阵 $A$ 分解成 $A = LU$
      • $L$ 是下三角矩阵，对角线全为 1
      • $U$ 是上三角矩阵
    • 消元法的每一步操作都对应着一个矩阵，$L$ 就是这些操作的累积
    • 分解的巧妙之处：$L$ 和 $U$ 可以拼接成一个矩阵存储
      • 左下部分存 $L$，右上部分存 $U$，对角线是 $U$ 的
      • 大小和原矩阵 $A$ 相同，所以可以在 $A$ 自己的存储空间上覆盖
      • $L$ 对角线的 1 和 $U$ 对角线以下的 0 完全不用存

  • 3.8 意外发生的情况：选主元

    • 如果消元过程中主元为 0，就无法顺利分解成 $A = LU$
    • 办法是交换行，对应左乘一个置换矩阵 $S$
    • 多次选主元的结果：$A = PLU$
      • $P = S_1 S_2 \dots S_k$
      • 每一行恰好一个 1、每一列恰好一个 1、其余为 0 的矩阵叫置换矩阵
      • 置换矩阵乘以向量，效果就是重新排列分量的顺序
      • 任何置换矩阵都能写成若干 $S$ 的乘积

  • 3.6 利用LU分解求解线性方程组

    • 把 $LU\mathbf{x} = \mathbf{y}$ 翻译一下：先乘以 $U$，再乘以 $L$
    • 于是分两步
      1. 求满足 $L\mathbf{z} = \mathbf{y}$ 的 $\mathbf{z}$
      2. 求满足 $U\mathbf{x} = \mathbf{z}$ 的 $\mathbf{x}$
    • 合起来 $A\mathbf{x} = LU\mathbf{x} = L(U\mathbf{x}) = L\mathbf{z} = \mathbf{y}$
    • 两步都是三角方程组，都很好解
    • 注意拼接布局的细节：对角线那一格存的是 $U$ 的对角元素，$L$ 的对角恒为 1 不存
      • 前代解 $L\mathbf{z} = \mathbf{y}$ 时对角系数按 1 算
      • 回代解 $U\mathbf{x} = \mathbf{z}$ 时才用 $U$ 的对角元素做除数

  • 3.6 顺带得到行列式

    • $\det A = \det(LU) = (\det L)(\det U)$

    • 下三角和上三角矩阵的行列式就是对角元素的乘积

    • 这里 $\det L$ 恰好等于 1（对角线全是 1），所以

      $\det A = (\text{U 的对角元素的乘积})$

    • 有选主元时：$\det A = (\det P)(\det U \text{ 的对角元素乘积})$

      • $\det P = +1$（选主元次数为偶数）或 $-1$（奇数）

  • 3.7 利用LU分解求逆矩阵

    • 把 $A^{-1} = X = (\mathbf{x}_1, \dots, \mathbf{x}_n)$ 按列分块
    • $AX = I$ 写成 $A(\mathbf{x}_i, \dots, \mathbf{x}_n) = (\mathbf{e}_1, \dots, \mathbf{e}_n)$
    • 每个等式 $A\mathbf{x}_i = \mathbf{e}_i$ 都是一个解线性方程的问题
    • 左边完全相同，右边不同 → 只需做一次 LU 分解，之后一劳永逸
    • 但实践中往往不需要整个 $A^{-1}$，只需要 $A^{-1}\mathbf{y}$
      • 那就直接解 $A\mathbf{x} = \mathbf{y}$，运算量更小，也避免累积误差

• 为什么 LU 分解值得单独一章

  • 一次分解，多次使用
    • 解 $n$ 个不同右端的方程组：消元只做一次
    • 求行列式：只算 $U$ 的对角乘积
    • 求逆矩阵：对每一列标准向量解一次
  • 这就是「把昂贵的部分提前做掉」的典型

用 Ruby 验证

# LU 分解（不选主元），结果覆盖写入原矩阵
def lu_decomp(mat)
  n = mat.size
  m = mat.map { |r| r.map(&:to_f) }
  (0...n).each do |i|
    (i + 1...n).each do |r|
      f = m[r][i] / m[i][i]
      m[r][i] = f                    # L 的对角线以下
      (i + 1...n).each { |c| m[r][c] -= f * m[i][c] }
    end
  end
  m
end

# 用 LU 分解解方程组
def solve_lu(a, y)
  n = a.size
  m = lu_decomp(a)

  # 前代：Lz = y，L 的对角线恒为 1（拼接时对角线存的是 U 的）
  z = Array.new(n, 0.0)
  n.times do |i|
    s = y[i] - (0...i).sum { |c| m[i][c] * z[c] }
    z[i] = s
  end

  # 回代：Ux = z
  x = Array.new(n, 0.0)
  (n - 1).downto(0) do |i|
    s = z[i] - ((i + 1)...n).sum { |c| m[i][c] * x[c] }
    x[i] = s / m[i][i]
  end
  x
end

a = [[2.0, 1.0, -1.0], [-3.0, -1.0, 2.0], [-2.0, 1.0, 2.0]]
y = [8.0, -11.0, -3.0]

x = solve_lu(a, y)
puts "解 = #{x.map { |v| v.round(6) }}"   # 解 = [2.0, 3.0, -1.0]

# 验证 A x = y
puts "A*x = #{a.each_index.map { |i| a[i].each_index.sum { |j| a[i][j] * x[j] }.round(6) }}"  # A*x = [8.0, -11.0, -3.0]

# 顺带算行列式：U 的对角元素乘积（L 的对角全是 1）
u = lu_decomp(a)
det = (0...3).inject(1.0) { |acc, i| acc * u[i][i] }
puts "det = #{det.round(6)}"            # det = -1.0

数值计算里的两个细节

选主元不是为了正确性，而是为了稳定性。

消元时如果主元很小（但不是 0），除法会把舍入误差放大。选绝对值最大的行做主元，能让后续计算的相对误差保持在可接受范围内。代价是要额外维护一个置换矩阵 $P$，以及把 $P^T\mathbf{y}$ 算一遍（只是重排分量，很便宜）。

存储技巧有实际价值。

$L$ 和 $U$ 不需要分开存。把它们拼在一个矩阵里：

[ L 的下三角 | U 的上三角 ]
     对角线是 U 的，L 的对角线恒为 1 不存

好处是内存占用和原矩阵完全一样，而且可以原地覆盖——不需要额外分配空间。对于大矩阵求解，这个差别很实际。

LU 分解在书中的位置

这一章在书里的标题是「计算机上的计算（1）」，说明作者把它定位成从数学到工程的转折点：

• 第 1、2 章是数学：矩阵是什么、映射是什么、解有什么性质

• 第 3 章开始是计算：怎么在有限精度下真的把它算出来

• 第 4 章（特征值、对角化）继续这条路：判断系统会不会「失控」

上一节我们看到，$A\mathbf{x} = \mathbf{y}$ 只要做一次 LU 分解就能解很多次。而当 $\mathbf{y}$ 会随时间变化时，问题就变成了「$A^t$ 会不会爆炸」，这就需要第 4 章的特征值工具。