Fast way to build a table that refers to most recent addition?

Question

I'm trying to construct a large table where new entries depend on the most recent entries. It's quite slow: takes 2 mins and I'm embedding this into a function that does it repeatedly so any reduction in time is a big savings. I think this issue is with "AppendTo" and I've seen repeated references on here to use Reap/Sow, but I'm having trouble implementing that because the next entries in the table are conditional on the previous. Here is an example of the code I'm working with right now.

SeedRandom[1234];

TransMark = Table[RandomReal[], {i, 7}, {j, 7}]; (*Just an example of a similar matrix to what I have*)

RandDraw = Table[RandomReal[], {i, 10000}, {j, 576}];

WorkerMat = Table[{{2, 1}}, {10000}]; 

Timing[Do[
   If[RandDraw[[i, j]] <= 
     TransMark[[All, 1]][[Last[WorkerMat[[i]]][[1]]]],
    If[Last[WorkerMat[[i]]][[1]] == 1, 
     Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
    AppendTo[WorkerMat[[i]], {1, Tenure}],
    If[RandDraw[[i, j]] <= 
      TransMark[[All, 2]][[Last[WorkerMat[[i]]][[1]]]],
     If[Last[WorkerMat[[i]]][[1]] == 2, 
      Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
     AppendTo[WorkerMat[[i]], {2, Tenure}],
     If[RandDraw[[i, j]] <= 
       TransMark[[All, 3]][[Last[WorkerMat[[i]]][[1]]]],
      If[Last[WorkerMat[[i]]][[1]] == 3, 
       Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
      AppendTo[WorkerMat[[i]], {3, Tenure}],
      If[RandDraw[[i, j]] <= 
        TransMark[[All, 4]][[Last[WorkerMat[[i]]][[1]]]],
       If[Last[WorkerMat[[i]]][[1]] == 4, 
        Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
       AppendTo[WorkerMat[[i]], {4, Tenure}],
       If[
        RandDraw[[i, j]] <= 
         TransMark[[All, 5]][[Last[WorkerMat[[i]]][[1]]]],
        If[Last[WorkerMat[[i]]][[1]] == 5, 
         Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
        AppendTo[WorkerMat[[i]], {5, Tenure}],
        If[
         RandDraw[[i, j]] <= 
          TransMark[[All, 6]][[Last[WorkerMat[[i]]][[1]]]],
         If[Last[WorkerMat[[i]]][[1]] == 6, 
          Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
         AppendTo[WorkerMat[[i]], {6, Tenure}],
         If[Last[WorkerMat[[i]]][[1]] == 7, 
          Tenure = Last[WorkerMat[[i]]][[2]] + 1, Tenure = 1];
         AppendTo[WorkerMat[[i]], {7, Tenure}];]]]]]] , {j, 1, 
    576}, {i, 1, 10000}];
 ];

Which takes over 2 mins on my machine. My apologies if this is a repeat question, but I've been looking for a similar post for many hours and have not discovered one. Thanks in advance!

Answer 1

There are two issues with your code: First, it is very procedural with lots of If statements. Mathematica is slow at processing procedural code. Good news: Procedural code in Mathematica can usually be compiled with Compile in order to speed up the execution.

The second and major flaw of your code however is that you use AppendTo quite often. Every time you use it, a list has to be copied which makes the code very memory bound and of computational complexity $O(N^2)$ for a list with $N$ elements (and thus slow).

Fortunately, the size of WorkingMat is known right from the start and Last[WorkerMat[[i]]] refers to WorkerMat[[i,j]] and Append[Last[WorkerMat[[i]]],value] can be replaced by WorkerMat[[i,j+1]] = value. Moreover, one can make the j-loop the outer loop and parallelize it. Afterwards, one obtains a code which is nicely compilable. Replacing the nested If construction by a lookup for the first entry in a certain list which is greater or equal a certain value, the result is the CompiledFunction cf below (I am skipping a significant amount of refactorization steps here.)

cf = Compile[{{rand, _Real, 1}, {TransMark, _Real, 2}},
  Block[{r, a, b, k},
   a = Table[0, {i, 1, Length[rand] + 1}, {j, 1, 2}];
   a[[1, 1]] = 2;
   a[[1, 2]] = 1;
   Do[
    r = Compile`GetElement[rand, j];
    b = Compile`GetElement[a, j, 1];
    k = 1;
    While[r > Compile`GetElement[TransMark, b, k] && k < 7, k++];
    a[[j + 1, 1]] = k;
    a[[j + 1, 2]] = If[b == k, Compile`GetElement[a, j, 2] + 1, 1];
    , {j, 1, Length[rand]}];
   a
   ],
  RuntimeAttributes -> {Listable},
  Parallelization -> True,
  CompilationTarget -> "C",
  RuntimeOptions -> "Speed"
  ]

Btw., the Compile`GetElement statements in the code above are read-only instructions that replace Part. I think it instructs the C code generator to skip certain bound checks. This can improve the performance considerably but makes it less "safe".

Generating some data...

SeedRandom[1234];
TransMark = RandomReal[{0, 1}, {7, 7}];
n = 10000;
m = 576;
RandDraw = RandomReal[{0, 1}, {n, m}];

... and running the code on the data

WorkerMat2 = cf[RandDraw, TransMark]; // AbsoluteTiming//First

0.05216

Your original code took about 86 seconds on my machine, so that's a speed-up of factor 1600 or so.

Answer 2

This is substantially faster and might be equivalent. If it is not, I may well have made mistakes in translating and, as was noted, a minimal working example is a useful commodity.

SeedRandom[1234];
TransMark = RandomReal[{0, 1}, {7, 7}];
maxi = 10000;
maxj = 576;
RandDraw = RandomReal[{0, 1}, {maxi, maxj}];
WorkerMat = Table[{{2, 1}}, {maxi}];

AbsoluteTiming[
 wm2 = Table[
   last = Last[WorkerMat[[i]]];
   rndi = RandDraw[[i]];
   ltable = Table[
     rndij = rndi[[j]];
     {l1, l2} = last;
     tm1 = TransMark[[l1]];
     tenure = 1;
     Which[rndij <= tm1[[1]], k = 1,
      rndij <= tm1[[2]], k = 2,
      rndij <= tm1[[3]], k = 3,
      rndij <= tm1[[4]], k = 4,
      rndij <= tm1[[5]], k = 5,
      rndij <= tm1[[6]], k = 6,
      True, k = 7
      ];
     If[l1 == k, tenure = l2 + 1];
     last = {k, tenure}
     , {j, maxj}];
   Developer`ToPackedArray[ltable]
   , {i, maxi}];]

(* Out[374]= {27.892029, Null} *)

wm3 = Join[WorkerMat, wm2, 2];

That last should give the original result.

I'm sure a further speed gain could be had using Compile.