1
$\begingroup$

I was trying to find out why the following function is rather slow.

i4p[a_, m1_, m2_, m3_, m4_] := (
If[Mod[m1 + m2 + m3 + m4, 2] == 1,
0,
z = 1/(2*(1 - a)); k = (m1 + m2 + m3 + m4)/2 ;
summation = 0;
Do[
 Do[
   Do[
     Do[
       Do[
         Do[

           If[idx21 + idx31 + idx41 + idx32 + idx42 + idx43 > k, 
            Continue[]];

           summation = 
            summation + (Pochhammer[-m1, idx21 + idx31 + idx41]*
                Pochhammer[-m2, idx21 + idx32 + idx42]*
                Pochhammer[-m3, idx31 + idx32 + idx43]*
                Pochhammer[-m4, idx41 + idx42 + idx43])/(idx21!*
                idx31!*idx41!*idx32!*idx42!*idx43!*
                Pochhammer[1/2 - k, 
                idx21 + idx31 + idx41 + idx32 + idx42 + idx43])*
              z^(idx21 + idx31 + idx41 + idx32 + idx42 + idx43)
           , {idx21, 0, k}];
         , {idx31, 0, k}];
       , {idx41, 0, k}];
     , {idx32, 0, k}];
   , {idx42, 0, k}];
 , {idx43, 0, k}];
Sqrt[a]*(-1/2*z)^(-k) Gamma[k + 1/2]*summation
]);

So I tried adding // RuntimeTools`Profile and it showed me as largest contributions:

Calls   Time    Code
1771561 56.557  If[idx21+idx31+idx41+idx32+idx42+idx43>k,Continue[]]
1771561 41.451  idx21+idx31+idx41+idx32+idx42+idx43>k
1771561 32.638  idx21+idx31+idx41+idx32+idx42+idx43
1763553 6.022   Continue[]

of a overall 69 seconds runtime. From this I take it that Mathematica is very efficient about caching the calls to Pochhammer, which is great, but seems to be very slow when adding up simple integers. How can this be? 2Million integer additions should be done in a hartbeat on modern processors. I appreciate any hints on how to speed this up.

$\endgroup$
5
  • 1
    $\begingroup$ Do is a very awkward construction in Mathematica and not in the functional programming spirit or style. Rewrite all your terms in the functional programming style and your code will be much faster. $\endgroup$ – David G. Stork Mar 5 '15 at 3:17
  • $\begingroup$ Can you give a set of valid {a, m1, m2, m3, m4} for test? $\endgroup$ – xzczd Mar 5 '15 at 3:19
  • 1
    $\begingroup$ This is not a performance problem. This is a coding problem. In addition, adding some description of what your after with this might prompt a much simpler construct than this monstrosity. As DS commented, coding this properly should result in greater speed, but I'd venture there's a more direct method vs just speeding up this. $\endgroup$ – ciao Mar 5 '15 at 3:35
  • $\begingroup$ you realize that Continue only exits the innermost loop. Using the procedural Do construct instead of that if, every level should be bounded by that criteria, ie, `{idx21, 0, k-(idx31 + idx41 + idx32 + idx42 + idx43)} $\endgroup$ – george2079 Mar 5 '15 at 4:36
  • $\begingroup$ The 6th root of 1771561 says k==11. Your profile shows that 99.5% of the time the bound is exceeded. You aren't immediately exiting any Do and skipping all the remaining iterations that it would Do. I had thought of pushing it all inside a Sum[expr, {idx43,1,11}, {idx42,1,11}, ...{idx21,1,11}] but with the 99.5% it seems much more likely to Break out of each loop as soon as the bound is exceeded. $\endgroup$ – Bill Mar 5 '15 at 4:40
3
$\begingroup$

First of all, we should be aware that given loop has number of operations proportional to $k^6$, so the time grows quite fast with $k$. I think the best way to implement

$i_{21} + i_{31} + i_{41} + i_{32} + i_{42} + i_{43} \le k$

is to sum over (you can insert it in your code instead of your arguments of Do)

{idx21, 0, k - idx43 - idx42 - idx32 - idx41 - idx31}];,
{idx31, 0, k - idx43 - idx42 - idx32 - idx41}];,
{idx41, 0, k - idx43 - idx42 - idx32}];,
{idx32, 0, k - idx43 - idx42}];,
{idx42, 0, k - idx43}];,
{idx43, 0, k}];

This is a major speedup in your loop: it reduces hypercube $k \times k \times k \times k \times k \times k$ to much smaller loop iterations (I'm a bit lazy to provide exact answer, but you can imagine it as geometric figure). The remaining part of my answer are minor speedup boosts.

Next thing you can do is remove next line from your procedure

If[idx21 + idx31 + idx41 + idx32 + idx42 + idx43 > k,  Continue[]];

(it will never trigger). Then you can wrap your code in

Module[{z, k, summation},
(*your procedure is here*)
]

To omit change of variables z, k, summation etc in a global context. Ok, it's not a speedup, but it can save your life if you have, for example, Global`z variable. There a lot minor improvements, but I think they are not worth it. To give you an answer, I have rewritten your function in functional style (mapping on tables instead of iterating by do-loops), and it worked faster (but not dramatically) before I've applied Major Trick, but in final version this didn't made any observable difference. To sum up, the point was in your procedure, and I believe same performance would be observed with different language). Still, time should grow like $C\, k^6$, and we only made $C$ considerably smaller.

Here is your code edited, for reference:

i4p[a_, m1_, m2_, m3_, m4_] /; (Mod[m1 + m2 + m3 + m4, 2] == 1) := 0
i4p[a_, m1_, m2_, m3_, m4_] :=
Module[{z, k},
  z = 1/(2*(1 - a));
  k = (m1 + m2 + m3 + m4)/2;
  summation = 0;
  Do[Do[Do[Do[Do[Do[
     summation =
         summation +
         (Pochhammer[-m1, idx21 + idx31 + idx41]*
         Pochhammer[-m2, idx21 + idx32 + idx42]*
         Pochhammer[-m3, idx31 + idx32 + idx43]*
         Pochhammer[-m4, idx41 + idx42 + idx43])/
         (idx21!*idx31!*idx41!*idx32!*idx42!*idx43!*
         Pochhammer[1/2 - k,idx21 + idx31 + idx41 + idx32 + idx42 + idx43])*
         z^(idx21 + idx31 + idx41 + idx32 + idx42 + idx43),
     {idx21, 0, k - idx43 - idx42 - idx32 - idx41 - idx31}];,
     {idx31, 0, k - idx43 - idx42 - idx32 - idx41}];,
     {idx41, 0, k - idx43 - idx42 - idx32}];,
     {idx32, 0, k - idx43 - idx42}];,
     {idx42, 0, k - idx43}];,
     {idx43, 0, k}];
  Sqrt[a]*(-1/2*z)^(-k) Gamma[k + 1/2]*summation
]
$\endgroup$
1
  • $\begingroup$ Slightly more compact, but no faster, is to use Sum (but not NSum) in place of the nested Do loops. In any case, effective solution. $\endgroup$ – bbgodfrey Mar 6 '15 at 13:09

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.