M.Hiroi's Home Page

Functional Programming

お気楽 Haskell プログラミング入門

[ PrevPage | Haskell | NextPage ]

線形代数編

●特殊な行列とその性質

Prelude Numeric.LinearAlgebra> a = matrix 2 [1,2,3,4]
Prelude Numeric.LinearAlgebra> a
(2><2)
 [ 1.0, 2.0
 , 3.0, 4.0 ]
Prelude Numeric.LinearAlgebra> tr a
(2><2)
 [ 1.0, 3.0
 , 2.0, 4.0 ]
Prelude Numeric.LinearAlgebra> mTm a
Herm (2><2)
 [ 10.0, 14.0
 , 14.0, 20.0 ]
Prelude Numeric.LinearAlgebra> (tr a) <> a
(2><2)
 [ 10.0, 14.0
 , 14.0, 20.0 ]

Prelude Numeric.LinearAlgebra> b = matrix 3 [1..9]
Prelude Numeric.LinearAlgebra> b
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> tr b
(3><3)
 [ 1.0, 4.0, 7.0
 , 2.0, 5.0, 8.0
 , 3.0, 6.0, 9.0 ]
Prelude Numeric.LinearAlgebra> mTm b
Herm (3><3)
 [ 66.0,  78.0,  90.0
 , 78.0,  93.0, 108.0
 , 90.0, 108.0, 126.0 ]
Prelude Numeric.LinearAlgebra> (tr b) <> b
(3><3)
 [ 66.0,  78.0,  90.0
 , 78.0,  93.0, 108.0
 , 90.0, 108.0, 126.0 ]

Prelude Numeric.LinearAlgebra> c = matrix 2 [1..6]
Prelude Numeric.LinearAlgebra> c
(3><2)
 [ 1.0, 2.0
 , 3.0, 4.0
 , 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> tr c
(2><3)
 [ 1.0, 3.0, 5.0
 , 2.0, 4.0, 6.0 ]
Prelude Numeric.LinearAlgebra> mTm c
Herm (2><2)
 [ 35.0, 44.0
 , 44.0, 56.0 ]
Prelude Numeric.LinearAlgebra> (tr c) <> c
(2><2)
 [ 35.0, 44.0
 , 44.0, 56.0 ]
Prelude Numeric.LinearAlgebra> a
(2><2)
 [ 1.0, 2.0
 , 3.0, 4.0 ]
Prelude Numeric.LinearAlgebra> trace a
5.0
Prelude Numeric.LinearAlgebra> trace $ tr a
5.0
Prelude Numeric.LinearAlgebra> det a
-2.0
Prelude Numeric.LinearAlgebra> det $ tr a
-2.0

Prelude Numeric.LinearAlgebra> d = matrix 2 [5,6,7,8]
Prelude Numeric.LinearAlgebra> d
(2><2)
 [ 5.0, 6.0
 , 7.0, 8.0 ]
Prelude Numeric.LinearAlgebra> tr (a <> d)
(2><2)
 [ 19.0, 43.0
 , 22.0, 50.0 ]
Prelude Numeric.LinearAlgebra> (tr d) <> (tr a)
(2><2)
 [ 19.0, 43.0
 , 22.0, 50.0 ]

Prelude Numeric.LinearAlgebra> inv (tr a)
(2><2)
 [ -2.0,  1.5
 ,  1.0, -0.5 ]
Prelude Numeric.LinearAlgebra> tr (inv a)
(2><2)
 [ -1.9999999999999998,   1.4999999999999998
 ,                 1.0, -0.49999999999999994 ]
Prelude Numeric.LinearAlgebra> a = matrix 2 [1,3,3,2]
Prelude Numeric.LinearAlgebra> a
(2><2)
 [ 1.0, 3.0
 , 3.0, 2.0 ]
Prelude Numeric.LinearAlgebra> tr a == a
True
Prelude Numeric.LinearAlgebra> b = matrix 3 [1, 4, 5, 4, 2, 6, 5, 6, 3]
Prelude Numeric.LinearAlgebra> b
(3><3)
 [ 1.0, 4.0, 5.0
 , 4.0, 2.0, 6.0
 , 5.0, 6.0, 3.0 ]
Prelude Numeric.LinearAlgebra> tr b == b
True
Prelude Numeric.LinearAlgebra> c = matrix 3 [1,0,0,0,2,0,0,0,3]
Prelude Numeric.LinearAlgebra> c
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 2.0, 0.0
 , 0.0, 0.0, 3.0 ]
Prelude Numeric.LinearAlgebra> tr c == c
True
Prelude Numeric.LinearAlgebra> a = matrix 3 [1..9]
Prelude Numeric.LinearAlgebra> a
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> sym a
Herm (3><3)
 [ 1.0, 3.0, 5.0
 , 3.0, 5.0, 7.0
 , 5.0, 7.0, 9.0 ]

Prelude Numeric.LinearAlgebra> takeDiag a
[1.0,5.0,9.0]
Prelude Numeric.LinearAlgebra> diag $ vector [1,5,9]
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 5.0, 0.0
 , 0.0, 0.0, 9.0 ]
Prelude Numeric.LinearAlgebra> diagl [1,5,9]
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 5.0, 0.0
 , 0.0, 0.0, 9.0 ]

Prelude Numeric.LinearAlgebra> diagRect 1 (vector [2, 3, 4]) 3 3
(3><3)
 [ 2.0, 1.0, 1.0
 , 1.0, 3.0, 1.0
 , 1.0, 1.0, 4.0 ]
Prelude Numeric.LinearAlgebra> diagRect 1 (vector [2, 3, 4]) 3 4
(3><4)
 [ 2.0, 1.0, 1.0, 1.0
 , 1.0, 3.0, 1.0, 1.0
 , 1.0, 1.0, 4.0, 1.0 ]
Prelude Numeric.LinearAlgebra> diagRect 1 (vector [2, 3, 4]) 4 3
(4><3)
 [ 2.0, 1.0, 1.0
 , 1.0, 3.0, 1.0
 , 1.0, 1.0, 4.0
 , 1.0, 1.0, 1.0 ]
Prelude Numeric.LinearAlgebra> a = diagl [1,2,3]
Prelude Numeric.LinearAlgebra> a
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 2.0, 0.0
 , 0.0, 0.0, 3.0 ]
Prelude Numeric.LinearAlgebra> det a
6.0
Prelude Numeric.LinearAlgebra> :t prodElements
prodElements :: Container c e => c e -> e
Prelude Numeric.LinearAlgebra> prodElements (takeDiag a :: Vector R)
6.0

Prelude Numeric.LinearAlgebra> a <> a
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 4.0, 0.0
 , 0.0, 0.0, 9.0 ]
Prelude Numeric.LinearAlgebra> a <> a <> a
(3><3)
 [ 1.0, 0.0,  0.0
 , 0.0, 8.0,  0.0
 , 0.0, 0.0, 27.0 ]

Prelude Numeric.LinearAlgebra> inv a
(3><3)
 [ 1.0, 0.0,                0.0
 , 0.0, 0.5,                0.0
 , 0.0, 0.0, 0.3333333333333333 ]
Prelude Numeric.LinearAlgebra> diag (1 / vector [1,2,3])
(3><3)
 [ 1.0, 0.0,                0.0
 , 0.0, 0.5,                0.0
 , 0.0, 0.0, 0.3333333333333333 ]
Prelude Numeric.LinearAlgebra> eigenvalues a
[1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0]
Prelude Numeric.LinearAlgebra> a = matrix 3 [1,2,3,4,5,6,7,8,10]
Prelude Numeric.LinearAlgebra> a
(3><3)
 [ 1.0, 2.0,  3.0
 , 4.0, 5.0,  6.0
 , 7.0, 8.0, 10.0 ]
Prelude Numeric.LinearAlgebra> la = matrix 3 [1,0,0,4,5,0,7,8,10]
Prelude Numeric.LinearAlgebra> la
(3><3)
 [ 1.0, 0.0,  0.0
 , 4.0, 5.0,  0.0
 , 7.0, 8.0, 10.0 ]
Prelude Numeric.LinearAlgebra> ua = matrix 3 [1,2,3,0,5,6,0,0,10]
Prelude Numeric.LinearAlgebra> ua
(3><3)
 [ 1.0, 2.0,  3.0
 , 0.0, 5.0,  6.0
 , 0.0, 0.0, 10.0 ]

Prelude Numeric.LinearAlgebra> tr la
(3><3)
 [ 1.0, 4.0,  7.0
 , 0.0, 5.0,  8.0
 , 0.0, 0.0, 10.0 ]
Prelude Numeric.LinearAlgebra> tr ua
(3><3)
 [ 1.0, 0.0,  0.0
 , 2.0, 5.0,  0.0
 , 3.0, 6.0, 10.0 ]

Prelude Numeric.LinearAlgebra> la <> la
(3><3)
 [   1.0,   0.0,   0.0
 ,  24.0,  25.0,   0.0
 , 109.0, 120.0, 100.0 ]
Prelude Numeric.LinearAlgebra> ua <> ua
(3><3)
 [ 1.0, 12.0,  45.0
 , 0.0, 25.0,  90.0
 , 0.0,  0.0, 100.0 ]

Prelude Numeric.LinearAlgebra> det la
49.99999999999999
Prelude Numeric.LinearAlgebra> det ua
50.0

Prelude Numeric.LinearAlgebra> inv la
(3><3)
 [     1.0000000000000002,                  0.0, -3.96508223080413e-18
 ,                   -0.8,  0.20000000000000004, -2.42861286636753e-17
 , -6.0000000000000046e-2, -0.16000000000000003,   0.10000000000000002 ]
Prelude Numeric.LinearAlgebra> inv ua
(3><3)
 [ 1.0, -0.4,              -6.0e-2
 , 0.0,  0.2, -0.12000000000000002
 , 0.0,  0.0,                  0.1 ]

Prelude Numeric.LinearAlgebra> eigenvalues la
[10.0 :+ 0.0,5.0 :+ 0.0,1.0 :+ 0.0]
Prelude Numeric.LinearAlgebra> eigenvalues ua
[1.0 :+ 0.0,5.0 :+ 0.0,10.0 :+ 0.0]
回転行列 R は直交行列

R = [ cos(x), -sin(x),
      sin(x),  cos(x) ]

転置行列
RT = [  cos(x), sin(x),
       -sin(x), cos(x) ]

行列式 |R| = cos(x) * cos(x) + sin(x) * sin(x) = 1

逆行列
R-1 = [  cos(x), sin(x),
        -sin(x), cos(x) ]          

逆行列の転置 (R に戻るので直交行列)
(R-1)T = [ cos(x), -sin(x),
           sin(x),  cos(x) ]
Prelude Numeric.LinearAlgebra> a = matrix 2 [0,1,1,0]
Prelude Numeric.LinearAlgebra> a
(2><2)
 [ 0.0, 1.0
 , 1.0, 0.0 ]
Prelude Numeric.LinearAlgebra> tr a
(2><2)
 [ 0.0, 1.0
 , 1.0, 0.0 ]
Prelude Numeric.LinearAlgebra> det a
-1.0
Prelude Numeric.LinearAlgebra> inv a
(2><2)
 [ 0.0, 1.0
 , 1.0, 0.0 ]

Prelude Numeric.LinearAlgebra> b = matrix 3 [0, 1, 0, 1, 0, 0, 0, 0, 1]
Prelude Numeric.LinearAlgebra> b
(3><3)
 [ 0.0, 1.0, 0.0
 , 1.0, 0.0, 0.0
 , 0.0, 0.0, 1.0 ]
Prelude Numeric.LinearAlgebra> tr b
(3><3)
 [ 0.0, 1.0, 0.0
 , 1.0, 0.0, 0.0
 , 0.0, 0.0, 1.0 ]
Prelude Numeric.LinearAlgebra> inv b
(3><3)
 [ 0.0, 1.0, 0.0
 , 1.0, 0.0, 0.0
 , 0.0, 0.0, 1.0 ]

Prelude Numeric.LinearAlgebra> c = matrix 3 [1..9]
Prelude Numeric.LinearAlgebra> c
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> c <> b
(3><3)
 [ 2.0, 1.0, 3.0
 , 5.0, 4.0, 6.0
 , 8.0, 7.0, 9.0 ]
Prelude Numeric.LinearAlgebra> b <> c
(3><3)
 [ 4.0, 5.0, 6.0
 , 1.0, 2.0, 3.0
 , 7.0, 8.0, 9.0 ]

●連立一次方程式の解法

今回は連立一次方程式を解く基本的なアルゴリズムを取り上げます。hmatrix には連立一次方程式 Ax = b の x を求める関数 (ソルバー) が用意されているので、実用的にはそちらを使えばいいのですが、Haskell (hmatrix) とアルゴリズムのお勉強ということで、あえてプログラムを作ってみましょう。

●ガウスの消去法

ガウスの消去法 (Gaussian elimination) は人が連立方程式を解くときの方法とほとんど同じです。次の図を見てください。

a1 * x + a2 * y + a3 * z = d1       a1 * x + a2  * y + a3  * z = d1        a1 * x + a2  * y + a3   * z = d1
b1 * x + b2 * y + b3 * z = d2  ==>           b2' * y + b3' * z = d2'  ==>           b2' * y + b3'  * z = d2'
c1 * x + c2 * y + c3 * z = d3                c2' * y + c3' * z = d3'                          c3'' * z = d3''

                                          図 : 前進消去

1 番目の式を b1 / a1 倍して 2 番目の式から引き算すると、x の項を消去することができます。同様の方法で 3 番目の式の x の項を消去することができます。次に、2 番目の式を c2' / b2' 倍して 3 番目の式から引き算すると y の項を消去することができます。この処理を「前進消去」といいます。

前進消去を行うと、最後の式には変数が z しかありません。z の値は d3'' / c3'' に決定することができます。次はひとつ前の式に戻り、 変数 z に d3'' / c3'' を代入すると変数 y の値を求めることができます。これを繰り返して、先頭の式に戻ると変数 x の値を求めることができます。この処理を「後退代入」といいます。

連立一次方程式は行列を使って Ax = b と表すことができます。A を係数行列といい、A と b を連結した行列を拡大係数行列といいます。

[ a1, a2, a3     [ a1, a2, a3, d1
, b1, b2, b3     , b1, b2, b3, d2
, c1, c2, c3 ]   , c1, c2, c3, d3 ]

   係数行列         拡大係数行列

ガウスの消去法は、拡大係数行列を使うと簡単にプログラムを作ることができます。Haskell でプログラムすると、次のようになります。

リスト : ガウスの消去法 (前進消去)

import Prelude hiding ((<>))
import Numeric.LinearAlgebra

gauss1 :: [Vector R] -> Int -> [Vector R]
gauss1 [v] _ = [v]
gauss1 (v:vs) n =
  v : gauss1 (map (\w -> let temp = (w ! n) / (v ! n)
                         in w - (v * scalar temp)) vs)
             (n + 1)

関数 gauss1 は前進消去を行います。第 1 引数は拡大係数行列を行単位で分解したリストで、要素は行を表すベクトルです。第 2 引数が係数を 0 にするインデックスを表します。リストの要素が一つしかない場合、一番最後の式に到達しました。0 から n - 1 番目までの係数は 0 になっているので、v をそのままリストに格納して返します。

そうでなければ、リストを v:vs で分解して、vs の n 番目の係数を 0 にします。これは map を使えば簡単です。ベクトル w の n 番目の要素を 0 にするので、wn / vn を v に乗算し、それを w から減算すれば、w の n 番目の係数を 0 にすることができます。vs に格納されたベクトルの n 番目の要素を 0 にしたら、gauss1 を再帰呼び出しして、n + 1 番目の係数を 0 にします。

簡単な実行例を示します。

 x +  y = 100
2x + 4y = 272
*Main> fromRows $ gauss1 [vector [1,1,100], vector [2,4,272]] 0
(2><3)
 [ 1.0, 1.0, 100.0
 , 0.0, 2.0,  72.0 ]
 x +  y +  z = 10
2x + 4y + 6z = 38
2x      + 4z = 14
*Main> fromRows $ gauss1 [vector [1,1,1,10], vector [2,4,6,38], vector [2,0,4,14]] 0
(3><4)
 [ 1.0, 1.0, 1.0, 10.0
 , 0.0, 2.0, 4.0, 18.0
 , 0.0, 0.0, 6.0, 12.0 ]
  a +  b +   c + d = -5
 -a +  b -   c + d = -7
 8a + 4b +  2c + d = -31
-8a + 4b + -2c + d = -35
*Main> fromRows $ gauss1 [vector[1,1,1,1,-5], vector[-1,1,-1,1,-7], vector[8,4,2,1,-31], vector[-8,4,-2,1,-35]] 0
(4><5)
 [ 1.0, 1.0,  1.0,  1.0,  -5.0
 , 0.0, 2.0,  0.0,  2.0, -12.0
 , 0.0, 0.0, -6.0, -3.0, -15.0
 , 0.0, 0.0,  0.0, -6.0, -18.0 ]

次は後退代入を行う関数 gauss2 を作ります。

リスト : ガウスの消去法 (後退代入)

gauss2 :: [Vector R] -> Int -> [R]
gauss2 [v] n = [(v ! (n + 1)) / (v ! n)]
gauss2 (v:vs) n =
  let m = n + 1
      x = gauss2 vs m
      y = subVector m (size v - 1 - m) v <.> vector x
      z = ((v ! (size v - 1)) - y)  / (v ! n)
  in z : x

第 1 引数のリストには gauss1 の返り値を渡します。第 2 引数が求める変数のインデックスを表します。返り値は変数の値を格納したリストです。リストの要素が一つの場合、ベクトル v には変数が一つしかないので、その値は簡単に求めることができます。v ! (n + 1) の値を v ! n で割り算するだけです。これをリストに格納して返します。

そうでなければ、リストを v と vs に分解し、vs に対して gauss2 を再帰呼び出しします。これで n + 1 (= m) 番目から最後までの変数の値を求めることができます。次に、v の m 番目から最後の変数に値を代入して計算し、変数 y にセットします。これは内積 <.> を使えば簡単ですね。subVector で部分ベクトルを取り出していることに注意してください。あとは、v の最後の値から y を引き算し、v ! n で割り算すれば n 番目の変数の値を求めることができます。

最後にガウスの消去法で連立方程式を解く関数 gauss を作ります。

リスト : ガウスの消去法

gauss :: Matrix R -> Vector R -> [R]
gauss xs ys =
  let zs = xs ||| asColumn ys
  in gauss2 (gauss1 (toRows zs) 0) 0

関数 gauss の引数 xs が係数行列、ys が右辺の定数を格納したベクトルです。まず最初に、拡大係数行列を生成して変数 zs にセットします。あとは toRows で zs を分解して、gauss1 と gauss2 を呼び出すだけです。

それでは実際に試してみましょう。

 x +  y = 100
2x + 4y = 272
*Main> gauss (matrix 2 [1,1,2,4]) (vector [100,272])
[64.0,36.0]
 x +  y +  z = 10
2x + 4y + 6z = 38
2x      + 4z = 14
*Main> gauss (matrix 3 [1,1,1,2,4,6,2,0,4]) (vector [10,38,14])
[3.0,5.0,2.0]
  a +  b +   c + d = -5
 -a +  b -   c + d = -7
 8a + 4b +  2c + d = -31
-8a + 4b + -2c + d = -35
*Main> gauss (matrix 4 [1,1,1,1,-1,1,-1,1,8,4,2,1,-8,4,-2,1]) (vector [-5,-7,-31,-35])
[0.0,-9.0,1.0,3.0]
  a -  b +  c - d + e = 1
12a - 6b + 2c         = 0 
  a +  b +  c + d + e = 8
12a + 6b + 2c         = 0
 4a + 3b + 2c + d     = 1
*Main> gauss (matrix 5 [1,-1,1,-1,1,12,-6,2,0,0,1,1,1,1,1,12,6,2,0,0,4,3,2,1,0]) (vector [1,0,8,0,1])
[0.3124999999999991,0.0,-1.874999999999997,3.5,6.062499999999997]

●ガウス・ジョルダン法

ガウスの消去法は、係数行列の部分を「上三角行列」に変換することで、連立一次方程式を解きました。ここで、係数行列を単位行列に変換できれば、後退代入することなく解を求めることができます。これを「ガウス・ジョルダン法 (Gauss - Jordan elimination)」といいます。簡単な動作例を下図に示します。

次の連立方程式をガウス・ジョルダン法で解く

  x +  y +  z = 10
 2x + 4y + 6z = 38
 2x      + 4z = 14

拡大係数行列を zs とする

(1) zs[1, 1] を 1 にする
[  1.   1.   1.  10.     # zs[1, 1:end] /= 1
   2.   4.   6.  38.
   2.   0.   4.  14.]

(2) zs[2, 1], zs[3, 1] を 0 にする
[  1.   1.   1.  10.      
   0.   2.   4.  18.      # zs[2, 1:end] -= zs[1, 1:end] * 2
   0.  -2.   2.  -6.]     # zs[3, 1:end] -= zs[1, 1:end] * 2

(3) zs[2, 2] を 1 にする
[  1.   1.   1.  10. 
   0.   1.   2.   9.      # zs[2, 2:end] /= 2
   0.  -2.   2.  -6.]

(4) zs[1, 2], zs[3, 2] を 0 にする
[  1.   0.  -1.   1.      # zs[1, 2:end] -= zs[2, 2:end] * 1
   0.   1.   2.   9. 
   0.   0.   6.  12.]     # zs[3, 2:end] -= zs[2, 2:end] * -2

(5) zs[3, 3] を 1 にする
[  1.   0.  -1.   1. 
   0.   1.   2.   9. 
   0.   0.   1.   2.]     # zs[3, 3:end] /= 6

(6) zs[1, 3], zs[2, 3] を 0 にする
[  1.   0.   0.   3.      # zs[1, 3:end] -= zs[3, 3:end] * -1
   0.   1.   0.   5.      # zs[2, 3:end] -= zs[3, 3:end] * 2
   0.   0.   1.   2.]

このように、ガウス・ジョルダン法の基本的な考え方は実にシンプルで、プログラムも簡単になるのですが、速度はガウスの消去法よりも少し遅くなるといわれています。以下にプログラムと実行例を示します。

リスト : ガウス・ジョルダン法

gj1 :: [Vector R] -> Int -> Int -> [Vector R]
gj1 zs n m
  | n == m    = zs
  | otherwise = gj1 (gj2 zs 0) (n + 1) m
      where
        v = zs !! n
        w = v / scalar (v ! n)
        gj2 [] _ = []
        gj2 (x:xs) i =
          (if i == n then w else x - (w * scalar (x ! n))) : gj2 xs (i + 1)

gaussJordan :: Matrix R -> Vector R -> [R]
gaussJordan xs ys =
  let zs = xs ||| asColumn ys
  in map (\v -> v ! (size v - 1)) $ gj1 (toRows zs) 0 (rows xs)

係数拡大行列の操作は関数 gj1 で行います。gj1 は係数拡大行列を返すので、それから行の末尾要素を map で取り出します。gj1 の実際の処理は局所関数 gj2 で行います。再帰呼び出しで n 行以外の行の n 番目の要素を 0 にしています。

 x +  y = 100
2x + 4y = 272
*Main> gaussJordan (matrix 2 [1,1,2,4]) (vector [100,272])
[64.0,36.0]
 x +  y +  z = 10
2x + 4y + 6z = 38
2x      + 4z = 14
*Main> gaussJordan (matrix 3 [1,1,1,2,4,6,2,0,4]) (vector [10,38,14])
[3.0,5.0,2.0]
  a +  b +   c + d = -5
 -a +  b -   c + d = -7
 8a + 4b +  2c + d = -31
-8a + 4b + -2c + d = -35
*Main> gaussJordan (matrix 4 [1,1,1,1,-1,1,-1,1,8,4,2,1,-8,4,-2,1]) (vector [-5,-7,-31,-35])
[0.0,-9.0,1.0,3.0]
  a -  b +  c - d + e = 1
12a - 6b + 2c         = 0 
  a +  b +  c + d + e = 8
12a + 6b + 2c         = 0
 4a + 3b + 2c + d     = 1
*Main> gaussJordan (matrix 5 [1,-1,1,-1,1,12,-6,2,0,0,1,1,1,1,1,12,6,2,0,0,4,3,2,1,0]) (vector [1,0,8,0,1])
[0.3125,8.881784197001252e-16,-1.875,3.5,6.0625]

●ピボット選択

ところで、ガウスの消去法とガウス・ジョルダン法は、拡大係数行列の対角成分 zs[i, i] が 0 になると計算できなくなります。たとえば、次の連立方程式は解くことができますが、zs[0, 0] が 0 になっているため、プログラムの実行中に値が無限大になってしまいます。

      2y + 4z = 14   [ 0.  2.  4. 14.
  x +  y +  z = 10     1.  1.  1. 10.
 4x + 2y + 6z = 38     4.  2.  6. 38. ]

    連立方程式        拡大係数行列 zs

このとき、拡大係数行列の行を交換すると連立一次方程式が解ける場合があります。また、zs[i, i] が 0 に等しくなくても 0 に近い値だと解の誤差が大きくなります。そこで、zs[i, i] の絶対値がなるべく大きくなるように行を交換します。これを「部分ピボット選択」といいます。なお、参考文献 1 によると、『行だけでなく列も交換する完全ピボット選択という方法もあるが、通常は部分ピボット選択で十分である』 とのことです。

プログラムは次のようになります。

リスト : 部分ピボット選択 (1)

findMaxRow :: [Vector R] -> Int -> Int
findMaxRow [] _ = 0
findMaxRow (z:zs) n =
  iter zs 1 (abs (z ! n)) 0
  where
    iter [] _ _ b = b
    iter (z:zs) i a b =
      let c = abs (z ! n)
      in if a < c
         then iter zs (i + 1) c i
         else iter zs (i + 1) a b

swapRow :: [a] -> Int -> Int -> [a]
swapRow zs m n =
  iter zs 0
  where
    x = zs !! m
    y = zs !! n
    iter [] _ = []
    iter (z:zs) i
      | i == m    = y : iter zs (i + 1)
      | i == n    = x : iter zs (i + 1)
      | otherwise = z : iter zs (i + 1)

selectPivot :: [Vector R] -> Int -> Int -> [Vector R]
selectPivot zs m n =
  let i = m + findMaxRow (drop m zs) n
  in if m == i
     then zs
     else swapRow zs m i

関数 findMaxRow は zs の中から n 番目の要素の絶対値が最大のベクトルを探します。実際の処理は局所関数 iter で行い、リスト zs を線形探索するだけです。関数 swapRow はリスト zs の m 番目と n 番目の要素を交換したリストを返します。実際の処理は局所関数 iter で行います。あらかじめ交換する要素を変数 x, y に求めておき、インデックス i が m と等しければ y を、n と等しければ x を、それ以外はそのままの値 (z) をコンス演算子 : で追加します。

ピボット選択は関数 selectPivot で行います。リスト zs の m 行以降で最大値の行を findMaxRow で選び、変数 i にセットします。i と m が同じであれば、m 行目が最大値なので何もせずに zs を返します。そうでなければ、swapRow で zs の m 行と i 行を交換します。

あとは selectPivot を呼び出すだけです。次のリストを見てください。

リスト : 部分ピボット選択 (ガウスの消去法)

-- 前進消去
gauss1' :: [Vector R] -> Int -> [Vector R]
gauss1' [v] _ = [v]
gauss1' (v:vs) n =
  v : gauss1' (selectPivot
               (map (\w -> let temp = (w ! n) / (v ! n)
                           in w - (v * scalar temp)) vs)
               0
               (n + 1))
              (n + 1)

gauss' :: Matrix R -> Vector R -> [R]
gauss' xs ys =
  let zs = xs ||| asColumn ys
  in gauss2 (gauss1' (selectPivot (toRows zs) 0 0) 0) 0

関数 gauss1' は前進消去でピボット選択を行います。第 1 引数にはピボット選択を行ったリストを渡します。gauss1' を再帰呼び出しするとき、map の返り値 (リスト) に selectPivot を適用します。この場合、リストの先頭と最大値の行を交換することになります。関数 gauss' は gauss1' を呼び出すとき selectPivot を適用して、先頭行 (0 番目) と最大値の行を交換します。

リスト : 部分ピボット選択 (ガウスのガウス・ジョルダン法)

gj1' :: [Vector R] -> Int -> Int -> [Vector R]
gj1' zs n m
  | n == m    = zs
  | otherwise = gj1' (selectPivot (gj2 zs 0) (n + 1) (n + 1)) (n + 1) m
      where
        v = zs !! n
        w = v / scalar (v ! n)
        gj2 [] _ = []
        gj2 (x:xs) i =
          (if i == n then w else x - (w * scalar (x ! n))) : gj2 xs (i + 1)

gaussJordan' :: Matrix R -> Vector R -> [R]
gaussJordan' xs ys =
  let zs = xs ||| asColumn ys
  in map (\v -> v ! (size v - 1)) $ gj1' (selectPivot (toRows zs) 0 0) 0 (rows xs)

関数 gaussJordan' も gauss' と同様に、関数 gj1' を呼び出すとき selectPivot を適用して、0 行目と最大値の行を交換します。関数 gj1' では、局所関数 gj2 の返り値に selectPivot を適用して、n + 1 行とそれ以降の最大値の行を交換します。gauss1' と違って、交換する行は先頭 (0 番目) ではなく n + 1 番目であることに注意してください。

それでは実行してみましょう。

*Main> gauss (matrix 3 [0,2,4,1,1,1,4,2,6]) (vector [14,10,38])
[NaN,NaN,NaN]
*Main> gauss' (matrix 3 [0,2,4,1,1,1,4,2,6]) (vector [14,10,38])
[5.0,3.0,2.0]
*Main> gauss' (matrix 4 [1,1,1,1,-1,1,-1,1,8,4,2,1,-8,4,-2,1]) (vector [-5,-7,-31,-35])
[0.0,-9.0,1.0,3.0]
*Main> gaussJordan (matrix 3 [0,2,4,1,1,1,4,2,6]) (vector [14,10,38])
[NaN,NaN,NaN]
*Main> gaussJordan' (matrix 3 [0,2,4,1,1,1,4,2,6]) (vector [14,10,38])
[5.0,3.0,2.0]
*Main> gaussJordan' (matrix 4 [1,1,1,1,-1,1,-1,1,8,4,2,1,-8,4,-2,1]) (vector [-5,-7,-31,-35])
[0.0,-9.0,1.0,3.0]

正しく解を求めることができました。

●逆行列

連立方程式 Ax = b を解く場合、A の逆行列 A-1 がわかれば、x = A-1b で解を求めることができます。実をいうと、逆行列はガウス・ジョルダン法を使って簡単に求めることができるのです。次の図を見てください。

逆行列の定義 : A * A-1 = I

[ a11 a12 a13        [ b11 b12 b13        [ 1 0 0 
  a21 a22 a23    *     b21 b22 b23     =    0 1 0 
  a31 a32 a33 ]        b31 b32 b33 ]        0 0 1 ]

   係数行列               逆行列         単位行列 I

上記の定義により逆行列 A-1 は次の方程式を解けば求めることができる

A * [b11, b21, b31] = [1, 0, 0]
A * [b12, b22, b32] = [0, 1, 0]
A * [b13, b23, b33] = [0, 0, 1]

これは次の拡大係数行列にガウス・ジョルダン法を適用すること同じ

[ a11 a12 a13  1  0  0 
  a21 a22 a23  0  1  0 
  a31 a32 a33  0  0  1 ]

拡大係数行列の右半分に逆行列が求まる

プログラムも簡単です。次のリストを見てください。

リスト : ガウス・ジョルダン法で逆行列を求める

matinv :: Matrix R -> Matrix R
matinv xs =
  let
    n = rows xs
    zs = xs ||| ident n
  in
    dropColumns n $ fromRows $ gj1' (selectPivot (toRows zs) 0 0) 0 n
end

単位行列は関数 ident n で作成することができます。あとはガウス・ジョルダン法で解を求め、最後に拡大係数行列の右半分を返します。それでは実行してみましょう。

*Main> a = matrix 2 [1,1,2,4]
*Main> a
(2><2)
 [ 1.0, 1.0
 , 2.0, 4.0 ]
*Main> ra = matinv a
*Main> ra
(2><2)
 [  2.0, -0.5
 , -1.0,  0.5 ]
*Main> a <> ra
(2><2)
 [ 1.0, 0.0
 , 0.0, 1.0 ]

*Main> b = matrix 3 [1,1,1,2,4,6,2,0,4]
*Main> b
(3><3)
 [ 1.0, 1.0, 1.0
 , 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0 ]
*Main> rb = matinv b
*Main> rb
(3><3)
 [  1.3333333333333333, -0.3333333333333333, 0.16666666666666669
 ,  0.3333333333333333, 0.16666666666666669, -0.3333333333333333
 , -0.6666666666666666, 0.16666666666666666, 0.16666666666666666 ]
*Main> disp 3 rb
3x3
 1.333  -0.333   0.167
 0.333   0.167  -0.333
-0.667   0.167   0.167
*Main> disp 3 $ b <> rb
3x3
1.000  0.000  0.000
0.000  1.000  0.000
0.000  0.000  1.000

*Main> c
(4><4)
 [  1.0, 1.0,  1.0, 1.0
 , -1.0, 1.0, -1.0, 1.0
 ,  8.0, 4.0,  2.0, 1.0
 , -8.0, 4.0, -2.0, 1.0 ]
*Main> rc = matinv c
*Main> disp 4 rc
4x4
-0.1667   0.1667   0.0833  -0.0833
-0.1667  -0.1667   0.1667   0.1667
 0.6667  -0.6667  -0.0833   0.0833
 0.6667   0.6667  -0.1667  -0.1667
*Main> disp 4 $ c <> rc
4x4
 1.0000   0.0000   0.0000  -0.0000
 0.0000   1.0000   0.0000   0.0000
-0.0000   0.0000   1.0000  -0.0000
 0.0000  -0.0000  -0.0000   1.0000

*Main> d = matrix 5 [1,-1,1,-1,1,12,-6,2,0,0,1,1,1,1,1,12,6,2,0,0,4,3,2,1,0]
*Main> d
(5><5)
 [  1.0, -1.0, 1.0, -1.0, 1.0
 , 12.0, -6.0, 2.0,  0.0, 0.0
 ,  1.0,  1.0, 1.0,  1.0, 1.0
 , 12.0,  6.0, 2.0,  0.0, 0.0
 ,  4.0,  3.0, 2.0,  1.0, 0.0 ]
*Main> rd = matinv d
*Main> disp 4 rd
5x5
-0.0625   0.0417   0.0625   0.0833  -0.1250
 0.0000  -0.0833   0.0000   0.0833   0.0000
 0.3750  -0.0000  -0.3750  -0.2500   0.7500
-0.5000   0.0833   0.5000  -0.0833  -0.0000
 0.1875  -0.0417   0.8125   0.1667  -0.6250
*Main> disp 4 $ d <> rd
5x5
 1.0000  -0.0000   0.0000   0.0000  0.0000
 0.0000   1.0000  -0.0000  -0.0000  0.0000
 0.0000   0.0000   1.0000   0.0000  0.0000
 0.0000  -0.0000  -0.0000   1.0000  0.0000
-0.0000  -0.0000   0.0000   0.0000  1.0000

*Main> e = matrix 3 [0,2,4,1,1,1,4,2,6]
*Main> e
(3><3)
 [ 0.0, 2.0, 4.0
 , 1.0, 1.0, 1.0
 , 4.0, 2.0, 6.0 ]
*Main> re = matinv e
*Main> disp 4 re
3x3
-0.3333   0.3333   0.1667
 0.1667   1.3333  -0.3333
 0.1667  -0.6667   0.1667
*Main> disp 4 $ e <> re
3x3
1.0000  0.0000  0.0000
0.0000  1.0000  0.0000
0.0000  0.0000  1.0000

disp n は行列 (Matrix Double) を表示する関数で、小数点数以下の桁数を引数 n で指定します。行列と逆行列の積を計算すると、浮動小数点数の計算による誤差がありますが、単位行列になることがわかります。

ところで、実際に連立一次方程式を解くとき、逆行列を使用することはほとんどありません。わざわざ逆行列を求めるよりも、ガウスの消去法を使ったほうがより速く解を求めることができるからです。同じ係数行列の連立方程式を何度も解く場合でも、逆行列を求めるのではなく、次に説明する LU 分解を使うことが多いようです。

●参考文献・URL

  1. 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. ガウスの消去法 - Wikipedia

●プログラムリスト

--
-- gauss.hs : ガウスの消去法
--
--            Copyright (C) 2021 Makoto Hiroi
--
import Prelude hiding ((<>))
import Numeric.LinearAlgebra

-- 前進消去
gauss1 :: [Vector R] -> Int -> [Vector R]
gauss1 [v] _ = [v]
gauss1 (v:vs) n =
  v : gauss1 (map (\w -> let temp = (w ! n) / (v ! n)
                         in w - (v * scalar temp)) vs)
             (n + 1)

-- 後退代入
gauss2 :: [Vector R] -> Int -> [R]
gauss2 [v] n = [(v ! (n + 1)) / (v ! n)]
gauss2 (v:vs) n =
  let m = n + 1
      x = gauss2 vs m
      y = subVector m (size v - 1 - m) v <.> vector x
      z = ((v ! (size v - 1)) - y)  / (v ! n)
  in z : x

-- ガウスの消去法
gauss :: Matrix R -> Vector R -> [R]
gauss xs ys =
  let zs = xs ||| asColumn ys
  in gauss2 (gauss1 (toRows zs) 0) 0

-- ガウス・ジョルダン法
gj1 :: [Vector R] -> Int -> Int -> [Vector R]
gj1 zs n m
  | n == m    = zs
  | otherwise = gj1 (gj2 zs 0) (n + 1) m
      where
        v = zs !! n
        w = v / scalar (v ! n)
        gj2 [] _ = []
        gj2 (x:xs) i =
          (if i == n then w else x - (w * scalar (x ! n))) : gj2 xs (i + 1)

gaussJordan :: Matrix R -> Vector R -> [R]
gaussJordan xs ys =
  let zs = xs ||| asColumn ys
  in map (\v -> v ! (size v - 1)) $ gj1 (toRows zs) 0 (rows xs)

--
-- ピボット選択
--

-- 最大値の行を選ぶ
findMaxRow :: [Vector R] -> Int -> Int
findMaxRow [] _ = 0
findMaxRow (z:zs) n =
  iter zs 1 (abs (z ! n)) 0
  where
    iter [] _ _ b = b
    iter (z:zs) i a b =
      let c = abs (z ! n)
      in if a < c
         then iter zs (i + 1) c i
         else iter zs (i + 1) a b

-- 行を交換する
swapRow :: [a] -> Int -> Int -> [a]
swapRow zs m n =
  iter zs 0
  where
    x = zs !! m
    y = zs !! n
    iter [] _ = []
    iter (z:zs) i
      | i == m    = y : iter zs (i + 1)
      | i == n    = x : iter zs (i + 1)
      | otherwise = z : iter zs (i + 1)

-- ピボット選択
selectPivot :: [Vector R] -> Int -> Int -> [Vector R]
selectPivot zs m n =
  let i = m + findMaxRow (drop m zs) n
  in if m == i
     then zs
     else swapRow zs m i

-- 前進消去 
gauss1' :: [Vector R] -> Int -> [Vector R]
gauss1' [v] _ = [v]
gauss1' (v:vs) n =
  v : gauss1' (selectPivot
               (map (\w -> let temp = (w ! n) / (v ! n)
                           in w - (v * scalar temp)) vs)
               0
               (n + 1))
              (n + 1)

-- ガウスの消去法
gauss' :: Matrix R -> Vector R -> [R]
gauss' xs ys =
  let zs = xs ||| asColumn ys
  in gauss2 (gauss1' (selectPivot (toRows zs) 0 0) 0) 0

-- ガウス・ジョルダン法
gj1' :: [Vector R] -> Int -> Int -> [Vector R]
gj1' zs n m
  | n == m    = zs
  | otherwise = gj1' (selectPivot (gj2 zs 0) (n + 1) (n + 1)) (n + 1) m
      where
        v = zs !! n
        w = v / scalar (v ! n)
        gj2 [] _ = []
        gj2 (x:xs) i =
          (if i == n then w else x - (w * scalar (x ! n))) : gj2 xs (i + 1)

gaussJordan' :: Matrix R -> Vector R -> [R]
gaussJordan' xs ys =
  let zs = xs ||| asColumn ys
  in map (\v -> v ! (size v - 1)) $ gj1' (selectPivot (toRows zs) 0 0) 0 (rows xs)

-- 逆行列
matinv :: Matrix R -> Matrix R
matinv xs =
  let
    n = rows xs
    zs = xs ||| ident n
  in
    dropColumns n $ fromRows $ gj1' (selectPivot (toRows zs) 0 0) 0 n

Copyright (C) 2021 Makoto Hiroi
All rights reserved.

[ PrevPage | Haskell | NextPage ]