JuliaTokyo #2 でモンテカルロ法について話してきました

9月27日の JuliaTokyo #2 で『Juliaで学ぶ Hamiltonian Monte Carlo (NUTS 入り)』と題してHamiltonian Monte Carlo (HMC) とStanで使われているNUTSについてお話してきました。

JuliaTokyo #2 - connpass

内容

Juliaで学ぶ Hamiltonian Monte Carlo (NUTS 入り) from Kenta Sato

サンプルコード: bicycle1885/JuliaTokyo2HMC · GitHub

MCMC自体については30分ではほとんど説明できませんでしたが、会場のほとんどの方が既に知っているようでしたのでちょうど良かったようです。

全体の流れとしては、Metropolis-Hastings、Hamiltonian Monte Carlo、No-U-Turn Samplerの順に3つサンプラーを実装コードと正規分布からのサンプリングを例に紹介し、どのように前のサンプラーの弱点を克服していってるかを説明しました。

Unicodeプログラミング

サンプルコードはJuliaの愉しみのひとつであるUnicode文字をふんだんに使っていたところが面白いところです。

function build_tree(L::Function, ∇L::Function, θ::Vector{Float64}, r::Vector{Float64}, u::Float64, v::Int, j::Int, ϵ::Float64)
    if j == 0
        θ′, r′ = leapfrog(∇L, θ, r, v * ϵ)
        C′ = u ≤ exp(L(θ′) - r′ ⋅ r′ / 2) ? Set([(θ′, r′)]) : Set([])
        s′ = int(L(θ′) - r′ ⋅ r′ / 2 > log(u) - Δmax)
        return θ′, r′, θ′, r′, C′, s′
    else
        θ⁻, r⁻, θ⁺, r⁺, C′, s′ = build_tree(L, ∇L, θ, r, u, v, j - 1, ϵ)
        if v == -1
            θ⁻, r⁻, _, _, C″, s″ = build_tree(L, ∇L, θ⁻, r⁻, u, v, j - 1, ϵ)
        else
            _, _, θ⁺, r⁺, C″, s″ = build_tree(L, ∇L, θ⁺, r⁺, u, v, j - 1, ϵ)
        end
        s′ = s′ * s″ * ((θ⁺ - θ⁻) ⋅ r⁻ ≥ 0) * ((θ⁺ - θ⁻) ⋅ r⁺ ≥ 0)
        C′ = C′ ∪ C″
        return θ⁻, r⁻, θ⁺, r⁺, C′, s′
    end
end

ちなみにjulia-vimはlatex記法サポートしてます https://t.co/iT48OQLi6r #JuliaTokyo

— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 9月 27

箴言

Knuth 「尚早な最適化は諸悪の根源」 M. Jackson 「最適化第一法則: 最適化するな, 第二法則: まだするな」 #JuliaTokyo

— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 9月 27

+α

俺のRandomForests.jlが某巨大企業で使われているらしい！凄い！ #JuliaTokyo

— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 9月 27

社内で改造しているようですが、手続き上のアレでプルリクは受けられないようです(´・ω・`)

コマンドを並列に実行するGNU parallelがとても便利

最近のコンピュータは複数のCPUコアを持っているので並列にコマンドを実行することができます。 たくさんの同じようなファイルに同じ処理を実行することは、私のやっているバイオインフォマティクスではよくあります。

しかし自分で並列に実行するスクリプトを書くことはそれほど簡単ではなく、ログや実行結果の確認など煩雑な処理を書かなければいけません。 この記事では、そうした処理を簡単にするGNU parallelというツールを紹介します。

GNU parallel

UNIX系のOSではインストールはとても簡単です。MacでしたらHomebrewを使って、Linuxでは各ディストリビューションのパッケージマネージャからインストールできます。 詳しくはGNU parallelのウェブページを参照して下さい(http://www.gnu.org/software/parallel/)。

Homebrew:

brew install parallel

パッケージマネージャが使えない環境でも、ソースコードからビルドすることができます。 ダウンロードページ(http://ftp.gnu.org/gnu/parallel/)から、最新のparallel(parallel-latest.tar.bz2)をダウンロードして展開し、ビルドしましょう。 必要ならばconfigureに--prefixも渡せます。

Source:

./configure && make && make install

なお、最初にparallelコマンドを実行する際に引用のお願いメッセージが出ますが、これはparallel --bibtexを一度実行すれば抑えられます。

Hello, GNU parallel !

インストールが終わったら、早速GNU parallelを実行してみましょう。

~/t/sample $ parallel echo ::: hello world !
hello
!
world

GNU parallelのparallelコマンドは、実行するコマンドと、そのコマンドに適用するいくつかの引数を引数としてとります。 上の例では、echoが実行するコマンドで、:::の後に続く3つの引数hello, world, !がそれぞれechoコマンドに適用する引数です。 echoされたメッセージが元の引数の順番と一致していないことから分かる通り、echoコマンドは並列に実行されるため、各引数の実行順は不定になります。

:::は複数つなげることで、引数の組み合わせの積をつくることもできます。

~/t/sample $ parallel echo ::: hello, bye, ::: Alice Bob Charlie ::: !
hello, Alice !
hello, Bob !
hello, Charlie !
bye, Alice !
bye, Bob !
bye, Charlie !

ファイルや標準入出力から引数を与える

以下のようなファイル名のリストを収めたファイル(list.txt)があるとします。

~/t/sample $ cat list.txt
foo.txt
bar.txt
baz.txt

-a <file>オプションで、ファイルの各行を引数としてコマンドを並列に走らせることができます。 list.txtの各ファイルを並列に圧縮するには以下のようにします。

~/t/sample $ ls
bar.txt  baz.txt  foo.txt  list.txt
~/t/sample $ parallel -a list.txt gzip
~/t/sample $ ls
bar.txt.gz baz.txt.gz foo.txt.gz list.txt

パイプで引数を渡すこともできます。その時は-a -とオプションを渡します。

~/t/sample $ cat list.txt | parallel -a - gzip
~/t/sample $ ls
bar.txt.gz baz.txt.gz foo.txt.gz list.txt

引数の区切り文字はデフォルトでは\n(newline)です。 -0(--null)オプションで区切り文字が\0(null)文字になります。

他にも、<や::::を使って引数の納められたファイルを渡すことができます。

~/t/sample $ parallel gzip < list.txt
~/t/sample $ parallel gzip :::: list.txt

引数の置き換え

コマンドに渡す引数の位置が最後でない場合や、特別な処理をしたい場合に引数の置換え位置を指定するをする必要があります。 その際には{}を使います。

~/t/sample $ parallel 'find {} -name "README*"' ::: ~/vendor/julia ~/vendor/vector
/Users/kenta/vendor/vector//old-testsuite/microsuite/README
/Users/kenta/vendor/vector//README.md
/Users/kenta/vendor/julia//contrib/mac/app/README
/Users/kenta/vendor/julia//contrib/README.ackrc.txt
...

さらに、{.}でファイル名の拡張子を除いた引数にしたり、{/}でファイル名だけ取り出したりできます。

~/t/sample $ parallel 'echo {.}' ::: tmp/foo.txt.gz tmp/bar.txt.gz
tmp/foo.txt
tmp/bar.txt
~/t/sample $ parallel 'echo {/}' ::: tmp/foo.txt.gz tmp/bar.txt.gz
foo.txt.gz
bar.txt.gz

実行されるコマンドの確認

--dry-runをつかえば、引数がどのようにコマンドに適用されてどのようなコマンドが実際に走るかを確認できます。 実際に実行する前に一度確認すると良いでしょう。

~/t/sample $ cat list.txt | parallel -a - --dry-run gzip
gzip foo.txt
gzip bar.txt
gzip baz.txt

実行結果の取得

標準出力と標準エラーに吐かれた結果をファイルとして取得することも容易です。 --results <outputdir>とすることで出力を各引数毎にディレクトリに構造化して保存できます。

~/t/sample $ parallel --results results 'perl -E "say STDOUT \"stdout\"; say STDERR \"stderr\""' ::: A B C
stdout
stderr
stdout
stderr
stdout
stderr
~/t/sample $ tree results/
results/
└── 1
    ├── A
    │   ├── stderr
    │   └── stdout
    ├── B
    │   ├── stderr
    │   └── stdout
    └── C
        ├── stderr
        └── stdout

4 directories, 6 files
~/t/sample $ cat results/1/A/stderr
stderr
~/t/sample $ cat results/1/A/stdout
stdout

実行結果の確認

時間のかかるコマンドや引数が多い場合など、実行結果の確認が大変な場合があります。 そこで、--joblog <logfile>を使えばコマンドのexit statusが確認しやすくなります。

~/t/sample $ parallel --joblog joblog.txt 'sleep {}; exit {}' ::: 0 1 2 3
~/t/sample $ cat joblog.txt
Seq     Host    Starttime       JobRuntime      Send    Receive Exitval Signal  Command
1       :       1407646480.280       0.072      0       0       0       0       sleep 0; exit 0
2       :       1407646480.284       1.144      0       0       1       0       sleep 1; exit 1
3       :       1407646480.289       2.125      0       0       2       0       sleep 2; exit 2
4       :       1407646480.294       3.153      0       0       3       0       sleep 3; exit 3

一目見て分かる通り、実行時間や実際に実行されたコマンドなども記録されるため、問題が発生した時のトラブルシューティングに非常に役に立ちます。 また、終わったコマンドから順に追記されていくため、進捗状況の確認もできます。

並列ジョブの制御

ファイルのダウンロードなどを並列に行いたいが過剰に負荷を掛けたくない場合や、逆にCPU boundでない処理をCPUのコア数より多く並列させたいときなどに同時実行するジョブ数を指定できると便利です。 デフォルトでは、マシンのコア数分だけ並列させるため、共有サーバーなどではCPUを占拠してしまい迷惑になることも考えられます。 そのようなときは--jobs <N>(-j <N>)オプションを使って、個分だけ並列化させる必要があります。

以下のように6つの引数を1, 2, 3, 6並列でそれぞれ実行して実行時間を見てみましょう。

~/t/sample $ time parallel --jobs 1 'sleep {}' ::: 1 1 1 1 1 1
        7.06 real         0.26 user         0.10 sys
~/t/sample $ time parallel --jobs 2 'sleep {}' ::: 1 1 1 1 1 1
        3.58 real         0.27 user         0.11 sys
~/t/sample $ time parallel --jobs 3 'sleep {}' ::: 1 1 1 1 1 1
        2.44 real         0.28 user         0.11 sys
~/t/sample $ time parallel --jobs 6 'sleep {}' ::: 1 1 1 1 1 1
        1.30 real         0.30 user         0.12 sys

--jobs 50%などと書くことで全体の半分のCPUコアを使う指定もできます。

~/t/sample $ time parallel --jobs 50% 'sleep {}' ::: 1 1 1 1 1 1
        3.59 real         0.26 user         0.11 sys

また、メモリを沢山必要とする計算などでは、--noswapオプションを指定することで、メモリのスワップが発生している時には新しいジョブを実行しないようにすることもできます。

他の便利そうな機能

ここで紹介した機能はGNU parallelのごく一部で、他にも有用なオプションなどがたくさんあります。 一度、man parallel_tutorialやman parallelに目を通すことをオススメします。 私はまだ実際に実行して確認していませんが、GNU parallelにはSSH越しに複数のノードで並列にジョブを実行させる機能やパイプライン処理で中間のボトルネックになっている処理だけ並列化するなどの機能もあるようです。 これらの機能は確認し次第、ブログに追記していこうと思います。

statisticsにt検定を実装しました

ちょっとHaskellのstatisticsパッケージを見ていたら、統計的仮説検定に何故かt検定が無かったのでチョイチョイと実装しました。

https://github.com/bos/statistics/pull/66

masterの一歩手前のconfidenceブランチに取り込まれましたので、そのうち使えるようになると思います。

この実装で以下の3つの2標本検定が使えるようになります。

• Student's t-test
• Welch's t-test
• paired samples t-test

詳細については、Wikipediaを参照下さい。

statisticsパッケージの検定のAPIがに統一感がなくてイケてない感じがするので、なにか思いついたらそちらも修正したいです。あとp値や統計量が得られないのはさすがに良くない気がします。

JuliaTokyo #1 を開催しました

日本で(多分)最初のJuliaの勉強会 JuliaTokyo #1 を開催しました。イベントページはこちらです。 http://juliatokyo.connpass.com/event/6891/

当日は40人ほどの参加者が集まり、いつの間にJuliaがそんなに一般化したんだという感じでしたが、Julia歴1週間とか当日の朝から始めましたというように、何処からかJuliaが熱いらしいとの噂を聞きつけてJuliaを使ってみようという新しいもの好きの人たちが多い印象でした。

私はパッケージの作り方についてお話をしてきた次第です。

以下は、それぞれの発表と私のひとこと感想です。

メインセッション

Juliaのこれまでとこれから - @sorami

Juliaの背景的な情報やどんなエコシステムになっているかを紹介していました。 JuliaはPythonやRからユーザを取り込みつつ、共に進化していく道を辿るようです。

Julia100本ノック - @chezou

Julia 100 exercises #JuliaTokyo from chezou

Juliaでnumpyの100 numpy exercisesを書いた話です。 Juliaのやり方とnumpyのやり方はベクトル化などの使い方がかなり違い、苦労があったようです。 あとJuliaのベクトル内包表記は速いのでガンガン使っていきましょう！

Juliaのパッケージを作ろう! - @bicycle1885

Juliaのパッケージをつくろう! from Kenta Sato

私の発表です。 Juliaのパッケージをどのようにつくるかを解説しました。 パッケージ名やディレクトリ構成などについて、ドキュメントに書いてない慣習などが大事だったりします。

ビジュアライズJulia - @nezuq

ビジュアライズJulia from nezuQ

"The Grammar of Graphics"のlayerやscaleの概念や、Gadfly.jlとの関係の話でした。 可視化関係のライブラリはまだそれほどなので、今後充実してくると良いですね。

メカ女子将棋 - @kimrin

メカ女子将棋Julia tokyo#1 from Takeshi Kimura

Juliaで将棋AIを作った話でした。やっぱりJuliaはCなどに比べ開発がずっとしやすいですね。 盤面の表現にはJuliaの128bit整数を使っているようです。

A Brief JuliaCon Report / Natural Language Processing with Julia - Pontus Stenetorpさん

https://github.com/JuliaCon/presentations/blob/master/JuliaNLP/JuliaNLP.pdf

先月Chicagoで開催されたJuliaConのレポートと、そのときに発表されたJuliaでの自然言語処理のお話をしていただきました。 JuliaConでは、アカデミックの人がメインだったものの、実際にサービスにJuliaを投入している会社もあったようです。 Juliaで自然言語のパーサallenを作って、パフォーマンスも満足いくものになったようです。

LT

Juliaでちゃんとしたコードを書く - 関さん

当日の朝からJuliaを始めて、どうするとJuliaが遅くなるかのお話でした。 マニュアルのperformance tipsマジ重要です。読みましょう!!

Julia0.3でランダムフォレスト - @gepuro

Julia0.3でランダムフォレスト from Atsushi Hayakawa

DecisionTree.jlを使ってみたようです。 DecisionTree.jlはID3を決定木構築のアルゴリズムに使っているので、CARTが使いたければ私のRandomForests.jlを使いましょう！

Plotly Julia API - Yasunobu Igarashiさん

Plotly Julia API from Yasunobu Igarashi

Plotlyというデータの可視化や共有のためのプラットフォームの紹介でした。 JuliaのAPIもあるんですね、すごい！

LightTableを使いましょう - @QuantixResearch

LightTable + Jewel でJuliaを開発しようという話でした。 インラインで変数の値が見れたり、補完もかなり強力なようです。あとで使ってみよう！

感想

新しいおもちゃを手にしたときの子供のような、Juliaに対する期待と興奮が垣間見れたように思います。 おぉ！Juliaすごい！いいじゃん！というような熱気がありました。 Twitterでしか見たことない人と実際にお話したり、情報交換をしたりも大変有益でした。 反省点としては、ちょっと初心者向けの話が少なすぎたかなと思います。次回は初心者向けセッションを設けて欲しいとの要望が多かったです。 会場を用意していただいた株式会社ブレインパッドの皆様、本当にありがとうございました。

あと、本家Juliaステッカーをもらえたのがすごい嬉しかった！

本場シカゴからやって来たステッカー貰った(((o(*ﾟ▽ﾟ*)o))) pic.twitter.com/1AJeQCp0ml

— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 7月 5

HaskellのFusionがあれば速度と抽象化を両立できる

データの分析をする際には配列やベクトルは欠かせないデータ構造です。 大体どの言語にも大体配列は用意されていて、そこにサンプルのデータ等を入れて統計量を計算したり関数に渡して回帰をしたりするわけです。 ベクトルという単位でデータの塊を扱うものの、実際のアルゴリズムではベクトル内の要素一つひとつを見ていって何か処理をしたり計算をすることが多いでしょう。 その際、命令型言語ではforループを陽に使って要素にアクセスすることになります。

簡単な例を見てみましょう。ベクトルv1とv2の内積を求める関数をC言語で書くと以下のようになります。

function dot(double* v1, double* v2, size_t n)
{
    double s = 0.0;
    for (int i = 0; i < n; i++) {
        s += v1[i] * v2[i];
    }
    return s;
}

特に難しいところはありませんね:)

ベクトル化された記法

しかし、大きなベクトルに対してこうして陽にループを書くことは、PythonやRではある種の「悪手」であると認識されています。 理由は単純で、PythonやRのインタープリタのループが遅すぎるため、NumPyのようなCやFortranなどで書かれた高速なライブラリに配列(実際にはPythonのオブジェクトが持つメモリ上のバッファへのポインタ)を渡して計算させるのが定石になっています*1。

NumPyやSciPyではこうしたベクトルに対する関数や演算子が大量に用意されており、これらを組み合わせれば、明示的なforループを避けて簡潔に高速な計算ができるようになっているわけです。こうした操作をベクトル化(vectorized)された操作と言われています。 例えば、先程と同様２つのベクトルの内積を計算するならNumPyにnumpy.dot関数が予め用意されています。

import bumpy as np

v1 = np.array([1.0, 2.0, 3.0])
v2 = np.array([2.0, 3.0, 4.0])

print(np.dot(v1, v2))

ベクトルの要素ごとの積(*)と総和を求める関数を使えば、内積計算を以下のようにも書くことができます。

def dot(v1, v2):
    return np.sum(v1 * v2)

分散の計算

こうしたベクトル化された関数を組み合わせて標本分散(不偏分散)を計算してみましょう。 数式での定義は以下のとおりです。

これを単純にNumPyを使ったPythonのコードに写してみましょう。(もちろん、NumPyには分散を計算するnp.var関数が用意されています！)

def var(x):
    return np.sum((x - np.mean(x))**2) / (np.size(x) - 1)

ベクトル化された関数を使えば、定義式をほとんどそのまま移すだけで、分散を計算する関数が作れました。

しかし、この方法には問題があります。それは、計算の途中に出てくる式の計算のために不要な一時配列が確保され、余計なメモリを食っていることです。 上の例では、x - np.mean(x)では、元の配列xと同じ大きさの配列が確保され、平均からの差分をそこに格納していく形になっています。(x - np.mean(x))**2ではさらに二乗の結果を格納する配列が確保されています。これらの一時的な配列は、本来分散を計算するためにはまったく不要なものです。 内積の例でも、np.sum(v1 * v2)のところでnp.sumに渡す一時ベクトルが新たに作られてしまっています。

def var(x):
    return np.sum((x - np.mean(x))**2) / (np.size(x) - 1)
#                  ~~~~~~~~~~~~~~
#                  ^
#                 ~~~~~~~~~~~~~~~~~~~
#                 ^
#                 余計なベクトルが2個作られている

Pythonのforループは遅くて使えないため、これを回避するためには、CやCythonで陽にループを書くハメになります。

#include <stdlib.h>
#include <stdio.h>

// Two-pass algorithm:
//   http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass_algorithm
double var(double* v, size_t n)
{
    double mean = 0.0;

    for (int i = 0; i < n; i++) {
        mean += v[i];
    }

    mean /= n;

    double ssd = 0.0;

    for (int i = 0; i < n; i++) {
        ssd += (v[i] - mean) * (v[i] - mean);
    }

    return ssd / (n - 1);
}

先ほどのベクトル化された計算と比べると、ループに展開したコードはかなり読みにくく、定義式との対応がわかりづらくなってしまいました。

ちなみに、Juliaでは、マクロを使ってベクトル化された記法からループへと展開するライブラリDevectorize.jlがあります。

fusionを使って宣言的に書き、効率的なコードを吐く

宣言的なベクトル化された記法をしつつ、ループと同等の効率的な実行を可能にするのがfusionです。 fusionでは、中間にできてしまう一時的なベクトルをライブラリによる自動的なコードの書き換え規則とコンパイラによる最適化で排除し、効率的なコードを生成する技術です。 Haskellのvectorパッケージでは、ほとんど意識する必要なくベクトル計算のfusionを行ってくれます。 unboxedな値を格納するベクトルを使うなら、Data.Vector.Unboxedモジュールをインポートして使います。ここでは見た目を簡潔にするため、Preludeの名前が衝突する関数を隠しています。

import Prelude hiding　(sum, zipWith)
import Data.Vector.Unboxed

dot :: Vector Double -> Vector Double -> Double
dot v1 v2 = sum $ zipWith (*) v1 v2

v1 * v2ではなくzipWithを使ってる点が若干分かりにくいかもしれませんが、ループを書かずに宣言的に内積が定義できています。それでいて、実はこのdot関数は計算の中間データを保持するベクトルを作っていません!

fusionがどうやって動いているのか

実際にfusionがどのような仕組みになっているか簡単な例で見てみましょう。 ここでは、Stream Fusion - From Lists to Streams to Nothing at Allのサンプルの一部ををGHCで試せるように書いてみます。 これは、stream fusionというfusionの一例です。 vectorパッケージではVector型に対してfusionを行っていますが、ここではリストを使っています。

minimal stream fusion

例えば、map f . map gというコードは関数合成の結合則とstream . unstream ≡ idという関係式から、

map f . map g ≡ (unstream . mapS f . stream) . (unstream . mapS g . stream)
              ≡ unstream . mapS f . (stream . unstream) . mapS g . stream
              ≡ ustream . mapS f . mapS g . stream

となります。ここで、streamとunstreamが打ち消しあうことで、リストとStreamの変換が消えますが、これはGHCのrewrite ruleによって実現しています。 さらに、GHCの最適化によりStepの生成も排除され、効率的なコードが吐かれるわけです。 リストに対するstream fusionのより詳しい実装は論文の著者等によるstream-fusionパッケージを参照してください。ただ、最近のGHCは既にリストののfusionをサポートしているため、上のコードやstream-fusionパッケージを使ったものは普通に書いたものと同等のパフォーマンスになりました*2。

分散の計算、再び

分散の計算をvectorパッケージを使って書いてみましょう。 fusionの効果を見るため、わざと関数を分けて定義しています。

import Prelude hiding (sum, length, zipWith, map)
import System.Environment (getArgs)
import Data.Vector.Unboxed as V
 
-- | pairwise product
(.*) :: Vector Double -> Vector Double -> Vector Double
x .* y = zipWith (*) x y
 
-- | sample mean
mean :: Vector Double -> Double
mean v = sum v / fromIntegral (length v)
 
-- | deviations from the mean
deviation :: Vector Double -> Vector Double
deviation v = map (\x -> x - mean v) v
 
-- | vectorized square
square :: Vector Double -> Vector Double
square v = v .* v
 
-- | sum of squared deviations from the mean
ssd :: Vector Double -> Double
ssd = sum . square . deviation
 
-- | unbiased sample variance
var :: Vector Double -> Double
var v = ssd v / fromIntegral (length v - 1)

このようにかなり宣言的にプログラムを書くことができ、それぞれの関数は単独で使うことも、他の場所で部品として使うこともできるなど再利用性も極めて高いです。Cの方だと可能なのはせいぜい平均値を計算するmean関数をくくりだすくらいでしょうか。

パフォーマンス比較

大きい配列に対してパフォーマンスを比較しみましょう。 長さ100,000,000の配列に対してCのループを使った実装とパフォーマンスを比べてみます。 使用したコードは以下のリンクにあります。 Gist - variance

Haskell: 0.49s user 0.35s system 99% cpu 0.839 total
C:       0.48s user 0.33s system 98% cpu 0.815 total

Haskellの実装とCの実装とでほとんど同じ速度がでていますね。Haskellのコンパイルには-fllvmでLLVMのコードを吐くようにすると何割か速くなりました。

本当に一時データが割り当てられていないかも見てみましょう。 Doubleは8byteなので、長さ100,000,000の入力ベクトルに対して800,000,000バイト程度が必要になります。 以下のヒープに割り当てられた領域の大きさからも、不要な配列が割り当てられていないのが分かります。

% ./variance-hs $(cat n) +RTS -s -RTS

     800,121,472 bytes allocated in the heap
          10,104 bytes copied during GC
          44,312 bytes maximum residency (1 sample(s))
          21,224 bytes maximum slop
             764 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0         1 colls,     0 par    0.00s    0.00s     0.0001s    0.0001s
  Gen  1         1 colls,     0 par    0.00s    0.06s     0.0562s    0.0562s

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.50s  (  0.72s elapsed)
  GC      time    0.00s  (  0.06s elapsed)
  EXIT    time    0.00s  (  0.06s elapsed)
  Total   time    0.50s  (  0.83s elapsed)

  %GC     time       0.0%  (6.8% elapsed)

  Alloc rate    1,610,816,342 bytes per MUT second

  Productivity 100.0% of total user, 60.0% of total elapsed

新しいfusion

このようにfusionは強力なツールですが、今までのstream fusionにはいくつか限界があります。 最後に、こうした制約を取っ払うために発案された2つのfusionを紹介しましょう。 実装の詳細に関しては、私の能力の限界を超えるため、興味のある方は元の論文を参照してください。

Generalized Stream Fusion

2つのベクトルを結合して新しいベクトルを作りたいことはよくあります。 しかし、既存のstream fusionでは一度に1つづつしか値を取り出せないため、ベクトルを塊として扱うmemcpyのような効率的なコピーができません。 また、同じ理由で複数の値に一気に操作をするSIMD命令も使うことができません。

しかしGeneralized Stream Fusionでは素朴なStreamをより柔軟なBundleという単位にまとめることでこれらの問題を解決してくれます。 2つのベクトルの内積を計算するベンチマークでは、SSEを使うようにしたCのコードと同等かそれ以上の速度を出しています*3。

Generalized Stream Fusionでは、既存のHaskellのコードをほとんど変更することなく、合成可能な再利用性を保ったまま高速化できるのが魅力です。 この機能はvectorパッケージ上に実装されています。 現在OpenBLASやEigen3などの数値計算ライブラリと同等のパフォーマンスは出せていませんが、開発が進めばこれらのライブラリに迫るものができるかもしれません。

Data Flow Fusion

もう一つのstream fusionの拡張が、Data Flow Fusionです。 元のstream fusionではベクトルの値を変換して別のベクトルなり集計値を計算する「消費者」は常にひとりだけに限られています。 しかし、1つのベクトルから複数の統計量を計算したりする分岐をしたいことが常です。 大きなベクトルを複数回ループすることはキャッシュの有効利用ができないため、1つのループにまとめるのに比べて不利になってしまいます。

以下のコードは入力ベクトル(vec1)の要素をすべてインクリメントし(vec2)、正数のみを集め(vec3)たものとその中の最大の数をタプルにして返す関数です。

filterMax :: Vector Int -> (Vector Int, Int)
filterMax vec1 =
  let vec2 = map    (+ 1) vec1
      vec3 = filter (> 0) vec2
      n    = foldl max 0 vec3
  in (vec3, n)

この処理は一度のループで実行可能なはずですが、stream fusionでは1つのfusionにできない処理のようです。また、複数のベクトルに対して同時にループを行うときにループカウンタが重複してしまいレジスタを無駄遣いしてしまうこともパフォーマンスを悪化させてしまいます。

こうした問題をベクトル長の変化やデータフローを解析する高度なfusionによって解決するのがData Flow Fusionです*4。 実装は、Repaの配列に対するGHCの最適化プラグイン(repa-plugin)として提供されているようです。 この最適化により通常のstream fusion(Stream)に比べてdata flow fusion(Flow)では大きく性能が向上し、人間が手でfusionしたCのコード(Hand-fused C)に迫るパフォーマンスの向上を実現しています。

fusionは既に実用的な技術になっており、より強力になったfusionが使えるようになるのもそう遠くないかもしれません。

*1:Pythonで数値計算のコツ：for文書いたら負けかなと思っている

*2:GHC User’s Guide - 7.21.4. List fusion

*3:Exploiting Vector Instructions with Generalized Stream Fusion

*4:Data Flow Fusion with Series Expressions in Haskell

Pythonのコマンドライン引数処理の決定版 docopt (あとJuliaに移植したよ)

Pythonをよく使う人にはよく知ってる人も多いのですが、docoptという便利ライブラリがあります。 docoptはargparseやoptparseのようなコマンドライン引数をパースするライブラリなのですが、その発想がコロンブスの卵なのです。

例えばPython標準のargparseだと、argparseのAPIを組み合わせてパーサを組み立てるわけです。するとパーサと一緒にヘルプも作ってくれて、"program --help"などとすると自動生成されたヘルプを表示してくれるようになります。 しかし、そのAPIを覚えるのが大変で、毎回ドキュメントを読まないと忘れちゃうわけです。

import argparse

parser = argparse.ArgumentParser(description='Process some integers.')
parser.add_argument('integers', metavar='N', type=int, nargs='+',
                   help='an integer for the accumulator')
parser.add_argument('--sum', dest='accumulate', action='store_const',
                   const=sum, default=max,
                   help='sum the integers (default: find the max)')

args = parser.parse_args()
print(args.accumulate(args.integers))

https://docs.python.org/dev/library/argparse.html

その問題をdocoptではスマートに解決してくれます。docoptでは、パーサからヘルプを生成するのではなくて、ヘルプからパーサを生成します。そしてコマンドライン引数は辞書(dict)を使って取得できるようになります。 docoptモジュールがエクスポートしているのはdocopt関数のみで、複雑なAPIを覚える必要もなく、馴染みのあるヘルプを文法に従って書けばパーサを生成してくれます。

先の例をdocoptで書き換えれば、以下のようになります。

"""Process some integers.

usage: prog.py [-h] [--sum] N [N...]

options:
    -h, --help  show this help message and exit
    --sum        sum the integers (default: find the max)
"""

from docopt import docopt

if __name__ == '__main__':
    args = docopt(__doc__)
    Ns = map(int, args["N"])
    if args["--sum"]:
        print(sum(Ns))
    else:
        print(max(Ns))

とても直感的ではありませんか！ 使用すると、以下のようになります。 引数が文法に合わなければ usage が表示されます。

/Users/kenta/tmp% python3 prog.py 4
4
/Users/kenta/tmp% python3 prog.py 4 3 5 2
5
/Users/kenta/tmp% python3 prog.py --sum 4 3 5 2
14
/Users/kenta/tmp% python3 prog.py
Usage: prog.py [-h] [--sum] N [N...]
/Users/kenta/tmp% python3 prog.py -h
Process some integers.

usage: prog.py [-h] [--sum] N [N...]

options:
    -h, --help  show this help message and exit
    --sum        sum the integers (default: find the max)

ヘルプの文法は標準的なもので、以下の様なものがあります。

要素	説明	補足
<argument>	位置引数	<x> <y>に引数1 2を与えれば、{"<x>": 1, "<y>": 2}と結び付けられる
--option, -o	オプション	ハイフンで始まる引数
(arguments)	必須の引数	引数はデフォルトで必須だが、グループ化するときに使える
[arguments]	必須でない引数	(...)と違って、なくてもエラーにはならない
arg1 | arg2	arg1もしくはarg2	3つ以上連ねることもできる (例. --up | --down | --left | --right)
args...	可変数引数	返り値はリストになる

これらを組み合わせて、以下のように複雑なものも書けます。これをargparseで書くのは至難の業でしょう。

Naval Fate.

Usage:
  naval_fate ship new <name>...
  naval_fate ship <name> move <x> <y> [--speed=<kn>]
  naval_fate ship shoot <x> <y>
  naval_fate mine (set|remove) <x> <y> [--moored|--drifting]
  naval_fate -h | --help
  naval_fate --version

Options:
  -h --help     Show this screen.
  --version     Show version.
  --speed=<kn>  Speed in knots [default: 10].
  --moored      Moored (anchored) mine.
  --drifting    Drifting mine.

(from docopt.org)

docoptを使うモチベーションやより詳しい文法については公式サイト(http://docopt.org/)や作者によるPyCon UK 2012の動画を御覧ください。

PyCon UK 2012: Create beautiful command-line ...

言語非依存なのも利点の1つです。一度使い方を覚えれば、他言語でもAPIを特別覚える必要がなく使えます。 オリジナルのPythonの他にも、Haskell, Go, C, Ruby, CoffeeScriptなどにも移植されています。

Julia版のdocopt - DocOpt.jl

個人的に最近Juliaをよく使っているのですが、docoptが無かったので移植したしました。 GitHubに上げたのは何ヶ月か前なのですが、最近やっと正式なdocoptファミリーの一員になり、Juliaの正式なパッケージにしたので改めて紹介します。

ソースコードはこちらです。https://github.com/docopt/DocOpt.jl

文法や使い方はPython版と同じですので、既にPython版を使ったことのある方は何も覚えることはありません。 using DocOptしてdocopt関数をインポートし、ヘルプを渡すだけです。

doc = """Naval Fate.

Usage:
  naval_fate.py ship new <name>...
  naval_fate.py ship <name> move <x> <y> [--speed=<kn>]
  naval_fate.py ship shoot <x> <y>
  naval_fate.py mine (set|remove) <x> <y> [--moored|--drifting]
  naval_fate.py -h | --help
  naval_fate.py --version

Options:
  -h --help     Show this screen.
  --version     Show version.
  --speed=<kn>  Speed in knots [default: 10].
  --moored      Moored (anchored) mine.
  --drifting    Drifting mine.

"""

using DocOpt  # import docopt function

arguments = docopt(doc, version=v"2.0.0")
dump(arguments)

パース結果はPython同様、辞書になります。

インストールはJuliaの対話環境から、Pkg.add("DocOpt")とすればインストールできます。 残念ながら今のところJuliaのv0.2系をサポートしていないので、v0.3のコンパイル済みprerelease版をインストールするか、ソースからビルドすることになります。 v0.2で使いたいという要望があれば、v0.2もサポートしようと思います。

Haskellでデータ処理がしたい！

この記事はIIJの@kazu_yamamotoさんの下で、アルバイトとして書いています。 最終的な目的はHaskellでの統計処理やデータ処理のツールを整備しようというものですが、その下調べをした内容をこの記事にまとめています。 ご意見などコメントいただけるとありがたいです。

---

今から3年前の2011年のICFPに、Emily G. Mitchell 氏のHaskellを使ったシミュレーションに関するExperience Reportが掲載されました。

[PDF] Functional Programming through Deep Time - Modeling the first complex ecosystems on Earth

このレポートでは筆者が行ったエディアカラ紀の生態系シミュレーションの実装に関してプログラミング言語の選択や、使った言語にどのような利点・欠点があったかが述べられています。 初めはRを用いていたようですが、実行速度とシミュレーションのコード修正をするのに問題があってHaskellを選択した経緯とそのときに感じたHaskellでデータ処理をする際の問題点はとても興味深いものでした。 筆者は元々Haskellの使い手だったわけではないようですが、夫があのNeil Michell氏であることもあってか、シュミレーションにはHaskellを使い、プロットなどを含むデータ処理にはRを使ったようです。

その中でRと比較したHaskellの利点として以下のようなことが挙げられています。

• 構文が一貫していること
• 純粋さと静的型付けのおかげでコードのリファクタリングが容易であったこと
• 型システムのおかげで、デバッグが容易であったこと

しかし、Haskellまわりのツールには不満点もあるようで、

• 異なるパッケージの似たようなデータ型で互いにデータを変換するのが困難であること
• Haskellパッケージのインストールが難しいこと (RはGUIでパッケージをインストールすることもできる)
• 標準ライブラリで統計量の計算やプロットができないこと

などが挙げられています。

現状Haskellの統計パッケージあたりがどういうのがあるか調べてみると、Haskell Wikiにあるのは、

• hstats
• hmatrix-gsl-stats
• hstatistics

の3つだけです。(hstatsは2008年で止まっている)

プロットに関しては、

• Chart
• plot
• easyplot

といった感じです。正直、充実しているとは言い難い状況ですね。

現状の機能では大変なことを実感するには、簡単な具体例を考えてみれば分かりやすいでしょう。以下の例は、2012年のHaskell cafeにあったスレッドの一部を翻訳したものです。

練習として、GHCiを使って適当な浮動小数点数の表になっているCSVファイルを読み込んで、ある列と別の列の間で線形回帰をして、回帰直線と一緒に散布図をプロットしてみてください。残差の正規性もチェックしてみるといいかも。大きなデータも扱えるようにするには、bytestringとかattoparsecなんかも使う必要があるでしょう。x86_64のGHCiからhmatrix/gslの関数を使うと、セグメンテーションフォルトとかバスエラーになる既知のバグがあることにも気をつけましょう*1。 何か当たり前なことを見落としてるかもしれませんが、(データの)コンテナ・永続化・パース・統計・プロットをするのにどのパッケージを選ぶのにすっっごく時間がかかるでしょう。

[Haskell-cafe] Mathematics and Statistics libraries

確かにこれをHaskellでやるのはちょっと考えてみただけどもやりたくない感じですね。 まずCSVのパースなのでcassavaあたりを使うといいのでしょうか。cassavaのData.CSVモジュールを見るとdecodeというのがあるのでこれを使えば良さそうですね。でも引数がByteString型なのでbytestringパッケージのData.ByteStringモジュールをインポートしてファイルを読む必要がありますね。返り値はVector型なのでそれから必要な列を取り出して線形回帰すればいいのですが、Haskellで線形回帰をするライブラリを探してぐぐってみるとstatistics-linregがありますね。よし、これをcabalでインストールしてとかやってられません。 しかしRやPythonなら、慣れている人なら5分とかからずにやってしまうでしょう。 試しにRでやってみるとたった6行です。

library(ggplot2)  # ggplot2ライブラリを読み込む
data = read.csv("data.csv")  # ヘッダーつきCSVファイルを読み込む
data.lm = lm(y ~ x, data)  # 線形回帰 
plot(data.lm)  # 残差の正規性をQ-Qプロットでチェック
ggplot(data, aes(x=x, y=y)) + geom_point() + stat_smooth(method="lm")  # 試しにプロット
ggsave("data.png")  # プロットを画像ファイルにセーブ

入力のCSVファイル:

id,x,y
0,0.125,2.140
1,0.065.567
2,-0.079,3.115
3,-0.193,1.765
...

プロット:

この簡潔さを可能にしているRの特徴として挙げられるのは、

• 表データを扱う組み込みのデータ構造 (データフレーム)
• データ処理のための組み込み関数 (read.csv(), lm() など)
• プロットのためのライブラリ (ggplot2, latticeなど)

などでしょう。 Pythonでも組み込みの機能はRほどではないにしろ、データフレームにはpandas、データ処理にはnumpy/scipy/scikit-learn、プロットにはmatplotlibがあり、IPythonを使えばRと同様インタラクティブに処理可能になります。

これを真似て処理するデータを収めるコンテナを統一的に扱えるようにし、データ処理によく使う関数をひとつのモジュールをインポートするだけで使えるようにすれば、結構Haskellでもいけるようになるのではないでしょうか。 加えて、快適に使うために欲しいものとしては、ドキュメントがサクサク読めて、データを眺めたりプロットをはめ込めるRStudioやIPython notebookのような統合的対話環境があると嬉しいです。

これらの機能を整備したとして、Haskellを使う利点は何でしょうか。 先の論文の要旨も踏まえて私の考えを3点述べると、

1. 静的な型付けによる安全性とメンテナンス性の向上
2. データ並列による処理の高速化
3. 入出力の抽象化とストリーミング処理

これらはRやPythonではそう簡単には実現できない特徴ですので、うまくデータ処理のニーズと合致すればとても魅力的なツールになるはずです。 最後に、データ解析者やHaskellerの方々の目線から、このような機能は必須だとか、Haskellのこのライブラリを使うといいなどの意見をコメントいただけるととてもありがたいです。

*1:GHCの7.8.xでは修正されているようです。 https://twitter.com/kazu_yamamoto/status/454741146053255168