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

上一篇我们拆了 Cholesky——30 行 C 代码把对称正定矩阵 $A$ 分解成 $A=R^{T}R$，$R$ 是上三角。当时只用半篇文章的篇幅扫过Cholesky 的求逆部分，留了个悬念：所谓"两步法求逆"到底是怎么把 $A^{-1}$ 算出来的？为什么它就比通用方法快一倍？求逆结果为什么是个"半成品"？

这一篇就把这件事讲透。聚焦点是Cholesky 的后半段（求逆部分），而不是 分解部分。这是教科书最含糊其辞、而工程代码最值得玩味的一段。读完你会明白：

• 为什么三角矩阵的逆只需要 $O(n^2)$ 回代而不是 $O(n^3)$ 消元；
• 为什么 求出来的逆只填了上三角，调用方还得自己"对称化"一次；
• 为什么这套两步法的总计算量是通用 Gauss-Jordan 的一半左右；
• 为什么行列式可以"顺路"从求逆循环里捞出来。

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

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

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

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

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

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

Cholesky 求逆的核心思想，就是把这份结构信息转化为计算量的节省。 先分解（上一篇讲过，$n^3\mathrm{/}3$ 量级），再分别处理两个三角矩阵（这一篇要讲，各 $n^3\mathrm{/}3$ 量级），合起来比 Gauss-Jordan 的 $2n^3\mathrm{/}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}$，满足 $RS=I$。展开 $RS=I$ 的第 $i$ 行第 $j$ 列元素（约定 $i\le j$，因为 $S$ 也是上三角，$i>j$ 时 $S_{ij}=0$）：

${\sum }_{k=i}^j{R}_{ik}{S}_{kj}={\delta }_{ij}$

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

把 $k=i$ 那一项拎出来：

$R_{ii}{S}_{ij}+{\sum }_{k=i+1}^j{R}_{ik}{S}_{kj}={\delta }_{ij}$

解出 $S_{ij}$：

$S_{ij}=\frac{{\delta }_{ij}-{\sum }_{k=i+1}^j{R}_{ik}{S}_{kj}}{R_{ii}}$

分两种情况看：

• 对角元（$i=j$）：${\delta }_{ii}=1$，求和项空集，得到 $S_{ii}=1\mathrm{/}{R}_{ii}$——对角元就是直接取倒数。
• 非对角元（）：${\delta }_{ij}=0$，得到 $S_{ij}=-\frac{{\sum }_{k=i+1}^j{R}_{ik}{S}_{kj}}{R_{ii}}$——一个内积再除以对角元。

注意求和 ${\sum }_{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^2)$，所以总代价 $O(n^3)$——但常数因子只是通用方法的一个零头（下面会精确算）。

合成 $A^{-1}$：三角乘三角

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

$B$ 的第 $i$ 行第 $j$ 列（$i\le j$）：

$B_{ij}={\sum }_k{S}_{ik}(S^T{)}_{kj}={\sum }_k{S}_{ik}{S}_{jk}={\sum }_{k=\max (i,j)}^{n-1}{S}_{ik}{S}_{jk}$

求和下界是 $\max (i,j)$——因为 $S$ 上三角，$S_{ik}=0$ 当 、$S_{jk}=0$ 当 ，非零项要求 $k\ge i$ 且 $k\ge j$，即 $k\ge \max (i,j)$；上界到 $n-1$。对 $i\le j$，$\max (i,j)=j$，所以：

$B_{ij}={\sum }_{k=j}^{n-1}{S}_{ik}{S}_{jk}$

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

整体计算量：和 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。 这就是"快一倍"的来源。

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

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

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

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

更重要的是：Cholesky 不平方放大条件数。 直接对 $A$ 求逆时，相对误差大约正比于 $\kappa (A)\cdot {\epsilon }_{\text{mach}}$；而如果我们走"最小二乘 → 构造 $A=X^{T}X$ → 求逆"这条路，构造 $X^{T}X$ 本身就把条件数平方了一次（$\kappa (X^{T}X)\approx \kappa (X{)}^2$），这是另一个故事，但至少在 $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\mathrm{/}{R}_{kk}$。这对应我们推导的 $S_{ii}=1\mathrm{/}{R}_{ii}$。

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

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

把 $t$ 设为 $-S_{kk}=-1\mathrm{/}{R}_{kk}$。然后——

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

dscal(n, alpha, x, incx) 是 BLAS 的"向量数乘"：$x\leftarrow \alpha \cdot x$。这里它把第 $k$ 列的前 $k$ 个元素（即 a[0 + k*lda] 到 a[k-1 + k*lda]，共 $k$ 个）整体乘以 $t=-S_{kk}$。

为什么？因为在求逆公式 $S_{ij}=-\frac{{\sum }_{k=i+1}^j{R}_{ik}{S}_{kj}}{R_{ii}}$ 里，分母 $R_{ii}$ 会被反复用到——把这一列的上半段预先乘以 $-1\mathrm{/}{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\leftarrow \alpha \cdot 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\mathrm{/}{R}_{kk}$。所以第 $k$ 列前 $k+1$ 个元素此时是：$(-S_{kk}\cdot R_{0k},-S_{kk}\cdot R_{1k},\dots ,-S_{kk}\cdot R_{k-1,k},S_{kk})$。

而第 $j$ 列（$j>k$）的前 $k+1$ 个元素，此时还是 $R$ 的内容：$(R_{0j},R_{1j},\dots ,R_{kj})$。

把它们对应位置相加（乘以 $t_j=R_{kj}$）后，第 $i$ 行（$0\le i\le k$）的新值是：

$R_{ij}+R_{kj}\cdot (-S_{kk}\cdot R_{ik}),$Rij​+Rkj​⋅(−Skk​⋅Rik​),i<k

$R_{kj}+R_{kj}\cdot S_{kk},$Rkj​+Rkj​⋅Skk​,i=k

注意 $R_{kk}\cdot S_{kk}=1$，所以上面第二式可以写成 $R_{kj}\cdot (1+S_{kk})$ 吗？不，这里要小心——当 $i=k$ 时，"第 $k$ 列的第 $k$ 个元素"已经是 $S_{kk}$，但 $R_{kk}\cdot S_{kk}=1$，所以"第 $k$ 列第 $k$ 行的位置"对 $t_j\cdot x_k$ 的贡献是 $t_j\cdot S_{kk}=R_{kj}\cdot S_{kk}$。然后 $y_k=R_{kj}$，更新后 $y_k\leftarrow R_{kj}+R_{kj}{S}_{kk}$？

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

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

对角元取倒数：$S_{kk}=1\mathrm{/}{R}_{kk}$

回到循环开头那一行：

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

这是公式 $S_{kk}=1\mathrm{/}{R}_{kk}$ 的直接翻译。一个除法，干脆利落。整个三角求逆里，除法只发生在对角元上——非对角元的"除以 $R_{ii}$"被合并到了 dscal 里，变成了乘法。这是性能的关键之一：除法比乘法慢 3~5 倍，能合并就合并。

内层循环把第 $k$ 行清零

注意这一行：

a[k + j * lda] = 0.0;

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

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

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

经过第一步，a 的上三角里存的是 $S=R^{-1}$（下三角还是原始 $A$ 的下三角，未被使用）。第二步要算 $B=S\cdot 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}={\sum }_{m=j}^{n-1}{S}_{im}{S}_{jm}$（$i\le j$）。这段代码的策略是：对每一列 $j$，把它的元素 $S_{kj}$ 当作系数去更新它左边的所有列 $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$（注意 ，所以这是上三角里的元素）。然后 daxpy(k+1, t, &a[0 + j*lda], 1, &a[0 + k*lda], 1) 按 BLAS 约定是 $y\leftarrow \alpha x+y$：这里 $x$ 是第 $j$ 列、$y$ 是第 $k$ 列，所以是把第 $j$ 列的前 $k+1$ 个元素乘以 $t$、累加到第 $k$ 列的前 $k+1$ 个元素上——被更新的是第 $k$ 列。

这一步在算什么？它把第 $k$ 列前 $k+1$ 个位置更新为 $S_{ik}+S_{kj}\cdot S_{ij}$（$i\le k$；此时第 $j$ 列仍存着原始 $S_{ij}$）。对照目标 $B_{ik}={\sum }_m{S}_{im}{S}_{km}$：每一项 $S_{im}{S}_{km}$，把 $S_{km}$ 当作标量去缩放第 $m$ 列、再累加到第 $k$ 列——而这里的 $m$ 正是 $j$（$j>k$），所以这次累加贡献的是 $S_{kj}\cdot S_{ij}$。

由于外层 $j$ 从 $k+1$ 跑到 $n-1$，位置 $(i,k)$（$i\le k$）最终累加到 ${\sum }_{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}\cdot S_{jj}$：目标 $B_{ij}={\sum }_{m=j}^{n-1}{S}_{im}{S}_{jm}$ 里，$m=j$ 这一项正是 $S_{ij}{S}_{jj}$（$i\le j$）。此时第 $j$ 列前 $j+1$ 个位置仍是原始 $S_{ij}$，dscal 整体乘以 $S_{jj}$ 后，它们便从 $S_{ij}$ 变成 $S_{ij}\cdot S_{jj}$。
2. 把第 $j$ 列定格为 $B$ 的第 $j$ 列：把这项和内层循环的累加合起来，位置 $(i,j)$（$i\le j$）最终等于 $S_{ij}\cdot S_{jj}+{\sum }_{m=j+1}^{n-1}{S}_{im}{S}_{jm}={\sum }_{m=j}^{n-1}{S}_{im}{S}_{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\times n$ 缓冲区。在 1978 年内存以 KB 计的年代，这是必须的；到今天依然有意义——一个 $10000\times 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 的例子。设：

​r00​00​r01​r11​0​r02​r12​r22​​

​,S=R−1=

​s00​00​s01​s11​0​s02​s12​s22​​

​

要求 $RS=I$。逐个元素算：

对角元（直接取倒数）：

$s_{00}=1\mathrm{/}{r}_{00},$s00​=1/r00​,s11​=1/r11​,s22​=1/r22​

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

$s_{12}$：从 $RS=I$ 的第 1 行第 2 列，${\sum }_k{r}_{1k}{s}_{k2}=0$，即 $r_{11}{s}_{12}+r_{12}{s}_{22}=0$，所以

$s_{12}=-r_{12}{s}_{22}\mathrm{/}{r}_{11}=-r_{12}\mathrm{/}(r_{11}{r}_{22})$

$s_{02}$：第 0 行第 2 列，$r_{00}{s}_{02}+r_{01}{s}_{12}+r_{02}{s}_{22}=0$，所以

$s_{02}=-(r_{01}{s}_{12}+r_{02}{s}_{22})\mathrm{/}{r}_{00}$

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

第 1 列的非对角元：

$s_{01}$：第 0 行第 1 列，$r_{00}{s}_{01}+r_{01}{s}_{11}=0$，所以

$s_{01}=-r_{01}{s}_{11}\mathrm{/}{r}_{00}=-r_{01}\mathrm{/}(r_{00}{r}_{11})$

最终：

​1/r00​00​−r01​/(r00​r11​)1/r11​0​(r01​r12​−r02​r11​)/(r00​r11​r22​)−r12​/(r11​r22​)1/r22​​

​

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

对比通用求逆：对一个 $3\times 3$ 的稠密矩阵做 Gauss-Jordan 消元，每个元素都要参与行变换，总 flop 约 $2n^3=54$。而三角求逆只要 $n^3\mathrm{/}6\approx 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\mathrm{/}{R}_{ii}$，到时候就再也读不到原始的 $R_{ii}$ 了。 这就是为什么行列式代码放在求逆代码之前。

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

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

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

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

（这部分逻辑上一篇 #1 已经详细讲过，这里不展开。本文要强调的是：行列式计算是求逆过程的"顺路产物"——调用方既能拿到 $A^{-1}$，又能几乎免费地拿到 $\det (A)$，一次调用两份收获。这在多元统计里特别有用：判别分析要算 $\mathrm{\mid }W\mathrm{\mid }$（组内协方差行列式）、贝叶斯模型选择要算 $\mathrm{\mid }\mathrm{\Sigma }\mathrm{\mid }$，都是和协方差求逆一起做的。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^{T}X{)}^{-1}$

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

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

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

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

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

优化：牛顿法的 Hessian 求逆

牛顿法更新 $\theta \leftarrow \theta -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^3\mathrm{/}6)$。
3. 真实代码不是公式的直接翻译： 把求逆重排成"对角元取倒数 + 列缩放 + rank-1 更新"的循环结构，每一步都落在 BLAS 的 dscal/daxpy 上，吃满向量化加速。
4. 求逆结果只在上三角，必须 symmetrize：函数只算"半个逆"——上三角是 $A^{-1}$ 的上三角，下三角还是原始数据。调用方要自己镜像填充，这是 LINPACK 风格的"半成品约定"，让"需要完整矩阵"和"只需要上三角"的两种调用方各取所需。
5. 行列式是求逆的副产品：$\det (A)=\prod R_{ii}^2$，几乎零成本（顺着对角线读一遍），还用手工科学计数法防溢出。job = 11 一次拿到逆和行列式，是多元统计的"组合套餐"。

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

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