Which Distributions can be Compiled using RandomVariate

Question

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)?

Answer 1

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.

Answer 2

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.

Answer 3

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

It looks fine compared to the built-in version:

Answer 4

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.