Haskell (hmatrix) には行列 \(A\) の行列式 (determinant) を求める関数 det \(A\) が用意されていますが、\(A\) を LU 分解するときに計算することもできます。次に示す行列式と三角行列の性質を使います。以下では \(A\) の行列式を \(|A|\) と表しています。
ここでは数学的な証明は行わずに簡単な例を示すだけにとどめます。
(1)
[ a, b , * [ e, f, = [ ae + bg, af + bh,
c, d ] g, h ] ce + dg, cf + dh ]
右辺の行列式の乗算 (ad - bc) * (eh - fg) = adeh - adfg - bceh + bcfg
左辺の行列式 (ae + bg)(cf + dh) - (af + bh)(ce + dg)
= acef + adeh + bcfg + bdgh - acef - adfg - bceh - bdgh
^^^^ ++++ ^^^^ ++++
= adeh - adfg - bceh + bcfg
(2)
[ a, b, [ c, d,
c, d ] a, b ]
ad - bc cb - da = - (ad - bc)
(3) [ a, b, [ a, 0, 0, d ] c, d ] ad - b0 = ad ad - 0c = ad [ a, b, c, d, e, f, => aei + bfg + chd - ceg - bdi - ahf g, h, i ], 上三角行列は d, h, g が 0 なので行列式は aei になる 下三角行列は b, c, f が 0 なので行列式は aei になる
前回作成したプログラムでは行列 \(A\) を LU に分解したとき、\(L\) の対角成分が 1 になるので、\(U\) の対角成分の積が \(A\) の行列式になります。ピボット選択を行う場合、行列式の符号を記憶しておいて、行を交換したときは符号を反転すればいいでしょう。
プログラムは次のようになります。
リスト : 行列式の計算 (lu.hs を修正)
-- ピボット選択
selectPivot' :: [Vector R] -> Int -> [Int] -> ([Vector R], [Int], R)
selectPivot' zs n idx =
let i = findMaxRow zs n
in if i == 0
then (zs, idx, 1)
else (swapRow zs i 0, swapRow idx (n + i) n, -1)
-- 上三角行列
makeUp' :: [Vector R] -> Int -> [Int] -> R -> ([Vector R], [Int], R)
makeUp' [v] _ idx sgn = ([v], idx, sgn)
makeUp' (v:vs) n idx sgn =
let
a = map (\w -> let temp = (w ! n) / (v ! n)
in w - (v * scalar temp))
vs
(b, d, s1) = selectPivot' a (n + 1) idx
(xs, ys, s2) = makeUp' b (n + 1) d (sgn * s1)
in
(v : xs, ys, s2)
-- LU 分解
makeLU' :: Matrix R -> (Matrix R, Matrix R, [Int], R)
makeLU' xs =
let
a = toRows xs
(b, d, s) = selectPivot' a 0 [0..((rows xs) - 1)]
(ys, d1, s1) = makeUp' b 0 d s
zs = makeLow (transfer a d1) 0
in
(fromColumns $ map fromList zs, fromRows ys, d1, s1)
-- 行列式を求める
det' :: Matrix R -> R
det' xs =
let (_, us, _, sgn) = makeLU' xs
in sgn * (prodElements $ takeDiag us)
selectPivot' は行列式の符号をタプルに格納して返すように変更します。行を交換したら -1 を、交換しなければ 1 を返します。makeUp' は引数 sgn に行列式の符号を受け取って、符号の結果をタプルに格納して返します。makeLU' も行列式の符号をタプルに格納して返します。関数 det' は makeLU' の返り値をタプル (_, us, _sgn) で受け取ります。そして、takeDiag で上三角行列の対角成分を取り出し、prodElements で要素を掛け算して、最後に符号 sgn を掛け算します。対角成分に 0 の要素があるならば行列式は 0 になります。
それでは実際に試してみましょう。
ghci> a = matrix 2 [1,1,2,4] ghci> makeLU' a ((2><2) [ 1.0, 0.0 , 0.5, 1.0 ],(2><2) [ 2.0, 4.0 , 0.0, -1.0 ],[1,0],-1.0) ghci> det' a 2.0 ghci> det a 2.0 ghci> b = matrix 3 [1,1,1,2,4,6,2,0,4] ghci> makeLU' b ((3><3) [ 1.0, 0.0, 0.0 , 1.0, 1.0, 0.0 , 0.5, 0.25, 1.0 ],(3><3) [ 2.0, 4.0, 6.0 , 0.0, -4.0, -2.0 , 0.0, 0.0, -1.5 ],[1,2,0],1.0) ghci> det' b 12.0 ghci> det b 12.0 ghci> c = matrix 4 [1,1,1,1,-1,1,-1,1,8,4,2,1,-8,4,-2,1] ghci> makeLU' c ((4><4) [ 1.0, 0.0, 0.0, 0.0 , -1.0, 1.0, 0.0, 0.0 , 0.125, 6.25e-2, 1.0, 0.0 , -0.125, 0.1875, -1.0, 1.0 ],(4><4) [ 8.0, 4.0, 2.0, 1.0 , 0.0, 8.0, 0.0, 2.0 , 0.0, 0.0, 0.75, 0.75 , 0.0, 0.0, 0.0, 1.5 ],[2,3,0,1],1.0) ghci> det' c 72.0 ghci> det c 72.0 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> makeLU' d ((5><5) [ 1.0, 0.0, 0.0, 0.0, 0.0 , 1.0, 1.0, 0.0, 0.0, 0.0 , 0.3333333333333333, 0.4166666666666667, 1.0, 0.0, 0.0 ,8.333333333333333e-2,-4.1666666666666664e-2,0.625, 1.0, 0.0 ,8.333333333333333e-2, 0.125,0.625,-0.23076923076923078, 1.0 ], (5><5) [ 12.0, -6.0, 2.0, 0.0, 0.0 , 0.0, 12.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 1.3333333333333335, 1.0, 0.0 , 0.0, 0.0, -1.1102230246251565e-16, -1.625, 1.0 , 0.0, 0.0, -1.3664283380001927e-16, 0.0, 1.2307692307692308 ], [1,3,4,0,2],-1.0) ghci> det' d 384.0000000000001 ghci> det d 384.0000000000001 ghci> e = matrix 4 [2, 4, 2, 2, 4, 10, 3, 3, 2, 6, 1, 1, 3, 7, 1, 4] ghci> makeLU' e ((4><4) [ 1.0, 0.0, 0.0, 0.0 , 0.5, 1.0, 0.0, 0.0 , 0.75, 0.5, 1.0, 0.0 , 0.5, -1.0, -0.0, 1.0 ],(4><4) [ 4.0, 10.0, 3.0, 3.0 , 0.0, -1.0, 0.5, 0.5 , 0.0, 0.0, -1.5, 1.5 , 0.0, 0.0, 0.0, 0.0 ],[1,0,3,2],1.0) ghci> det' e 0.0 ghci> det e 0.0
最後は行列式が 0 になる場合です。この場合、解を一意的に定めることはできません。解は無数にあるか存在しないかのどちらかになります。
連立一次方程式の解の存在条件は行列の「階数 (rank)」を使って判定することができます。階数の定義はいくつか方法がありますが、ここでは行列の基本変形を使うことにしましょう。以下に示す操作を行の基本変形といいます。
同様に列の基本変形を定義することができます。これらの基本変形は、ガウスの消去法やガウス・ジョルダン法で用いる操作と同じです。一般に、任意の行列は基本変形を何回か適用すると「階段行列」に変形することが可能です。階段行列は左下半分に 0 が階段状に並んだ行列のことです。
[ a, b, c, d, [ a, b, c, d, [ a, b, c, d,
0, e, f, g, 0, e, f, g, 0, e, f, g,
0, 0, h, i ] 0, 0, 0, h ] 0, 0, 0, 0 ]
(1) rank = 3 (2) rank = 3 (3) rank = 2
図 : 階段行列,
階段行列に変形したあと、零ベクトル (要素がすべて 0 のベクトル) を除いた行ベクトルの本数が階数になります。上図でいうと、(1) と (2) の行列には零ベクトルはないので階数は 3 になりますが、(3) の行列は零ベクトルが一つあるので階数は 2 になります。
ここで、上図の行列が拡大係数行列を表しているとしましょう。(1) の場合、係数行列と拡大係数行列の階数はどちらも 3 になり、拡大係数行列の行数と一致します。この場合、連立一次方程式の解は一意的に定まります。
(2) の場合、3 行 3 列の係数行列としてみると、その階数は 2 になり、拡大係数行列の階数よりも少なくなります。この場合は解がありません。係数がすべて 0 なのに定数が残っているので、連立方程式として矛盾しているわけです。
(3) の場合、係数行列と拡大係数行列の階数はともに 2 になりますが、拡大係数行列の行数とは一致しません。このとき、零ベクトルに対応する変数に適当な値を代入すると、残りの変数を決定することができます。つまり、連立一次方程式の解は存在するが、一意的には定めることができないわけです。まとめると次のようになります。
簡単な例を示しましょう。
(1) 2a + 4b + 2c + 2d = 8 4a + 10b + 3c + 3d = 17 2a + 6b + c + d = 9 3a + 7b + c + 4d = 11
行の基本変形 [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 1. -0.5 -0.5 0.5 0. -0.5 -1.25 1.75 -1.75 ] [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 0. 0. 0. 0. 0. 0. -1.5 1.5 -1.5 ] [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 0. -1.5 1.5 -1.5 0. 0. 0. 0. 0. ] 係数行列と拡大係数行列の階数は 3 d = k とおくと c = (-1.5 - 1.5d) / -1.5 = k + 1 b = (-0.5 - 0.5c - 0.5d) / -1 = k + 1 a = (17 - 10b - 3c - 3d) / 4 = -4k + 1 k は任意の定数なので解は無数にある
(2) 2a + 4b + 2c + 2d = 8 4a + 10b + 3c + 3d = 17 2a + 6b + c + d = 10 3a + 7b + c + 4d = 11 行の基本変形 [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 1. -0.5 -0.5 1.5 0. -0.5 -1.25 1.75 -1.75 ] [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 0. 0. 0. 1. 0. 0. -1.5 1.5 -1.5 ] [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 0. -1.5 1.5 -1.5 0. 0. 0. 0. 1. ] 係数行列の階数は 3 で拡大係数行列の階数は 4 なので解はない
(3) 2a + 4b + 2c + 2d = 8 4a + 10b + 3c + 3d = 17 2a + 6b + c + 2d = 9 3a + 7b + c + 4d = 11
行の基本変形 [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 1. -0.5 0.5 0.5 0. -0.5 -1.25 1.75 -1.75 ] [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 0. 0. 1. 0. 0. 0. -1.5 1.5 -1.5 ] [ 4. 10. 3. 3. 17. 0. -1. 0.5 0.5 -0.5 0. 0. -1.5 1.5 -1.5 0. 0. 0. 1. 0. ] 係数行列と拡大係数行列の階数は 4 d = 0 / 1 = 0 c = (-1.5 - 1.5*0) / -1.5 = 1 b = (-0.5 - 0.5*1 - 0.5*0) / -1 = 1 a = (17 - 10*1 - 3*1 - 3*0) / 4 = 1 解は一意的に定まる
拙作のページで説明した「ガウスの消去法」は、拡大係数行列を階段行列に変形することで解を求めています。このような方法を「直接法 (Direct Method)」といいます。これに対し、適当な初期解から始めて、繰り返し計算することで真の解に収束させていく方法が考えられます。これを「反復法 (Iterative Method)」といいます。
直接法は厳密解を求めることができますが、係数行列が大きくなると解くのに時間がかかるようになります。厳密解を求めるのが難しい場合、反復法を使うと現実的な時間で近似解を求めることが可能です。今回は基本的な反復法である「ヤコビ法」と「ガウス・ザイデル法」について簡単に説明します。
連立方程式 \(Ax = b\) の係数行列 \(A\) を対角行列 \(D\) とそれ以外の要素を持つ行列 \(A'\) に分解します。すると、方程式は次のように変形することができます。
ここで、最後の式を漸化式 \(x_{i+1} = D^{-1}(b - A'x_i)\) と考えて、反復処理により解を求めるのが「ヤコビ法 (Jacobi Method)」です。たとえば、三元連立方程式の場合、漸化式は下図のようになります。
a1 * x + a2 * y + a3 * z = d1
b1 * x + b2 * y + b3 * z = d2
c1 * x + c2 * y + c3 * z = d3
A' = [ 0, a2, a3, D = [ a1, 0, 0, D-1 = [ 1/a1, 0, 0,
b1, 0, b3, 0, b2, 0, 0, 1/b2, 0,
c1, c2, 0 ] 0, 0, c3 ] 0, 0, 1/c3 ]
xi+1 = (d1 - a2 * yi - a3 * zi) / a1
yi+1 = (d2 - b1 * xi - b3 * zi) / b2
zi+1 = (d3 - c1 * xi - c2 * yi) / c3
ヤコビ法で解が収束した場合、それが元の連立方程式の解となります。なお、行列 \(A\) によっては解が収束しない (発散する) 場合があります。収束の十分条件ですが、『係数行列 \(A\) が対角優位な行列である場合に収束する』ことが知られています。行列 \(A\) の対角成分 \(a_{ii}\) の絶対値が、他の行の成分 \(a_{ij}\) の絶対値の合計よりも大きいことを「対角優位」といいます。
上図の三元連立方程式でいえば、|a1| > |b1| + |c1|, |b2| > |a2| + |c2|, |c3| > |a3| + |b3| を満たすとき、ヤコビ法の解は収束することになります。
収束の判定ですが、一般的には \(x_{i+1}\) と \(x_i\) の差分が許容誤差 \(\varepsilon\) に収まったときに収束と判定します。具体的な方法を以下に示します。
1 は差分の合計値、2 は差分を解で割った値の合計値、3 は差分を解で割った値の最大値が ε 以下になったとき収束と判定します。今回は 2 の方法を使うことにします。
それではプログラムを作りましょう。Haskell (hmatrix) を使うとヤコビ法は簡単にプログラムすることができます。
リスト : ヤコビ法
jacobi :: Matrix R -> Matrix R -> Maybe (Matrix R, Int)
jacobi a b =
iter 0 (col $ replicate (rows a) 0.0)
where
d = asColumn $ takeDiag a
iter i x
| i >= 512 = Nothing
| otherwise = let y = (b - (a <> x - d * x)) / d
in if sumElements (cmap abs ((y - x) / y)) < 1e-6
then Just (y, i + 1)
else iter (i + 1) y
引数 a が係数行列、b が右辺値を格納した配列です。最初に a の対角成分を takeDiag で取り出し、列ベクトルに変換して変数 d にセットします。実際の処理は局所関数 iter で行います。引数 x には解を格納する列ベクトルをセットします。x は 0 で初期化します。あとは再帰呼び出しで漸化式 y = D-1(b - A'x) を繰り返し計算します。実際には、a <> x を計算して、そこから d * x を引き算しています。hmatrix の場合、演算子 / は要素同士の割り算になるので、対角行列 d の逆行列を求める必要はありません。
収束条件のチェックも簡単です。sumElements (cmap abs ((y - x) / y)) で合計値を求め、それが 1e-6 よりも小さくなったならば解は収束したので Just (y, i) を返します。繰り返しが 512 回以上になった場合、解は収束していないので Nothing を返します。
それでは実行してみましょう。
9x + y + 2z = 9 x + 9y + z = 18 2x + y + 9z = -5
ghci> jacobi (matrix 3 [9,1,2,1,9,1,2,1,9]) (col [9,18,-5]) Just ((3><1) [ 0.9999999667086195 , 1.9999999761533003 , -1.0000000318590117 ],14)
厳密解は [1, -2, -1] ですが、ヤコビ法だと 14 回で収束しています。次のように Debug.Trace の関数を使うと、収束の様子を表示することができます。
let y = traceShowId $ (b - (a <> x - d * x)) / d
興味のある方は試してみてください。ガウス・ザイデル法を使うと、これよりも少ない回数で解を求めることができます。
ヤコビ法で漸化式を計算するとき、右辺式は xi, yi zi の値を使いますが、1 行目から順番に計算していくと、2 行目の yi+1 を計算するときには xi+1 の値がすでに求まっています。3 行目の zi+1 を計算するときは、xi+1 と yi+1 の値が求まっています。これらの値を使って漸化式を計算する方法を「ガウス・ザイデル法 (Gauss-Siedel Method)」といいます。
たとえば、三元連立方程式の場合、ガウス・ザイデル法の漸化式は下図のようになります。
a1 * x + a2 * y + a3 * z = d1
b1 * x + b2 * y + b3 * z = d2
c1 * x + c2 * y + c3 * z = d3
A' = [ 0, a2, a3, D = [[a1, 0, 0, D-1 = [ 1/a1, 0, 0,
b1, 0, b3, [ 0, b2, 0, 0, 1/b2, 0,
c1, c2, 0 ] [ 0, 0, c3 ] 0, 0, 1/c3 ]
xi+1 = (d1 - a2 * yi - a3 * zi) / a1
yi+1 = (d2 - b1 * xi+1 - b3 * zi) / b2
zi+1 = (d3 - c1 * xi+1 - c2 * yi+1) / c3
ガウス・ザイデル法の解の収束条件や収束の判定方法はヤコビ法と同じです。
それではプログラムを作りましょう。次のリストを見てください。
リスト : ガウス・ザイデル法
gaussSeidel :: Matrix R -> Matrix R -> Maybe (Matrix R, Int)
gaussSeidel a b =
iter 0 (col $ replicate (rows a) 0.0)
where
d = takeDiag a
zs = toRows a
iter i x
| i >= 512 = Nothing
| otherwise = let y = iter2 zs 0 (toList $ flatten x) []
in if sumElements (cmap abs ((y - x) / y)) < 1e-6
then Just (y, i + 1)
else iter (i + 1) y
iter2 [] _ _ ys = col ys
iter2 (z:zs) i (_:xs) ys =
let y = ((b ! i ! 0) - (z <.> vector (ys ++ [0] ++ xs))) / (d ! i)
in iter2 zs (i + 1) xs (ys ++ [y])
ポイントは局所関数 iter2 で変数の値を逐次的に計算していくところです。第 1 引数には係数行列を toRows でリストに変換したものを渡します。第 2 引数が処理する変数の位置、第 3 引数が一つ前の変数の値を格納したリスト、第 4 引数が求めた変数の値を格納したリストです。
vector (ys ++ [0] ++ xs) で求めた変数の値と、ひとつ前の変数の値を結合してベクトルに変換し、それと z との内積を計算します。その値を b ! i ! 0 から引き算して d ! i で割り算すれば、i 番目の変数の値を求めることができます。iter2 を再帰呼び出しするときは、第 3 引数の先頭要素を取り除き、求めた変数の値 y を ys の末尾に追加します。
それでは実際に試してみましょう。
9x + y + 2z = 9 x + 9y + z = 18 2x + y + 9z = -5
ghci> gaussSeidel (matrix 3 [9,1,2,1,9,1,2,1,9]) (col [9,18,-5]) Just ((3><1) [ 0.9999999187463209 , 2.000000011961912 , -0.9999999832727282 ],6)
ヤコビ法だと収束するまで 14 回かかっていたのが、ガウス・ザイデル法を使うと 6 回ですみました。ただし、逐次的に処理を行うため、プログラムがちょっと複雑になってしまいました。hmatrix を無理やり使うよりも、mutable な配列を使ったほうがよかったのかもしれません。
それでは簡単なテストとして、ガウスの消去法、ヤコビ法、ガウス・ザイデル法の実行時間を比較してみましょう。次のリストを見てください。
リスト : 簡単なテスト -- テストデータの生成 makeData :: Int -> IO (Matrix R, Matrix R, Matrix R) makeData n = do a <- rand n n x <- rand n 1 let a1 = a + (diag $ vector $ map (\v -> 2 * sumElements v) (toRows a)) return (a1, x, a1 <> x) -- 答え合わせ checkAnswer :: (Matrix R, Matrix R, Matrix R) -> Maybe (Matrix R, Int) -> Maybe (R, Int) checkAnswer _ Nothing = Nothing checkAnswer (a, x, b) (Just (y, n)) = Just (sumElements (cmap abs (y - x)), n)
乱数でテストデータを生成します。これを関数 makeData で行っています。連立方程式 Ax = b の係数行列 A と解 x を乱数で生成して、右辺値の b を計算します。係数行列は対角優位になるよう修正しています。checkAnswer は makeData と解法プログラムの結果を受け取り、x との誤差を計算します。
実行結果は次のようになりました。
ghci> :set +s ghci> (a, x, b) <- makeData 500 (0.00 secs, 35,280 bytes) ghci> checkAnswer (a, x, b) (jacobi a b) Just (2.6354544861652458e-8,33) (0.06 secs, 25,852,096 bytes) ghci> checkAnswer (a, x, b) (gaussSeidel a b) Just (9.66711508260977e-10,12) (0.15 secs, 236,512,152 bytes) ghci> checkAnswer (a, x, b) (Just ((gauss' a b), 0)) Just (3.150071349600214e-13,0) (0.86 secs, 1,792,777,824 bytes) ghci> (a, x, b) <- makeData 750 (0.00 secs, 35,208 bytes) ghci> checkAnswer (a, x, b) (jacobi a b) Just (1.989642568560844e-8,34) (0.06 secs, 55,191,544 bytes) ghci> checkAnswer (a, x, b) (gaussSeidel a b) Just (1.3893328199053936e-9,12) (0.43 secs, 554,765,352 bytes) ghci> checkAnswer (a, x, b) (Just ((gauss' a b), 0)) Just (5.717837227997569e-13,0) (3.13 secs, 5,720,292,096 bytes) ghci> (a, x, b) <- makeData 1000 (0.00 secs, 35,224 bytes) ghci> checkAnswer (a, x, b) (jacobi a b) Just (2.7103529010334346e-8,34) (0.11 secs, 95,416,912 bytes) ghci> checkAnswer (a, x, b) (gaussSeidel a b) Just (2.004639381033256e-9,12) (0.74 secs, 909,521,752 bytes) ghci> checkAnswer (a, x, b) (Just ((gauss' a b), 0)) Just (8.207556813008554e-13,0) (5.63 secs, 13,166,904,016 bytes) 実行環境 : Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
変数の個数が増えるにしたがい、ガウスの消去法の実行時間は大幅に増加します。反復法でも遅くはなりますが、ガウスの消去法よりはずっと高速です。ヤコビ法とガウス・ザイデル法を比較した場合、ヤコビ法のほうが速くなりました。反復回数はガウス・ザイデル法のほうが少ないので、逐次的に計算する処理で時間がかかっていると思われます。hmatrix を使うのであれば、なるべく行列 (ベクトル) で計算したほうがよさそうです。
もちろん、コンパイルすれば実行速度はもっと速くなると思います。興味のある方はいろいろ試してみてください。
--
-- jacobi.hs : 反復法
--
-- Copyright (C) 2021 Makoto Hiroi
--
import Prelude hiding ((<>))
import Numeric.LinearAlgebra
import Debug.Trace
--
-- ガウスの消去法
--
-- ピボット選択
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)
-- 後退代入
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 -> Matrix R -> Matrix R
gauss' xs ys =
let zs = xs ||| ys
in col $ gauss2 (gauss1' (selectPivot (toRows zs) 0 0) 0) 0
--
-- 反復法
--
-- ヤコビ法
jacobi :: Matrix R -> Matrix R -> Maybe (Matrix R, Int)
jacobi a b =
iter 0 (col $ replicate (rows a) 0.0)
where
d = asColumn $ takeDiag a
iter i x
| i >= 512 = Nothing
| otherwise = let y = (b - (a <> x - d * x)) / d
in if sumElements (cmap abs ((y - x) / y)) < 1e-6
then Just (y, i + 1)
else iter (i + 1) y
-- ガウス・ザイデル法
gaussSeidel :: Matrix R -> Matrix R -> Maybe (Matrix R, Int)
gaussSeidel a b =
iter 0 (col $ replicate (rows a) 0.0)
where
d = takeDiag a
zs = toRows a
iter i x
| i >= 512 = Nothing
| otherwise = let y = iter2 zs 0 (toList $ flatten x) []
in if sumElements (cmap abs ((y - x) / y)) < 1e-6
then Just (y, i + 1)
else iter (i + 1) y
iter2 [] _ _ ys = col ys
iter2 (z:zs) i (_:xs) ys =
let y = ((b ! i ! 0) - (z <.> vector (ys ++ [0] ++ xs))) / (d ! i)
in iter2 zs (i + 1) xs (ys ++ [y])
-- テストデータの生成
makeData :: Int -> IO (Matrix R, Matrix R, Matrix R)
makeData n = do
a <- rand n n
x <- rand n 1
let a1 = a + (diag $ vector $ map (\v -> 2 * sumElements v) (toRows a))
return (a1, x, a1 <> x)
-- 答え合わせ
checkAnswer :: (Matrix R, Matrix R, Matrix R) -> Maybe (Matrix R, Int)
-> Maybe (R, Int)
checkAnswer _ Nothing = Nothing
checkAnswer (a, x, b) (Just (y, n)) = Just (sumElements (cmap abs (y - x)), n)