手把手教您用Cholesky分解求解线性方程组：比高斯消元快一倍的方法

您正在求解一个形如 $Ax = b$ 的线性方程组。如果系数矩阵 $A$ 是对称正定矩阵，有一种名为 Cholesky分解 的特殊技巧，能让计算速度直接翻倍，同时保持极高的数值稳定性。本指南将带您从零掌握这个方法。

1. 识别问题与准备工具

确保 您的系数矩阵满足两个关键条件。对称正定矩阵是 Cholesky分解 的唯一适用对象。

• 对称性：矩阵的转置等于自身，即 $A = A^T$。例如，矩阵 [[4, 1], [1, 3]] 是对称的。
• 正定性：对于任意非零向量 $x$，满足 $x^T A x > 0$。这等价于矩阵的所有特征值大于零。

检查 矩阵是否正定。常用方法是尝试进行一次Cholesky分解，如果中途遇到负数开平方（即负的平方根），则矩阵不是正定的。许多数值计算库（如 numpy.linalg.cholesky）会在此情况下报错。

2. 核心原理：将矩阵分解为“下三角 × 其转置”

Cholesky分解的核心公式是：

$$A = L \cdot L^T$$

其中，$L$ 是一个下三角矩阵（对角线以上的元素全部为零）。求解 $Ax = b$ 的过程被巧妙地分解为两个连续的三步，每一步的计算量都远小于原始的高斯消元。

2.1 为什么这比其他方法快？

• 计算量减半：普通高斯消元法求解一个 $n \times n$ 的方程组大约需要 $\frac{2}{3}n^3$ 次浮点运算。而Cholesky分解仅需要 $\frac{1}{3}n^3$ 次运算，正好快了一倍。
• 内存占用更低：因为 $A$ 是对称的，您只需要存储矩阵的下三角部分（包括对角线）。分解得到的 $L$ 可以直接覆盖存储 $A$ 的数组，无需额外内存。
• 数值稳定性极高：分解过程中不会出现像高斯消元那样需要通过“部分主元”来避免除零的麻烦。只要矩阵是正定的，分解过程就天然稳定，不会产生巨大误差。

3. 操作步骤：用Cholesky分解求解方程组

假设您有一个正定对称矩阵 $A$ 和方程组 $Ax = b$。我们通过手算示例和Python代码两个视角演示。

3.1 手算一个3x3的例子

给定 矩阵 $A$ 和向量 $b$：

$$A = \begin{bmatrix} 4 & 2 & 1 \\ 2 & 5 & 3 \\ 1 & 3 & 6 \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}$$

步骤 1：计算Cholesky分解 $A = L L^T$。

• 计算 $L_{11} = \sqrt{A_{11}} = \sqrt{4} = 2$
• 计算 $L_{21} = A_{21} / L_{11} = 2 / 2 = 1$
• 计算 $L_{31} = A_{31} / L_{11} = 1 / 2 = 0.5$
• 计算 $L_{22} = \sqrt{A_{22} - L_{21}^2} = \sqrt{5 - 1^2} = \sqrt{4} = 2$
• 计算 $L_{32} = (A_{32} - L_{21} \cdot L_{31}) / L_{22} = (3 - 1 \cdot 0.5) / 2 = 2.5 / 2 = 1.25$
• 计算 $L_{33} = \sqrt{A_{33} - (L_{31}^2 + L_{32}^2)} = \sqrt{6 - (0.5^2 + 1.25^2)} = \sqrt{6 - (0.25 + 1.5625)} = \sqrt{4.1875} \approx 2.046$

得到 下三角矩阵 $L$：

$$L = \begin{bmatrix} 2 & 0 & 0 \\ 1 & 2 & 0 \\ 0.5 & 1.25 & 2.046 \end{bmatrix}$$

步骤 2：解下三角方程组 $Ly = b$。

• 从第一行开始：$2y_1 = 1 \Rightarrow y_1 = 0.5$
• 第二行：$1 \cdot 0.5 + 2y_2 = 2 \Rightarrow y_2 = (2 - 0.5) / 2 = 0.75$
• 第三行：$0.5 \cdot 0.5 + 1.25 \cdot 0.75 + 2.046 y_3 = 3 \Rightarrow y_3 = (3 - 0.25 - 0.9375) / 2.046 \approx 0.886$

得到 中间向量 $y = [0.5, 0.75, 0.886]^T$

步骤 3：解上三角方程组 $L^T x = y$。

• 从第三行开始（因为 $L^T$ 的上三角部分）：$2.046 x_3 = 0.886 \Rightarrow x_3 \approx 0.433$
• 第二行：$2x_2 + 1.25 \cdot 0.433 = 0.75 \Rightarrow x_2 = (0.75 - 0.541) / 2 = 0.1045$
• 第一行：$2x_1 + 1 \cdot 0.1045 + 0.5 \cdot 0.433 = 0.5 \Rightarrow x_1 = (0.5 - 0.1045 - 0.2165) / 2 = 0.0895$

最终得到 解向量 $x \approx [0.0895, 0.1045, 0.433]^T$。您可以用原方程 $Ax$ 验证是否等于 $b$，误差应在可接受范围内。

3.2 用Python代码一步到位

如果手算繁琐，您可以直接用 numpy 库完成全部工作。请确保已安装 numpy。

import numpy as np

# 定义矩阵A和向量b
A = np.array([[4, 2, 1],
              [2, 5, 3],
              [1, 3, 6]], dtype=float)

b = np.array([1, 2, 3], dtype=float)

# 步骤1：进行Cholesky分解
L = np.linalg.cholesky(A)
# 验证：A = L @ L.T 应近似成立
print("L矩阵:\n", L)
print("重构的A = L @ L.T:\n", L @ L.T)

# 步骤2：解 Ly = b（下三角）
y = np.linalg.solve(L, b)

# 步骤3：解 L^T x = y（上三角）
x = np.linalg.solve(L.T, y)

print("解向量 x:\n", x)

执行 上述代码，您会得到与手算近似的解。如果您想一步完成所有求解，可以直接使用 scipy.linalg.cho_solve，它专为Cholesky分解后的方程组优化。

from scipy.linalg import cho_factor, cho_solve

# 进行Cholesky分解并返回一个紧凑的结构
L_and_lower = cho_factor(A)

# 直接求解 x = A^{-1} b
x = cho_solve(L_and_lower, b)
print("用 scipy 一步求解的 x:\n", x)

4. 应用场景与性能验证

Cholesky分解在实际工程中广泛应用，尤其是在大型稀疏矩阵（如有限元分析、图论中的拉普拉斯矩阵）和统计建模（如多元正态分布的协方差矩阵）中。它的数值优势在以下场景体现得淋漓尽致：

• 科学计算：求解由偏微分方程离散化产生的线性系统。
• 金融模型：计算投资组合的VaR（风险价值）时，需要反复求解协方差矩阵。
• 机器学习：线性回归中的正规方程 $X^T X \theta = X^T y$，其中 $X^T X$ 就是对称正定矩阵。

验证 其性能优势：尝试用一个大矩阵（例如1000×1000）对比Cholesky分解与高斯消元的求解速度。

import time
import numpy as np

n = 1000
# 生成一个对称正定矩阵
A = np.random.rand(n, n)
A = A @ A.T + np.eye(n) * n  # 保证正定
b = np.random.rand(n)

# 使用Cholesky分解
start = time.time()
L = np.linalg.cholesky(A)
y = np.linalg.solve(L, b)
x_chol = np.linalg.solve(L.T, y)
time_chol = time.time() - start

# 使用普通高斯消元（numpy自动选择）
start = time.time()
x_gauss = np.linalg.solve(A, b)
time_gauss = time.time() - start

print(f"Cholesky耗时: {time_chol:.4f} 秒")
print(f"高斯消元耗时: {time_gauss:.4f} 秒")
print(f"速度提升比例: {time_gauss / time_chol:.2f} 倍")

运行 这段代码，您会看到Cholesky方法通常比默认的高斯消元快1.5到2倍，且矩阵越大优势越明显。核心原因是它利用对称性减少了约一半的运算，并避免了主元选取的额外开销。

5. 处理边界情况：当矩阵不是正定时怎么办？

如果您的矩阵是对称的但不是正定的，Cholesky分解会失败（出现负的平方根）。此时有不同的处理策略：

• 改用LDL分解：如果矩阵是对称不定的（即某些特征值为负），使用 $A = L \cdot D \cdot L^T$ 的分解，其中 $D$ 是对角矩阵。在 scipy.linalg 中可以用 ldl 函数。
• 采用主元策略：如果问题来自接触力学或约束优化，可以使用部分主元Cholesky，或降级使用LU分解（scipy.linalg.lu）。
• 近似正定化：如果矩阵在数值上接近正定（比如包含微小噪音），可以尝试做一个小量的正则化：$A' = A + \epsilon I$，其中 $\epsilon$ 取一个很小的正数（如 1e-8）。这在统计中称为“岭回归”技巧。