Cholesky分解

Cholesky 分解

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

在求解线性方程组、矩阵求逆等核心线性代数运算中，算法的效率至关重要。当矩阵具有特殊的“正定性”时，我们可以利用其优良的数学性质，采用一种比通用LU分解快约一倍的分解方法—Cholesky分解。本节将深入讲解Cholesky分解的定义、唯一性、构造性证明、高效算法及其在检验矩阵正定性、求解线性方程组和矩阵求逆中的关键应用，帮助你掌握这一处理正定矩阵的利器。

💡 核心导读

本节你将了解到：

1. 定义与唯一性：Cholesky分解是正定矩阵的“平方根”分解，形式为 $A = LL^T$，其中 $L$ 是对角元为正的下三角矩阵，且分解唯一。
2. 构造与算法：通过数学归纳法和对称的初等行列变换，可以构造性地证明分解的存在性，并推导出时间复杂度为 $O(\frac{1}{3}n^3)$ 的实用算法。
3. 核心应用一：检验正定性：尝试对矩阵进行Cholesky分解是检验其是否正定的最有效数值方法之一，算法失败即表明矩阵非正定。
4. 核心应用二：高效解方程：对于正定线性方程组 $Ax=b$，利用Cholesky分解 $A=LL^T$ 并求解两个三角方程组，其复杂度优于通用LU分解。
5. 核心应用三：稳定求逆：结合Cholesky分解和下三角矩阵求逆的平方复杂度算法，可以稳定、高效地计算正定矩阵的逆，但仍需注意求逆本身的高计算成本。

1. 定义、性质与直观理解

Cholesky分解（Cholesky decomposition）是专门针对正定矩阵的一种特殊矩阵分解。

定义（Cholesky分解）：设 $A \in \mathbb{R}^{n \times n}$ 是一个对称正定矩阵。则存在一个唯一的对角元全为正实数的下三角矩阵 $L \in \mathbb{R}^{n \times n}$，使得 $$A = L L^T.$$ 我们称 $L$ 为矩阵 $A$ 的 Cholesky因子。

直观理解：

• 矩阵的“平方根”：正如一个正数 $a$ 可以写成 $a = (\sqrt{a})^2$，正定矩阵 $A$ 可以写成 $A = L L^T$。因此，$L$ 可以视为 $A$ 的某种“平方根”。
• 唯一性：对于正定矩阵 $A$，满足上述条件的下三角矩阵 $L$ 是唯一的。如果 $A$ 是半正定矩阵，则分解存在但不唯一。
• 不考虑非正定情况：对于非正定（不定或负定）矩阵，我们不讨论其Cholesky分解。

Cholesky分解主要有两个实用出发点：

1. 已知正定，进行计算：当我们已知矩阵 $A$ 是正定的，并需要对其进行某些操作（如解方程、求逆）时，先进行Cholesky分解往往是高效的第一步。
2. 检验是否正定：当我们不确定一个对称矩阵是否正定时，可以尝试对其执行Cholesky分解算法。如果算法顺利完成，则 $A$ 是正定的；如果算法在中途因出现零或负数（如试图对非正数开平方）而失败，则 $A$ 不是正定矩阵。

2. 存在性证明与计算思想

Cholesky分解的存在性可以通过数学归纳法构造性地证明，这个证明过程也揭示了其核心计算思想。

证明思路（数学归纳法）：

1. 归纳基础：当 $n=1$ 时，$A = [a_{11}]$ 是一个正数。取 $L = [\sqrt{a_{11}}]$，显然有 $A = L L^T$，结论成立。
2. 归纳假设：假设对于所有 $n-1$ 阶正定矩阵，Cholesky分解存在。
3. 归纳步骤：考虑 $n$ 阶正定矩阵 $A$。因其正定且对称，可进行如下分块： $$A = \begin{pmatrix} a_{11} & \mathbf{w}^T \ \mathbf{w} & K \end{pmatrix}$$ 其中 $a_{11} > 0$，$\mathbf{w}$ 是 $n-1$ 维列向量，$K$ 是 $n-1$ 阶矩阵。
4. 消去变换：我们的目标是通过对称的初等行列变换，消去分块中的非对角块 $\mathbf{w}$ 和 $\mathbf{w}^T$。构造下三角矩阵： $$L_1 = \begin{pmatrix} 1 & \mathbf{0}^T \ -\frac{1}{a_{11}}\mathbf{w} & I_{n-1} \end{pmatrix}$$ 对 $A$ 进行合同变换： $$L_1 A L_1^T = \begin{pmatrix} 1 & \mathbf{0}^T \ -\frac{1}{a_{11}}\mathbf{w} & I_{n-1} \end{pmatrix} \begin{pmatrix} a_{11} & \mathbf{w}^T \ \mathbf{w} & K \end{pmatrix} \begin{pmatrix} 1 & -\frac{1}{a_{11}}\mathbf{w}^T \ \mathbf{0} & I_{n-1} \end{pmatrix} = \begin{pmatrix} a_{11} & \mathbf{0}^T \ \mathbf{0} & K - \frac{1}{a_{11}}\mathbf{w}\mathbf{w}^T \end{pmatrix}.$$ 记 $A_2 = K - \frac{1}{a_{11}}\mathbf{w}\mathbf{w}^T$。
5. 子矩阵的正定性：可以证明，$A_2$ 仍然是正定矩阵。因为合同变换不改变矩阵的正定性，且 $L_1 A L_1^T$ 的左上角 $a_{11}>0$，根据正定矩阵的判定定理（各阶顺序主子式大于零），其右下角的 $n-1$ 阶子矩阵 $A_2$ 也必为正定。
6. 递归分解：由归纳假设，$n-1$ 阶正定矩阵 $A_2$ 存在Cholesky分解：$A_2 = L_2 L_2^T$，其中 $L_2$ 是下三角矩阵。
7. 组合因子：至此我们已有 $$L_1 A L_1^T = \begin{pmatrix} a_{11} & \mathbf{0}^T \ \mathbf{0} & A_2 \end{pmatrix} = \begin{pmatrix} \sqrt{a_{11}} & \mathbf{0}^T \ \mathbf{0} & L_2 \end{pmatrix} \begin{pmatrix} \sqrt{a_{11}} & \mathbf{0}^T \ \mathbf{0} & L_2 \end{pmatrix}^T.$$ 对等式两边左乘 $L_1^{-1}$、右乘 $(L_1^T)^{-1}$，由 $L_1^{-1} = \begin{pmatrix} 1 & \mathbf{0}^T \ \frac{1}{a_{11}}\mathbf{w} & I_{n-1} \end{pmatrix}$，可得 $$A = \left(L_1^{-1} \begin{pmatrix} \sqrt{a_{11}} & \mathbf{0}^T \ \mathbf{0} & L_2 \end{pmatrix}\right) \left(L_1^{-1} \begin{pmatrix} \sqrt{a_{11}} & \mathbf{0}^T \ \mathbf{0} & L_2 \end{pmatrix}\right)^T.$$ 因此 $$L = L_1^{-1} \begin{pmatrix} \sqrt{a_{11}} & \mathbf{0}^T \ \mathbf{0} & L_2 \end{pmatrix} = \begin{pmatrix} 1 & \mathbf{0}^T \ \frac{1}{a_{11}}\mathbf{w} & I_{n-1} \end{pmatrix} \begin{pmatrix} \sqrt{a_{11}} & \mathbf{0}^T \ \mathbf{0} & L_2 \end{pmatrix} = \begin{pmatrix} \sqrt{a_{11}} & \mathbf{0}^T \ \frac{1}{\sqrt{a_{11}}}\mathbf{w} & L_2 \end{pmatrix}.$$ 它是对角元为正的下三角矩阵（$\sqrt{a_{11}} > 0$，$L_2$ 的对角元由归纳假设为正），且满足 $A = L L^T$，归纳完成。对于 $n$ 阶矩阵，上述“消元—分解”过程可递归执行，最终 $L$ 等于所有初等下三角矩阵逆的乘积，仍为下三角矩阵。

上述证明过程本质上是高斯消元法在正定矩阵情形下的一个变体，但利用了对称性，只需对矩阵的一半元素进行操作。

算例：以矩阵 $A = \begin{pmatrix}4 & 2 \ 2 & 5\end{pmatrix}$ 为例。

1. 第一步：消去 $a_{21}=2$。变换矩阵 $L_1 = \begin{pmatrix}1 & 0 \ -\frac{2}{4} & 1\end{pmatrix} = \begin{pmatrix}1 & 0 \ -0.5 & 1\end{pmatrix}$。 $$L_1 A L_1^T = \begin{pmatrix}1 & 0 \ -0.5 & 1\end{pmatrix} \begin{pmatrix}4 & 2 \ 2 & 5\end{pmatrix} \begin{pmatrix}1 & -0.5 \ 0 & 1\end{pmatrix} = \begin{pmatrix}4 & 0 \ 0 & 4\end{pmatrix}.$$ 这里 $A_2 = [4]$。
2. 对 $A_2$ 分解：$4 = 2 \times 2$。
3. 组合：根据构造，$L = \begin{pmatrix}\sqrt{4} & 0 \ \frac{2}{\sqrt{4}} & 2\end{pmatrix} = \begin{pmatrix}2 & 0 \ 1 & 2\end{pmatrix}$。 验证：$L L^T = \begin{pmatrix}2 & 0 \ 1 & 2\end{pmatrix} \begin{pmatrix}2 & 1 \ 0 & 2\end{pmatrix} = \begin{pmatrix}4 & 2 \ 2 & 5\end{pmatrix} = A$。

3. 算法、复杂度与算例

基于上述构造思想，可以得到Cholesky分解的数值算法（Cholesky–Banachiewicz算法）。

算法描述： 输入：$n$ 阶对称正定矩阵 $A$。 输出：下三角矩阵 $L$，使得 $A = L L^T$。

1. 初始化 $L$ 为零矩阵。
2. For $j = 1$ to $n$: a. 计算 $L_{jj} = \sqrt{A_{jj} - \sum_{k=1}^{j-1} L_{jk}^2}$。 b. For $i = j+1$ to $n$: $L_{ij} = \frac{1}{L_{jj}} (A_{ij} - \sum_{k=1}^{j-1} L_{ik} L_{jk})$。
3. 返回 $L$。

算法复杂度分析：该算法的主要计算量集中在三重循环的内层求和上。总运算量约为 $\frac{1}{3}n^3$ 次浮点运算（FLOPs）。作为对比，通用的LU分解需要约 $\frac{2}{3}n^3$ FLOPs。因此，Cholesky分解的计算量约为LU分解的一半。

稳定性：由于正定矩阵的顺序主子式均大于零，算法中 $L_{jj}$ 的计算不会出现对零或负数开根号的情况，因此不需要选主元（pivoting）。这使得Cholesky分解在数值上非常稳定。

下面我们用NumPy来验证Cholesky分解及其应用。

import numpy as np

# 定义一个对称正定矩阵 A
A = np.array([[4.0, 2.0, 1.0],
              [2.0, 5.0, 3.0],
              [1.0, 3.0, 6.0]])

print("原矩阵 A:")
print(A)

# 应用1: 使用numpy的Cholesky分解，返回下三角矩阵 L，满足 A = L L^T
L = np.linalg.cholesky(A)  # 下三角Cholesky因子
L_T = L.T  # 上三角因子 L^T
print("\nCholesky因子 L (下三角):")
print(L)
print("\n验证 A = L L^T:")
print(L @ L_T)

# 应用2: 解线性方程组 Ax = b
b = np.array([1.0, 3.0, 2.0])
print(f"\n右端项 b: {b}")
# 方法：先解下三角方程组 Ly = b (前向代入)
# 注：此处用通用solve演示，生产环境可使用scipy.linalg.solve_triangular(L, b, lower=True)优化
y = np.linalg.solve(L, b)
# 再解上三角方程组 L^T x = y (后向代入)
x = np.linalg.solve(L_T, y)
print(f"解线性方程组 Ax=b 的解 x: {x}")
print(f"验证 Ax 是否等于 b: {np.allclose(A @ x, b)}")

# 应用3: 利用分解求逆 A^{-1} = (L^{-1})^T (L^{-1})
# 首先求下三角矩阵L的逆。下三角矩阵的逆也是下三角矩阵，可以通过前向代入法高效求解。
# 这里我们使用一个简单实现（对于教学演示），实际中可使用scipy.linalg.solve_triangular。
def invert_lower_triangular(L):
    """求下三角矩阵L的逆。假设L非奇异且为下三角。"""
    n = L.shape[0]
    invL = np.zeros_like(L)
    for i in range(n):
        invL[i, i] = 1.0 / L[i, i]  # 对角线元素
        for j in range(i):
            s = 0.0
            for k in range(j, i+1): # 注意求和范围
                # 更精确的写法是：对于i>j， sum(L[i,k]*invL[k,j]) for k=j to i-1
                if k < i: # 当k=i时，是L[i,i]*invL[i,j]，此项是待求的，移到等式左边
                    s += L[i, k] * invL[k, j]
            # L[i,i]*invL[i,j] + s = 0 => invL[i,j] = -s / L[i,i]
            invL[i, j] = -s / L[i, i]
    return invL

invL = invert_lower_triangular(L)
invA_via_cholesky = invL.T @ invL
print("\n通过Cholesky分解求得的逆矩阵 A_inv:")
print(invA_via_cholesky)

# 与numpy直接求逆的结果对比
invA_np = np.linalg.inv(A)
print("\nNumPy直接求逆的结果:")
print(invA_np)
print(f"\n两种方法结果是否接近: {np.allclose(invA_via_cholesky, invA_np)}")

4. 三大核心应用场景

4.1 检验矩阵的正定性

这是Cholesky分解一个非常独特且重要的应用。对于一个对称矩阵 $A$，要判断其是否正定，最直接的数值方法就是尝试进行Cholesky分解。

• 如果分解成功完成，则 $A$ 是正定矩阵，并且同时得到了其Cholesky因子 $L$。
• 如果分解过程中失败（例如，在计算 $L_{jj}$ 时遇到 $A_{jj} - \sum_{k=1}^{j-1} L_{jk}^2 \le 0$），则 $A$ 不是正定矩阵（可能是半正定或不定）。 这种方法的时间复杂度为 $O(\frac{1}{3}n^3)$，是检验正定性的高效数值手段。

4.2 求解正定线性方程组

对于线性方程组 $A\mathbf{x} = \mathbf{b}$，其中 $A$ 是 $n \times n$ 正定矩阵。

1. 首先计算Cholesky分解：$A = L L^T$。
2. 原方程化为 $L L^T \mathbf{x} = \mathbf{b}$。
3. 令 $\mathbf{y} = L^T \mathbf{x}$，则先解 下三角方程组 $L \mathbf{y} = \mathbf{b}$（前向代入，复杂度 $O(n^2)$）。
4. 再解 上三角方程组 $L^T \mathbf{x} = \mathbf{y}$（后向代入，复杂度 $O(n^2)$）。

总复杂度主要由Cholesky分解的 $O(\frac{1}{3}n^3)$ 主导，优于使用LU分解的 $O(\frac{2}{3}n^3)$。在科学计算中，对于大规模正定系统（如来自有限元法、优化问题），Cholesky分解是首选方法。

4.3 计算正定矩阵的逆

虽然直接求逆通常应避免，但在某些情况下不可避免。对于正定矩阵 $A$，利用Cholesky分解可以稳定地求逆： $$A^{-1} = (L L^T)^{-1} = (L^T)^{-1} L^{-1} = (L^{-1})^T (L^{-1}).$$ 关键在于求下三角矩阵 $L$ 的逆 $L^{-1}$。由于 $L$ 是下三角矩阵，其逆也是下三角矩阵，并且可以通过简单的前向代入法以 $O(n^2)$ 的复杂度求得（如前面代码示例所示）。因此，总的求逆复杂度仍由Cholesky分解的 $O(\frac{1}{3}n^3)$ 决定。

重要提醒：尽管这种方法比通用的高斯消元法求逆（$O(n^3)$）在常数因子上有优势，但矩阵求逆本身始终是 $O(n^3)$ 复杂度的操作，计算成本很高。在机器学习和优化中，许多涉及逆矩阵的公式（如正态方程 $(X^TX)^{-1}X^Ty$）都可以通过Cholesky分解求解线性方程组来规避显式求逆，这是更优的实践。

📝 动手练一练

1. 手动计算：对以下正定矩阵 $A = \begin{pmatrix} 9 & 3 \ 3 & 5 \end{pmatrix}$ 进行Cholesky分解，即求下三角矩阵 $L$ 使得 $A = L L^T$。
2. 编程验证：使用NumPy的 np.linalg.cholesky 函数验证你手算的 $L$ 是否正确，并利用该分解求解方程组 $A\mathbf{x} = \begin{pmatrix} 12 \ 11 \end{pmatrix}$。

参考答案：

1. 手算过程： 设 $L = \begin{pmatrix} l_{11} & 0 \ l_{21} & l_{22} \end{pmatrix}$。 由 $A = L L^T$： $\begin{pmatrix} 9 & 3 \ 3 & 5 \end{pmatrix} = \begin{pmatrix} l_{11} & 0 \ l_{21} & l_{22} \end{pmatrix} \begin{pmatrix} l_{11} & l_{21} \ 0 & l_{22} \end{pmatrix} = \begin{pmatrix} l_{11}^2 & l_{11}l_{21} \ l_{11}l_{21} & l_{21}^2 + l_{22}^2 \end{pmatrix}$。 比较元素：

   • $l_{11}^2 = 9 \Rightarrow l_{11} = 3$（取正）。
   • $l_{11}l_{21} = 3 \Rightarrow 3 \cdot l_{21} = 3 \Rightarrow l_{21} = 1$。
   • $l_{21}^2 + l_{22}^2 = 5 \Rightarrow 1 + l_{22}^2 = 5 \Rightarrow l_{22}^2 = 4 \Rightarrow l_{22} = 2$（取正）。 因此，$L = \begin{pmatrix} 3 & 0 \ 1 & 2 \end{pmatrix}$。

2. 编程验证：

   import numpy as np
   A = np.array([[9.0, 3.0], [3.0, 5.0]])
   b = np.array([12.0, 11.0])
   
   # 使用numpy的Cholesky分解，返回下三角因子 L
   L_np = np.linalg.cholesky(A)
   L_T_np = L_np.T  # 上三角因子 L^T
   print("NumPy计算的 L:")
   print(L_np)
   
   # 验证分解
   print("\n验证 A = L L^T:")
   print(L_np @ L_T_np)
   
   # 解方程 Ly = b
   y = np.linalg.solve(L_np, b)
   # 解方程 L^T x = y
   x = np.linalg.solve(L_T_np, y)
   print(f"\n方程的解 x: {x}")
   print(f"验证 Ax = b: {np.allclose(A @ x, b)}")

本章小结

本节深入探讨了针对对称正定矩阵的Cholesky分解。

要点回顾：

• 定义：正定矩阵 $A$ 可唯一分解为 $A = LL^T$，$L$ 是对角元为正的下三角矩阵。
• 算法与复杂度：基于高斯消元思想的Cholesky算法复杂度约为 $\frac{1}{3}n^3$，且因正定性保证而无需选主元，数值稳定。
• 核心应用：
  1. 高效检验正定性：尝试分解是最实用的数值判定方法。
  2. 快速求解线性方程组：对于 $Ax=b$，分解后求解两个三角方程组，速度快于通用LU分解。
  3. 稳定求逆：结合下三角矩阵求逆，可稳定计算 $A^{-1}$，但应优先考虑避免显式求逆的方案。

行动清单：

1. 掌握手算：对于2维或3维简单正定矩阵，能够手动完成Cholesky分解。
2. 熟练调用：在Python中，学会使用 np.linalg.cholesky 进行分解，并利用其结果解方程。
3. 建立条件反射：在后续学习机器学习模型（如高斯过程、优化中的牛顿法）时，遇到对称正定矩阵，优先考虑使用Cholesky分解来提升计算效率和数值稳定性。

— 小象教研组