# Optimizing Mathematica Code

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?

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.

• About the code-review part: Avoid names with capital letters to prevent conflicts with builtins. N is a builtin. Do would usually be more readable than For. Apr 19, 2017 at 18:01
• This code looks like an obvious candidate for compilation, but it is hard to be more specific without seeing some actual working examples of how you typically use it, together with your timings. Apr 19, 2017 at 18:03
• As for the performance tuning part: please show an example usage with actual parameters, and sufficiently long timing so that is is usable as a benchmark. Apr 19, 2017 at 18:03
• Hi, thanks, edited with some examples of how it will be used, basically this is a part of training algorithm so it will run on sequences of length 100, for at least 100000 examples. The tm and em parameters are in Length and form are almost exactly what I will be using IRL. Apr 19, 2017 at 18:25
• Instead of Table[P[[i]]*EM[[i, Part[Seq, 1]]], {i, 1, N}] use P*EM[[All, Part[Seq, 1]]] Apr 19, 2017 at 18:43

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]


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}
]

• Try to use Compile, like @Leonid Shifring suggested. Apr 19, 2017 at 19:36
• My original version of the function would not compile easily, but the functional style approach certainly does. Apr 19, 2017 at 20:15
• @mmeent Current compiled version will use MainEvaluate for norm, which can slow things down. See <<CompiledFunctionTools  and CompilePrint[ForwardProcedure] for virtual machine code. Apr 19, 2017 at 20:24
• You can inline norm into your code i.e. substitute norm[#1.TM #2] & with ( (#1.TM #2)/Total[#1.TM #2]) &. This will make compiled code even faster. Apr 19, 2017 at 20:32
• @Shadowray indeed Apr 19, 2017 at 20:32