行列式の計算

　行列式を効率的に計算する方法を考える。

まず、n*n行列Aの行列式は、

と定義される。ε()は、与えられた順列が奇順列の場合は-1、偶順列の場合は1となるような関数である。

　上の定義より、n*n行列の行列式をそのまま計算すると、O(n*n!)の計算量が必要となってしまう。大きなサイズ(n>20)の行列の行列式を簡単に計算するために、以下の性質を用いる。

行列Aのi行から、j行の定数倍を引いても、行列式の値は変わらない。

上記の性質を使って行列Aを上三角行列に変更する。行列式の定義を見ると、上三角行列の行列式は、対角成分の積と等しいことに気付く。

三角行列への変換はO(n^3)、対角成分の積を求める処理はO(n)なので、全体でO(n^3)の計算量で行列式を求めることができる。

以下サンプルソース。

const int N = 3;

double determinant(const double x[][N]) {
    double tmp[N][N];

    memcpy(tmp, x, sizeof(tmp));

    // make pivots non-zero numbers
    for (int i = 0; i < N; i++) {
        if (tmp[i][i] != 0)
            continue;
        int k;
        for (k = 0; k < N; k++)
            if (tmp[k][i] != 0)
                break;

        if (k == N)          // all numbers in a column is zero
            return 0;

        for (int j = 0; j < N; j++)
            tmp[i][j] += tmp[k][j];
    }

    // make a triangular matrix
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (i == j)
                continue;
            double c = tmp[j][i] / tmp[i][i];
            for (int k = 0; k < N; k++)
                tmp[j][k] -=  c * tmp[i][k];
        }
    }

    double det = 1;
    for (int i = 0; i < N; i++)
        det *= tmp[i][i];

    return det;
}

int main() {
    double x[N][N];

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            x[i][j] = -50 + rand() % 100;
    }

    cout << determinant(x) << endl;
    return 0;
}