Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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).


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


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.

share|improve this question
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? –  blochwave Jul 30 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 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 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 at 12:33

1 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 = 
    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.

share|improve this answer
Nice work untangling those loops! Any particular reason not to use IdentityMatrix[{J, n}] ? –  Simon Woods Jul 31 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 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. –  blochwave Jul 31 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 at 7:57

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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