March 14, 2014 at 00:53Quickly generate random number arrays for simulation, statistical estimation and re-sampling
- Tutorial
Simulation using Monte Carlo methods nowadays is used in almost all areas of operational activity where multiple decisions are required based on the analysis of data received from the outside world. At the same time, the quality, performance and accessibility of random number generators, which are used to give the abstract method the features of a real problem solved by a specialist, begins to play an important role. As I recently found out, this issue begins to play a decisive role in the transition to parallel programming ... You also encountered this problem and want to know how to quickly get random number arrays with the necessary distribution on Windows?
/ * The article uses a significant number of citations in English, I assume a sufficient level of language training for the reader * /
One of my duty service applications runs simulation models on a long sequence of factors from a database. Without going into details, a task arose recently of studying the probabilistic characteristics of a repeating binary (two-source) process with a strong element of randomness from the history of its actual implementations in order to reduce uncertainty in predicting future outcomes. Implementations of the process are considered independent of each other and equally (normally) distributed, the process is considered stationary. The goal - the most accurate prediction of the results of N successive process implementations as possible - implies an estimate of the average and lowest expected frequencies of positive process implementations over a period of length N.
Here you can use the standard statistical approach with the calculation of the confidence level at which the realizations of a random variable do not deviate from its expectation by more than k values of its standard deviation. This approach is based on the idea of frequency convergence to probabilities ( Bernoulli theorem ). Of course, it is very tempting to calculate the average values based on several hundred thousand process implementations recorded in the database, but since I have a restriction on the length of the subsample N (of the order of several tens), the random nature of the process on a limited sample will surely bring me surprises in the form of a deviation from average. That is why I need to know with a given degree of certainty what maximum possible deviations from the average expected frequency should be "laid down" for a given N.
I decided to get the required estimates by the method of statistical bootstrapping (bootstrapping) - multiple generation of samples from N realizations using the Monte Carlo method from the main set of examples. The essence of the method is to form a sufficiently large number of pseudo-samples from the available sample, consisting of random combinations of the original set of elements (as a result, some initial elements may occur several times in one pseudo-sample, while others may be absent altogether), and for each pseudo-sample obtained determine the values of the analyzed statistical characteristics in order to study their scatter.
According to preliminary estimates, I will need to generate ~ 10 million pseudo-samples for each of the thousands of parameter subgroups. Given the amount of data and other factors, we are talking about tens of hours of server computing time. In Visual Basic 6, a prototype application simulator was written and tested. However, after compiling the multi-threaded combat version, it suddenly turned out that when the simulation was launched in several threads, the load per CPU core was only 5% instead of the expected 100%.
Basic Hotspots Analysis in the Intel Vtune Amplifier profiler showed nothing unusual. I was sure that the simulator should fully load all threads with work, since the code did not contain anything other than arithmetic operations and memory access by random indices:
Is this code so tied to memory bandwidth ?! And why then do we need multi-megabyte processor caches, why can't they cope ?! The problem baffled me for a few days. Belatedly launched Concurrency Analysis revealed a huge drawdown somewhere in the bowels of the VB system library, but could not show me where these costly system functions are called from:
Finally, experimentally, I found out that the problem was with the Rnd operator! Apparently, parallel use of Rnd is not provided by developers; perhaps the Microsoft implementation uses some kind of internal pending locks. In fact, parallelizing Rnd instead of speeding up the application slowed down THOUSANDS! I think a similar behavior will be typical for products from the MS Office package that use VBA. Never before has the issue of quickly generating random integer indices been so relevant to me! I decided to look for ready-made libraries for generating arrays of random numbers.
The choice fell on the Intel® Math Kernel Library with its Statistical Functions package .
MKL is available in test mode. Currently (version 11.02), 11 basic random variable generators, 12 continuous, 11 discrete output distributions are available. What turned out to be very important for me, the developer gives a guarantee of the thread safety of the library, moreover, it points out the advantages of its product when used in multi-threaded mode ... In the MKL documentation, my attention was drawn to the following section:
Working examples of using random numbers from MKL VSL for simulation using the Xeon Phi coprocessor architecture, for example, Shuo Li's Case Study: Achieving High Performance on Monte Carlo European Option Using Stepwise Optimization Framework , are already available. ". Perhaps someone wants to translate this article? In it, by the way, when switching from sequential generation of random numbers to" streaming ", a good speed increase is recorded:
Recompiling for the MIC architecture provides even greater growth compared to the two 8-core Sandy Bridge @ 2.6 Hz Xeons:
However, from the article I still did not understand which model Phi the author used, apparently, 7120X (16GB, 1.238 GHz, 61 core) for $ 4129. It's time, probably, to take Phi, who have not had time before, not waiting for Knights Landing, guys ;-) It will be interesting to try in battle and "say your phi."
In any case, we are now interested in MKL's approach to generating random number blocks. First, by calling vslNewStream, a random number stream generator with specified characteristics is created, for example, the name of the base generator VSL_BRNG_SFMT19937 (A SIMD-oriented Fast Mersenne Twister pseudorandom number generator) sounds attractive. Then this generator can be used (and from different streams) to directly fill the allocated memory block with a sequence of random numbers with the desired distribution and data type. A call, for example, virnguniform (method, stream, n, r, a, b) fills the array r based on the stream n generator with uniformly distributed random values from a to b. After work, the stream should be deleted by calling vslDeleteStream.
There are important aspects in the parallel use / generation of random numbers, due to the finite period of the pseudo-random generators, beautifully painted by our compatriots Sergei Maydanov and Andrei Naraikin in the article " Making the Monte Carlo Approach Even Easier and Faster ". For the curious, I quote an excerpt (without translation):
It must be said that in their MKL library, smart guys from Intel implemented both skip-ahead and leapfrog approaches to break up subsequences. However, in the framework of my task, although several streams from one basic generator are used, there are no requirements for the absolute statistical independence of streams from each other, since they are used independently and on different subgroups of data. Therefore, the main technical difficulty for me is connecting the MKL interfaces sharpened by C and Fortran to the Visual Basic project. The article Use Intel® MKL from Microsoft * Office Excel comes to the rescue , which describes how to generate a wrapper dll using MKL itself, and gives examples of calls to MLK functions (stdcall convention, of course) from Excel.
The concept is this: you can ask MKL Builder to automatically generate a dll that exports an arbitrary list of MKL functions you need with a given calling convention. Of course, you can create your own dll project in VC ++, but it seemed to me that automatic generation is a more universal way, and even helps to get the required file faster. It quickly became clear that the article was out of date, and even the test command nmake libia32 threading = sequential interface = stdcall export = user_example_list name = mkl_custom crashed immediately for many reasons. I solved one of the problems by adding the msvcrt.lib file to the mkl \ tools \ builder \ folder. Separately enjoyed the multiple messages unresolved external symbol ___security_cookie .
Using the recommendationsMicrosoft and having killed in search of half a day, I came to the conclusion that you need to modify the makefile in the builder folder by assigning bufferoverflowU.lib parameter to the linker command line after out: $ (CB_NAME) .dll. Well, then the mechanism with a scratch worked, and by issuing nmake libia32 threading = parallel interface = stdcall export = vml_vsl_stdcall_example_list name = mkl_custom , I became the lucky owner of a 38-megabyte dynamic library containing all the functions of the Vector Statistical Library. Overcoming the healthy capitalist greed in myself, I left only 3 necessary functions in the library, after which the size decreased to a few megabytes.
Now about the calling code. In VB6 and VBA, the required functions should be declared like this:
As a test, let's try to generate an array of 5 million random integer indices from 1 to 10,000 (not including 10,000) 100 times using the base generator VSL_BRNG_MCG31 (A 31-bit multiplicative congruential generator):
Runtime - 1.19 seconds. By the way, the VSL_BRNG_SFMT19937 generator that knows how to use SSE can do the same task in 0.99 seconds. Well, we add the calls of the necessary operations to the main code, and monitor the processor load:
Bingo, full load! There was an idea to make a GPGPU implementation using the cuRAND package, whose advantages are clearly indicated by Nvidia:
Blazing Fast cuRAND Performance compared to Intel MKL
However, the initial goal was achieved using MKL, and there is no free time to test the Cuda version yet. :-)
To summarize: we learned about some problems in the parallel generation of random numbers, and received practical recommendations for solving them using Intel MKL. You can use this powerful library from any language of modern Visual Studio - C ++, C #, VB, as well as Visual Fortran, and, using the voiced technique, from any environment that can call external dlls, for example, VB6, VBA, LabView, Matlab and many others.
Thanks for attention!
/ * The article uses a significant number of citations in English, I assume a sufficient level of language training for the reader * /
One of my duty service applications runs simulation models on a long sequence of factors from a database. Without going into details, a task arose recently of studying the probabilistic characteristics of a repeating binary (two-source) process with a strong element of randomness from the history of its actual implementations in order to reduce uncertainty in predicting future outcomes. Implementations of the process are considered independent of each other and equally (normally) distributed, the process is considered stationary. The goal - the most accurate prediction of the results of N successive process implementations as possible - implies an estimate of the average and lowest expected frequencies of positive process implementations over a period of length N.
Here you can use the standard statistical approach with the calculation of the confidence level at which the realizations of a random variable do not deviate from its expectation by more than k values of its standard deviation. This approach is based on the idea of frequency convergence to probabilities ( Bernoulli theorem ). Of course, it is very tempting to calculate the average values based on several hundred thousand process implementations recorded in the database, but since I have a restriction on the length of the subsample N (of the order of several tens), the random nature of the process on a limited sample will surely bring me surprises in the form of a deviation from average. That is why I need to know with a given degree of certainty what maximum possible deviations from the average expected frequency should be "laid down" for a given N.
I decided to get the required estimates by the method of statistical bootstrapping (bootstrapping) - multiple generation of samples from N realizations using the Monte Carlo method from the main set of examples. The essence of the method is to form a sufficiently large number of pseudo-samples from the available sample, consisting of random combinations of the original set of elements (as a result, some initial elements may occur several times in one pseudo-sample, while others may be absent altogether), and for each pseudo-sample obtained determine the values of the analyzed statistical characteristics in order to study their scatter.
According to preliminary estimates, I will need to generate ~ 10 million pseudo-samples for each of the thousands of parameter subgroups. Given the amount of data and other factors, we are talking about tens of hours of server computing time. In Visual Basic 6, a prototype application simulator was written and tested. However, after compiling the multi-threaded combat version, it suddenly turned out that when the simulation was launched in several threads, the load per CPU core was only 5% instead of the expected 100%.
Basic Hotspots Analysis in the Intel Vtune Amplifier profiler showed nothing unusual. I was sure that the simulator should fully load all threads with work, since the code did not contain anything other than arithmetic operations and memory access by random indices:
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
Is this code so tied to memory bandwidth ?! And why then do we need multi-megabyte processor caches, why can't they cope ?! The problem baffled me for a few days. Belatedly launched Concurrency Analysis revealed a huge drawdown somewhere in the bowels of the VB system library, but could not show me where these costly system functions are called from:
Finally, experimentally, I found out that the problem was with the Rnd operator! Apparently, parallel use of Rnd is not provided by developers; perhaps the Microsoft implementation uses some kind of internal pending locks. In fact, parallelizing Rnd instead of speeding up the application slowed down THOUSANDS! I think a similar behavior will be typical for products from the MS Office package that use VBA. Never before has the issue of quickly generating random integer indices been so relevant to me! I decided to look for ready-made libraries for generating arrays of random numbers.
The choice fell on the Intel® Math Kernel Library with its Statistical Functions package .
MKL is available in test mode. Currently (version 11.02), 11 basic random variable generators, 12 continuous, 11 discrete output distributions are available. What turned out to be very important for me, the developer gives a guarantee of the thread safety of the library, moreover, it points out the advantages of its product when used in multi-threaded mode ... In the MKL documentation, my attention was drawn to the following section:
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].I have long wanted to try this “iron” feature, but, unfortunately, it turned out that my processors are not equipped with a hardware random number generator. In addition to the above, an important feature of the MKL package is the paradigm of block generation of random number streams, in contrast to the usual practice of obtaining one random value per operator call in most programming languages.
Working examples of using random numbers from MKL VSL for simulation using the Xeon Phi coprocessor architecture, for example, Shuo Li's Case Study: Achieving High Performance on Monte Carlo European Option Using Stepwise Optimization Framework , are already available. ". Perhaps someone wants to translate this article? In it, by the way, when switching from sequential generation of random numbers to" streaming ", a good speed increase is recorded:
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 ...
Recompiling for the MIC architecture provides even greater growth compared to the two 8-core Sandy Bridge @ 2.6 Hz Xeons:
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.
However, from the article I still did not understand which model Phi the author used, apparently, 7120X (16GB, 1.238 GHz, 61 core) for $ 4129. It's time, probably, to take Phi, who have not had time before, not waiting for Knights Landing, guys ;-) It will be interesting to try in battle and "say your phi."
In any case, we are now interested in MKL's approach to generating random number blocks. First, by calling vslNewStream, a random number stream generator with specified characteristics is created, for example, the name of the base generator VSL_BRNG_SFMT19937 (A SIMD-oriented Fast Mersenne Twister pseudorandom number generator) sounds attractive. Then this generator can be used (and from different streams) to directly fill the allocated memory block with a sequence of random numbers with the desired distribution and data type. A call, for example, virnguniform (method, stream, n, r, a, b) fills the array r based on the stream n generator with uniformly distributed random values from a to b. After work, the stream should be deleted by calling vslDeleteStream.
There are important aspects in the parallel use / generation of random numbers, due to the finite period of the pseudo-random generators, beautifully painted by our compatriots Sergei Maydanov and Andrei Naraikin in the article " Making the Monte Carlo Approach Even Easier and Faster ". For the curious, I quote an excerpt (without translation):
Hidden text
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.
It must be said that in their MKL library, smart guys from Intel implemented both skip-ahead and leapfrog approaches to break up subsequences. However, in the framework of my task, although several streams from one basic generator are used, there are no requirements for the absolute statistical independence of streams from each other, since they are used independently and on different subgroups of data. Therefore, the main technical difficulty for me is connecting the MKL interfaces sharpened by C and Fortran to the Visual Basic project. The article Use Intel® MKL from Microsoft * Office Excel comes to the rescue , which describes how to generate a wrapper dll using MKL itself, and gives examples of calls to MLK functions (stdcall convention, of course) from Excel.
The concept is this: you can ask MKL Builder to automatically generate a dll that exports an arbitrary list of MKL functions you need with a given calling convention. Of course, you can create your own dll project in VC ++, but it seemed to me that automatic generation is a more universal way, and even helps to get the required file faster. It quickly became clear that the article was out of date, and even the test command nmake libia32 threading = sequential interface = stdcall export = user_example_list name = mkl_custom crashed immediately for many reasons. I solved one of the problems by adding the msvcrt.lib file to the mkl \ tools \ builder \ folder. Separately enjoyed the multiple messages unresolved external symbol ___security_cookie .
Using the recommendationsMicrosoft and having killed in search of half a day, I came to the conclusion that you need to modify the makefile in the builder folder by assigning bufferoverflowU.lib parameter to the linker command line after out: $ (CB_NAME) .dll. Well, then the mechanism with a scratch worked, and by issuing nmake libia32 threading = parallel interface = stdcall export = vml_vsl_stdcall_example_list name = mkl_custom , I became the lucky owner of a 38-megabyte dynamic library containing all the functions of the Vector Statistical Library. Overcoming the healthy capitalist greed in myself, I left only 3 necessary functions in the library, after which the size decreased to a few megabytes.
Now about the calling code. In VB6 and VBA, the required functions should be declared like this:
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
As a test, let's try to generate an array of 5 million random integer indices from 1 to 10,000 (not including 10,000) 100 times using the base generator 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)
Runtime - 1.19 seconds. By the way, the VSL_BRNG_SFMT19937 generator that knows how to use SSE can do the same task in 0.99 seconds. Well, we add the calls of the necessary operations to the main code, and monitor the processor load:
Bingo, full load! There was an idea to make a GPGPU implementation using the cuRAND package, whose advantages are clearly indicated by Nvidia:
Blazing Fast cuRAND Performance compared to Intel MKL
However, the initial goal was achieved using MKL, and there is no free time to test the Cuda version yet. :-)
To summarize: we learned about some problems in the parallel generation of random numbers, and received practical recommendations for solving them using Intel MKL. You can use this powerful library from any language of modern Visual Studio - C ++, C #, VB, as well as Visual Fortran, and, using the voiced technique, from any environment that can call external dlls, for example, VB6, VBA, LabView, Matlab and many others.
Thanks for attention!
