2
$\begingroup$

J. M. suggested that I pose the following as a question, a matter arising in relation to the function probit:

SpecialFunctions`Probit; (* force autoload *)
probit = Compile[{{u, _Real}},
                 With[{d = 17/40},
                      If[Abs[u - 1/2] <= d,
                         System`StatisticalFunctionsDump`CompiledProbitCentralMinimax[u], 
                         System`StatisticalFunctionsDump`CompiledProbitAsymptotic[u]]],
                 CompilationOptions -> {"InlineCompiledFunctions" -> True,
                                        "InlineExternalDefinitions" -> True}, 
                 RuntimeAttributes -> {Listable}];

that he gave in his answer to my question Can I use Compile to speed up InverseCDF? . In a comment, he suggested that in order to compare performance times with the code of Henrik Schumacher (also given in an earlier answer of Schumacher to the same question), one should "insert CompilationTarget -> "C", Parallelization -> True after the RuntimeAttributes -> {Listable} part" into the probit code, that is

probit = Compile[{{u, _Real}},
                 With[{d = 17/40},
                      If[Abs[u - 1/2] <= d,
                         System`StatisticalFunctionsDump`CompiledProbitCentralMinimax[u], 
                         System`StatisticalFunctionsDump`CompiledProbitAsymptotic[u]]],
                 CompilationOptions -> {"InlineCompiledFunctions" -> True,
                                        "InlineExternalDefinitions" -> True}, 
                 RuntimeAttributes -> {Listable}, CompilationTarget -> "C",
                 Parallelization -> True];

But when I use a list (as I need to)--say {1/3,1/4}--as the input to the so-augmented probit code I get the error message

CompiledFunction::pext: Instruction 12 in
CompiledFunction[{10,11.3,7900},{_Real},{{3,0,0},{3,0,5}},<<4>>,Evaluate,
LibraryFunction[/Users/Paul/Library/Mathematica/ApplicationData/CCompilerDriver/BuildFolder/
slater-71924/compiledFunction0.dylib,compiledFunction0,{{Real,0,Constant}},Real]]
calls ordinary code that can be evaluated on only one thread at a time.

Without Parallelization -> True, I don't get the error message, but the program then runs more slowly than I would hope. (I'm rather befuddled, since this specific problem did not seem to arise in some related analyses of mine yesterday.)

$\endgroup$
2
$\begingroup$

Good ol' inlining probem. Enforcing correct inlining with With should work. The problem was that the uninlined CompiledProbitCentralMinimax and CompiledProbitAsymptotic enforced calls to the infamous MainEvaluate.

SpecialFunctions`Probit;
probit = With[{
    d = N[17/40],
    cf1 = System`StatisticalFunctionsDump`CompiledProbitCentralMinimax,
    cf2 = System`StatisticalFunctionsDump`CompiledProbitAsymptotic
    },
   Compile[{{u, _Real}},
    If[Abs[u - 0.5] <= d,
     cf1[u],
     cf2[u]
     ],
    RuntimeAttributes -> {Listable},
    CompilationTarget -> "C",
    Parallelization -> True,
    CompilationOptions -> {
      "InlineCompiledFunctions" -> True,
      "InlineExternalDefinitions" -> True
      },
    RuntimeOptions -> "Speed"
    ]
   ];

Here is also J.M.s function AcklamQuantile:

AcklamQuantile = Block[{u}, Compile[{{u, _Real}}, #,
      CompilationTarget -> "C",
      RuntimeAttributes -> {Listable},
      Parallelization -> True,
      RuntimeOptions -> "Speed"
      ] & @@ 
    Hold[With[{a = {-39.69683028665376, 
         220.9460984245205, -275.9285104469687, 
         138.3577518672690, -30.66479806614716, 2.506628277459239}, 
       b = {-54.47609879822406, 161.5858368580409, -155.6989798598866,
          66.80131188771972, -13.28068155288572, 1}, 
       c = {-0.007784894002430293, -0.3223964580411365, \
-2.400758277161838, -2.549732539343734, 4.374664141464968, 
         2.938163982698783}, 
       d = {0.007784695709041462, 0.3224671290700398, 
         2.445134137142996, 3.754408661907416, 1.}}, 
      Which[0.02435 <= u <= 0.97575, 
       With[{v = u - 1/2}, 
         v Fold[(#1 v^2 + #2) &, 0, a]/Fold[(#1 v^2 + #2) &, 0, b]] //
         Evaluate, u > 0.97575, 
       With[{q = Sqrt[-2 Log[1 - u]]}, -Fold[(#1 q + #2) &, 0, c]/
          Fold[(#1 q + #2) &, 0, d]] // Evaluate, True, 
       With[{q = Sqrt[-2 Log[u]]}, 
         Fold[(#1 q + #2) &, 0, c]/Fold[(#1 q + #2) &, 0, d]] // 
        Evaluate]
      ]
     ]
   ];

Speed test against my function cfinv from the linked post (running on a Haswell Quad Core):

T = RandomReal[{0., 1.}, {1000000}];
a = cfinv[T]; // RepeatedTiming // First
b = probit[T]; // RepeatedTiming // First
c = AcklamQuantile[T]; // RepeatedTiming // First
Max[Abs[a - b]]
Max[Abs[a - c]]

0.0525

0.0348

0.0342

2.65592*10^-11

5.59094*10^-9

Well, J.M.'s code is indeed quite a lot faster, in particular at the tails of the distribution (where Newton's method converges slowly).

$\endgroup$
  • $\begingroup$ Interesting, that needs to be injected too... I guess I should edit my answer in the other question. $\endgroup$ – J. M. is in limbo Sep 26 '18 at 18:31
  • $\begingroup$ Misunderstanding--I (stupidly) thought the comment was of Henrik Schumacher. $\endgroup$ – Paul B. Slater Sep 26 '18 at 18:35
  • $\begingroup$ Well, the "correct inlining" code above does now work on lists in my context--but it only seems to give a rather minor (1-2%) improvement over the original cfinv function of Schumacher I've been using. The AcklamQuantile function, I believe (I have to check further) gives a more substantial improvement, but its precision is lower, so maybe I won't convert to that option either. $\endgroup$ – Paul B. Slater Sep 26 '18 at 18:47
  • $\begingroup$ The AcklamQuantile function (with the CompilationTarget -> "C", Parallelization -> True insertion [no "inlining" apparently needed]) gives about a 5% improvement in run time. (The improvement times I give in this comment and the previous one apply to my overall program--and not specifically to the InverseCDF-related commands. So, their improvement times for these commands themselves should be somewhat superior.) $\endgroup$ – Paul B. Slater Sep 26 '18 at 19:04
  • $\begingroup$ Ah, I see. Makes sense: Now other parts of your code are becoming bottlenecks rather than the inverse CDF... $\endgroup$ – Henrik Schumacher Sep 26 '18 at 19:06

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.