Help us understand the problem. What is going on with this article?

OCaml で数値計算

More than 3 years have passed since last update.

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 は行列同士の積だけ並列化が効くらしい

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

対話環境での利用

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

$ 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 から使う話が投稿されていました。

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/Cublasmake && make install からの spoc, spoc_cublas を依存パッケージに加え、cclib オプションで libcublas.dylib をリンクしたら行けるかとこんな感じで

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 を動かそうと思ったんですが、 CUDA 初心者すぎて動かせませんでした、雑魚すぎてごめんなさい。

付録

今回使用したコード

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

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

# 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

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away