# My Baum-Welch algorithm is very slow. Is it due to Mathematica?

I've also posted this question here, since I'm not sure which community is more adequate for this type of question.

I've programmed imperatively the Baum-Welch algorithm for a 2 state hidden Markov model (HMM), as is explained in this book. When I try it with a size sample of 100, I usually get a decent estimate after 18 seconds. However, when I use a sample size of 500, it takes me between 290-350 seconds to get an improved estimation of the HMM parameters.

I get the impression the algorithm should not perform this badly. However I'm not sure why it does.

One thing I've discovered is that Mathematica is more functional oriented, so when I program iteratively I'm losing some computation 'power'. But is this really enough for this bad performance??

What is your experience with this algorithm?

I'm using as stopping criterion when the norm of the difference between the current estimate and the previous one is less than $10^{-3}$. This criterion just seemed simple enough. I'm also trying to use other stopping criteria, but it doesn't seems to alter much. The problem is each iteration of the algorithm, in a sample size of 500, takes almost 20 seconds more or less... This worries me.

Any help would be appreciated.

Here is my code.

EM1[y_, ite_] :=
Module[{n, listainicial, teta0, teta1, teta2, likeham, likehamtotal,
pimattotal, sigmatotal, tetatotal, θk1, θk2,
listak, listbk, fikyj, iterations},
n = Length[y] - 1;
listainicial = N[Table[0, {j, 1, n}]];
pimattotal = listainicial;
tetatotal = listainicial;
sigmatotal = listainicial;
likehamtotal = listainicial;

likeham[pimat_, fiyj_] := Module[{pizero, Pmat, counter, mult},
pizero = pimat[[1]];
Pmat = {pimat[[2]], pimat[[3]]};
counter = 1;
mult = IdentityMatrix[2];
While[counter <= n,
mult =
Dot[mult,
Dot[Pmat, {{fiyj[[1, counter]], 0}, {0, fiyj[[2, counter]]}}]];
counter = counter + 1;
];
Dot[Dot[pizero, mult], {1, 1}]
];

iterations = 1;
While[iterations <= ite,
(*initialization*)

pim0 = RandomVariate[UniformDistribution[{0, 1}]];
pim1 = RandomVariate[UniformDistribution[{0, 1}]];
pim2 = RandomVariate[UniformDistribution[{0, 1}]];

pimatk = N[{{pim0, 1 - pim0}, {pim1, 1 - pim1}, {pim2, 1 - pim2}}];
fikyj = N[Table[Table[1/n, {j, 1, n}], {i, 1, 2}]];

ak = N[Table[Table[0.1, {j, 1, n}], {i, 1, 2}]];
bk = N[Table[Table[1, {j, 1, n}], {i, 1, 2}]];
θk1 = {0.1, 0.1, 0.1};
θk2 = {0.5, 0.5, 0.5};

listak = listainicial;
listbk = listainicial;

While[Norm[θk1 - θk2] > 10^-3,

θk1 = θk2;

ak[[1, 1]] = pimatk[[1, 1]]*fikyj[[1, 1]];
ak[[2, 1]] = pimatk[[1, 2]]*fikyj[[2, 1]];

ak[[1, 1]] = ak[[1, 1]];
ak[[2, 1]] = ak[[2, 1]];

(*Induction*)
counterE = 1;
While [counterE <= n - 1,
ak[[1,
counterE +
1]] = (pimatk[[1 + 1, 1]]*ak[[1, counterE]] +
pimatk[[2 + 1, 1]]*ak[[2, counterE]])*
fikyj[[1, counterE + 1]];
ak[[2,
counterE +
1]] = (pimatk[[1 + 1, 2]]*ak[[1, counterE]] +
pimatk[[2 + 1, 2]]*ak[[2, counterE]])*
fikyj[[2, counterE + 1]];

ak[[1, counterE + 1]] = ak[[1, counterE + 1]];
ak[[2, counterE + 1]] = ak[[2, counterE + 1]];

bk[[1, n - counterE]] =
pimatk[[1 + 1, 1]]*fikyj[[1, n - counterE + 1]]*
bk[[1, n - counterE + 1]] +
pimatk[[1 + 1, 2]]*fikyj[[2, n - counterE + 1]]*
bk[[2, n - counterE + 1]];
bk[[2, n - counterE]] =
pimatk[[2 + 1, 1]]*fikyj[[1, n - counterE + 1]]*
bk[[1, n - counterE + 1]] +
pimatk[[2 + 1, 2]]*fikyj[[2, n - counterE + 1]]*
bk[[2, n - counterE + 1]];

bk[[1, n - counterE]] = bk[[1, n - counterE]];
bk[[2, n - counterE]] = bk[[2, n - counterE]];

counterE = counterE + 1;
];

τhijk =
Table[Table[
Table[(ak[[h, j]]*pimatk[[h + 1, i]]*fikyj[[i, j + 1]]*
bk[[i, j + 1]])/(Sum[
Sum[ak[[hh, j]]*pimatk[[hh + 1, ii]]*fikyj[[ii, j + 1]]*
bk[[ii, j + 1]], {ii, 1, 2}], {hh, 1, 2}]), {j, 1,
n - 1}], {i, 1, 2}], {h, 1, 2}];
τijk =
Table[Table[
ak[[i, j]]*
bk[[i, j]]/(Sum[
ak[[h, j]]*10^listak[[j]]*bk[[h, j]], {h, 1, 2}]), {j, 1,
n}], {i, 1, 2}];

(*M-step*)

teta0 = (Sum[
Sum[(Sum[
y[[j + 1]]*τijk[[i, j]], {j, 1, n}])*τijk[[i,
k]]*y[[k]]/(Sum[τijk[[i, j]], {j, 1, n}]) -
y[[k + 1]]*y[[k]]*τijk[[i, k]], {i, 1, 2}], {k, 1,
n}])/(Sum[
Sum[(Sum[
y[[j]]*τijk[[i, j]], {j, 1,
n}]/(Sum[τijk[[i, j]], {j, 1, n}]) -
y[[k]])*τijk[[i, k]]*y[[k]], {i, 1, 2}], {k, 1, n}]);

teta1 =
Sum[τijk[[1, j]]*(y[[j + 1]] - teta0*y[[j]]), {j, 1, n}]/
Sum[τijk[[1, j]], {j, 1, n}];

teta2 =
Sum[τijk[[2, j]]*(y[[j + 1]] - teta0*y[[j]]), {j, 1, n}]/
Sum[τijk[[2, j]], {j, 1, n}];

sigmaest =
Sqrt[Sum[((y[[j + 1]] - teta1 - teta0*y[[j]])^2*τijk[[1,
j]] + (y[[j + 1]] - teta2 - teta0*y[[j]])^2*τijk[[2,
j]]), {j, 1, n}]/(n)];

θk2 = {teta0, teta1, teta2};

(*Q[θ_,σ_]:=(τijk[[1,1]]*Log[pimatk[[1,
1]]]+τijk[[2,1]]*Log[pimatk[[1,2]]])+Sum[Sum[
Sum[τhijk[[h,i,j]]*Log[pimatk[[h+1,i]]],{j,1,n-1}],{i,1,
2}],{h,1,2}]+Sum[Sum[τijk[[i,j]]*Log[1/(Sqrt[2*
Pi]*σ)]*(-(y[[j+1]]-θ[[i+1]]-θ[[1]]*y[[
j]])^2/(2*σ^2)),{j,1,n}],{i,1,2}];*)

pimatk[[1, 1]] = τijk[[1, 1]];
pimatk[[1, 2]] = τijk[[2, 1]];(*π_ 0i=τijk[[i,
1]]*)

pimatk[[2, 1]] =
Sum[τhijk[[1, 1, j]], {j, 1, n - 1}]/
Sum[τijk[[1, j]], {j, 2, n}];
pimatk[[2, 2]] =
Sum[τhijk[[1, 2, j]], {j, 1, n - 1}]/
Sum[τijk[[1, j]], {j, 2, n}];
pimatk[[3, 1]] =
Sum[τhijk[[2, 1, j]], {j, 1, n - 1}]/
Sum[τijk[[2, j]], {j, 2, n}];
pimatk[[3, 2]] =
Sum[τhijk[[2, 2, j]], {j, 1, n - 1}]/
Sum[τijk[[2, j]], {j, 2, n}];(*π_hi=(\!$$\*UnderoverscriptBox[\(∑$$, $$j = 2$$, $$n$$]
\*SubscriptBox[$$τ$$, $$hij$$]\))/\!$$\*UnderoverscriptBox[\(∑$$, $$j = 2$$, $$n$$]
\*SubscriptBox[$$τ$$, $$hj$$]\)*)
(*Print["pimatk =",pimatk];*)

counterM = 1;
While[counterM <= n,
fikyj[[1, counterM]] = (1/(Sqrt[2*Pi]*sigmaest))*
E^(-(y[[counterM + 1]] - θk2[[2]] - θk2[[1]]*
y[[counterM]])^2/(2*sigmaest^2));

fikyj[[2, counterM]] = (1/(Sqrt[2*Pi]*sigmaest))*
E^(-(y[[counterM + 1]] - θk2[[3]] - θk2[[1]]*
y[[counterM]])^2/(2*sigmaest^2));
counterM = counterM + 1;

];

];

likehamtotal[[iterations]] = likeham[pimatk, fikyj];
pimattotal[[iterations]] = pimatk;
tetatotal[[iterations]] = θk2;
sigmatotal[[iterations]] = sigmaest;

iterations = iterations + 1;
];

iterations = 1;
pos = 1;
While[iterations <= ite ,
If[likehamtotal[[iterations]] > likehamtotal[[pos]],
pos = iterations];
iterations = iterations + 1
];
{pimattotal[[pos]], tetatotal[[pos]], sigmatotal[[pos]]}
];

• Can you post your code? We can't really help you otherwise! – dr.blochwave Oct 1 '14 at 9:14
• @blochwave Thanks for your interest. I've just posted my code. – An old man in the sea. Oct 1 '14 at 10:47
• That is a lot of code for anyone on this site to have time to read/critique/analyze it fully and most likely, no one is going to. However from a quick glance at the code, I can most certainly say that the slowness is due to the procedural programming style. All those Sum[Sum[Sum[Sum... terms in nested Whiles are bound to give terrible performance. Writing the algorithm in a functional manner would certainly speed things up a lot. Unfortunately, it is hard (for me at least) to recommend specific changes to make to your code. If I were you, I'd start afresh with functional (use Igor's link) – rm -rf Oct 1 '14 at 14:37
• I was about to comment on the nested Sum[Sum[Sum[...s! – dr.blochwave Oct 1 '14 at 14:40
• I would really like to see a Compiled version of this though. I think it could be pretty fast. I don't see anything that can't be compiled to C. – C. E. Oct 1 '14 at 15:14