OCaml で数値計算



## 1.はじめに

本当は、最近 OCaml 書きはじめた人なので開発環境を紹介するつもりでした。

### ビルドツール

http://omake-japanese.sourceforge.jp/

OMake賢いです。何よりビルド時間の短縮につながります。

### テストツール

http://ounit.forge.ocamlcore.org/api-ounit

OUnit便利。Eclipse でも動くのですね。

### ビルド&テスト自動化

```
omake check -P # check は OUnit テスト時のオプション
```

OMakeが常駐し、自動で変更検出 -> ビルド -> テスト実行します。何がミスったかすぐわかって最高に楽。ただし、今回はこれらを活用する話ではありません。。。


## 2. OCaml で数値計算

BLAS(+Lapack) は古来から伝わる超高速線形代数ライブラリです。その性能と認知度から多くの言語や環境に共通APIが存在します。ATLAS 等をインストールする際のエゲつないパラメータチューニングを目の当たりにすると、伝統の重みを感じざるを得ません。もちろん OCaml にもありました。

http://mmottl.github.io/lacaml/


それでは、OCaml BLASライブラリ・Lacamlと、 C++ のモダンな行列演算 Template ライブラリ Eigen, Boost::numeric::ublas と比較してみます。

### 4096 * 4096 同士の単精度行列積 (sgemm)

|言語|ライブラリ|計算時間|
|:---|:---|---:|
|C|na|510,728,861[ms]|
|C++ |Boost|533,703,624[ms]|
|C++ |Eigen|4,291,730[ms]|
|OCaml|Lacaml + BLAS|942,867[ms]|


[Eigen は行列同士の積だけ並列化が効くらしい](http://eigen.tuxfamily.org/index.php?title=FAQ#How_does_Eigen_compare_to_BLAS.2FLAPACK.3F) 

上記ページで Eigen の人も言ってますが、素の BLAS は使い勝手が悪いです。とはいえ速度は流石です。今回は OSX 標準の ``/System/Library/Frameworks/Accelerate.framework`` が手軽だったので使用しました。探せばもっと早い BLAS 実装も存在すると思います。 http://ibisforest.org/index.php?BLAS

### 対話環境での利用

加えて OCaml なら Lacaml で巧妙にラップされています。C++ にはない対話環境での利用は Python(+NumPy) のように手軽です。

``` ocaml
$ opam install lacaml
$ utop

> #require "lacaml";;

> open Lacaml.S;; (* 単精度演算のお知らせ *)

> let a = Mat.of_array [|[|2.0; 3.0|]; [|1.0; 0.5|]|];;

val a : Lacaml_float32.mat =    C1  C2
                             R1  2   3
                             R2  1 0.5

> let inverse m = let mi = lacpy m in getri mi; mi;; (* 引数に代入してくる関数なので一旦コピー *)

val inverse : Lacaml_float32.mat -> Lacaml_float32.mat = <fun>

> gemm a (inverse a);; (* 逆行列との行列積 *)

- : Lacaml_float32.mat =    C1 C2
                         R1  1  0
                         R2  0  1

> let s, v, d = gesvd (Mat.random 1024 3);; (* 特異値分解 *)
```

タイムリーにも昨日 [C++ を Python から使う話](http://qiita.com/termoshtt/items/81eeb0467d9087958f7f)が投稿されていました。


## 3. おわりに

他の方々とは違う方向で OCaml の表現力について紹介・検証しました。今回はちょっと動かして性能を検証するに留まりましたが、次はこれらを利用して ML (Machine Learning) のプログラムを作れたらいいですね。


## 4. 補足 (挫折したこと)


###cuBLAS で行列演算

ついでに GPGPU です。参考までに nVidia社の cuBLAS 登場です。


+ 4096 * 4096 同士の単精度行列積


|言語|ライブラリ|計算時間|GPUから結果を取り出す時間|
|:---|:---|---:|---:|
|C|na|510,728,861[ms]|na
|C++ |Boost|533,703,624[ms]|na
|C++ |Eigen|4,291,730[ms]|na
|OCaml|Lacaml| 942,867[ms]|na
|C++|cuBLAS| 72[ms]| 146,085[ms]


これはいいですね。一応 OCaml にも GPGPU 系のライブラリがあって動かしてみたのですが...

http://mathiasbourgoin.github.io/SPOC/


``SPOC/SpocLib/Cublas`` で ``make && make install`` からの ``spoc, spoc_cublas`` を依存パッケージに加え、``cclib`` オプションで ``libcublas.dylib`` をリンクしたら行けるかとこんな感じで

``` ocaml:main.ml
open Core.Std
open Core_bench.Std

open Lacaml.S


open Spoc
open Cublas
       
let devices = Spoc.Devices.init ()
let dev = ref devices.(0)
let x = 1
let auto_transfers = ref false  
  
let () =
  Cublas.init ();
  Random.self_init();
  let a = Spoc.Vector.create Spoc.Vector.float32 (x * x)
  and b = Spoc.Vector.create Spoc.Vector.float32 (x * x)
  and res = Spoc.Vector.create Spoc.Vector.float32 (x * x) in
  for i = 0 to (Spoc.Vector.length a - 1) do
    Spoc.Mem.set a i (Random.float 32.);
    Spoc.Mem.set b i (Random.float 32.);
    Spoc.Mem.set res i (Random.float 32.);
  done;

  Spoc.Mem.auto_transfers !auto_transfers;
  if (not !auto_transfers) then
	begin
	  Printf.printf "Transfering Vectors (on Device memory)\n";
	  flush stdout;
	  Spoc.Mem.to_device a !dev;
	  Spoc.Mem.to_device b !dev;
	  Spoc.Mem.to_device res !dev;	    
	end;
  
  begin
    Printf.printf "Computing\n";
	flush stdout;
    Cublas.cublasSgemm 'n' 'n' x x x 1.0 a x b x 0.0 res x !dev;
  end;

  if (not !auto_transfers) then
	begin
	  Printf.printf "Transfering Vectors (on CPU memory)\n";
	  Pervasives.flush stdout;
	  Spoc.Mem.to_cpu res ();
	end;
      
  Cublas.shutdown ();
  Spoc.Devices.flush !dev ();
  
```

``` :実行結果
Transfering Vectors (on Device memory)
Computing
[1]    29511 segmentation fault  ./src/main
```

...。

SPOC の自前 kernel 使ったサンプルは全部動いてたので(古い9600GT搭載マシンだと、こっちもセグフォしちゃったけど...)、nVidia の [matrixMul_kernel.cu](http://docs.nvidia.com/cuda/cuda-samples/index.html#matrix-multiplication--cuda-driver-api-version-) を動かそうと思ったんですが、 CUDA 初心者すぎて動かせませんでした、雑魚すぎてごめんなさい。


## 付録

### 今回使用したコード

https://bitbucket.org/shigekikarita/ocaml-ac2014/src

CUDA は、こんな感じで PATH 通ってる前提です

``` sh
# CUDA
export PATH=/usr/local/cuda/bin:$PATH
export DYLD_LIBRARY_PATH=/usr/local/cuda/lib:$DYLD_LIBRARY_PATH
export LD_LIBRARY_PATH=/usr/local/cuda/lib:$LD_LIBRARY_PATH
export CPLUS_INCLUDE_PATH=/usr/local/cuda/include:/usr/local/cuda/samples/common/inc:$CPLUS_INCLUDE_PATH
export C_INCLUDE_PATH=/usr/local/cuda/include:/usr/local/cuda/samples/common/inc:$C_INCLUDE_PATH
```

Boost と Eigen も入れてください。

### 一応スペック


CPU Intel i7-4770K

GPU Geforce GTX760