35
$\begingroup$

Recently, Oleksandr kindly showed a list of Mathematica commands that can be compiled. RandomVariate was part of that list. However, whether this can be compiled depends upon the distribution that is being sampled.

Needs["CompiledFunctionTools`"]

cf1 = Compile[{{m, _Real}, {s, _Real}},
  Module[{v1, v2, v3, v4, v5, v6},
   v1 = RandomVariate[NormalDistribution[m, s]];
   v2 = RandomVariate[UniformDistribution[{m, s}]];
   v3 = RandomVariate[GammaDistribution[m, s]];
   v4 = RandomVariate[PoissonDistribution[m]];
   v5 = RandomVariate[ChiSquareDistribution[m]];
   v6 = RandomVariate[ExponentialDistribution[m]];
   {v1, v2, v3, v4, v5, v6}
   ]
  ]

Using CompilePrint shows that RandomVariate can be compiled for the Normal Distribution or the Uniform Distribution and not with some others.

CompilePrint[cf1]

2 arguments
        4 Integer registers
        8 Real registers
        1 Tensor register
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        R0 = A1
        R1 = A2
        Result = T(R1)0

1   R2 = RandomNormal[ R0, R1]]
2   R3 = RandomReal[ R0, R1]]
3   I0 = MainEvaluate[ Function[{m, s}, 
             RandomVariate[GammaDistribution[m, s]]][ R0, R1]]
4   I1 = MainEvaluate[ Function[{m, s}, 
             RandomVariate[PoissonDistribution[m]]][ R0, R1]]
5   I2 = MainEvaluate[ Function[{m, s}, 
             RandomVariate[ChiSquareDistribution[m]]][ R0, R1]]
6   I3 = MainEvaluate[ Function[{m, s}, 
             RandomVariate[ExponentialDistribution[m]]][ R0, R1]]
7   R4 = I0
8   R5 = I1
9   R6 = I2
10  R7 = I3
11  T(R1)0 ={ R2, R3, R4, R5, R6, R7 }
12  Return

Does anyone have a list of all the distributions that can be compiled (including, PDF, CDF and RandomVariate functionality)?

$\endgroup$
4
  • $\begingroup$ Well, you can check each of them, or semi-automate checking. When one is not compiled, and you try to use it, you'll easily find out anyway: see my answer to that post asking about compilable functions. $\endgroup$
    – Szabolcs
    Feb 1, 2012 at 19:33
  • $\begingroup$ @szabolcs There are too many combinations to try out, unfortunately. $\endgroup$
    – asim
    Feb 2, 2012 at 1:43
  • 3
    $\begingroup$ Beta distribution can be compiled too.. $\endgroup$
    – asim
    Feb 10, 2012 at 14:09
  • 1
    $\begingroup$ Since the results of RandomVariate are ALWAYS going to be Reals (for a continuous distribution), and Integers (for most inbuilt discrete distributions), I can't see why there would be any advantage to compiling them, that could not have automatically been built into the function. Does anyone have some timing tests where manual compilation yields substantive benefits to something like: RandomVariate[ dist, {10^6}] ? $\endgroup$
    – wolfies
    Aug 5, 2014 at 15:37

4 Answers 4

26
$\begingroup$

To my knowledge UniformDistribution and NormalDistribution are the only distributions that are directly compilable for RandomVariate.

Consider that sampling from a UniformDistribution is what RandomReal was originally designed to do. This code is likely written deep down in C and so compiles without any special effort. In order to hook up RandomVariate for uniforms Compile just needs to recognize that this is really just a call to RandomReal.

Now, sampling from a NormalDistribution is so common that it was considered worth the time investment to make it compilable. Notice that the call to RandomVariate actually produces a call to RandomNormal which was almost certainly written for this purpose.

As for other distributions, special code would need to be written for each one in a similar fashion to RandomNormal for them to be "supported" by Compile. Since there are well over 100 of these, it would be a huge undertaking. An argument could be made for doing this for a few distributions but who is to decide which ones are most important?

There is a sunny side. Most distributions have their own dedicated and highly optimized methods for random number generation. Often Compile is used under the hood when machine precision numbers are requested.

Because of this, even if they were directly compilable you probably wouldn't see much of a speed boost since the code is already optimized.

Fortunately Compile can happily handle arrays of numbers. I typically just rely on the optimized code used by RandomVariate to generate the numbers and subsequently pass them in as an argument to the compiled function.

Incidentally, everything I just said about RandomVariate is also true of distribution functions like PDF, CDF, etc. Obviously these are just pure functions (in the univariate case) and unless they are built with some exotic components they should compile assuming you evaluate them before putting them into your compiled function.

$\endgroup$
4
  • $\begingroup$ Also of note is that Erf[] is in the list of compilable functions in Olek's answer. So yes, I think the normal distribution is pretty well-covered here. $\endgroup$ Feb 2, 2012 at 6:29
  • 1
    $\begingroup$ I just found out that RandomVariate[BetaDistribution[a,b]] can also be compiled. By the way, it is not always that simple to generate the random numbers before passing them to a compiled function, especially, if the number of draws are endogenously determined within the function, as part of some while loop. $\endgroup$
    – asim
    Feb 10, 2012 at 14:05
  • $\begingroup$ @asim this is definitely true. I was unaware the BetaDistribution was compilable, thanks! $\endgroup$
    – Andy Ross
    Feb 10, 2012 at 16:09
  • $\begingroup$ @asim I have had situations like this. Typically I generate far more numbers than I think I will need and pass them in. It is faster to generate a bunch in one go than to repeatedly call RandomVariate. The downside is having to leave Compile and generate more if I didn't generate enough numbers in the first go. $\endgroup$
    – Andy Ross
    Feb 10, 2012 at 23:22
20
$\begingroup$

Another option to the answer posted by Andy Ross cropped up in a recent question of mine about corrupting an image with Poisson noise. In my own answer, I made use of LibraryLink to utilise the distributions built into C++.

This was especially useful in my case because Poisson noise in an image otherwise relies on a call to RandomVariate for each pixel (although the built-in Mathematica function for doing this via ImageEffect works differently, but is still slower than my answer).

Repeated calls to RandomVariate rapidly become time-consuming. Here I can get a ~200x speed-up!

(* Generate random list *)
listofnums1d = RandomReal[{0, 100}, {10^6}];

(* Calling RandomVariate each time takes 146.78 seconds *)
RandomVariate[PoissonDistribution[#]] & /@ listofnums1d; // AbsoluteTiming

(* Loading my C++ function via LibraryLink (code below) takes 0.77 seconds *) 
mypoissondistribution[listofnums1d]; // AbsoluteTiming

C++ has a whole host of distributions built into <random>, including Gamma, Exponential, ChiSquare, Cauchy and Weibull.

If repeated calls to RandomVariate are a problem, this could be a way to go.


C++ and LibraryLink code

#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 poissondistribution(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

    int err; // error code

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

    mint const* dims; // dimensions of the tensor

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

    mint i; // bean counter

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

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

    for(i=0;i<dims[0];i++) {
        if(data1[i]==0.)
        {
            data2[i] = 0.; 
            // If value is 0, return 0. 
            // You can avoid this by not passing a value of 0, of course.
        }
        else
        {
            std::poisson_distribution<int> poissondistribution(data1[i]); // Poisson distribution
            data2[i] = poissondistribution(generator); // If value is >0, sample from the Poisson distribution
        }
    }

    MArgument_setMTensor(Res, m2);
    return LIBRARY_NO_ERROR;
}

Put the C++ code into a file (e.g. in the notebook directory) and you can compile it like this.

Needs["CCompilerDriver`"]

distributionlib = 
  CreateLibrary[{FileNameJoin[{NotebookDirectory[], "filename"}]},
  "libraryname", "Debug" -> False, 
  "CompileOptions" -> "-O2 -std=c++11"];

Once the library has been compiled and you stop making changes to the C++ code, you can save time by taking it from the default location Mathematica places it ($LibraryPath, found using FindLibrary[]), and move it into your notebook directory and call it. Remember that the file extension is .dll on Windows, .so on Linux, and .dylib on Mac.

This works because LibraryFunctionLoad[] uses FindLibrary internally.

mypoissondistribution = 
  LibraryFunctionLoad[
   FileNameJoin[{NotebookDirectory[], "filename"}], 
   "functionname", {{Real, 1}}, {Real, 1}];

Note that this code takes and returns a 1D list of values. The code in my linked question about images takes a 2D list.

I believe that you can then include mypoissondistribution in a compiled function without a call to MainEvaluate. Perhaps someone can clarify? I don't know if you need to set "InlineExternalDefinitions" and/or "InlineCompiledFunctions" to True for this to work.

$\endgroup$
1
  • 1
    $\begingroup$ This is very neat!! +1 $\endgroup$
    – brama
    Aug 6, 2014 at 10:28
9
$\begingroup$

Here is a quick (and dirty?) implementation of the Laplace distribution. Relationship to the uniform distribution as given on Wikipedia:

randomLaplace = With[{$MachineEpsilon = $MachineEpsilon},
  Compile[{{µ, _Real, 0}, {b, _Real, 0}, {n, _Integer, 0}},
   With[{u = RandomReal[{-1/2, 1/2} (1 - $MachineEpsilon), n]},
    µ - b Sign[u] Log[1 - 2 Abs[u]]
   ]
  ]
 ];

Histogram of compiled Laplace distribution results

It looks fine compared to the built-in version:

Histogram of built-in Laplace distribution results

$\endgroup$
0
3
$\begingroup$

Another distribution which is missing is MultinormalDistribution[v,S]. However, variates from this distribution are just affine transformations of vectors of iid normal variates. Indeed, if y=RandomVariate[NormalDistribution[],Length@s], then x=v+A.y has the distribution of interest, where A is any decomposition of S such that S==A.Transpose[A].

Since RandomVariate[NormalDistribution[],{d1,d2,..}] can be compiled, it suffices to implement a Cholesky decomposition inside Compile (the builtin one cannot be compiled). For example, we can steal the procedural algorithm from this post. We slightly modify it to find the transpose of the Cholesky decomposition, so that the matrix multiplication is more natural with our array shape, and we also only allow real matrices.

multinormalVariate = Compile[{{v, _Real, 1}, {S, _Real, 2}, {n, _Integer}},
  Module[{d, A, iidnorm, x},
    dimens = Length@S;
    A = Table[0.0, {j, d}, {i, d}];
    Do[
      Do[
        If[i == j,
          A[[i, i]] = Sqrt[S[[i, i]] - Sum[A[[k, i]]^2, {k, 1, i - 1}]];,
          A[[j, i]] = (S[[i, j]] - Sum[A[[k, i]]*A[[k, j]], {k, 1, j - 1}])/A[[j, j]];
        ];,
        {i, j, d}];,
      {j, d}
    ];
    (*A is triangular, but Dot is faster than taking advantage of this in a loop*)
    (v + #) & /@ (RandomVariate[NormalDistribution[], {n, dimens}].A)
  ],
  CompilationTarget->"C"
];

In my tests, this function scales about the same as the builtin function with respect to n, where most of the work is in matrix multiplication and choosing samples. However, the builtin one scales better for large d. A more efficient decomposer would be needed. For d that aren't too large (<~500), this still may be useful in a large compiled function where you want to limit calls to MainEvaluate.

Module[{d = 10, n = 10000, S, A, v, data1, data2, t},
  A = RandomReal[{0, 1}, {d, d}];
  S = A.Transpose[A];
  v = RandomReal[{-2, 2}, d];
  {t, data1} = AbsoluteTiming@multinormalVariate[v, S, n]; Print@t;
  {t, data2} = AbsoluteTiming@RandomVariate[MultinormalDistribution[v, S], n]; Print@t;
  Print[Norm /@ {v - Mean@data1, v - Mean@data2}];
  Print[Norm /@ {Covariance[data1] - S, Covariance[data2] - S}]
]

Edit: this decomposition will fail for singular or poorly conditioned S.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.