概要
- 余因子展開で行列式と逆行列を求める。
計算量が非常に多いため、大きな行列の計算はメモリの容量や計算時間がネックになるかもしれない。MATLABやPythonのライブラリで計算した方が良い。 LU分解とかで求められるみたいだが..つまり今回のプログラムに有用性はない。 - C言語では関数への2次元配列の引き渡しが面倒であるため、逆行列の計算はmain関数に書いた。行列式の計算は再帰を実行する必要があるため関数にまとめた。
計算過程
- 元のn次正方行列、その逆行列などを格納する2次元配列は動的確保した。なお、関数への引き渡しが上手くいくように少し工夫をして動的確保をした(これに関しては本題からそれるため省略、別記事で書く予定)
- フローチャート (全体の概要)
- ソースコード
- このフローチャートと対応している部分は以下になる。(「行列の要素を具体的に定める(数値を代入)」過程は省略)
int i, j, _i, _j;
// n次正方行列
float **matrixA;
// Determinant of matrix A
float determinantMatrixA;
// matrix dimension
int n;
// n次正方行列の余因子展開で生じる(n-1)次正方行列の要素番号(indexRow:行, indexColumn:列)
int indexRow, indexColumn;
// n次正方行列の余因子展開で生じる正方行列の次元(=n-1)
int m;
// n次正方行列の余因子展開で生じる(n-1)次正方行列
float **matrixB;
// (i, j)余因子の符号
int sign;
// 余因子
float cofactor;
// 余因子行列
float **cofactorMatrixA;
// 逆行列
float **inverseMatrixA;
// 単位行列
float **identityMatrixE;
// 行列式の計算
determinantMatrixA = calcDeterminant(&matrixA[0][0], n);
printf("-- The determinant of matrix A = %.2f --\n\n", determinantMatrixA);
/* 2次元配列引き渡しによるメモリの扱いが大変なため
逆行列の計算はmain関数内に記述する */
if(fabs(determinantMatrixA) < 1.0e-10){
// (1) detA=|A|=0のときは逆行列が存在しないため計算を実行せずに終了とする
printf("Determinant = 0, Cannot find inverse matrix.\n");
}else{
// (2) detA=|A|=0でなければ逆行列を計算する
if(n == 1){
// (2-1) 1次正方行列の逆行列は定義通りに計算する
inverseMatrixA[0][0] = 1.0/determinantMatrixA;
}else if(n == 2){
// (2-2) 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][0]/determinantMatrixA;
}else{
// (2-3) 以下より行列式が0ではない3次以上の正方行列の逆行列を計算する
// (インデントがすごいことになっている)
// 余因子行列を求める
for(i = 0; i < n; i++){
for(j = 0; j < n; j++){
// i行j列(a(i,j))における余因子の計算
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 = (float)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;
}
}
// 逆行列の計算ここまで
}
// 3次以上の逆行列の計算はここまで
printf("-- The inverse matrix of A --\n");
printMatrix(&inverseMatrixA[0][0], n);
// 行列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("-- The identity matrix, E=A*A(-1) --\n");
printMatrix(&identityMatrixE[0][0], n);
}
- 行列式計算(calcDeterminant関数)
- 行列式も余因子展開で計算しているが、
逆行列はすべての(i, j)の余因子を求めているのに対して、
determinamtはn列目(プログラム的にはj=n-1)のみ余因子展開をしているため、main関数内での逆行列の計算よりは少し簡単となる。
- 行列式も余因子展開で計算しているが、
- calcDeterminant関数のソースコード
- (本当はフローチャートも書きたいのだが、再帰処理や余因子展開の細かい流れを書くのが非常に大変そうで断念した..)
- 下記ソースコードにある2次元配列matrixBはn次正方行列matrixAを余因子展開するときに生じる(n-1)次正方行列を格納するためのものである。n次正方行列を余因子展開するときは(n-1)次正方行列の行列式を求める必要がある
int i, _i, j;
int m;
float **matrixA;
float **matrixB;
int indexRow;
int sign;
float cofactor;
float determinant;
if(n == 1){
// 1次正方行列(n=1)ならば値をそのまま返す
determinant = matrixA[0][0];
}else if(n == 2){
// 2次正方行列(n=2)ならばdeterminantを公式の通りに計算する
determinant = matrixA[0][0]*matrixA[1][1] - matrixA[0][1]*matrixA[1][0];
}else{
// 3次以上の正方行列の行列式を計算する
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 = (float)sign * calcDeterminant(&matrixB[0][0], m);
determinant += matrixA[i][n-1] * cofactor;
}
}
実際の計算
A=\begin{pmatrix}
2.0 & -1.0 & 2.0 & 1.0 \\
4.0 & -1.0 & 6.0 & 3.0 \\
-2.0 & 2.0 & 4.0 & 2.0 \\
-6.0 & 5.0 & 3.0 & 9.0
\end{pmatrix}
上記の4次正方行列の行列式と逆行列を求める。
- 計算結果
- 特に問題はなさそう
-- 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,
\begin{vmatrix}
A
\end{vmatrix}=60.00
A^{-1}=\begin{pmatrix}
-2.00 & 1.00 & -0.50 & -0.00 \\
-3.50 & 1.50 & -0.50 & 0.00 \\
0.53 & -0.20 & 0.33 & -0.07 \\
0.43 & -0.10 & -0.17 & 0.13
\end{pmatrix}
あとがき
冒頭でも申し上げた通り、このプログラムを書くのにすごく時間をかけたが、見返りは特にない。もともと余因子展開で行列式や逆行列をもとめる作業は4次あたりで手計算をしたことはあるし、とくに新しい知見を得られたわけではない。非常に要領の悪いことをしたことになる。