Original

Kraut 방법을 이용한 $O(N^3)$ 행렬식 계산

이 글에서는 $O(N^3)$에 동작하는 Kraut 방법을 이용하여 행렬의 행렬식을 구하는 방법을 설명합니다.

Kraut 알고리즘은 행렬 $A$를 $A = L U$ 꼴로 분해하며, 여기서 $L$은 하삼각 행렬이고 $U$는 상삼각 행렬입니다. 일반성을 잃지 않고 $L$의 대각 원소가 모두 1이라고 가정할 수 있습니다. 이 행렬들을 알고 나면 $A$의 행렬식은 쉽게 계산할 수 있습니다. 행렬식은 행렬 $U$의 주대각선 위에 있는 모든 원소의 곱과 같습니다.

가역 행렬은 모든 주소행렬식이 0이 아닐 때, 그리고 오직 그때에만 LU 분해를 가지며 그 분해가 유일하다는 정리가 있습니다. 여기서는 행렬 $L$의 대각선이 모두 1로 이루어진 분해만을 고려합니다.

$A$를 행렬, $N$을 그 크기라고 합시다. 다음 단계에 따라 행렬 $L$과 $U$의 원소를 구합니다:

1. $i = 1, 2, ..., N$에 대해 $L_{i i} = 1$로 둡니다.

2. 각 $j = 1, 2, ..., N$에 대해 다음을 수행합니다:

   • $i = 1, 2, ..., j$에 대해 다음 값을 구합니다

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

   • 다음으로, $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;
    }
}