New Compiler Technology in V12

Question

I am wondering whether the speed of compilation for the new compiler can be improved.

Here are the results from my test.

Consider the following function for generating a random draw from a gamma distribution

ablRanGamma[alpha_, beta_] := 
    Module[
      {d, c, x = 0.0, v = 0.0, u, cond},

        d = alpha - 1.0/3.0;
        c = 1.0/Sqrt[9.0*d];
        cond = 1;

        While[cond == 1,
            v = -1.0;
            While[v <= 0.0,
                    x = RandomVariate[NormalDistribution[0, 1]];
                    v = 1.0 + c*x;
                ];
            v = v^3;
            u = RandomReal[];
            If[u > 1. - 0.0331*x^4 && 
         Log[u] > 0.5*x^2 + d*(1.0 - v + Log[v]), cond = 1, cond = 0];
        ];
        beta*d*v
        ]

The old compiler technology can be used to compile the function to C, as follows

    ablRanGammaC = Compile[{{alpha, _Real}, {beta, _Real}},
   Module[
        {d, c, x = 0.0, v = 0.0, u, cond},
        d = alpha - 1.0/3.0; c = 1.0/Sqrt[9.0*d]; cond = 1;
        While[cond == 1,
            v = -1.0;
            While[v <= 0.0,
                x = RandomVariate[NormalDistribution[0, 1]];
                v = 1.0 + c*x;
            ];
            v = v^3;
            u = RandomReal[];
            If[
      u > 1. - 0.0331*x^4 && Log[u] > 0.5*x^2 + d*(1.0 - v + Log[v]), 
      cond = 1, cond = 0];
        ];
        beta*d*v
        ], 
    CompilationTarget -> "C"
   ];

This results in a significant improvement is speed, as expected.

In[10]:= AbsoluteTiming[Mean@Table[ablRanGamma[2.1, 3.2], {500000}]]

Out[10]= {11.7353, 6.71994}

In[7]:= AbsoluteTiming[Mean@Table[ablRanGammaC[2.1, 3.2], {500000}]]

Out[7]= {0.142951, 6.72236}

Now moving to the new compiler. First, RandomVariate[NormalDistribution[0,1]] cannot be compiled directly, so we need to use a KernelFunction. We create two functions to get random normal draws, one compiled and the other not.

f1[m_, s_] := RandomVariate[NormalDistribution[m, s]]

f1C = Compile[{{m, _Real}, {s, _Real}}, 
   RandomVariate[NormalDistribution[m, s]],
   CompilationTarget -> "C"
   ];

Then we can use these as KernelFunctions using the new compiler. Consider the function which uses f1.

ablRanGammaN1 = FunctionCompile[
   Function[{Typed[alpha, "Real64"], Typed[beta, "Real64"]},
    Module[
     {d, c, x = 0.0, v = 0.0, u, cond, 
      f1 = Typed[
        KernelFunction[f1], {"Real64", "Real64"} -> "Real64"]},
        d = alpha - 1.0/3.0; c = 1.0/Sqrt[9.0*d]; cond = 1;
        While[cond == 1,
            v = -1.0;
            While[v <= 0.0,
                x = f1[0.0, 1.0];
                v = 1.0 + c*x;
            ];
            v = v^3;
            u = RandomReal[];
            If[
       u > 1. - 0.0331*x^4 && Log[u] > 0.5*x^2 + d*(1.0 - v + Log[v]),
        cond = 1, cond = 0];
        ];
        beta*d*v
        ]
    ]
   ];

Its performance is

In[8]:= AbsoluteTiming[Mean@Table[ablRanGammaN1[2.1, 3.2], {500000}]]

Out[8]= {2.46856, 6.72463}

This can be improved by using f1C instead of f1 as a kernel function.

ablRanGammaN2 = FunctionCompile[
   Function[{Typed[alpha, "Real64"], Typed[beta, "Real64"]},
    Module[
     {d, c, x = 0.0, v = 0.0, u, cond, 
      f1 = Typed[
        KernelFunction[f1C], {"Real64", "Real64"} -> "Real64"]},
        d = alpha - 1.0/3.0; c = 1.0/Sqrt[9.0*d]; cond = 1;
        While[cond == 1,
            v = -1.0;
            While[v <= 0.0,
                x = f1[0.0, 1.0];
                v = 1.0 + c*x;
            ];
            v = v^3;
            u = RandomReal[];
            If[
       u > 1. - 0.0331*x^4 && Log[u] > 0.5*x^2 + d*(1.0 - v + Log[v]),
        cond = 1, cond = 0];
        ];
        beta*d*v
        ]
    ]
   ];

We now have

In[9]:= AbsoluteTiming[Mean@Table[ablRanGammaN2[2.1, 3.2], {500000}]]

Out[9]= {0.81675, 6.72306}

This is clearly better than the uncompiled function, but still is significantly slower than when we compile to C using the old technology. Can anybody explain why this is the case. In my experiments, this does not seem to be due to the overhead of using a kernel function. The C compiled code appears to be faster than the machine compiled code for the new technology.

Answer 1

I wanted to look into this for a long time, and it took only 2 months to come back to this question. To answer your main question: Yes, it appears that the new compiler only supports uniformly distributed random numbers. There are three native functions: Native`Random, Native`RandomGeneratorBitFunction, and Native`RandomRuntime. If you call RandomReal it is compiled to Native`Random. I'll come back to the issue, why it's so hard to find information on how to use them.

To give you a solution to your problem first: Your function basically uses random numbers from a normal- and a uniform-distribution to compute a random number which is Gamma-distributed. I believe this is Marsaglia's simple transformation-rejection method. What you haven't considered is that you can also compute a normally-distributed random number from uniformly distributed random numbers. An easy way is the Box-Muller method which only uses elementary functions and two random reals.

You can replace your f1 entirely by this simple and calculate and assign a random number which comes from a (0,1)-normal-distribution. The code would then look like this:

f = Function[{Typed[alpha, "Real64"], Typed[beta, "Real64"]},
   Module[
    {d, c, x = 0.0, v = 0.0, u, cond},
    d = alpha - 1.0/3.0;
    c = 1.0/Sqrt[9.0*d];
    cond = 1;
    While[cond == 1,
     v = -1.0;
     While[v <= 0.0,
      x = Sqrt[-2.0*Log[Native`Random[0.0, 1.0]]]*Cos[2.0*Pi*Native`Random[0.0, 1.0]];
      v = 1.0 + c*x
      ];
     v = v^3;
     u = Native`Random[0., 1.];
     If[
      u > 1. - 0.0331*x^4 && Log[u] > 0.5*x^2 + d*(1.0 - v + Log[v]),
      cond = 1,
      cond = 0
      ];
     ];
    beta*d*v]];

ablRanGammaN3 = FunctionCompile[f];

The timings on my machine are:

• uncompiled: 10.6206
• compiled C: 0.126193
• ablRanGammaN1: 2.27945
• ablRanGammaN2: 0.80765
• ablRanGammaN3: 0.369398

You see, we stole another 50% of the runtime but we still have 3x the runtime of the old compiler. And the histograms give hope that my implementation is not incorrect

Gain more insight

I have followed the rise of the new compiler pretty much from the start but I admit that I haven't used it very often. The big disadvantage is that the most reliable documentation resources are still the videos on Youtube (video 1, video 2, video 3, video 4).

If you want to look a bit deeper into it, my first advice is that you use the functions under

<<Compile`

and in particular CompileToCodeFunction which seems to do the same as FunctionCompile but has many options and you can use to as an entry point to inspect the underlying code using PrintDefinitions. Other than that, you have CompileToAST, CompileToIR, etc which really help because you can call ["toString"] on their output.

Another thing to look into is Names["Native`*"] which seems to be all the functions that are used internally and which you can also use directly in your code. The downside is that I'm currently not aware of how we can find out how to use them. If you watch the videos I linked, you'll see that they are heavily used for performance tweaking and to represent things the compiler cannot do at the moment. Pay close attention especially to the sort example of Tom-Wickham Jones at the end of his video.

Finally, the two major advantages of the new technology for me personally are:

• we can link functions written in C/C++, Rust, ... in an instance. There are a lot of examples in one of the videos
• we can create code for different platforms like web-assembly

It's just that the new compiler isn't there yet when it comes to documentation. At the moment it's more a playground for the adventurous heroes among us.

Answer 2

It seems that there has been a little bit improvement as of v13.1 but I still find the documentation rather sparse.

Adding to the examples given above I have implemented the Normal Deviates by Ratio-of-Uniforms method given in Numerical Recipes, 3rd ed. (deviates.h) using FunctionDeclaration so that normaldev can then be used within ablRanGammaN4 without further declaration.

$CompilerEnvironment = CreatCompilerEnvironment[];
CompilerEnvironmentAppendTo[
 FunctionDeclaration[
  normaldev,
  Typed[{"Real64", "Real64"} -> "Real64" ] @ Function[
    {
     mu, sig
     }
    ,
    Module[
     {
      (* typed declarations need to be initialized *)
      Typed[u, "Real64"] = 1.,
      Typed[v, "Real64"] = 1.,
      Typed[x, "Real64"] = 1.,
      Typed[y, "Real64"] = 1.,
      Typed[q, "Real64"] = 1.,
      Typed[cond, "Boolean"] = True,
     }
     ,
     While[
      cond,
      (
        u = RandomReal[];
        v = 1.7156 * (RandomReal[] - 0.5);
        x = u - 0.449871;
        y = Abs[v] + 0.386595;
        q = (x*x) + y *(0.19600*y - 0.25472 * x);
        If[ q > 0.27597 && ( q > 0.27846 || (v*v) > -4. Log[u] * (u*u)), 
            (* then *) cond = True, 
            (* else *) cond = False 
        ];
      )
     ];
     mu + sig * v/u 
    ]
  ]
 ]
]

The implementation for ablRanGamma is only slightly modified. Note, that I tend to get error messages like Compile::err: TypeError. The value of '...' has not been assigned on all branches (types Real64 and Uninitialized are not equal) when I am not used Typed in the scoping construct. If Typed is used, then it seems an initial value has to be assigned as well.

ablRanGammaN4 = FunctionCompile[
  Function[
   {
    Typed[alpha, "Real64"], 
    Typed[beta, "Real64"]
    }, 
   Module[
    {
     Typed[d, "Real64"] = alpha - 1./3., 
     Typed[c, "Real64"] = 1., 
     Typed[x, "Real64"] = 0., 
     Typed[v, "Real64"] = 0., 
     Typed[u, "Real64"] = 1.,
     Typed[cond, "Boolean"] = True
     },
    c = c /Sqrt[9.* d];
    While[ 
     cond, 
     (
      v = -1.; 
      While[ v <= 0.
       , 
       (
         x = normaldev[0., 1.];
         v = 1. + c * x;
       )
      ];
      v = v^3; 
      u = RandomReal[];
      If[ u > 1. - 0.0331*x^4 && Log[u] > 0.5*x^2 + d*(1. - v + Log[v]), 
       cond = True, 
       cond = False 
      ];
     )
    ]; 
    beta * d * v
   ]
  ],
  CompilerOptions -> {
    "InvocationMode" -> "Standalone",
    "AbortHandling" -> False,
    "LLVMOptimization" -> "ClangOptimization"[3]
  }
]

The timings I get with AbsoluteTiming on my machine are:

• uncompiled: 12.674 s
• compiled C: 0.086 s
• ablRanGammaN4: 0.114 s

So FunctionCompile is taking less than 1.5x the old compiler times. Unfortunately, CompilerOptions are not really documented, so I have been experimenting with what is stated here or here.