F#で数値・線形代数計算をするためのライブラリ紹介（F# PowerPack, F# MathProvider）
http://d.hatena.ne.jp/teramonagi/20111215/1323874810#20111215fn6
を見ると、F# PowerPack と F# MathProvider を使うと、逆行列が簡単に解けます。
#r "FSharp.PowerPack.dll";
#r @"D:toolsMathProviderNet 4.5MathProvider.dll";
let A = matrix [[3.0; 1.0;];
[2.0; 5.0;]]
let det = MathProvider.LinearAlgebra.det A
let inv = MathProvider.LinearAlgebra.inv A
こんな風にすると、
val A : matrix = matrix [[3.0; 1.0]
[2.0; 5.0]]
val det : float = 13.0
val inv : matrix = matrix [[0.3846153846; -0.07692307692]
[-0.1538461538; 0.2307692308]]
と出ます。めでたしめでたし…なのは、後から知ったわけで、実は F# MathProvider は必ず Blas.dll と lapack.dll が必要なのかと思って躊躇してたんですよね。できることならば、F# だけで作りたい、というか、ストアアプリとか自前で高速化（有限要素法で1万節点とか）を考えると、中身を理解しておきたい、という願望があったので…車輪の再発明をしてしまいました。
let inv ( A' : matrix ) =
let A = A'.Copy()
let n = A.NumCols-1
let B = Matrix.identity A.NumCols
for k in 0..n do
// akkを1にする
let p = A.[k,k]
for j in 0..n do
A.[k,j] <- A.[k,j] / p
B.[k,j] <- B.[k,j] / p
// k行以外から引く
for i in 0..n do
if i <> k then
let d = A.[i,k]
for j in 0..n do
A.[i,j] <- A.[i,j] - A.[k,j] * d
B.[i,j] <- B.[i,j] - B.[k,j] * d
B
逆行列の計算
http://www.asahi-net.or.jp/~uc3k-ymd/Lesson/Section03/invmat.html
を参考にしてべたべたな手計算の方法をF#に直しています。結果が0になるところも計算してしまって無駄が多いので、高速化のし甲斐がある…というか、Fortran とか C言語とかのコードもあちこちにあるので、それらをコンバートしたほうが良いかも。
ただし、有限要素法で使う逆行列の場合には、対称性をうまく利用したり計算量を減らすことができるので、一般的な逆行列よりも手を加えたほうが実用的かもしれません。が、ADVENTURE のように大規模な構造（１億節点など）とは違って、節点数のオーダーが違えば、最近のPCならばさして高速化せずとも結構なところまでいけるのでは？と考えています。これは先行き試してみたいところです。
実はこのF#のコード、半日ほどかかったんですよね。なかなか数値が合わなくて苦労したのは、F#の文法に馴れていなかったためもあるのですが、既存の Fortran コードとは別の形で実装して手計算そのままを移したためです。まあ、勉強になったからいいか。正月だし。
これが「できたーっ!!!」後に気付いたのが、F# MathProvider の実装です。Blas.dll を使わない場合には、自前の F# コードを使うようにできているんですね。なるほど。
以下が抜粋、LU分解と上三角を使っています。
/// Given A[n,n] find it's inverse.
/// This call may fail.
let Inverse a =
if HaveService() then LinearAlgebraService.inverse a
else LinearAlgebraManaged.Inverse a
let Inverse A =
let (n,m) = matrixDims A
if n <> m then invalid_arg "Matrix must be square when computing its inverse."
let P,L,U = LU A
(SolveTriangularLinearSystems U (SolveTriangularLinearSystems L ((Matrix.identity n).PermuteRows P) true) false)
let LU A =
let (nA,mA) = matrixDims A
if nA<>mA then invalid_arg "lu: not square";
let L = Matrix.zero nA nA
let U = Matrix.copy A
let P = [| 0 .. nA-1 |]
let abs (x:float) = System.Math.Abs x
let swapR X i j = // REVIEW should we make this a general method?
let (nX,mX) = matrixDims X
let t = X.[i .. i,0 .. ]
for k=0 to mX-1 do
X.[i,k] <- X.[j,k]
X.[j,k] <- t.[0,k]
done
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
(((*P.Length,*)Permutation.ofArray P),L + Matrix.identity nA,U)
let SolveTriangularLinearSystems K B isLower =
if isLower then
let (nK,mK) = matrixDims K
let (nB,mB) = matrixDims B
if nK<>mK || nB<> nK then invalid_arg "Cannot do backward substitution on non-square matrices."
let X = Matrix.zero nK mK
for i=0 to nK-1 do
for k=0 to mB-1 do
let s = ref B.[i,k]
for j=0 to i-1 do
s := !s - K.[i,j] * X.[j,k]
done
X.[i,k] <- !s / K.[i,i]
done
done
X
else
let (nK,mK) = matrixDims K
let (nB,mB) = matrixDims B
if nK<>mK || nB<> nK then invalid_arg "Cannot do backward substitution on non-square matrices."
let X = Matrix.zero nK mK
for i=0 to nK-1 do
for k=0 to mB-1 do
let s = ref B.[nK-i-1,k]
for j=0 to i-1 do
s := !s - K.[nK-i-1,nK-j-1] * X.[nK-j-1,k]
done
X.[nK-i-1,k] <- !s / K.[nK-i-1,nK-i-1]
done
done
X
まあ、数式は理解して使ったほうがいいので、しばらくは自前で実装ということで。
