逆行列を計算するEXCELマクロ
逆行列は非常に適用範囲の広いものであるだけでなく、その計算原理は他の計算問題に応用可能であ
り、一度はマスターしておくことが望ましいものである。
まず、行列に不慣れな学生の便宜のために、逆行列の定義を説明しておこう。行列とは、縦方向と横
方向とにデータを長方形に並べた表のようなものであり、その位置も意味を持つものである。したがっ
て、EXCELで使用される２次元配列と全く同じものとなる。行列における縦方向の行数と横方向の列数
は全く任意であり、行数がｍ、列数がｎの行列は簡単に「ｍ×ｎ次の行列」と呼ばる。また、行数と列
数が等しい行列は特に「正方行列」と呼ばれ、ｍ×ｍ次の行列は「ｍ次の正方行列」と呼ばれる。さら
に、上からｉ行目で左からｊ番目の要素は、その行列の「i-j 要素」と呼ばれる。
行列同士の演算も定義されているので、和、差、積について簡単に説明しておく。２つの行列同士の
和ないし差とは、両行列の対応する要素同士の和ないし差を求めることである。したがって、和ないし
差が計算できるためには、両行列の行数および列数が一致していなければならない。また、行列ＡとＢ
の和は A+B と記述され、行列ＣとＤの差は C-D と記述される。
行列同士の積は若干複雑であるので、例によって説明することにする。行列ＸとＹの積計算の結果を
行列Ｚとすることにしよう。すなわち、計算により行列Ｚの各要素を求めることになる。行列Ｚのi-j
要素は、行列Ｘの「i番目の行」の全要素と、行列Ｙの「j番目の列」の全要素を取り出し、それぞれの
対応する位置にある要素同士を掛け、それらを総和したものとなる。したがって、行列Ｘの列数と行列
Ｙの行数とが一致していなければならない。すなわち、Ｘがｍ×ｎ次の行列であり、Ｙがｐ×ｑ次の行
列であるとき、ＸとＹの積が計算可能であるためには、n=p でなければならない。また、このような２
行列の積計算の結果であるＺはｍ×ｑ次の行列となる。行列Ｚの一つの要素を求めるためにｎ（ないし
ｐ）回のかけ算とその結果の総和計算が必要であり、行列Ｚの全要素を求めるためには、かけ算だけで
もｐ×ｍ×ｑ回行わなければならない。小さな行列同士の積であれば手計算でも可能であろうが、少し
でも大きくなればパソコンを含めた計算機の力を借りた方がよいことになる。そこで、行列同士の積計
算を行うプログラムを考えてみよう。
上で説明に使用した例を用いることとし、行列Ｘの各要素の値が配列Xに代入されており、行列Yの各
要素の値が配列Ｙに代入されているものとし、計算結果を配列Ｚに代入させるものとする。先の積計算
の説明から明らかなように、Z(i, j) すなわち行列Ｚのi-j要素は
X(i, 1)*Y(1, j)+X(i, 2)*Y(2, j)+ ･･･+X(i, n)*Y(n, j)
と計算されるものである。この計算をするためのプログラムが
s = 0
For k = 1 To n
s = s + X(i, k) * Y(k, j)
Next k
Z(j, j) = s
となることは明らかであろう。さらに、このプログラムが Z(i, j) を求めるためのものであり、配列
ないし行列Ｚを求めるためには、ｉを1,2,･･･,m、ｊを1,2,･･･,q と変化させなければならないので、
行列（配列）の積を求めるプログラムは次のようになる。
For i = 1 To m
For j = 1 To q
s = 0
For k = 1 To n
s = s + X(i, k) * Y(k, j)
Next k
Z(i, j) = s
Next j
Next i
ここで、本節の主題である逆行列に話を戻そう。逆行列とは普通の数値における逆数に対応するもの
であり、２つの正方行列ＡとＢの積が単位行列（ＩないしＥと記述される）になるとき、「ＡはＢの逆
行列」ないし「ＢはＡの逆行列」と呼ばれる。また、行列Ａの逆行列は A-1 と記述される。なお、単
位行列とは、その対角要素がすべて１で、他の要素がすべて０であるような正方行列のことである。
このような逆行列を求める方法として、色々なものが提案されているが、ここでは最も基本的な「掃
き出し法」と呼ばれる方法について説明する。掃き出し法の原理は、行列の「基本変形」を繰り返すこ
とで、与えられた（逆行列を求めるべき）行列を単位行列にすることである。行列の基本変形とは、対
象の行列にある種の行列をかけることに相当し、基本変形を繰り返すことは順次行列をかけていくこと
に相当する。そして、このような繰り返しの結果として単位行列が得られたとすれば、そのような基本
変形の合成、すなわち、基本変形に対応した行列をすべてかけ合わせたものが逆行列となる。例えば
Hr Hr-1 ･･･H2 H1 A = I
となるものとしよう。ただし、A は与えられた行列であり、H1、H2、･･･、Hr は繰り返し行われた基本
変形を行列表現したものである。すなわち、ｒ回の基本変形で行列Ａを単位行列Ｉに変換できたものと
する。このとき、ｒ個の行列の積 HrHr-1･･･H2H1 をＢで表せば、BA=I となり、行列Ｂが行列Ａの逆行
列となる。
つぎに、行列の基本変形について説明しよう。基本変形には三つのタイプのものがある。基本変形は
行と列とのいずれに対しても可能であるが、ここでは行に対する基本変形に説明を限定する。すなわち
① ある行の全要素にある数値をかける。
┌
│
└
|
1 2 3
4 5 6
7 8 9
|
┐
│
┘
|
|
→
|
┌
│
└
|
3 6 9
4 5 6
7 8 9
|
┐
│
┘
|
|
｛第１行に３をかける｝
|
② 任意の２つの行を入れ換える。
┌
│
└
|
1 2 3
4 5 6
7 8 9
|
┐
│
┘
|
|
→
|
┌
│
└
|
1 2 3
7 8 9
4 5 6
|
┐
│
┘
|
|
｛第２行と第３行を入れ換える｝
|
③ ある行の全要素にある数値をかけたものを別の行に加える。
┌
│
└
|
1 2 3
4 5 6
7 8 9
|
┐
│
┘
|
|
→
|
┌
│
└
|
1 2 3
2 1 0
7 8 9
|
┐
│
┘
|
|
｛第１行に -2 をかけたものを第２行に加える｝
|
である。また、これらの基本変形に対応した行列は
①
|
┌
│
└
|
3 0 0
0 1 0
0 0 1
|
┐
│
┘
|
|
②
|
┌
│
└
|
1 0 0
0 0 1
0 1 0
|
┐
│
┘
|
|
③
|
┌
│
└
|
1 0 0
0 1 0
-2 0 1
|
┐
│
┘
|
|
となる（対象に行列に、これらの行列を左からかけた結果が、先の基本変形と結果と同じになることを
確認せよ）。
掃き出し法とは、このような基本変形を繰り返して、与えられた行列を単位行列に変換するものであ
るが、いっぺんに単位行列に変換するのではなく、１列ずつ処理していく方法である。すなわち、まず
最初に、1-1 要素を１にし、第１列の他の要素を０にする。次に、2-2 要素を１にし、第２列の他の要
素を０にする。などの処理を繰り返して、最後に、n-n 要素を１にし、第ｎ列の他の要素を０にして、
作業を完了するものである（ただし、対象の行列がｎ次の正方行列とする）。
具体的な方法を説明するために、ここでは３次の正方行列Ａ（２次元配列Ａに代入されているものと
する）
|
┌
│
└
|
2 3 4
5 6 7
8 9 0
|
┐
│
┘
|
|
の逆行列を求めること考えてみよう。まず、1-1 要素を１にしなければならないが、このために基本変
形①を使用する。すなわち、元の 1-1 要素の値 2 の逆数 0.5 を第１行にかける。
|
┌
│
└
|
1 1.5 2
5 6 7
8 9 0
|
┐
│
┘
|
|
次に、第１列の他の要素を０にするために、基本変形③を２回行う。すなわち、第１行に -5 をかけた
ものを第２行に加え、さらに、第１行に -8 をかけたものを第３行に加える。
|
┌
│
└
|
1 1.5 2
0 -1.5 -3
0 -3 -16
|
┐
│
┘
|
|
以下、同様な処理を第２列と第３列について行えば、次のようになり、対象の行列を単位行列に変換す
ることができる。
┌
│
└
|
1 1.5 2
0 1 2
0 -3 -16
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 0 -1
0 1 2
0 0 -10
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 0 -1
0 1 2
0 0 1
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 0 0
0 1 0
0 0 1
|
┐
│
┘
|
|
それでは、逆行列はどこに求まるのであろうか。逆行列を求めるためには、ここで用いられた基本変形
を行列表現に直し、それらをかけ合わせなければならないが、これは大変な作業となる。前述の表現を
借りれば、逆行列を求める作業は
HrHr-1･･･H2H1A = I
となるような基本変形の並び H1、H2、･･･、Hr を求めることであったが、[A I] なるn×2n次の行列に
これらの基本変形を行ったとすれば
HrHr-1･･･H2H1[A I] = [I HrHr-1･･･H2H1]
となり、右辺の行列の右半分の部分に逆行列が求まることになる（ただし、行列Hiをかけることは、対
応する基本変形を行うことと等価であることに注意せよ）。したがって、ここでの例においては、次の
ようにすればよいことになる。
┌
│
└
|
2 3 4 1 0 0
5 6 7 0 1 0
8 9 0 0 0 1
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 1.5 2 0.5 0 0
5 6 7 0 1 0
8 9 0 0 0 1
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 1.5 2 0.5 0 0
0 -1.5 -3 -2.5 1 0
0 -3 -16 -4 0 1
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 1.5 2 0.5 0 0
0 1 2 5/3 -2/3 0
0 -3 -16 -4 0 1
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 0 -1 -2 1 0
0 1 2 5/3 -2/3 0
0 0 -10 1 -2 1
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 0 -1 -2 1 0
0 1 2 5/3 -2/3 0
0 0 1 -0.1 0.2 -0.1
|
┐
│
┘
|
|
→ |
┌
│
└
|
1 0 0 -2.1 1.2 -0.1
0 1 0 5.6/3 -3.2/3 0.2
0 0 1 -0.1 0.2 -0.1
|
┐
│
┘
|
|
以上の説明で掃き出し法について理解できたことと思うが、ここでの例のような小さな行列の場合で
も、逆行列を求めるためにはかなりの手間が必要である。しかるに、このような単純作業の繰り返しこ
そが計算機の得意とするものであるので、この種の計算は計算機に任せるべきであろう。
それでは、掃き出し法により逆行列を計算するためのプログラムを考えてみよう。いま、大きさ n×
2n の配列Ａの左半分に逆行列を求めるべき行列が代入されており、右半分には単位行列が代入されて
いるものとしよう。このとき最初にすることは、基本変形によって第１列を[1, 0,･･･, 0] とすること
である。そのためには、まず、A(1, 1) が１となるようにしなければならないが、そのためには第１行
の全要素を A(1, 1) で割ればよい（A(1, 1)の逆数をかけることと同じ）。プログラムは
x = A(1, 1)
For j = 1 To n + n
A(1, j) = A(1, j) / x
Next j
となる。ここで注意すべきことは、A(1, 1) の値を変数名Xに取っておいてから割り算を行わねばなら
ないことである。よくやる誤りは
For j = 1 To n + n
A(1, j) = A(1, j) / A(1, 1)
Next j
としてしまうことである。このプログラムでは、最初にｊが１のときに A(1, 1)=A(1, 1)/A(1, 1) と
なり、A(1, 1) に１が代入されてしまう。さらに、ｊが２以上となったとき、この A(1, 1)、すなわ
ち、１で割られることになり、当初の目的が達成されない。また、終端値がn+nになっていることに注
意せよ。
次に、第１列の他の要素を０とするような基本変形を行わねばならない。例えば、A(2, 1) を０と
するためには、A(1, 1) が既に１となっているので、第１行に -A(2, 1) をかけたものを第２行に加
える、ないし、第１行に A(2, 1) をかけたものを第２行から引けばよい。すなわち
x = A(2, 1)
For j = 1 To n + n
A(2, j) = A(2, j) - x * A(1, j)
Next j
このときも、A(2, 1) の値を取っておいてから計算しなければならないことは明らかであろう。この
プログラムは、A(2, 1) の部分を０にするためのものであるが、A(3, 1)、･･･、A(n, 1) を０にする
ためのプログラムも同型であり、2の部分を 3、･･･、n と変えるだけでよい。したがって、A(2, 1)、
A(3, 1)、･･･、A(n, 1) を0にするためのプログラムは次のようになる。
For i = 2 To n
x = A(i, 1)
For j = 1 To n + n
A(i, j) = A(i, j) - x * A(1, j)
Next j
Next i
以上のプログラムで、第1 列に対する処理は完了する。
第2 列に対する処理も同様にして
x = A(2, 2)
For j = 2 To n + n
A(2, j) = A(2, j) / x
Next j
For i = 1 To n
If i <> 2 THEN
x = A(i, 2)
For j = 2 To n + n
A(i, j) = A(i, j) - x * A(2, j)
Next j
EndIf
Next i
となる。このプログラムは第１列に対するものとほぼ同型であるが、いくつかの点で異なっている。
その第１のものは、0にすべき場所が A(1, 2)、A(3, 2)、･･･、A(n, 2)と飛んでいるので、そのため
の処理がIf文として加えられていることである。したがって、第１列に対するプログラムをこちらに
合わせれば
For i = 1 To n
If i <> 1 THEN
x = A(i, 1)
For j = 1 To n + n
A(i, j) = A(i, j) - x * A(1, j)
Next j
EndIf
Next i
となることは明らかであろう。第二の差異は、変数名ｊを用いた２カ所の For 文の初期値が２になっ
ていることである。これは、掃き出し法の原理に戻って考えれば明らかであるが、第２列を処理する段
階では、A(2, 1) は既に０となっており、最初の For 文において、それをｘで割ったとしても０であ
ることには変わりがないし、第二の For 文において、それにｘをかけて第１列の他の要素から引いた
としても変化はない。すなわち、無駄な計算を省くために、このように初期値が２とされている。
さらに、第３列以降を処理するプログラムは、この第２列に対するプログラムで、２となっている場
所をすべて変えることでよいことも明らかであろう。したがって、各列を処理するためのプログラムは
同型となり、全体を For～Next 文で処理できることになる。すなわち
For k = 1 To n
x = A(k, k)
For j = k To n + n
A(k, j) = A(k, j) / x
Next j
For i = 1 To n
If i = k THEN
x = A(i, k)
For j = k To n + n
A(i, j) = A(i, j) - x * A(k, j)
Next j
EndIf
Next i
Next k
賢明な学生は気が付いているかもしれないが、このプログラムの３行目において、ｘすなわちA(k,k)
で割り算を行っているが、A(k, k) の値が０となることはないのか、という疑問が生じるであろう。実
際、そのようなことが起こる可能性は多いにある。そこで、このような問題を解決するための処理を付
け加えなければならない。この処理は「ピボット処理」と呼ばれるものであり、これまで使用してこな
かった基本変形②を用いるものである。すなわち、第k列を処理するときに A(k, k) が０であったとし
ても、その下の部分に０でない値が存在すれば、その行とｋ行とを入れ換えることで、A(k, k) の部分
を０以外の値にすることができる。ここで、下の部分としたのは、上の部分に０以外の値があったとし
てその行とｋ行を入れ換えると、せっかく処理の済んでいる第ｋ列以前の列の形を崩すことになり、困
ったことになるからである。次に示すプログラムでは、「ピボット処理」として、A(k, k) 以下の部分
（A(k, k)、A(k+1, k)、･･･、A(n, k)の中）で絶対値が最大なものをみつけ、その最大値のあった行と
ｋ行とを入れ換えるようになっている。また、A(k, k) 以下の部分がすべて０のときは、与えられた行
列が、逆行列の存在しないものであることを意味するので、そのことをメッセージで表示して、実行を
終了するようにしてある。
n = ActiveCell
ActiveCell.Offset(1, 0).Range("A1").Select
nn = n + n
ReDim A(n, nn)
|
|
For i = 1 To n
For j = 1 To n
A(i, j) = ActiveCell
ActiveCell.Offset(0, 1).Range("A1").Select
Next j
ActiveCell.Offset(1, -n).Range("A1").Select
Next i
|
① |
ActiveCell.Offset(1, 0).Range("A1").Select
|
|
For i = 1 To n
For j = n + 1 To nn: A(i, j) = 0: Next j
A(i, i + n) = 1
Next i
|
② |
'
For k = 1 To n
|
|
p = 0.0000000001: c = 0
For i = k To n
If Abs(A(i, k)) > p Then p = Abs(A(i, k)): c = i
Next i
|
③ |
If c = 0 Then ActiveCell = "逆行列が存在しない": GoTo xyz
|
④ |
If c <> k Then
|
⑤ |
For j = k To nn
y = A(k, j): A(k, j) = A(c, j): A(c, j) = y
Next j
|
⑥ |
End If
|
|
x = A(k, k)
For j = k To nn
A(k, j) = A(k, j) / x
Next j
|
⑦ |
For i = 1 To n
If i <> k Then
x = A(i, k)
For j = k To nn
A(i, j) = A(i, j) - x * A(k, j)
Next j
End If
Next i
|
⑧ |
Next k
'
|
|
For i = 1 To n
For j = n + 1 To nn
ActiveCell = A(i, j)
ActiveCell.Offset(0, 1).Range("A1").Select
Next j
ActiveCell.Offset(1, -n).Range("A1").Select
Next i
|
⑨ |
xyz:
|
⑩ |
このプログラムの①は、逆行列を求めるべき行列を読み込み、配列Aの左半分に代入するためのもの
であり、②は右半分に単位行列を代入するためのものである。
③は、A(k, k) 以下の部分で絶対値が最大なものの場所（行）をみつけるためのものである。
④は、A(k, k) 以下の部分に絶対値が０より大きいものが存在しない、すなわち、その部分の値がす
べて０であったとき、計算の継続が不可能、すなわち、与えられた行列には逆行列が存在しないので、
その旨を表示して実行を終了すもためのものである。
⑤の命令は、A(k, k) 以下の部分での最大値が A(k, k) 自身であったときに、自分自身を自分自身
と入れ換えるという無駄な計算をさせないためのものである。しかし、この命令がなかったとしても、
多少余分な計算をするだけで問題はない。
⑥は、ｋ行とｃ行の要素を入れ替えるためのものである。このとき、ｋ列より左側は共に０となって
いるので、入れ替える必要はなく、For j = k となっている。
⑦は、ｋ行の各要素を A(k, k) で割り、A(k, k) 要素を１にするためのものである。ここでもｋ列
より左側は０となっているので、割る必要がなく、For j = k となっている。
⑧は、⑦で１となった A(k, k) を用いてｋ列の各要素を０にするためのものである。ただし、i =
kのときは、自分自身で自分を引いてしまい、誤った計算をするため、i = k のときは計算しないよう
に、すなわち i<>k のときのみ計算するようになっている。
⑨は、求められた逆行列を表示するためのものである。このプログラムは、一般に行列を表示すると
きに使用可能なものであるので憶えておくと便利であろう。
⑩は、④でのジャンプ先を示すラベルである。
|