Can Mathematica's random number generation be improved?

Question

Performance on Random Number generation is intolerable. Whenever you need many deviates from the same distribution but with different parameters it takes forever.

• Here is a thread on Poisson deviates being slow. The comments and answers are all aimed at Poisson Deviates and do not address the fundamental shortcoming Mathematica has generating random numbers.

• Here I asked why MultiNomial Deviates were so slow, and didn't get an answer, although one person suggested there was no solution.

• Elsewhere I asked how to generate BinomialDeviates and got a helpful Binomial-specific answer.

The problem is I don't want to spend hours trying to work around Mathematica's intolerable performance each time I need some deviates from a standard distribution.

---

This is essentially a suggestion that WRI needs to add function calls for people who need large numbers of random deviates in which all the parameters are different. Below I show the syntax for such a call for binomial variates in R, Octave, and scipy.

R:

num<-10^5    
n<-sample(1:10,num,TRUE)
p<-runif(num)

ptm <- proc.time()
out=rbinom(num,n,p)
proc.time()-ptm

Here is the same thing in Octave:

num=10^7;
n=randi([1,10],num,1);
p=rand(num,1);

t=cputime;
out=binornd(n,p);
cputime-t

and scipy:

num=10**7
n=stats.randint.rvs(1,10,size=num)
p=stats.uniform.rvs(0,1,num)
out=stats.binom.rvs(n,p)

NOTE: The point is not to generate the distribution of deviates that are obtained for uniformly distributed p and n. The point is to generate a bunch of deviates when each deviate is drawn for a different parameter value. The manner in which the parameters are obtained should be considered as unknown. I used uniform for p and n because they had to be set to something, not because I am interested in the resulting distribution.

As you can see the calls in R, Octave, aand Python (scipy) where the actual work desired is performed is a single line of self-evident code in each case. The performance of all these languages blows Mathematica away. Mathematica has Solve and NSolve, DSolve & NDSolve, Integrate & NIntegrate, Sum and NSum, Product & NProduct... Something similar (to the numerical versions) should be implemented for random number generation.

Below there are 3 answers. blochwave's solution has the best performance and makes calls to C++ routines which the user must write. The routine poissonvariate.cpp has 48 lines. In addition there are two lines of Mathematica which must be written as well. The 48 lines of cpp could obviously be copied from the working poissonvariate and then morphed into the code needed for the desired generator but one must still perform substantially more work to implement this solution for new generators than one must in R, Octave, and python. My solution was to use Rlink and generate the numbers in R. THis gave a performance improvement of about 100 over native Mathematica which for my purposes was adequate. The amount of code required to implement the R solution is 2 + n where n is the number of distinct random variate types you need. (BY "type" I mean, exponential, binomial, poisson, multinomial, beta, chisquare, etc. etc. etc.) . I believe that most users value both code simplicity and performance. The 3rd answer did not provide a solution to the problem that was posed.

---

Below I compare computing 10,000 random deviates in Mathematica against 10,000,000 in C. I draw from 5 distributions in Mathematica and C: Binomial, Exponential, Chi-Squared, Poisson and Gamma. In C I generated about 50 million deviates in about 10 seconds. In Mathematica I generated about 50,000 deviates in about 8 seconds. If I generate 10,000 deviates of each distribution in C the code runs in about 0.03 seconds which is roughly 250 times faster than Mathematica. Is there really no way to improve this?

Here is the comparison between Mathematica and C. The upshot is that there is no comparison. C is about 100-1000 times faster.

num = 10000

ints  = RandomInteger[{1, 10}, num];
reals = RandomReal[{0, 1}, num];
mus   = RandomReal[{0, 1}, num];

mf    = Thread[BinomialDistribution[ints, reals]];    
gf    = Thread[GammaDistribution[mus, reals]];

RandomVariate /@ mf; // Timing     
RandomVariate /@ gf; // Timing

RandomVariate[PoissonDistribution[#]] & /@ mus; // Timing
RandomVariate[ChiSquareDistribution[#]] & /@ ints; // Timing
RandomVariate[ExponentialDistribution[#]] & /@ mus; // Timing

The 10,000 Binomial deviates take 6 seconds (on a MacBook Air) to generate. The Poisson deviates take over 1 second.

The C code below uses the library ranlib, which is documented and virtually bullet-proof as far as I know. You can get ranlib.c here - it's the 7th library down.

To run the code below save it to a file called ran_test.c. Get the source code from the library and put linpack.c, ranlib.c, com.c and ranlib.h into the same directory. On a Unix machine (e.g. Linux or Mac) that has GCC installed type the following to compile.

gcc -O2 ran_test.c linpack.c ranlib.c com.c -o tst  

To time and execute the code enter: time ./tst on the command line and hit Return. On my MacBook Air the C generates 10 million deviates from each of the following distributions: Gamma, Binomial, Chi-Squared, Exponential and Poisson. The parameters for the distributions are different on each call. The code runs in about 10 seconds. Thus it is roughly 1,000 times faster than Mathematica.

#include <stdio.h>
#include <stdlib.h>
#include "ranlib.h"

#define MILLION 1000000
#define TEN_MILLION 10000000

main(){
  long int i,j,k,dev,count, ix,n;
  long iseed1=100,iseed2=1000;
  long int ncat, rint;
  float p, x,pp;
  float mu;

   setall(iseed1,iseed2);

  rint = ignuin(1,TEN_MILLION);
  printf("random int=%ld\n",rint);

  for(k=0;k<TEN_MILLION;k++){
    p=genunf(0,1);
    ncat=ignuin(1,10);
    n=ignbin(ncat, p);
    mu=genunf(0,10000);
    x=gengam(genunf(0,100),genunf(0,10000));

    /* print a few of the deviates */

    if(k%rint==0)printf("gamma deviate=%g\n",x);
    if(k%rint==0)printf("binomial deviate=%ld\n",n);
    n=ignpoi(mu);
    if(k%rint==0)printf("poisson variate n=%ld\n",n);
    p=genexp(mu);
    if(k%rint==0)printf("exp variate =%g\n",p);
    p=genchi(mu);
    if(k%rint==0)printf("chi square =%g\n",p);
    if(k%rint==0)printf("\n");
  } 
}

Answer 1

Update 06/10/2016

Now I'm using Mathematica 11.0, things seem to have improved in some cases:

list1 = RandomInteger[{1, 10}, {10^5}];
list2 = RandomReal[{0, 1}, {10^5}];
MapThread[
    RandomVariate[BinomialDistribution[#1, #2]] &, {list1, 
     list2}]; // AbsoluteTiming

(* 38.59 seconds on MMA10 *)
(*  2.94 seconds on MMA11, same machine *)

Still not as fast as my compiled version (0.02 seconds!), but a definite improvement. Looking at the other distributions:

list = RandomReal[{0, 100}, {10^5}];
RandomVariate[PoissonDistribution[#]] & /@ list; // AbsoluteTiming

(* 13.16 seconds on MMA10 *)
(*  3.24 seconds on MMA11, same machine *)

list = RandomReal[{0, 100}, {10^6}];
RandomVariate[ExponentialDistribution[#]] & /@ list; // AbsoluteTiming

(*  5.06 seconds on MMA10 *)
(* 12.83 seconds on MMA11, same machine *)

list = RandomReal[{0, 100}, {10^6}];
RandomVariate[ChiSquareDistribution[#]] & /@ list; // AbsoluteTiming

(* 11.10 seconds on MMA10 *)
(*  6.32 seconds on MMA11, same machine *)

Oh dear, not good news for ExponentialDistribution[]!

---

Update 20/02/2015

I've now adapted my code to work with a binomial distribution, since that seemed to be the most problematic case for the OP. The relevant C++ code is at the bottom of this answer.

Needs["CCompilerDriver`"]

BinomialVariateLib = 
  CreateLibrary[{ToString[NotebookDirectory[]] <> 
     "binomialvariate.cpp"}, "BinomialVariateLib", "Debug" -> False, 
   "CompileOptions" -> "-std=c++11 -O3"];
BinomialVariate = 
  LibraryFunctionLoad[BinomialVariateLib, 
   "BinomialVariate", {{Real, 1}, {Real, 1}}, {Real, 1}];

list1 = RandomInteger[{1, 10}, {10^5}];
list2 = RandomReal[{0, 1}, {10^5}];

l1 = BinomialVariate[list1, list2]; // AbsoluteTiming
l2 = MapThread[
    RandomVariate[BinomialDistribution[#1, #2]] &, {list1, 
     list2}]; // AbsoluteTiming

(* 0.020 seconds *)
(* 38.59 seconds *)

I make that a speed-up of 2000x, with no multi-thread/parallel computing invoked. Trying with a larger list of $10^{6}$ as the OP does in another answer, my method takes 0.16 seconds (I didn't bother timing the MMA version, apparently it takes minutes). This appears to be a similar performance to the R code given by the OP in a separate answer.

It's even better on Windows 8.1 with the Visual Studio compiler, I get a 2700x speed-up. Perhaps the Intel compiler would be even quicker. Unfortunately the OP couldn't get my code to work on a Mac. I guess it has something to do with the use of std::random and std::chrono, but these could be overcome by using a non-C++11 solution. Alternatively, I believe it is possible to install GCC on Mac OS? The above example is compiled with GCC.

Note: @ciao's answer using a parametric mixture distribution is very quick, even more so than my unoptimized C++ below. It takes about 0.13 seconds on my machine to generate the distribution and then draw $10^{6}$ samples from it. At a guess, putting together a similar method in C++ will be faster still, but perhaps not worth the effort. It may also not be appropriate for the OP's case if the distribution needs to be changed each time.

---

Original answer

One example where drawing random numbers from distributions with different parameters is for adding Poisson noise to an image, where the level of noise depends on the underlying signal.

I encountered this very problem, and because only certain distributions can be compiled in Mathematica, I asked a question about it: Speed up adding Poisson noise to an image.

Eventually I came up with a solution using C++ to deal with the distributions, and pass the data back to MMA using LibraryLink. Adapting that solution, the following improvement is possible for the Poisson distribution. The relevant C++ code is at the bottom of this answer.

Needs["CCompilerDriver`"]

PoissonVariateLib = 
  CreateLibrary[{ToString[NotebookDirectory[]] <> 
     "poissonvariate.cpp"}, "PoissonVariateLib", "Debug" -> False, 
   "CompileOptions" -> "-std=c++11 -O3"];
PoissonVariate = 
  LibraryFunctionLoad[PoissonVariateLib, 
   "PoissonVariate", {{Real, 1}}, {Real, 1}];

list = RandomReal[{0, 100}, {10^5}];

PoissonVariate@list; // AbsoluteTiming
RandomVariate[PoissonDistribution[#]] & /@ list; // AbsoluteTiming

(* 0.08 seconds *)
(* 13.16 seconds *)

This is 150x speed-up over MMA! I've not looked into optimizing the C++ at all, and nor am I making much use of multiple CPU cores by parallelizing any of my code at either the C++ or MMA level.

However, this answer should demonstrate that if time is a critical factor for you, then LibraryLink or MathLink to C/C++ is a pretty good option to consider.

The same process can be done for other distributions in your example, by using what's available in std::random in C++. First the exponential distribution (remember to adapt the C++ code!):

(* Trying it with a bigger list this time *)
list = RandomReal[{0, 100}, {10^6}];

ExponentialVariate@list; // AbsoluteTiming
RandomVariate[ExponentialDistribution[#]] & /@ list; // AbsoluteTiming

(* 0.22 seconds *)
(* 5.06 seconds *)

Another impressive speed-up. Changing the code again for the Chi-Squared distribution, we find another good speed-up:

(* Trying it with a bigger list this time *)
list = RandomReal[{0, 100}, {10^6}];

ChiSquaredVariate@list; // AbsoluteTiming
RandomVariate[ChiSquareDistribution[#]] & /@ list; // AbsoluteTiming

(* 0.11 seconds *)
(* 11.1 seconds *)

---

C++ code

This code is for the file binomialvariate.cpp, as per my updated answer. The code for the Poisson, exponential and chi-squared distributions (and indeed many other distributions) should be straightforward to implement from this, or from my answer here.

#include "WolframLibrary.h"
#include "stdio.h"
#include "stdlib.h"
#include <random>
#include <chrono>

EXTERN_C DLLEXPORT mint WolframLibrary_getVersion(){return WolframLibraryVersion;}
EXTERN_C DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {return 0;}
EXTERN_C DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {}

EXTERN_C DLLEXPORT int BinomialVariate(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

    int err; // error code

    MTensor m1; // input tensor 2
    MTensor m2; // input tensor 2
    MTensor m3; // output tensor

    mint const* dims1; // dimensions of the tensor

    double* data1; // actual data of the input tensor
    double* data2;
    double* data3; // data for the output tensor

    mint i; // bean counters

    m1 = MArgument_getMTensor(Args[0]);
    dims1 = libData->MTensor_getDimensions(m1);
    m2 = MArgument_getMTensor(Args[1]);
    err = libData->MTensor_new(MType_Real, 1, dims1, &m3);
    data1 = libData->MTensor_getRealData(m1);
    data2 = libData->MTensor_getRealData(m2);
    data3 = libData->MTensor_getRealData(m3);

    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::mt19937 generator (seed); // RNG

    for(i=0;i<dims1[0];i++) {
        std::binomial_distribution<int> distribution(data1[i], data2[i]); 
        data3[i] = distribution(generator); 
    }

    MArgument_setMTensor(Res, m3);
    return LIBRARY_NO_ERROR;
}

Answer 2

This is a poor use of the random functions in Mathematica. As clearly stated in the documentation, generation of variates one at a time has significant overhead, and generating them en masse has significant benefits, particularly with statistical distributions:

For statistical distributions, the speed advantage of generating many numbers at once can be even greater. In addition to the efficiency benefit inherited from the uniform number generators used, many statistical distributions also benefit from vectorized evaluation of elementary and special functions. For instance, WeibullDistribution benefits from vector evaluations of the elementary functions Power, Times, and Log.

This is mentioned in a tech note here.

E.g., generating 100K Weibull RV has a speed difference of over 300X.

If you need such variates generated individually, or because they cannot be generated en masse, explore other options. Using a ParameterMixture distribution can speed things dramatically when the parameters themselves are RV (as in your examples), you can build custom generators within Mathematica for specific use cases/custom distributions that can exhibit marked performance benefits, or if you have external facilities that already generate them in the form and with the speed you desire, just call them from within Mathematica.

Per Blochwave's comment, an example of the technique I use for such things. This is the binomial case:

dist = ParameterMixtureDistribution[
   BinomialDistribution[a, b], {a \[Distributed] DiscreteUniformDistribution[{1, 10}], 
                                b \[Distributed] UniformDistribution[]}];

pdf = PDF[dist, #] & /@ Range[0, 10];

variates = RandomChoice[pdf -> Range[0, 10], 10^6]; // Timing

(* {0.202801,Null} *)

That's on an old netbook. The R example on the same netbook takes ~1.6 seconds...

Similar techniques can be used for continuous derived distributions, substituting inversion (or technique of choice) for the RandomChoice.

As with any language, one must think about the tools the language makes available, and how to use them efficiently (want the exact symbolic PDF from R to do such things? Oops, out of luck...). Any pair of languages lend themselves to arbitrary efficiency drag races where one might appear to be clearly superior, until one uses the tools intelligently. One does not expect the Ferrari to be an efficient garbage scow, every software tool has its place and its own idiosyncrasies and coding techniques that can be opaque to many.

Answer 3

Here is my solution. The first 3 lines of code were what I needed to make this work on my Mac. It is not in the documentation and you'd need something else on a PC.

Needs["RLink`"]
SetEnvironment["DYLD_LIBRARY_PATH" -> "/Library/Frameworks/R.framework/Resources/lib"];
InstallR["RHomeLocation" -> "/Library/Frameworks/R.framework/Resources"];
REvaluate["R.version.string"]

num = 10^6;
ints = RandomInteger[{1, 10}, num];
reals = RandomReal[{0, 1}, num];
binom = REvaluate["rbinom"];
pois = REvaluate["rpois"];
binom[num, ints, reals]; // Timing
pois[num, reals]; // Timing

I gather that comparing Mathematica to C is considered invalid although I don't understand why. Is it also invalid to compare the performance of Mathematica to R on the same problem? Are all comparisons of the performance of other languages to Mathematica invalid or are there some valid ones?

In the hope that R is considered a more appropriate comparison, The following R is faster than my compiled C by a factor of 5 or 10.

num=10^6; 
ptm <- proc.time(); 
out=rbinom(num,sample(1:10,num,TRUE),runif(num)); 
proc.time()-ptm

The above line of R runs in about 150 milliseconds. The equivalent line of Mathematica takes over 10 minutes.

num = 10^6; 
MapThread[
    RandomVariate[BinomialDistribution[#1, #2]] &,  
        {RandomInteger[{1, 10}, num], RandomReal[{0, 1}, num]}
   ]; // Timing

(* 629.910404, Null *)

The fact that the lousy performance the random number generation is documented is of little help. While blochwave's solution will prove useful it will require some effort to make it work on my system. It took about 5 minutes to read the R documentation and write the code.

Here is the Octave code:

$ cat rnumcall.m    
t=cputime;
num=10^7;
n=randi([1,5],num,1);
p=rand(num,1);
out=binornd(n,p);
cputime-t

octave rnumcall.m
ans =  2.8404

or with num=10^6,

$ octave rnumcall.m    
ans =  0.29141

which is of course faster than Mathematica but now only 2000x as fast rather than R's 4000x.

Answer 4

I don't know why Mathematica's RandomVariate is so slow. But I have an easy method to generate lots of random values for exponential distributions.

Let $\varrho$ be a propabiblity distribution with continuous probability density function $f$ (PDF). The idea is that the commulative distribution function $F$ (CDF) of $\varrho$ pushes $f$ forward to the uniform probability distribution on the unit interval. So, we can generate random values for $\varrho$ by transforming uniformly distributed pseudorandom numbers on the unit interval with the inverse of F. For that, we need only a single line of code:

getExponentialDistribution[λ_] := -Log[RandomReal[{0, 1}, Dimensions[λ]]]/λ

The following is to convince you that generated random numbers really follow the exponential distribution (and that the generation is reasonably fast):

n = 10^6;
a = getExponentialDistribution[ConstantArray[1., {n}]]; // AbsoluteTiming // First
EstimatedDistribution[a, ExponentialDistribution[x]]
Show[
 Histogram[a, {1/5}],
 Plot[Evaluate[n/5 PDF[ExponentialDistribution[1]][x]], {x, 0, 10}, 
  PlotRange -> All],
 PlotRange -> All
 ]

0.016884

ExponentialDistribution[0.999564]

And here is a comparison of this with RandomVariate when applied to exponential distributions with randomized parameters:

λ = RandomReal[{0, 100}, {n}];
RandomVariate[ExponentialDistribution[#]] & /@ λ; // AbsoluteTiming // First
getExponentialDistribution[λ]; // AbsoluteTiming // First

2.79397

0.016625