＊ここに掲載するのは、エジンバラ大学のMichal Piotrowski博士によるHECToRレポート「Optimisation of R bootstrapping on HECToR with SPRINT, Michal Piotrowski, 08/05/2012」を要約したものです。
概要
このプロジェクトでは、サンプル推定の正確性の計算と統計的推論に広く用いられ、計算的に負荷の大きなリサンプリング手法であるブートストラップ法の並列版の実装を行います。SPRINTプロジェクトにおいてユーザ要求に関するサーベイが行われた結果、並列化対象としてブートストラップ関数が選ばれました。
The Simple Parallel R Interface (SPRINT) [1]は、エジンバラ大学のDivision of Pathway MedicineおよびEPCCの共同プロジェクトとして開始され、特別な知識を必要とせずにRユーザがHECToRスーパーコンピューティング資源や他のHPCシステムを利用できるようにすることを目的としています。Rはフリーの統計学計算のための言語および計算環境です(http://www.r-project.org/)。
このプロジェクトの結果、SPRINTフレームワークの新しいバージョン1.0がThe Comprehensive R Archive Network (CRAN)にリリースされました。
ブートストラップ法
ブートストラップ法は70年代にEfronにより紹介された方法で、コンピュータの処理能力の発展と深く結びついています。これは、母集団の推定における不確かさの決定に用いられるコンピュータシミュレーションの一般的な技法です。元の観測データの復元抽出によるリサンプリングが行われ、元のデータの分布特性とは異なる同じ大きさの新しいサンプル群が作成されます。このサンプル数を増やすと結果の推定精度が向上しますが、多くの処理能力が要求されます。嘗てはブートストラップ標本数は計算資源により制約されるのが通常でした。ブートストラップ法は幅広く応用が可能な良く知られた方法で、計算集約的かつSPRINTに適した並列化のし易い方法です。
ブートストラップ法の動作を理解するために例題を考えます。復元抽出により未知の母集団Pから平均Mを推定することを考えます。サンプル平均mを計算し、推定における不確かさの度合いを知るために標準誤差も推定したいのですが、これは未知の母集団Pの分散に依存します。ブートストラップ法が基にするアイデアは、与えられたサンプルそのものだけを用いて母集団全体の分布がシミュレーション可能であるというものです。ここで、サンプルと母集団の分布は似ており、復元抽出によるリサンプリングで新しいサンプルを生成できることを仮定しています。これは測定に統計的ばらつきを導入するものです。こうして、ブートストラップ法により得られた多数のm値の標準偏差を計算することで、Pの分散を推定することが可能になります。
Rでは、ブートストラップ標本を生成する`boot`関数がbootパッケージに実装されています。これはRで表現される如何なる関数からもブートストラップを生成可能です。また、こうしたサンプルからバイアスとブートストラップの信頼区間の推定が可能です。`boot`コマンドは、データセットのリサンプリングやこれらサンプルの統計計算を実行します。
実装
最初の並列ブートストラップ法である`pboot`はEPCCのMScプロジェクトで実装されていました。この初期バージョンの性能テスト結果は良好でしたが、そこ際用いた古いSPRINTでは、同時に全てのプロセスがRインタープリターへアクセスできない状態でした。結果として、その他のSPRINT関数はこれと適合しなく、並列処理内での利用が不可能でした。さらにデータ分散が不十分だったため、HECToRフェーズ3のマシン上では8プロセッサまでしかスケールしていませんでした。また既存のブートストラップ機能の一部が実装されたのみで、ノンパラメトリックケースのみに限定されていました。
SPRINTの新しいバージョンは、全てのノードに個別のRプロセスとそのランタイム環境を用意することが出来ます。これは、R又は部分的にRで記述されたより複雑な関数の並列化を大幅に簡略化します。関数はC言語で書き換える必要はなく、シリアライズ後に全ノードへ転送されます。その後デシリアライズ化されてRインタープリターへ返却されて計算が実行されます。並列化はデータ分散で行われ、各ノードでは同じ関数がデータセットの一部を処理します。この手法はデータ依存性がない場合にのみ実行可能です。
このモデルは、より複雑なRオブジェクトをノード間で送信することや、ユーザ定義関数を引数として渡すことを可能にします。これはブートストラップ法実装に極めて適しており、多くの異なるサンプルデータを用いて任意の関数のコピーを多くのプロセスで実行することが可能です。ここではこれら関数の実行にデータ依存性はありません。さらにインターフェイスについてはシリアル版と同じで良く、全引数を含めた既存の関数呼出しをシリアライズして、全ノード上のリサンプリング数を用いて実行することになります。
統計関数はブートストラップ複製数のループで計算されます。その前にブートストラップ複製インデックスをシリアルに生成しています。これは性能に影響を与えますが、簡潔なコード作成とシリアル版との再現性確保のために残しました。
各繰り返しの最後に、一つのブートストラップの統計量が各プロセスで計算されローカルリストに追加されます。全計算が完了した後に、このリストはSPRINTの並列リダクション関数で結合されます。これは当初、並列ランダムフォレスト実装においてツリーリダクションを高速化するために実装されたものです。この結合関数は並列リダクション関数へ引数として渡されます。このやり方を用いれば、並列結合関数は線形から対数的な複雑性へ減少されて性能を大きく改善します。
並列ブートストラップインターフェイスはシリアル実行の場合と同じです。必要な変更は、SPRINTライブラリのロードと、boot関数からpbootへの変更、およびスクリプトの最後にMPI環境を終了する、という最小限で済みます。
結果
2816ノード、各ノードに2つの16コアのAMD Opteron 2.3Ghz Interlagosプロセッサ、総コア数が90,112コアの現状のフェーズ3システムを用いました。ベンチマーク全てでGolub et al.のデータを用いました。これにはトレーニングサンプル群とテストサンプルが含まれます。これは、急性リンパ性白血病（ALL）の47患者と急性骨髄性白血病（AML）の25患者、および7129遺伝子発現数から成ります。3種の統計関数に対してブートストラップ法を適用し、関数の複雑性と複製数に対するスケール性の依存性を調査しました。
比較的に複雑性の小さな関数例として中央偏差と標準偏差を選びました。両者の実行時条件の違いはブートストラップ複製数です。中央偏差では24999、標準偏差では9999個の複製を用いました。シリアル実行される複製インデックス生成の影響により、4プロセスおよび16プロセス以降に性能は大きく劣化します。
最後のベンチマークはmedoids関数です。ここで、サンプルサイズを2500遺伝子、複製数を999個にして、6 medoids のクラスタリングを行いました。PAM計算の複雑性は以前のものに比べて極めて高いため、スケール性を持つプロセッサー数は伸びますが、32プロセッサーを超えると劣化します。
謝辞
このプロジェクトは、NAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびNAG Ltd.のCSEサポートサービスにより管理運営されています。
文献
| [1] | J. Hill, M. Hambley, T. Forster, M. Mewissen, T. Sloan, F. Scharinger, A. Trew, and P. Ghazal. SPRINT: A new parallel framework for R. BMC Bioinformatics, 9(1):558, 2008 |