写在前面：从“会分解”到“会求逆”

上一篇我们详细拆解了 Cholesky 分解——用 30 行 C 代码将对称正定矩阵 A 分解为 A = R^T R，其中 R 是上三角矩阵。当时只用半篇文章的篇幅简单介绍了 Cholesky 求逆部分，留下了一个悬念：所谓“两步法求逆”究竟如何高效计算出 A^-1？为什么它比通用方法快一倍？求逆结果为何是个“半成品”？

本文就彻底讲透这个问题。我们将聚焦于 Cholesky 的后半段（即求逆部分），而非分解部分。这是教科书中最含糊其辞、而工程代码中最值得玩味的一段内容。读完你会明白：

• 为什么三角矩阵的逆只需 O(n²) 回代而非 O(n³) 消元；
• 为什么求出的逆只填充了上三角，调用方还需自己“对称化”一次；
• 为什么这套两步法的总计算量约为通用 Gauss-Jordan 的一半；
• 为什么行列式可以“顺路”从求逆循环中提取出来。

一句话概括本篇主旨：通用的求逆算法把矩阵当作“一坨数”无差别地处理；而 Cholesky 求逆则尊重矩阵的结构（对称、正定、已分解），把一个大问题拆成两个能利用结构的小问题。结构信息用得越充分，节省的计算量就越多。

钩子：教科书默认的求逆方法其实是“暴力法”

随便翻开一本数值分析或线性代数教材，“求矩阵的逆”那一节八成是这样教的：

然后给你一个 3×3 的手算例子，让你在草稿纸上写满半个黑板。考试也考这个，于是你记住了：“求逆 = 增广 + 全选主元 + 消元”。

这个方法有一个致命问题：它无视矩阵的一切结构。不管 A 是对称的、稀疏的、正定的、还是带状的，Gauss-Jordan 都把它当作 n × n 的“一坨数”，从左上角消到右下角，每个元素都参与消元。代价是实打实的 O(n³) 浮点运算，且过程中要做主元选取、行交换，数值上还得提防大主元相消。

但如果你已经知道 A 是对称正定的呢？这时候“暴力法”就是暴殄天物——你手握一份关于矩阵结构的强信息（它有 Cholesky 分解、可以分解成两个三角矩阵的乘积），却完全没用上。

Cholesky 求逆的核心思想，就是把这份结构信息转化为计算量的节省。先分解（上一篇讲过，n³/3 量级），再分别处理两个三角矩阵（这一篇要讲，各 n³/3 量级），合起来比 Gauss-Jordan 的 2n³/3 少了一截，而且更稳定。

下面正式拆解程序。

数学原理：把求逆拆成两个“小一号”的问题

核心恒等式

A^-1 = (R^T R)^-1 = R^-1 (R^-1)^T

证明只需两行：(R^T R)^-1 = R^-1 (R^T)^-1 = R^-1 (R^-1)^T，最后一步利用了“转置的逆等于逆的转置”。这条恒等式看似不起眼，却是整个两步法的根基。

它告诉我们：求 A^-1 不必碰 A，只需求 R^-1 然后做一次三角乘法即可。而 R 是上三角矩阵，三角矩阵的求逆有专门的、便宜得多的算法——回代。

三角矩阵求逆的本质：回代，每个元素 O(n)

为什么三角矩阵求逆便宜？因为三角矩阵的逆“结构上”还是三角矩阵（上三角的逆仍是上三角），而且每个元素都可以用一次回代得到。

把 R 的逆记作 S = R^-1，满足 R S = I。展开 R S = I 的第 i 行第 j 列元素（约定 i ≤ j，因为 S 也是上三角，i > j 时 S_ij = 0）：

∑_k=i^j R_ik S_kj = δ_ij

注意求和从 k = i 开始（因为 R 是上三角，R_ik = 0 当 k < i），到 k = j 结束（因为 S 是上三角，S_kj = 0 当 k > j）。这是一个只跨越 [i, j] 区间的小型求和，而不是全跨度 [1, n]。

把 k = i 那一项拎出来：

R_ii S_ij + ∑_k=i+1^j R_ik S_kj = δ_ij

解出 S_ij：

S_ij = (δ_ij - ∑_k=i+1^j R_ik S_kj) / R_ii

分两种情况看：

• 对角元（i = j）：δ_ii = 1，求和项为空集，得到 S_ii = 1 / R_ii——对角元就是直接取倒数。
• 非对角元（i < j）：δ_ij = 0，得到 S_ij = -(∑_k=i+1^j R_ik S_kj) / R_ii——一个内积再除以对角元。

注意求和 ∑_k=i+1^j 里：R_ik 是 R 第 i 行从 i+1 到 j 的元素，S_kj 是 S 第 j 列从 i+1 到 j 的元素。这两个片段在算法开始计算 S_ij 时必须已经算好。这就规定了计算顺序：对每一列 j，从下往上（i 从 j-1 递减到 0）算非对角元，最后处理对角元 S_jj。

这就是“三角矩阵求逆 = 回代”。每个 S_ij 的代价是 O(j - i)，平均 O(n)；总元素数 O(n²)，所以总代价 O(n³)——但常数因子只是通用方法的一个零头（下面会精确计算）。

合成 A^-1：三角乘三角

第二步：拿到 S = R^-1 后，计算 B = S · S^T。这也是一个三角矩阵乘以三角矩阵的乘法——S 上三角、S^T 下三角，乘积 B 是对称的，而且只需计算上三角就够了（下三角由对称性可得，这就是 dpodi 只填上三角的伏笔）。

B 的第 i 行第 j 列（i ≤ j）：

B_ij = ∑_k S_ik (S^T)_kj = ∑_k S_ik S_jk = ∑_k=max(i,j)^n-1 S_ik S_jk

求和下界是 max(i, j)——因为 S 上三角，S_ik = 0 当 k < i、S_jk = 0 当 k < j，非零项要求 k ≥ i 且 k ≥ j，即 k ≥ max(i, j)；上界到 n-1。对 i ≤ j，max(i,j) = j，所以：

B_ij = ∑_k=j^n-1 S_ik S_jk

这又是一个 O(n-j) 的内积。计算 B 上三角的所有元素总代价也是 O(n³) 量级（带一个比较小的常数）。

整体计算量：和 Gauss-Jordan 精确对比

把所有阶段加起来（以下“flop”指一次浮点加法或乘法，是经典的计数单位）：

阶段	flop 数（近似）
Cholesky 分解	n³/3
三角求逆（反演 R）	n³/6（实际带常数更接近此值的一半）
三角乘法合成 A⁻¹	n³/6
Cholesky 求逆合计	≈ n³/3 + n³/3 = 2n³/3 的下界
通用 Gauss-Jordan 求逆	≈ 2n³

精确数字在不同教材里略有出入（取决于 flop 的定义、是否计入行交换），但结论稳定：Cholesky 两步法求逆的总计算量大约是 Gauss-Jordan 的 1/2 到 1/3。这就是“快一倍”的来源。

更妙的是：如果调用方本来就需要分解（比如回归里既算 β 又算 (X^TX)^-1），那 n³/3 的分解成本是无论如何都要付的，求逆的边际成本只剩 n³/3——和“再跑一遍 Gauss-Jordan 求 (X^TX)^-1”相比，省了不止一半。

数值稳定性：条件数不平方放大

计算量只是一个维度，另一个维度是数值稳定性。Gauss-Jordan 在消元时经常引入大主元相消（一个大数减一个接近的大数，丢掉有效数字），数值误差会显著累积。

Cholesky 两步法则不然。Cholesky 分解本身不需要选主元——正定性天然保证了 R_ii > 0，每一步的除法分母都是正数，没有符号抵消。而三角求逆里的除法也永远除以 R_ii 或 S_ii，分母都不小。

更重要的是：Cholesky 不平方放大条件数。直接对 A 求逆时，相对误差大约正比于 κ(A) · ε_mach；而如果我们走“最小二乘 → 构造 A = X^TX → 求逆”这条路，构造 X^TX 本身就把条件数平方了一次（κ(X^TX) ≈ κ(X)²），这是另一个故事，但至少在 A 已经给定的情况下，Cholesky 求逆不会进一步放大误差。

一句话：Cholesky 求逆是“用结构换精度”的典范。你告诉它“这是个正定矩阵”，它就回报你一个又快又稳定的逆。

逐段拆解：真实的求逆代码

铺垫完了数学，现在对着源码看。函数分为两段：算行列式（上一篇已讲，这里略），以及求逆。求逆部分又分两步：

• 第一步：反演 R 得到 R^-1（原地存上三角）
• 第二步：合成 A^-1 = R^-1 (R^-1)^T（填满上三角）

我们一段一段看。

第一步：反演上三角 R

/* ---- 求逆：第一步，反演 R 得到 R⁻¹（仍存上三角）---- */
if (job % 10 != 0) {
    for (int k = 0; k < n; k++) {
        a[k + k * lda] = 1.0 / a[k + k * lda];  /* R⁻¹ 对角 */
        double t = -a[k + k * lda];
        dscal(k, t, &a[0 + k * lda], 1);         /* 缩放本列上方 */
        for (int j = k + 1; j < n; j++) {
            double tj = a[k + j * lda];
            a[k + j * lda] = 0.0;
            daxpy(k + 1, tj, &a[0 + k * lda], 1, &a[0 + j * lda], 1);
        }
    }
}

这段代码的循环顺序看起来有点奇怪——不是“对每一列 j 从下往上算”，而是“对每个对角元 k 做一次结构化的更新”。这是工程实现和教科书推导之间的一道鸿沟。我们把它对到刚才推导的公式上。

外层循环：逐个对角元 k 处理

for (int k = 0; k < n; k++) {
    a[k + k * lda] = 1.0 / a[k + k * lda];  /* R⁻¹ 对角 */

第 k 次迭代，先把 R 的第 k 个对角元 R_kk 原地替换成它的倒数 S_kk = 1 / R_kk。这对应我们推导的 S_ii = 1/R_ii。

但这里有个微妙的“时机”问题。注意此时的 a[k + k*lda] 已经是 S_kk（取过倒数了），下面这一行：

double t = -a[k + k * lda];

把 t 设为 -S_kk = -1/R_kk。然后——

dscal(k, t, &a[0 + k * lda], 1);  /* 缩放本列上方 */

dscal(n, alpha, x, incx) 是 BLAS 的“向量数乘”：x ← α · x。这里它把第 k 列的前 k 个元素（即 a[0 + k*lda] 到 a[k-1 + k*lda]，共 k 个）整体乘以 t = -S_kk。

为什么？因为在求逆公式 S_ij = -(∑_k=i+1^j R_ik S_kj) / R_ii 里，分母 R_ii 会被反复用到——把这一列的上半段预先乘以 -1/R_ii，相当于把除法“折叠”进了乘法。这是一种典型的循环不变量外提：与其在每个内层迭代里都做一次除法，不如一次性把整列的公共因子提出来。

不过这里有个小细节需要解释——为什么缩放的是第 k 列的前 k 个元素？这就要看 dscal 之后那一段 daxpy 的循环了。

内层循环：用第 k 列更新它右边的所有列

for (int j = k + 1; j < n; j++) {
    double tj = a[k + j * lda];
    a[k + j * lda] = 0.0;
    daxpy(k + 1, tj, &a[0 + k * lda], 1, &a[0 + j * lda], 1);
}

daxpy(n, alpha, x, incx, y, incy) 是 BLAS 的“axpy”操作：y ← α · x + y。这里它把第 j 列的前 k+1 个元素（y = a[0 + j*lda] 起，共 k+1 个）更新为：原值加上 t_j 倍的第 k 列前 k+1 个元素（x = a[0 + k*lda] 起）。

注意循环开始前，第 k 列前 k 个元素已经被 dscal 乘了 -S_kk，第 k 个位置（即对角元位置）现在是 S_kk = 1/R_kk。所以第 k 列前 k+1 个元素此时是：(-S_kk · R_0k, -S_kk · R_1k, …, -S_kk · R_k-1,k, S_kk)。

而第 j 列（j > k）的前 k+1 个元素，此时还是 R 的内容：(R_0j, R_1j, …, R_kj)。

把它们对应位置相加（乘以 t_j = R_kj）后，第 i 行（0 ≤ i ≤ k）的新值是：

R_ij + R_kj · (-S_kk · R_ik), i < k

R_kj + R_kj · S_kk, i = k

注意 R_kk · S_kk = 1，所以上面第二式可以写成 R_kj · (1 + S_kk) 吗？不，这里要小心——当 i = k 时，“第 k 列的第 k 个元素”已经是 S_kk，但 R_kk · S_kk = 1，所以“第 k 列第 k 行的位置”对 t_j · x_k 的贡献是 t_j · S_kk = R_kj · S_kk。然后 y_k = R_kj，更新后 y_k ← R_kj + R_kj S_kk？

读者如果跟着推到这里会卡住——真实的 LINPACK 风格代码在这里做了一件事：它不是直接维护 R 或 S，而是维护一个“中间矩阵”，通过一系列 rank-1 更新逐步把 R 变成 S。每一轮 k 循环都完成了一部分反演。这种“原地 rank-1 累积”的实现，比“显式按公式 S_ij = -∑/diag 一个一个算”更高效，因为所有运算都打包成了 BLAS 调用，能直接吃到向量化加速。

要害在于：这段代码不是数学公式的逐字翻译，而是等价但更高效的算法重排。数学公式告诉你“每个 S_ij 怎么算”，代码把它重排成“每个 R_kj 如何传播到其他位置”。重排前后数学等价，但内存访问模式天差地别——重排后每一步都是连续内存的 BLAS 调用，对 CPU 缓存友好。

对角元取倒数：S_kk = 1 / R_kk

回到循环开头那一行：

a[k + k * lda] = 1.0 / a[k + k * lda];

这是公式 S_kk = 1 / R_kk 的直接翻译。一个除法，干脆利落。整个三角求逆里，除法只发生在对角元上——非对角元的“除以 R_ii”被合并到了 dscal 里，变成了乘法。这是性能的关键之一：除法比乘法慢 3~5 倍，能合并就合并。

内层循环把第 k 行清零

注意这一行：

a[k + j * lda] = 0.0;

把第 j 列的第 k 行清零，然后才调用 daxpy。为什么？因为 daxpy 是 y ← αx + y，会把 y 的原值加进来。但这里我们不想保留原值 R_kj（它已经被存到 tj 里了），所以要先把那个位置清零，让 daxpy 干净地写入新值。这是 LINPACK 代码里反复出现的模式：用 daxpy 做 rank-1 更新之前，先把要覆盖的位置清零。

读完这一段，你应该能体会到 dpodi 第一步的精妙之处：它把“三角矩阵求逆”这个看似需要 n² 次独立回代的操作，重构成了 n 次“对角元取倒数 + 一列缩放 + 一连串 rank-1 更新”，每次更新都是 BLAS 的连续内存操作。这就是为什么实测性能远高于教科书伪代码。

第二步：合成 A⁻¹ = R⁻¹ · (R⁻¹)ᵀ

经过第一步，a 的上三角里存的是 S = R^-1（下三角还是原始 A 的下三角，未被使用）。第二步要算 B = S · S^T：

/* ---- 第二步，合成 A⁻¹ = R⁻¹·(R⁻¹)ᵀ（填满上三角）---- */
for (int j = 0; j < n; j++) {
    for (int k = 0; k < j; k++) {
        double t = a[k + j * lda];
        daxpy(k + 1, t, &a[0 + j * lda], 1, &a[0 + k * lda], 1);
    }
    double t = a[j + j * lda];
    dscal(j + 1, t, &a[0 + j * lda], 1);
}

对应数学公式 B_ij = ∑_m=j^n-1 S_im S_jm（i ≤ j）。这段代码的策略是：对每一列 j，把它的元素 S_kj 当作系数去更新它左边的所有列 k（k），最后再缩放第 j 列本身。

内层循环：用第 j 列更新它左边的所有列

for (int k = 0; k < j; k++) {
    double t = a[k + j * lda];
    daxpy(k + 1, t, &a[0 + j * lda], 1, &a[0 + k * lda], 1);
}

第 j 列第 k 行的元素 S_kj 取出来作为标量 t（注意 k < j，所以这是上三角里的元素）。然后 daxpy(k+1, t, &a[0 + j*lda], 1, &a[0 + k*lda], 1) 按 BLAS 约定是 y ← αx + y：这里 x 是第 j 列、y 是第 k 列，所以是把第 j 列的前 k+1 个元素乘以 t、累加到第 k 列的前 k+1 个元素上——被更新的是第 k 列。

这一步在算什么？它把第 k 列前 k+1 个位置更新为 S_ik + S_kj·S_ij（i ≤ k；此时第 j 列仍存着原始 S_ij）。对照目标 B_ik = ∑_m S_im S_km：每一项 S_imS_km，把 S_km 当作标量去缩放第 m 列、再累加到第 k 列——而这里的 m 正是 j（j>k），所以这次累加贡献的是 S_kj · S_ij。

由于外层 j 从 k+1 跑到 n-1，位置 (i,k)（i ≤ k）最终累加到 ∑_j=k+1^n-1 S_kj S_ij——只差 m=k 的对角项没补上。

收尾：处理对角贡献并缩放

double t = a[j + j * lda];
dscal(j + 1, t, &a[0 + j * lda], 1);

取出第 j 列的对角元 S_jj 作为标量，把第 j 列前 j+1 个元素整体乘以 S_jj。

这一步同时完成了两件事：

1. 补上对角贡献 S_ij · S_jj：目标 B_ij = ∑_m=j^n-1 S_im S_jm 里，m = j 这一项正是 S_ij S_jj（i ≤ j）。此时第 j 列前 j+1 个位置仍是原始 S_ij，dscal 整体乘以 S_jj 后，它们便从 S_ij 变成 S_ij · S_jj。
2. 把第 j 列定格为 B 的第 j 列：把这项和内层循环的累加合起来，位置 (i,j)（i ≤ j）最终等于 S_ij · S_jj + ∑_m=j+1^n-1 S_imS_jm = ∑_m=j^n-1 S_imS_jm = B_ij。核心思想是：用一次 dscal 同时完成“补对角贡献”和“乘上 S_jj”两件事，又一次体现了循环不变量外提的工程智慧。

读完这两步，你大概会感慨：真实的数值代码不是数学公式的直接翻译，而是一连串巧妙的算法重排，让每一步都落在 BLAS 这种高度优化的基本操作上。这是教科书伪代码和 LINPACK 生产代码之间最深的一道鸿沟。

工程细节一：求逆结果只在上三角，必须 symmetrize

这是最容易被忽略、却最坑调用方的一点。看函数头注释：

/* job%10 != 0 : 原地把 A 的上三角替换为 A⁻¹ 的上三角（下三角由对称性可得）
 * A⁻¹ = R⁻¹ · (R⁻¹)ᵀ—— 先反演上三角 R，再做三角乘法合成 */

注意那一句“下三角由对称性可得”。意思是：跑完之后，a 的下三角里还是原始 A 的下三角（或者留下的“未被使用”的下三角），而不是 A^-1 的下三角。调用方必须自己把上三角的值镜像填到下三角，才能得到完整的 A^-1。

为什么这样设计？两个原因：

1. 对称矩阵的逆仍是对称矩阵，下三角完全是上三角的镜像，没必要浪费计算量去算两遍。
2. 内存零浪费：原地存储意味着不需要额外的 n × n 缓冲区。在 1978 年内存以 KB 计的年代，这是必须的；到今天依然有意义——一个 10000 × 10000 的 double 矩阵占 800 MB，能省一份就省一份。

但代价是：调用方必须知道这个约定。如果直接拿程序的输出当 A^-1 用、不做对称化，下三角就全是错的（还是原始 A 的数据）。这是一个典型的“接口陷阱”——函数名说“求逆”，但实际只给了“半个逆”。

在工业软件底层代码里，调用之后的标准模式长这样（伪代码）：

dpodi(a, lda, n, det, 11);  /* job=11 表示同时算行列式和求逆 */
/* 现在的上三角是 A⁻¹ 的上三角 */
for (int j = 0; j < n; j++) {
    for (int i = j + 1; i < n; i++) {
        a[i + j * lda] = a[j + i * lda];  /* 镜像：把上三角的值复制到下三角 */
    }
}
/* 现在的 a 才是完整的 A⁻¹ */

这一步叫 symmetrize（对称化）。它不在内部做，而是留给调用方——为什么？因为有些调用方根本不需要下三角。比如只算 A^-1 b（解方程组），用上三角的 A^-1 配合 BLAS 的 dtrmv（三角矩阵-向量乘法）就够了，根本不用对称化。把对称化放到调用方，让“需要完整矩阵”和“只需要上三角”的两种调用方各取所需，是更灵活的设计。

这就是“半成品约定”的工程智慧：函数只做最小必要的工作，把可选的后续留给调用方。教科书从来不讲这种约定，但每个用过 LINPACK/LAPACK 的人都踩过这个坑。

工程细节二：三角矩阵求逆的逐步代数（3×3 例子）

为了让你彻底信服“三角求逆 = 回代”，我们手算一个 3×3 的例子。设：

R = ( r₀₀ r₀₁ r₀₂ ; 0 r₁₁ r₁₂ ; 0 0 r₂₂ ) , S = R^-1 = ( s₀₀ s₀₁ s₀₂ ; 0 s₁₁ s₁₂ ; 0 0 s₂₂ )

要求 R S = I。逐个元素算：

对角元（直接取倒数）：

s₀₀ = 1/r₀₀, s₁₁ = 1/r₁₁, s₂₂ = 1/r₂₂

第 2 列的非对角元（从下往上）：

s₁₂：从 R S = I 的第 1 行第 2 列，∑_k r_1k s_k2 = 0，即 r₁₁ s₁₂ + r₁₂ s₂₂ = 0，所以

s₁₂ = -r₁₂ s₂₂ / r₁₁ = -r₁₂ / (r₁₁ r₂₂)

s₀₂：第 0 行第 2 列，r₀₀ s₀₂ + r₀₁ s₁₂ + r₀₂ s₂₂ = 0，所以

s₀₂ = -(r₀₁ s₁₂ + r₀₂ s₂₂) / r₀₀

注意 s₁₂ 已经在上一步算好了，这里直接代入——这就是“回代”的含义。每个元素的计算依赖同列更靠下的元素，所以必须从下往上算。

第 1 列的非对角元：

s₀₁：第 0 行第 1 列，r₀₀ s₀₁ + r₀₁ s₁₁ = 0，所以

s₀₁ = -r₀₁ s₁₁ / r₀₀ = -r₀₁ / (r₀₀ r₁₁)

最终：

R^-1 = ( 1/r₀₀ , -r₀₁/(r₀₀ r₁₁) , (r₀₁ r₁₂ - r₀₂ r₁₁)/(r₀₀ r₁₁ r₂₂) ; 0 , 1/r₁₁ , -r₁₂/(r₁₁ r₂₂) ; 0 , 0 , 1/r₂₂ )

数一下运算量：3 个对角元（3 次除法）+ 3 个非对角元（每个 1 次内积 + 1 次除法）。对一般的 n，三角求逆大约 ½ n(n-1) 个非对角元，每个 O(n) 内积——总 flop 约 n³/6 量级。

对比通用求逆：对一个 3 × 3 的稠密矩阵做 Gauss-Jordan 消元，每个元素都要参与行变换，总 flop 约 2n³ = 54。而三角求逆只要 n³/6 ≈ 4.5。差距随 n 增大而放大。这就是结构带来的红利。

工程细节三：行列式作为求逆循环的副产品

现在看函数最前面那段算行列式的代码。注意它的位置——在求逆之前，且单独由 job / 10 != 0 控制：

/* ---- 行列式 ---- */
if (job / 10 != 0) {
    det[0] = 1.0; det[1] = 0.0;
    const double ten = 10.0;
    for (int i = 0; i < n; i++) {
        det[0] *= a[i + i * lda] * a[i + i * lda];  /* Π R(i,i)² */
        if (det[0] == 0.0) break;
        while (det[0] < 1.0) { det[0] *= ten; det[1] -= 1.0; }
        while (det[0] >= ten) { det[0] /= ten; det[1] += 1.0; }
    }
}

这里的关键观察：这段代码读的是 a[i + i*lda]，即 R 的对角元 R_ii。它必须在求逆之前跑——因为求逆的第一步就会把 R_ii 替换成 1/R_ii，到时候就再也读不到原始的 R_ii 了。这就是为什么行列式代码放在求逆代码之前。

数学上 det(A) = det(R^T R) = det(R)² = ∏_i R_ii²，所以代码就是读 R 的对角元，平方，累乘。几乎零额外成本——分解已经把 R_ii 都算好了，这里只是顺着对角线读一遍。

手工科学计数法防溢出（简述）

det[0] 存尾数，永远保持在 [1, 10) 区间；det[1] 存 10 的指数。每乘一次就立刻归一化：小于 1 就乘 10、指数减 1；大于等于 10 就除以 10、指数加 1。于是 det(A) = det[0] × 10^det[1]，永远不溢出。

为什么必须这么做？一个 100 × 100 的矩阵，∏ R_ii² 轻松达到 10²⁰⁰，中间过程还可能更高。double 的上限约 1.8 × 10³⁰⁸，但只要矩阵规模再大一点或者 R_ii 略大于 1，就会爆。手工科学计数法是“教科书数学”和“生产代码”之间最典型的鸿沟——数学上 ∏ R_ii² 是个干净的公式，计算机上直接乘会炸。工业软件用两行 while 解决了它。

（这部分逻辑上一篇已经详细讲过，这里不展开。本文要强调的是：行列式计算是求逆过程的“顺路产物”——调用方既能拿到 A^-1，又能几乎免费地拿到 det(A)，一次调用两份收获。这在多元统计里特别有用：判别分析要算 |W|（组内协方差行列式）、贝叶斯模型选择要算 |Σ|，都是和协方差求逆一起做的。job = 11 就是一次拿到逆和行列式的“组合套餐”。）

job 参数的设计：位运算复用一个 int

dpodi 的接口用一个 int 的两个十进制位同时控制“算不算行列式”和“算不算逆”：

job 值	job/10（算行列式？）	job%10（算逆？）	行为
0	0	0	啥也不算（合法但无用）
1	0	1	只算逆
10	1	0	只算行列式
11	1	1	都算

这是 1970 年代 Fortran 接口的典型风格——用一个 int 的十进制位打包多个布尔开关，省下函数参数。现代 C 代码会用位掩码（JOB_DET | JOB_INV），但 LINPACK 选择了十进制位。读这种老代码时要注意：job / 10 和 job % 10 不是在做数学除法，是在拆开关。

把整条链串起来：求逆在统计计算里出现在哪里

理解了两步法，回头看统计软件的底层，到处都是它的影子：

线性回归：算 (X^TX)^-1

最小二乘回归 β̂ = (X^TX)^-1 X^T y，要算 (X^TX)^-1。X^TX 是对称正定，走 Cholesky：分解、求逆。回归系数的标准误、置信区间、t 统计量，全都依赖 (X^TX)^-1 的对角元。这个矩阵的精度直接决定了统计推断的可靠性——所以“又快又稳”的 Cholesky 求逆在这里是刚需。

协方差矩阵求逆：多元统计的常客

判别分析、多元回归、PCA、因子分析、结构方程……都要算协方差矩阵 Σ 的逆。Σ 对称正定（假设非奇异），又是 Cholesky 求逆的舞台。判别分析里的马氏距离 d² = (x - μ)^T Σ^-1 (x - μ)，每判一个样本都要算一次 Σ^-1 x，预计算 Σ^-1 用 dpodi，运行时只做矩阵-向量乘法。

贝叶斯推断：从多元正态后验采样

贝叶斯回归的后验往往是多元正态 N(μ, Σ)，采样需要算 μ + L z，其中 L L^T = Σ。这里的 L 就是 Cholesky 因子，采样过程不需要显式求 Σ^-1，但如果要做后验预测、算协方差调整，就得让求逆出马。MCMC 的每一步迭代都可能触发一次分解+求逆，性能至关重要。

优化：牛顿法的 Hessian 求逆

牛顿法更新 θ ← θ - H^-1 g，其中 H 是 Hessian。如果 H 正定（凸优化场景），又是 Cholesky 求逆。信任域方法、内点法、L-BFGS 的初始化，底层都藏着 Cholesky。

与 Gauss-Jordan 的全面对比

最后做一个总账对比，把性能、稳定性、适用性都摆出来：

维度	Cholesky 两步法求逆	Gauss-Jordan 消元求逆
适用矩阵	对称正定	任意非奇异
计算量（flop）	≈ n³/3 + n³/3 ≈ 2n³/3	≈ 2n³
实测速度	快一倍左右	基准
主元选取	不需要（正定性保证）	必须（部分主元或全主元）
行交换	无	有（破坏原矩阵布局）
数值稳定性	极好（无大主元相消）	中等（依赖主元策略）
内存	原地，无额外分配	原地，但需要 pivot 数组
副产品	行列式（几乎免费）	行列式（要单独算）
结果存储	只填上三角，需 symmetrize	填满全矩阵

一句话总结：Cholesky 求逆是“用适用性换性能”的典范——你告诉它矩阵是正定的，它就回报你 2 倍速度和更高的精度。如果你的矩阵恰好是对称正定（而统计计算里一半的矩阵都是），没有理由用 Gauss-Jordan。

总结：教科书不讲的 5 个点

1. 核心恒等式：A^-1 = (R^T R)^-1 = R^-1 (R^-1)^T。把稠密求逆拆成两个三角求逆+一次三角乘法，是两步法的数学根基。
2. 三角矩阵求逆 = 回代：每个元素用一个内积加一次除法得到，对角元直接取倒数，非对角元从下往上回代，总代价 O(n³/6)。
3. 真实代码不是公式的直接翻译：把求逆重排成“对角元取倒数 + 列缩放 + rank-1 更新”的循环结构，每一步都落在 BLAS 的 dscal/daxpy 上，吃满向量化加速。
4. 求逆结果只在上三角，必须 symmetrize：函数只算“半个逆”——上三角是 A^-1 的上三角，下三角还是原始数据。调用方要自己镜像填充，这是 LINPACK 风格的“半成品约定”，让“需要完整矩阵”和“只需要上三角”的两种调用方各取所需。
5. 行列式是求逆的副产品：det(A) = ∏ R_ii²，几乎零成本（顺着对角线读一遍），还用手工科学计数法防溢出。job = 11 一次拿到逆和行列式，是多元统计的“组合套餐”。

如果只记一句话，那就是：通用的求逆算法把矩阵当成一坨无差别的数；Cholesky 求逆尊重矩阵的结构，用结构信息换计算量和精度。工程上凡是对称正定，第一选择永远是 Cholesky，没有例外。

下一篇，我们会拆 dtrsl——三角方程组求解，也就是 R x = b 怎么解。那是个看似平凡、但藏着“前向代入 vs 后向代入”和“右端向量批处理”两个工程细节的故事。同样是教科书一句带过、生产代码精妙绝伦。