# Why is my compiled code slower?

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.

• Your compiled version doesn't run for me - looks like some syntax errors with the brackets? e.g. (p/Length[NRR[[j]]) should be (p/Length[NRR[[j]]]). Also can you post some example inputs X,J,n,p and NRR? – dr.blochwave Jul 30 '14 at 10:53
• The example code is missing braces and is syntactically not correct. In its current form the question cannot be answered. – halirutan Jul 30 '14 at 11:19
• way before messing with Compile learn to utilize mathematica object oriented programming and wealth of built in functions. For a start look up InentityMatrix.. – george2079 Jul 30 '14 at 12:18
• @george2079...that is, IdentityMatrix... :) @OP: Note that K has a special use in Mathematica; evaluate ?K to find out. – Michael E2 Jul 30 '14 at 12:33

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:= 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[]]
(* 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.

• Nice work untangling those loops! Any particular reason not to use IdentityMatrix[{J, n}] ? – Simon Woods Jul 31 '14 at 11:31
• @SimonWoods The particular reason that "I fer-got it existed...". Feel free to edit that in if it works. – Verbeia Jul 31 '14 at 11:54
• Good work, +1! Should all be compilable too, according to the list here: mathematica.stackexchange.com/questions/1096/…, though not necessarily without calls to MainEvaluate. – dr.blochwave Jul 31 '14 at 18:44
• Thanks a bunch for the improvements. I will incorporate them into my program straight away. Just to clarify what I should probably have included before, J>n (indeed much greater since I do not need the full tensor hKer for my purposes), Length[NRR]=J, and Length[NRR[[j]]]=4 for most but not all j. – Joshua Cooperman Aug 1 '14 at 7:57