Implementing EM and Viterbi algorithms for Hidden Markov Model in linear memory

Abstract

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 http://logos.cs.uno.edu/~achurban. 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–12] implements the Baum-Welch algorithm in $O(\sqrt{T}N)$ 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 $\sqrt{T}$ blocks of $\sqrt{T}$ 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.

Table 1 The computational expense of different algorithm implementations running on HMM.

Methods and Results

HMM definition, EM learning and Viterbi decoding

The following parameters describe the conventional HMM implementation according to [14]:

• A set of states S = {S₁,..., S _N } with q _t being the state visited at time t,

• A set of PDFs B = {b₁(o),..., b _N (o)}, describing the emission probabilities b _j (o _t ) = p(o _t |q _t = S _j ) for 1 ≤ j ≤ N, where o _t is the observation at time-point t from the sequence of observations O = {o₁,..., o _T },

• The state-transition probability matrix A = {a_{i, j}} for 1 ≤ i, j ≤ N, where a_{i, j}= p(q_{t + 1}= S _j |q _t = S _i ),

• The initial state distribution vector Π = {π₁,..., π _N }.

A set of parameters λ = (Π, A, B) completely specifies an HMM. Here we describe the HMM parameter update rules for the EM learning algorithm rigorously derived in [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.

Table 2 The Viterbi decoding, forward and backward procedures.

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

$${\xi }_t(i,j)=p(q_t=S_i,q_{t+1}=S_j|O,\lambda )=\frac{{\alpha }_t(i)a_{i,j}{b}_j(o_{t+1}){\beta }_{t+1}(j)}{p(O|\lambda )},$$

(1)

and γ _t (i) as the probability of being in state i at time t, given the observation sequence and the model

$${\gamma }_t(i)=p(q_t=S_i|O,\lambda )=\frac{{\alpha }_t(i){\beta }_t(i)}{{\displaystyle {\sum }_{i=1}^N{\alpha }_t(i){\beta }_t(i)}}={\displaystyle \sum _{j=1}^N{\xi }_t(i,j)}.$$

(2)

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

Table 3 The maximization step in HMM learning. states.

Forward sweep strategy explained

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:

Figure 1

Time trellis for simple model where possible emissions of 0 and 1 are shown above and below trellis. Probabilities of emissions that happen after each transition are shown in bold and transitions of interest taken at certain time-point are underlined.

p(Making transition 1 → 2 at time 1–2) = α₁(1) × a_1,2 × b₂(o₂) × β₂(2), p(Making transition 1 → 2 at time 2–3) = α₂(1) × a_1,2 × b₂(o₃) × β₃(2).

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)). Then we multiply this probability of being in state 1 by the transition a_1,2 and emission b₂(o₂) 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)), 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).

According to [13] t_{i, j}(t, m) is the weighted sum of probabilities of all possible state paths that emit subsequence o₁,..., o _t and finish in state m, taking an i → j transition at least once where the weight of each state path is the number of i → j transitions taken. Processing of the entire t_{i, j}(t, m) recurrence takes memory proportional to O(NQ) and processor time O(TNQQ _max ). Initially we have t_{i, j}(1, m) = 0 since no transitions have been made. After initialization, we perform the following recurrence operations:

$$t_{i,j}(t,m)={\alpha }_{t-1}(i)a_{i,m}{b}_m(o_t)\delta (m=j)$$

(3)

$$+{\displaystyle \sum _{n=1}^N{t}_{i,j}(t-1,n)a_{n,m}{b}_m(o_t)},$$

(4)

where $\delta (m=j)=\{\begin{array}{ll}1,\hfill & \begin{array}{cc}\text{if}& m=j\end{array}\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}$. 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 subsequence o₁,..., o_t-1and finish in state i (forward probability α_t-1(i)), multiplied by transition a_{i, j}and emission b _j (o _t ) probabilities upon arrival to o _t . We extend weighted paths containing evidence of i → j transitions made at previous time-points 1,..., t - 1 further down the trellis in subequation (4). Finally, by the end of the recurrence we marginalize the final state m out of probability t_{i, j}(T, m) to get a weighted sum of state paths taking transition i → j at various time-points that is equivalent to the numerator in the transition probability estimate shown in Table 3. Thus, we estimate transition utilization using:

$$\begin{array}{cc}{t}_{i,j}^{END}={\displaystyle \sum _{m=1}^N{t}_{i,j}(T,m),}& a_{i,j}=\frac{t_{i,j}^{END}}{{\displaystyle {\sum }_{j\in out(S_i)}{t}_{i,j}^{END}}},\end{array}$$

where out(S _i ) is the set of nodes connected by edges from S _i .

According to [13] e _i (γ, t, m) is the weighted sum of probabilities of all possible state paths that emit subsequence o₁,..., o _t and finish in state m, for which state i emits observation γ at least once where the weight of each state path is the number of γ emissions that it makes from state i.

The following algorithm updates parameters for the set of discrete symbol probability distributions.

Initialization step e _i (γ, 1, m) = α₁(m) δ (i = m) δ (γ = o₁). After initialization we make the recurrence steps, where we correct the emission recurrence presented in [13] [see Additional File 1]:

$$e_i(\gamma ,t,m)={\alpha }_t(m)\delta (i=m)\delta (\gamma =o_t)$$

(5)

$$={\displaystyle \sum _{n=1}^N{e}_i(\gamma ,t-1,n)a_{n,m}{b}_m(o_t)}.$$

(6)

Finally, by the end of the recurrence we marginalize the final state m out of e _i (γ, T, m) and estimate the emission parameters through normalization

$$\begin{array}{cc}{e}_i^{END}(\gamma )={\displaystyle \sum _{m=1}^N{e}_i(\gamma ,T,m),}& {\widehat{b}}_j(\gamma )=\frac{e_j^{END}(\gamma )}{{\displaystyle {\sum }_{\gamma =1}^D{e}_j^{END}(\gamma )}}.\end{array}$$

The algorithm for discrete emission parameters estimate B = {b₁(o),..., b _N (o)} takes in O(NED) memory and O(TNEDQ _max ) time. As discussed [see Subsection HMM definition, EM learning and Viterbi decoding] the forward sweep takes O(TNQ _max ) time, where only the values of α_t-1(i) for 1 ≤ i ≤ N are needed to evaluate αt(i), thus reducing the memory requirement to O(N) for the forward algorithm. Computing e _i (γ, t, m) takes O(NED) previous probabilities of e _i (γ, t - 1, m) for 1 ≤ m ≤ N, 1 ≤ i ≤ E, 1≤ γ ≤ D. Recurrent updating of each e _i (γ, t, m) probability elements takes O(Q _max ) summations, totalling O(TNEDQ _max ).

Theorem 1 e _i (γ, t, m) is the weighted sum of probabilities of all possible state paths that emit subsequence o₁,..., o _t and finish in state m, for which state i emits observation γ at least once.

Proof

Initialization The only chance for a path at a time-point 1 to emit symbol γ at least once from state i and finish in state m is α₁(m) in case i = m and γ = o₁. Therefore, initialization conditions satisfy definition of e _i (γ, t, m).

Induction Suppose e _i (γ, t - 1, m) represents correct weighted sum of probabilities of all possible state paths that emit subsequence o₁,..., o_t-1and finish in state m, for which state i emits observation γ at least once. We need to prove the above holds for time-point t. Following equation (1) in recurrence part (5) we score the evidence of symbol o _t emission from state i at time-point t, which will be further propagated down the trellis in recurrence part (6). Chances of such event is α _t (m), i.e. sum of probabilities of all possible state paths finishing in state m at time-point t under conditions i = m and γ = o _t . The second part of the recurrence (6) extends the weighted paths containing evidence of γ symbol emission from state i at previous time-points 1,..., t - 1 and finishing in state n further down 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 time-points 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 HMMs with silent Start/End states, and automatically handles the prior probabilities estimates for the states as standard transitions connecting Start with other non-silent states. The prior transition estimates a_{Start, i}are naturally involved within recurrent updates of t_{i, j}(t, m), which takes an additional O(N²) memory if all N non-silent states have non-zero priors with time cost O(TN²Q _max ). In order to compute the prior estimates in the conventional HMM formulation we need to know the backward probability at time-point 1, which requires calculation of the entire backward table. Therefore, in the next section we present a linear memory Baum-Welch algorithm modification built around a backward sweep with scaling, which only involves calculation of α₁(i) for 1 ≤ i ≤ N to estimate priors in O(N) time and O(N) memory.

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:

E _i (γ, t, m) – the weighted sum of probabilities of all possible state paths that emit subsequence o _t ,..., o _T and finish in state m, for which state i emits observation γ at least once, where the weight of each state path is the number of γ emissions that it makes from state i.

T_{i, j}(t, m) – the weighted sum of probabilities of all possible state paths that emit subsequence o _t ,..., o _T and finish in state m, taking i → j transition at least once, where the weight of each state path is the number of i → j transitions that it takes.

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 $\mathcal{E}$ is a set of nodes with free emission parameters and $\mathcal{T}$ 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

Figure 2

The linear memory implementation of Baum-Welch learning algorithm for HMM. This algorithm takes set of HMM parameters λ and sequence of symbols O. Expected HMM parameters are calculated according to formulas [see Subsection Parameters update].

$$\log \left({\displaystyle \sum _{t=0}^{N-1}{p}_i}\right)=\log p_0+\log \left(1+{\displaystyle \sum _{i=1}^{N-1}{e}^{\log p_i-\log p_0}}\right).$$

(7)

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.

Table 4 The scoring functions for discrete and continuous emissions.

Parameters update

We estimate the initial probability according to equations presented in Table 3, where D₁ is defined in Appendix A

$${\widehat{\pi }}_i=\frac{{\alpha }_1(i){\tilde{\beta }}_i(1)}{{\displaystyle {\sum }_{i=1}^N{\alpha }_1(i)}{\tilde{\beta }}_i(1)}=\frac{{\alpha }_1(i)D_1{\beta }_1(i)}{{\displaystyle {\sum }_{i=1}^N{\alpha }_1(i)D_1}{\beta }_1(i)}=\frac{{\alpha }_1(i){\beta }_1(i)}{{\displaystyle {\sum }_{i=1}^N{\alpha }_1(i)}{\beta }_1(i)}.$$

The emissions estimate for the discrete case are

$${\widehat{b}}_j(\gamma )=\frac{{\tilde{E}}_j^{END}(\gamma )}{{\displaystyle {\sum }_{\gamma =1}^D{\tilde{E}}_j^{END}(\gamma )}}=\frac{D_1{E}_j^{END}(\gamma )}{D_1{\displaystyle {\sum }_{\gamma =1}^D{E}_j^{END}(\gamma )}}=\frac{E_j^{END}(\gamma )}{{\displaystyle {\sum }_{\gamma =1}^D{E}_j^{END}(\gamma )}}.$$

(8)

For normally distributed continuous observation PDF

$${\widehat{b}}_j(o)\to \mu =\frac{{\tilde{E}}_j^{END}(1)}{{\tilde{E}}_j^{END}(3)}=\frac{D_1{E}_j^{END}(1)}{D_1{E}_j^{END}(3)}=\frac{E_j^{END}(1)}{E_j^{END}(3)},$$

(9)

$${\widehat{b}}_j(o)\to {\sigma }^2=\frac{{\tilde{E}}_j^{END}(2)}{{\tilde{E}}_j^{END}(3)}=\frac{D_1{E}_j^{END}(2)}{D_1{E}_j^{END}(3)}=\frac{E_j^{END}(2)}{E_j^{END}(3)}.$$

(10)

The transition estimate is calculated as following

$$a_{i,j}=\frac{{\tilde{T}}_{i,j}^{END}}{{\displaystyle {\sum }_{j\in out(S_i)}{\tilde{T}}_{i,j}^{END}}}=\frac{D_1{T}_{i,j}^{END}}{D_1{\displaystyle {\sum }_{j\in out(S_i)}{T}_{i,j}^{END}}}=\frac{T_{i,j}^{END}}{{\displaystyle {\sum }_{j\in out(S_i)}{T}_{i,j}^{END}}}\text{ for }i\in \mathcal{T}.$$

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 $\sqrt{T}$ blocks of $\sqrt{T}$ 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 $\sqrt{T}$ states of the most likely state path and so on, until the entire sequence decoding is reconstructed. The algorithm requires memory space proportional to $O(N\sqrt{T}+T)$ 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 time-point 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. It runs in the same O(TNQ _max ) time as a conventional Viterbi decoder, but takes the amount of memory O(T) as has been demonstrated by our simulations [see Section Computational performance].

Figure 3

Viterbi algorithm implementation with linked list. Here ${\psi }_{t+1}(q_{t+1}^{\ast })\to prev$ is reference to the previous node.

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, the 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.

Figure 4

Explicit DHMM topology. Here the first aggregate state emits 0 with probability 0.75 and 1 with probability 0.25 and the second aggregate state emits 0 with probability 0.25 and 1 with probability 0.75. Duration histograms are shown for each aggregate state.

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.

Figure 6

Spike detection loop topology.

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.

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 values of D, the checkpointing algorithm claims the least memory. This is because the table sizes in the linear memory EM implementation grow quickly as the number of states and transitions in the model increases. Garbage collection for large D is the lowest for the checkpointing EM compared to other implementations.

Figure 8

In subfigures 8(a) – 8(c) performance of Baum-Welch used on DHMM topology with two aggregate states of various maximum duration D . In subfigures 8(d) – 8(f) performance of Baum-Welch algorithm used on spike topology for various ionic flow durations is shown.

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 compared 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$\sqrt{T}$ > N(Q + ED), which reduces to 2$\sqrt{T}$ > Q + ED condition – quite realistic for sufficiently large values of T. The memory advantage is the key attribute since efforts are underway (Z. Jiang and S. Winters-Hilt) to perform segmented linear HMM processing on a distributed platform, where the speed deficiency is offset by M-fold speedup when using M nodes.

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 table ξ of size N² × T, although we do it in N × T through progressive summation of forward-backward tables), where both modified EM algorithm implementations have a very compact memory signature. An HMM algorithm implemented based on forward sweep strategy with silent Start/End states runs slower and takes more memory compared to the proposed backward sweep strategy in case of DHMM topology. This is because prior probabilities of the states are estimated as regular transitions from the Start state, thus substantially increasing t_{i, j}(t, m) table size and time required for a recurrent step.

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.

Table 5 Memory use for Viterbi decoding on spike topology with loop sizes 6 and 11.

Table 6 Memory use for Viterbi decoding on explicit DHMM with D = 60 and two aggregate

Discussion and Conclusion

We have discussed implementation of Baum-Welch and Viterbi algorithms in linear memory. We successfully used 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 "real-time" 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).

Availability & requirements

University of New Orleans bioinformatics group: http://logos.cs.uno.edu/~achurban

Appendices

Appendix A – Proofs of scaling relationships

The scaling steps in Figure 2 need additional rigorous justification. Our proofs are partially inspired by recurrences presented in [24] with further clarifications.

Theorem 2 $\tilde{\beta }$ _t (m) = d _t $\overline{\beta }$ _t (m)

Proof Let us define $D_t=\frac{1}{{\displaystyle {\sum }_{i=1}^N{\beta }_t(i)}},d_t=\frac{1}{{\displaystyle {\sum }_{i=1}^N{\beta }_t(i)}},{\tilde{\beta }}_t(m)=D_t{\beta }_t(m)$,

$$\begin{array}{c}{\overline{\beta }}_t(m)={\displaystyle \sum _{j=1}^N{a}_{m,j}{b}_j(o_{t+1}){\tilde{\beta }}_{t+1}(j)}\\ =D_{t+1}{\displaystyle \sum _{j=1}^N{a}_{m,j}{b}_j(o_{t+1}){\beta }_{t+1}(j)}\\ =D_{t+1}{\beta }_t(m),\\ d_t=\frac{1}{{\displaystyle {\sum }_{i=1}^N{\overline{\beta }}_t(i)}}=\frac{1}{D_{t+1}{\displaystyle {\sum }_{i=1}^N{\beta }_t(i)}},\\ {\tilde{\beta }}_t(m)=d_t{\overline{\beta }}_t(m)=d_t{D}_{t+1}{\beta }_t(m)=\frac{1}{D_{t+1}{\displaystyle {\sum }_{i=1}^N{\beta }_t(i)}}{D}_{t+1}{\beta }_t(m)=D_t{\beta }_t(m).\end{array}$$

Here we observe useful relationships D _t = d _t D_{t + 1}and $\overline{\beta }$ _t (m) = D_{t + 1}β _t (m) necessary in follow-up proves.   □

Theorem 3 $\tilde{T}$_{i, j}(t, m) = d _t $\overline{T}$_{i, j}(t, m)

Proof Let us define $\tilde{T}$_{i, j}(t, m) = D _t T_{i, j}(t, m),

$$\begin{array}{c}{\tilde{T}}_{i,j}(t,m)=d_t{\overline{T}}_{i,j}(t,m)\\ =d_t\left[{\tilde{\beta }}_{t+1}(j)a_{m,j}{b}_j(o_{t+1})\delta (i=m)+{\displaystyle \sum _{n=1}^N{a}_{m,n}{\tilde{T}}_{i,j}(t+1,n)b_n(o_{t+1})}\right]\\ =d_t\left[D_{t+1}{\beta }_{t+1}(j)a_{m,j}{b}_j(o_{t+1})\delta (i=m)+D_{t+1}{\displaystyle \sum _{n=1}^N{a}_{m,n}{T}_{i,j}(t+1,n)b_n(o_{t+1})}\right]\\ =d_t{D}_{t+1}\left[{\beta }_{t+1}(j)a_{m,j}{b}_j(o_{t+1})\delta (i=m)+{\displaystyle \sum _{n=1}^N{a}_{m,n}{T}_{i,j}(t+1,n)b_n(o_{t+1})}\right]\\ =d_t{D}_{t+1}{T}_{i,j}(t,m)\\ =D_t{T}_{i,j}(t,m).\end{array}$$

Theorem 4 $\tilde{E}$ _i (γ, t, m) = d _t $\overline{E}$ _i (γ, t, m)

Proof Let us define $\tilde{E}$ _i (γ, t, m) = D _t E _i (γ, t, m),

$$\begin{array}{c}{\tilde{E}}_i(\gamma ,t,m)=d_t{\overline{E}}_i(\gamma ,t,m)\\ =d_t\left[{\displaystyle \sum _{n=1}^N{b}_n(o_{t+1})a_{m,n}{\tilde{E}}_i(\gamma ,t+1,n)+{\overline{\beta }}_t(m)\text{ SCORE}(o_t,\gamma )\delta (m=i)}\right]\\ =d_t\left[D_{t+1}{\displaystyle \sum _{n=1}^N{b}_n(o_{t+1})a_{m,n}{E}_i(\gamma ,t+1,n)+D_{t+1}{\beta }_t(m)\text{ SCORE}(o_t,\gamma )\delta (m=i)}\right]\\ =d_t{D}_{t+1}\left[{\displaystyle \sum _{n=1}^N{b}_n(o_{t+1})a_{m,n}{E}_i(\gamma ,t+1,n)+{\beta }_t(m)\text{ SCORE}(o_t,\gamma )\delta (m=i)}\right]\\ =d_t{D}_{t+1}{E}_i(\gamma ,t,m)\\ =D_t{E}_i(\gamma ,t,m).\end{array}$$