📚 课程目录第4章 Cholesky分解与对称正定系统

4. Cholesky分解与对称正定系统

在方程组求解的众多结构中，对称正定（Symmetric Positive Definite，SPD）矩阵占据特殊地位：它既保证了分解的存在性与唯一性，又允许我们使用比LU分解少一半计算量、且无需选主元的Cholesky分解。本节围绕SPD矩阵的性质、Cholesky分解的算法推导、平方根法的实现，以及分解存在条件与数值稳定性展开。

一、对称正定矩阵的性质

设 A ∈ R^(n×n) 满足 A = Aᵀ，且对任意非零向量 x 都有 xᵀAx > 0，则称 A 为对称正定矩阵。这一条件等价于以下几条常用判据：

特征值判据：所有特征值 λᵢ > 0；
顺序主子式判据：各阶顺序主子式 Δₖ > 0（k = 1,…,n）；
合同判据：存在非奇异矩阵 B 使 A = BᵀB。

SPD矩阵还带来若干工程上极有价值的推论：对角元恒为正（aᵢᵢ > 0）；绝对值最大的元素必出现在对角线上；矩阵非奇异且条件数 κ₂(A) = λmax/λmin ≥ 1。这些性质正是Cholesky分解能够稳定进行、且无需行交换的根本原因。

二、Cholesky分解算法

Cholesky分解断言：若 A 对称正定，则存在唯一的对角元全为正的下三角矩阵 L，使得

A = L·Lᵀ

其中 L 的对角元 lᵢᵢ > 0。将 A = LLᵀ 按元素展开，对 i ≥ j 有

aᵢⱼ = Σₖ₌₁^j lᵢₖ·lⱼₖ

由此可逐列（或逐行）递推得到计算公式。按列计算时，第 j 列的元素为：

lⱼⱼ = sqrt( aⱼⱼ − Σₖ₌₁^(j−1) lⱼₖ² )
lᵢⱼ = ( aᵢⱼ − Σₖ₌₁^(j−1) lᵢₖ·lⱼₖ ) / lⱼⱼ ,  i = j+1,…,n

该递推的关键在于根号内的量 aⱼⱼ − Σ lⱼₖ²。对SPD矩阵，可以证明它恒为正，因此开方总有实数解；一旦该量出现零或负数，就说明矩阵并非正定（或已严重病态），这是实现中最重要的“健康检查点”。

三、平方根法的实现

所谓“平方根法”，即利用 A = LLᵀ 将原方程组 Ax = b 转化为两个三角方程组：先解 Ly = b（前代），再解 Lᵀx = y（回代）。完整实现如下（Python 示例，仅示意算法逻辑）：

import math

def cholesky(A):
    n = len(A)
    L = [[0.0]*n for _ in range(n)]
    for j in range(n):
        s = A[j][j] - sum(L[j][k]**2 for k in range(j))
        if s <= 0:
            raise ValueError("矩阵非正定或已病态：主元非正")
        L[j][j] = math.sqrt(s)
        for i in range(j+1, n):
            t = A[i][j] - sum(L[i][k]*L[j][k] for k in range(j))
            L[i][j] = t / L[j][j]
    return L

def solve_cholesky(A, b):
    L = cholesky(A)
    n = len(A)
    y = [0.0]*n
    for i in range(n):                       # 前代 Ly = b
        y[i] = (b[i] - sum(L[i][k]*y[k] for k in range(i))) / L[i][i]
    x = [0.0]*n
    for i in range(n-1, -1, -1):             # 回代 Lᵀx = y
        x[i] = (y[i] - sum(L[k][i]*x[k] for k in range(i+1, n))) / L[i][i]
    return x

计算量方面，Cholesky分解约需 n³/6 次乘加与 n 次开方，约为LU分解（n³/3）的一半；存储上只需保存 L 的下三角部分，约为 n²/2 个元素。对于对称正定的大型稀疏系统（如有限元刚度矩阵），这一优势尤为明显。

四、分解存在条件与数值稳定性

Cholesky分解存在的充要条件正是矩阵对称正定。若矩阵仅半正定（存在零特征值），则某个主元会恰为零，分解在实数域内失败；若矩阵对称但不定，则中途会出现负的主元开方。因此实现中必须对根号内量做非正检测，而不能像LU那样依赖选主元“补救”——Cholesky分解本身不允许行交换，否则会破坏对称结构。

数值稳定性方面，Cholesky分解是无主元分解中稳定性最好的一种：可以证明其元素满足 |lᵢⱼ| ≤ √aⱼⱼ，即 L 的元素不会相对 A 的对角元发生爆炸式增长，前向误差可控。但“稳定”不等于“精确”：当 A 的条件数很大时，κ₂(A) = κ₂(L)²，即 L 的条件数是原矩阵条件数的平方根，这实际上比直接对 A 操作更温和；然而求解误差仍随 κ₂(A) 放大。若 A 接近奇异，主元 s 会非常小，开方后 lⱼⱼ 也很小，后续除法会放大舍入误差，此时应结合条件数估计或改用带对角修正（如 A + δI）的策略。

工程建议：对SPD系统，优先使用Cholesky而非通用LU；在分解循环中记录最小主元 min(s)，若其与 max(aᵢᵢ) 之比低于约 1e-12，应视为病态预警，提示检查矩阵构造或改用更高精度/正则化。

上图把分解、主元健康检查与两段三角求解串成一条完整链路：左侧给出 L 的下三角结构，中部菱形是唯一需要做非正判断的位置，右侧失败分支提示非正定或病态时的处理方向。实际编码时，把主元检测放在分解循环内，比事后检查 A 的特征值更经济，也更能定位问题出现的列。

最后需要强调：Cholesky分解的“稳定”是相对于无主元LU而言的，它并不能消除问题本身的病态性。当物理模型或数据导致矩阵接近奇异时，正确的做法是识别病态（观察最小主元、估计条件数），再通过正则化、缩放或改用QR/奇异值分解等手段处理，而不是指望Cholesky分解自动“救回”精度。