2012年2月11日土曜日

Rのちょっと速いコードの書き方

このエントリーをはてなブックマークに追加
Pocket

Rはループ速度が段違いに遅いと言われる。確かにループとメソッド呼び出しで構成したマイクロベンチマークを実行すると、Javaが6.32秒、C++で6.33秒で終わる処理が、87時間18分16.0秒(推定値)かかったりする。S-PLUSやMatlabなどの他の同種の言語よりは高速か同等と指摘されているが、汎用言語に比べると断然遅い(Benchmark 2)。

もちろん大半の計算は問題ない。標本数1万ぐらいのサンプルでプロビット分析を行っても、1秒もかからず計算が終わる。コマンドを打っている時間の方が圧倒的に長い。しかし人間はどのような環境でも速度に憧れるものだ。そしてRでも短時間に処理を終わらせる為のコツはある。

1. パッケージや内部関数を使う
大抵の著名パッケージはC言語で実装されているので、内部的な処理は高速だ。Rではなるべくコードを書かない方が良い。スクリプト言語やインタープリッタ全般に言える特性で、PHP等にも言える特性だ。関数やパッケージを覚えていく事が、まず大事な事になる。
2. 行列演算≿apply関数≿ループ演算
Rと言えばベクトルや行列演算。そこだけ実行速度が速い。ある対数尤度関数をループ処理で記述すると1.11秒かかったが、apply()関数でまとめると0.48秒に、配列演算にまとめると0.03秒になった(Rは行列演算≿apply関数≿ループ演算)。線形代数の問題集を買ってきてトレーニングすると、良いRプログラマになれるかも知れない。
対数尤度関数の記述方法と処理時間
記述方法 経過時間 行列演算比
行列演算 0.03秒 1倍
apply関数 0.48秒 16倍
ループ演算 1.11秒 37倍
一般的にはベクトルにしないとループ演算が増えるので、誤解される事は無いであろうが、ベクトルや行列演算を使うと言う事は、X[i]のような添字操作を行わず、sumやmean等のベクトル操作を行う関数を用いる事だ。

追記(2012/02/15 02:00):別所でifelse(x>=0,y,z)とするよりは、yとzが必ず値を返すのであれば(x>=0)*y+(x<0)*zの方が速いとツッコミを受けたが、これもベクター演算化と言う事になる(ifelse()と論理値の掛け算)。

3. 馬鹿とリストは使いよう
ベクターと比較してリストは遅いと言う指摘がある。なるほど、同じ処理をすれば遅い。しかしベクターで出来ることをリストにしようと言う人はいないであろう。データの紐付けなどをRで行う場合は、subset関数でデータフレームの検索をいちいち行うより、ハッシュ的にリストを使う方が高速になる(関連ページ:Rでミクロデータにマクロデータを紐付ける)。
4. 同じ計算は二度しない
x/2 + y/2 + z/2と(x + y + z)/2を比較すると、コンパイルを行う言語とは違いRは後者の方が速い(加算後の除算≿除算後の加算)。実行時に最適化が行われないからだ。速度を出したい場合は、演算回数を減らす工夫も有効になる。ただし数値演算時にはオーバーフローによる誤差を避けるためにテイラー展開を行い、速度を犠牲にしているケースもある。無闇な最適化は誤計算の元だ。
5. 配列やベクターのサイズは固定する
Rではx <- numeric(0)とベクターを初期化していても、x[2] <- 2123などと代入が出来る。自動的にベクターが増加するためだ。しかし、配列やベクター・サイズの変更はコストがかかる。同じ50,000のサンプル・サイズでも、1回の初期化で済ませれば0.16秒で済む処理が、50,000回のサイズ変更を伴う処理にすると3.46秒かかったりするので、大量のデータ処理を行う場合は注意をした方が良いかも知れない(配列の拡張とパフォーマンス)。
6. 高速な数値演算アルゴリズムを選択する
ある1,000サンプルのProbit Modelで、ほぼ同一の結果を返す計算方法を何種類も試してみた。大雑把に言って、一般化線形モデル(GLM)が最も速く、最尤法が続き、MCMC(M-Hアルゴリズム)が遅い。MCMCは分布も計算するので計算結果が異なるわけだが、パラメーターの推定量を出すだけならばGLMが望ましい事になる。
数値演算アルゴリズムと処理時間
アルゴリズム 関数名 経過時間
GLM glm() 0.02秒
Newton-Raphson nlm() 0.03秒
BFGS optim() 0.04秒
L-BFGS-B optim() 0.04秒
Nelder-Mead optim() 0.06秒
CG optim() 0.45秒
SANN optim() 4.70秒
MCMC(Metropolis-Hastings) MCMCmetrop1R 9.61秒
7. whileよりもfor
ループ構造を作るときに、n回処理を繰り返すような大抵のケースではwhileよりもforを使う方が高速になる。古いCプログラマだとRのforはメモリーを無駄に使っている気がして避けてしまうのだが、Rでは連続配列の初期化が速く、演算代入処理が遅いために何割かの差がつくようだ。最も単純な演算で30%強の経過時間の削減ができた。
8. 高速化パッケージを用いる
ビルトインされたcompilerパッケージにあるcompile/cmpfun関数でバイトコードを生成したり、Snow/ffパッケージで複数プロセッサによる処理を行ったりできる(Rの高速化)。適用は容易であるが、ループをベクトル演算に展開したり、外部ネイティブ・コードを呼び出す方が効果的なようだ。FizzBuzz問題で素朴なRのコード、最適化されたRのコード、Cで書かれた外部コードを比較すると、11.08秒、0.58秒、0.40秒であるが、素朴なRのコードにcmpfun関数を適用しても4.72秒にしかならなかった。また最適化されたRのコードはほとんど変化が無い(Rで解くFizzBuzz問題Rで解いていないFizzBuzz問題)。
9. 処理時間を計測する
system.time({ ... })で処理時間を計測できるので、少数サンプルで試行錯誤を行う事で、効率の良いコードを書く事ができる。

MCMCで何時間もかかる演算をする人は、引数として与える確率分布関数を工夫する価値はあるようだが、特にテクニカルな事をしない人は余り速度は気にしなくて良いようだ。

究極的にはC言語で書いたライブラリを呼び出すしか無いのであろうが、その場合は限りなくC言語になるため、C++/Boostでも使う方が手っ取り早くなる。筆者のようなライト・ユーザーには少々荷が重く、Rでちょっと高速なコードを書いて自己満足を得るのも悪く無いと思う。

0 コメント:

コメントを投稿