プログラムの高速化・並列化サービスの事例

＊ここに掲載するのは、エジンバラ大学EPCCのSavvas Petrou博士によるHECToRレポート「SPRINT D3.2 Final project report, Savvas Petrou」を要約したものです。

概要

遺伝情報の解析には膨大な計算演算能力とメモリーが必要になります。この数年間に、生物研究分野においてハイスループットおよび高並列の実験技術が広く導入されました。マイクロアレイ技術はその際立った例の一つで、数千から数百万の遺伝子や数10から数千の異なるサンプルの遺伝子配列の同時測定を可能にしました。こうした研究は前例のない量のデータを生成し、既存のバイオインフォマティクス計算設備が試されています。新たな全ゲノム関連研究や臨床プロジェクトは、数百から数千のマイクロアレイ実験を必要とします。

GPX-MEA [3]やArrayExpress [4]といったマイクロアレイ実験の蓄積が成長し続けていますが、 その研究にはこれらデータベースに含まれるデータから有用な新規の科学的知識を抽出し、解析するための計算資源が不足しています。このトレンドは続くことが予想されており、その間の技術の進歩は解析対象のデータ量増大に寄与し続けていきます。網羅性の強化は、より多くの遺伝子配列を一つのアレイ上で解析することを可能にします。その結果、さらに大きなデータ量が創出されます。著名なフリーソフトのデータ解析パッケージとしてR[5]が有ります。Rは最近のBBSRC[6]の検証の中で、生物生命科学で利用される重要なコードの一つとして認識されています。

最近になってEPCC[7]はDivision of Pathway Medicine (DPM) [8]と共に、SPRINT[9]と呼ぶRを用いたプロトタイプパッケージを開発しました。これは遺伝子解析に利用される主要な統計機能を並列化したものです[10]。このプロトタイプは、並列統計関数がRとインターフェイスされて、生物研究者がHPC分野に参入できる簡便なやり方を提供可能であることを実証しました。SPRINTの目的は、最小限のHPC知識を持って、Rスクリプトの修正を最小限にしつつ最大の性能を引き出すことです。

SPRINTはCとMPI[11]で記述されているため、HECToR[12]を含めた多くのプラットフォームに対してポータブルです。筆者らはこれまでSPRINTに関する多くの発表を行ってきました[13,14]。

このプロジェクトではSPRINTをHECToRへ移植し、重要なルーチンを最適化し、付加的な機能を追加します。プロジェクトの目的は、SPRINTに実装した2つの並列処理関数の最適化と並列化の結果を示すことに在ります。第一の作業[15]では、既存の並列の相関計算関数のスケール性の調査と、512プロセスまでスケールする新しいバージョンの実装を行いました。第二の作業[16]では、同じく512プロセスまでスケールする並べ替え検定関数の並列版の実装を行いました。

この仕事は、Rユーザミーティング[28]とData-Intensive Research Workshop[29]で口頭発表されました。また、ACM International Symposium on High Performance Distributed Computing (HPDC) [31]のEmerging Computational Methods for the Life Sciences Workshop [30]において論文を発表しました。この概要もuseR! 2010 International Conference [32]で公開されました。

このプロジェクトのソースコードはGNUライセンス化でWebサイトから入手可能です。CRAN [33]リポジトリで公開準備中です。

背景

·R統計言語

Rは統計処理環境の興味深いプログラミング言語です。これはオープンソースのGNUパッケージであり、S言語と似たものです。この数年に渡ってR利用者は増加し、毎月のように多くの新しいライブラリが開発され公開されています。現状のリリース版では、Rは並列処理サポートはされておらず、HPC環境を用いるには不都合です。近年、幾つかのRコミュニティグループは、Rで並列処理を可能にする開発を行っています[17,18,19,20,21] 。SPRINTはそうしたパッケージの内の一つです。

·SPRINTパッケージ

SPRINT(The Simple Parallel R INTerface)はRパッケージとして、並列関数ライブラリおよびRへ並列関数を追加するインターフェイス両方を備えます。SPRINTは2007年に3年間のプロジェクトとして開始され、任意のパッケージがR言語を用いてHPC資源に容易にアクセスが出来るような機能を提供することを目的とされました。 現在のプロジェクトは2009年4月に開始されて、ウェルカム・トラスト[22]によりファンドされた2年間のプロジェクトです。主たる目的はSPRINTの機能拡張です。このプロジェクトはエジンバラ大学の、EPCC[7]およびDivision of Pathway Medicine (DPM) [8]のコラボプロジェクトです。
SPRINTパッケージは、R統計処理言語と、C言語およびMPIへのインターフェイスの3種の部品から成ります。既存パッケージの拡張にはこれらの知識が必要です。

·HECToRシステム

この新しい並列関数のテストとベンチマークには英国スーパーコンピューティングサービスのHECToR Cray XTシステム[12]を用いました。このシステム（フェーズ2a）は1416ブレードで構成され、各ブレードは4つのクアッドコアプロセッサーを搭載します。このCPUはAMD 2.3 GHz Opteronでメモリーは8GBです。全体で、22,656コア、45.3TBメモリーであり、理論ピーク性能は208TFlopsです。
このシステムはCray X2ベクトル演算システムを持ちます。これは28ベクトルユニットで構成され、各々は4つのCrayベクトルプロセッサを搭載し、32GBのメモリーを共有します。全体で、112ベクトルプロセッサで、896GBメモリーです。このプロセッサの理論ピーク性能は25GFlopsであり、全体では2.87TFlopsです。

HECToR上でRパッケージのインストールの詳細はSPRINTのWebサイト[9]に在ります。

並列相関関数：pcor

以前のSPRINTプロジェクト[10]でHillらは、ピアソンのペアワイズ相関関数の並列化を試みました。Edinburgh Compute and Data Facility [23]におけるベンチマークの結果から、4プロセッサを超えるとスケールしないという結果が示されました。このプロジェクトでは最初に、Hillらのpcor実装をHECToR上で同じパラメータを用いてベンチマークを行いました。

ベンチマークは11.001遺伝子の320サンプルを含むデータセットで実行しました。この結果から以下のことが観察されます：
・計算カーネルのスケール性は64ワーカープロセスを超えると劣化する。
・256プロセスを超えるとコードは動作しない。
・出力ステップはシリアル実行のためスケールせず、膨大な時間が掛かる。

また、実装上の特性として以下の点があります：
・全係数配列サイズは入力データセットよりも何倍も大きい。これはマスタープロセスで膨大なメモリーを要求し、扱いが可能な入力データ量を制限している。
・作業割り振り後のマスタープロセスへの結果返却処理では、マスタープロセスの処理能力を超えて多くのリクエストが待ち合わされて、膨大な時間が掛かる。これによりワーカープロセスはブロックされて無駄な時間が浪費される。
・ピアソン相関計算に必要な複数のパラメータは、保存や再利用ではなく再計算されている。

·新しい並列アルゴリズム

既存の出力実装はマスタープロセスにしかありません。ワーカープロセスは割り当てられた作業が終わったら、その結果をマスタープロセスへ送信します。この実装の問題点は以下の通りです：
・マスタープロセスは全ての結果を保存する配列があるための十分大きなメモリーを必要とします。この配列は入力データの列数に等しいランク数を持つ対称な2次元配列です。これにより、処理可能なデータサイズは制限されます。
・計算カーネル間で膨大な量のデータ交換が発生し、処理の応答が非効率となります。

この解決にはMPIの並列I/Oの利用が有効です。出力処理を並列化すればこうした問題を回避することが出来ます。この時には、相関係数はワーカープロセスにローカルに保存され、並列にファイル出力されます。結果はワーカープロセス全体に分散することになるため、マスタープロセスはもはや全結果配列を必要としません。この並列出力はHECToRのLustreファイルシステムの高速性を生かして、出力時間を削減するでしょう。さらにこれは、ファイルシステムの最大I/Oスループットに達するまで、プロセス数の増加に伴いスケールすると考えられます。

しかしながらMPI-I/Oはバイナリー形式のファイルしか扱えません。バイナリー形式であればファイルの読み書きや、ワーカープロセス間で重なりの無いようにファイルを分割して、効率的な速度とバンド幅で集団的に読み書きが可能です。C言語とR環境は共に同じバイナリー形式を取るため、MPI-I/Oを用いてCで書かれたファイルを、Rのバイナリーファイル関数で扱うことが可能です。

1ファイルへの結果の出力処理に関する考察から、膨大な結果を生成する入力配列の処理に新たな実装を考えることが出来ます。ワーカープロセスを増やせば利用できるメモリーは増えるので、結果の保存により多くのメモリーを用いる事が可能になります。しかしながら入力データサイズは事前に決まっているため、ワーカープロセスに保存可能な相関係数用のメモリーサイズは、入力データサイズで決まります。

RのFFパッケージ[24]は、R環境においてメモリーマップトファイルが利用できます。これはバイナリーファイルのopen,read,writeにC言語ライブラリを用います。FFが操作可能なバイナリーファイル構造は、並列実装で生成された構造と同一です。これは、MPI-I/Oで生成された結果ファイルをFFオブジェクトを用いてR環境で簡単に扱うことが出来ることを意味します。最終結果の書き込みが終了して制御がR環境に戻った時点で、FFオブジェクトが生成されてハンドルがユーザに戻ります。このハンドルを用いればユーザは、あたかもそれがメモリー内の通常のR配列であるかのようにファイルの内容を読み書きすることができます。

·結果

並列I/Oによりプロファイルは劇的に変化しました。計算カーネルは512プロセスまでスケールし、出力時間は200倍以上高速化しました。

さらに、入力データサイズを変えて測定しました。この実行は以前はメモリー制限で不可能だったデータサイズです。実行は全て256プロセスで行いました。

この新しい並列実装版は、以前には実行不可能だったデータの計算が可能になり、妥当な時間で計算を実行可能です。

並べ替え検定関数：pmaxの並列化

並べ替え検定は、データセットの統計的有意性を計るために使われる統計分析法の1つです。並列化対象はmulttest[26]パッケージの関数mt.maxT[25]です。これはステップダウン法の多重検定での調整済みp値の計算です[27]。

並べ替え検定の制約に、並べ替えの回数があります。回数が増えると、関数計算時間が増加します。統計的有意性の精度は並べ替え回数に依存します。理想的には正確な値を得るには全並べ替えを実行することが必要です。しかしながら全並べ替え数は小さなデータにおいても膨大になります。こうした大きな回数の場合には、シリアル実行の計算は数日かかります。

·シリアルコード実装

シリアルコードは2種の並べ替え生成器を実装しています。それぞれは、ランダム並べ替えと完全並べ替えです。これら生成器は6種の統計手法で用いられています。さらにユーザは計算前に並べ替え結果をメモリーに保存するか、一度に計算するかの選択が可能です。

完全計算では並べ替え結果はメモリーに保存されません。Block-Fにおいても同様です。

·新しい並列手法

シリアル実行と同じ結果を生成するには、各プロセスで実行される並べ替えの選択には注意が必要です。最初の並べ替えでは最初の行ラベルが用いられ、マスタープロセスで1度だけ考慮が必要です。他の全てのプロセスはこの並べ替えをスキップします。さらに、各プロセスの生成器は適切な並べ替え値に対応しなくてはなりません。初期化関数に追加の変数が渡され、その値によって生成器は数サイクル空回りし、適切な並べ替え値へその生成器を初期化します。

全プロセスは、全データセットのコピーとユーザ指定オプションを持つ必要があります。計算に入る前に、データセットとユーザオプションが全プロセスにブロードキャストされます。さらに計算カーネルの最後では、部分結果がマスタープロセスへ集約されます。これらの結果から最終的なp値が計算されます。

この並列実装でのプロセスは、計算の最初と最後のみ通信します。計算カーネル実行時には各プロセスは依存関係はありません。並列版は前処理と計算カーネルの間に、追加のステップがあります。この新しいステップは、全プロセスにデータをブロードキャストするもので、データ収集と最終結果計算の役割も果たします。

第二の大きな変更は、生成器が正しく初期化されることを確認する部分です。生成器のインタフェイスを変更し、必要に応じて空回りさせる修正を加えました。また、マスタープロセスのみが最初の並べ替えを実行することをチェックするコードも加えられました。

pmaxTのインターフェイスはmt.maxTと同じです。全機能が並列化されました。

·結果

ベンチマークは、6102列の遺伝子と76行のサンプルによるデータセットを用い、150,000並べ替えを実行しました。

計算カーネルのみを見ればほぼ理想的なスケール性を示しました。パラメータのブロードキャストとp値計算の集団通信では、プロセスの増加に伴い時間は増加します。しかしながらHECToRにおいては集団通信が最適化されているため、プロセス数を倍にした際に増える時間は極めて小さくなっています。

全体としてみれば、プロセス数と計算時間の関係は線形と言えます。この関係は列数と計算時間についても同様のことが言えます。

謝辞

このプロジェクトは、NAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびNAG Ltd.のCSEサポートサービスにより管理運営されています。

文献

[1]	dCSE website, “Distributed CSE Support.” Available at : http://www.hector.ac.uk/cse/distributedcse/, accessed 10 March 2010.
[2]	T. Sloan and M. Mewissen, “SPRINTing with HECToR dCSE Proposal.”
[3]	G. Grimes, S. Moodie, J. Beattie, M. Craigon, P. Dickinson, T. Forster, A. Livingston, M. Mewissen, K. Robertson, A. Ross, G. Sing, and P. Ghazal, “GPXMacrophage Expression Atlas: A database for expression profiles of macrophages challenged with a variety of pro-inflammatory, anti-inflammatory, benign and pathogen insults,” BMC Genomics, vol.6, p.178, 2005.
[4]	A. Brazma, H. Parkinson, U. Sarkans, M. Shojatalab, J. Vilo, N. Abeygunawardena, E. Holloway, M. Kapushesky, P. Kemmeren, G. G. Lara, A. Oezcimen, P. Rocca-Serra, and S.-A. Sansone, “ArrayExpress? a public repository for microarray gene expression data at the EBI,” Nucl. Acids Res., vol.31, pp.68-71, 2003.
[5]	The R package, “The R Project for Statistical Computing.” Available at : http://www.r-project.org/, accessed 10 March 2010.
[6]	BBSRC web site, “Biotechnology and Biological Sciences Research Council.” Available at : http://www.bbsrc.ac.uk/.
[7]	EPCC website, “Edinburgh Parallel Computing Centre.” Available at : http://www.epcc.ed.ac.uk/, accessed 10 March 2010.
[8]	DPM website, “Division of Pathway Medicine.” Available at : http://www.ed.ac.uk/schools-departments/pathway-medicine, accessed 10 March 2010.
[9]	SPRINT’s website, “SPRINT: A new parallel framework for R.” Available at : http://www.r-sprint.org/, accessed 10 March 2010.
[10]	Jon Hill, Matthew Hambley, Thorsten Forster, Muriel Mewissen, Terence M Sloan, Florian Scharinger, Arthur Trew and Peter Ghazal, “SPRINT: A new parallel framework for R,” BMC Bioinformatics, December 2008.
[11]	The MPI parallel library, “The Message Passing Interface (MPI) standard.” Available at : http://www.mcs.anl.gov/research/projects/mpi/, accessed 10 March 2010.
[12]	HECToR’s user support, “HECToR : Hardware overview.” Available at : http://www.hector.ac.uk/service/hardware/, accessed 10 March 2010.
[13]	J. Hill, “Parallel statistics with R,” EPCC News, vol.62, p.7, 2008
[14]	M. Dublin, “New software tool lets researchers user R across a computer cluster.” Available at : http://www.genomeweb.com/new-software-tool-lets-researchers-use-r-across-compute-cluster, accessed 29 March 2010, 2009.
[15]	Savvas Petrou, “D1.1 Work Package 1 deliverable. Performance analysis of parallel correlation function pcor.”
[16]	Savvas Petrou, “D2.1 Work Package 2 deliverable. Performance analysis of parallel permutation testing function pmaxT.”
[17]	NWS, “NetWorkSpaces for R.” Available at : http://nws-r.sourceforge.net/, accessed 3 March 2010.
[18]	Rmpi, “Rmpi for R.” Available at : http://www.stats.uwo.ca/faculty/yu/Rmpi/, accessed 3 March 2010.
[19]	R/parallel, “R/parallel framework.” Available at : http://www.rparallel.org/, accessed 3 March 2010.
[20]	papply, “papply - R parallel apply function.” Available at : http://math.acadiau.ca/ACMMaC/software/papply.html, accessed 3 March 2010.
[21]	SNOW, “Simple Network of Workstations for R.” Available at : http://www.stat.uiowa.edu/~luke/R/cluster/cluster.html, accessed 3 March 2010.
[22]	Wellcome Trust web site, “The Wellcome Trust.” Available at : http://www.wellcome.ac.uk/.
[23]	ECDF, “The Edinburgh Compute and Data Facility.” Available at : http://www.ecdf.ed.ac.uk/, accessed 10 March 2010.
[24]	The FF R package, “ff: memory-efficient storage of large data on disk and fast access functions.” Available at : http://cran.r-project.org/web/packages/ff/index.html, accessed 10 March 2010.
[25]	Y. Ge, S. Dudoit and T. P. Speed, “Resampling-based multiple testing for microarray data hypothesis,” TEST, vol.12, pp.1-77, June 2003.
[26]	The multtest R package, “multtest: Resampling-based multiple hypothesis testing.” Available at : http://cran.r-project.org/web/packages/multtest/index.html, accessed 10 March 2010.
[27]	Westfall,P.H. and Young, S.S., “Resampling-Based Multiple Testing: Examples and Methods for P-Value Adjustment.,” p.?340, 1993. Wiley, New York.
[28]	R User Meeting, “SPRINT Project.” Available at : http://www.nesc.ac.uk/action/esi/contribution.cfm?Title=1044, accessed 30 March 2010.
[29]	Data-Intensive Research Workshop, “SPRINT Project.” Available at : http://www.nesc.ac.uk/esi/events/1047/, accessed 30 March 2010.
[30]	HPDC 2010 Workshop, “Emerging Computational Methods for the Life Sciences (ECML).” Available at : http://salsahpc.indiana.edu/ECMLS2010/ Accessed 25 April 2010.
[31]	ACM International Symposium, “High Performance Distributed Computing (HPDC).” Available at : http://hpdc2010.eecs.northwestern.edu/ Accessed 25 April 2010.
[32]	useR! 2010, “The R User Conference 2010.” Available at : http://user2010.org/ Accessed 25 April 2010.
[33]	CRAN, “The Comprehensive R Archive Network.” Available at : http://cran.r-project.org/, accessed 30 March 2010.