Optimizing Mathematica Code

Question

I have implemented the whole Baum-Welch Algorithm for training Hidden Markov Models as described in Rabiner's paper Here is my code that calculates the Forward Probabilities`

ForwardProcedure[TM_, EM_, P_, N_, M_, Seq_] := 
 Module[{DP, i, j, DPPart},
  DP = ConstantArray[0, {N, Length@Seq}];

  DP[[1, All]] = Table[P[[i]]*EM[[i, Part[Seq, 1]]], {i, 1, N}];

  For[j = 2, j <= Length@Seq, j++,

   DPPart = DP[[All, j - 1]];

   For[i = 1, i <= N, i++,

    DP[[i, j]] = Dot[ DPPart, TM[[All, i]]]*EM[[i, Part[Seq, j]]];
    ];

   DP[[All, j]] = DP[[All, j]]*1/Total[DP[[All, j]]];
   ];
  DP
  ]

For all who are interested, and/or know the procedure TM = transition Matrix, EM = Emission Matrix, P = initial Probability Distribution, N = Number of states, M = Number of Emissions, Seq = A list containing the observations.

The code runs slow because of a lot of Part and the For statements. Taking under consideration the this fantastic answer

I tried

1. Preallocate large lists and use ConstantArray instead of Table.
2. Used Dot instead of Total.
3. Extracted the previous column only once for each time step only once.

However my ${Mathematica}$ knowledge is lacking, to use some of the advice given in Leonid Shifrin's answer. I was wondering if any of you have ideas on how to make this code run faster -- maybe paralelization? maybe using more kernels?

Additional Information

There are several similar matrix calculating functions that are implemented, I can post their codes as well. I have code that has Modules inside Modules is that good practice?

Update

tm = Table[0.1, {10}, {10}];
em = Table[0.05, {10}, {20}];
p = Table[0.1, {10}];
n = 10;
m = 20;
In[16]:= DP = 
   Do[ForwardProcedure[tm, em, p, n, m, 
     Table[RandomInteger[{1, 10}], 100]], 1]; // RepeatedTiming

Out[16]= {0.0062, Null}

So one iteration costs about 0.0062,a 100000 iteration will cost like 11 minutes.

After all there is a function called TrainHMM Which looks like

TrainHMM[Pi_, TM_, EM_, n_, m_, trainData_] := 
  Module[{newTM = TM, newEM = EM, newPi = Pi, tmpVar, 
    i, α, β, γ, ξ},
   For[i = 1, i <= Length@trainData, i++,
    α = 
     ForwardProcedure[ newTM, newEM, newPi, n, m, trainData[[i]]];
    β = 
     BackwardProcedure[newTM, newEM, newPi, n, m, trainData[[i]]];
    γ = createGammaMatrix[α, β];
    ξ = 
     CreateTransitionMatrixGivenSeq[α, β, newTM, newEM, 
      trainData[[i]]];
    tmpVar = 
     UpdateProbabilities[newPi, newTM, newEM, γ, ξ, 
      trainData[[i]]]; 
    newPi = tmpVar[[1]];
    newTM = tmpVar[[2]];
    newEM = tmpVar[[3]];
    ];
   tmpVar
   ];

Beta and Gamma matrices are calculated similarly to ForwardProcedure in terms of computational complexity they are the same. UpdateProbabilites Has the time complexity $O(n*Length@Seq)$

So overall training would take a 1-1.5 hrs which is well, maybe not so bad, but still, there is always room for improvement. The trainData consists of 100000 lists of length 100.

Answer 1

A version of mmeent's code in a more functional style:

norm[x_] := x / Total[x]

ForwardProcedure[TM_, EM_, P_, N_, M_, Seq_] := Module[{DP, m},
  m = Transpose[EM][[Seq]];
  DP = FoldList[norm[#1.TM #2] &, P First[m], Rest[m]];
  Transpose @ DP]

Answer 2

Here is an optimized version of your function

ForwardProcedure[TM_, EM_, P_, N_, M_, Seq_] := Module[
 {DP, DPPart},
  DP = ConstantArray[0, {N, Length@Seq}];
  DP[[All, 1]] = P*EM[[All, Part[Seq, 1]]];
  Do[
   DPPart = DP[[All, j - 1]].TM;
   DP[[All, j]] = DPPart*EM[[All, Seq[[j]]]];
   DP[[All, j]] = DP[[All, j]]*1/Total[DP[[All, j]]];
   , {j, 2, Length@Seq}
   ];
 DP]

The general idea is to avoid explicit looping constructs, and try to have operations work on lists/arrays as a whole. On my machine this accelerated evaluation by about a factor 3. Feel free to ask about some of the changes.

Update: A version of Simon Woods's code using compile:

norm = Compile[{{x, _Real, 1}}, x/Total[x]]
ForwardProcedure = 
 Compile[
  {
   {TM, _Real, 2},
   {EM, _Real, 2},
   {P, _Real, 1},
   {n, _Integer}, 
   {m2, _Integer},
   {Seq, _Integer, 1}
  },
  With[
   {m = Transpose[EM][[Seq]]},
   Transpose@FoldList[norm[#1.TM #2] &, P First[m], Rest[m]]
  ],
  CompilationOptions -> {"InlineCompiledFunctions" -> True}
 ]