[F#] LU分解を作る

投稿日: 2014/01/09 作成者: masuda

旧タイプなので、連立一次方程式を解くときに逆行列を使って解けば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 と計算していくんですね。もう一度読み直したらそう書いてあるし…

出来上がったコードがこんな感じ。 

/// 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.&#91;i,j&#93;
            for k in 0..i-1 do
                U.&#91;i,j&#93; <- U.&#91;i,j&#93; - L.&#91;i,k&#93;*U.&#91;k,j&#93;
        // Lを計算
        for i in j+1..n-1 do
            L.&#91;i,j&#93; <- A.&#91;i,j&#93;
            for k in 0..j-1 do
                L.&#91;i,j&#93; <- L.&#91;i,j&#93; - L.&#91;i,k&#93;*U.&#91;k,j&#93;
            L.&#91;i,j&#93; <- L.&#91;i,j&#93; / U.&#91;j,j&#93;
    // 結果を返す
    (L,U)
&#91;/code&#93;
<p>
次のように行列を設定しておくと
</p>
[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

こんな感じで、L と U にして返してくれます。 

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

先のページにも書いてあるのですが、Ujj で割ってしまうので、ピボット選択は必須…なのか？（ちょっとよくわからず）。浮動小数点だったら大丈夫な気もするのですが、これはあとで確認してみます。まあ、割り算のところにチェックを入れるのは数値計算の基本なので、なるべく割り算で誤差が出ないようにする方向でよいかと。

ちなみに F# MathProvider を使うと 

let (p,L,U) = MathProvider.LinearAlgebra.lu A

として、次なる結果を得られます。 

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

MathProvider の LU 分解のコードはもっとシンプルで、こんな感じになっています。うまく交互に計算すると L のほうはほとんど計算せずに済みみたいです。行列数が多くなると swapR のパフォーマンスが問題になりますが、有限要素法以外だったら大丈夫かと。

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# らしく書き換えられないかな、と思案中です。コードが短くなるしΣを使ったほうが、数式を同じになるので。

カテゴリー: F# パーマリンク