行列計算ライブラリの簡易ベンチマーク

ClojureとF#は行列計算をネイティブコードで実行するライブラリがあるのでそれを使うとどのくらい速いのかをテストしてみました。

3000x3000で要素は[0,1]区間の一様乱数で埋め尽くした行列Mを作ってMとMの逆行列を掛けてトレースを計算する(答えは3000になるはず)という内容です。

条件はMacBookPro late 2013 16GB Memory Win7 Pro x64です。

すべてのライブラリをカバーしてないのと、実行環境によって変化が大きそうなので冒頭での結論提示はしません。参考情報ということで以下でご参照ください。

R (3.0.3 x64)

まずベースラインとしてR言語で。

N <- 3000
 
pt <- proc.time()
m <- matrix(runif(N*N), nrow=N, ncol=N)
ans <- sum(diag(m %*% solve(m)))
pt <- proc.time()-pt
 
print(paste("ans =",ans))
print(paste("time = ",pt[3]))

何回かやってみて、だいたい平均が43秒でした。ちなみに以下のF#とClojureではもっと速くなる結果が出てますが、F#やClojueの言語そのもので書くとRよりはるかに遅くなります(実際にその言語で書かれたライブラリでもやりました)。

F# (Math.NET Numerics + MKL x64)

open MathNet.Numerics.LinearAlgebra
open MathNet.Numerics.LinearAlgebra.Double
open MathNet.Numerics.Distributions
 
open System
open System.Diagnostics
 
MathNet.Numerics.Control.LinearAlgebraProvider <-
   new MathNet.Numerics.Algorithms.LinearAlgebra.Mkl.MklLinearAlgebraProvider()
 
let test n =
    let m : Matrix = 
        upcast DenseMatrix.randomCreate n n (ContinuousUniform (0.0,1.0))
    (m * m.Inverse()).Trace()
 
[<EntryPoint>]
let main argv = 
    let sw = new Stopwatch()
    sw.Start()
    let ans = test 3000
    sw.Stop()
    printfn "%A" ans
    printfn "elapsed = %d ms" <| sw.ElapsedMilliseconds
    0

結果はだいたい10秒程度でした。Rよりだいぶ速い。

.NETは商用のものも多いですがフリーのライブラリはMath.NET Numericsという共通インタフェースにMKLというネイティブライブラリを入れたものがDonSyme著「Expert F# 3.0」によると良いらしいのでそれで試しました。

Clojure (JBLASベースのClatrix)

(ns clojure-matrix-bench.core
  (:use clatrix.core))
 
(defn- my-main []
  (let [n 3000
        m (rand n n)]
    (trace (* m (i m)))))
 
(defn -main []
  (print (time my-main)))

結果はだいたい3〜5秒でした(計測ミス)47〜50秒でした。

Clojureもcore.matrixという共通インタフェースが定義されててそれを各ネイティブラッパーが実装提供する構造になってます。

JVM系はフリーのいいライブラリがいっぱいあるのが魅力です。ここに一覧が。あと最近でたっぽいこれとかも。

まとめ

僕のPCはハイパースレッドつき2コアなので物理スレッドは4個ですが、上記のプログラムのいずれも実行中はWindowsのCPUモニタを見るとどうも逆行列計算では1/4しか使っていないっぽかった。(一部、traceのところで50%に行く時があったので単純な和の場合は並列演算を利かすのでしょう。逆行列計算は途中がモノイドにならないということなのかな)。あとメモリは1億足らずの要素数だからか、全然余裕でした。

というわけでJBlasベースのclatrixがIntel MKLより速かったのがちょっと意外。(計測ミス) Intel MKLは速いですね。

ソースは(こちら)。後でほかのライブラリ(Parallel ColtとEJML)や、ほかの言語のPythonのNumPyとかJulia(たしかBLASのはず)とかでも試してみます。

なにか意味のある結果を出しました的なエントリじゃなくて恐縮ですが、「関数型言語+データ分析関連ライブラリ」での使える武器を増やしまくりたいという趣旨でもあります。

あと、乱数生成(特にR)が遅い可能性あるという指摘をいただいたので行列演算のみで計測しなおしましたがいずれも1%も変化しませんでした。ただそれは一様乱数だからで、正規乱数とかはパフォーマンスに気を付けたほうがいいかもです。