余因子展開による行列式と逆行列を求めるコード

余因子展開で行列式と逆行列を求めるコードを書いた。
行列式の計算はそれ専用の関数を実装したが、逆行列はmain関数に書いた。C/C++だと2次元配列の引き渡しが面倒だからである。

• 逆行列を求める

int main(void)

    int i, j, _i, _j;
    // n次正方行列
    double **matrixA;
    // Determinant of matrix A
    double determinantMatrixA;
    // matrix dimension
    int n = 4;
    // n次正方行列の余因子展開で生じる(n-1)次正方行列の要素番号(indexRow:行, indexColumn:列)
    int indexRow, indexColumn;
    // n次正方行列の余因子展開で生じる(n-1)次正方行列の(n-1)を格納する変数
    int m;
    // n次正方行列の余因子展開で生じる(n-1)次正方行列
    double **matrixB;
    // (i, j)余因子の符号
    int sign;
    // 余因子
    double cofactor;
    // 余因子行列
    double **cofactorMatrixA;
    // 逆行列
    double **inverseMatrixA;
    // 単位行列
    double **identityMatrixE;
    double *baseAddress, *baseAddress2, *baseAddress3, *baseAddress4, *baseAddress5;

    // n次正方行列のメモリを確保する
    matrixA = new double*[n];
	baseAddress = new double[n*n];
    for(i = 0; i < n; i++){
        matrixA[i] = baseAddress + i*n;
    }

    // 余因子展開により生じるn個の行列を順に格納する(n-1)次正方行列のメモリを確保する
    m = n - 1;
    matrixB = new double*[m];
    baseAddress2 = new double[m*m];
    for(i = 0; i < m; i++){
        matrixB[i] = baseAddress2 + i*m;
    }

    // 余因子行列のメモリを確保する
    cofactorMatrixA = new double*[n];
    baseAddress3 = new double[n*n];
	for(i = 0; i < n; i++){
        cofactorMatrixA[i] = baseAddress3 + i*n;
    }

    // 逆行列のメモリを確保する
    inverseMatrixA = new double*[n];
    baseAddress4 = new double[n*n];
	for(i = 0; i < n; i++){
        inverseMatrixA[i] = baseAddress4 + i*n;
    }

    // 単位行列のメモリを確保する
    identityMatrixE = new double*[n];
    baseAddress5 = new double[n*n];
	for(i = 0; i < n; i++){
        identityMatrixE[i] = baseAddress5 + i*n;
    }

↑それぞれが何の配列かはコメントに書いた。

    matrixA[0][0] = 2.0;
    matrixA[1][0] = 4.0;
    matrixA[2][0] = -2.0;
    matrixA[3][0] = -6.0;
    matrixA[0][1] = -1.0;
    matrixA[1][1] = -1.0;
    matrixA[2][1] = 2.0;
    matrixA[3][1] = 5.0;
    matrixA[0][2] = 2.0;
    matrixA[1][2] = 6.0;
    matrixA[2][2] = 4.0;
    matrixA[3][2] = 3.0;
    matrixA[0][3] = 1.0;
    matrixA[1][3] = 3.0;
    matrixA[2][3] = 2.0;
    matrixA[3][3] = 9.0;

今回は4次正方行列の行列式(determinant)と逆行列を求める。
各要素の数字の格納の仕方はなんでもよい。上記ではハードコーディングで代入している。

    // 行列式の計算
    determinantMatrixA = calcDeterminant(&matrixA[0][0], n);

calcDeterminant関数で行列式を計算し、それを戻り値としている。
C/C++での2次元配列の引き渡しは面倒である。

    /* 2次元配列引き渡しによるメモリの扱いが大変なため
        逆行列の計算はmain関数内に記述する */
    // detA=|A|=0のときは逆行列が存在しないため計算を実行せずに終了とする
    if(-1.0e-10 < determinantMatrixA && determinantMatrixA < 1.0e-10){
        printf("cannot find inverse matrix\n");
        return(0);
    }

    // 2次正方行列は簡単に計算できる
    if(n == 2){
        inverseMatrixA[0][0] =  matrixA[1][1]/determinantMatrixA; inverseMatrixA[0][1] = -matrixA[0][1]/determinantMatrixA;
        inverseMatrixA[1][0] = -matrixA[1][0]/determinantMatrixA; inverseMatrixA[1][1] =  matrixA[0][1]/determinantMatrixA;
        return(0);
    }

    // 余因子行列の計算
    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){
            // i行j列(aij)のおける余因子の計算
            indexRow = 0;
            for(_i = 0; _i < n; _i++){
                if(_i == i){
                    continue;
                }
                indexColumn = 0;
                for(_j = 0; _j < n; _j++){
                    if(_j == j){
                        continue;
                    }
                    // i行j列目を除く
                    matrixB[indexRow][indexColumn] = matrixA[_i][_j]; 
                    indexColumn++;
                }
                indexRow++;
            }
            // (i, j)余因子(cofactor)を求める
            sign = ( (((i+1) + (j+1))%2 == 0) ? 1 : -1 );
            cofactor = (double)sign * calcDeterminant(&matrixB[0][0], m);
            // (i, j)余因子を余因子行列のj行i列目に代入する(転置)
            cofactorMatrixA[j][i] = cofactor;
        }
    }

    // 逆行列を計算する
    for(i = 0; i < n; i++){
        for(j = 0 ; j < n; j++){
            inverseMatrixA[i][j] = cofactorMatrixA[i][j]/determinantMatrixA;
        }
    }

上記の流れとしては

1. 行列式が0ならば逆行列の計算をしないでプログラムを終了する。
2. 2次正方行列は公式で解く。
3. n次正方行列(n>2)は余因子行列を計算する。
4. 余因子行列から逆行列を研鑽する。

// 行列Aとその逆行列A(-1)の積を計算し単位行列Eになることを確認する
    initializeMatrix(&identityMatrixE[0][0], n);
    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){
            for(int k = 0; k < n; k++){
                identityMatrixE[i][j] += matrixA[i][k] * inverseMatrixA[k][j];
            }
        }
    }

念のため、もとの行列とその逆行列を掛けて単位行列になるかを確かめる。
なお、このコードは見やすさためにprintfで行列の各要素を表示する部分は省略した。
↓実行結果
行列式も逆行列も計算できている。

-- Find the inverse of matrixA below. --
a(1,1)=2.00, a(1,2)=-1.00, a(1,3)=2.00, a(1,4)=1.00, 
a(2,1)=4.00, a(2,2)=-1.00, a(2,3)=6.00, a(2,4)=3.00, 
a(3,1)=-2.00, a(3,2)=2.00, a(3,3)=4.00, a(3,4)=2.00, 
a(4,1)=-6.00, a(4,2)=5.00, a(4,3)=3.00, a(4,4)=9.00, 

-- The determinant of matrix A = 60.00 --

-- The inverse matrix of A --
a(1,1)=-2.00, a(1,2)=1.00, a(1,3)=-0.50, a(1,4)=-0.00, 
a(2,1)=-3.50, a(2,2)=1.50, a(2,3)=-0.50, a(2,4)=0.00, 
a(3,1)=0.53, a(3,2)=-0.20, a(3,3)=0.33, a(3,4)=-0.07, 
a(4,1)=0.43, a(4,2)=-0.10, a(4,3)=-0.17, a(4,4)=0.13, 

-- The identity matrix, E=A*A(-1) --
a(1,1)=1.00, a(1,2)=-0.00, a(1,3)=-0.00, a(1,4)=0.00, 
a(2,1)=0.00, a(2,2)=1.00, a(2,3)=0.00, a(2,4)=0.00, 
a(3,1)=0.00, a(3,2)=-0.00, a(3,3)=1.00, a(3,4)=0.00, 
a(4,1)=0.00, a(4,2)=-0.00, a(4,3)=0.00, a(4,4)=1.00,

• 行列式を計算するcalcDeterminant関数
  コードので宣言されている変数の役割は上記のmain関数での同じ変数名と同じ。

double calcDeterminant(double *, int)

double calcDeterminant(double *matrixTmp, int n){

    int i, _i, j;
    int m;
    double **matrixA;
    double **matrixB;
    int indexRow;
    int sign;
    double cofactor;
    double determinant;
    double *baseAddress, *baseAddress2;

    // n次正方行列のメモリを確保する
    matrixA = new double*[n];
    baseAddress = new double[n*n];
	for(i = 0; i < n; i++){
        matrixA[i] = baseAddress + i*n;
    }
    // 1次元配列として引き渡された配列を2次元配列(n次正方行列)に代入する
    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){
            matrixA[i][j] = matrixTmp[i*n + j];
        }
    }

    // 余因子展開により生じるn個の行列を順に格納する(n-1)次正方行列のメモリを確保する。
    m = n - 1;
    matrixB = new double*[m];
    baseAddress2 = new double[m*m];
    for(i = 0; i < m; i++){
        matrixB[i] = baseAddress2 + i*m;
    }
    
    // 2次正方行列(n=2)ならばdeterminantを公式の通りに計算する
    if(n == 2){
        determinant = matrixA[0][0]*matrixA[1][1] - matrixA[0][1]*matrixA[1][0];
        return(determinant);
    }
    
    determinant = 0.0;
    // 第n列(j=n-1)での余因子展開をする
    for(i = 0; i < n; i++){
        initializeMatrix(&matrixB[0][0], m);
        // i行n列(ain)のおける余因子の計算
        indexRow = 0;
        for(_i = 0; _i < n; _i++){
            if(_i == i){
                continue;
            }
            for(j = 0; j < n-1; j++){
                // i行n列目を除く
                matrixB[indexRow][j] = matrixA[_i][j];
            }
            indexRow++;
        }
        // 余因子を求める
        sign = ( (((i+1) + n)%2 == 0) ? 1 : -1 );
        cofactor = (double)sign * calcDeterminant(&matrixB[0][0], m);
        // 余因子と(i, j)から行列式を求める
        determinant += matrixA[i][n-1] * cofactor;
    }

    // メモリ解放
    delete[] matrixA;
    delete[] matrixB;
    delete[] baseAddress;
    delete[] baseAddress2;

    return(determinant);

}

上記の流れとしては

1. 2次正方行列の行列式は公式で解く。
2. n次正方行列に対して第n列(j=n-1)で余因子展開をする。
3. 余因子展開のときに(n-1)正方行列が出てくるため、これをまた同じ関数に引き渡す(再帰)。簡単な公式で求められる2次正方行列になるまで再帰を繰り返す。
4. i行n列(0<=i<n, j=n-1)の余因子を計算する。
5. 求めた余因子から行列式を計算する。

• 実装の方法に関して
  冒頭でも述べたが、C/C++は2次元配列の引き渡しが面倒なので、逆行列の計算に関してはmain関数で実装した。
  行列式の計算と逆行列の計算は一部が似ている(余因子展開)ため、二つの分けずに一緒に求められないかと考えた。実装は出来るのかもしれないが、試したところややこしいと感じたため断念した。
  このコードはC言語だけでも書ける。

• その他
  逆行列の求め方は掃き出し法によりもっと短い行数で書けるが、そのサンプルコードを見てもよく分からなかった。
  そしてこの行列計算はMATLABやPythonで簡単にできる。
  LU分解とかで求められるみたいだが..

(随時更新予定)