3
$\begingroup$

For some table tab having N rows and n columns, I need to go over all its rows $j$ and columns $f$ and sum some expression depending on the tabs elements:

$$ \text{count} = \sum_{j = 1}^{N}\sum_{f_{1} = 1}^{n-1}\sum_{f_{2} = f_{1}+1}^{n}\text{expression}(\text{tab}_{j,f_{1}},\text{tab}_{j,f_{2}}) $$

Here is my code assuming that the expression is just some boolean condition:

ncols = 2;
tab = RandomReal[{0, 1}, {10^7, ncols}];
comp1 = Hold@Compile[{{tab, _Real, 2}}, Module[{el1, el2, count},
      count = 0;
      (*Cycle over rows*)
      Do[
       (*Cycle over f1*)
       Do[
        el1 = Compile`GetElement[tab, j, f1];
        (*Cycle over f2*)
        Do[
         el2 = Compile`GetElement[tab, j, f2];
         (*adding to count*)
         count += Boole[el1*el2 > 0.4 && el1^2 + el2^2 < 0.9]
         , {f2, f1 + 1, ncols, 1}]
        , {f1, 1, ncols - 1, 1}]
       , {j, 1, Length[tab], 1}];
      count], CompilationTarget -> "C", RuntimeOptions -> "Speed"] /. 
   OwnValues@ncols // ReleaseHold

The code structure is convenient for my real task (see below), so I want to preserve it.

Another version uses For instead of Do:

comp3 = Hold@
    Compile[{{tab, _Real, 2}}, Module[{el1, el2, count, j, f1, f2},
      count = 0;
      (*Cycle over rows*)
      For[j = 1, j <= Length[tab], j++,
       For[f1 = 1, f1 <= ncols - 1, f1++,
        el1 = Compile`GetElement[tab, j, f1];
        (*Cycle over f2*)
        For[f2 = f1 + 1, f2 <= ncols, f2++,
         el2 = Compile`GetElement[tab, j, f2];
         (*adding to count*)
         count += Boole[el1*el2 > 0.4 && el1^2 + el2^2 < 0.9]]]];
      count], CompilationTarget -> "C", RuntimeOptions -> "Speed"] /. 
      OwnValues@ncols // ReleaseHold

For n = 2, the same summation may be done using this simple code:

comp2 = Hold@Compile[{{tab, _Real, 2}}, Module[{el1, el2, count},
      count = 0;
      Do[
       el1 = Compile`GetElement[tab, j, 1];
       el2 = Compile`GetElement[tab, j, 2];
       count += Boole[el1*el2 > 0.4 && el1^2 + el2^2 < 0.9]
       , {j, 1, Length[tab], 1}];
      count], CompilationTarget -> "C", RuntimeOptions -> "Speed"] /. 
   OwnValues@ncols // ReleaseHold

comp2 is much faster than comp1, faster than comp3, while comp3 is faster than comp1

comp1[tab] // AbsoluteTiming
comp2[tab] // AbsoluteTiming
comp3[tab] // AbsoluteTiming

{0.0745858,162135}

{0.0325867,162135}

{0.048013,162135}

The same timing hierarchy between comp1, comp3 holds if increasing the number of columns in tab.

Could you please tell me what is the reason for the observed slowdown of comp3/comp1 compared to comp2, and of comp1 compared to comp3 (the only difference is the replacement Do->For), and how to improve comp1/comp3 if preserving a similar structure with the cycles?

Edit

After the changes from the answer by Michael E2, the timings are

{0.0548117,162135}

{0.0360111,162135}

{0.04748,162135}

Still, the Do version is slower on my machine.

P.S.

I need this structure because in my real problem (illustrated in a simplified way by this question):

  1. I do not just sum the elements but rather some values obtained using expressions involving these elements,

  2. the expressions are computationally expensive,

  3. and there are some additional conditions inside the summation (say, I may want to break the summation over the columns because of some reason).

Update Consider the following code:

comp4 = Compile[{{tab, _Real, 2}}, 
  Module[{el1, el2, count, ncols, f1max},
   count = 0;
   ncols = (Length[Compile`GetElement[tab, 1]]*2)/2;
   f1max = ncols - 1;
   (*Cycle over rows*)
   Do[(*Cycle over f1*)
    Do[el1 = Compile`GetElement[tab, j, f1];
     (*Cycle over f2*)
     Do[el2 = Compile`GetElement[tab, j, f1 + f2];(*!*)
      (*adding to count*)
      count += Boole[el1*el2 > 0.4 && el1^2 + el2^2 < 0.9]
      , {f2, ncols - f1}],(*!*)
     {f1, f1max}],(*!*)
    {j, Length[tab]}];
   count], CompilationTarget -> "C", RuntimeOptions -> "Speed"]
comp1[tab] // AbsoluteTiming
comp4[tab] // AbsoluteTiming

Note the string (Length[CompileGetElement[tab, 1]]*2)/2;. Normally, it is just Length[CompileGetElement[tab, 1]]. But the stupid code probably treats it without the simplification. As a result, comp4 is evaluated 3 times slower than comp1. Maybe some analog of Simplify would work, but Simplify itself is not compilable.

$\endgroup$
6
  • 1
    $\begingroup$ Why don't you use Sum and Map and Plus/Total? $\endgroup$
    – lericr
    Jul 14 at 12:46
  • $\begingroup$ @lericr : because in my real problem I do not just sum the elements of the table but some values obtained using expressions which use these elements. I will add the clarification. $\endgroup$ Jul 14 at 12:48
  • 1
    $\begingroup$ Sum can manage the indices. $\endgroup$
    – lericr
    Jul 14 at 13:00
  • $\begingroup$ @lericr : currently, I use this Do structure for the code from this question: mathematica.stackexchange.com/questions/287399/… (the actual relevant code is acceptanceOptimized). I have doubts about whether Sum may be used there. That's why I prefer (in the lack of other ways) to keep Do in this example. $\endgroup$ Jul 14 at 13:03
  • $\begingroup$ (1) In compiling to WVM, Mma constructs have limitations that facilitate translation to (optimized) C-like code. (2) In my own head-to-head test of triply nested Do[] v. For[], For[] still lags behind Do[] by quite a bit (5X slower in WVM, 30K X slower in C!). A moderately slow calculation in the body of the loop might make that overhead insignificant. My conclusion is that compiled For[] is still slower than Do[]. $\endgroup$
    – Michael E2
    Jul 14 at 16:17

1 Answer 1

2
$\begingroup$

The culprit seems to be an IteratorCountI[] instruction when Do[] has a variable starting point. In this case in the inner loop {f2, f1+1, ncols}. I have no idea why this should be so slow.

ncols = 2;
tab = RandomReal[{0, 1}, {10^7, ncols}];
comp1 = Hold@Compile[{{tab, _Real, 2}},
      Module[{el1, el2, count},
       count = 0;
       (*Cycle over rows*)
       Do[
        (*Cycle over f1*)
        Do[
         el1 = Compile`GetElement[tab, j, f1];
         (*Cycle over f2*)
         Do[
          el2 = Compile`GetElement[tab, j, f1 + f2]; (* ! *)
          (*adding to count*)
          count += Boole[el1*el2 > 0.4 && el1^2 + el2^2 < 0.9],
          {f2, ncols - f1}],                         (* ! *)
         {f1, "ncols-1"}],                           (* ! *)
        {j, Length[tab]}];
       count
       ],
      CompilationTarget -> "C", RuntimeOptions -> "Speed"] /.
     OwnValues@ncols /.
   "ncols-1" -> ncols - 1 // (* a slight efficiency *)
  ReleaseHold

Timings:

comp1[tab] // RepeatedTiming
comp2[tab] // RepeatedTiming
comp3[tab] // RepeatedTiming
(*
  {0.010221,   162488}
  {0.00751259, 162488}
  {0.0145402,  162488}
*)
$\endgroup$
9
  • $\begingroup$ Thanks! And without this fix, you also had a significant slowdown of comp1? $\endgroup$ Jul 14 at 18:04
  • 1
    $\begingroup$ @JohnTaylor For your original comp1, it was about 6X times slower than comp2 (0.044 s) and 3X times slower than comp3. I think the "ncols-1" trick sped it up to 5X slower than comp2 (0.036 s). For some reason ncols - 1 is optimized in For[] but not in Do[]. (It's computed each iteration in Do[] as 2 - 1, instead of compiled as 1). $\endgroup$
    – Michael E2
    Jul 14 at 18:21
  • $\begingroup$ And the CPU is M1? $\endgroup$ Jul 14 at 19:19
  • 1
    $\begingroup$ Apparently division is parsed as returning a Real, which makes ncols be typed as a Real. You can declare types with initialization in Module: Module[{el1 = 0., el2 = 0., count = 0, ncols = 0, f1max = 0},...]. This adds minimal initialization overhead if done outside loops. You also need to explicitly cast the division to an integer: ncols = Floor[(Length[Compile`GetElement[tab, 1]]*2)/2];. BTW, I don't know much about Compile. It's all trial and error. We used to have a few users who seemed to know some of the internals. Like optimization... $\endgroup$
    – Michael E2
    Jul 15 at 13:53
  • 2
    $\begingroup$ ....This is not optimized the way you might expect: Experimental`OptimizeExpression[(Length[Compile`GetElement[tab, 1]]*2)/2 // Unevaluated, "OptimizationLevel" -> 2]. This might be why, or it might be that Length[..]*2/2 is type Real and Length[..] is type Integer; so no simplification. -- Passing a scalar argument to the CompiledFunction has minimal overhead. Why not pass ncols = Length[tab[[1]]]/x? I'm not sure the ncols-1 trick is that robust. $\endgroup$
    – Michael E2
    Jul 15 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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