行列式の計算

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

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

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

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

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

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

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

以下サンプルソース。

view plainprint?

1. const int N = 3;  
2.   
3. double determinant(const double x[][N]) {  
4.     double tmp[N][N];  
5.   
6.     memcpy(tmp, x, sizeof(tmp));  
7.   
8.     // make pivots non-zero numbers  
9.     for (int i = 0; i < N; i++) {  
10.         if (tmp[i][i] != 0)  
11.             continue;  
12.         int k;  
13.         for (k = 0; k < N; k++)  
14.             if (tmp[k][i] != 0)  
15.                 break;  
16.   
17.         if (k == N)          // all numbers in a column is zero  
18.             return 0;  
19.   
20.         for (int j = 0; j < N; j++)  
21.             tmp[i][j] += tmp[k][j];  
22.     }  
23.   
24.     // make a triangular matrix  
25.     for (int i = 0; i < N; i++) {  
26.         for (int j = 0; j < N; j++) {  
27.             if (i == j)  
28.                 continue;  
29.             double c = tmp[j][i] / tmp[i][i];  
30.             for (int k = 0; k < N; k++)  
31.                 tmp[j][k] -=  c * tmp[i][k];  
32.         }  
33.     }  
34.   
35.     double det = 1;  
36.     for (int i = 0; i < N; i++)  
37.         det *= tmp[i][i];  
38.   
39.     return det;  
40. }  
41.   
42. int main() {  
43.     double x[N][N];  
44.   
45.     for (int i = 0; i < N; i++) {  
46.         for (int j = 0; j < N; j++)  
47.             x[i][j] = -50 + rand() % 100;  
48.     }  
49.   
50.     cout << determinant(x) << endl;  
51.     return 0;  
52. }  

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