行列のapplicationとして、逆行列の計算を行うmrevを考えます。アルゴリズムとしては、ガウスの消去法を利用します。また、問題を単純にするために引数で行列の次元を与えること、また逆行列が存在しない場合はnilを返すものとします。
ガウスの消去法は、逆行列を計算したい行列と、単位行列を並べ、元の行列を単位行列に変化させるという作業(掃出し)を行う際に、同じ操作を、並べてある単位行列に適用することによって実現します。
ここで、単位行列にするための操作は、原則として、次の三種類があります。
これらの操作は、実は、いずれも、可逆な作業であり、かつ、同時に、実は、ある種の行列の掛け算(可逆ということは、すなわち、逆行列を持つ)であることが、示されるわけですが、これに関しては、線形代数学に譲るとして、ここでは、アルゴリズムとして、これらの操作を考えます。
消去法では、まず、前進の掃出しと、後退の掃出しがあります。前進が終了すれば、後退は、無条件に行えます。
前進消去は、行の順番に行います。ここで、処理する行の対角要素が0の場合は、不都合があるので、他の行と交換する必要があります(枢軸選び)。
もし、都合の良い行が見つからなかった場合は、その行列は、正則でなかったことがわかり、そこでアルゴリズムは終了します。
全体の流れとしては、次のような手順になるでしょう。
逆行列の計算を行う関数mrevをいきなり作るのは大変なので簡単な関数から次第にくみ上げて行く形にします。このために、まずは、行列操作の基本関数から作ることにしましょう。
行列aのn行目横ベクトルを返す関数(nraw a n)
(nraw '((1 2 3) (3 5 6) (7 8 9)) 2)の結果は(3 5 6)
(nraw '((1 2 3) (3 5 6) (7 8 9)) 1)の結果は(1 2 3)
行列aのn行、m列の要素を返す関数(nmelm a n m)
(nmelm '((1 2 3) (3 5 6) (7 8 9)) 2 1)の結果は3
(nmelm '((1 2 3) (3 5 6) (7 8 9)) 1 2)の結果は2
行列aのn行目横ベクトルをvに変更した結果を返す関数(nreplace a n v)
(nreplace '((1 2 3) (3 5 6) (7 8 9)) 2 '(0 0 0))の結果は((1 2 3) (0 0 0) (7 8 9))
(nreplace '((1 2 3) (3 5 6) (7 8 9)) 3 '(0 0 0) )の結果は((1 2 3) (3 5 6) (0 0 0))
行列aのn行目をc倍する関数(ncprod a n c)
(ncprod '( (1 2 3) (4 5 6) (7 8 9) ) 2 3)の結果は((1 2 3) (12 15 18) (7 8 9))
行列aのn行目からm行目のc倍を引く(nmcsub a n m c )
(nmcsub '( (1 2 3) (4 5 6) (7 8 9) ) 3 1 7 )の結果は((1 2 3) (4 5 6) (0 -6 -12))
行列aのn行目とm行目を交換する(nmchange a n m)
(nmchange '((1 2 3)(4 5 6)(7 8 9)) 1 2)の結果は((4 5 6) (1 2 3) (7 8 9))
(nmchange '((1 2 3)(4 5 6)(7 8 9)) 3 2)の結果は((1 2 3) (7 8 9) (4 5 6))
以下の処理では、色々なところで、単位行列が必要になります。そこで、n次元の正方行列を求める関数があると便利です。ここでは、そのn次元の正方行列を求める関数(nunit n)を考えることにします。さて、正方行列の計算ですが、これは、よく考えると、単位ベクトルが並んでいるものと考えることができます。したがって、最初に、n次元の、i番目の単位ベクトルを返すを考えることにします。(nunitv n i)さて、n次元のi番目の単位ベクトルというのは、ようするにn個の要素を持つlistで、i番目が1で、その他の要素は、0となるようなベクトルです。従って、n個の要素を持つlistを、後ろから要素を決めながら、作って行くことになります。従って、その定義は次のようになります。
(defun nunitv (n i) (nunitvsub n i 1) ) (defun nunitvsub ( n i j ) (cond ( (< n j) nil ) ( t (cons (cond ( (= i j) 1 ) ( t 0 ) ) ( nunitvsub n i (+ j 1) ) )) ))
このnunitvがあれば、nunitも同様に実現することができます。
(defun nunit ( n ) (nunitsub n 1) ) (defun nunitsub ( n i ) (cond ( (< n i) nil ) ( t (cons (nunitv n i) (nunitsub n (+ i 1))) ) ))
さて、ほぼ、逆行列の計算の本体部分に手を付ける準備が殆どできあがってきたのですが、その前に、もう一つだけ、考えます。この逆行列の計算では、与えられたn次元の正方行列の要素を消去しながら、その操作を、同時に、元々単位行列だった行列に適用します。もし、元の行列が単位行列になったら、今度は、その単位行列に操作を加えたものが元の行列の逆行列になっているという点が、このアルゴリズムの基本です。従って、以下の関数では常に、元の行列と、元単位行列だった行列の二つの行列を対にして扱います。そこで、準備の最後として、二つのデータの対を扱う次の様な関数を作成します。
計算の順番としては、最後になりますが、簡単な方から、片づけることにして、後退消去から考えます。後退消去は、すでに、上三角行列になっており、更に、対角要素は、全て1となっている場合の処理です。原則としては、下の列の要素から、順に0に消去して行きます。その消去は、nmcsubで行うので、その処理は機械的です。今、n行の正方行列で、m列目の消去を行う関数を、
(gbcnm n m p)で表すとするとします。ただし、pは、元の正方行列と単位行列に、これまでの操作を加えた結果を要素とするデータ対です。これの定義は、次のようになります。
(defun gbcnm ( n m p ) ( gbcnmsub n m (+ m 1) p ) ) (defun gbcnmsub ( n m i p ) (cond ( (< n i) p ) ( t (gbcnmsub n m (+ i 1) (makepair (nmcsub (pair1 p) m i (nmelm (pair1 p) m i )) (nmcsub (pair2 p) m i (nmelm (pair1 p) m i )) )))))
この結果、例えば、
> (gbcnm 3 2 (makepair '( (1 2 3) (0 1 2) (0 0 1) ) (nunit 3) ) ) (((1 2 3) (0 1 0) (0 0 1)) ((1 0 0) (0 1 -2) (0 0 1))) > (gbcnm 3 1 (gbcnm 3 2 (makepair '((1 2 3) (0 1 2) (0 0 1)) (nunit 3)) )) (((1 0 0) (0 1 0) (0 0 1)) ((1 -2 1) (0 1 -2) (0 0 1))) >
となり、確かに、一つ目の要素が、単位行列になっていることがわかります。最後に、この計算結果の二つ目の要素と、最初の上三角行列の積を計算すると、確かに、単位行列になることがわかります。
> (mprod '( (1 2 3) (0 1 2) (0 0 1) ) '((1 -2 1) (0 1 -2) (0 0 1))) ((1 0 0) (0 1 0) (0 0 1)) >
これを利用して、n次元の行列のガウスの後退消去を行う(gbc n p)は、次のようになります(もちろん、pは、データ対です)。
(defun gbc ( n p ) ( gbcsub n (- n 1) p ) ) (defun gbcsub ( n i p ) (cond ( (< i 1) p ) ( t (gbcsub n (- i 1) (gbcnm n i p)) ) ))
これを使えば、
> (gbc 3 (makepair '( (1 2 3) (0 1 2) (0 0 1) ) (nunit 3) ) ) (((1 0 0) (0 1 0) (0 0 1)) ((1 -2 1) (0 1 -2) (0 0 1))) >
となり、後退消去がうまく行くことがわかります。
前進消去は、基本的に、後退消去と同じ操作を、今度は、上から下へと行い、上三角行列を作る作業です。したがって、ここでも、活躍するのは、nmcsubです。
しかし、前進消去が、後退消去と大きく異なる点は、後退消去が、無条件に行えるのに、前進消去は、場合によっては、行えなかったり(正則でないと、途中でできなくなります)。あるいは、操作を加える行を交換したりする(枢軸選び.. )が必要になるという点です。
この点を注意しながら、順番に作って行くことにしましょう。
まず、考えるのは、後退消去と同様、m行目の処理を行う関数
(gfcnm n m p)
です。
引数のn, m, pは、gbcと同じです。また、ここではpの一つ目に入っている消去の対象となる行列は、m - 1行目までは、すでに、上三角化されているものとして、処理します。
(defun gfcnm ( n m p ) ( gfcnmsub n m 1 p ) ) (defun gfcnmsub ( n m i p ) (cond ( (< n i) p ) ( t (gfcnmsub n m (+ i 1) (makepair (nmcsub (pair1 p) m i (nmelm (pair1 p) m i )) (nmcsub (pair2 p) m i (nmelm (pair1 p) m i )) )))))
これを利用すれば、例えば、次のような結果が得られます。
> (gfcnm 3 3 (makepair '( (1 2 3) (0 1 6) (7 8 9) ) (nunit 3) ) ) (((1 2 3) (0 1 6) (0 0 -552)) ((1 0 0) (0 1 0) (161 -138 -23))) >
次に考えるのは、正規化と、枢軸選びの問題です。
正規化は単に、対角要素を1にする操作で、これは、単に、対角要素を含む行を、その対角要素で割ることで実現します。但し、問題は、たまたま、その対角要素が0の場合です。この場合は、当然のことながら、0で割るということはできないので、困ってしまうわけです。もし、対角要素が0であれば、対角要素が0にならないように、他の行と交換を行うことにします。これが枢軸選びです。そこで、与えられた行列から、必要に応じて、枢軸を選び、正規化して、返すような関数( gfcnormal n m p )を考えることにします。ここで、もし、枢軸が選べない場合は、nilを返すことにします。
(defun gfcnormal ( n i p ) (cond ( (zerop (nmelm (pair1 p) i i)) (gfcnormalsub n i (gfcfind n i p)) ) ( t (gfcnormalsub n i p) ) ) ) (defun gfcnormalsub ( n i p ) (cond ( (null p) nil) ( t (makepair (ncprod (pair1 p) i (/ 1 (nmelm (pair1 p) i i))) (ncprod (pair2 p) i (/ 1 (nmelm (pair1 p) i i))) ) ) ) ) (defun gfcfind (n i p) (gfcfindsub n i (+ i 1) p) ) (defun gfcfindsub (n i j p) (cond ( (< n j) nil ) ( (zerop (nmelm (pair1 p) j i ) ) (gfcfindsub n i (+ j 1) p) ) ( t (makepair (nmchange (pair1 p) i j) (nmchange (pair2 p) i j) ) ) ))
ここで、gfcfindは、枢軸を探して、見つかれば、行を交換したものを返します。もし、枢軸にふさわしい行がなければ、nilが帰ります。gfcnormalsubは、枢軸が与えられていれば、正規化を行う関数です。
さて、いよいよ、前進消去が行える準備ができました。前進消去と、後退消去の違いは、単に、途中で、正規化に失敗した場合に、値としてnilを返すというだけのことです。定義は、以下のようになります。
(defun gfc ( n p ) ( gfcsub n 2 (gfcnormal n 1 p) ) ) (defun gfcsub ( n i p ) (cond ( (null p) nil ) ( (< n i) p ) ( t (gfcsub n (+ i 1) (gfcnormal n i (gfcnm n i p) ) ) ) ))
以上で準備はおしまいです。
逆行列の計算を行う関数mrevは、次のような感じになります。
(defun mrev ( n m ) (mrevsub n (makepair m (nunit n))) ) (defun mrevsub ( n p ) (mrevsub2 n (gfc n p) ) ) (defun mrevsub2 ( n p ) (cond ( (null p) nil ) ( t (pair2 (gbc n p)) ) ) )
以上の結果をまとめたものがmatrix.txtです。