为什么 Matlab 比 C++ 快 11 倍
Why is Matlab 11 times faster than C++
我正在比较 Matlab 和 C++ 之间普通看涨期权的 Monte Carlo 定价算法的速度。这与 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 class 中,如下所示
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 倍的加速。
就目前而言,您正在生成单线程代码。猜测,Matlab 使用的是多线程代码。这允许它 运行 快大约 N 倍,其中 N = CPU.
中的核心数
不过,故事的内容远不止于此。出现的另一个问题是您正在使用 rand()
,它使用隐藏的全局状态。因此,如果您对代码进行简单的重写以使其成为多线程的,那么 rand()
内部状态的冲突很有可能会阻止您获得很大的速度提升(并且可能很容易 运行 更慢——也许慢了很多)。
要解决这个问题,您可以考虑(例如)使用 C++11 中添加的新 运行dom 数字生成(以及可能的分发)类。有了这些,您可以为每个线程创建一个单独的 运行dom 编号生成器实例,防止它们的内部状态发生冲突。
我稍微重写了您的代码以使用这些代码,并调用了该函数以得到:
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 处理器上,这个 运行 在 ~.16 秒内。
当然,如果您的 CPU 具有更多内核,那么您很有可能 运行 宁得更快。
就目前而言,我的代码需要 VC++ 2015 而不是 2013,但我怀疑它对性能的影响很大。这主要是为了方便,比如使用 10'000'000
而不是 10000000
(但我不会在这台机器上安装 2013 的副本只是为了弄清楚我需要更改什么以适应它)。
另请注意,与最近的英特尔处理器相比,您可能(或可能不会)通过将 arch:AVX
更改为 arch:AVX2
来获得一些改进。
快速检查单线程代码表明您的 Box-Muller 分发代码可能比标准库的正常分发代码快一点,因此切换到线程友好的版本可能会获得更快的速度(优化后的版本应该大约是原来的两倍)。
我正在比较 Matlab 和 C++ 之间普通看涨期权的 Monte Carlo 定价算法的速度。这与 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 class 中,如下所示
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 倍的加速。
就目前而言,您正在生成单线程代码。猜测,Matlab 使用的是多线程代码。这允许它 运行 快大约 N 倍,其中 N = CPU.
中的核心数不过,故事的内容远不止于此。出现的另一个问题是您正在使用 rand()
,它使用隐藏的全局状态。因此,如果您对代码进行简单的重写以使其成为多线程的,那么 rand()
内部状态的冲突很有可能会阻止您获得很大的速度提升(并且可能很容易 运行 更慢——也许慢了很多)。
要解决这个问题,您可以考虑(例如)使用 C++11 中添加的新 运行dom 数字生成(以及可能的分发)类。有了这些,您可以为每个线程创建一个单独的 运行dom 编号生成器实例,防止它们的内部状态发生冲突。
我稍微重写了您的代码以使用这些代码,并调用了该函数以得到:
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 处理器上,这个 运行 在 ~.16 秒内。
当然,如果您的 CPU 具有更多内核,那么您很有可能 运行 宁得更快。
就目前而言,我的代码需要 VC++ 2015 而不是 2013,但我怀疑它对性能的影响很大。这主要是为了方便,比如使用 10'000'000
而不是 10000000
(但我不会在这台机器上安装 2013 的副本只是为了弄清楚我需要更改什么以适应它)。
另请注意,与最近的英特尔处理器相比,您可能(或可能不会)通过将 arch:AVX
更改为 arch:AVX2
来获得一些改进。
快速检查单线程代码表明您的 Box-Muller 分发代码可能比标准库的正常分发代码快一点,因此切换到线程友好的版本可能会获得更快的速度(优化后的版本应该大约是原来的两倍)。