Can Mathematica's random number generation be improved?

23

12

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

JEP

Posted 2015-02-18T06:25:31.647

Reputation: 1 238

11What'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 – 2015-02-18T06:58:17.417

1I 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 – 2015-02-18T10:57:49.587

1I was not aware that symbolic calculations have anything to do with random number generation? ?? – JEP – 2015-02-18T12:53:08.307

f should not have been there. Other than that ints, reals, and mus were all used as arguments to the distributions.

I am detecting hostility here which I do not understand. I don't understand the point of comments that don't address the question.

My point is that if you can generate random variates all at once the performance is comparable to C so it is false that "of course compiled C is faster than Mathematica." There are many simulation problems in which random numbers with different parameters are needed as in the examples I provided. – JEP – 2015-02-18T13:04:32.223

8Yes - 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 – 2015-02-18T15:13:26.533

3If 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 – 2015-02-20T15:50:56.817

1Sorry about that. I wasn't trying to trick anyone into reading anything. The first line in my post is: "Performance on Random Number generation is intolerable." I'm not sure why you kept reading if you don't care about performance. WRI can fix it or not. – JEP – 2015-02-20T17:10:50.620

Seeing as intelligent use of Mathematica provided faster generation than the C and R examples below, I don't think it's WRI that needs to do any fixing. Perhaps the Experimental`AngerManagement package needs a spin... – ciao – 2015-02-20T19:27:39.293

@AlbertRetey well said. – ciao – 2015-02-20T19:28:59.517

@rasher out of interest, what was your code for the parametric generation of 10^6 binomial variates? 200ms is close to the performance of the C++ code (150ms) – dr.blochwave – 2015-02-20T21:47:55.223

2@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 – 2015-02-20T22:00:11.203

@rasher look forward to seeing it, especially if it beats the C++! – dr.blochwave – 2015-02-20T22:01:42.453

@blochwave: Added... – ciao – 2015-02-20T22:12:47.340

If you are assuming a distribution for the input parameters on the parametric generation then you are not generating the random numbers that I need. – JEP – 2015-02-20T22:23:05.647

@JEP why not? The example provided by rasher below gives the same histogram as MMA or C++, and is very quick to generate and draw from (especially for very large draws of $>10^{7}$) – dr.blochwave – 2015-02-20T22:31:20.310

My input parameters come from a simulation. They don't have a known parametric distribution. – JEP – 2015-02-20T22:56:23.970

Ok, I understand now, thanks for clarifiying – dr.blochwave – 2015-02-20T23:07:07.727

1@jep Then you use the same technique where the parameter RVs are defined by your inputs. Pretty basic stuff. – ciao – 2015-02-20T23:45:02.397

Answers

7

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.

JEP

Posted 2015-02-18T06:25:31.647

Reputation: 1 238

6I 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. – 2015-02-19T00:09:58.683

2Thanks. 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 – 2015-02-19T19:06:31.973

1I don't really agree with rasher's answer, but I do think it's fair to point out that Mathematica's focus is usually on generality over run-time efficiency. The per-call overhead is often very high, especially for those functions implemented in terms of top-level code. This shouldn't necessarily be excused, but it is a fact one needs to know. One could read the R or Octave source code and reproduce the same implementation in Mathematica--clearly R is doing something clever if it can beat C++. Or, if this isn't very appealing, RLink and MATLink would allow you to call R or MATLAB directly. – Oleksandr R. – 2015-02-20T00:19:10.553

RLink huh? I've spent the last hour turning my R code into a script w/ command line arguments that go into rbinom. 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++. – JEP – 2015-02-20T00:51:47.597

3Using 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 – 2015-02-20T02:15:04.793

Thanks for the tip Oleksandr R. – JEP – 2015-02-20T02:23:17.843

If all I'm after is the raw generation of variates??? WTF? Why on earth would I be after the raw generation of variates? Where I come from we call "solutions" such as the one that yuou posted that received 14 upvotes "vaporware." I'd vote it down 14 times if they'd let me. – JEP – 2015-02-20T02:25:12.683

Oleksandr R. : I had never used Mathematica as a language for simulation and I guess you are saying that it isn't a good language for simulation. A year ago I wrote some simulation code in MMA because it was fairly easy to write. THe code was pretty slow, but fast enough that I could write the code and do the simulations within my deadlines. Now I am trying to determine whether I can continue using that code. The removal of this bottle neck might be enough to that I'll continue using MMA. – JEP – 2015-02-20T06:23:43.310

10@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 – 2015-02-20T08:12:52.140

4"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 – 2015-02-20T08:26:03.163

1I wouldn't say that Mathematica isn't good for writing simulations in general, but when you don't want to do any of the numerical implementation yourself, you are of course at the mercy of the standard library. It's to be expected that different packages have different strengths and weaknesses here. It also depends very strongly on the specific algorithm (although Mathematica is probably worse than most languages in terms of its extreme sensitivity to coding style--things like Julia may be better). I like Mathematica and I use it a lot, but it simply can't be the best tool for every job. – Oleksandr R. – 2015-02-20T08:49:27.350

1I agree blochwave. I did make a faulty comparison. Afterwards I tried to point that out in the comments but probably wasn't clear. As you point out it appears that C++ and R give comparable performance. WRI could easily fix this (for example by implementing a performance tuned version of RandomVariate, call it NRandomVariate or something) but based on the voting a substantial number of users prefer the poor performance. – JEP – 2015-02-20T15:23:05.353

37

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

dr.blochwave

Posted 2015-02-18T06:25:31.647

Reputation: 8 483

1Very nice answer, +1 – ciao – 2015-02-18T08:44:11.283

This is interesting. It doesn't work on my Mac. I put the c++ code in a file called poissonvariate.cpp which is in the working directory. It gives the following error message: Compile error: /Users/xxx/poissonvariate.cpp:4:10: fatal error:
'random' file not found >>
– JEP – 2015-02-18T13:45:05.193

Not having a Mac I can't check this for you, but perhaps you can check the examples in http://reference.wolfram.com/language/guide/LibraryLink.html and see if they work first?

– dr.blochwave – 2015-02-18T13:50:47.990

@JEP Perhaps also look at http://reference.wolfram.com/language/CCompilerDriver/ref/CreateLibrary.html

– dr.blochwave – 2015-02-18T13:52:18.937

1Notice that it allows you to define the C++ code inside Mathematica, rather than as a separate file. – dr.blochwave – 2015-02-18T13:53:11.430

I'm following your suggestions to look at the examples but in the mean time: what are the files "random" and "chrono" (from the 4th and 5th lines of your C++ ?). – JEP – 2015-02-18T14:04:57.083

<random> is to include std::random, and <chrono> is only to get the time for the random seed - it can be left out if you want to use the default random seed.

– dr.blochwave – 2015-02-18T14:28:43.567

Or you can use a given seed to compare the different methods. – dr.blochwave – 2015-02-18T17:01:57.570

1@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 – 2015-02-18T21:22:16.380

Thanks celtschk. I'll check on that. – JEP – 2015-02-18T22:34:51.253

OK. I ran the example here: http://stackoverflow.com/questions/14228856/how-to-compile-c-with-c11-support-in-mac-terminal

and added the line #include <random> and it compiled and ran so I think it works. Now to figure out how to tell Mathematica to use those options.

– JEP – 2015-02-18T22:42:08.913

@JEP Note that CreateLibrary has plenty of options including CompilerName and CompileOptions if you look under the Details tab in the documentation.

– dr.blochwave – 2015-02-19T08:16:53.387

Thanks. That will probably help. I'm getting some weird bugs right now. I've tried it on 3 version of Mac OS X: 10.7.4, 10.8.5, and 10.9.4 . I keep getting an error message: Compile error: clang: error: invalid deployment target for -stdlib=libc++ (requires OS X 10.7 or later) >>

INTERNAL SELF-TEST ERROR: Lex|c|3138 Click here to find out if this problem is known, and to help improve the Wolfram System by reporting it to Wolfram Research. "Click here" takes me here: http://support.wolfram.com/error-report/?Version=10.0.1&SystemID=MacOSX-x86&LicenseID=L2953-7873&location=Lex-c-3138

– JEP – 2015-02-19T19:12:57.483

Setting the flag: "CompilerName" -> "clang++" results in the error message reported in my previous comment: CreateLibrary::cmperr: Compile error: clang: error: invalid deployment target for -stdlib=libc++ (requires OS X 10.7 or later) >>

This is on OS X 10.8.5 but:

cat t.cpp #include <iostream> #include <random>

int main() { int* p = nullptr; std::cout << p << std::endl; std::cout << "Hello world"<<std::endl; }

clang++ -std=c++11 -stdlib=libc++ -O3 t.cpp

$ ./a.out 0x0 Hello world – JEP – 2015-02-19T21:30:37.110

blochwave and celttschk: I appreciate your efforts on this. You might be interested in my RLink solution below. With 4 lines of code you are set up. You can then define whatever generators you want with a single line of code for each generator. You can call them with lists as arguments and as far as I know it just works. Less code to maintain and I believe it is substantially faster than your C++ version. And it works on my Mac. :-) – JEP – 2015-02-20T03:13:29.950

I spoke to soon. On the Poisson deviate with calling R from Mathematica I only get a speed up of 85. The speed up on the binomial is about a factor of 260. – JEP – 2015-02-20T05:54:23.963

@JEP you'll see that the performance on the binomial distribution in C++ is the same as R... – dr.blochwave – 2015-02-20T08:39:24.697

Yes. When I said I "spoke to soon" I was trying to point out that your C++ Poisson solution is nearly twice as fast as my R one. They're both far faster than native MMA. For me the wins with R are (1)from the moment Oleksandr told me to try RLink it took 15 or 20 minutes and I had working code, while I spent the better part of two days failing to get your working C++ to work on my Macs & (2)the R sol'n takes substantially less code. Your C++ wins as you can performance tune the C++, if need be. It might not help because the performance might be bound by the overhead calling R or C++ from MMA. – JEP – 2015-02-20T15:06:41.117

I've checked, the overhead calling the C++ from MMA is negligible. Regarding the code not working, I think it must be Mac specific I'm afraid, this works for me on Windows and Linux. – dr.blochwave – 2015-02-20T15:23:12.283

2

@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 – 2015-02-21T08:56:26.347

1@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 – 2016-02-03T18:51:10.560

@RolfMertig thanks - I'll check this out when I'm back at a PC with MMA tomorrow! – dr.blochwave – 2016-02-03T19:10:06.493

@blochwave Sorry, the comparison should be AbsoluteTiming[ l2 = Flatten@ MapThread[BinomialVariate[{#1}, {#2}] &, {list1, list2}];] Then I get about 23 fold time difference (on windows). – Rolf Mertig – 2016-02-03T23:43:30.320

1I get a little over 21s on 10.4.1, and 18.5 on 11.0.1. So, it's "improved." :) – rcollyer – 2016-10-06T12:17:07.470

34

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.

ciao

Posted 2015-02-18T06:25:31.647

Reputation: 23 752

4An example would have been helpful. – JEP – 2015-02-18T13:07:18.463

1The 300x speedup of weibull is impressive. On this problem Mathematica is only 123 times slower than Octave. cat wbl.m t=cputime;num=10^5;n=50rand(num,1);p=100rand(num,1);out=wblrnd(n,p);cputime-t $ octave wbl.m

ans = 0.018561 – JEP – 2015-02-19T23:23:58.160

7@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 – 2015-02-19T23:32:35.863

1It would be stupid not to complain about a slow buggy $1500 software package that claims about itself: "built to provide industrial-strength capabilities—with robust, efficient algorithms across all areas" & "draws on its algorithmic power—as well as the careful design of the Wolfram Language—to create a system that's uniquely easy to use" & this "With its intuitive English-like function names and coherent design, the Wolfram Language is uniquely easy to read, write, and learn" and which in addition has been unresponsive to bug reports. It is stupid to defend it. – JEP – 2015-02-20T19:23:50.283

7@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 – 2015-02-20T19:44:01.853

@rasher - Cheers for adding the code, it's very quick - 0.03 seconds on my machine, compared to the 0.15 seconds for C++, although a fairer comparison including the generation of the mixture distribution bumps that up to 0.13 seconds. – dr.blochwave – 2015-02-20T22:21:28.427

@blochwave: Since the distribution is defined for the generation, that cost is amortized- there's no extra cost generating 10 or 10 Billion variates. Only if the distribution's parameter RVs distributions themselves changed would that come into play, and then there's other ways to skin that cat... so I think it more than "fair" as is, not to mention shows techniques that cannot even be done with purely numeric s/w. Thanks for timing on your box! – ciao – 2015-02-20T22:59:30.590

1Mr Wizard, I am not confused. If you are going to start chastising me for my tone you ought to start at the beginning to this thread and read the links. I was told I was brain dead, then lazy, then stupid. All for expecting MMA to be able to generate random variates with performance similar to that of other languages. – JEP – 2015-02-20T23:03:32.923

@JEP I fail to see where "brain dead" originated. However calling people stupid, even by implication, is not appropriate. rasher: it should be possible to argue a position without stooping to ad hominem attacks, which you are flirting with above; please do not cross that line. – Mr.Wizard – 2015-02-20T23:28:40.873

2

@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 – 2015-02-20T23:32:53.557

6To @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 – 2015-02-20T23:37:48.397

1I'm glad to see that you've added code. THis is not a solution to the general problem in which the input parameters are not parametrically distributed. – JEP – 2015-02-20T23:40:22.163

Mr Wizard, I didn't really want you to go through and read all these links! That is a waste of time. The "brain dead" was rasher's response the first time I raised this issue. Yes. I thought people came here to to help and be helped. I don't know how to make it clearer that I was looking for a solution than posting the solution as I did. blochwave posted a fine solution which I couldn't get to work on my Macs. I imagine my solution will work on any OSX Mac. It received two up votes whereas rasher's solution which does not solve the general problem received 17. Go figure. – JEP – 2015-02-20T23:47:43.333

2@Mr.Wizard: I get your point, but a "That car is crap and needs fixin' 'cause I crashed" post when it turns out the driver didn't know how to drive is... well, I'll leave it to you to judge. Now the question appears to have changed - nowhere in the OP is the fact that the distribution of the parameters themselves is some undefinable construct. My vote to close stands - this feels more like a marketing rant from a competitor. I shall refrain from commenting further on this post, the hackles are at attention. – ciao – 2015-02-20T23:50:40.803

2

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]

enter image description here

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

Henrik Schumacher

Posted 2015-02-18T06:25:31.647

Reputation: 85 430