BMC Bioinformatics BioMed Central Methodology article Implementing EM and Viterbi algorithms for Hidden Markov Model in linear memory

Background The Baum-Welch learning procedure for Hidden Markov Models (HMMs) provides a powerful tool for tailoring HMM topologies to data for use in knowledge discovery and clustering. A linear memory procedure recently proposed by Miklós, I. and Meyer, I.M. describes a memory sparse version of the Baum-Welch algorithm with modifications to the original probabilistic table topologies to make memory use independent of sequence length (and linearly dependent on state number). The original description of the technique has some errors that we amend. We then compare the corrected implementation on a variety of data sets with conventional and checkpointing implementations. Results We provide a correct recurrence relation for the emission parameter estimate and extend it to parameter estimates of the Normal distribution. To accelerate estimation of the prior state probabilities, and decrease memory use, we reverse the originally proposed forward sweep. We describe different scaling strategies necessary in all real implementations of the algorithm to prevent underflow. In this paper we also describe our approach to a linear memory implementation of the Viterbi decoding algorithm (with linearity in the sequence length, while memory use is approximately independent of state number). We demonstrate the use of the linear memory implementation on an extended Duration Hidden Markov Model (DHMM) and on an HMM with a spike detection topology. Comparing the various implementations of the Baum-Welch procedure we find that the checkpointing algorithm produces the best overall tradeoff between memory use and speed. In cases where sequence length is very large (for Baum-Welch), or state number is very large (for Viterbi), the linear memory methods outlined may offer some utility. Conclusion Our performance-optimized Java implementations of Baum-Welch algorithm are available at . The described method and implementations will aid sequence alignment, gene structure prediction, HMM profile training, nanopore ionic flow blockades analysis and many other domains that require efficient HMM training with EM.


Background
Hidden Markov Models (HMMs) are a widely accepted modeling tool [1] used in various domains, such as speech recognition [2] and bioinformatics [3]. An HMM can be described as a stochastic finite state machine where each transition between hidden states ends with a symbol emission. The HMM can be represented as a directed graph with N states where each state can emit either a discrete character or a continuous value drawn from a Probability Density Function (PDF).
We are interested in a distributed HMM analysis of the channel current blockade signal caused by a single DNA hairpin molecule held in a nanopore detector [4,5]. The molecules examined frequently produce toggles with stationary statistical profiles for thousands of milliseconds.
With a sampling rate of 20 μs, processing even a modest blockade signal of 200 ms duration (10,000 sample points) becomes problematic, mostly because of the size of the dynamic programming tables required in the conventional implementations of the HMM's Baum-Welch and Viterbi decoding algorithms. Since we are also trying to model durations [6] and spike phenomena [7], by increasing the number of HMM states, conventional HMM implementations are found to be prohibitively expensive in terms of memory use.
The Baum-Welch algorithm is an Expectation Maximization (EM) algorithm invented by Leonard E. Baum and Lloyd R. Welch, and first appears in [8]. A later refinement, Hirschberg's algorithm for an HMM [9], reduces the memory footprint by recursively halving the pairwise alignment dynamic programming table for sequences of comparable size. In our application domain, the length of the observed emission sequence (in the case of nanopore ionic flow blockade analysis or gene structure prediction) is prohibitively long compared to the number of HMM states. Further, Baum-Welch requires multiple paths, instead of the most likely one, making this strategy less than optimal.
The checkpointing algorithm [10][11][12] implements the Baum-Welch algorithm in memory and in O(TNQ max + T(Q + E)) processor time, where T is the length of the observed sequence, Q max is the maximum HMM node out-degree, E is the number of free emission parameters and Q is the number of free transition parameters. It divides the input sequence into blocks of symbols each, and, during the forward pass, retains the first column of the forward probability table for each block. When the reverse sweep starts, the forward values for each block are sequentially re-evaluated, beginning with their corresponding checkpoints, to update the parameter estimates.
Further refinement to the algorithm, as described in [13] and amended here, has rendered the memory demands independent of the observed sequence length, with O(N(Q + ED)) memory usage and O(TNQ max (Q + ED)) running time, where D is the dimensionality of a vector that stores statistics on the emission PDF parameter estimates. Performance of the various algorithms is summarized in Table 1. In this work, we also present a modification of one of the key HMM algorithms, the Viterbi algorithm, improving the memory profile without affecting the execution time.

HMM definition, EM learning and Viterbi decoding
The following parameters describe the conventional HMM implementation according to [14]: • A set of states S = {S 1 ,..., S N } with q t being the state visited at time t,  [15]. The Viterbi algorithm, as shown in Table 2, is a dynamic programming algorithm that runs an HMM to find the most likely sequence of hidden states, called the Viterbi path, that result in an observed sequence. When training the HMM using the Baum-Welch algorithm (an Expectation Maximization procedure), first we need to find the expected probabilities of being at a certain state at a certain time-point using the forward-backward procedure as shown in Table 2. The forward, backward, and Viterbi algorithms take O(TNQ max ) time to execute.
Let us define ξ t (i, j) as the probability of being in state i at time t, and state j at time t + 1, given the model and the observation sequence and γ t (i) as the probability of being in state i at time t, given the observation sequence and the model The HMM maximization step using these probabilities is shown in Table 3. The conventional EM procedure for HMM [14] takes O(TN) memory and O(TNQ max + T (Q + E)) processor time. An HMM containing empty internal states (see for example [3]) and Hierarchical HMM (HHMM) could be converted into canonical HMM form through stack transformation as discussed in [16]. Figure 1 outlines initial, simple transition probability calculations for all possible paths through a "toy" HMM. In Figure 1, to estimate the probability of transition from state 1 to state 2 (1 → 2), we calculate the probability of transition utilization at time intervals 1-2 and 2-3 as:

Forward sweep strategy explained
In particular, to estimate the probability of transition 1 → 2 at time interval 1-2 we calculate the sum of probabilities of all possible paths that lead to state 1 at time-point 1 (forward probability α 1 (1)). Then we multiply this probability of being in state 1 by the transition a 1,2 and emission b 2 (o 2 ) probabilities.
Further multiplication by the sum of probabilities of all possible paths from state 2 at time 2 until the end of the emission sequence (backward probability β 2 (2)), is the expected probability of transition use. The sum of these estimates at time-points 1-2 and 2-3 is equivalent to the transition probability estimate in Table 3 (prior to normalization).
• Finally is the sequence likelihood . Following equation (1), at a certain time-point t we need to score the evidence supporting transition between nodes i and j, which is the sum of probabilities of all possible state paths that emit the transition probability estimate shown in Table 3. Thus, we estimate transition utilization using: where out(S i ) is the set of nodes connected by edges from S i .
According to [13]  The following algorithm updates parameters for the set of discrete symbol probability distributions.
After initialization we make the recurrence steps, where we correct the emission recurrence presented in [13]

Initial probability estimate Transition probability estimate
Emission parameters estimate , for 1 ≤ j ≤ N, The linear memory implementation of Baum-Welch learning algorithm for HMM 1 Initialization Recurrence the trellis through available transitions a n, m . Thus the definition of e i (γ, t, m) holds true for the time-point t.
At the end of recurrence we marginalize the final state m out of probability e i (γ, T, m) to get the weighted sum of all state paths making observation γ in state i at various timepoints equivalent to the numerator of the discrete emission parameter estimate in Table 3, which is a weighted sum of all possible paths that score emissions evidence at certain time-points. By normalizing these scores we estimate the emission parameters.
The forward sweep strategy was originally formulated in [13] for

Linear memory Baum-Welch using a backward sweep with scaling
The objective of the algorithm presented in this section is equivalent to that discussed previously [see Section Forward sweep strategy explained] based on forward probabilities of state occupation. However, by using the backward probabilities of state occupation we are able to estimate initial HMM state probabilities much more quickly. In the description that follows we introduce a new set of probabilities:  All calculations are based on backward probability β t (i), but there is inevitably insufficient precision to directly represent these values for significantly long emission sequences. Therefore we scale the backward probability during the recursion.
The linear memory Baum-Welch implementation is shown in Figure 2, where is a set of nodes with free emission parameters and is a set of nodes with free emanating transitions. Scaling relationships used at every iteration are rigorously proven [see Appendix A]. An alternative to scaling is relation (7) presented in [17] showing how to efficiently add log probabilities The scoring functions used for the emissions updates are shown in Table 4. With discrete emission we sum all the backward probabilities of state occupation given discrete symbol emission at certain time-points and later we normalize these counts in (8). In the case of a normally distributed continuous PDF we sum emissions and emission deviation from state i mean squared. These sums are scaled by probability of state occupation. We use these counts to estimate the emission mean (9) and variance (10) for a given state.

Parameters update
We estimate the initial probability according to equations presented in Table 3

Discrete emission Continuous Gaussian emission
The transition estimate is calculated as following

Viterbi decoding in linear memory
In this section we describe results when using a "linear memory" modification of the original Viterbi algorithm that was introduced in [18] by Andrew J. Viterbi. As mentioned previously, the Viterbi algorithm is a dynamic programming algorithm for finding the most likely sequence of hidden states, called the "Viterbi path", corresponding to the sequence of observed events in the context of an HMM. Viterbi checkpointing implementation introduced in [11] divides input sequence into blocks of symbols each and during the first Viterbi pass keeps only the first column of the d table for each block. The reconstruction of the most probable state path starts with the last block, where we use the last checkpointing column to initialize recovery of the last states of the most likely state path and so on, until the entire sequence decoding is reconstructed. The algorithm requires memory space pro-portional to and takes computational time O(TNQ max ) twice as much as conventional implementation would take because of multiple sweeps. Additional computations could be easily justified by the lower memory use, which in practice leads to substantial speedup.
Only two columns are needed for the δ array that stores maximum probability scores for a state at a given timepoint for the algorithm to run (referring to the relationship shown in Table 2). We replace the array, needed to store the dynamic programming backtrack pointers, by a linked list. Our approximately linear memory implementation follows the observation that the backtrack paths typically converge to the most likely state path and travel all together to the beginning of the decoding table.
Although the approximate linearity depends on model structure and emission sequence decoded, and is not guaranteed, this behavior is typical for the HMM topologies we use. The possibility of using O(N log(T)) space (in case we write to disk the most likely path before the coalescence point, i.e. the first state on the backtrack path where only a single candidate is left for the initial segment of the most probable state path) has only been rigorously proven for simple symmetric two-state HMM [19].
The modified Viterbi algorithm is shown in Figure 3. This approach has proven to be useful in decoding the explicit Duration Hidden Markov Model (DHMM) topology introduced in [6], as can be seen in Figure 4. Although this implementation closely follows the originally proposed non-parametric duration density [20] design, thê , ( Figure 3 Viterbi algorithm implementation with linked list. Here is reference to the previous node.

Viterbi algorithm implementation with linked list
with optimal decoding Q * = {q * 1 , q * 2 , . . . , q * T }. advantage is that we no longer have to use highly modified Forward-Backward and Viterbi algorithms [21]. This topology directly corresponds to the Generalized Hidden Markov Model (GHMM) used in GENSCAN [22], one of the most accurate gene structure prediction tools. The potential for a very large number of nodes in our proposed topology is compensated by introducing the linear memory Viterbi and Baum-Welch implementations with unaltered running time O(STM), where M is the maximum duration of an aggregate state and S is the number of aggregate states. An example of backtracking path compression for such a topology with discrete emissions can be seen in Figure 5, where the most likely trail of states quickly combines with all the alternative trails. Another interesting topology used by our laboratory for "spike" detection is shown in Figure 6, where the spikes are modelled as a mixture of two trajectories interconnected with an underlying set of ionic flow blockade states. The Viterbi decoding trail for this topology, detecting two sequential spikes in samples from the real continuous ionic flow blockade, is shown in Figure 7. Again, the backtracks quickly converge to the most likely state sequence.
Our particular implementation relies on the Java Garbage Collector (GC), which periodically deletes all the linked list nodes allocated in the heap that are no longer referenced by the active program stack as shaded in light gray color in Figures 5 and 7. On multiple core machines the GC could run concurrently with the main computational thread thus not obstructing execution of the method. In the lower level languages (like C/C++) "smart pointers" could be used to deallocate unused links when their reference count drops to zero, which is in some ways even more efficient than Java's garbage collection.

Computational performance
We conducted experiments on the HMM topologies mentioned above [see Section Viterbi decoding in linear memory] with both artificial and real data sets, and evaluated performance of the various implementations of the Viterbi and EM learning algorithms. We describe the performance of the Java Virtual Machive (JVM) after the HotSpot Java byte code optimizer burns-in, i.e. after it optimizes both memory use and execution time within EM cycles. The linear memory, checkpointing and conventional algorithm implementations are thereby streamlined to avoid an unfair comparison due to obvious performance bottlenecks.
For the DHMM topology shown in Figure 4 we have chosen to systematically alter the size of two aggregate states from 50 to 500 when learning on an artificially generated sequence of 10,000 discrete symbols to see how the number of free learning parameters affects the performance of the EM learning algorithms. In Subfigures 8(a) -8(c) we see that the running time of the linear implementation grows dramatically compared to conventional and checkpointing implementations, making it a very slow alternative for such a scenario. Although the linear implementation memory profile is low, as expected, for high Explicit Duration HMM trellis for the observation string shown below Figure 5 Explicit Duration HMM trellis for the observation string shown below. The most likely sequence of states for the observation string shown below is shaded. The lightly grayed states will be deallocated by garbage collector. In experiments on EM learning on a spike detection HMM topology, shown in Figure 6, we have systematically varied the ionic flow duration from 1,000 ms to 64,000 ms. Although in Subfigures 8(d) -8(f) duration of the time cycle of the linear memory implementation is not so inflated in this situation, it is still many times higher than for conventional and checkpointing algorithms. Please note that conventional and checkpointing algorithms spend almost identical time per cycle. The conventional implementation still takes the largest amount of memory and once again checkpointing takes less memory com-pared to the linear memory implementation in the case of long signal duration. Garbage collection in the case of the conventional implementation starts taking a substantial fraction of the CPU time for maximum signal duration, which advocates against using the conventional implementation for the analysis of long signals.
Theoretically, for the linear memory algorithm to run faster than checkpointing alternative the following condition should hold true 2TNQ max + T (Q + E) > TNQ max (Q + ED) which reduces to the condition, unrealistic for any practical model 2 > Q + ED. Thus, the linear memory algorithm will always run slower than checkpointing. The memory condition for the linear memory EM algorithm implementation to run in less space is 2N > N(Q + ED), which reduces to 2 > Q + ED condition -quite realistic for sufficiently In both test scenarios shown in Figure 8 we see that conventional implementation of Baum-Welch aggressively claims very large amounts of heap space, even for modestly sized problems (in some applications, such as the JAHMM package [23], it allocates the probabilistic In Tables 5 and 6 we list the ratio of memory used by the linked list nodes referenced from the active program stack to the sequence length T. As could be seen, it quickly becomes proportional to 1.0 in both spike detection and the explicit DHMM topologies as the decoded sequence length grows.

Discussion and Conclusion
We have discussed implementation of Baum-Welch and Viterbi algorithms in linear memory. We successfully used Trellis for loopy topology used for spikes detection where shallow spike (states 1-6) and deep spike (states 7-17) are conse-quently decoded Figure 7 Trellis for loopy topology used for spikes detection where shallow spike (states 1-6) and deep spike (states 7-17) are consequently decoded. The most likely sequence of states for the sequence of observed ionic flow current blockades (in pA) shown below is shaded. The lightly grayed states will be deallocated by garbage collector.   In subfigures 8(a) -8(c) performance of Baum-Welch used on DHMM topology with two aggregate states of various maximum duration D Figure 8 In these implementations in analysis of nanopore blockade signals with very large sample sizes (more than 3,000 ms) where the main limitation becomes processing time rather than memory overflow. We are currently working on efficient distributed implementations of the linear memory algorithms to facilitate quicker, potentially "realtime" applications.
In both of the test scenarios considered, the linear memory modification of the EM algorithm takes significantly more time to execute compared to conventional and checkpointing implementations. Despite being the fastest in many realistic scenarios, the conventional implementation of the EM learning algorithm suffers from substantial memory use. The checkpointing algorithm alleviates both of these extremes, sometimes running even faster than the conventional algorithm due to lower memory management overhead. The checkpointing algorithm seems to provide an excellent tradeoff between memory use and speed. We are trying to understand if the running time of our linear memory EM algorithm implementation can be constrained in a way similar to the checkpointing algorithm. A demo program featuring the canonical, checkpointing and linear memory implementations of the EM learning and the Viterbi decoding algorithms is available on our web site (see Availability & requirements section below).