`parfor` を使用したモンテカルロコードの高速化

ライブスクリプトを開く

この例では、parfor ループを使用してモンテカルロコードを高速化する方法を説明します。モンテカルロ法は、物理学、数学、生物学、金融学など、さまざまな分野で使用されています。モンテカルロ法では、ランダムに分散された入力によって 1 つの関数が何度も実行されます。Parallel Computing Toolbox を使用すると、for ループを parfor ループに置き換えて、コードを簡単に高速化できます。

この例では、ドルオークションに基づいて単純な確率的シミュレーションを実行します。モンテカルロ法を使用し、複数のシミュレーションを実行して 1 ドル紙幣の市場価値を求めます。この例では、ドルオークションは、ランダム過程に依存する出力を生成するブラックボックス関数として扱われます。このモデルの詳細については、ドルオークションを参照してください。モンテカルロコードを高速化するための一般的な方法を確認するには、parfor ループを使用した市場価値の推定を参照してください。

ドルオークション

ドルオークションとは、1971 年に Martin Shubik によって紹介されたノンゼロサムゲームです。ゲームの中で、参加者は 1 ドル紙幣に入札します。1 人の参加者が入札後、他の参加者はそれぞれ前の入札者より高額での入札を選ぶことができます。入札を決断する参加者がいなくなると、オークションは終了します。最高額入札者が 1 ドル紙幣を受け取りますが、通常のオークションとは異なり、最高額の入札者と 2 番目に高い金額を提示した入札者の両方が自分の入札額をオークション主催者に渡します。

確率的モデル

確率的モデルを使用すると、ドルオークションに似たゲームをモデル化できます。状態 (現在の入札額と参加者数) は、マルコフ過程を使用してモデル化できるため、結果 (市場価値) は適切に定義された統計量をもつと予測できます。結果は条件付き分布から導き出されるため、ドルオークションはモンテカルロ解析には理想的です。市場価値は以下の要因に影響されます。

参加者数 (nPlayers)
参加者の行動

この例では、以下のアルゴリズムで、状態次第で参加者がどのような行動を取るか (入札するか離脱するか) を判定します。

入札額を前の入札額プラス incr に設定します。
前の入札者以外の参加者の中から、1 人の参加者をランダムに選択します。
まだ入札が行われていない場合は、8 に進みます。
前の入札額が 1 未満の場合は、0 から 1 の間の乱数を生成します。この乱数が dropoutRate 未満の場合は、7 に進みます。
参加者が落札した場合に稼ぐことのできる金額 (gain) を計算します。
参加者が 2 番目に高い金額の入札者となった場合に失う金額 (loss) を計算します。gain が loss より大きい場合は、8 に進みます。
参加者が離脱します。その参加者を参加者グループから削除して、9 に進みます。
参加者が入札します。
2 人以上の参加者が残っている場合は、手順 1 に進みます。

サポート関数 dollarAuction がドルオークションをシミュレートします。コードを表示するには、dollarAuction.m を参照してください。この関数は、nPlayers、incr および dropoutRate の 3 つの入力を取ります。それぞれの値を設定します。

nPlayers = 20;
incr = 0.05;
dropoutRate = 0.01;

関数 dollarAuction を実行することにより、ランダムなシナリオを実行します。出力 bids および dropouts を保存します。

[bids,dropouts] = dollarAuction(nPlayers,incr,dropoutRate);

ゲームが進むにつれ、入札する参加者と離脱する参加者が出てきます。入札額が 1 を超えた場合、残り 1 名になるまで、参加者たちは "入札合戦" に陥ります。

table dropouts には、2 つの変数が含まれており、Player は各参加者に割り当てられている一意な番号、Epoch は、Player が何周目の入札時に離脱したかを示します。findgroups を使用して dropouts.Epoch をグループ化し、splitapply を使用して dropouts.Epoch のそれぞれの周回で離脱した参加者の数を取得します。

[G,epochs] = findgroups(dropouts.Epoch);
numberDropouts = splitapply(@numel,dropouts.Epoch,G);

最初は、離脱者はいません。この情報を、epochs と numberDropouts の前に 1 と 0 を付けることによって追加します。

epochs = [1;epochs];
numberDropouts = [0;numberDropouts];

nPlayers と cumsum を使用して、残りの参加者数を numberDropouts から計算します。incr と epochs を使用して入札額を計算します。stairs を使用して、numberDropouts の累積和に対して入札額をプロットします。

playersRemaining = nPlayers - cumsum(numberDropouts);
stairs(incr*epochs,playersRemaining)
xlabel('Bid')
ylabel('Number of players remaining')

モンテカルロ法を使用した市場価値の推定

モンテカルロ法を使用することにより、紙幣の市場価値を値 origValue で推定できます。ここでは、モンテカルロモデルを生成し、Parallel Computing Toolbox を使用する場合と使用しない場合の速度を比較します。結果をランダムにサンプリングするために使用する試行回数 nTrials を設定します。

nTrials = 10000;

サポート関数 dollarAuction を複数回実行すると、起こり得る結果をサンプリングできます。for ループを使用して nTrials 個のサンプルを生成し、各試行の最後の入札額を B に保存します。関数 dollarAuction を実行するたびに、異なる結果が得られます。しかし、関数の実行回数が多くなると、すべての実行から得られる結果は適切に定義された統計量をもちます。

nTrials 回のシミュレーションの計算が完了するまでの所要時間を記録します。経過時間の統計ノイズを削減するために、この処理を 5 回繰り返して、最短の経過時間を取ります。

t = zeros(1,5);
for j = 1:5
    tic
    B = zeros(1,nTrials);
    for i = 1:nTrials
        bids = dollarAuction(nPlayers,incr,dropoutRate);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
forTime = min(t)

forTime = 21.4323

histogram を使用して、最終入札額 B のヒストグラムをプロットします。xline を使用して、元の値 (1 ドル) と mean によって得られた平均市場価値をプロットに重ね合わせます。

histogram(B);
origLine = xline(1,'k','LineWidth',3);
marketLine = xline(mean(B),'k--','LineWidth',3);
xlabel('Market value')
ylabel('Frequency')
legend([origLine, marketLine],{'Original value','Market value'},'Location','NorthEast')

与えられたアルゴリズムと入力パラメーターでは、平均市場価値は元の価値より大きくなりました。

`parfor` ループを使用した市場価値の推定

Parallel Computing Toolbox を使用すると、モンテカルロコードを簡単に高速化できます。まず、'Processes' プロファイルを使用して、4 つのワーカーをもつ並列プールを作成します。

p = parpool('Processes',4);

Starting parallel pool (parpool) using the 'Processes' profile ...
Connected to parallel pool with 4 workers.

for ループを parfor ループに置き換えます。nTrials 回のシミュレーションの計算が完了するまでの所要時間を記録します。経過時間の統計ノイズを削減するために、この処理を 5 回繰り返して、最短の経過時間を取ります。

t = zeros(1,5);
for j = 1:5
    tic
    parfor i = 1:nTrials
        bids = dollarAuction(nPlayers,incr,dropoutRate);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
parforTime = min(t)

parforTime = 5.9174

ワーカーが 4 つの場合、parfor ループを使用するとコードを 3 倍以上高速で実行できることが結果によって示されています。

`parfor` ループでの乱数による再現可能な結果の生成

parfor ループで乱数を生成すると、ループを実行するたびに異なる結果が生成される可能性があります。再現可能な結果を作成するには、ループの各反復で、乱数発生器の状態が確定的でなければなりません。詳細については、parfor ループでの乱数の繰り返しを参照してください。

サポート関数 dollarAuctionStream は、4 番目の引数 s を取ります。このサポート関数は、指定されたストリームを使用して乱数を生成します。コードを表示するには、dollarAuctionStream.m を参照してください。

ストリームを作成するとき、そのストリームのサブストリームは統計的に独立しています。詳細については、RandStreamを参照してください。コードで生成される結果の分布を毎回同じにするには、ループの反復ごとに乱数発生器ストリームを作成し、Substream プロパティをループインデックスに設定します。dollarAuction を dollarAuctionStream に置き換え、s を使用してワーカー上で dollarAuctionStream を実行します。

t = zeros(1,5);
for j = 1:5
    tic
    parfor i = 1:nTrials
        s = RandStream('Threefry');
        s.Substream = i;
        bids = dollarAuctionStream(nPlayers,incr,dropoutRate,s);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
parforTime = min(t)

parforTime = 8.7355

デスクトップからクラスターへのスケールアップ

デスクトップから、より多くのワーカーをもつクラスターへ、コードをスケールアップすることができます。デスクトップからクラスターへのスケールアップの詳細については、デスクトップからクラスターへのスケールアップを参照してください。

delete を使用して既存の並列プールをシャットダウンします。

delete(p);

サポート関数 dollarAuctionStream を parfor ループ内で計算します。同じ parfor ループを異なる数のワーカーで実行し、経過時間を記録します。経過時間の統計ノイズを削減するために、parfor ループを 5 回実行して、最短の経過時間を取ります。最短時間を配列 elapsedTimes に記録します。次のコードの MyCluster をクラスタープロファイルの名前に置き換えます。

workers = [1 2 4 8 16 32];
elapsedTimes = zeros(1,numel(workers));

% Create a pool using the 'MyCluster' cluster profile
p = parpool('MyCluster', 32);

Starting parallel pool (parpool) using the 'MyCluster' profile ...
Connected to the parallel pool (number of workers: 32).

for k = 1:numel(workers)    
    t = zeros(1,5);
    for j = 1:5
        tic
        parfor (i = 1:nTrials, workers(k))
            s = RandStream('Threefry');
            s.Substream = i;
            bids = dollarAuctionStream(nPlayers,incr,dropoutRate,s);
            B(i) = bids.Bid(end);
        end
        t(j) = toc;
    end
    
    elapsedTimes(k) = min(t);
end

Analyzing and transferring files to the workers ...done.

elapsedTimes(1) を elapsedTimes の時間で除算して、計算の高速化を計算します。ワーカー数に対する高速化をプロットすることにより、ストロングスケーリングを調べます。

speedup = elapsedTimes(1) ./ elapsedTimes;
plot(workers,speedup)
xlabel('Number of workers')
ylabel('Computational speedup')

計算の高速化はワーカー数に伴って上昇します。

parfor を使用したモンテカルロ コードの高速化

ドル オークション