I am obtaining slow gains in compiling a function, and would like to ask if there is something I am obviously missing.

The function I intend to speed up has the form:

        PDF[NormalDistribution[.2, 1], r]  / CDF[NormalDistribution[.2, 1],r],
        {r, 1, 100000}]

where in large evaluations the number of elements in the table, and the parameters of the Normal distribution, may change.

This way of evaluating will return a result with a timing of around 8.1 seconds.

I am compiling this function as:

t1 = Compile[{{r, _Real, 0}},
      PDF[NormalDistribution[.2, 1], r]  / CDF[NormalDistribution[.2, 1], r], 
      RuntimeAttributes -> {Listable}];

t2 = Compile[{{rr, _Real, 1}},

and then invoking

t2[Table[r, {r, 1, 100000}]

In this case, I get a timing of around 6 seconds - an improvement, but not what I expected. I don't see any MainEvaluate calls. Is there anything obvious I am missing?

  • 1
    $\begingroup$ Wait a minute, you don't see MainEvaluate in CompilePrint[t1]? Did these functions recently become compilable...? I'm on 10.0.2 and CompilePrint[t1] has 5 lines for me, two of which are MainEvaluate :p $\endgroup$ – Marius Ladegård Meyer Dec 5 '16 at 22:03
  • 1
    $\begingroup$ @MariusLadegårdMeyer, same result on 11.0.1 - I think fred is mistaken. $\endgroup$ – Simon Woods Dec 5 '16 at 22:06
  • $\begingroup$ @SimonWoods, thanks for verifying. fred, that explains the lack of speedup of course. $\endgroup$ – Marius Ladegård Meyer Dec 5 '16 at 22:10
  • $\begingroup$ Hi all - you are actually both right, my apologies. If I type "CompiledFunctionTools``CompilePrint[t1]" I only see the instruction; If I call Needs["CompiledFunctionTools"] first, I see them. My inexperience on Compile hits again. [the inverted accent is not showing well in the in-line code above] $\endgroup$ – Fred Dec 5 '16 at 22:11
  • 1
    $\begingroup$ You should also note that compiled code can only deal with machine numbers and your result is 72 trillion orders of magnitude smaller than the smallest machine real. $\endgroup$ – Simon Woods Dec 5 '16 at 22:26

Since the PDF and CDF of a normal distribution can be expressed analytically as a function of Erfc, you can already speed up your calculation by avoiding the PDF and CDF calculation / lookup at every step:

f[mu_, sigma_] = 
 PDF[NormalDistribution[mu, sigma], r] / CDF[NormalDistribution[mu, sigma], r]


Times @@ Table[f[.2, 1.], {r, 1, 100000}]

(*Out: 3.3929765235*10^-72383065132319 *)

This is already significantly faster than the compiled version:

Times @@ Table[f[.2, 1.], {r, 1, 100000}]; // RepeatedTiming

(* Out: {2.7, Null} *)

An alternative formulation, using the Listable properties of the functions involved, would be:

f2[mu_, sigma_][r_] = 
 PDF[NormalDistribution[mu, sigma], r] / CDF[NormalDistribution[mu, sigma], r]

Times @@ f2[0.2, 1][Range[100000]]; // RepeatedTiming

(* Out: {1.164, Null} *)

Of course, the results are all the same:

Times @@ Table[
  PDF[NormalDistribution[.2, 1], r]/CDF[NormalDistribution[.2, 1], r],
  {r, 1, 100000}
Times @@ f2[0.2, 1][Range[100000]];
Times @@ Table[f[.2, 1.], {r, 1, 100000}];

%%% == %% == %

(* Out: True *)
  • $\begingroup$ This - plus the comments above on the MainEvaluate calls - is very helpful. I have embedded everything into a single compiled function that makes is very fast now. Will try to update the code above. $\endgroup$ – Fred Dec 5 '16 at 22:35
  • $\begingroup$ @fred Something probably got lost in your code when copy-pasting. Try to enclose the code in comments inside backticks like this. I would also urge you to write your own answer combining the insight you gained. $\endgroup$ – MarcoB Dec 5 '16 at 22:41

Following the suggestions above, here's my improved version.

I have found that coding Erfc directly eliminates the MainEvaluate and makes everything remarkably fast. This my improved code:

t3 = Compile[{{rr, _Real, 1}},
        Times @@ Table[
              (Exp[-((.2 + rr[[f]])/1)^2] ((2/\[Pi])^(1/2))  )/
                      (1 Erfc[(.2 - rr[[f]])/(1 2^0.5)]),
          {f, 1, 100000}]
  • $\begingroup$ Fred, Thank you for sharing your solution. Unfortunately, however, I can't seem to get your code to work. Can you show a use example for t3? $\endgroup$ – MarcoB Dec 5 '16 at 23:34
  • $\begingroup$ Sure, I should have been more clear. In my context, the tensor rr has the same length as the table inside the Compile function. As an example, the instructions x=Table[0.2,{i,1,100000}] and t3[x] work. $\endgroup$ – Fred Dec 6 '16 at 2:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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