Why is my compiled code slower?

Question

I am trying to make a Mathematica program run faster. Here is part of the original version. It takes as input the integers X (of order 100), J (of order 10000), and n (of order 100), the real number p (of order 1), and the rank 2 list NRR of integers (or order 10000 by 4).

 hKer=Table[0,{k,1,X},{j,1,J},{i,1,n}];
    For[j=1,j<n+1,j++,
    For[i=1,i<n+1,i++,
    hKer[[1,j,i]]=KroneckerDelta[j,i]]];
    For[k=2,k<X+1,k++,
    For[j=1,j<J+1,j++,
    For[i=1,i<n+1,i++,
    hKer[[k,j,i]]=N[(1-p)hKer[[k-1,j,i]]+(p/Length[NRR[[j]]])Sum[hKer[[k-1,NRR[[j,l]],i]],{l,1,Length[NRR[[j]]]}]]]]];

Searching the internet for advice led me to try compiling this program. Copying an example on stackoverflow, I arrived at the following version.

fhKer=Compile[{{X,_Integer},{J,_Integer},{n,_Integer},{p,_Real},{NRR,_Integer,2}},
Module[{hKer},
hKer=Table[0.,{k,1,X},{j,1,J},{i,1,n}];
Do[hKer[[1,j,i]]=If[j==i,1,0],{j,1,n},{i,1,n}];
Do[hKer[[k,j,i]]=N[(1-p)K[[k-1,j,i]]+(p/Length[NRR[[j]]])Sum[hKer[[k-1,NRR[[j,l]],i]],{l,1,Length[NRR[[j]]]}],{k,2,X},{j,1,J},{i,1,n}];
hKer]]

The compiled program runs more slowly than the original version. I checked that the compiled program does not contain a MainEvaluate, which was put forth as a potential problem on stackexchange. I do not currently have any good ideas either for making the compiled program run more quickly or for modifying the original version to run more quickly. I would greatly appreciate any suggestions.

Answer 1

In this instance, compiling your triply nested For loops is unlikely to be the most effective way to speed up your code. I may have misunderstood what you are trying to do, being allergic to nested loops. However, two obvious speedups occur to me:

1. Define your starting value in a more efficient way
2. Use functional programming (specifically NestList) to build your tensor.

Here is a toy example showing what I mean. First, let's set the sizes. I assume that this is what you mean. Note that NRR needs to have elements that are never greater than J. I also believe that J and n must be equal, but perhaps I've misunderstood.

 J = 100; X = 100; n = 100; p = 0.6;    
 NRR = RandomInteger[{1, 10}, {100, 4}];

Here is your code. It is pretty slow on my machine.

In[229]:= hKer = Table[0, {k, 1, X}, {j, 1, J}, {i, 1, n}];
For[j = 1, j < n + 1, j++, 
  For[i = 1, i < n + 1, i++, 
   hKer[[1, j, i]] = KroneckerDelta[j, i]]];
For[k = 2, k < X + 1, k++, 
   For[j = 1, j < J + 1, j++, 
    For[i = 1, i < n + 1, i++, 
     hKer[[k, j, i]] = 
      N[(1 - p) hKer[[k - 1, j, i]] + (p/Length[NRR[[j]]]) Sum[
          hKer[[k - 1, NRR[[j, l]], i]], {l, 1, 
           Length[NRR[[j]]]}]]]]]; // Timing

(*{14.724984, Null}*)

Now let's try a functional version. First, I'm pretty sure that NRR is a matrix, so Length[NRR[[j]] ] is the same regardless of j. (But again, perhaps I have misunderstood.) If that's true, you can save yourself some hassle by defining it once.

lenNRR = Length[NRR[[1]]]
(* 4 *)

Now, here is my version, which builds up the tensor matrix by matrix recursively. The main things to note are that (1) Join[IdentityMatrix[n], ConstantArray[0, {J - n, n}]] stands in for your first set of (doubly) nested loops setting up the starting value, and (2) that you can extract multiple parts of the previous result (hKer[[k-1]] in your notation) by taking each row of NRR at once into the Part ([[]]) function. I use the With construct to avoid getting my pure function Slot (#) notation confused.

betterhKer = 
   NestList[
    With[{prev = #1}, (1 - p) prev + 
       Total /@ ((p/lenNRR)*(prev[[#]] & /@ NRR))] &, 
    Join[IdentityMatrix[n], ConstantArray[0, {J - n, n}]] , 
    X]; // Timing

(* {0.086085, Null} *)

171 times faster, and the results are the same.

This is an example where learning the idiom of Mathematica and its functional and matrix-based features will save you time and hassle.