c++ - 为什么 Matlab 比 C++ 快 11 倍

Question

我正在比较 Matlab 和 C++ 之间普通看涨期权的蒙特卡罗定价算法的速度。这与 Why is MATLAB so fast in matrix multiplication? 不同因为加速不是由于矩阵乘法(只有一个快速完成的点积)，而是由于其高效的高斯随机数生成器。

在 Matlab 中，代码已经被矢量化，代码如下所示

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

如果我按如下方式使用 1000 万条路径运行它

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

现在我在 VS2013 中用 C++ 做同样的事情

我的代码在 OptionMC 类中，如下所示

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;
}

BM代码如下

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;
}

我已将优化选项设置为优化 Visual Studio 中的速度。

对于 10m 的路径，它需要 4.124 秒。

这比 Matlab 慢 11 倍。

谁能解释一下区别？

编辑:在进一步测试中，减速似乎是对 GaussianRVByBoxMuller 的调用。 Matlab 似乎有一个非常有效的实现 - Ziggurat 方法。请注意，BM 在这里不是最优的，因为它生成 2 个 RV，而我只使用 1 个。仅修复此问题将得到 2 倍的加速。

Answer 1

就目前而言，您正在生成单线程代码。猜测，Matlab 使用的是多线程代码。这使它的运行速度提高了大约 N 倍，其中 N = CPU 中的内核数。

不过，故事的内容远不止于此。出现的另一个问题是您正在使用 rand()，它使用隐藏的全局状态。因此，如果您对代码进行简单的重写以使其成为多线程的，rand() 内部状态的冲突很有可能会阻止您获得很大的速度提升(并且可能很容易运行得更慢——也许慢很多)。

要解决这个问题，您可以考虑(例如)使用 C++11 中添加的新随机数生成(和可能的分布)类。有了这些，您可以为每个线程创建一个单独的随机数生成器实例，防止它们的内部状态发生冲突。

我稍微重写了您的代码以使用这些代码，并调用该函数以获取此信息:

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";
}

我使用 VC++ 2015 编译了这个，使用以下命令行:

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

在 AMD A8-7600 上，它运行了 ~.31 秒。
在 Intel i7 处理器上，运行时间约为 0.16 秒。

当然，如果您的 CPU 有更多的内核，您很有可能会运行得更快一些。

就目前而言，我的代码需要 VC++ 2015 而不是 2013，但我怀疑它对性能的影响很大。这主要是为了方便，比如使用 10'000'000 而不是 10000000 (但我不会在这台机器上安装 2013 的拷贝只是为了弄清楚我需要更改以适应它)。

另请注意，与最近的 Intel 处理器相比，您可能(或可能不会)通过将 arch:AVX 更改为 arch:AVX2 来获得一些改进。

快速检查单线程代码表明您的 Box-Muller 分发代码可能比标准库的正常分发代码快一点，因此切换到线程友好的版本可能会获得更快的速度(并且优化版本应该也大约是它的两倍)。