克劳特算法计算行列式（Kraut）

使用克劳特算法以 $O(N^3)$ 的复杂度计算行列式（Calculating the determinant using Kraut method in $O(N^3)$）

在本文中，我们将描述如何使用 Kraut 方法找到矩阵的行列式，该方法的时间复杂度为 $O(N^3)$ 。

克劳特算法发现矩阵 $A$ 的分解为 $A=LU$ ，其中 $L$ 是下三角矩阵， $U$ 是上三角矩阵。在不丧失一般性的情况下，我们可以假设 $L$ 的所有对角线元素都等于 $1$ 。一旦我们知道这些矩阵，就很容易计算出 $A$ 的行列式：它等于矩阵 $U$ 的主对角线上所有元素的乘积。

有一个定理表明，任何可逆矩阵都有 $LU$ 分解，并且它是唯一的，当且仅当它的所有主余子式都是非零的。我们只考虑这样的分解，其中矩阵 $L$ 的对角线是由矩阵的对角线构成的。

设 $A$ 为矩阵， $N$ 为矩阵的大小。我们将使用以下步骤找到矩阵的 $L$ 和 $U$ 元素：

1. 令 $L_{ii}=1$ for $i=1,2,...,N$
2. 对于每个 $j=1,2,...,N$ ，执行：
   • for $i=1,2,...,j$ 寻找值
     $U_{ij}=A_{ij}-{\sum }_{k=1}^{i-1}{L}_{ik}\cdot U_{kj}$
   • 接下来，for $i=j+1,j+2,...,N$ 寻找值
     $L_{ij}=\frac{1}{U_{jj}}\left(A_{ij}-{\sum }_{k=1}^{j-1}{L}_{ik}\cdot U_{kj}\right)$

---

实现

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;
	}
}