# Can Mathematica's random number generation be improved?

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.

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];

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");
}
}

• What's your point? Of course compiled C is faster than Mathematica. So call your C from it if that kind of speed is critical for this simplistic stuff. On the other hand, I'll wait patiently while you use your C code to evaluate complex probabilistic equations symbolically... – ciao Feb 18 '15 at 6:58
• I don't get the point either. If you want to do a comparison, get it down to 2 lines of code. Do you really need to define num, ints, reals, mus, f, mf, rvs, gf ... just to do a timing test? Talking of efficient, who has time to look through all that? – wolfies Feb 18 '15 at 10:57
• Yes - but your question is messy. A good question simplifies the problem down to 1 or 2 lines (if possible), and it is certainly possible to do so here. By contrast, your question appears lazy, because it just dumps lines and lines of code from whatever you were doing, without honing it down to the heart of the matter. 2 lines: no more. Further, your title asks if it can be FIXED, which infers that it is broken. But your question is not about quality, nor do you compare the quality of pseudorandom drawings generated, nor does your post suggest anything is broken. – wolfies Feb 18 '15 at 15:13
• If you feel hostility I think that might have to do with the fact that you chose to use a quite aggressive title, which also isn't in line with your actual question. You might realize that many readers would be seriously concerned about broken random number generators while they are neither surprised nor worried to hear that Mathematica is slow for this kind of task. They might (as I did) feel being tricked into reading your (lengthy) question in the first place... – Albert Retey Feb 20 '15 at 15:50
• @blochwave: That 200ms was on a netbook, 6x faster than the R example on same - but it's basic MMA thinking: Use ParametricMixture to get dist, PDF (or CDF depending on dist type), use PDF to drive RandomChoice (or CDF for inversion). Will add binomial example to my answer when time permits. – ciao Feb 20 '15 at 22:00

## 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}];
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 =
"BinomialVariate", {{Real, 1}, {Real, 1}}, {Real, 1}];

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

l1 = BinomialVariate[list1, list2]; // AbsoluteTiming
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.

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 =
"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;
}

• Very nice answer, +1 – ciao Feb 18 '15 at 8:44
• Notice that it allows you to define the C++ code inside Mathematica, rather than as a separate file. – dr.blochwave Feb 18 '15 at 13:53
• @JEP: The headers random and chrono were added to C++ in the C++11 standard. It might be that your C++ compiler doesn't support C++11, or that you have to give a special flag to compile C++11 code. – celtschk Feb 18 '15 at 21:22
• @JEP and others: I got it to work on OS X 10.10.2 with g++ installed with homebrew and adding "CompilerInstallation" -> "/usr/local/bin/g++" to CreateLibrary options. Timings are 0.012565 and 51.336809 for the compiled and native method, respectively. – shrx Feb 21 '15 at 8:56
• @blochwave May I suggest to change BinomialVariate = LibraryFunctionLoad[BinomialVariateLib, "BinomialVariate", {{Real, 1}, {Real, 1}}, {Integer, 1}]; and in the C++ code mint* data3; and  data3 = libData->MTensor_getIntegerData(m3); and data3[i] = (int) distribution(generator); ? Also I think this AbsoluteTiming[ l2 = MapThread[BinomialVariate[{#1, #2}] &, {list1, list2}];] is a better comparison. – Rolf Mertig Feb 3 '16 at 18:51

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.

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.

• An example would have been helpful. – JEP Feb 18 '15 at 13:07
• @jep Again, what's your point? Mathematica has no automatic detection of stupid coding and coders. If you prefer some other system, use it, complain on their site... – ciao Feb 19 '15 at 23:32
• @JEP You seem to be confused. This site is not run by or directly affiliated with the developers or publishers of the software Mathematica. This is not the place for such complaints, at least not in that tone. Please desist from making them. Constructive participation is expected here. – Mr.Wizard Feb 20 '15 at 19:44
• @JEP There are plenty of things in any system to complain about; I've certainly pointed out my share of problems. However this community should be about helping each other. We have had people come to this site who seem to be here to do nothing but tear down Mathematica and its users, and I think you struck a raw nerve. It is entirely valid to publicly wonder why Mathematica has certain failings but you must find a way to make it apparent that you are working toward solutions not simply complaining. – Mr.Wizard Feb 20 '15 at 23:32
• To @rasher and others I ask you to take a step back and consider how you might assist a user to better understand the software, accept its limitations, and participate constructively. Rejoinders like "What's your point? Of course..." sound confrontational even when they are true. Issues that are "of course" understood by experienced users may be new to others. – Mr.Wizard Feb 20 '15 at 23:37

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;
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.

• I think that if blochwave's answer is not that appealing to you, then you probably came to the wrong place to raise this concern. Although I agree with you that these functions seem to be unfeasibly slow (especially in comparison to R, which is not known for its speed), the fact is that there's really nothing we as users can do about it. So, if this still seriously concerns you even though workarounds are available, I'd advocate making the request to WRI. If enough people ask, they do try to spend the time to improve features. But if nobody asks, it probably won't get done. – Oleksandr R. Feb 19 '15 at 0:09
• Thanks. I intend to implement blochwave's solution and I truly appreciate that he took the time to share the code. But appealing? Of course it isn't appealing. With my "lazy" "brain dead" programming style I can accomplish the same "simplistic" task in 4 or 5 lines in both R and Octave. R & Octave (which are both free and open source) are of course both about 1,000 times faster than MMA. The code blochwave provided doesn't run on any of my 3 Macs. Based on the upvotes of rasher's "answer" a large number of users here here think MMA's lousy performance is a feature not a bug. – JEP Feb 19 '15 at 19:06
• Using Mathematica sanely, using the R example, a parametric mixture CDF and random choice generates the 10^6 variates in ~200 milliseconds. On an old friggin' netbook... So again, using the proper tool and using it properly matters. OTOH, if all you're after is raw generation of such variates, you can better R by orders of magnitude with other tools. – ciao Feb 20 '15 at 2:15
• @JEP You are completely right: performance matters and WRI managenent should take it much much more seriouly. My vote is more speed and less Twitter! – Rolf Mertig Feb 20 '15 at 8:12
• "It occurred to me that since R is 4,000 times faster than MMA while C++ is only 150 times faster that the thinking man, even one that is brain dead and lazy, would probably choose R over C++" - you aren't comparing like-with-like. You're comparing a Binomial Distribution in R with a Poisson Distribution in C++. You'll see from my update that the actual difference between R and C++ when comparing the same thing is negligible. – dr.blochwave Feb 20 '15 at 8:26

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[λ]]]/λ


This is you 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