F#で数値・線形代数計算をするためのライブラリ紹介(F# PowerPack, F# MathProvider)

に参戦中でして、その12月15日分の記事です。

日本だと未だに「関数型って仕事じゃ使えない」的な話をよく目にしますが、外資系企業、特に私のいる金融業界だと既に関数型での開発が行われていて、F#の使用例としてはクレディスイスという外資系金融機関がF#を使用しているという話は結構有名です*1。私自身はまだF#自体仕事で使えていないのですが、ゆくゆくはF#を使って数値計算モデリング的な事をやっていきたいな考えていて、今日はそれに関して調べたことを書こうと思います。

F#で使える数値計算用のライブラリ

当然「F#を使ってゼロからライブラリ開発!」なんて元気はあまりない&.NETのものがそのまま使えるというメリットを活かさない理由がないということで、既存ライブラリをまずは調査。この辺は結構豊富で、以下のようなものがみつかりました。

上を書きあげてから気がついたのですが、Overview: Numerical Libraries for F# and .NET Framework | Microsoft Docsでは*2、以下のように各種数値計算用のライブラリとその特徴を表としてまとめてくれていました。

Library License Linear Algebra Optimization FFT Statistics Designed for F#
Alglib Dual (GPL, Commercial) Managed Yes Yes Yes No
Math.NET Numerics MIT Managed, Native No No Yes Yes
DotNumerics GPL Managed Yes No No No
F# PowerPack and F# MathProvider Apache 2.0 Managed, Native No No No Yes
Microsoft Sho Noncommercial Native Yes Yes Yes No
Microsoft Solver Foundation Free, Commercial, Academic N/A Yes No No No
Extreme Optimization Commercial Managed, Native Yes Yes No Samples
F# for Numerics, F# for Visualization Commercial Managed No Yes Yes Yes
NMath Commercial Native Yes Yes No No

数値計算に使えそうなF#のライブラリ、結構たくさんありますね!
今回はその中でもF# PowerPackF# MathProviderを紹介したいと思います。

F# PowerPack

F# PowerPackは皆さん既に使っていそうですがの中にある機能をご紹介したいと思います。
インストール自体は上記のLINKから行ってダウンロードをした後、「インストーラー起動」⇒「Nextを続けて選択」でOKでした。
F# PowerPackにはチュートリアルがありますので、まずはこれを参考にするのが最も良いと思います。

では早速スクリプトを書いていきましょう。まずは準備としてライブラリへのリファレンス(参照)を張っておきます。

//PowerPackへのリファレンスを追加
#r "FSharp.PowerPack.dll"
複素数

complex関数を使うと複素数型の変数束縛が可能になります。ここではシンプルに"-1"を複素数型として定義して、それのルートを計算させてみます。

  • 1のルートなんでこれはi(虚数単位)になることが期待されますが、実際ちゃんとそうなっていますよと、そして最後に虚数単位を二乗して-1に戻るかを確認しました。
//-1を複素数型として定義
let x = complex -1.0 0.0;;
//xの平方根(虚数単位)
let unit_imaginary = sqrt x;;
unit_imaginary * unit_imaginary;;

実行結果は

val x : complex = -1r+0i
val unit_imaginary : Complex = 6.12303176911189e-17r+1i
val it : Complex =
  -1r+1.22460635382238e-16i {Conjugate = -1r-1.22460635382238e-16i;
                             ImaginaryPart = 1.224606354e-16;
                             Magnitude = 1.0;
                             Phase = 3.141592654;
                             RealPart = -1.0;
                             i = 1.224606354e-16;
                             r = -1.0;}

となりバッチリ、最終的に虚数単位を2乗して-1に答えが戻ってくる様がわかります。

行列

PowerPackを使うと行列も定義できて、単位行列や零行列等基本的なものに関してはMatrixモジュールを使うことで簡単に作成できるようになっています。

//零行列
Matrix.zero 3 3;;
//単位行列
Matrix.identity 3;;
//行列を定義
let sample_matrix = matrix [ [ 1.; 2.];[3.; 4.]];;

実行結果はそれぞれ

val it : matrix = matrix [[0.0; 0.0; 0.0]
                          [0.0; 0.0; 0.0]
                          [0.0; 0.0; 0.0]]
val it : matrix = matrix [[1.0; 0.0; 0.0]
                          [0.0; 1.0; 0.0]
                          [0.0; 0.0; 1.0]]
val sample_matrix : matrix = matrix [[1.0; 2.0]
                                     [3.0; 4.0]]

となります。

行列の初期化に関してもう少し言うと、Matrix.init関数は更に強力で行列の”行”・”列”を引数にしたlamda式を使って色々な行列を作成することができます。ここではサンプルとして下三角行列を作成してみます。

> Matrix.init 3 3 (fun row col -> if row >= col then 1.0 else 0.0);;
val it : matrix = matrix [[1.0; 0.0; 0.0]
                          [1.0; 1.0; 0.0]
                          [1.0; 1.0; 1.0]]

以降、上で作った”sample_matrix”行列を例に色々な行列処理を施してみます。

まず個別の要素にアクセスするには

> sample_matrix.[0, 0];;
val it : float = 1.0
> sample_matrix.[0, 1];;
val it : float = 2.0
> sample_matrix.[0 .. 0, 0 .. 1];;
val it : Matrix<float> = matrix [[1.0; 2.0]]

と取得したい行と列をそれぞれ指定する".[row, col]"のようなスタイルで書きます。例にあるように行列のインデックスは0始りで、".."を使い複数行・列に渡って同時にアクセス可能です。

行列の特定の各行・列のみを抽出したいときにはRow・Columnメンバ関数をそれぞれ使ってもOKです*3

> sample_matrix.Row(0);;
val it : RowVector<float> = rowvec [|1.0; 2.0|]
> sample_matrix.Column(0);;
val it : Vector<float> = vector [|1.0; 3.0|]

ただ型がMatrix型からVector型に変わってしまう点が問題でしょうか・・・
Vector型の有難味がまだよく解っていないので、どちらが良いのかはもうちょい実践経験を積まないとなんとも言えませんorz

行列の要素の値を変更するには

> sample_matrix.[0, 0] <- 100.0;;
val it : unit = ()
> sample_matrix;;
val it : matrix = matrix [[100.0; 2.0]
                          [3.0; 4.0]]
> sample_matrix.[0, 0] <- 1.0;;
val it : unit = ()
> sample_matrix;;
val it : matrix = matrix [[1.0; 2.0]
                          [3.0; 4.0]]

のように"<-"を使用して変更します。mutableな変数と同じですね。

行列同士の基本的な四則演算の機能も当然存在して、

> sample_matrix + sample_matrix;;
val it : Matrix<float> = matrix [[2.0; 4.0]
                                 [6.0; 8.0]]
> sample_matrix - sample_matrix;;
val it : Matrix<float> = matrix [[0.0; 0.0]
                                 [0.0; 0.0]]
> sample_matrix * sample_matrix;;
val it : Matrix<float> = matrix [[7.0; 10.0]
                                 [15.0; 22.0]]
> sample_matrix .* sample_matrix;;
val it : Matrix<float> = matrix [[1.0; 4.0]
                                 [9.0; 16.0]]

のように、上から順に足し算(+)・引き算(-)・行列の掛け算(*)・行列の要素ごとの掛け算(.*)となります。

行列の転置・総和・対角和を計算するための関数もあります。*4

> Matrix.transpose sample_matrix;;
val it : matrix = matrix [[1.0; 3.0]
                          [2.0; 4.0]]
> Matrix.sum sample_matrix;;
val it : float = 10.0
> Matrix.trace sample_matrix;;
val it : float = 5.0

行列の次元数・行サイズ・列サイズを抽出するにはそれぞれメンバ関数Dimensions・NumRows・NumColsを使います。

> sample_matrix.Dimensions;;
val it : int * int = (2, 2)
> sample_matrix.NumRows;;
val it : int = 2
> sample_matrix.NumCols;;
val it : int = 2

F#ではおなじみのmap・mapi・fold*5処理を実行させることもできます。以下の例で、map処理の方では行列の全要素を2倍してます*6。mapi処理の方では行列を上三角の要素だけ残す処理、foldでは全要素の積を計算してます。

> sample_matrix |> Matrix.map (fun x -> x * 2.0);;
val it : matrix = matrix [[2.0; 4.0]
                          [6.0; 8.0]]
> sample_matrix |> Matrix.mapi (fun i j value -> if i > j then 0.0 else value) ;;
val it : matrix = matrix [[1.0; 2.0]
                          [0.0; 4.0]]
> sample_matrix |> Matrix.fold (fun x y -> x * y) 1.0;;
val it : float = 24.0

行列の要素がある特定の条件を満たしているかどうかはMatrix.forall,Matrix.exists関数を使用します。
これは行列の個々の要素に対して指定したlambda式がtrueを返すか否かを判定して、

  • forallの場合:全部
  • exitstの場合:少なくとも1つ

の要素でtrueが帰ってくればtrueを返すような関数です。他言語にあるall,anyみたいな奴ということです。

> sample_matrix |> Matrix.forall (fun x -> if x > 1.0 then true else false);;
val it : bool = false
> sample_matrix |> Matrix.exists (fun x -> if x > 1.0 then true else false);;
val it : bool = true
> sample_matrix |> Matrix.exists (fun x -> if x > 1.0 then true else false);;

F# MathProvider

次にF# MathProviderを紹介します。F# MathProviderは線形代数計算ライブラリであるBLASLAPACKへのラッパーです。このライブラリで線形代数計算を行う際の行列型はPowerPackのものが使われるのでそちらも予めインストールしておく必要があります。これのインストール自体は簡単でF# MathProviderからダウンロードしたzipファイルを解凍し、適当なフォルダにおいておけばOKです。私はCドライブ直下に保存しました。

これを使うにはまずライブラリへのリファレンスを張る必要があるので、以下のようにPowerPackも含めてリファレンスを張っておきます。

> #r "FSharp.PowerPack.dll";;

--> Referenced 'C:\Program Files\FSharpPowerPack-2.0.0.0\bin\FSharp.PowerPack.dll'

> #r "C:\MathProvider\Net 4.0\MathProvider.dll";;

--> Referenced 'C:\MathProvider\Net 4.0\MathProvider.dll'

次にMathProviderについてくるlapack.dllとblas.dllを使いたいのでパスの通ってるフォルダにそれらをコピーしておくか、以下のようにカレントの作業フォルダをdllが入っているフォルダに変更しておきます。

> System.Environment.CurrentDirectory <- @"C:\MathProvider\LapackRuntime";;
val it : unit = ()

このままでも以下に記述したサンプルコードは動きますが、固有値計算・SVD分解等、関数によっては

System.Exception: Not yet implemented, managed fallback linear algebra ops coming soon

のように「未実装だよ!」という旨のエラーがでます。これはlapackblasへのProviderを起動させないとManaged F#での実装が使用されるためで、そちらでの実装がまだ終わっていないことに起因しているようです。Nativeなlapack実装で計算させるためには以下のようにstartProviderメソッドを実行しておきます。

> MathProvider.LinearAlgebra.startProvider();;
Fatlab: dense linear algebra service = "Native Lapack"
val it : unit = ()

以下の使用例では以下のような行列を例に、MathProviderから提供される行列処理を行ってみます。

> let A = matrix [ [3.; 1.0]; [2.; 5.] ];;

val A : matrix = matrix [[3.0; 1.0]
                         [2.0; 5.0]]
行列式
> let det = MathProvider.LinearAlgebra.det A;;

val det : float = 13.0
逆行列
> let inv = MathProvider.LinearAlgebra.inv A;;

val inv : matrix = matrix [[0.3846153846; -0.07692307692]
                           [-0.1538461538; 0.2307692308]]
コレスキー分解
> let ch  = MathProvider.LinearAlgebra.chol A;;

val ch : matrix = matrix [[1.732050808; 1.154700538]
                          [0.0; 1.914854216]]
LU分解
> let p, l, u  = MathProvider.LinearAlgebra.lu  A;;

val u : matrix = matrix [[3.0; 1.0]
                         [0.0; 4.333333333]]
val p : (int -> int)
val l : matrix = matrix [[1.0; 0.0]
                         [0.6666666667; 1.0]]
QR分解
> let q, r  = MathProvider.LinearAlgebra.qr  A;;

val r : matrix = matrix [[-3.605551275; -3.605551275]
                         [4.440892099e-16; -3.605551275]]
val q : matrix = matrix [[-0.8320502943; 0.5547001962]
                         [-0.5547001962; -0.8320502943]]
SVD分解
> let v, s, ut = MathProvider.LinearAlgebra.svd A;;

val v : matrix = matrix [[-0.4161611044; -0.9092908969]
                         [-0.9092908969; 0.4161611044]]
val ut : matrix = matrix [[-0.5257311121; -0.8506508084]
                          [-0.8506508084; 0.5257311121]]
val s : Vector<float> = vector [|5.833904512; 2.228353236|]
固有値固有ベクトル計算
> let a, b = MathProvider.LinearAlgebra.eigenSym A;;

val b : Matrix<float> = matrix [[-0.9238795325; 0.3826834324]
                                [0.3826834324; 0.9238795325]]
val a : vector = vector [|2.585786438; 5.414213562|]

まとめ

今回はF# PowerPack・F# MathProviderに的を絞って特に行列の扱いについて重点的に機能を紹介しました。これは数値計算の基礎の基礎的な部分なので、今後は上で紹介したShoやMath.NETを使ってより実践的な内容をやっていきたいと思ってます。ここで紹介仕切れなかった行列型の機能については参考LINKの「Type Microsoft.FSharp.Math.Matrix」や「Module Microsoft.FSharp.Math.Matrix」を見てください。スパース行列等も作れるようです。

*1:参考:http://www.microsoft.com/casestudies/Microsoft-Visual-Studio-2010/Financial-Services-Firm/Banking-Firm-Uses-Functional-Language-to-Speed-Development-by-50-Percent/4000006794

*2:ちょい古めかもです

*3:上で見たように".."を使ってもできます

*4:転置に関してはsample_matrix.Transposeでも計算できた

*5:foldiもあり

*6:単純に2.0を行列に掛け算してもOKです、そっちの方が早いし楽…