Compile nested Sums

Question

I want to compile an expressions that contains nested Sum-expressions. A simple example that gives me problems is:

sumsum = Compile[{{n, _Integer}}, Sum[1/Sum[i + j, {j, n}], {i, n}]]

Note: This is not the problem I want to solve, it is just a simple example code, where the problem, which is described below, occurs.

While the expression compiles, using CompilePrint shows that the whole evaluation does not use compiled code, because there clearly is a call to MainEvaluate when the first call to Sum appears:

Needs["CompiledFunctionTools`"]
CompilePrint[sumsum]

A quick fix would be to use Plus@@Table instead of Sum, but this combination is in general slower than using Sum. Indeed

Plus @@ Table[1/Sum[i + j, {j, n}], {i, n}]

Does not show any call to MainEvaluate and runs faster than the uncompiled code.

Any suggestions on how to get nested sums to compile correctly? (once again, I am interested in a general solution to compiling nested sums, not a specific solution to the example stated above)

Answer 1

Let me elaborate my comment into an answer. To make the nested Sum compiled, let's first have a close look at the compiling result of code containing one Sum:

sum = Compile[{{n, _Integer}}, 1/Sum[3.141 + j, {j, n}]];
Needs["CompiledFunctionTools`"]
p1 = CompilePrint@sum

It's not hard to notice that Sum is actually translated into a loop by Compile. It's reasonable to guess, if we express the summation with a loop instead of Sum, it will be finally translated into the same result, and it's true:

sum2 = Compile[{{n, _Integer}}, Module[{k = 0.}, 1/(Do[k += 3.141 + j, {j, n}]; k)]];
p2 = CompilePrint@sum2;
Grid[{{p1, p2}}, Frame -> All]

You can see the compiling results are almost same.

Manually replacing Sum with Module[…, Do[…]] is tedious, let's make use of code generation technique to facilitate the replacement:

sum[expr_, iter_] := Module[{k = 0.}, Do[k += expr, iter]; k];

sum2test = Hold@Compile[{{n, _Integer}}, 1/sum[3.141 + j, {j, n}]] //. DownValues@sum // 
   ReleaseHold;

sum2test@1000 == sum2@1000
(* True *)

Then since the Module[… Do[… ]] algorithm is closer to the the one adopted inside Compile, it's reasonable to guess that Compile can translate it more easily i.e. it may be compiled successfully in a situation that Sum fails, for example, nested Sum, and again, it's true:

sumsum2 = Hold@Compile[{{n, _Integer}}, sum[1/Sum[i + j, {j, n}], {i, n}]] //. 
DownValues@sum // ReleaseHold;

CompilePrint@sumsum2    

Pictured by Simon Wood's shadow.

Now the code is completely compiled and I believe it had the same performance as the compiled nested Sum (if it was correctly compiled).

Another acceptable solution is

sumsum3 = Compile[{{n, _Integer}}, Total@Table[1/Sum[i + j, {j, n}], {i, n}]];

It's a little faster than sumsum2:

sumsum3[10^4] // AbsoluteTiming
sumsum2[10^4] // AbsoluteTiming

{4.073000, 0.000109848}
{4.156000, 0.000109848}

However, what I really want to say is, I believe a nested Sum is not the fastest solution for your original problem at all, considering the facts that most of the arithmetic functions in Mathematica are Listable and taking advantage of the vectorization is one of the most important approaches to speed up codes even inside Compile. Here's a vectorized version of your toy code:

sumsum0 = Compile[{{n, _Integer}}, Total[1/(n Range@n + Total@Range@n)]];
sumsum0[10^4] // AbsoluteTiming

{0.001000, 0.000109848}