Почему Matlab в 11 раз быстрее, чем Stack Overflow
Я сравниваю скорость алгоритма ценообразования Монте-Карло для опциона колл ваниль между Matlab и C ++. Это не то же самое, что Почему MATLAB так быстр в умножении матриц? поскольку ускорение происходит не из-за умножения матриц (есть только точечное произведение, которое выполняется быстро), а из-за его высокоэффективного генератора гауссовых случайных чисел.
В 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
Если я запускаю это с 10 миллионами путей следующим образом
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
Теперь я делаю то же самое в C ++ в VS2013
Мой код находится в классе 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.
Для 10-метровых дорожек это занимает 4,124 секунды.
Это в 11 раз медленнее, чем Matlab.
Кто-нибудь может объяснить разницу?
РЕДАКТИРОВАТЬ: При дальнейшем тестировании замедление, похоже, вызов GaussianRVByBoxMuller. Кажется, у Matlab очень эффективная реализация — метод Ziggurat. Обратите внимание, что BM здесь неоптимальный, поскольку он генерирует 2 RV, и я использую только 1. Исправление только это даст коэффициент ускорения 2.
Решение
В настоящий момент вы генерируете однопоточный код. По-видимому, Matlab использует многопоточный код. Это позволяет ему работать быстрее примерно в N раз, где N = количество ядер в вашем процессоре.
В этой истории есть немного больше, чем это. Еще одна проблема, которая возникает в том, что вы используете 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 это длилось ~ 0,31 секунды.
На процессоре Intel i7 это длилось ~ .16 секунд.
Конечно, если у вас процессор с гораздо большим количеством ядер, у вас есть неплохие шансы, что он будет работать еще немного быстрее.
На данный момент мой код требует VC ++ 2015 вместо 2013, но я сомневаюсь, что это сильно влияет на производительность. Это в основном просто удобство, как использование 10'000'000 вместо 10000000 (но я не собираюсь устанавливать копию 2013 года на эту машину, чтобы просто выяснить, что мне нужно изменить, чтобы приспособиться к ней).
Также обратите внимание, что на последнем процессоре Intel вы можете (или не можете) получить некоторое улучшение, изменив arch:AVX в arch:AVX2 вместо.
Быстрая проверка однопоточного кода показывает, что ваш дистрибутивный код Box-Muller может быть немного быстрее, чем обычный дистрибутивный код стандартной библиотеки, поэтому переход на ориентированную на многопотоковую версию версию может получить чуть большую скорость (и оптимизированную версию). должен примерно вдвое больше).
Другие решения
Других решений пока нет …