\documentstyle{csart}

\begin{document}

\def\diag{\mbox {\bf diag}\,}
\newcommand{\boldmat}[1]{\mbox{\boldmath ${#1}$}}

\sloppy

\pagestyle{empty}


\title{Comments on Parallelism in Hidden Markov Model Algorithms}

\author{Bryan Carpenter}
%\address{Department of Electronics and Computer Science, \\
%         Southampton University, \\
%         Southampton, S09 5NH} 

\maketitle

\pagestyle{plain}

\section*{Hidden Markov Models}

     A Hidden Markov Model (HMM) is characterised by a set of $N$ states and
a set of $M$ possible observations on these states.
It is also parametrised by a vector and
two matrices: $\pi_i$ is the probability of the system starting in state $i$;
$a_{ij}$ is the single-step transition probability from state $i$ to $j$;
and $b_i(O)$ is the probability of making observation $O$ when the
system is in state $i$.  Reference \cite{hmm} is a simple
introduction to some of the basic computations involving HMMs.
The three problems discussed there are
\begin{enumerate}
\item{``Evaluation''.} For a given HMM, what is the probability of
of making the set of observations $O_1, \ldots , O_T$ in time-steps
$1, \ldots, T$?
\item{``Decoding''.} For a given HMM and a given sequence of observations
$O_1, \ldots , O_T$, what is the best guess at the underlying sequence of
states that yielded these observations?
\item{``Estimation''.} Given a training sequence of $T$ observations
(presumably $T$ is quite large), how do you choose $\pi$, $a$ and $b$
to maximise the probability of that observation sequence?
\end{enumerate}
Problem 1 can be solved by the ``forward-backward'' algorithm, problem
2 by the ``Viterbi'' algorithm, and 3 by the ``Baum-Welch'' re-estimation
formulae (which improve the target probability: they can be applied
iteratively to approach an optimum).  These algorithms are all presented
in \cite{hmm}, and I'm not going to reproduce the details here.

     An important ingredient is the forward-backward scheme.
If you consider problem 1, the probability of the system passing
through the sequence of states $i_1, \ldots , i_T$ in the first
$T$ steps is
\begin{equation}
\pi_{i_1} a_{i_1i_2} a_{i_2i_3} \ldots a_{i_{T - 1}i_T}
\end{equation}
{If} the system passed through this sequence of states the
probability of the observation sequence would be
\begin{equation}
b_{i_1}(O_1) b_{i_2}(O_2) \ldots b_{i_T}(O_T)
\end{equation}
Multiplying and summing over {\em all} possible state sequences gives the
probability
\begin{eqnarray*}
&& \sum_{i_1, i_2, \ldots, i_T}
\pi_{i_1} b_{i_1}(O_1) a_{i_1i_2} b_{i_2}(O_2) a_{i_2i_3} \ldots
a_{i_{T - 1}i_T} b_{i_T}(O_T) \\
&=& \boldmat{\pi}^\dagger \cdot \diag (\boldmat{b} (O_1)) \cdot \boldmat{a} \ldots
\diag (\boldmat{b} (O_{T-1})) \cdot \boldmat{a} \cdot \boldmat{b} (O_T)
\end{eqnarray*}
Two obvious orders of evaluation for the matrix product
are to start at the left and go
forwards, or on the right and go backwards (this sounds like a trivial
observation, but the intermediate results from the two schemes
are both needed in the Baum-Welch re-estimation formulae).

\section*{Remarks on Complexity}

     I'll briefly discuss ``fast'' parallel versions of the algorithms.
Mostly, these schemes aren't particularly practical, but they're a starting
point.

     The evaluation problem, and calculation of the forward and
backward intermediate products, are ``in NC''---that is, they can be computed
in (poly-)logarithmic time with a polynomially bounded number of processors.
Using the standard parallel prefix technique \cite{gibryt},
applied here to multiplication of $N$ by $N$ matrices, the partial results for
\begin{equation}
\prod_{t = 1}^{T - 1} \diag (\boldmat{b} (O_t)) \cdot \boldmat{a}
\end{equation}
can all be computed on $\sigma_m T / \log T$ processors in $\tau_m \log T$
time, where $\tau_m$, $\sigma_m$ are the time and processor requirements
for a single matrix multiplication.  Using the obvious tree-based scheme
for the sums in the matrix multiplication we can easily get parallel time
$\log N \log T$ on $N^3 T / (\log N \log T)$ processors.  The problem here is
that the sequential algorithm only needs to work with vector intermediate
results and matrix-vector products, so its operation count is
only $N^2 T$.  This means the parallel algorithm just described isn't
optimal, in the sense of using processor cycles efficiently.
I don't know how to solve this problem and maintain logarithmic
time.  Is there such an algorithm?

     The Baum-Welch re-estimation formulae just involve taking products
and sums of the intermediate vectors calculated from the
forward-backward schemes.  Given the output of the forward-backward
algorithm, calculating all the $\gamma$ and $\xi$ coefficients
described in \cite{hmm} takes $N^2 T$ operations.  Calculating the new
$\pi$, $a$, $b$ coefficients from these takes the same again.  There
aren't any complicated dependencies other than sums, and it should be
straightforward to do this part of the calculation in $\log T$ time on
an optimal number of processors.  Unfortunately, there's no more work
here than in the forward-backward phase, so optimal parallelisation of
this stage leaves the first one as a bottleneck.

     Parallelising the Viterbi algorithm looks a bit more interesting,
and I don't know exactly how to do it.  It's a sort of dynamic
programming problem.  Fast algorithms for a different class of dynamic
programming problems are given in \cite{gibryt}.  Also, I think the
Viterbi algorithm probably comes under the class of straight-line problems
for which fast algorithms are given in the same place,
but I haven't checked that in detail.
I'd guess that there probably is a logarithmic-time
parallelisation, but if it's similar to the algorithms just mentioned
it probably won't be optimal in processor usage.

\section*{The Real World}

     As pointed out above, a maximally parallel version of the
important forward-backward algorithm seems to involve
more operations by a factor of $N$ than the sequential version.
This {\em may} be acceptable if $N$ is very small.  For
example, in \cite{hmm}, Rabiner and Juang describe an isolated word
recognition example with 5 states.

     If, on the other hand, $N$ is large (some of Sharman's
examples in \cite{tag} involve hundreds or even thousands of
states---I'm not sure how practical these examples are meant
to be) it may be sufficient to extract parallelism purely from
the individual matrix-vector products.  For vectors of size
a hundred or so (perhaps less) it should be possible to
get healthy speedups on existing parallel machines.
In this regime of $N$ it should also be
possible to accelerate the Viterbi algorithm by parallelising
the individual maximisation stages.

     Applications described by Rabiner and Juang involve
building separate HMMs for each known word.
Each word to be recognised is checked against the models
in turn.  This usage is obviously amenable to straightforward,
highly efficient parallelisation, by farming.

\begin{thebibliography}{9}

\bibitem{hmm} Rabiner, L.R.~and B.H.~Juang, {\em An Introduction
to Hidden Markov Models}, IEEE ASSP Mag. vol.~3, no.~1 Jan 1986.

\bibitem{gibryt} Gibbons, A.~ and W.~Rytter, {\em Efficient
Parallel Algorithms}, Cambridge University Press, 1988.

\bibitem{tag} Sharman, R.A., {\em Hidden Markov Model Methods
for Word Tagging}, IBM UKSC Report 214.

\end{thebibliography}

\end{document}



