なぜMatlabはC++よりも11倍高速です

Question

MatlabとC++の間のバニラコールオプションのMonte Carlo価格設定アルゴリズムの速度を比較しています。高速化は行列乗算によるものではなく（高速に実行されるドットプロダクトのみであるため）、非常に効率的なガウス乱数ジェネレータによるものと思われるため、これはWhy is MATLAB so fast in matrix multiplication?と同じではありません。コードがベクトル化されており、

function [ value ] = OptionMCValue(yearsToExpiry, spot, strike, riskFreeRate, dividendYield, volatility, numPaths ) 

    sd = volatility*sqrt(yearsToExpiry); 
    sAdjusted = spot * exp((riskFreeRate - dividendYield - 0.5*volatility*volatility) * yearsToExpiry); 

    g = randn(1,numPaths); 
    sT = sAdjusted * exp(g * sd); 
    values = max(sT-strike,0);` 
    value = mean(values); 
    value = value * exp(-riskFreeRate * yearsToExpiry); 

end 

strike = 100.0; 
yearsToExpiry = 2.16563; 
spot = 100.0; 
volatility = 0.20; 
dividendYield = 0.03; 
riskFreeRate = 0.05; 
oneMillion = 1000000; 
numPaths = 10*oneMillion; 

tic 
value = OptionMCValue(yearsToExpiry, spot, strike, riskFreeRate, dividendYield, volatility, numPaths ); 
toc 

を次のように私は千万のパスでこれを実行すると、私は

Elapsed time is 0.359304 seconds. 
    12.8311 

を取得するには、次のようにコードが見えMATLABで

今私はVS2013のC++で同じことを行います

私のコードはOptionMCクラスであると私はVisual Studioで速度を最適化するために、最適化オプションを設定している

double GaussianRVByBoxMuller() 
{ 
double result; 
double x; double y;; 
double w; 

do 
{ 
    x = 2.0*rand()/static_cast<double>(RAND_MAX)-1; 
    y = 2.0*rand()/static_cast<double>(RAND_MAX)-1; 
    w = x*x + y*y; 
} while (w >= 1.0); 

w = sqrt(-2.0 * log(w)/w); 
result = x*w; 

return result; 
} 

を次のようにBMコードがある

double OptionMC::value(double yearsToExpiry, 
        double spot, 
        double riskFreeRate, 
        double dividendYield, 
        double volatility, 
        unsigned long numPaths) 
{ 
    double sd = volatility*sqrt(yearsToExpiry); 
    double sAdjusted = spot * exp((riskFreeRate - dividendYield - 0.5*volatility*volatility) * yearsToExpiry); 
    double value = 0.0; 
    double g, sT; 

    for (unsigned long i = 0; i < numPaths; i++) 
    { 
     g = GaussianRVByBoxMuller(); 
     sT = sAdjusted * exp(g * sd); 
     value += Max(sT - m_strike, 0.0); 
    } 

    value = value * exp(-riskFreeRate * yearsToExpiry); 
    value /= (double) numPaths; 
    return value; 
} 

を以下のようです。

10mパスの場合、4.124秒かかります。

これはMatlabの11倍です。

誰でもその違いを説明できますか？

EDIT：さらなるテストでは、スローダウンはGaussianRVByBoxMullerの呼び出しであるようです。 Matlabは非常に効率的な実装、つまりZigguratメソッドを持っているようです。ここでBMは2つのRVを生成するので、ここでは最適ではないことに注意してください。私は1だけを使用します。これを固定すると2倍のスピードアップが得られます。

Answer 1

今のところ、シングルスレッドコードを生成しています。推測すると、Matlabはマルチスレッドコードを使用しています。これにより、約N倍の速さで実行できます。ここで、NはCPUのコア数です。

これ以上の話があります。もう1つの問題は、rand()を使用していることです。これは、隠されたグローバルな状態を使用します。したがって、コードをマルチスレッド化するためにコードを書き直すと、rand()の内部状態との競合がスピードの向上を妨げる可能性があります。少し遅い）。

これを解決するには、C++ 11で追加された新しい乱数生成（および可能であれば配布）クラスを使用することなどが考えられます。これらを使用すると、スレッドごとに乱数ジェネレータのインスタンスを個別に作成し、内部状態に対する競合を防ぐことができます。

私はそれらを使用するために少しをあなたのコードを書き直し、この取得するには、関数を呼び出す：

cl -openmp -Qpar -arch:AVX -O2b2 -GL test.cpp 

：私は、次のコマンドラインを使用して、VC++ 2015でこれをコンパイル

double m_strike = 100.0; 

class generator { 
    std::normal_distribution<double> dis; 
    std::mt19937_64 gen; 
public: 
    generator(double lower = 0.0, double upper = 1.0) 
     : gen(std::random_device()()), dis(lower, upper) {} 

    double operator()() { 
     return dis(gen); 
    } 
}; 

double value(double yearsToExpiry, 
    double spot, 
    double riskFreeRate, 
    double dividendYield, 
    double volatility, 
    unsigned long numPaths) 
{ 
    double sd = volatility*sqrt(yearsToExpiry); 
    double sAdjusted = spot * exp((riskFreeRate - dividendYield - 0.5*volatility*volatility) * yearsToExpiry); 
    double value = 0.0; 
    double g, sT; 

    generator gen; 

// run iterations in parallel, with a private random number generator for each thread: 
#pragma omp parallel for reduction(+:value) private(gen) 
    for (long i = 0; i < numPaths; i++) 
    { 
     g = gen(); // GaussianRVByBoxMuller(); 
     sT = sAdjusted * exp(g * sd); 
     value += std::max(sT - m_strike, 0.0); 
    } 

    value = value * exp(-riskFreeRate * yearsToExpiry); 
    value /= (double)numPaths; 
    return value; 
} 

int main() { 
    std::cout << "value: " << value(2.16563, 100.0, 0.05, 0.03, 0.2, 10'000'000) << "\n"; 
} 

を

AMD A8-7600では、これは約0.31秒で実行されました。
インテルi7プロセッサでは、これは〜16秒で実行されました。

もちろん、コアがたくさんあるCPUを使用している場合は、これをかなり高速に実行することができます。

私のコードでは、2013年の代わりにVC++ 2015が必要ですが、パフォーマンスに大きな影響を与えるかどうかは疑問です。 10000000の代わりに10'000'000を使用するような便利さがほとんどです（ただし、このマシンに2013のコピーをインストールして、それに合わせて変更する必要があるかどうかはわかりません）。

最新のIntelプロセッサーよりも、arch:AVXをarch:AVX2に変更すると、改善される可能性があります（またはそうでないかもしれません）。

シングルスレッドコードのクイックチェックは、Box-Mullerディストリビューションコードが標準ライブラリの通常のディストリビューションコードよりも速いことを示しています。そのため、スレッドに適したバージョンに切り替えると、最適化されたバージョンはの約2倍です。