旧タイプなので、連立一次方程式を解くときに逆行列を使って解けばOK…と頭から思い込んでいたのですが、Gaussの消去法を使うと逆行列は出てこなくて「？」となるわけです。上三角化を使うのだから、LU分解でいいんですよね。
LU分解で行列式と逆行列の計算 その３: メモブログ
http://sssiii.seesaa.net/article/379560261.html
のところから、
連立１次方程式 I
http://nalab.mind.meiji.ac.jp/~mk/labo/text/linear-eq-1.pdf
のPDFを読んで、なるほど、LU分解で十分ですね。ってことで自前で実装してみます。
4 LU分解
http://akita-nct.jp/yamamoto/lecture/2004/5E/linear_equations/text/html/node4.html#eq:LU____________
を参考にして式を流用します。最初、手順がうまくわからなくて、自分で 4×4 の行列を作って手計算をした後に、数式に直して、元の式と一致することを確認。なるほど、行単位じゃなくて「列単位」で計算するので、計算する順番は、β11,α21,α31,α41 で計算した後に、β12,β22,α32,α42 と計算していくんですね。もう一度読み直したらそう書いてあるし…
出来上がったコードがこんな感じ。
[code lang=”csharp”]
/// LU分解
let LU ( A : matrix ) =
let n = A.NumCols
let L = Matrix.identity n
let U = Matrix.zero n n
for j in 0..n-1 do
// Uを計算
for i in 0..j do
U.[i,j] <- A.[i,j]
for k in 0..i-1 do
U.[i,j] <- U.[i,j] - L.[i,k]*U.[k,j]
// Lを計算
for i in j+1..n-1 do
L.[i,j] <- A.[i,j]
for k in 0..j-1 do
L.[i,j] <- L.[i,j] - L.[i,k]*U.[k,j]
L.[i,j] <- L.[i,j] / U.[j,j]
// 結果を返す
(L,U)
[/code]
次のように行列を設定しておくと
[code lang=”csharp”]
let A = matrix [[1.0;2.0;1.0];
[2.0;1.0;0.0];
[1.0;1.0;2.0]]
LU A
[/code]
こんな感じで、L と U にして返してくれます。
[code lang=”csharp”]
val it : matrix * matrix =
(matrix [[1.0; 0.0; 0.0]
[2.0; 1.0; 0.0]
[1.0; 0.3333333333; 1.0]], matrix [[1.0; 2.0; 1.0]
[0.0; -3.0; -2.0]
[0.0; 0.0; 1.666666667]])
[/code]
先のページにも書いてあるのですが、Ujj で割ってしまうので、ピボット選択は必須…なのか？（ちょっとよくわからず）。浮動小数点だったら大丈夫な気もするのですが、これはあとで確認してみます。まあ、割り算のところにチェックを入れるのは数値計算の基本なので、なるべく割り算で誤差が出ないようにする方向でよいかと。
ちなみに F# MathProvider を使うと
[code lang=”csharp”]
let (p,L,U) = MathProvider.LinearAlgebra.lu A
[/code]
として、次なる結果を得られます。
[code lang=”csharp”]
val p : (int -> int)
val U : matrix = matrix [[2.0; 1.0; 0.0]
[0.0; 1.5; 1.0]
[0.0; 0.0; 1.666666667]]
val L : matrix = matrix [[1.0; 0.0; 0.0]
[0.5; 1.0; 0.0]
[0.5; 0.3333333333; 1.0]]
[/code]
MathProvider の LU 分解のコードはもっとシンプルで、こんな感じになっています。うまく交互に計算すると L のほうはほとんど計算せずに済みみたいです。行列数が多くなると swapR のパフォーマンスが問題になりますが、有限要素法以外だったら大丈夫かと。
[code lang=”csharp”]
for i=0 to nA-2 do
let mutable maxi = i // Find the biggest pivot element.
for k=i+1 to nA-1 do
if abs(U.[maxi,i]) < abs(U.[k,i]) then maxi <- k
done
//let maxi = maxIndex (i+1) (nA-1) (fun k -> abs(U.[k,i]))
if maxi <> i then
swapR U i maxi // Swap rows if needed.
swapR L i maxi // REVIEW can be made more performant.
let t = P.[i]
P.[i] <- P.[maxi]
P.[maxi] <- t
for j=i+1 to nA-1 do
L.[j,i] <- U.[j,i] / U.[i,i]
for k=i+1 to nA-1 do
U.[j,k] <- U.[j,k] - L.[j,i] * U.[i,k]
done
U.[j,i] <- 0.0
done
done
[/code]
sum のところは、もっと F# らしく書き換えられないかな、と思案中です。コードが短くなるしΣを使ったほうが、数式を同じになるので。