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
- Preallocate large lists and use
ConstantArray
instead ofTable
. - Used
Dot
instead ofTotal
. - 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 Module
s inside Module
s 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.
code-review
part: Avoid names with capital letters to prevent conflicts with builtins.N
is a builtin.Do
would usually be more readable thanFor
. $\endgroup$Table[P[[i]]*EM[[i, Part[Seq, 1]]], {i, 1, N}]
useP*EM[[All, Part[Seq, 1]]]
$\endgroup$