逆行列を求める際のオーソドックスな手法に掃き出し法というものがあります。やり方の解説やサンプルコードはすぐ見つかると思います。今回はコードと合わせてなぜそれで答えが出るのか試行錯誤してみます。
掃き出し法
掃き出し法とか、ガウスの消去法とか、ガウス・ジョルダンの消去法とか色々な呼ばれ方をしているこの手法ですが、基本的には連立1次方程式を解くための手法です。
具体的には、行列を構成する各行に対して、
行を入れ替える
行に0でない定数をかける
行に他の行の定数倍を加える
この3つの操作を繰り返していく事で行列を所望の形に変形させていく方法です。
この変形で単位行列を作っていく事で逆行列を求めるのが掃き出し法を使った逆行列の算出になります。
具体的なサンプルで確かめてみます。
$$ A=\left(\begin{array}{c}x_1 & x_2 & x_3 \\ x_4 & x_5 & x_6 \\ x_7 & x_8 & x_9\end{array}\right) $$
な行列を例にするとその行列の横に単位行列を置いた
$$ \left(\begin{array}{c}x_1 & x_2 & x_3 & | & 1 & 0 & 0\\ x_4 & x_5 & x_6 & | & 0 & 1 & 0\\ x_7 & x_8 & x_9 & | & 0 & 0 & 1 \end{array}\right) $$
な行列を作ります。この3行6列の行列に対して上の3つの操作を繰り返していき、左半分を単位行列に変形していきます。
まずは1行目左上を1にします。これは1行目を\(1/x_1\)することで求まります。
$$ \left(\begin{array}{c}1 & x_2/x_1 & x_3/x_1 & | & 1/x_1 & 0 & 0\\ x_4 & x_5 & x_6 & | & 0 & 1 & 0\\ x_7 & x_8 & x_9 & | & 0 & 0 & 1 \end{array}\right) $$
続いて2行目と3行目の左を0にします。これは1行目をそれぞれ\(x_4\)倍、\(x_7\)倍して各行から減じれば求まります。
$$ \left(\begin{array}{c}1 & x_2/x_1 & x_3/x_1 & | & 1/x_1 & 0 & 0\\ 0 & x_5-x_4*x_2/x_1 & x_6-x_4*x_3/x_1 & | & -x_4/x_1 & 1 & 0\\ 0 & x_8-x_7*x_2/x_1 & x_9-x_7*x_3/x_1 & | & -x_7/x_1 & 0 & 1 \end{array}\right) $$
続いて2行目の2列目を1にします。これも同じように2行目を\(\frac{1}{x_5-x_4x_2/x_1}\)することで求まります。
$$ \left(\begin{array}{c}1 & x_2/x_1 & x_3/x_1 & | & 1/x_1 & 0 & 0\\ 0 & 1 & \frac{x_6-x_4*x_3/x_1}{x_5-x_4*x_2/x_1} & | & \frac{-x_4/x_1}{x_5-x_4*x_2/x_1} & \frac{1}{x_5-x_4*x_2/x_1} & 0\\ 0 & x_8-x_7*x_2/x_1 & x_9-x_7*x_3/x_1 & | & -x_7/x_1 & 0 & 1 \end{array}\right) $$
とこんな感じで左から順に単位行列な形を作っていくと、最終的には右半分が逆行列になる。そんな具合です。
$$ \left(\begin{array}{c}1 & 0 & 0 & | & x’_1 & x’_2 & x’_3\\ 0 & 1 & 0 & | & x’_4 & x’_5 & x’_6\\ 0 & 0 & 1 & | & x’_7 & x’_8 & x’_9 \end{array}\right) $$
いわゆるこの手法的な解説は色々なサイトでされていると思います。確かにそのように変換すれば「なぜか」逆行列が求まります。
今回はこの「なぜか」をもう少し突っ込んで考えてみます。
3つの操作の意味
上述した3つの操作ですが、この3つの操作はすべて行列の掛け算で表現できる点がこの変換のポイントになります。
ある3×3の行列
$$ A = \left(\begin{array}{c}a & b & c \\ d & e & f \\ g & h & i \end{array}\right) $$
をサンプルに変換を行ってみます。
行を入れ替える
例えば1行目と2行目の入れ替えを考えます。その変換を実現させるためには左から
$$ X_1 = \left(\begin{array}{c}0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right) $$
を乗じる事で入れ替わります。試してみればすぐわかると思います。
$$ X_1A = \left(\begin{array}{c}0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right)\left(\begin{array}{c}a & b & c \\ d & e & f \\ g & h & i \end{array}\right)=\left(\begin{array}{c}d & e & f \\ a & b & c \\ g & h & i \end{array}\right) $$
一見回りくどい計算になりますが、この入れ替えの操作が行列の掛け算になっています。
行に0でない定数をかける
これも同様です。例えば3行目をj倍してみます。その変換は同様に
$$ X_2 = \left(\begin{array}{c}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & j \end{array}\right) $$
を左から乗じます。
$$ X_2A = \left(\begin{array}{c}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & j \end{array}\right)\left(\begin{array}{c}a & b & c \\ d & e & f \\ g & h & i \end{array}\right)=\left(\begin{array}{c}a & b & c \\ d & e & f \\ jg & jh & ji \end{array}\right) $$
単位行列をちょいといじるだけで行の定数倍が簡単に作れます。
行に他の行の定数倍を加える
最後にある行の定数倍を他の行に加える（減じる）操作です。これは結局上記2つの操作の組み合わせで実現できる操作で、例えば2行目をj倍して、3行目に足す。といった操作であれば、
$$ X_3 = \left(\begin{array}{c}1 & 0 & 0 \\ 0 & 1 & 0 \\0 & j & 1 \end{array}\right) $$
な行列を左から乗じる事で実現できてしまいます。
$$ X_3A = \left(\begin{array}{c}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & j & 1 \end{array}\right)\left(\begin{array}{c}a & b & c \\ d & e & f \\ g & h & i \end{array}\right)=\left(\begin{array}{c}a & b & c \\ d & e & f \\ jd+g & je+h & jf+i \end{array}\right) $$
手計算で簡単に確かめられます。
このようにこの3つの基本操作はすべて行列の掛け算で表すことができるということがわかりました。
掃き出し法による逆行列算出
先ほどのサンプルAに対してその逆行列を求めよ。といった場合になぜこの基本操作を繰り返していくと求まるのか？少し考えてみます。
掃き出し法による逆行列の算出は、基本操作となる行列Xの掛け算を繰り返していき、最終的にその形を単位行列にしていく手法です。この操作が意味することは、
$$ X_1X_2X_3X_4..X_nA = I $$
と変換を続けていく事になります。この連続操作\(X_1X_2X_3X_4..X_n\)を1つの行列とみなすと
$$ XA = I $$
となり、XがAの逆行列のなっていることがわかります。3つの操作の連続が逆行列を求める操作そのものになっているというわけです。
これをシステマチックに行っているのがこの掃き出し法になります。\(X_1X_2X_3X_4..X_n\)を求めるために、同じ操作を単位行列に対しても行って行きます。
$$ X_1X_2X_3X_4..X_nI = XI = X $$
となりXすなわちAの逆行列が導き出されます。
ひとまず手段と数学的な意味はこんなところで、続いてコードにしてみようと思います。
コード
Pythonを使います。検算が簡単だから。ニーズはないでしょうね。Pythonなら1行なので。C言語に移植しやすいようにCっぽく書いてみます。3行3列で書いていますが、基本的にはn行n列でも動くはずです。
上述したように、とある行列Aと単位行列Iを用意します。
n = 3
A = np.array([[random.random(), random.random(), random.random()],
[random.random(), random.random(), random.random()],
[random.random(), random.random(), random.random()]])
I = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]).astype(np.float64)
適当な乱数で作ります。
まずは対角成分を1にするための操作を行います。yが0の時は0,0の場所。1になると1,1といったようにgainには対角成分の逆数が入り、横方向にその対象となるy行の値をgain倍します。
for y in range(n): # 操作中の行がy
gain = 1/A[y, y] # 対角成分
for x in range(n):
A[y, x] = A[y, x] * gain # 対角成分を1へ
I[y, x] = I[y, x] * gain
yが0の時を例にすると、
$$ \left(\begin{array}{c}1 & x_2/x_1 & x_3/x_1 & | & 1/x_1 & 0 & 0\\ x_4 & x_5 & x_6 & | & 0 & 1 & 0\\ x_7 & x_8 & x_9 & | & 0 & 0 & 1 \end{array}\right) $$
この計算です。同じ計算を単位行列Iに対しても行います。
続いてy行目とは異なる行に対して、y行目の定数倍を引き算して0にします。
for yy in range(n): # 操作される行（消される行）
if y != yy: # 自分自身の行でないときに
gain = A[yy, y] # 対角成分
for x in range(n):
A[yy, x] = A[yy, x] - A[y, x] * gain
I[yy, x] = I[yy, x] - I[y, x] * gain
これもyを0を例にすると、0行目にそれぞれ\(gain=x_4 or x_7\)を掛けて、1行目と2行目から引きます。
$$ \left(\begin{array}{c}1 & x_2/x_1 & x_3/x_1 & | & 1/x_1 & 0 & 0\\ 0 & x_5-x_2/x_1*x_4 & x_6-x_3/x_1*x_4 & | & -1/x_1*x_4 & 1 & 0\\ 0 & x_8-x_2/x_1*x_7 & x_9-x_3/x_1*x_7 & | & -1/x_1*x_7 & 0 & 1 \end{array}\right) $$
この計算をしています。これをすべてのyに対して行うことでAが単位行列になり、Iに答えとなる逆行列が入ってきます。
ピボット選択
残念ながらこの計算には一部問題があります。ゼロ割が入ってしまう可能性があります。この
gain = 1/A[y, y] # 対角成分
を計算するときにA[y, y]がゼロである可能性があります。ゼロじゃ無いにしてもとても小さな値の時にgainが無限大に発散してしまい解の精度が悪くなります。
その場合にyより下にある行とy行を入れ替え、ゼロ割りを回避するのがピボットの選択とか入れ替えとか交換とかいう手法です。せっかく作ったyまでの行は残しておきたいですし。
入れ替えは通常その列の絶対値が最も大きいものと入れ替えるのがセオリーの様です。
例えば2行目の計算中に3行目の値の値の方が大きい\(x_4 < x_7\)となったら、
$$ \left(\begin{array}{c}1 & x_1 & x_2 & | & x_3 & 0 & 0\\ 0 & x_4 & x_5 & | & x_6 & 1 & 0\\ 0 & x_7 & x_8 & | & x_9 & 0 & 1 \end{array}\right) \Rightarrow \left(\begin{array}{c}1 & x_1 & x_2 & | & x_3 & 0 & 0\\ 0 & x_7 & x_8 & | & x_9 & 0 & 1\\ 0 & x_4 & x_5 & | & x_6 & 1 & 0 \end{array}\right) $$
2行目と3行目を入れ替えてやる。（その時1行目はそっとしておく）という計算になります。入れ替えるも基本操作の1つです。
泥臭く、Cっぽいコードにすると行の入れ替えは
max = abs(A[y, y])
indx = y
for yy in range(y + 1, n):
if max < abs(A[yy, y]):
max = abs(A[yy, y])
indx = yy
if indx != y:
for x in range(n):
tmp = A[indx, x]
A[indx, x] = A[y, x]
A[y, x] = tmp
tmp = I[indx, x]
I[indx, x] = I[y, x]
I[y, x] = tmp
こんな感じ。本来は入れ替えた上で最大値がゼロだったら逆行列が存在しない。ということになってエラーを返した方がいいかもしれません。
検算付の全コードです。
import numpy as np
import random
n = 3
A = np.array([[random.random(), random.random(), random.random()],
[random.random(), random.random(), random.random()],
[random.random(), random.random(), random.random()]])
I = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]).astype(np.float64)
invA = np.linalg.inv(A)
for y in range(n): # 操作中の行がy
max = abs(A[y, y])
indx = y
for yy in range(y + 1, n):
if max < abs(A[yy, y]):
max = abs(A[yy, y])
indx = yy
if indx != y:
for x in range(n):
tmp = A[indx, x]
A[indx, x] = A[y, x]
A[y, x] = tmp
tmp = I[indx, x]
I[indx, x] = I[y, x]
I[y, x] = tmp
gain = 1/A[y, y] # 対角成分
for x in range(n):
A[y, x] = A[y, x] * gain # 対角成分を1へ
I[y, x] = I[y, x] * gain
for yy in range(n): # 操作される行（消される行）
if y != yy: # 自分自身の行でないときに
gain = A[yy, y] # 対角成分
for x in range(n):
A[yy, x] = A[yy, x] - A[y, x] * gain
I[yy, x] = I[yy, x] - I[y, x] * gain
めでたくinvAとIが一致しました。
まとめ
掃き出し法による逆行列の計算を行ってみました。手法としては知ってはいたのですが、なぜ？の部分がイマイチ残っていたので今回自分用のメモとしてまとめてみました。
解が一致したので計算は間違っていないと思います。理解の足しになれば。
今回C言語への移植を考慮してべたべたなコードを書いています。Pythonで書いてあっても誰もうれしく無いだろうし。
おまけ（C言語で書くと…）
ちなみにCで書くと
for(int y=0;y<n;y++){ //操作中の行がy
double max = abs(A[y][y]);
int indx = y;
for(int yy=y+1;yy<n;yy++){
if(max < abs(A[yy][y])){
max = abs(A[yy][y]);
indx = yy;
}
}
if(indx != y){
for(int x=0;x<n;x++){
double tmp = A[indx][x];
A[indx][x] = A[y][x];
A[y][x] = tmp;
tmp = I[indx][x];
I[indx][x] = I[y][x];
I[y][x] = tmp;
}
}
double gain = 1.0/A[y][y]; //対角成分
for(int x=0;x<n;x++){
A[y][x] *= gain; //対角成分を1へ
I[y][x] *= gain;
}
for(int yy=0;yy<n;yy++){ //操作される行（消される行）
if(y != yy){ //自分自身の行で無いときに
gain = A[yy][y];
for(int x=0;x<n;x++){
A[yy][x] -= A[y][x] * gain;
I[yy][x] -= I[y][x] * gain;
}
}
}
}
こんな感じですかね。こっちの方がニーズがありますかね。ほとんど同じコードです。コンパイラが無いので試してませんが、多分こんな感じ。
コメント