GPU 上での同じ処理の繰り返しからなる確率的アルゴリズムの解析とその応用飯田正成 / 仁木直人 (東京理科大学工学研究科 / 工学部) スレッド内ループの繰返し数が確率変動する算法におけるこの結果から，以下のことがいえる． warp (SIMD group of cores) の漸近的寿命分布を導出する．ま (1) アドレス方式の影響は機種によって異なる．Titan ではどのた，モンテカルロ計算に必要な棄却法による乱数発生を実装し条件でも大きな差がない．一方，Tesla では，k が小さい時はあて，計算時間分布を調べた実験結果を報告する．まり差がないが，p, k がともに大きい時には影響が見える． (2) 標準的アドレス方式が確率的算法においても推奨される． 1 はじめに対象とするのは，多数のスレッドが同じ処理を行い，かつ，各ス 4 実証実験例 −−− 棄却法による乱数発生レッドに含まれるループの実行回数 (試行回数と呼ぶ) や条件分 ∃ 密度関数 g(x) , h(x) 間に p ∈(0, 1), p g(x) ≤ h(x) なる関係があ岐による実行パスが確率的に変動する算法である. このような，るとき，棄却法による g(x) に従う k 個の乱数発生とは，以下の確率的算法は，様々なモンテカルロ計算や統計計算に現れる． GPU の SIMD 処理方式では，warp 内の全 m スレッド中，最ような手続きをいう． (1) h(x) に従う乱数 X を生成する. も遅いものの計算終了が warp における処理時間を決める． (2) 区間 (0,1) 上の一様乱数 U を生成する. 2 Warp での処理時間モデルと漸近的性質 ([2]) (3) U ≤ p g(X)/h(X) なら X を出力．そうでなければ (1) へ. Warp の生存時間 (処理終了時間) S は (4) k 個の乱数が出力されたら終了，そうでなければ (1) へ . Y ∑ この時 , 試行毎に乱数を得る確率 ( 受容率 ) は p である．また， (1)S = R0 + Rj , Y = max {X1, . . . , Xm} (3) ループ中の条件分岐による実行パスは (一様乱数発生に関す j=1 る部分を除き ) 少なくとも 3 ある．と書けるものとする．ここに，Xi は第 i スレッド (i = 1, . . . , m) ここで扱うのは, 正規乱数 [3] (可能な実行パス 3) と形状母数が予め与えられた同じ計算を k 回終えるまでの試行回数，Rj 3 (j = 1, 2, . . . ) は j 回目の試行に要する時間，R0 は試行回数に 0.5 のガンマ乱数 [1] (同 3 ) であり，それぞれ実線と破線で表す．無関係な処理時間の合計をそれぞれ表す確率変数である．SIMD 4.1 共通メモリへの出力や資源競合がない場合方式から，R0, Rj は各スレッドに共通であり，条件分岐パスのモデル検証実験として，単一ブロック単一 warp のみ ( n = 32 ) 選択確率がすでに k 回の計算を終えたスレッド数 c に依存するを動作させ，かつ，書込みをレジスタに限った場合に， 10,000 回ことから，Rj も c に依存した分布に従う．の乱数発生に要した平均時間 (単位 clock) を以下に示す．ここで，誤りのない算法で成立すると思われる仮定を置く． clock clock Tesla Titan (1) 1 成功までの試行回数は互いに独立同分布に従い，その r 次積率 µr は全て有限である． (2) Rj は有限区間内に分布し，E[R1] ≥ E[R2] ≥ . . . ．すると，各スレッドに与える計算回数 k → ∞ のとき，1 個当た k k りの計算時間 S/k の分布の特性関数 λ(t) は d 図の印象から，モデルは warp での処理時間を概ねうまく説明し λ(t) → ν(−i log ρ(t)) ているように思える．で漸近近似される．ただし，ν(t), ρ(t) はそれぞれ Y /k の極限 4.2 現実的な状況における平均処理時間分布 (m 個の正規分布の最大値)，R1 が従う分布の特性関数であ各スレッドへの負担 k が大き過ぎると，資源の有効利用が図れり，i は虚数単位である．また，さらに仮定ず，総処理時間が長くなることは容易に想像できる．また，共通 (3) R1 は 1 点分布である．メモリへのアクセスにも影響することも考えられるため，乱数を追加すれば，S/k は近似的に (m = 32 のとき，形状母数が 13 総数を 14 × 2ℓ (ℓ = 20, 21, . . . , 24)，ブロック内スレッド数をあるいはやや大きい) 3 母数ガンマ分布に従う．そして，この結 256 として, 発生した乱数を順次標準的アドレス方式により共通論は，実際に見られる S の分布をよく説明している．メモリに書込む実験を行った結果 (1 乱数当りの実計算時間，単算法設計における意味は，以下のようにまとめられよう．位 ns) を以下に示す． p ns ns (1) 平均試行回数は E[Y /k] → µ1 (k → ∞) である． Tesla Titan p (2) 平均処理時間は E[S/k] → µ1 E[R1] (k → ∞) である． (3) E[Xi/k] = µ1 であることを考慮すれば, 各スレッドに与える負荷 k をある程度以上大きくすることにより, 対象とする確率的算法は GPU 上で十分効率的に動作する. k k æ 20 000 æ 20 000 æ æ æ 15 000 æ 15 000 æ æ æ æ æ æ æ æ æ æ 10 000 æ æ æ æ æ æ à à à à 500 1000 10 000 à à à à 5000 5000 à 5 à à 10 à à 50 à 100 à à 500 1000 æ 10 æ à æ à ì à æ æ à ô ò ì æ ì ò ô à æ à ô ì ò æ ì ò à ô æ à ô ì ò æ à ì æ à ì ò æ à ò ì ô æ à ò ô ì æ ì à ô ò æ à ò ì ô æ ô ò à ì æ à ì ô ò æ à ô ò ì à æ ò ô ì æ ì ô à ò à ò ô ì à ô ì ò ì ò ô ò ì ô ò ô ò ô à à 10 à 50 100 æ l = 20 4 10 100 1000 104 à l = 21 à æ à æ ô ô 2 1 書込みアドレスの決定方法 5 6 2 3 à 8æ æ à æ ì ò 4 à 10 8 æ 6 à 1 ì l = 22 æ ì à ò à ô æ æ à ì ì æ à ì ò æ à ì æ æ ò ô à ì à ô ò à ì ò ô ò æ ì æ ì æ ò æ æ à ì ô à à ì à ì ò ò ô ô ò ô 10 100 æ à ò ì ô æ à ô ò ì æ à ì ò ô ô æ à ò ì æ à ì ò ô ô à æ ò ì 1000 ò l = 23 ò à ì ô ò ì ò ô ô ò ô ô ô ò à ì ô ì ò ô ô l = 24 104 なお, 共通メモリでなくレジスタに書込む場合, どちらの乱数生メモリアクセス速度は総処理速度に大きく影響し，その隠匿方成法でも k の増加による処理時間増加はない. 法も GPGPU 設計に不可欠である．中でも，低速な共通メモリ他の実装実験結果を含めて考えれば，次のことが言えよう． (global memory) への出力は，注意を払う必要がある． (1) k を必要以上に大きくすると，全体の計算効率が下がる可グリッド内通し番号 0 ≤ t ≤ n − 1 を持つスレッドからの第 i 出能性があり, k の値は数十 (16 ∼ 64) 位がよさそうである．力の書込みアドレス a は，a = n (i − 1) + t とすることが推奨 (2) 総スレッド数による影響は大きく無い. されている．しかし，この方式は同時に連続したメモリへの出力を行う前提があり，確率的算法にも適するか検証が必要である．なお，Intel core i7 2600K での平均乱数発生時間が，正規乱数: 107 ns, ガンマ乱数 : 123 ns であることから， Tesla で 20 ∼ 30 倍，ただし，メモリ関連ハードウェアが機種によって異なるため，検 Titan で 40 ∼ 50 倍程度の高速化ができている．証は実験によるしかない．以下の図は，独立に確率 p = 1/µ1 で書込みが起きる実験を参考文献 Tesla (Tesla C2070, 倍精度), Titan (GTX Titan, 倍精度) 上で行った [1] J. H. Arenas and U. Dieter: Computer methods for sampling from 20 ときの平均書込み時間を示す．両者とも，n = 14 × 2 /k , ブロッ gamma, beta, Poisson and binominal distributions, Computing, ク内スレッド数は 256 である．また，実線は a = n (i − 1)+ t, 12, 223-246 (1974). 破線は a = k t + i − 1 のアドレス方式に対応し，各線は下から [2] M. Iida and N. Niki: Lifespan distribution of SIMD groups on p = 1, 0.8, 0.5, 0.2 のときである． GPU engaged in a class of probabilistic computation (to appear ns ns in J. Jpn. Soc. Comput. Statist. ). Tesla Titan [3] A. J. Kinderman and J. F. Monahan: Computer generation of random variables using the ratio of uniform deviates, ACM Trans. Math. Soft., 3, 257-260 (1977). 8ò 8 ò ì à 6æ ì ì à à æ æ 4 ì à ì æ à ò 6 ò ò ò ì ì à ì à ì æ æ à æ à æ ò ò ì à æ ò ò ì æ ò ò à ì ì ì à à æ æ à ò ò ì æ ì à æ ì à æ ì ì à à à æ æ ò ò ì à æ à ò ò ì æ à à ò 4 æ æ ì à ì æ à ò ò ì à æ à æ ò ì à æ à æ 2 ò ì à æ à æ 2 1 5 10 50 100 500 1000 k ò ì æ à ì à ò ì æ à ì à æ æ ò æ ì à ì à ò æ ì à à ò æ ì à à ò æ ì à æ æ æ æ 500 1000 k 5 10 50 100