Why doesn't R use Inverse Transform Sampling to sample from the Exponential Distribution?

Question

I was reading this question about the algorithm that R uses to sample from the Exponential($\lambda$) distribution. It looks like R uses the Ahrens-Dieter algorithm to sample from the exponential distribution. The reference for it is a paper from 1972 that makes a surprising claim that this algorithm is a little faster (in terms of microseconds, not in terms of Big O computational cycles) than just delivering $\frac{-\ln(u)}{\lambda}$ (where U~Uniform(0,1); as the Inverse Transform would have us do).

I find this very surprising. First of all, I am not sure the results from 1972 carry much weight in 2022. Secondly, why the need for an algorithm in the first place as opposed to just delivering $\frac{-\ln(u)}{\lambda}$ (as the Inverse Transform would have us do)? How slow can calculating $\frac{-\ln(u)}{\lambda}$ be? Maybe it was slow in 1972. Or maybe there is another reason?

Answer 1

I have written rather a lot of code to test this on two machines, the first of which is a Windows machine that has an Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz processor. I am using a recent MinGW compiler via the CodeBlocks IDE, with optimizations set at O3 and expensive_optimizations enabled, but no chip-specific optimizations.

I have run the same code on an Intel(R) Xeon CPU E5-2650L v4 @ 1.70GHz processor, Ubuntu OS, and GCC 7.5 compiler, with the same compiler options.

I am using the pcg32 random number generator to generate the required uniform variates inside the code helpfully provided by @Alex. pcg32 is extremely fast - my implementation, copied from who knows where, takes only 18% longer than the C library rand() function while returning an unsigned integer between 0 and 4294967296$ = 2^{32}$, whereas rand() returns a signed integer between 0 and 32767 (much poorer granularity), and has excellent properties. See https://www.pcg-random.org/ for more.

Fast Windows machine: Generating 10 million variates using the pcg32 version of the code that R uses took 568,095 microseconds, including the overhead induced by the for loop.

Generating 10 million variates using the inverse probability transform took 527,346 microseconds, including the overhead induced by the for loop. This is roughly 93% of the time that the A-D algorithm uses. Roughly 40% of the time for the inverse probability transform algorithm appears to be loop overhead and the uniform RNG.

Slower Linux machine: the A-D algorithm took 666,358 microseconds, the inverse probability transform took only 555,192 microseconds. This is roughly 83% of the time that the A-D algorithm uses.

These results certainly validate the OP's suspicion that things may have changed since the 1970s. Regardless of the algorithm choice, being able to generate 15-20 million exponential random numbers in one second on one thread is certainly a nice capability to have!

One interesting finding is that the runtimes don't scale with the CPU speed; they are a little greater on the Linux box than on the Windows box, but not nearly the 2x+ difference in speed. Of course, the compilers are different, and GHz is by no means the sole influence on runtime.

I'd be happy to post the code, but there's about 130 lines of it, including some comments and blank lines. Thoughts?

Answer 2

Just for the review sake, the source code for this function in R (written in the C language):

#include "nmath.h"

double exp_rand(void)
{
    /* q[k-1] = sum(log(2)^k / k!)  k=1,..,n,
       The highest n (here 16) is determined by q[n-1] = 1.0
       within standard precision */
    const static double q[] = {
        0.6931471805599453,
        0.9333736875190459,
        0.9888777961838675,
        0.9984959252914960,
        0.9998292811061389,
        0.9999833164100727,
        0.9999985691438767,
        0.9999998906925558,
        0.9999999924734159,
        0.9999999995283275,
        0.9999999999728814,
        0.9999999999985598,
        0.9999999999999289,
        0.9999999999999968,
        0.9999999999999999,
        1.0000000000000000
    };

    double a = 0.0;
    double u = unif_rand();    /* precaution if u = 0 is ever returned */
    while(u <= 0.0 || u >= 1.0) u = unif_rand();
    for (;;) {
        u += u;
        if (u > 1.0) break;
        a += q[0];
    }
    u -= 1.0;

    if (u <= q[0]) return a + u;

    int i = 0;
    double ustar = unif_rand(), umin = ustar;
    do {
        ustar = unif_rand();
        if (umin > ustar) umin = ustar;
        i++;
    } while (u > q[i]);
    return a + umin * q[0];
}
```