このページでは、C言語での「逆行列の求め方」について解説していきます。
逆行列を求める方法はいくつかあると思いますが、今回は「掃き出し法」（より正確にはガウスジョルダン法）により逆行列を求めていきたいと思います。
この方法を利用すれば、4×4 や 4×4 といった高次の逆行列も比較的簡単に求めることが可能です！
C言語での行列の扱い方については下記ページで解説していますので、そもそもC言語で行列を扱う方法をご存知ない方は、事前に下記ページを読んでいただくことをオススメします。
【C言語】行列の扱い方
Contents
掃き出し法での逆行列の求め方
最初に掃き出し法で逆行列を求める際の考え方をおさらいし、続いて手順について解説していきたいと思います。
基本的には実装に必要になる手順のみを解説していきますので、「掃き出し法で逆行列を求められる理由」等を知りたい方は、別途「掃き出し法 逆行列 証明」などでググって調べていただければと思います。
掃き出し法による逆行列の算出
まずは、掃き出し法で下記の３次の正方行列 \(A\) の逆行列 \(A^{-1}\) を求めることを例に考えていきましょう。
$$ A = \left ( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right ) $$
掃き出し法においては、まず最初に逆行列を求めたい正方行列の右側に単位行列を追加して拡張した行列を用意します。
$$ \left ( \begin{array}{cccccc} a_{11} & a_{12} & a_{13} & 1 & 0 & 0 \\ a_{21} & a_{22} & a_{23} & 0 & 1 & 0 \\ a_{31} & a_{32} & a_{33} & 0 & 0 & 1 \end{array} \right ) $$
さらに、掃き出し法においては、上記の拡張した行列に次の３つの操作（行基本変形）を繰り返し行うことで、
- （１）：２つの行を入れ替える
- （２）：ある行を定数倍する（０倍以外）
- （３）：ある行に他の行を定数倍したものを加える
下記のように、行列の左側が単位行列になるように行列の成分を変化させていきます。
$$ \left ( \begin{array}{cccccc} 1 & 0 & 0 & x_{11} & x_{12} & x_{13} \\ 0 & 1 & 0 & x_{21} & x_{22} & x_{23} \\ 0 & 0 & 1 & x_{31} & x_{32} & x_{33} \end{array} \right ) $$
そして、この状態になった際の、元々単位行列だった右側の行列部分が、元々行列の左側にあった行列の逆行列となります。
つまり、上記の場合であれば、次の行列が正方行列 \(A\) の逆行列 \(A^{-1}\) ということになりますね！
$$ A^{-1} = \left ( \begin{array}{ccc} x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{array} \right ) $$
掃き出し法で逆行列を求める場合、どちらかというと上記の逆行列 \(A^{-1}\) における成分 \(x_{ij}\) を求めるというよりかは、用意した行列の左側が単位行列になるまで３つの操作を繰り返せば、自然に \(x_{ij}\) が求まっているという感覚に近いと思います。
スポンサーリンク
掃き出し法で逆行列を求める手順
では、続いては、行列の左側を単位行列にするための手順について解説していきます。
前述の通り、掃き出し法で逆行列を求める際に行う操作は下記の３つのみとなります。
- （１）：２つの行を入れ替える
- （２）：ある行を定数倍する（０倍以外）
- （３）：ある行に他の行を定数倍したものを加える
で、掃き出し法では、この３つの操作を繰り返し行なって「行列の左側部分を単位行列に変形する」ことで逆行列を求めていくことになります。
上記のような行列の変形は、対角成分に注目し、その対角成分が存在する列ごとに単位行列に近づけていくことで実現することができます。
対角成分に注目する
もう少し具体的に手順を解説していきます！
まずは、１つの対角成分に注目します。最初に注目するのは１番左の列の対角成分になります。
注目している対角成分を 1 にする
次に、この注目している対角成分を 1 にすることを考えます。
この時一番手っ取り早いのは、対角成分自身で割ってやることになりますね！
ただ、１つの成分に対してのみ計算を行う操作は上記の３つの操作にはありません。そのため、ここでは「注目している対角成分が存在する行」を定数倍することで対角成分を 1 にします。より具体的には、その行を 1 / 対角成分 倍します。
これは３つの操作のうちの（２）の操作にあたります。
これにより必ず注目している対角成分が 1 になります。
上記の３つの操作は全て、行全体に対する操作であることに注意してください。逆行列を求めたい左側の行列に対してのみ行うのではなく、追加した右側の行列に対しても同じ操作を行う必要があります。
残りの列の成分を 0 にする
さらに、次は「注目している対角成分が存在する列」の残りの成分を 0 にしていくことを考えます（その列の対角成分以外を 0 にする）。現状は、一番左の列の対角成分に注目しているのですから、一番左の列を 0 にしていきます（ここからは図中の行列の成分を簡略化して書いていきます）。
ここでポイントになる操作が３つの操作の（３）です。掃き出し法においては「ある行に他の行を定数倍したものを加える」操作が可能です。
ここでは、この「注目している対角成分が存在する行」を定数倍したものを他の行に加えていくことで、他の行の「注目している対角成分が存在する列」の成分を 0 にしていきます。
注目している対角成分は操作（２）によって 1 になっていますので、定数倍する時に用いる倍数として「注目している対角成分が存在する列」の成分の符号を逆転した値を用いれば、２つの行を足し合わせた時に「注目している対角成分が存在する列」の成分を 0 にすることができます。
で、これを「注目している対角成分が存在する行」以外の全ての行に対して繰り返し行えば、「注目している対角成分が存在する列」の各成分は下記のように変化することになります。
- 対角成分：
1 - それ以外：
0
つまり、「現在注目している対角成分が存在する列」に関しては、対角成分のみを 1 に、それ以外を 0 にすることができたことになります。
対各成分を移動しながら同様の処理を繰り返す
あとは、注目する対角成分を移動しながら同様の処理を行えば良いことになります。
まず、注目する対角成分を次の列のものに移し、
次に（２）の操作により注目している対角成分を 1 にし、
さらに（３）の操作によって「注目している対角成分が存在する列」の成分を全て 0 にします（対角成分を除いて）。
そうすれば、「今まで注目していた対角成分が存在する列」全てを、対角成分のみを 1 に、それ以外を 0 にすることができます。
そして、この一連の処理を、全ての列の対角成分に対して繰り返し行えば、いずれは行列の左側が全て下記の状態を満たすことになり、行列の左側が単位行列になったことになります。
- 対角成分：
1 - それ以外：
0
最初にも説明した通り、この状態になった時の行列の右側部分が、元々行列の左側にあった行列の逆行列ということになります。ですので、この時点で逆行列を求めるための処理は完了です。
ここまで図としては「３次の正方行列の逆行列を求める例」のものを示してきましたので、対角成分の個数は３つになります。したがって、３回分上記のような処理を繰り返してやれば、逆行列が求まることになります。
もちろん他の次数の正方行列であっても同様の処理で逆行列を求めることができ、次数に応じた回数分上記のような処理を繰り返すことになります。
ここまでの操作のみではゼロ除算が発生する
ただ、１点注意点があって、ここまで解説してきた手順で逆行列が求められるのは、「注目する対角成分が全て 0 以外である場合のみ」となります。
これは、注目している対角成分が 0 の時には 1 / 対角成分 の計算時にゼロ除算が発生し、そのままでは（２）の操作による対角成分を 1 にする計算が上手く行えないからです。
これを解決するために、ここまで使用しなかった（１）の操作「２つの行を入れ替える」が活躍することになります。
ただ、ここまでの解説内容のみを実装したとしても対角成分が 0 にならなければ逆行列を求めることは可能なので、一旦ここまでの流れを整理するためにも、まずはここまで解説してきた内容の実装について解説していきたいと思います。
そしてその後に、対角成分が 0 になった時に「ゼロ除算を回避する方法」について解説していきたいと思います。
掃き出し法の実装
では、ここまで解説してきた内容をC言語で実装する流れを解説していきます。
２次元配列の変数宣言と成分の格納
まず行列を扱うので、下記のページで解説している通り２次元配列の変数宣言を行います。
【C言語】行列の扱い方
掃き出し法で逆行列を求める際に必要な２次元配列
ここで変数宣言を行う必要のある２次元配列は下記の３つになります。
mat：逆行列を求める対象となる正方行列を扱う２次元配列inv：逆行列を扱う２次元配列sweep：掃き出し法の操作を行う２次元配列
２次元配列の変数宣言時のポイント
これらを変数宣言する際のポイントもいくつかあって、まず mat の逆行列を求めるわけですから、mat は正方行列である必要があります。つまり、行列の行数と列数が同じようになるように変数宣言してやる必要があります。
また、inv に関しては mat の逆行列となるわけですから、mat と同じ行数・列数である必要があります。
さらに、sweep は初期状態では、掃き出し法での逆行列を求める際の考え方 で紹介した下記のような、左側に逆行列を求めたい正方行列の成分、右側に単位行列の成分を持った行列を扱う必要があります。
$$ \left ( \begin{array}{cccccc} a_{11} & a_{12} & a_{13} & 1 & 0 & 0 \\ a_{21} & a_{22} & a_{23} & 0 & 1 & 0 \\ a_{31} & a_{32} & a_{33} & 0 & 0 & 1 \end{array} \right ) $$
で、この sweep に 掃き出し法での逆行列を求める際の考え方 で紹介した３つの操作を行うことで逆行列を求めていくことになります。
今回逆行列を求めたい正方行列は mat であり、sweep には mat に加えて右側の単位行列の成分も格納できるようにする必要があるため、行数は mat と同じ、列数は mat の２倍となるように変数宣言する必要があります。
また、逆行列を求める際には小数点以下の値も扱う必要があるため、sweep と inv に関しては小数点以下も扱える型（float や double 等）として宣言する必要があります（mat に関しては扱う値による）。
以上を考慮すると、掃き出し法で逆行列を求めていく際には、まず下記のような変数宣言を行う必要があることになります。
/* 逆行列を求める行列用の２次元配列 */
double mat[N][N];
/* 逆行列用の２次元配列 */
double inv[N][N];
/* 掃き出し法を行う行列 */
double sweep[N][N * 2];N は #define で定義してやる必要があり、N を変更することでさまざまな次数の正方行列の逆行列を求めることができます。
また、mat は小数点以下の値も扱えるよう、double 型として宣言しています。
sweep への成分の格納
続いて、掃き出し法を行なっていくための準備として、sweep への行列の成分の格納を行なっていきます。
前述の通り、sweep の左側には mat の成分が、右側には単位行列の成分が格納されるようにしていく必要があります。
これは、逆行列を求める対象となる行列 mat は行数列数ともに N なので、sweep の第 0 列から第 N - 1 列に mat の成分を、sweep の第 N 列から第 2 * N - 1 列に単位行列の成分を格納していくことで実現することができます。
具体的には、下記の処理により、掃き出し法を行なっていく上で必要になる行列 sweep を準備することができます。
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
/* sweepの左側に逆行列を求める行列をセット */
sweep[i][j] = mat[i][j];
/* sweepの右側に単位行列をセット */
sweep[i][N + j] = (i == j) ? 1 : 0;
}
}スポンサーリンク
操作（２）で対角成分を 1 にする
ここまでで掃き出し法を行うための準備が整ったことになりますので、次は掃き出し法で逆行列を求めていきます。
基本的に実装すべきことは 掃き出し法での逆行列を求める際の考え方 で書いた通りで、まず対角成分を 1 にするために操作（２）を行います。
ここで、現在注目している対角成分が存在する列を k とすれば、この操作は、1 / sweep[k][k]を求め、k 行目の全成分に 1 / sweep[k][k] を掛けることで実現することができます。
/* sweep[k][k]に掛けると1になる値を求める */
a = 1 / sweep[k][k];
/* 操作（２）：k行目をa倍する */
for (j = 0; j < N * 2; j++) {
/* これによりsweep[k][k]が1になる */
sweep[k][j] *= a;
}これにより、sweep の k 行目の成分、すなわち sweep[k][0] 〜 sweep[k][2 * N - 1] の値が変化し、さらに sweep[k][k] が 1 となります。
操作（３）で対角成分の列を 0 にする
続いて、k 行目以外の k 列目の成分が全て 0 になるように、操作（３）を行なっていきます。
操作（２）により、sweep[k][k] が 1 となっていますので、例えば i 行目の k 行目の成分、すなわち sweep[i][k] を 0 にするためには、k 行目を -sweep[i][k]倍し、その結果を i 行目に足してやれば良いことになります。
したがって、下記のような処理により、i 行目の k 行目の成分を 0 にすることができます。
/* k行目に掛ける値を求める */
a = -sweep[i][k];
for (j = 0; j < N * 2; j++) {
/* i行目にk行目をa倍した行を足す */
/* これによりsweep[i][k]が0になる */
sweep[i][j] += sweep[k][j] * a;
}で、ここで行いたいのは「k 行目以外の k 列目の成分を全て 0 にする」ことですので、行の添字となる i を 0 から N - 1 まで変化させながら上記の処理を繰り返す必要があります。
ただし、i が k の場合に同様のことを行うと対角成分 sweep[k][k] まで 0 になってしまうため、i が k の場合は上記の処理を行わないように制御する必要があります。
具体的には、次のように i に対するループの中で上記の処理を行うようにします。
/* 操作（３）によりk行目以外の行のk列目を0にする */
for (i = 0; i < N; i++) {
if (i == k) {
/* k行目はそのまま */
continue;
}
/* k行目に掛ける値を求める */
a = -sweep[i][k];
for (j = 0; j < N * 2; j++) {
/* i行目にk行目をa倍した行を足す */
/* これによりsweep[i][k]が0になる */
sweep[i][j] += sweep[k][j] * a;
}
}ここまでの処理により、sweep の k 列目の値が、対角成分の sweep[k][k] のみ 1 となり、それ以外の成分は 0 となることになります。
操作（２）と操作（３）を全ての対角成分に対して繰り返す
あとは、ここまで説明してきた処理を対角成分を移動させながら繰り返し行なっていけば良いだけです。
ここまで説明してきた処理は k 列目の対角成分に注目した処理ですので、注目する対角成分の列を表す k を 0 から N - 1 まで変化させながら、ここまで説明してきた処理を繰り返し実行してやれば良いです。
/* 全ての列の対角成分に対する繰り返し */
for (k = 0; k < N; k++) {
/* sweep[k][k]に掛けると1になる値を求める */
a = 1 / sweep[k][k];
/* 操作（２）：k行目をa倍する */
for (j = 0; j < N * 2; j++) {
/* これによりsweep[k][k]が1になる */
sweep[k][j] *= a;
}
/* 操作（３）によりk行目以外の行のk列目を0にする */
for (i = 0; i < N; i++) {
if (i == k) {
/* k行目はそのまま */
continue;
}
/* k行目に掛ける値を求める */
a = -sweep[i][k];
for (j = 0; j < N * 2; j++) {
/* i行目にk行目をa倍した行を足す */
/* これによりsweep[i][k]が0になる */
sweep[i][j] += sweep[k][j] * a;
}
}
}これにより、sweep の 0 列目から N - 1 列目までが単位行列となり、N 列目から 2 * N - 1 列目までが行列 mat の逆行列となります。
スポンサーリンク
逆行列部分のみ抽出する
最後に sweep から逆行列の部分のみを抽出して行列 inv に格納してやれば、逆行列 inv が完成したことになります。
具体的には、下記の処理により逆行列 inv が完成します。
/* sweepの右半分がmatの逆行列 */
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
inv[i][j] = sweep[i][N + j];
}
}掃き出し法のプログラム例（ゼロ除算対策なし）
ここまで紹介してきたソースコードを１つにまとめた例が下記となります。
#include <stdio.h>
#include <math.h>
#define N 3 /* 逆行列を求める行列の行数・列数 */
#define MAX_ERR 1e-10 /* 許容する誤差 */
int check(double mat[N][N], double inv[N][N]) {
double inner_product;
int i, j, k;
double ans;
double diff;
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
inner_product = 0;
for (k = 0; k < N; k++) {
inner_product += mat[i][k] * inv[k][j];
}
ans = (i == j) ? 1 : 0;
diff = ans - inner_product;
if (fabs(diff) > MAX_ERR) {
return 0;
}
}
}
return 1;
}
int main(void) {
/* 逆行列を求める行列用の２次元配列 */
double mat[N][N];
/* 逆行列用の２次元配列 */
double inv[N][N];
/* 掃き出し法を行う行列 */
double sweep[N][N * 2];
int i; /* 行 */
int j; /* 列 */
int k; /* 注目対角成分が存在する列 */
double a; /* 定数倍用 */
/* 正方行列の各成分をセット */
mat[0][0] = 2; mat[0][1] = 1; mat[0][2] = 4;
mat[1][0] = 4; mat[1][1] = 3; mat[1][2] = 4;
mat[2][0] = 1; mat[2][1] = 0; mat[2][2] = 2;
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
/* sweepの左側に逆行列を求める行列をセット */
sweep[i][j] = mat[i][j];
/* sweepの右側に単位行列をセット */
sweep[i][N + j] = (i == j) ? 1 : 0;
}
}
/* 全ての列の対角成分に対する繰り返し */
for (k = 0; k < N; k++) {
/* sweep[k][k]に掛けると1になる値を求める */
a = 1 / sweep[k][k];
/* 操作（２）：k行目をa倍する */
for (j = 0; j < N * 2; j++) {
/* これによりsweep[k][k]が1になる */
sweep[k][j] *= a;
}
/* 操作（３）によりk行目以外の行のk列目を0にする */
for (i = 0; i < N; i++) {
if (i == k) {
/* k行目はそのまま */
continue;
}
/* k行目に掛ける値を求める */
a = -sweep[i][k];
for (j = 0; j < N * 2; j++) {
/* i行目にk行目をa倍した行を足す */
/* これによりsweep[i][k]が0になる */
sweep[i][j] += sweep[k][j] * a;
}
}
}
/* sweepの右半分がmatの逆行列 */
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
inv[i][j] = sweep[i][N + j];
}
}
/* 逆行列invを表示 */
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
printf("%f, ", inv[i][j]);
}
printf("\n");
}
/* 検算 */
if (check(mat, inv)) {
printf("invはmatの逆行列です！！\n");
} else {
printf("invはmatの逆行列になってません...\n");
}
return 0;
}実行すれば下記のように inv の各成分、すなわち mat の逆行列が表示されます。
-1.500000, 0.500000, 2.000000, 1.000000, 0.000000, -2.000000, 0.750000, -0.250000, -0.500000, invはmatの逆行列です！！
mat の各成分の値は私がてきとうに決めていますので、逆行列を求めたい行列に合わせて適宜変更していただければと思います。N を変更すれば扱う正方行列の次数も変更できます。
また、check 関数は、mat と inv の積が単位行列になるかどうかを確認するための関数になります。ただし、掃き出し法の操作を行なっている際にゼロ除算が発生した場合はうまく動作しないので注意してください。
check 関数では、単純に２つの行列の積を計算しながら、各成分の値が確定した際に単位行列の成分と一致するかどうかを判断しています。全ての成分が単位行列と一致すれば、inv は mat の逆行列になっているため 1 を返却し、１つでも成分が異なるものがあれば 0 を返却します。
ただし、扱っているのが浮動小数点数で誤差が発生するため、一致するかどうかは比較対象の２つの値の差の絶対値が MAX_ERR 以下であるかどうかで判断しています。
C言語での行列の積の求め方については下記ページで解説していますので、積の求め方を知りたい方は下記を参照していただければと思います。
【C言語】行列の積の求め方
で、見出しにもつけていますが、上記のソースコードの場合、mat の成分によってはゼロ除算が発生し、逆行列がうまく求められない or 実行時にエラーが発生する、といった現象が起こります。
例えば mat[0][0] = 2; の部分を mat[0][0] = 0; に変更すれば、k が 0 の時の a = 1 / sweep[k][k]; でゼロ除算が発生し、その結果として表示される逆行列は私の環境では下記のようになりました（前述の通り、ゼロ除算が発生すると check 関数がうまく動作しないのでその点もご注意ください）。
nan, nan, nan, nan, nan, nan, nan, nan, nan, invはmatの逆行列です！！
環境によって結果が異なる可能性はありますが、いずれにせようまく逆行列は求められないはずです。
もちろん、全ての正方行列で逆行列を求めることが可能というわけではないですが（正方行列が特異行列の場合は正方行列は求められない）、上記のソースコードの mat においては、mat[0][0] = 0; に変更した場合でも理論上は逆行列を求めることが可能です。
なので、単に上記のソースコードの実装が不十分なだけで、ゼロ除算が行われないように対策すれば逆行列を求めることができるようになります。
次は、このゼロ除算対策を行う方法について解説していきます。
ゼロ除算対策
前述の解説の通り、注目している対角成分が 0 の場合は操作（２）を行う際の倍数の計算時にゼロ除算が発生します。
スポンサーリンク
操作（１）を行なってゼロ除算を防ぐ
このゼロ除算は、操作（１）、すなわち「２つの行を入れ替える」ことで回避することができます。
具体的には、対角成分を 1 にするために行う操作（２）の前に下記のような処理を行います。
まず注目している対角成分以降の行において、その対角成分が存在する列の中から最大の絶対値を持つ成分を探索します。
そして、その「最大の絶対値持つ成分が存在する行」と、「注目している対角成分が存在する行」を入れ替えます。すなわち、この２つの行に対して操作（１）を行います。
そうすれば、もし注目した対角成分が 0 の場合であっても、行の入れ替えにより対角成分の絶対値が 0 よりも大きなものに入れ替わることになります。そしてその後に操作（２）を行うようにすることで、ゼロ除算を防ぐことができます。
このように除算を行う対角成分を選ぶことを「ピボット選択」と呼びます（部分ピボット選択）
ただし、最大の絶対値が 0 である場合もあり得ます。
この場合は、逆行列を求めようとしている行列が「逆行列を持たない」ことを意味します。つまり、元々逆行列が求められないケースということになるので、逆行列を求める処理を終了してしまえば良いです
また、絶対値が最大の成分を探すのは、あくまでも注目している対角成分以降の行のみです。
注目している対角成分よりも上側の行に関しては、ここまでの処理によって既に対角成分が 1 になっています。その行を他の行に入れ替えてしまうと、せっかく 1 にした対角成分が他の値に変化し、逆行列がうまく求められなくなってしまいます。
最大の絶対値の成分を持つ行と入れ替える理由
ゼロ除算を防ぐためであれば「最大の絶対値の成分を持つ行」と入れ替える必要もなく、「0 以外の成分を持つ行」であればどの行と入れ替えても良さそうに感じるかもしれません。
「最大の絶対値の成分を持つ行」と入れ替えるのには、ゼロ除算を防ぐためだけでなく「浮動小数点数を扱う際に発生する誤差を出来るだけ抑えるため」という理由もあります。
絶対値の小さな値で割り算を行うと誤差が大きくなる可能性がありますので、行の入れ替えを行なって絶対値の大きな値で割り算が行えるようにすることで、誤差を極力小さくするようにしているというわけです（ただ場合によっては誤差が大きくなることもあるようです）。
ゼロ除算対策の実装
では、ここまで解説してきたゼロ除算対策を実装していきましょう！
このゼロ除算対策は、操作（２）を行うための「倍数を計算する前」に行う必要があります。
したがって、掃き出し法のプログラム例（ゼロ除算対策なし） のソースコードにおける a = 1 / sweep[k][k]; の行よりも前、すなわち、k に対する forループ内の最初の処理として、このゼロ除算対策の処理を追記します。
具体的には、下記のゼロ除算対策を追記します。
/* 最大の絶対値を注目対角成分の絶対値と仮定 */
double max = fabs(sweep[k][k]);
int max_i = k;
/* k列目が最大の絶対値となる行を探す */
for (i = k + 1; i < N; i++) {
if (fabs(sweep[i][k]) > max) {
max = fabs(sweep[i][k]);
max_i = i;
}
}
if (fabs(sweep[max_i][k]) <= MAX_ERR) {
/* 逆行列は求められない */
printf("逆行列は求められません...\n");
return 0;
}
/* 操作（１）：k行目とmax_i行目を入れ替える */
if (k != max_i) {
for (j = 0; j < N * 2; j++) {
double tmp = sweep[max_i][j];
sweep[max_i][j] = sweep[k][j];
sweep[k][j] = tmp;
}
}まず、下記部分で k 行目と入れ替える行の選択を行なっています（k 行目は現在注目している対角成分が存在する行）。
具体的には、k 行目以降の行から k 列目の成分の絶対値が最大となる行を探索し、その行を max_i としています。
/* 最大の絶対値を注目対角成分の絶対値と仮定 */
double max = fabs(sweep[k][k]);
int max_i = k;
/* k列目が最大の絶対値となる行を探す */
for (i = k + 1; i < N; i++) {
if (fabs(sweep[i][k]) > max) {
max = fabs(sweep[i][k]);
max_i = i;
}
}max_i 行目・k 列目の成分の絶対値、すなわち fabs(sweep[max_i][k] が 0 の場合があるので、その場合は下記で「逆行列を求められない」と判断してプログラムを終了するようにしています。
ただし、sweep[max_i][k] には浮動小数点数同士の演算による誤差が含まれている可能性があるため、単に 0 と比較してもうまく判断できないことがあります。
そのため、ある程度誤差を考慮して比較が行えるよう、sweep[max_i][k] の絶対値が MAX_ERR 以下であるかどうかで判断するようにしています。
if (fabs(sweep[max_i][k]) <= MAX_ERR) {
/* 逆行列は求められない */
printf("逆行列は求められません...\n");
return 0;
}
逆行列が求められる場合は、下記で実際に sweep の k 行目と max_i 行目の入れ替えを行っています。これが操作（１）となります。
/* 操作（１）：k行目とmax_i行目を入れ替える */
if (k != max_i) {
for (j = 0; j < N * 2; j++) {
double tmp = sweep[max_i][j];
sweep[max_i][j] = sweep[k][j];
sweep[k][j] = tmp;
}
}この入れ替えを行なった時点で、sweep[k][k] は 0 以外の値になっているため、1 / sweep[k][k] 実行時のゼロ除算を防ぐことができたことになります。
スポンサーリンク
掃き出し法で逆行列を求めるプログラム例
以上が掃き出し法での逆行列の求め方の解説となります。
最後に、ここまでの解説全てを含めた「掃き出し法で逆行列を求めるプログラム」のソースコードの例を紹介しておきます。
#include <stdio.h>
#include <math.h>
#define N 3 /* 逆行列を求める行列の行数・列数 */
#define MAX_ERR 1e-10 /* 許容する誤差 */
int check(double mat[N][N], double inv[N][N]) {
double inner_product;
int i, j, k;
double ans;
double diff;
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
inner_product = 0;
for (k = 0; k < N; k++) {
inner_product += mat[i][k] * inv[k][j];
}
ans = (i == j) ? 1 : 0;
diff = ans - inner_product;
if (fabs(diff) > MAX_ERR) {
return 0;
}
}
}
return 1;
}
int main(void) {
/* 逆行列を求める行列用の２次元配列 */
double mat[N][N];
/* 逆行列用の２次元配列 */
double inv[N][N];
/* 掃き出し法を行う行列 */
double sweep[N][N * 2];
int i; /* 行 */
int j; /* 列 */
int k; /* 注目対角成分が存在する列 */
double a; /* 定数倍用 */
/* 正方行列の各成分をセット */
mat[0][0] = 0; mat[0][1] = 1; mat[0][2] = 4;
mat[1][0] = 4; mat[1][1] = 3; mat[1][2] = 4;
mat[2][0] = 1; mat[2][1] = 0; mat[2][2] = 2;
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
/* sweepの左側に逆行列を求める行列をセット */
sweep[i][j] = mat[i][j];
/* sweepの右側に単位行列をセット */
sweep[i][N + j] = (i == j) ? 1 : 0;
}
}
/* 全ての列の対角成分に対する繰り返し */
for (k = 0; k < N; k++) {
/* 最大の絶対値を注目対角成分の絶対値と仮定 */
double max = fabs(sweep[k][k]);
int max_i = k;
/* k列目が最大の絶対値となる行を探す */
for (i = k + 1; i < N; i++) {
if (fabs(sweep[i][k]) > max) {
max = fabs(sweep[i][k]);
max_i = i;
}
}
if (fabs(sweep[max_i][k]) <= MAX_ERR) {
/* 逆行列は求められない */
printf("逆行列は求められません...\n");
return 0;
}
/* 操作（１）：k行目とmax_i行目を入れ替える */
if (k != max_i) {
for (j = 0; j < N * 2; j++) {
double tmp = sweep[max_i][j];
sweep[max_i][j] = sweep[k][j];
sweep[k][j] = tmp;
}
}
/* sweep[k][k]に掛けると1になる値を求める */
a = 1 / sweep[k][k];
/* 操作（２）：k行目をa倍する */
for (j = 0; j < N * 2; j++) {
/* これによりsweep[k][k]が1になる */
sweep[k][j] *= a;
}
/* 操作（３）によりk行目以外の行のk列目を0にする */
for (i = 0; i < N; i++) {
if (i == k) {
/* k行目はそのまま */
continue;
}
/* k行目に掛ける値を求める */
a = -sweep[i][k];
for (j = 0; j < N * 2; j++) {
/* i行目にk行目をa倍した行を足す */
/* これによりsweep[i][k]が0になる */
sweep[i][j] += sweep[k][j] * a;
}
}
}
/* sweepの右半分がmatの逆行列 */
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
inv[i][j] = sweep[i][N + j];
}
}
/* 逆行列invを表示 */
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
printf("%f, ", inv[i][j]);
}
printf("\n");
}
/* 検算 */
if (check(mat, inv)) {
printf("invはmatの逆行列です！！\n");
} else {
printf("invはmatの逆行列になってません...\n");
}
return 0;
}掃き出し法のプログラム例（ゼロ除算対策なし） では mat[0][0] が 0 の場合は逆行列が上手く求められませんでしたが、上記のソースコードのプログラムであれば逆行列をうまく求めることができ、下記のように結果が表示されます。
-0.375000, 0.125000, 0.500000, 0.250000, 0.250000, -1.000000, 0.187500, -0.062500, 0.250000, invはmatの逆行列です！！
もちろん N を変更することでもっと大きな行数・列数の正方行列の逆行列を求めることも可能です。
例えば N の define を下記のように変更し、
#define N 4 /* 逆行列を求める行列の行数・列数 */さらに、mat の各成分を下記のように変更した場合でも逆行列を求めることができ、
/* 正方行列の各成分をセット */
mat[0][0] = 0; mat[0][1] = 1; mat[0][2] = 4; mat[0][3] = 5;
mat[1][0] = 4; mat[1][1] = 3; mat[1][2] = 4; mat[1][3] = 9;
mat[2][0] = 1; mat[2][1] = 0; mat[2][2] = 2; mat[2][3] = 7;
mat[3][0] = 8; mat[3][1] = 4; mat[3][2] = 1; mat[3][3] = 5;結果としては下記のように表示されることになります。
0.597561, -0.914634, 0.365854, 0.536585, -0.987805, 1.573171, -0.829268, -0.682927, 0.939024, -0.865854, 0.146341, 0.414634, -0.353659, 0.378049, 0.048780, -0.195122, invはmatの逆行列です！！
ただし、mat の成分によっては、mat が逆行列を持たない場合（特異行列の場合）もあるので注意してください。例えば matの各成分を下記のように設定した場合（N は 3）、mat は逆行列を持たないので逆行列を求めることができません。
/* 正方行列の各成分をセット */
mat[0][0] = 2; mat[0][1] = 1; mat[0][2] = 4;
mat[1][0] = 4; mat[1][1] = 2; mat[1][2] = 8;
mat[2][0] = 1; mat[2][1] = 0; mat[2][2] = 2;このような場合は、下記のようにメッセージが表示されるようなっています。
逆行列は求められません...
あと、扱う値が浮動小数点数のため誤差が発生する可能性もありますので、この点についてもご注意ください。特に次数が大きくなると誤差も大きくなります。
まとめ
このページでは、C言語での「掃き出し法を用いた逆行列の求め方」について解説しました！
逆行列を求める方法としては他の方法も存在しますが、掃き出し法は処理が直感的に理解しやすい方法であり、割と実装もしやすいのではないかと思います。
対角成分に注目しながら、列ごとに単位行列に近づけていくところがポイントだと思います。ただし、ゼロ除算が発生するとうまく逆行列が求められないため、行の入れ替えを行うことでゼロ除算の対策を行う必要があるので注意してください。
逆行列を求められるようになればプログラミングで出来ることの幅が広がります。もちろん実現可能な行列演算も増えますし、画像処理などでも逆行列が活躍します。
今後役に立つ可能性が高いですので、是非この逆行列の求め方は覚えておいてください！
