Fast generation of random variates as test case
The earlier question (195435) about performance tuning with regard to fast generation of random variates for a gamma distribution has stirred my curiosity and below I am using slightly simpler code from the 3rd ed. of Numerical Recipes (Normaldev) to test different versions of a numerical routine to generate random variates for a normal distribution.
Different Implementations
The following implementations use the same core body, but differ with regard to the way a function is declared in Mathematica:
normaldev
is an uncompiled implementationnormaldevC
usesCompile
withCompilationTarget → "C"
normaldevCLP
usesCompile
as above, but adds listablility and parallelizationnormaldevLLVM
finally uses the new compiler andFunctionCompile
Code (`body`
has been escaped as\`body\`
):
{ normaldev, normaldevC, normaldevCLP, normaldevLLVM } = With[
{
funcBody = "
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.25482 * x);
If[ q > 0.27597 && ( q > 0.27846 || v*v > -4. * Log[u] * u * u ), cond = True, cond = False ]
];
mu + sig * v/u"
},
Module[
{ normaldev, normaldevC, normaldevCLP, normaldevLLVM }
,
normaldev = StringTemplate["
Function[ {mu, sig },
Module[
{
u = 1.0,
v = 1.0,
x = 1.0,
y = 1.0,
q = 1.0,
cond = True
}
,
\`body\`
]
]"
];
normaldevC = StringTemplate["
Compile[ { { mu, _Real }, { sig, _Real } },
Module[
{
u = 1.0,
v = 1.0,
x = 1.0,
y = 1.0,
q = 1.0,
cond = True
}
,
\`body\`
],
CompilationTarget -> \"C\"
]"
];
normaldevCLP = StringTemplate["
Compile[ { { mu, _Real }, { sig, _Real } },
Module[
{
u = 1.0,
v = 1.0,
x = 1.0,
y = 1.0,
q = 1.0,
cond = True
}
,
\`body\`
],
CompilationTarget -> \"C\",
RuntimeAttributes -> {Listable},
Parallelization -> True
]"
];
normaldevLLVM = StringTemplate["
FunctionCompile[
Function[
{
Typed[ mu, \"Real64\" ],
Typed[ sig, \"Real64\" ]
}
,
Module[
{
Typed[ u, \"Real64\" ] = 1.0,
Typed[ v, \"Real64\" ] = 1.0,
Typed[ x, \"Real64\" ] = 1.0,
Typed[ y, \"Real64\" ] = 1.0,
Typed[ q, \"Real64\" ] = 1.0,
Typed[ cond, \"Boolean\" ] = True
}
,
\`body\`
]
],
CompilerOptions -> {
\"AbortHandling\" -> False,
\"LLVMOptimization\" -> \"ClangOptimization\"[3]
}
]"
];
Map[ ToExpression @ TemplateApply[ #, <| "body" -> funcBody |>]&, { normaldev, normaldevC, normaldevCLP, normaldevLLVM } ]
]
];
Performance evaluation
Using RepeatedTiming
we can compare run times for these implementations, which I have sorted below from slowest to fastest.
RepeatedTiming @ Through[ {Mean, StandardDeviation}@ Table[ normaldev[0.,1.], {1000000}] ]
(* {0.535271,{-0.00048065,1.00115}} *)
RepeatedTiming @ Through[ {Mean, StandardDeviation}@ Table[ normaldevLLVM[0., 1.], {1000000}] ]
(* {0.199849,{0.00011759,1.00026}} *)
RepeatedTiming @ Through[ {Mean, StandardDeviation}@ Table[ normaldevC[0.,1.], {1000000}] ]
(* {0.18677,{-0.00244043,0.998699}} *)
With[ { args = Transpose @ ConstantArray[{0.0,1.0},1000000]},
RepeatedTiming @ Through[ {Mean, StandardDeviation}@ ( normaldevCLP @@ args )]
]
(* {0.119437,{0.000838859,0.999507}} *)
Questions
- What would be the recommended implementation for the
LLVM
version if speed is the main goal, e.g., setting options, usingNative`Random
instead ofRandomReal
? - How can the benefits of listability and parallelization be obtained for
Compiled CodeFunction
?