並列計算
Wolframシステムコンパイラは計算を並列で行うことができる.並列計算はコンパイルされた関数をデータのリストに並列に縫い込むことで行う.まず,Listable属性のある関数をコンパイルする.
cS = Compile[{{x}}, x ^ 2, RuntimeAttributes -> {Listable}]入力がコンパイルされた関数のタイプ指定と一致する場合は通常の動作をする.次の例題では,実数入力がコンパイルされた関数のタイプ指定と一致するので関数が実行される.
cS[2.5]ここでは,コンパイルされた関数がリスト入力を受け取る.これは入力よりも高い階級なので,入力に縫い込まれる.
cS[{1, 2, 3, 4}]リスト可能なコンパイルされた関数も並列で実行することができる.これはParallelizationオプションをTrueにして行う.
cP = Compile[{{x}}, x ^ 2, RuntimeAttributes -> {Listable}, Parallelization -> True]入力がコンパイル関数のタイプ指定と一致する場合,その入力は正常に働く.
cP[2.5]次では,コンパイルされた関数が入力のリストを受け取っている.これは,コンパイルされた関数が縫い込まれた入力よりも階数が高いためである.これもまた並列で実行される.
cP[{1, 2, 3, 4}]計算の高速化
並列化を使う主な理由は計算の高速化である.これは,多くの場合うまくいく.
arg = Range[ -50., 50, 0.02];
cS = Compile[{{x}},
Module[{sum = 1.0, inc = 1.0}, Do[inc = inc * x / i;sum = sum + inc, {i, 10000}];sum],
RuntimeAttributes -> {Listable}, Parallelization -> False];
cS[arg];//AbsoluteTiming次は,上と等価の並列で実行される関数をコンパイルする.筆者のマシン上では速度が2倍になる.
cP = Compile[{{x}},
Module[{sum = 1.0, inc = 1.0}, Do[inc = inc * x / i;sum = sum + inc, {i, 10000}];sum],
RuntimeAttributes -> {Listable}, Parallelization -> True];
cP[arg];//AbsoluteTimingコンパイルされた関数を実行する並列スレッドのデフォルト数は$ProcessorCountの設定値である.
$ProcessorCountコンパイルのターゲット
CompilationTargetをCにすると並列にコンパイルされた関数をさらに高速化することができる.
cP = Compile[{{x}},
Module[{sum = 1.0, inc = 1.0}, Do[inc = inc * x / i;sum = sum + inc, {i, 10000}];sum],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True];
cP[arg];//AbsoluteTiming外部呼出し
並列計算で重要な問題に同じことをしようとしている2つの異なるスレッドをどのように扱うかがある.コンパイルされた関数を使った並列計算の場合これは自動的に処理される.
fFun[x_] := (cnt++;x ^ 2)
cnt = 0;cP = Compile[{{x}}, fFun[x],
RuntimeAttributes -> {Listable}, Parallelization -> True];次は入力データのリスト上でコンパイルされた関数を並列に実行する.
cP[{1, 2, 3, 4, 5}]cnt事実,並列実行されているコンパイルされた関数が外部呼出しを行うときは.実際に呼出しを行うのが1度に1つのスレッドだけになるように常に同期されたプリミティブによって行われる.これは,並列にコンパイルされた関数が多くの外部呼出しを行う場合は,並列による高速化がうまく行かないことを示している.
乱数
モンテカルロ法等の乱数を使う計算の多くは,並列で行うと計算を速めることができる.しかし,乱数を速く効率的に並列で使用するには,独立で作動し,他のスレッドで生成される数とは統計的に独立している乱数を生成する別々の実行スレッドに生成器がなくてはならない.このため,Wolfram言語がデフォルトで並列計算に使用する乱数生成器は,連続計算に使われるものとは異なり,必然的に実際の乱数も異なったものになる.さらに,任意の並列計算について,計算の特定部分に割り当てられる実行スレッドは,実行ごとに違うものになることがあり,このため,最初は同じ乱数の状態であっても,結果は実行ごとに異なることがある.
以下では,格子上の n ステップのランダムウォーク(複数のシミュレーションについて1つは連続して実行され,もう1つは並列で実行される)のシミュレーションを行う2つのCompiledFunctionオブジェクトを作成する.
{rws, rwp} = Map[Compile[{{n, _Integer}}, Total[RandomChoice[{{1, 0}, {-1, 0}, {0, 1}, {0, -1}}, n]], RuntimeAttributes -> Listable, Parallelization -> #]&, {False, True}];連続実行するものは,BlockRandom内で実行するたびに毎回同じ結果を返す.
BlockRandom[rws[ConstantArray[100, 20]]]BlockRandom[rws[ConstantArray[100, 20]]]p1 = BlockRandom[rwp[ConstantArray[100, 20]]]p2 = BlockRandom[rwp[ConstantArray[100, 20]]]違いの大部分は,結果の順番の違いにある.Intersectionを使うと,多くは実際には同じであることが分かる.
Intersection[p1, p2]同じであるものは,毎回同じ実行スレッドで実行されたものである.
通常BlockRandomあるいはSeedRandomを使う場合,並列計算外部でこれを使うことが望ましい.これらのコマンドを並列評価内で使いたい場合は,現行の実行スレッドだけに作動する.このことについては,「乱数生成」 の 「並列計算のSeedRandomおよびBlockRandom」のセクションに説明されており,心に留めておいたほうがよい.
SeedRandom[seed,Method->"ParallelGenerator"]を使って,並列計算用の乱数生成器を変更することができる.並列計算のデフォルトの乱数生成器は,高品質の乱数を生成する一連の1024のメルセンヌツイスタ生成器である.Wolfram言語は,並列計算のそれぞれの実行スレッドに毎回異なるこれらの生成器のいずれかを使う.これらの重要な機能は,それぞれが他から独立して乱数を生成するということで,このためそれぞれの実行スレッドで行われる計算間で相関関係はないということである.このことを示すには,ブロック検定と呼ばれる乱数の標準検定を行えばよい.
ブロック検定は,任意分布から生成された数の標本平均が中心極限定理で収束すべきように,正規分布に収束することを確認する.収束しない場合は,分布に問題のある可能性がある.
以下は,等しい確率で0か1になる n 個の整数からの標本平均を得て,リスト引数を与えられると並列実行を行うCompiledFunctionを定義する.
smean = Compile[{{n, _Integer}}, Total[RandomInteger[1, n]] / N[n], RuntimeAttributes -> Listable, Parallelization -> True];これを使って,連続と並列で n ビットの m 個の標本平均を生成する関数を定義する.
serial[n_, m_] := Table[smean[n], {m}];
parallel[n_, m_] := smean[ConstantArray[n, m]];中心極限定理から,標本の合計は同じ平均と標準偏差
を持つ正規分布に従わなくてはならない.
sumdist[n_] := Module[{u = BernoulliDistribution[1 / 2]}, NormalDistribution[Mean[u], StandardDeviation[u] / Sqrt[n]]];n = 10 ^ 4; m = 10 ^ 3;
{AbsoluteTiming[sdata = serial[n, m];], AbsoluteTiming[pdata = parallel[n, m];]}予想される分布の確率分布関数と比べた,連続と並列のデータのヒストグラムを示す.
pdf = Plot[PDF[sumdist[n], x], {x, .45, .55}, PlotRange -> All];Row[{Show[Histogram[sdata, Automatic, "ProbabilityDensity"], pdf], Show[Histogram[pdata, Automatic, "ProbabilityDensity"], pdf]}]視覚的にはどちらがいいかを判断することは難しい.DistributionFitTestを使って,データが予測される分布と同じ分布に従う帰無仮説
と同じ分布に従わない対立仮説
を持つ適合度仮説検定の
値を得る方が望ましい方法である.
{DistributionFitTest[sdata, sumdist[n]], DistributionFitTest[pdata, sumdist[n]]}DistributionFitTestの検定統計は,それ自体が(連続)一様分布でなければならないので,一番よい情報は,何度も連続実行と並列実行を行って比較することから得られるものである.
{First[AbsoluteTiming[pdata = Table[DistributionFitTest[parallel[n, m], sumdist[n]], {10 ^ 3}];]], DistributionFitTest[pdata, UniformDistribution[]], Histogram[pdata, Automatic, "ProbabilityDensity"]}{First[AbsoluteTiming[pdata = Table[DistributionFitTest[serial[n, m], sumdist[n]], {10 ^ 3}];]], DistributionFitTest[pdata, UniformDistribution[]], Histogram[pdata, Automatic, "ProbabilityDensity"]}並列制御
コンパイルされた関数は複数の実行スレッドを使って並列に実行する.スレッド数はまず$ProcessorCountによって設定される.
$ProcessorCount実際の設定はSystemOptionsで"ParallelOptions"の"ParallelThreadNumber" サブオプションを使って変更することができる.
SystemOptions["ParallelOptions" -> "ParallelThreadNumber"]次は並列にコンパイルされた関数のデモンストレーションである.
arg = Range[ -50., 50, 0.02];
cP = Compile[{{x}},
Module[{sum = 1.0, inc = 1.0}, Do[inc = inc * x / i;sum = sum + inc, {i, 10000}];sum],
RuntimeAttributes -> {Listable}, Parallelization -> True];
cP[arg];//AbsoluteTimingここでは,スレッド数が1に設定されている.このため,コンパイルされた関数が強制的に連続実行されることになる.
SetSystemOptions["ParallelOptions" -> "ParallelThreadNumber" -> 1]実行時間はParallelizationがFalseに設定されたかのようのに等しくなる.
cP[arg];//AbsoluteTiming