Имитационное моделирование с использованием методов Монте-Карло в наше время используется практически во всех областях операционной деятельности, где требуется многократное принятие решений по итогам анализа поступающих из внешнего мира данных. При этом важную роль начинает играть качество, производительность и доступность генераторов случайных чисел, использующихся для придания абстрактному методу черт реальной задачи, решаемой специалистом. Как я недавно выяснил, этот вопрос начинает играть решающее значение при переходе к параллельному программированию… Вы тоже столкнулись с этой проблемой, и хотите знать, как в Windows можно быстро получить массивы случайных чисел с нужным распределением?
/* В статье используется значительное количество цитат на английском языке, я предполагаю достаточный уровень языковой подготовки читателя */
Одно из моих приложений по долгу службы выполняет прогон имитационных моделей на длинной последовательности факторов из базы данных. Не вдаваясь в подробности — недавно возникла задача изучения вероятностных характеристик повторяющегося бинарного (двухисходового) процесса, обладающего сильным элементом случайности, по истории его фактических реализаций с целью снижения неопределённости в предсказании будущих исходов. Реализации процесса считаются независимыми друг от друга и одинаково (нормально) распределёнными, процесс считается стационарным. Цель — по возможности наиболее точное предсказание результатов N подряд идущих реализаций процесса, — подразумевает оценку средней и наименьшей ожидаемых частот положительных реализаций процесса за период длины N.
Тут можно использовать стандартный статистический подход с расчётом доверительного уровня, на котором реализации случайной величины не отклонятся от её матожидания больше чем на k величин её стандартного отклонения. Данный подход основывается на идее сходимости частот к вероятностям (теореме Бернулли). Конечно, очень заманчиво рассчитать средние значения на основе нескольких сот тысяч зафиксированных в базе данных реализаций процесса, но, поскольку у меня есть ограничение на длину подвыборки N (порядка нескольких десятков), случайная природа процесса на ограниченной выборке наверняка преподнесёт мне сюрпризы в виде отклонения от среднего. Именно поэтому мне нужно знать с заданной долей уверенности, на какие максимально возможные отклонения от среднеожидаемой частоты следует «закладываться» при заданном N.
Требуемые оценки я решил получать методом статистического бутстреппинга (bootstrapping) — многократной генерации выборок по N реализаций методом Монте-Карло из основной совокупности примеров. Суть метода состоит в том, чтобы из имеющейся выборки сформировать достаточно большое количество псевдовыборок, состоящих из случайных комбинаций исходного набора элементов (в результате в одной псевдовыборке некоторые исходные элементы могут встретиться несколько раз, тогда как другие могут вообще отсутствовать), и для каждой полученной псевдовыборки определить значения анализируемых статистических характеристик с целью изучения их разброса.
По предварительным оценкам, мне потребуется генерация ~ 10 млн псевдовыборок для каждой из тысяч подгрупп параметров. С учётом количества данных и других факторов, речь идёт о десятках часов вычислительного времени сервера. На Visual Basic 6 был написан и протестирован прототип приложения-симулятора. Однако после компиляции мультипоточной боевой версии внезапно обнаружилось, что при запуске симуляции в несколько потоков загрузка в расчёте на ядро CPU составляет всего 5% вместо ожидаемых 100%.
Basic Hotspots Analysis в профилировщике Intel Vtune Amplifier не показал ничего необычного. Я был уверен, что симулятор должен полностью загружать все потоки работой, поскольку ничего кроме арифметических операций и доступа к памяти по случайным индексам код не содержал:
For k = 1 To NumIters
Randomize
For i = 1 To NumPockets
NumFired = 0
For j = 1 To chainLength
ind = d + Int(B * Rnd)
NumFired = NumFired + FiredOrNot(ind)
Next j
Tmp = NumFired / chainLength
If Tmp < MinProb Then MinProb = Tmp
If Tmp > MaxProb Then MaxProb = Tmp
Values(i) = Tmp
TotalFired = TotalFired + NumFired
Next i
AvgProb = TotalFired / (k * NumPockets) / chainLength
Next k
Неужели этот код настолько завязан на пропускную способность памяти (memory-bound)?! А зачем же тогда нужны многомегабайтные процессорные кэши, почему они не справляются?! Проблема поставила меня в тупик на несколько дней. Запоздало запущеный Concurrency Analysis выявил огромную просадку где-то в недрах системной библиотеки VB, но не смог показать мне, из какого места моей программы вызываются эти затратные системные функции:
Наконец, опытным путём я выяснил, что проблема заключалась в операторе Rnd! По всей видимости, параллельное использование Rnd не предусмотрено разработчиками; возможно, реализация Microsoft использует какие-то внутренние ожидающие блокировки. По факту, распараллеливание Rnd вместо ускорения замедляло работу приложения в ТЫСЯЧИ раз! Думаю, аналогичное поведение будет свойственно и продуктам из пакета MS Office, использующих VBA. Никогда ещё вопрос быстрой генерации случайных целочисленных индексов не был для меня так актуален! Я решил поискать готовые библиотеки для генерации массивов случайных чисел.
Выбор пал на Intel® Math Kernel Library с её пакетом Statistical Functions.
MKL доступна в испытательном режиме. В настоящее время (версия 11.02) доступны 11 базовых генераторов случайных величин, 12 непрерывных, 11 дискретных выходных распределений. Что для меня оказалось очень важным, разработчик даёт гарантию потокобезопасности библиотеки, более того, указывает на преимущества своего продукта при использовании в мультипоточном режиме… В документации MKL моё внимание привлёк следующий раздел:
Non-deterministic random number generator (RDRAND-based generators only) [AVX], [IntelSWMan]. You can use a non-deterministic random number generator only if the underlying hardware supports it. For instructions on how to detect if an Intel CPU supports a non-deterministic random number generator see, for example, Chapter 8: Post-32nm Processor Instructions in [AVX] or Chapter 4: RdRand Instruction Usage in [BMT].
Давно хотелось попробовать эту «железную» возможность, но, к сожалению, оказалось, что мои процессоры не оснащены аппаратным генератором случайных чисел. Помимо перечисленного, важной особенностью пакета MKL является парадигма блочной генерации потоков случайных чисел, в отличие от привычной практики получения одного случайного значения за один вызов оператора в большинстве языков программирования.
В сети уже доступны рабочие примеры использования случайных чисел из MKL VSL для имитационного моделирования под сопроцессорную архитектуру Xeon Phi, например, статья Shuo Li "Case Study: Achieving High Performance on Monte Carlo European Option Using Stepwise Optimization Framework ". Возможно, кто-то захочет сделать перевод этой статьи? В ней, кстати, при переходе от последовательной генерации случайных чисел к «потоковому» фиксируется хороший прирост скорости:
Changing the RNG from C++ TR1 to Intel MKL not only delivered 5.53X improvement on the original code, it also enabled the compiler to do inline function calls and vectorized the code in the inner loop...
Перекомпиляция под архитектуру MIC обеспечивает ещё больший прирост по сравнению с двумя 8-ядерными Ксеонами Sandy Bridge @2.6 Гц:
The same program that took 44 seconds to complete on the host system, now completes in just less than 5 seconds, an instant 8.82X improvement.
Правда, из статьи я так и не понял, какую именно модель Phi использовал автор, видимо, 7120X (16GB, 1.238 GHz, 61 core) за $4129. Пора, наверное, брать Phi, кто ещё не успел, не дожидаясь Knights Landing, ребята ;-) Интересно будет попробовать в бою и «сказать своё фи».
В любом случае, нас сейчас интересует подход MKL к генерации блоков случайных чисел. Сначала вызовом vslNewStream создаётся генератор потока случайных числе с заданными характеристиками, например, привлекательно звучит название базового генератора VSL_BRNG_SFMT19937 (A SIMD-oriented Fast Mersenne Twister pseudorandom number generator). Затем этот генератор можно использовать (причём из разных потоков) для непосредственного заполнения выделенного блока памяти последовательность случайных чисел, обладающих нужным распределением и типом данных. Вызов, например, virnguniform (method, stream, n, r, a, b) заполняет массив r на основе генератора stream n равномерно распределёнными случайными значениями от a до b. После работы поток следует удалить вызовом vslDeleteStream.
Есть важные аспекты при параллельном использовании/генерации случайных чисел, обусловленные конечностью периода псевдослучайных генераторов, прекрасно расписанные нашими соотечественниками Сергеем Майдановым и Андреем Нарайкиным в статье "Making the Monte Carlo Approach Even Easier and Faster". Для любознательных привожу отрывок (без перевода):
Due to that trend and the computational cost of Monte Carlo methods, the question of their use in parallel environments is an important one. Many Monte Carlo methods allow quite natural parallelization. In this case, the choice of random number generators should be made according to their ability to work in parallel. One important aspect that should be considered is
independence between sequences generated in parallel. There are a number of methods to generate independent sequences. We consider three of them below.The first is simultaneous use of several random-number generators, parameters of which are chosen so that sequences produced by different generators are independent in some sense (for example, in terms of the spectral test). The other two use splitting of the original pseudorandom sequence into k non-overlapping subsequences, where k is the number of threads/processes so that different threads/processes use random numbers from the corresponding subsequence only.
One of them is known as the block-splitting or skip-ahead method. In this case, the original sequence is being split into k non-overlapping blocks of subsequent elements, and each block maps to a corresponding subsequence. The other method is known as the leapfrog method. This differs from the skip-ahead method in that the leapfrog method splits the original sequence x1, x2,x3, … into k subsequences so that the first subsequence generates random numbers x1, xk+1, x2k+1, x3k+1, …, the second generates random numbers x2, xk+2, x2k+2, x3k+2, …, and, finally, the kth generates random numbers xk, x2k, x3k, …
There are advantages and disadvantages to these approaches in various parallel Monte Carlo methods. Using the first method, the maximal number of independent streams is limited by the number of suitable parameters chosen. With the skip-ahead method, a high correlation between blocks is possible, although the randomness properties of subsequences themselves are the same as with the original sequence.
When the leapfrog method is used, randomness properties of subsequences degrade dramatically when the number of subsequences increases. Finally, for some generators, the implementation of the latter two methods could be as inefficient as generation of the whole original sequence to pick out a required subsequence. Hence, the most appropriate generators should be selected while taking into account all these considerations.
Нужно сказать, что в своей библиотеке MKL умные парни из Intel реализовали и skip-ahead, и leapfrog подходы к разбиению на подпоследовательности. Однако в рамках моей задачи, хотя и используется несколько потоков от одного базового генератора, нет требований по абсолютной статистической независимости потоков друг от друга, поскольку они используются независимо и на разных подгруппах данных. Поэтому основная техническая сложность для меня — подключение заточенных под C и Fortran интерфейсов MKL к проекту на Visual Basic. На помощь приходит статья Use Intel® MKL from Microsoft* Office Excel, в которой расписано, как сгенерировать обёрточную dll силами самой MKL, и даны примеры вызовов функций MLK (конвенция stdcall, разумеется) из Excel.
Концепция такова: можно попросить MKL Builder автоматически сгенерировать dll, экспортирующую произвольный список нужных Вам функций MKL, с заданной конвенцией вызова. Конечно, можно создать свой собственный проект dll в VC++, но мне показалось, что автоматическая генерация — более универсальный способ, да ещё и поможет быстрее получить требуемый файл. Быстро выяснилось, что статья устарела, и даже тестовая команда nmake libia32 threading=sequential interface=stdcall export=user_example_list name=mkl_custom сбоила сразу по многим причинам. Одну из проблем решил добавлением файла msvcrt.lib в папку mkltoolsbuilder. Отдельно доставили удовольствие множественные сообщения unresolved external symbol ___security_cookie.
Пользуясь рекомендациями Microsoft и убив на поиски полдня, пришёл к выводу, что нужно модифицировать файл makefile в папке builder, приписав в командную строку линкера после out:$(CB_NAME).dll параметр bufferoverflowU.lib. Ну, тогда уж механизм со скрипом заработал, и, выдав nmake libia32 threading=parallel interface=stdcall export=vml_vsl_stdcall_example_list name=mkl_custom, я стал счастливым обладателем 38-мегабайтной динамической библиотеки, содержащей все функции Vector Statistical Library. Поборов в себе здоровую капиталистическую жадность, я оставил в библиотеке лишь 3 нужные функции, после чего размер снизился до нескольких мегабайт.
Теперь о вызывающем коде. В VB6 и VBA нужные функции следует объявлять так:
Public Declare Function vslNewStream Lib "mkl_custom.dll"
(ByRef StreamPtr As Long, ByVal brng As Long, ByVal seed As Long) As Long
Public Declare Function vslDeleteStream Lib "mkl_custom.dll"
(ByRef StreamPtr As Long) As Long
Public Declare Function viRngUniform Lib "mkl_custom.dll"
(ByVal method As Long, ByVal StreamPtr As Long, ByVal n As Long, ByRef Data As Long,
ByVal A As Long, ByVal B As Long) As Long
В качестве теста попробуем 100 раз сгенерировать массив из 5 млн случайных целочисленных индексов от 1 до 10000 (не включая 10000), используя базовый генератор VSL_BRNG_MCG31 (A 31-bit multiplicative congruential generator):
Dim i&, theStream&, v&(), res&, t!
Dim theStream As Long
ReDim v(1 To 5000000):
res = vslNewStream(theStream, VSL_BRNG_MCG31, 777)
t = Timer
For i = 1 To 100
res = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, theStream, 5000000, v(1), 1, 10000)
Next i
Debug.Print Timer - t
res = vslDeleteStream(theStream)
Время выполнения — 1.19 сек. Кстати, умеющий использовать SSE генератор VSL_BRNG_SFMT19937 с этой же задачей справляется за 0.99 сек. Ну что же, добавляем вызовы нужных операций в основной код, и мониторим загрузку процессора:
Бинго, полная загрузка! Была мысль сделать GPGPU-реализацию с использованием пакета cuRAND, на чьи преимущества недвусмысленно указывает Nvidia:
Blazing Fast cuRAND Performance compared to Intel MKL
Однако начальная цель была достигнута с помощью MKL, а на тестирование Cuda-версии пока нет свободного времени. :-)
Подведём итог: мы узнали о некоторых проблемах при параллельной генерации случайных чисел, и получили практические рекомендации по их решению с помошью Intel MKL. Использовать эту мощную библиотеку можно из любого языка современной Visual Studio — C++, C#, VB, а также Visual Fortran, и, с помощью озвученной методики, из любой среды, способной вызывать внешние dll, например, VB6, VBA, LabView, Matlab и множество других.
Спасибо за внимание!
Автор: fingoldo
