矩阵分解：LU分解与Cholesky分解

引言

在数值分析中，矩阵分解是解决线性方程组的重要方法。通过将复杂的矩阵分解为结构简单的矩阵乘积，我们可以更高效地求解线性方程组、计算矩阵的逆以及求解特征值问题。本文将详细介绍两种重要的矩阵分解方法：LU分解和Cholesky分解。

1. LU分解

1.1 定义

LU分解（LU Decomposition）是将一个矩阵分解为一个下三角矩阵L（Lower triangular matrix）和一个上三角矩阵U（Upper triangular matrix）的乘积。

对于 $n \times n$ 矩阵 $A$，如果存在下三角矩阵 $L$ 和上三角矩阵 $U$，使得：

$$A = LU $$

其中：

  • $L$ 是下三角矩阵：$L = \begin{pmatrix} l_{11} & 0 & \cdots & 0 \\ l_{21} & l_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{pmatrix}$

  • $U$ 是上三角矩阵：$U = \begin{pmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{pmatrix}$

1.2 存在性和唯一性

定理：设 $A$ 是 $n \times n$ 矩阵，如果 $A$ 的所有顺序主子式都非零，即：

$$\det(A_k) \neq 0, \quad k = 1, 2, \ldots, n-1 $$

其中 $A_k$ 是 $A$ 的前 $k$ 行前 $k$ 列构成的子矩阵，则 $A$ 存在唯一的LU分解，且可以规定 $L$ 的对角元素全为1（称为Doolittle分解）或 $U$ 的对角元素全为1（称为Crout分解）。

1.3 高斯消元法与LU分解

LU分解实际上是高斯消元法的矩阵表示。在高斯消元过程中：

  • $L$ 记录了消元过程中的乘数
  • $U$ 是消元后得到的上三角矩阵

1.4 Doolittle分解算法

Doolittle分解假设 $L$ 的对角元素为1，算法步骤如下：

算法流程：

  1. 初始化：$l_{ii} = 1$ for $i = 1, 2, \ldots, n$
  2. 对于 $k = 1, 2, \ldots, n$：
    • 计算 $U$ 的第 $k$ 行：$u_{kj} = a_{kj} - \sum_{s=1}^{k-1} l_{ks}u_{sj}$ for $j = k, k+1, \ldots, n$
    • 计算 $L$ 的第 $k$ 列：$l_{ik} = \frac{a_{ik} - \sum_{s=1}^{k-1} l_{is}u_{sk}}{u_{kk}}$ for $i = k+1, k+2, \ldots, n$

伪代码：

1
2
3
4
5
6
7
8
for k = 1 to n do
    for j = k to n do
        u[k][j] = a[k][j] - sum(l[k][s] * u[s][j] for s = 1 to k-1)
    end for
    for i = k+1 to n do
        l[i][k] = (a[i][k] - sum(l[i][s] * u[s][k] for s = 1 to k-1)) / u[k][k]
    end for
end for

1.5 时间复杂度

LU分解的时间复杂度为 $O(n^3)$。具体分析：

  • 外层循环执行 $n$ 次
  • 对于第 $k$ 步，计算 $U$ 的第 $k$ 行需要 $O(n-k+1) \times O(k) = O(k(n-k+1))$ 次运算
  • 计算 $L$ 的第 $k$ 列需要 $O(n-k) \times O(k) = O(k(n-k))$ 次运算
  • 总运算次数：$\sum_{k=1}^{n} O(k(n-k+1)) + O(k(n-k)) = O(n^3)$

精确地说，需要约 $\frac{2n^3}{3}$ 次浮点运算。

1.6 实际例子

例1：对矩阵 $A = \begin{pmatrix} 2 & 1 & 3 \\ 4 & 3 & 1 \\ 6 & 5 & 2 \end{pmatrix}$ 进行LU分解。

解：使用Doolittle分解（$L$ 的对角元素为1）

第1步 ($k=1$)：

  • 计算 $U$ 的第1行：
    • $u_{11} = a_{11} = 2$
    • $u_{12} = a_{12} = 1$
    • $u_{13} = a_{13} = 3$
  • 计算 $L$ 的第1列：
    • $l_{21} = \frac{a_{21}}{u_{11}} = \frac{4}{2} = 2$
    • $l_{31} = \frac{a_{31}}{u_{11}} = \frac{6}{2} = 3$

第2步 ($k=2$)：

  • 计算 $U$ 的第2行：
    • $u_{22} = a_{22} - l_{21}u_{12} = 3 - 2 \times 1 = 1$
    • $u_{23} = a_{23} - l_{21}u_{13} = 1 - 2 \times 3 = -5$
  • 计算 $L$ 的第2列：
    • $l_{32} = \frac{a_{32} - l_{31}u_{12}}{u_{22}} = \frac{5 - 3 \times 1}{1} = 2$

第3步 ($k=3$)：

  • 计算 $U$ 的第3行：
    • $u_{33} = a_{33} - l_{31}u_{13} - l_{32}u_{23} = 2 - 3 \times 3 - 2 \times (-5) = 1$

结果：

$$L = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 2 & 1 \end{pmatrix}, \quad U = \begin{pmatrix} 2 & 1 & 3 \\ 0 & 1 & -5 \\ 0 & 0 & 1 \end{pmatrix} $$

验证：$LU = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 2 & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & 3 \\ 0 & 1 & -5 \\ 0 & 0 & 1 \end{pmatrix} = \begin{pmatrix} 2 & 1 & 3 \\ 4 & 3 & 1 \\ 6 & 5 & 2 \end{pmatrix} = A$ ✓

1.7 LU分解求解线性方程组

对于线性方程组 $Ax = b$，如果 $A = LU$，则：

$$LUx = b $$

可以分两步求解：

  1. 前向替换：求解 $Ly = b$
  2. 后向替换：求解 $Ux = y$

例2：利用上例的LU分解求解 $Ax = \begin{pmatrix} 10 \\ 14 \\ 20 \end{pmatrix}$

第1步：前向替换求解 $Ly = b$

$$\begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 2 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} 10 \\ 14 \\ 20 \end{pmatrix} $$

  • $y_1 = 10$
  • $y_2 = 14 - 2y_1 = 14 - 20 = -6$
  • $y_3 = 20 - 3y_1 - 2y_2 = 20 - 30 + 12 = 2$

第2步：后向替换求解 $Ux = y$

$$\begin{pmatrix} 2 & 1 & 3 \\ 0 & 1 & -5 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 10 \\ -6 \\ 2 \end{pmatrix} $$

  • $x_3 = 2$
  • $x_2 = -6 + 5x_3 = -6 + 10 = 4$
  • $x_1 = \frac{10 - x_2 - 3x_3}{2} = \frac{10 - 4 - 6}{2} = 0$

解：$x = \begin{pmatrix} 0 \\ 4 \\ 2 \end{pmatrix}$

1.8 选主元LU分解

当矩阵的某个主元素很小时，会导致数值不稳定。选主元LU分解通过行交换来选择较大的主元素：

$$PA = LU $$

其中 $P$ 是置换矩阵。

2. Cholesky分解

2.1 定义

Cholesky分解是针对正定矩阵的一种特殊LU分解。对于对称正定矩阵 $A$，存在唯一的下三角矩阵 $L$（对角元素为正），使得：

$$A = LL^T $$

其中 $L$ 称为Cholesky因子。

2.2 存在性定理

定理：设 $A$ 是 $n \times n$ 实对称矩阵，则 $A$ 可以进行Cholesky分解当且仅当 $A$ 是正定矩阵。

证明要点：

  • 必要性：如果 $A = LL^T$，则对任意非零向量 $x$，有 $x^TAx = x^TLL^Tx = ||L^Tx||^2 > 0$
  • 充分性：可以通过构造性证明，利用正定矩阵的性质逐步构造 $L$

2.3 Cholesky分解算法

算法步骤：
对于 $k = 1, 2, \ldots, n$：

  1. 计算对角元素：$l_{kk} = \sqrt{a_{kk} - \sum_{j=1}^{k-1} l_{kj}^2}$
  2. 计算第 $k$ 列的其他元素：$l_{ik} = \frac{a_{ik} - \sum_{j=1}^{k-1} l_{ij}l_{kj}}{l_{kk}}$ for $i = k+1, \ldots, n$

伪代码：

1
2
3
4
5
6
for k = 1 to n do
    l[k][k] = sqrt(a[k][k] - sum(l[k][j]^2 for j = 1 to k-1))
    for i = k+1 to n do
        l[i][k] = (a[i][k] - sum(l[i][j] * l[k][j] for j = 1 to k-1)) / l[k][k]
    end for
end for

2.4 时间复杂度

Cholesky分解的时间复杂度为 $O(n^3)$，但常数因子比LU分解小约一半。具体需要约 $\frac{n^3}{3}$ 次浮点运算和 $n$ 次开方运算。

由于利用了矩阵的对称性，Cholesky分解比一般的LU分解更高效：

  • 存储空间：只需存储下三角部分
  • 计算量：约为LU分解的一半
  • 数值稳定性：更好，不需要选主元

2.5 实际例子

例3：对矩阵 $A = \begin{pmatrix} 4 & 2 & 1 \\ 2 & 5 & 3 \\ 1 & 3 & 6 \end{pmatrix}$ 进行Cholesky分解。

首先验证 $A$ 是正定矩阵（所有顺序主子式为正）：

  • $M_1 = 4 > 0$
  • $M_2 = \det\begin{pmatrix} 4 & 2 \\ 2 & 5 \end{pmatrix} = 20 - 4 = 16 > 0$
  • $M_3 = \det(A) = 4(5 \times 6 - 3 \times 3) - 2(2 \times 6 - 1 \times 3) + 1(2 \times 3 - 1 \times 5) = 4 \times 21 - 2 \times 9 + 1 \times 1 = 67 > 0$

第1步 ($k=1$)：

  • $l_{11} = \sqrt{a_{11}} = \sqrt{4} = 2$
  • $l_{21} = \frac{a_{21}}{l_{11}} = \frac{2}{2} = 1$
  • $l_{31} = \frac{a_{31}}{l_{11}} = \frac{1}{2}$

第2步 ($k=2$)：

  • $l_{22} = \sqrt{a_{22} - l_{21}^2} = \sqrt{5 - 1^2} = \sqrt{4} = 2$
  • $l_{32} = \frac{a_{32} - l_{31}l_{21}}{l_{22}} = \frac{3 - \frac{1}{2} \times 1}{2} = \frac{5/2}{2} = \frac{5}{4}$

第3步 ($k=3$)：

  • $l_{33} = \sqrt{a_{33} - l_{31}^2 - l_{32}^2} = \sqrt{6 - \frac{1}{4} - \frac{25}{16}} = \sqrt{6 - \frac{4 + 25}{16}} = \sqrt{\frac{96 - 29}{16}} = \sqrt{\frac{67}{16}} = \frac{\sqrt{67}}{4}$

结果：

$$L = \begin{pmatrix} 2 & 0 & 0 \\ 1 & 2 & 0 \\ \frac{1}{2} & \frac{5}{4} & \frac{\sqrt{67}}{4} \end{pmatrix} $$

验证：

$$LL^T = \begin{pmatrix} 2 & 0 & 0 \\ 1 & 2 & 0 \\ \frac{1}{2} & \frac{5}{4} & \frac{\sqrt{67}}{4} \end{pmatrix} \begin{pmatrix} 2 & 1 & \frac{1}{2} \\ 0 & 2 & \frac{5}{4} \\ 0 & 0 & \frac{\sqrt{67}}{4} \end{pmatrix} = \begin{pmatrix} 4 & 2 & 1 \\ 2 & 5 & 3 \\ 1 & 3 & 6 \end{pmatrix} = A$$ ✓ ### 2.6 Cholesky分解求解线性方程组 对于对称正定线性方程组 $Ax = b$，如果 $A = LL^T$，则： $$LL^Tx = b$$

求解步骤：

  1. 前向替换：求解 $Ly = b$
  2. 后向替换：求解 $L^Tx = y$

例4：利用上例的Cholesky分解求解 $Ax = \begin{pmatrix} 9 \\ 13 \\ 16 \end{pmatrix}$

第1步：前向替换求解 $Ly = b$

  • $y_1 = \frac{9}{2}$
  • $y_2 = \frac{13 - 1 \times \frac{9}{2}}{2} = \frac{13 - 4.5}{2} = \frac{17}{4}$
  • $y_3 = \frac{16 - \frac{1}{2} \times \frac{9}{2} - \frac{5}{4} \times \frac{17}{4}}{\frac{\sqrt{67}}{4}} = \frac{\sqrt{67}}{4}$

第2步：后向替换求解 $L^Tx = y$（具体计算略）

3. 应用与比较

3.1 应用场景

LU分解的应用：

  • 求解一般线性方程组
  • 计算矩阵的行列式：$\det(A) = \det(L)\det(U) = \prod_{i=1}^n u_{ii}$
  • 计算矩阵的逆
  • 求解多个右端向量的线性方程组

Cholesky分解的应用：

  • 求解对称正定线性方程组
  • 最小二乘问题：$A^TAx = A^Tb$，其中 $A^TA$ 是对称正定的
  • 蒙特卡罗模拟中生成多元正态分布随机数
  • 优化问题中的海塞矩阵分解

3.2 算法比较

特性 LU分解 Cholesky分解
适用矩阵 一般方阵 对称正定矩阵
时间复杂度 $O(n^3)$ ($\frac{2n^3}{3}$ 次运算) $O(n^3)$ ($\frac{n^3}{3}$ 次运算)
存储需求 $n^2$ $\frac{n(n+1)}{2}$
数值稳定性 可能需要选主元 天然稳定
唯一性 需要额外约束 唯一（对角元为正）

3.3 数值稳定性考虑

LU分解的数值问题：

  • 当主元素很小时，会放大舍入误差
  • 解决方案：部分选主元或完全选主元

Cholesky分解的优势：

  • 正定矩阵保证了所有主元为正
  • 不需要选主元策略
  • 数值稳定性好

3.4 实现优化

内存访问优化：

  • 分块算法：减少Cache miss
  • 原地分解：节省存储空间

并行化：

  • LU分解可以进行列并行化
  • Cholesky分解可以利用对称性进行并行化

4. 总结

  1. LU分解是矩阵分解的基础，适用于一般方阵，通过将矩阵分解为下三角矩阵和上三角矩阵的乘积来简化线性方程组的求解。

  2. Cholesky分解是LU分解在对称正定矩阵上的特化，具有更好的数值稳定性和计算效率。

  3. 两种分解都有 $O(n^3)$ 的时间复杂度，但Cholesky分解的常数因子更小。

  4. 在实际应用中，应根据矩阵的性质选择合适的分解方法：

    • 一般矩阵选择LU分解（可能需要选主元）
    • 对称正定矩阵优先选择Cholesky分解

  5. 这些分解方法是数值线性代数的基础，在科学计算、工程应用和机器学习中有广泛的应用。

通过掌握这两种矩阵分解方法，我们可以更有效地解决各种线性代数问题，为进一步学习更高级的数值方法奠定坚实的基础。