cp-algorithms 指南：Kraut 方法 O(N³) 计算矩阵行列式（LU 分解法）

本篇技术指南基于 cp-algorithms 开源仓库中的 determinant-kraut.md 文章，深入讲解如何利用 Kraut 方法在 O(N³) 时间内计算方阵的行列式。Kraut 方法的核心思想是先将矩阵分解为下三角矩阵 L 与上三角矩阵 U 的乘积（LU 分解），再利用「行列式等于 U 矩阵主对角线元素之积」这一性质快速得出结果。读完本文，你将掌握 LU 分解的完整推导过程、原地计算的实现细节，以及一套带主元选取与高精度 BigDecimal 支持的 Java 参考实现，可直接用于竞赛编程与工程实践。

方法总览：从 LU 分解到行列式

给定一个大小为 N×N 的方阵 A，Kraut 方法并不直接枚举排列计算行列式（那将是指数级复杂度），而是先将 A 分解为：

$$A = L \cdot U$$

其中：

• L 是下三角矩阵（lower triangular），对角线以下有值、对角线以上全为零；
• U 是上三角矩阵（upper triangular），对角线以上有值、对角线以下全为零。

不失一般性，可以约定 L 的所有主对角线元素都等于 1。这一约定消除了分解的非唯一性，使 L 与 U 能被唯一确定。一旦得到这两个矩阵，行列式的计算就变得极为简单：

$$\det(A) = \prod_{i=1}^{N} U_{ii}$$

即 A 的行列式等于 U 矩阵主对角线全部元素的乘积。这是因为 $\det(A) = \det(L) \cdot \det(U)$，而 L 是对角线全为 1 的下三角矩阵，其行列式恒等于 1。

存在性与唯一性定理

关于 LU 分解，线性代数中有如下定理：任何一个可逆（invertible）矩阵都存在 LU 分解，并且当且仅当矩阵的所有主子式（principle minors）均非零时，该分解是唯一的。本文讨论的正是这种「L 对角线全为 1」意义下的唯一分解。

需要特别强调的是：当矩阵退化（奇异）时，Kraut 过程会遇到除零或无法继续的环节，因此实现中必须显式检测这一情况并返回 0。这一点与仓库中同目录的 Gauss 求行列式方法 的处理方式一致——该方法同样规定「若在当前列找不到非零元素，则算法立即停止并返回 0」。

算法推导：逐步计算 L 与 U 的元素

设 A 为 N×N 矩阵。Kraut 算法按下述步骤逐列求出 L 与 U 的全部元素：

第一步：对 $i = 1, 2, \dots, N$，令

$$L_{ii} = 1$$

第二步：对每一列 $j = 1, 2, \dots, N$，依次执行两个子步骤：

1. 对 $i = 1, 2, \dots, j$，计算 U 的第 j 列元素：

$$U_{ij} = A_{ij} - \sum_{k=1}^{i-1} L_{ik} \cdot U_{kj}$$

2. 对 $i = j+1, j+2, \dots, N$，计算 L 的第 j 列元素：

$$L_{ij} = \frac{1}{U_{jj}} \left( A_{ij} - \sum_{k=1}^{j-1} L_{ik} \cdot U_{kj} \right)$$

正确性直观

观察以上公式可以发现，计算 $U_{ij}$ 时用到的求和项只涉及「当前行左侧」已求出的 L 元素和「当前列上方」已求出的 U 元素；计算 $L_{ij}$ 时同理。这正是按列主序推进的原因：处理第 j 列时，L 的前 j 列与 U 的第 1..j 行均已就绪，所有需要的乘积项 $L_{ik} \cdot U_{kj}$ 都已具备，因此可以无回溯地一次性完成整列的求解。整个算法的总复杂度为：

$$\underbrace{O(N^2)}_{\text{每列计算}} \times \underbrace{O(N)}_{\text{每项求和}} = O(N^3)$$

Java 参考实现解析

原文档给出了一份基于 java.math.BigDecimal 的高精度 Java 实现。该实现不仅完成 LU 分解，还额外处理了三件对数值稳定性至关重要的事情：零行检测、隐式主元选取（scaled partial pivoting）、行交换的符号跟踪。下面结合源码逐一剖析。

完整代码

static BigInteger det (BigDecimal a [][], int n) {
	try {

	for (int i=0; i<n; i++) {
		boolean nonzero = false;
		for (int j=0; j<n; j++)
			if (a[i][j].compareTo (new BigDecimal (BigInteger.ZERO)) > 0)
				nonzero = true;
		if (!nonzero)
			return BigInteger.ZERO;
	}

	BigDecimal scaling [] = new BigDecimal [n];
	for (int i=0; i<n; i++) {
		BigDecimal big = new BigDecimal (BigInteger.ZERO);
		for (int j=0; j<n; j++)
			if (a[i][j].abs().compareTo (big) > 0)
				big = a[i][j].abs();
		scaling[i] = (new BigDecimal (BigInteger.ONE)) .divide
			(big, 100, BigDecimal.ROUND_HALF_EVEN);
	}

	int sign = 1;

	for (int j=0; j<n; j++) {
		for (int i=0; i<j; i++) {
			BigDecimal sum = a[i][j];
			for (int k=0; k<i; k++)
				sum = sum.subtract (a[i][k].multiply (a[k][j]));
			a[i][j] = sum;
		}

		BigDecimal big = new BigDecimal (BigInteger.ZERO);
		int imax = -1;
		for (int i=j; i<n; i++) {
			BigDecimal sum = a[i][j];
			for (int k=0; k<j; k++)
				sum = sum.subtract (a[i][k].multiply (a[k][j]));
			a[i][j] = sum;
			BigDecimal cur = sum.abs();
			cur = cur.multiply (scaling[i]);
			if (cur.compareTo (big) >= 0) {
				big = cur;
				imax = i;
			}
		}

		if (j != imax) {
			for (int k=0; k<n; k++) {
				BigDecimal t = a[j][k];
				a[j][k] = a[imax][k];
				a[imax][k] = t;
			}

			BigDecimal t = scaling[imax];
			scaling[imax] = scaling[j];
			scaling[j] = t;

			sign = -sign;
		}

		if (j != n-1)
			for (int i=j+1; i<n; i++)
				a[i][j] = a[i][j].divide
					(a[j][j], 100, BigDecimal.ROUND_HALF_EVEN);

	}

	BigDecimal result = new BigDecimal (1);
	if (sign == -1)
		result = result.negate();
	for (int i=0; i<n; i++)
		result = result.multiply (a[i][i]);

	return result.divide
		(BigDecimal.valueOf(1), 0, BigDecimal.ROUND_HALF_EVEN).toBigInteger();
	}
	catch (Exception e) {
		return BigInteger.ZERO;
	}
}

三个关键设计点

1. 零行预检测（防止错误早退）

函数开头对每一行做检查：若某一行全部为零（compareTo 均不大于 0），则矩阵必然奇异，直接返回 BigInteger.ZERO。注意该判断使用的是「> 0」而非「!= 0」——对于采用浮点或高精度十进制的矩阵，这一写法会连带判定「全为负数的行」为非零，属于原实现的既定语义，使用时需按输入数据的符号分布理解其行为。

2. 隐式主元选取（scaled partial pivoting）

为提高数值稳定性，实现先为每一行计算缩放因子 scaling<a href="https://link.gitcode.com/i/96cfa1a400234f038cd768186236b49d" target="_blank">i]（等于该行最大绝对值的倒数，除法精度设为 100 位小数），然后在第 j 列的主元搜索阶段，对每个候选元素按「绝对值 × 所在行缩放因子」的比较结果选取最大值所在行 imax。这种「先按行内量级归一化再比较」的策略与 [Gauss-Jordan 消元 中关于主元选取启发式的论述一致——该文指出：若某一行的方程被整体乘以 $10^6$，未归一化的朴素选主元几乎必然选中它，导致数值稳定性下降，因此隐式主元选取（implicit pivoting）是更稳妥的选择。

3. 行交换与符号跟踪

选完主元后，若 j != imax，代码将两行完整交换（同时交换缩放因子），并将 sign 取反。行交换会改变行列式的符号——每交换一次，行列式乘上 −1。最终结果计算时先累乘主对角线 $a[i][i]$（分解完成后对角线即 U 的主对角元素），再根据 sign 决定是否取负。

4. 原地存储与精度控制

整个分解**原地（in-place）**进行：L 的下三角部分（含对角线上的 1）与 U 的上三角部分（含对角线）共用矩阵 a，每列算完后，属于 L 的 $L_{ij}$ 立即除以 $U_{jj}$ 完成归一化。所有除法都使用 ROUND_HALF_EVEN（银行家舍入）且小数位设为 100，最后把乘积结果舍入到整数再转成 BigInteger，从而把浮点误差控制在可接受范围。异常（如除零、BigDecimal 溢出）统一被 catch 捕获并返回 0。

与 Gauss 方法的对比及选用建议

仓库的 线性代数模块 中提供了两条求行列式的路线，它们在 navigation.md 中并列展示：

维度	Kraut 方法（本文）	Gauss 消元法
源文档	determinant-kraut.md	determinant-gauss.md
核心思想	LU 分解，行列式 = U 对角线乘积	初等行变换化为对角矩阵，行列式 = 对角线乘积（跟踪符号）
复杂度	O(N³)	O(N³)
数值策略	隐式（按行缩放）选主元，BigDecimal 高精度	朴素列选主元 + 行交换，double + EPS 判断
输出类型	BigInteger（需精确整数结果时更稳）	double（浮点）
适用场景	需要高精度整数行列式、数值较大易溢出时	竞赛中常规浮点矩阵、对精度要求不高时

从代码结构看，Kraut 实现的一大优势在于：它把 LU 分解与行列式计算合二为一，若后续还需用同一分解求解线性方程组，可直接复用分解结果（即 L 与 U 因子），这是纯 Gauss 消元版本不具备的附加值。

在仓库中的定位与验证方式

该文章属于仓库「Linear Algebra → Matrices」栏目下的第三篇文章，与 Gauss 求行列式、高斯解线性方程组、矩阵秩 同属矩阵运算专题。仓库本身是 cp-algorithms.com 站的 Markdown 源码库（基于 e-maxx.ru 经典资源翻译并扩展），所有文章统一以 src/<分类>/<主题>.md 形式存放。

仓库的 test/ 目录下配套了完整的代码验证体系（extract_snippets.py 负责从文章中抽取代码块，test.sh 使用 g++ 与 C++17 标准编译并运行各用例）。虽然当前 determinant-kraut.md 的 Java 代码块未被打包进该 C++ 测试流水线，但相邻的 Gauss 方法有对应的 test_gauss.cpp 覆盖「唯一解、无穷多解、无解」三类场景。若读者需要自行验证本文 Java 实现，建议构造已知行列式的小矩阵（如 2×2 的 [[1,2],[3,4]]，行列式为 −2）逐例比对结果。

总结

Kraut 方法通过 LU 分解把 O(N!) 的直接行列式枚举降到 O(N³)，其关键在于：

• 分解结构：$A = LU$，L 对角线全为 1，分解唯一性由主子式非零保证；
• 计算顺序：按列推进，第 j 列先算 U 前 j 行、再算 L 下方各行，全部依赖都已就绪；
• 求值收尾：$\det(A) = \prod U_{ii}$，配合行交换符号跟踪与隐式选主元，兼顾正确性与数值稳定性。

对于需要在高精度场景下（如结果必须是精确整数、元素数量级巨大）求行列式的竞赛题目，本文给出的 BigDecimal 版本是可靠的参考实现；对于常规浮点矩阵，可选用仓库中更简洁的 Gauss 版本。两种方法共同构成了 cp-algorithms 在矩阵行列式计算上的完整工具箱。

延伸阅读

• Gauss 消元法求行列式：O(N³) 浮点实现的另一种选择
• 高斯消元解线性方程组：主元选取启发式的理论背景与隐式主元的由来
• 矩阵的秩：与行列式紧密相关的矩阵基本量
• 测试框架：仓库如何编译并验证文章中抽取的代码片段