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

●特殊な行列とその性質

転置行列 (transposed matrix)

行列 A の行と列を入れ替えたもので、\(A^{\mathrm{T}}\) と表記する
Haskell (hmatrix) の場合、関数 tr A で転置行列を求めることができる
転置行列との積 \(AA^{\mathrm{T}}\) を求める関数 mTm A もある

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

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

ghci> c = matrix 2 [1..6]
ghci> c
(3><2)
 [ 1.0, 2.0
 , 3.0, 4.0
 , 5.0, 6.0 ]
ghci> tr c
(2><3)
 [ 1.0, 3.0, 5.0
 , 2.0, 4.0, 6.0 ]
ghci> mTm c
Herm (2><2)
 [ 35.0, 44.0
 , 44.0, 56.0 ]
ghci> (tr c) <> c
(2><2)
 [ 35.0, 44.0
 , 44.0, 56.0 ]

転置行列の基本的な性質 (A, B は正方行列)

\(trace \ A^{\mathrm{T}} = trace \ A\)
\(|A^{\mathrm{T}}| = |A|\), (\(|A|\) は \(A\) の行列式)
\((AB)^{\mathrm{T}} = B^{\mathrm{T}}A^{\mathrm{T}}\)
\((A^{\mathrm{T}})^{-1} = (A^{-1})^{\mathrm{T}}\)

ghci> a
(2><2)
 [ 1.0, 2.0
 , 3.0, 4.0 ]
ghci> trace a
5.0
ghci> trace $ tr a
5.0
ghci> det a
-2.0
ghci> det $ tr a
-2.0

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

ghci> inv (tr a)
(2><2)
 [ -2.0,  1.5
 ,  1.0, -0.5 ]
ghci> tr (inv a)
(2><2)
 [ -1.9999999999999998,   1.4999999999999998
 ,                 1.0, -0.49999999999999994 ]

対称行列と対角行列
- 正方行列 \(A\) とその転置行列 \(A^{\mathrm{T}}\) が等しい行列のことを「対称行列 (symmetric matrix)」という
- その中で対角成分 \(a_{ii}\) 以外の要素が 0 の行列のことを「対角行列 (diagonal matrix)」という
- 単位行列は対角行列の一種

ghci> a = matrix 2 [1,3,3,2]
ghci> a
(2><2)
 [ 1.0, 3.0
 , 3.0, 2.0 ]
ghci> tr a == a
True
ghci> b = matrix 3 [1, 4, 5, 4, 2, 6, 5, 6, 3]
ghci> b
(3><3)
 [ 1.0, 4.0, 5.0
 , 4.0, 2.0, 6.0
 , 5.0, 6.0, 3.0 ]
ghci> tr b == b
True
ghci> c = matrix 3 [1,0,0,0,2,0,0,0,3]
ghci> c
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 2.0, 0.0
 , 0.0, 0.0, 3.0 ]
ghci> tr c == c
True

操作関数
sym A, 行列 \(A\) から対称行列を生成する
- 計算式は (A + tr A) / 2
- データ型は対称行列を表すものになる
diag v, diagl xs, 対角行列の生成

diag :: (Num a, Element a) => Vector a -> Matrix a
diagl :: [Double] -> Matrix Double

ベクトル v, リスト xs の要素が対角成分になる

takeDiag A, 行列 A の対角線上の要素を取得する
diag, diagl は対角線以外の要素が 0 になるが、値を指定したい場合は関数 diagRect を使う

diagRect :: Foreign.Storable.Storable t => t -> Vector t -> Int -> Int -> Matrix t

第 1 引数が指定する値、第 2 引数が対角線上の要素、第 3, 4 引数が行列の大きさ

ghci> a = matrix 3 [1..9]
ghci> a
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
ghci> sym a
Herm (3><3)
 [ 1.0, 3.0, 5.0
 , 3.0, 5.0, 7.0
 , 5.0, 7.0, 9.0 ]

ghci> takeDiag a
[1.0,5.0,9.0]
ghci> 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 ]
ghci> diagl [1,5,9]
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 5.0, 0.0
 , 0.0, 0.0, 9.0 ]

ghci> 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 ]
ghci> 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 ]
ghci> 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 ]

対角行列の基本的な性質 (\(A\) は対角行列)
1. \(|A|\) は対角成分 \(a_{ii}\) を乗算した値になる
2. \(A^n\) は対角成分 \(a_{ii}\) を n 乗したものになる
3. \(A^{-1}\) は対角成分 \(a_{ii}\) を逆数にしたものになる
4. 対角行列の固有値は対角成分 \(a_{ii}\) になる (固有値は後で説明する)

ghci> a = diagl [1,2,3]
ghci> a
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 2.0, 0.0
 , 0.0, 0.0, 3.0 ]
ghci> det a
6.0
ghci> :t prodElements
prodElements :: Container c e => c e -> e
ghci> prodElements (takeDiag a :: Vector R)
6.0

ghci> a <> a
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 4.0, 0.0
 , 0.0, 0.0, 9.0 ]
ghci> a <> a <> a
(3><3)
 [ 1.0, 0.0,  0.0
 , 0.0, 8.0,  0.0
 , 0.0, 0.0, 27.0 ]

ghci> inv a
(3><3)
 [ 1.0, 0.0,                0.0
 , 0.0, 0.5,                0.0
 , 0.0, 0.0, 0.3333333333333333 ]
ghci> 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 ]
ghci> eigenvalues a
[1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0]

eigenvalues A は行列 \(A\) の固有値を求める関数

eigenvalues :: Field t => Matrix t -> Vector (Complex Double)

Complex Double は実部と虚部が Double 型の複素数
複素数は実部 :+ 虚部で表す

三角行列

上三角行列 (Upper triangular matrix) は対角成分よりも下側の成分が 0 である行列のこと
下三角行列 (lower triangular matrix) は対角成分よりも上側の成分が 0 である行列のこと

\(\begin{array}{llll} \begin{pmatrix} a & b & c \\ 0 & d & e \\ 0 & 0 & f \end{pmatrix} & \begin{pmatrix} a & b & c & d \\ 0 & e & f & g \\ 0 & 0 & h & i \\ 0 & 0 & 0 & j \end{pmatrix} & \begin{pmatrix} 0 & 0 & a \\ 0 & b & c \\ d & e & f \end{pmatrix} & \begin{pmatrix} 0 & 0 & 0 & a \\ 0 & 0 & b & c \\ 0 & d & e & f \\ g & h & i & j \end{pmatrix} \end{array}\)

三角行列の主な性質 (\(A\) は三角行列, \(L\) は下三角行列, \(U\) は上三角行列)
1. \(L^{\mathrm{T}}\) は \(U\) に, \(U^{\mathrm{T}}\) は \(L\) になる
2. \(U \times U\) は \(U\) に、\(L \times L\) は \(L\) になる
3. \(|A|\) は対角成分 \(a_{ii}\) を乗算した値になる
4. \(U^{-1}\) は \(U\) に、\(L^{-1}\) は \(L\) になる
5. \(A\) の固有値は対角成分 \(a_{ii}\) になる (固有値は後で説明する)

ghci> a = matrix 3 [1,2,3,4,5,6,7,8,10]
ghci> a
(3><3)
 [ 1.0, 2.0,  3.0
 , 4.0, 5.0,  6.0
 , 7.0, 8.0, 10.0 ]
ghci> la = matrix 3 [1,0,0,4,5,0,7,8,10]
ghci> la
(3><3)
 [ 1.0, 0.0,  0.0
 , 4.0, 5.0,  0.0
 , 7.0, 8.0, 10.0 ]
ghci> ua = matrix 3 [1,2,3,0,5,6,0,0,10]
ghci> ua
(3><3)
 [ 1.0, 2.0,  3.0
 , 0.0, 5.0,  6.0
 , 0.0, 0.0, 10.0 ]

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

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

ghci> det la
49.99999999999999
ghci> det ua
50.0

ghci> 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 ]
ghci> inv ua
(3><3)
 [ 1.0, -0.4,              -6.0e-2
 , 0.0,  0.2, -0.12000000000000002
 , 0.0,  0.0,                  0.1 ]

ghci> eigenvalues la
[10.0 :+ 0.0,5.0 :+ 0.0,1.0 :+ 0.0]
ghci> eigenvalues ua
[1.0 :+ 0.0,5.0 :+ 0.0,10.0 :+ 0.0]

直交行列 (orthogonal matrix)
- \(A^{\mathrm{T}}\) と \(A^{-1}\) が等しい行列のこと
- 直交行列の簡単な例として、回転行列や置換行列などがある

直交行列の基本的な性質
1. 直交行列の行列式は \(\pm1\)
2. 直交行列の逆行列も直交行列
3. 対称行列は直交行列で対角化できる (あとで説明する)

回転行列 \(R\) は直交行列
\( R = \begin{pmatrix} \cos x & -\sin x \\ \sin x & \cos x \end{pmatrix} \)

転置行列
\( R^{\mathrm{T}} = \begin{pmatrix} \cos x & \sin x \\ -\sin x & \cos x \end{pmatrix} \)

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

逆行列
\( R^{-1} = \begin{pmatrix} \cos x & \sin x \\ -\sin x & \cos x \end{pmatrix} \)

逆行列の転置 (\(R\) に戻るので直交行列)
\( (R^{-1})^{\mathrm{T}} = \begin{pmatrix} \cos x & -\sin x \\ \sin x & \cos x \end{pmatrix} \)

ghci> a = matrix 2 [0,1,1,0]
ghci> a
(2><2)
 [ 0.0, 1.0
 , 1.0, 0.0 ]
ghci> tr a
(2><2)
 [ 0.0, 1.0
 , 1.0, 0.0 ]
ghci> det a
-1.0
ghci> inv a
(2><2)
 [ 0.0, 1.0
 , 1.0, 0.0 ]

ghci> b = matrix 3 [0, 1, 0, 1, 0, 0, 0, 0, 1]
ghci> b
(3><3)
 [ 0.0, 1.0, 0.0
 , 1.0, 0.0, 0.0
 , 0.0, 0.0, 1.0 ]
ghci> tr b
(3><3)
 [ 0.0, 1.0, 0.0
 , 1.0, 0.0, 0.0
 , 0.0, 0.0, 1.0 ]
ghci> inv b
(3><3)
 [ 0.0, 1.0, 0.0
 , 1.0, 0.0, 0.0
 , 0.0, 0.0, 1.0 ]

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

各行各列にちょうど一つだけ 1 の要素を持ち、それ以外は全て 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 にするので、w_n / v_n を v に乗算し、それを w から減算すれば、w の n 番目の係数を 0 にすることができます。vs に格納されたベクトルの n 番目の要素を 0 にしたら、gauss1 を再帰呼び出しして、n + 1 番目の係数を 0 にします。

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

 x +  y = 100
2x + 4y = 272

ghci> 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

ghci> 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

ghci> 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

ghci> 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

ghci> 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

ghci> 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

ghci> 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

ghci> 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

ghci> 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

ghci> 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

ghci> 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)

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 番目であることに注意してください。

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

ghci> gauss (matrix 3 [0,2,4,1,1,1,4,2,6]) (vector [14,10,38])
[NaN,NaN,NaN]
ghci> gauss' (matrix 3 [0,2,4,1,1,1,4,2,6]) (vector [14,10,38])
[5.0,3.0,2.0]
ghci> 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]
ghci> gaussJordan (matrix 3 [0,2,4,1,1,1,4,2,6]) (vector [14,10,38])
[NaN,NaN,NaN]
ghci> gaussJordan' (matrix 3 [0,2,4,1,1,1,4,2,6]) (vector [14,10,38])
[5.0,3.0,2.0]
ghci> 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 で作成することができます。あとはガウス・ジョルダン法で解を求め、最後に拡大係数行列の右半分を返します。それでは実行してみましょう。

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

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

ghci> 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 ]
ghci> rc = matinv c
ghci> 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
ghci> 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

ghci> 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]
ghci> 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 ]
ghci> rd = matinv d
ghci> 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
ghci> 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

ghci> e = matrix 3 [0,2,4,1,1,1,4,2,6]
ghci> e
(3><3)
 [ 0.0, 2.0, 4.0
 , 1.0, 1.0, 1.0
 , 4.0, 2.0, 6.0 ]
ghci> re = matinv e
ghci> disp 4 re
3x3
-0.3333   0.3333   0.1667
 0.1667   1.3333  -0.3333
 0.1667  -0.6667   0.1667
ghci> 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

奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
ガウスの消去法 - 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