对称不定分解（Bunch-Kaufman）：为何 Cholesky 不够用

写在前面：上一篇留下的一个洞

前面讲 Cholesky 时，我们把那 30 行代码捧上了天——线性回归的正规方程、ANOVA 的二次型、PCA 的预条件、Kalman 滤波的初始化、贝叶斯采样的每一步 MCMC——都建立在对称正定矩阵之上，都跑那 30 行 dpofa。

可现实里，你拿到的对称矩阵，未必正定。

打开任何一份真实的协方差矩阵，你都有可能撞上：

• 共线性：两个自变量高度相关，$X^{T}X$ 接近奇异，最小特征值趋近于 0；
• 奇异矩阵：虚拟变量陷阱、自由度耗尽，矩阵根本不可逆；
• 负特征值：在约束优化、Hessian 矩阵、二阶导数分析里，对称矩阵出现负特征值是家常便饭——只要目标函数在当前点不是凸的，Hessian 就不是半正定。

这时候如果你还硬上 Cholesky，会发生什么？回头看 dpofa 那行唯一的错误处理：

if (s <= 0.0) return j + 1; // 非正定：返回出错列号

$s$ 是 $R(j,j{)}^2$。一旦矩阵不正定，平方根里面的值就为负，s 为负，函数直接返回——分解失败，整个计算链断裂。你点了"回归"，屏幕报"matrix is singular"或者"system is computationally singular"，然后呢？然后没有然后了。

可问题是：负特征值的矩阵，难道就不能分解了吗？

数学上当然可以。对称矩阵永远可以对角化（谱定理），永远可以写成 $A=Q\mathrm{\Lambda }{Q}^T$。但是谱分解是 $O(n^3)$ 的迭代算法，不稳定也不高效。我们需要一个像 Cholesky 那样直接、$O(n^3\mathrm{/}3)$、数值稳定的分解，但它得能吃下不定矩阵。

1971 年，Bunch 和 Kaufman 在 Numerische Mathematik 上给出了答案：允许 2×2 块作主元，把对称矩阵分解成 $A=UD{U}^T$，其中 $D$ 是由 1×1 和 2×2 块组成的块对角阵。这篇论文的精妙之处不在"允许 2×2"——这事谁都想到——而在于每一步到底该选 1×1 还是 2×2，有一个解析最优判据，判据里那个古怪的常数 $\alpha =(1+\sqrt{17})\mathrm{/}8\approx 1.65$

​)/8≈1.65，是论文花了大半篇幅推导出来的。

今天我们就用 工业级算法代码的简化版 sym.c，把这个算法一行行拆开看。读懂它，你才真正读懂了统计软件底层对"对称矩阵"的全部处理——不只是正定那一半。

数学原理：从 Cholesky 到 LDLᵀ，再到块对角 D

第一步：把平方根拿掉——LDLᵀ 分解

Cholesky 分解 $A=R^{T}R$ 里有一个开方。开方就是不定矩阵的死穴——负数开不了。所以第一步，把开方拿出来。

任何对称矩阵（不必正定）都可以写成：

$A=LD{L}^T$

其中 $L$ 是单位下三角（对角元全是 1），$D$ 是对角阵。如果 $A$ 正定，那么 $D$ 全正，再对 $D$ 开方得到 $D^{1\mathrm{/}2}$，就有 $A=(L{D}^{1\mathrm{/}2})(L{D}^{1\mathrm{/}2}{)}^T$，回到了 Cholesky。

LDLᵀ 的好处是：$D$ 的元素可以是负的，不需要开方。只要 $A$ 非奇异，$D$ 的对角元就非零，分解照样进行。

第二步：单对角 $D$ 还不够——为什么需要 2×2 块

但单对角 $D$ 有一个隐患。考虑这个对称矩阵：

它是对称、非奇异（行列式 $-1$），但两个对角元都是 0。如果你硬要写 $A=LD{L}^T$，第一步就要除以 $A(1,1)=0$——直接爆炸。可这个矩阵本身完全正常，它的特征值是 $\pm 1$，一正一负。

问题出在哪？ 出在"主元"选错了。如果允许 $D$ 的第一个"块"是一个 2×2 块而不是 1×1 标量：

$D=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),$D=(01​10​),L=I,U=I

分解直接完成，没有除以零，没有开方。这个 2×2 块的特征值一正一负，完美吃下了原矩阵的不定性。

这就是 Bunch-Kaufman 的核心思想：$D$ 不再是单纯的对角阵，而是一个由 1×1 和 2×2 块拼起来的"块对角阵"。

$A=UD{U}^T,$A=UDUT,D=diag(D1​,D2​,…,Dp​),Di​是1×1或2×2

每个 2×2 块负责"吞掉"一个负特征值附近的奇异性。整个分解仍然是 $O(n^3\mathrm{/}3)$ 量级的直接方法，不迭代、不收敛、不依赖初值。

第三步：每步怎么选？——α 常数的来历

最难的问题来了：每一步，到底用 1×1 主元还是 2×2 主元？

朴素的想法是"谁大用谁"，但这没有稳定性保证。高斯消元里我们学过"部分主元法"——选列中最大的元素作主元——来控制元素增长（element growth），把舍入误差压在可控范围内。Bunch-Kaufman 把这个思想推广到了对称情形，但对称性带来了一个微妙的约束：选一个主元，就要对称地交换对应的行和列（不能只换行不换列，否则对称性破坏）。

Bunch 在 1971 年的论文里做了一件非常硬核的事：他把"元素增长上界"写成主元阈值的函数，然后求这个函数的最小值，得到解析解。

具体地，每一步考察当前活跃子矩阵的第 $k$ 列。记：

• ${\alpha }_{kk}=\mathrm{\mid }A(k,k)\mathrm{\mid }$（对角元的绝对值）
• ${\alpha }_{ik}={\max }_{i\mathrm{\ne }k}\mathrm{\mid }A(i,k)\mathrm{\mid }$（这一列里最大的非对角元绝对值）

候选方案有两个：

1. 用 $A(k,k)$ 作 1×1 主元：稳定性要求 ${\alpha }_{kk}$ 不要比 ${\alpha }_{ik}$ 小太多；
2. 用 $A(k,k)$ 和 $A(i_{\max },k)$ 配对作 2×2 块主元：稳定性要求这个 2×2 块的行列式不能太小。

Bunch 证明：存在一个常数 $\alpha \in (0,1)$，使得只要在每步满足 $\mathrm{\mid }A(k,k)\mathrm{\mid }\ge \alpha \cdot {\max }_i\mathrm{\mid }A(i,k)\mathrm{\mid }$ 时选 1×1，否则考虑 2×2，就能让元素增长因子的上界达到最小。这个最优的 $\alpha$，就是 $(1+\sqrt{17})\mathrm{/}8\approx 1.65$

​)/8≈1.65 的倒数在某种意义下的最优权衡点——精确推导涉及对 2×2 块消元后元素增长率的二次分析，论文里是几页不等式放缩。

之所以是 $\sqrt{17}$

​，不是 $\sqrt{16}$

​ 或 $\sqrt{25}$

​，是因为增长率函数的极值条件恰好满足 $8\alpha -1=\sqrt{17}$

​，化简得 $\alpha =(1+\sqrt{17})\mathrm{/}8$

​)/8。这不是拍脑袋，是 1971 年那篇论文花了大半篇幅算出来的解析最优解。

记住这个数字，等下你会看到它在 sym.c 第 27 行的精确出现。

逐段拆解 dsifa：一个 goto 编织的主元选择机器

打开 sym.c，第一眼你会看到文件头那段郑重其事的注释：

/*============================================================================* *sym.c —— 对称（不定）矩阵的 Bunch-Kaufman 分解与求解 C 移植 *对应 Fortran：dsifa.for / dsisl.for * *数学：任意对称矩阵 A（不必正定）可分解为A = U · D · Uᵀ *其中 U 是单位上三角，D 是由 1×1 与 2×2 块组成的块对角阵。 *每步在第 K 列用阈值准则（常数 α=(1+√17)/8）在「1×1 主元」与 *「2×2 主元块」之间权衡，以保证数值稳定。这是 Cholesky 在 *不定矩阵上的推广。 *============================================================================*/

注意"任意对称矩阵（不必正定）"——这是它和 Cholesky 最本质的区别。cholesky.c 文件头一句是"对称正定矩阵"，这里去掉了"正定"二字，整篇文章的故事就在这两个字里。

第 27 行：那个神秘常数精确出现

int dsifa(double *a, int lda, int n, int *kpvt) {#define A(i,j) a[((i)-1) + ((j)-1)*lda]#define KP(k)kpvt[(k)-1]const double alpha = (1.0 + sqrt(17.0)) / 8.0;int info = 0;int k = n;...

(1.0 + sqrt(17.0)) / 8.0 —— 1971 年 Bunch 论文的最优阈值常数，原封不动地写在这里。整个算法的稳定性，就压在这一行 const double alpha 上。 改掉这个数字，分解仍然能跑（数学上不会错），但元素增长可能失控，舍入误差在某个病态矩阵上会爆炸。SPSS、R、SciPy、MATLAB——全世界所有用 Bunch-Kaufman 的软件，这一行都长一模一样，因为这个值是数学上最优的，没有优化空间。

另一个值得注意的细节是 int k = n;——分解从最后一列往前做，而不是从前往后。这是 LINPACK dsifa 的约定，和 dpofa（Cholesky）从前往后不一样。从后往前的好处是：活跃子矩阵始终在数组前 $k$ 个位置，索引更整齐，便于调用 BLAS。这种"从尾巴向前剥"的写法在 70 年代 Fortran 代码里非常普遍。

L20-L70：1×1 与 2×2 主元的四级判据

这是整段代码最精彩的部分。一个主元的选择，要经过四级判据，对应代码里的四个 goto 跳转。我们逐段看：

L20:km1 = k - 1;absakk = fabs(A(k,k));imax = idamax(k - 1, &A(1,k), 1) + 1;/* +1：C 版 idamax 返回 0-based */colmax = fabs(A(imax,k));if (absakk < alpha * colmax) goto L30;kstep = 1; swap = 0;goto L90;

第一级判据：看当前列的对角元够不够大。absakk 是 $\mathrm{\mid }A(k,k)\mathrm{\mid }$；imax 是这一列（前 $k-1$ 个元素里）绝对值最大的行号，由 BLAS 函数 idamax 给出；colmax 是这个最大绝对值。

判据是：如果 $\mathrm{\mid }A(k,k)\mathrm{\mid }\ge \alpha \cdot {\max }_i\mathrm{\mid }A(i,k)\mathrm{\mid }$，直接用 $A(k,k)$ 作 1×1 主元，不交换（swap = 0）。

这就是我们前面讲的 Bunch 阈值准则的直接翻译：对角元只要不比列里最大元小太多（差距小于 $\alpha \approx 1.65$ 倍），就足够稳定，直接用。这一步命中率最高——大多数良态矩阵的绝大多数列，都在这一步被处理掉，代码极其快。

否则跳到 L30，进入第二级判据：

L30:rowmax = 0.0;imaxp1 = imax + 1;for (j = imaxp1; j <= k; j++)rowmax = (rowmax > fabs(A(imax,j))) ? rowmax : fabs(A(imax,j));if (imax != 1) {jmax = idamax(imax - 1, &A(1,imax), 1) + 1;rowmax = (rowmax > fabs(A(jmax,imax))) ? rowmax : fabs(A(jmax,imax));}if (fabs(A(imax,imax)) < alpha * rowmax) goto L60;kstep = 1; swap = 1;goto L80;

对角元不够大，那就考虑"换一个更大的对角元来作 1×1 主元"。候选是 imax（刚才找到的列最大元所在行）。但换之前要确认：imax 自己作为对角元，够不够大？

为此要算 rowmax：imax 行（在活跃子矩阵里）的最大绝对值。注意代码分两段处理——imax+1 到 k 这一段直接 for 循环扫，而 1 到 imax-1 这一段再次调 idamax。为什么分两段？因为对称矩阵只存了一半（这里存在下三角），所以 imax 行的前半段和后半段在内存里位置不同，得分开扫。

判据：如果 $\mathrm{\mid }A(imax,imax)\mathrm{\mid }\ge \alpha \cdot \text{rowmax}$，那就把 imax 行/列换到 $k$ 位置，作 1×1 主元（swap = 1）。这也是部分主元思想——找一个更稳的对角元。

否则继续到 L60，第三级判据：

L60:if (absakk < alpha * colmax * (colmax / rowmax)) goto L70;kstep = 1; swap = 0;goto L80;

imax 自己也不够格。现在有两个候选：老老实实回到 $A(k,k)$ 作 1×1（不交换），或者干脆上 2×2 块。这一级判据是 $\mathrm{\mid }A(k,k)\mathrm{\mid }$ 和 $\alpha \cdot {\text{colmax}}^2\mathrm{/}\text{rowmax}$ 的比较——为什么是这个古怪表达式？因为这是 2×2 块消元后元素增长率的精确上界，Bunch 论文里推导出来的。如果原对角元比这个上界还大，那 1×1 也够稳，不必上 2×2。

否则跳到 L70，第四级——决定用 2×2 块：

L70:kstep = 2;swap = (imax != km1);

kstep = 2 表示这一步消去两列（一个 2×2 块），swap 决定要不要把 imax 行/列换到 $k-1$ 位置（让 2×2 块由 $A(k-1,k)$ 和 $A(k,k)$ 配对组成）。

四级判据，层层下探，最终一定有一个分支被选中。 这就是 Bunch-Kaufman 的"最优权衡"在代码里的真实样子——不是一句"比较一下大小"，而是四个层次的不等式，每一个都对应一种稳定性分析。

L100-L120：1×1 主元的消元

判据选定 1×1 主元后（无论换不换），进入 L100。先做必要的行列对称交换：

L100:if (kstep == 2) goto L140;/* ===== 1×1 主元 ===== */if (!swap) goto L120;dswap(imax, &A(1,imax), 1, &A(1,k), 1);for (jj = imax; jj <= k; jj++) {j = k + imax - jj; /* 对称地交换行/列 imax 与 k */t = A(j,k);A(j,k) = A(imax,j);A(imax,j) = t;}

注意 dswap 只换了行，后面那个 for 循环才换列——这就是对称主元法的关键：行换完必须列也换，否则对称性破坏。循环里那个 j = k + imax - jj 是反向遍历，巧妙地避免了交换时元素被覆盖。这种细节在教科书里几乎从不出现，但在写实际代码时如果不小心，交换顺序错一步就会把矩阵搞乱。

交换完成后是消元主循环：

L120:for (jj = 1; jj <= km1; jj++) {j = k - jj;mulk = -A(j,k) / A(k,k);daxpy(j, mulk, &A(1,k), 1, &A(1,j), 1);A(j,k) = mulk;}KP(k) = k;if (swap) KP(k) = imax;

对前 $k-1$ 列每一列 $j$，算出乘子 mulk = -A(j,k)/A(k,k)，用 daxpy 把第 $k$ 列的 mulk 倍加到第 $j$ 列上去——这就是高斯消元的对称版本。和 Cholesky 不同的是，这里要除以 $A(k,k)$，而不是开方——这就是为什么不定矩阵也能算（$A(k,k)$ 可以是负的，只要不是 0）。

L140-L160：2×2 主元块的就地消元

最精彩的是 2×2 块的处理。判据选定 2×2 后：

L140:/* ===== 2×2 主元块 ===== */if (!swap) goto L160;dswap(imax, &A(1,imax), 1, &A(1,k-1), 1);for (jj = imax; jj <= km1; jj++) {j = km1 + imax - jj;t = A(j,k-1);A(j,k-1) = A(imax,j);A(imax,j) = t;}t = A(k-1,k);A(k-1,k) = A(imax,k);A(imax,k) = t;

同样是行列对称交换，把 imax 换到 $k-1$ 位置。交换完成后，$A(k-1,k)$ 是连接两个主元列的"桥梁"元素——它就是 2×2 块的非对角元。

然后是核心的 2×2 块消元：

L160:km2 = k - 2;if (km2 != 0) {ak = A(k,k) / A(k-1,k);akm1 = A(k-1,k-1) / A(k-1,k);denom = 1.0 - ak * akm1;for (jj = 1; jj <= km2; jj++) {j = km1 - jj;bk = A(j,k) / A(k-1,k);bkm1 = A(j,k-1) / A(k-1,k);mulk = (akm1 * bk - bkm1) / denom;mulkm1 = (ak * bkm1 - bk) / denom;daxpy(j, mulk, &A(1,k), 1, &A(1,j), 1);daxpy(j, mulkm1, &A(1,k-1), 1, &A(1,j), 1);A(j,k) = mulk;A(j,k-1) = mulkm1;}}

注意这里的数学：所有计算都除以 $A(k-1,k)$，而不是除以某个对角元。这是因为 2×2 块的"逆"是通过它的非对角元来表达的——前面讲过，2×2 块 的行列式是 $ac-b^2$，对于不定矩阵这个行列式是负的（特征值一正一负），但 $b$（也就是 $A(k-1,k)$）非零，所以可以拿 $b$ 当"桥梁"做消元。

denom = 1.0 - ak * akm1 就是归一化后的 2×2 块行列式。如果这个 denom 接近 0，说明 2×2 块本身奇异——这正是前面四级判据要极力避免的情况，所以走到这里 denom 一般不会出问题。

daxpy 被调用了两次——一次对第 $k$ 列、一次对第 $k-1$ 列。这是因为 2×2 块主元同时消去两列，每一列都要把它的倍数加到前面的所有列上去。两个 mulk 和 mulkm1 就是 2×2 块"逆"作用在前 $k-2$ 列上的结果。

整个 2×2 块处理的精妙之处在于：它不需要真的去求 2×2 块的特征值，不需要真的把矩阵对角化，只需要做一次"以非对角元为桥梁"的对称消元。这是 Bunch-Kaufman 相对于谱分解最大的工程优势——同样是处理不定矩阵，谱分解是 $O(n^3)$ 的迭代算法，Bunch-Kaufman 是 $O(n^3\mathrm{/}3)$ 的直接算法。

kpvt 数组：主元信息的"账本"

整个分解结束后，因子 $U$ 存在 $A$ 的下三角里（覆盖了原矩阵），而主元信息全部记在 kpvt[] 数组里。看这两段：

KP(k) = k;if (swap) KP(k) = imax;

（1×1 主元分支）

KP(k) = 1 - k;if (swap) KP(k) = -imax;KP(k-1) = KP(k);

（2×2 主元分支）

kpvt 的编码规则（这是 LINPACK 的 1-based 约定，sym.c 文件头明确保留了）：

• 正值：表示 1×1 主元。值就是真正的主元行号。KP(k) = k 表示没换过（用原对角元），KP(k) = imax 表示把 imax 行/列换到了 $k$ 位置。
• 负值：表示 2×2 块主元。绝对值是块的另一行号。KP(k) = 1 - k 表示块由 $k-1$ 和 $k$ 组成（没换），KP(k) = -imax 表示把 imax 换到了 $k-1$ 位置。注意 KP(k-1) = KP(k)——2×2 块占两列，两列的 kpvt 值相同，这是后面求解时识别"哪两列一起处理"的标志。

这个 kpvt 数组是后续求解（dsisl）的关键——求解时必须"反向应用"这些主元交换，才能正确地把右端项 $b$ 变换回去。没有 kpvt，分解因子 $U$ 就是一堆没有意义的数字。

工程细节：教科书从来不讲的几件事

细节一：Bunch 常数不是"经验值"，是解析最优

很多教材写到 $\alpha =(1+\sqrt{17})\mathrm{/}8$

​)/8 就打住了，最多加一句"由 Bunch 1971 论文给出"。但这个数到底怎么来的？

它是元素增长上界的极小值点。 Bunch 在论文里把每一步消元后矩阵元素的最大增长率写成 $\alpha$ 的函数 $g(\alpha )$，然后求 $g^{\mathrm{\prime }}(\alpha )=0$。这个方程化简后是 $\alpha$ 的一个二次关系，解出来正好是 $(1+\sqrt{17})\mathrm{/}8$

​)/8。

为什么是 17？因为 $17=1^2+16=1+4^2$，这是 2×2 块消元后增长率公式里系数的几何结构决定的。不是凑出来的，是算出来的。 任何一个学过数值线性代数的人，如果只记一个"工程常数"，那它一定是这个 $\alpha$——它是数值稳定性理论里少有的、能写出闭式最优解的例子。

细节二：每一步四个分支，对应四种稳定性情形

回头看 L20 到 L70 那段判据。四级判据不是冗余，每一种都对应一种稳定性情形：

1. 对角元本身就够大（L20 直接通过）：良态矩阵的常态，最快的路径。
2. 换一个对角元作主元（L30 通过）：类似部分主元法，找一个更稳的对角元。
3. 回到原对角元，1×1 凑合用（L60 通过）：原对角元比 2×2 块消元更稳，不必上 2×2。
4. 必须用 2×2 块（L70）：前面三种都不行，才上 2×2，吃下负特征值。

这四种情形的判据，全部源自 1971 年那篇论文的元素增长分析。 不是凭经验，不是启发式，是数学上的最优分类。这就是为什么 Bunch-Kaufman 50 多年来基本没被改动——它已经是最优的了。

细节三：2×2 块的就地处理——不开方、不对角化

看 L160 那段 2×2 块消元。注意一个细节：代码从头到尾没有对任何东西开方，也没有调用任何特征值函数。整个 2×2 块的处理，是通过"除以非对角元 $A(k-1,k)$"完成的。

为什么能这样？因为 2×2 对称块 的"作用"（把它从矩阵里消去）等价于做一个变换，这个变换的解析表达只涉及 $a,b,c$ 和 $b^2-ac$（也就是 $-$ 行列式）。只要 $b\mathrm{\ne }0$，就可以消元——不需要特征值，不需要开方。这正是前面四级判据保证的：走到 2×2 这一步，$b=A(k-1,k)$ 一定是非零的（否则判据会让它走 1×1 分支）。

这个设计极其精妙：它把"处理不定矩阵"这件事简化成了"多记一个非对角元"。整个算法的控制结构和 Cholesky 几乎一样，只是多了一个分支判断和 kpvt 记录。这就是为什么 LINPACK 当年选择 Bunch-Kaufman 而不是谱分解——它在工程上几乎是 Cholesky 的"无缝升级"。

细节四：kpvt 的反向应用——dsisl 里的两个分支

分解存在 $A$ 里、主元记在 kpvt 里，求解时怎么用？看 dsisl 的前代阶段：

L10:if (k == 0) goto L80;if (KP(k) < 0) goto L40;/* 2×2 块 */if (k != 1) {kp = KP(k);if (kp != k) { temp = B(k); B(k) = B(kp); B(kp) = temp; }daxpy(k - 1, B(k), &A(1,k), 1, &B(1), 1);}B(k) = B(k) / A(k,k);k = k - 1;goto L70;L40:if (k != 2) {kp = (KP(k) >= 0) ? KP(k) : -KP(k); /* IABS(KPVT(K)) */if (kp != k - 1) { temp = B(k-1); B(k-1) = B(kp); B(kp) = temp; }daxpy(k - 2, B(k), &A(1,k), 1, &B(1), 1);daxpy(k - 2, B(k-1), &A(1,k-1), 1, &B(1), 1);}...

if (KP(k) < 0) goto L40——通过 kpvt 的正负号区分 1×1 和 2×2 块。如果是 1×1，就单列处理；如果是 2×2，就两列一起处理，并且先把右端项 $b$ 的对应分量按主元交换的逆序换回去（temp = B(k-1); B(k-1) = B(kp); ...）。

这就是"反向应用主元"——分解时做了哪些行/列交换，求解时必须一模一样地在 $b$ 上做一遍，否则解就是错的。kpvt 数组就是这张"交换清单"，没有它，分解因子 $U$ 完全没法用。

注意 kp = (KP(k) >= 0) ? KP(k) : -KP(k)——这是手动取绝对值，因为 kpvt 的负值表示 2×2 块，但求解时我们需要的是"另一个主元行号"（正数）。代码用三元运算符实现了 Fortran 的 IABS，简洁但容易看漏。

与 Cholesky 的正面对比

把 sym.c 和 cholesky.c 放在一起看，差异一目了然：

最关键的一行对比：Cholesky 遇负就死，Bunch-Kaufman 把负吃下去。

// Cholesky (dpofa)if (s <= 0.0) return j + 1;// 非正定：直接失败// Bunch-Kaufman (dsifa)if (fabs(absakk) > 0.0 || fabs(colmax) > 0.0) goto L100;// 只要有非零元就继续KP(k) = k;info = k;// 只有"完全零列"才算奇异

dsifa 的失败条件是整列全零（absakk == 0 && colmax == 0），这比 Cholesky 的"对角元开方为负"要严格得多——只有真正奇异的矩阵才会失败。绝大多数"不正定"但"非奇异"的矩阵，Bunch-Kaufman 都能分解。这就是它存在的意义。

代价是计算量翻倍（$O(n^3\mathrm{/}3)$ vs $O(n^3\mathrm{/}6)$）和多了 kpvt 的开销。但比起"直接报错无法计算"，这个代价完全值得。

把整条链串起来：Bunch-Kaufman 在哪里默默运行

读懂 sym.c 这 200 行，你就读懂了统计软件处理"对称矩阵"的另一半底层：

Bunch-Kaufman 分解 (dsifa + dsisl) ├── 广义线性模型（GLM）：Hessian 矩阵往往不定，牛顿法每步都要分解 ├── 协方差矩阵求逆：实际数据里协方差矩阵未必正定（共线性、小样本） ├── 岭回归前的检测：判断矩阵是否需要正则化，先尝试 Bunch-Kaufman ├── 约束优化：拉格朗日 Hessian 在 saddle point 附近不定 ├── 因子分析：协方差结构分解，遇到病态样本时 Cholesky 失败 └── 时间序列：ARMA 模型的信息矩阵在某些参数区域不定

只要对称矩阵有可能不正定，用的就是 Bunch-Kaufman，不是 Cholesky。 这就是为什么 R 的 solve() 在面对对称矩阵时默认调 LAPACK 的 dsysv（Bunch-Kaufman 的现代版本），而不是 Cholesky——它更"鲁棒"，多付一倍计算量换"几乎不会失败"。

学完这一篇，你再看 R 的 glm()、Python 的 scipy.linalg.ldl()（注意，这就是 LDLᵀ，Bunch-Kaufman 的标准接口）、SPSS 的"非线性回归"对话框，会明白它们点下去的那一瞬间，跑的就是这 200 行——准确地说，是那个 (1.0 + sqrt(17.0)) / 8.0 和它后面那一长串 goto 跳转。

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

1. Cholesky 只能处理对称正定，Bunch-Kaufman 能处理任意对称。前提的放宽来自"允许 2×2 块主元"——这是教科书一句带过、但在代码里占了一半篇幅的核心机制。

2. Bunch 常数 $(1+\sqrt{17})\mathrm{/}8$

   ​)/8 是 1971 年论文的解析最优解，不是经验值。它最小化元素增长上界，化简后正好是 $\sqrt{17}$

   ​。全世界所有 Bunch-Kaufman 实现，这一行都长一样。

3. 每步四级判据，对应四种稳定性情形。从"对角元够大"到"必须上 2×2 块"，层层下探，每一级的判据都源自元素增长分析。不是启发式，是最优分类。

4. 2×2 块消元靠"非对角元"完成——不开方、不求特征值。只要 $A(k-1,k)\mathrm{\ne }0$，就能消元。这就是为什么不定矩阵也能算。

5. kpvt 数组记录主元信息，正值=1×1、负值=2×2。求解时必须"反向应用"这些交换，否则解全错。它是分解因子 $U$ 的"使用说明书"。

6. 失败条件是"整列全零"，不是"对角元为负"。这就是 Bunch-Kaufman 比 Cholesky 鲁棒的根本原因——只有真正奇异才失败，绝大多数"不正定但非奇异"的矩阵都能分解。

如果只记一句话，那就是：Cholesky 要求正定、用 $R^{T}R$；Bunch-Kaufman 接受不定、用 $UD{U}^T$，靠 2×2 块主元吃下负特征值。 那个看似古怪的 $(1+\sqrt{17})\mathrm{/}8$

​)/8，是 50 年前一篇论文算出来的最优阈值，今天每一份统计软件都还在用它，一字未改。

下一篇，我们会把镜头转向奇异值分解（SVD）——当矩阵不仅不定，甚至不是方阵、严重病态时，Bunch-Kaufman 也救不了，这时候就要上 SVD 了。那是另一套精妙绝伦的故事，涉及 Householder 双对角化和 Golub-Kahan 迭代，敬请期待。