Parameter estimation of hidden Markov models: comparison of EM and quasi-Newton methods with a new hybrid algorithm
Sidonie Foulon Thérèse Truong Anne-Louise Leutenegger Hervé Perdry


[
    [
    3 September 2024
====================







§ ABSTRACT


Hidden Markov Models (HMM) model a sequence of observations that are dependent on a hidden (or latent) state that follow a Markov chain. These models are widely used in diverse fields including ecology, speech recognition, and genetics.

Parameter estimation in HMM is typically performed using the Baum-Welch algorithm, a special case of the Expectation-Maximisation (EM) algorithm. While this method guarantee the convergence to a local maximum, its convergence rates is usually slow.

Alternative methods, such as the direct maximisation of the likelihood using quasi-Newton methods (such as L-BFGS-B) can offer faster convergence but can be more complicated to implement due to challenges to deal with the presence of bounds on the space of parameters.

We propose a novel hybrid algorithm, QNEM, that combines the Baum-Welch and the quasi-Newton algorithms. QNEM aims to leverage the strength of both algorithms by switching from one method to the other based on the convexity of the likelihood function.

We conducted a comparative analysis between QNEM, the Baum-Welch algorithm, an EM acceleration algorithm called SQUAREM (Varadhan, 2008, Scand J Statist), and the L-BFGS-B quasi-Newton method by applying these algorithms to four examples built on different models. We estimated the parameters of each model using the different algorithms and evaluated their performances.

Our results show that the best-performing algorithm depends on the model considered. QNEM performs well overall, always being faster or equivalent to L-BFGS-B. The Baum-Welch and SQUAREM algorithms are faster than the quasi-Newton and QNEM algorithms in certain scenarios with multiple optimum. In conclusion, QNEM offers a promising alternative to existing algorithms.




§ KEYWORDS


Hidden Markov models, Baum-Welch, quasi-Newton, L-BFGS-B, SQUAREM, optimisation, computational statistics



§ INTRODUCTION


Hidden Markov Models (HMMs) are widely used in time series analysis, with application across diverse fields, including ecology, speech recognition and many others. A particular interest is found in biology, especially in genetics: HMM are often used for genome annotation. For example, CpG islands identification (regions of the genome with a high frequency of cytosine-guanine dinucleotides) <cit.>, can be performed through the use of HMM. It can also be used to differentiate exons (protein-coding sequences) and introns (non-coding sequences) within genes <cit.>. 
. 

To compute the maximum likelihood estimators of the parameters of HMM, two approaches are naturally employed: the Baum-Welch algorithm <cit.>, which is a special case of the Expectation-Maximisation (EM) algorithm <cit.> (cf section <ref> below) and direct maximisation of the likelihood using a quasi-Newton method (section <ref>), especially BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm <cit.>. Both methods have advantages and drawbacks. In particular, the Baum-Welch algorithm accommodates naturally the presence of bounds on the space of parameters, which can be a hurdle for the quasi-Newton, in particular when some parameters live in an open interval. On the other hand, the Baum-Welch algorithm usually converges more slowly than the quasi-Newton method.

In this paper we propose a new algorithm, QNEM, for HMM parameter estimation. This algorithm starts with one or more EM steps, and switches to a quasi-Newton algorithm, the BFGS, as soon as a criterion on local convexity of the log-likelihood is met. Quasi-Newton iterations are then performed, until convergence or until the aforementioned criterion on local convexity is no longer satisfied, in which case it switches back to the EM algorithm, etc.

We compare its performances with several algorithms: the Baum-Welch algorithm, an accelerated EM algorithm called SQUAREM (section <ref>) and the L-BFGS-B (Limited-memory BFGS with Bound constraints) algorithm, the “off-the-shelf” quasi-Newton method optimized for problems with bounded parameters.

We first review the theory of HMMs, and then expose the four algorithms of interest as well as the examples on which the comparisons are performed. Finally, we present the results obtained and discuss them.


    

§ HIDDEN MARKOV MODELS


    
A sequence of random variables measured at successive moments (X_i)_i ∈ℕ is a Markov chain if it satisfies the Markov property: for any time i, the distribution of X_i+1 conditional to (X_j)_j≤ i is equal to the distribution of X_i+1 conditional to the single value X_i. 
The model is classically illustrated by a directed acyclic graph (DAG) in Figure <ref>.




One of the property of a Markov chain is the homogeneity of the sequence. A sequence is homogeneous if the transition probabilities ℙ(X_i+1 | X_i) are identical in every point of the sequence. This property is sometimes contradicted by a given sequence of observations, making the Markov model inappropriate. A more complex model, in which there is an additional hidden layer (S_i)_i ∈ℕ which is a Markov chain, and the distribution of X_i is assumed to depend only on the “hidden state” S_i, can be more appropriate. This model is called Hidden Markov Model (HMM) and is illustrated by the DAG in Figure <ref>.



In some situations, the Markov model is known to be adequate, but it can be impossible to directly access the Markov chain, but only variables depending on the current state of the Markov chain. In this case, the HMM provides a natural framework for analysis. HMMs allow to retrieve the states S_i of the latent Markov chain, using the observations X_i whose distribution depends on the hidden states. 


        

 §.§ Direct maximisation of the likelihood with a quasi-Newton type method
 

Direct maximisation of the likelihood (DML) with quasi-Newton methods is a possible way to estimate the HMM parameters. 

Quasi-Newton methods are particularly suited for maximisation of differentiable scalar functions, such as likelihood functions. They have the advantage of avoiding the calculation of the Hessian matrix by approximating the Hessian matrix (or its inverse) at each iteration. 
The off-the-shelf methods are BFGS or L-BFGS-B when some parameters are constrained, standing respectively for Broyden-Fletcher-Goldfarb-Shanno <cit.>, and for Limited-memory BFGS with Bound constraints. 

This method relies on a forward step, computing probabilities of the hidden states, either conditionally to the X_i, or jointly with the X_i. These probabilities allow to compute the likelihood. The most common choice for the DML is using joint probabilities going forward <cit.>. A log-transformation is often employed to get rid of the numerical underflow implied by the joint probabilities <cit.>. In this paper, we chose to use the conditional version of the algorithm for the forward step (cf below section <ref>), which avoids the numerical underflow and therefore going through the log-transformation. 



        

 §.§ Baum-Welch algorithm, an EM algorithm
 





The Baum-Welch algorithm <cit.> was developed to estimate the HMM parameters.

The respective contributions of Baum and Welch in the development of the Baum-Welch algorithm is explained in Welch's Shannon lecture, in 2003 <cit.>. After working separately in the calculation of a posteriori probabilities of the hidden states, Baum and Welch worked together on the re-estimation of the parameters from these a posteriori probabilities. 

The Baum-Welch algorithm was later recognised as a particular case of the more general Expectation-Maximisation (EM) algorithm formalised by Dempster, Laird and Rubin in 1977 <cit.>. 

The Baum-Welch algorithm relies on a forward-backward algorithm for calculating the a posteriori probabilities of the hidden states (E-step). This algorithm is constructed with a forward step as sketched before. We chose to use the forward algorithm based on conditional probabilities, which is more commonly used for the Baum-Welch algorithm. It is also possible to use the forward algorithm based on the joint probabilities, together with an adequate version of the backward step (cf <cit.>, chapter 6).
Then the parameters are re-estimated from the a posteriori probabilities (M-step). This step depends on the model considered.







§ MATERIAL AND METHODS





For the following, let X = (X_1, ..., X_L) = X_1^L the observations and S = (S_1, ..., S_L) = S_1^L the hidden states along the chain of length L.

Let (s, t) = ℙ_θ(S_i+1 = t| S_i = s) be the transition probability from state s to state t and (x,s) = ℙ_θ(X_i = x| S_i = s) the emission probability of observation x from state s.

We consider four maximum likelihood methods to estimate the parameter θ.



 §.§ Direct maximisation of the likelihood


As seen in for example Zucchini <cit.> and Turner <cit.>, a variant of the forward algorithm based on joint probabilities can be used to compute the likelihood. If a_i(s) = ℙ(S_i = s, X_1^i-1 = x_1^i-1) and b_i(s) = ℙ_θ(S_i = s,X_1^i = x_1^i), the quantities a_i(s) and b_i(s) can be computed recursively as follows:


    
  * Initialisation : 

    
    a_1(s) = π(s)

    
    b_1(s) = a_1(s) (x_1, s)

    with π(s) the initial distribution of state s.
    
  * Recursion : 

    
    a_i(s) = ∑_t b_i-1(t) (t, s)

    
    b_i(s) = a_i(s) (x_i, s)

    for i ∈2, ..., L


Finally the likelihood can be computed as
    
    L(θ ; X) = ℙ_θ(X_1^L = x_1^L) = ∑_s a_L(s)


As the values of the a_i(s) and b_i(s) become very small when i increases, this algorithm is prone to numerical underflow. The usual solution is to compute the values on the logarithmic scale <cit.>, which involves rewriting the above equations accordingly. In this paper, we prefer to use the version of the forward algorithm based on conditional probabilities.

This forward algorithm computes classically two conditional probabilities alternatively: 

    α_i(s) = ℙ_θ(S_i=s | X_1^i-1 = x_1^i-1)

the “forecast” probability and

    β_i(s) = ℙ_θ(S_i=s | X_1^i = x_1^i)

the “filtering” probability. 

To compute the log-likelihood we need to compute also for each i
γ_i = log(ℙ_θ(X_1^i = x_1^i) ). These quantities can be computed recursively as follows:


    
  * Initialisation : 


    
    α_1(s) = π_s
 
    
    β_1(s) = (x_1,s) ×α_1(s)/∑_i (x_1,i) ×α_1(i)
 
    
    γ_1 = log(∑_s α_1(s) (x_1, s)).


    
  * Recursion : 

 
    
    α_i(s) = ∑_t (t,s) ×β_i-1(t)
 
    
    β_i(s) = (x_i,s) ×α_i(s)/∑_t (x_i,t) ×α_i(t)
 
    
    γ_i = γ_i-1 + log(∑_s (x_i, s) α_i(s))

    for i ∈2, ..., L

Finally, the log-likelihood ℓ(θ ; 𝐗) is equal to γ_L. The likelihood is still computed on the logarithmic scale, but the α_i(s) and β_i(s) can be computed without rewriting all the equations on the logarithmic scale, which make the computation simpler.

The L-BFGS-B algorithm is then performed on the likelihood computed in the last step using the R  function. The gradient of the likelihood along θ can be computed using automatic differentiation (our implementation relies on the R package salad <cit.>). 

    

 §.§ Baum-Welch algorithm
 

The Baum-Welch EM algorithm is composed by two steps:

    
  * E step: knowing θ, compute the a posteriori probabilities of the hidden states along the chain (forward-backward algorithm)
    
  * M step: re-estimate θ using these probabilities

The Baum-Welch algorithm alternates these 2 steps until convergence. As mentioned in section <ref>, we will use here the forward-backward algorithm presented in <cit.> or <cit.>. 

        

  §.§.§ E step : Forward-Backward algorithm
 

The forward-backward algorithm runs through the observations twice: the first time forward, ie. from first to last observations ; the second backward, ie. from last to first observation. 

After computing the forward quantities shown in equations (<ref>) and (<ref>), the following “backward quantities” will be computed:

    δ_i(s,t) = ℙ_θ(S_i-1=s, S_i=t | X_1^L = x_1^L)

the “smoothing" probability and

    φ_i(s) = ℙ_θ(S_i=s | X_1^L = x_1^L)

the “marginal" probability. 


    
  * Initialisation: 

    
    φ_L(s) = β_L(s)

    
  * Recursion:

    
    δ_i(s,t) = (s,t) ×β_i-1(s)/α_i(t)×φ_i(t)

    
    φ_i-1(s) = ∑_t δ_i(s,t)

for i ∈L, ..., 2.


At the end of the forward-backward, we have computed the probabilities φ_i(s) of the hidden states, and the probability δ_i(s,t) of two consecutive hidden states, given all the observations. These probabilities will be used in the M step.


It is worth to note that the backward algorithm has a higher computational cost than the forward algorithm, due to the necessity to calculate the δ_i(s,t) for all pairs of states s and t, at each time point i.
As mentioned previously, there exists a variant of the backward algorithm based the forward quantities a_i and b_i involving joint probabilities, as seen in <cit.>.

        

  §.§.§ M step


The aim in M step is to re-estimate the parameter θ that will be used in the following E step.

Note θ^(k) the parameter at the k-th iteration. For each iteration w of the Baum-Welch algorithm, in M step we are going to estimate the parameter θ^(k+1) for the next iteration.

The complete likelihood of θ, ie the likelihood assuming the hidden states are observed, is: 

    L(θ ; S = s, X = x) =  ℙ_θ(S_1 = s_1)    ·∏_i=2^L ℙ_θ(S_i = s_i | S_i-1 = s_i-1) 
       ·∏_i=1^L ℙ_θ(X_i = x_i | S_i = s_i)


Now the log-likelihood of θ is: 

    ℓ(θ ; S = s, X = x) =  logℙ_θ(S_1 = s_1)     + ∑_i=2^L logℙ_θ(S_i = s_i | S_i-1 = s_i-1) 
        + ∑_i=1^L logℙ_θ(X_i = x_i | S_i = s_i)



The M step relies on the expected value of the log-likelihood: 

    Q(θ ; θ^(k))     = 𝔼( ℓ(θ ; S = s, X = x) |  θ^(k)) 
       = ∑_s ( logℙ_θ(S_1 = s) ·φ_1(s) ) + ∑_i=2^L[ ∑_s,t( logℙ_θ(S_i = t | S_i-1 = s) ·δ_i(s,t) ) ] 
        + ∑_i=1^L [ ∑_s ( logℙ_θ(X_i = x_i | S_i = s) ·φ_i(s) ) ]


The M step ends with finding the parameters that maximise Q(θ ; θ^(k)) : 

    θ^(k+1) = _θ Q(θ ; θ^(k)).


In concrete examples, the δ_i(s,t) will be used to compute the transition parameters of the Markov chain, and the ϕ_i(s) the parameters relatives to the emission probabilities.

    

 §.§ EM acceleration
 

The SQUAREM method was presented in Varadhan, 2008 <cit.>. It is a method created to accelerate the EM algorithm. This method seems interesting when we aim to apply EM to high-dimensional data or data based on complex statistical models. In a nutshell, the SQUAREM makes two EM iterations, and uses these values to extrapolate the trajectory of the algorithm, and estimate a future point several iterations ahead. It then iterates again from this extrapolated point. The constraints on the parameters, if any, are simply dealt with by constraining the extrapolated point to be in the feasible space. We applied this method to our examples to test its performances as well.

    

 §.§ QNEM: a mix of the two algorithms

    
We propose to mix the Baum-Welch algorithm and the quasi-Newton algorithm BFGS as presented in the algorithm <ref>.  



If the current value of the parameter is θ^(k), the BFGS algorithm finds θ^(k+1) by performing a linear search in a search direction p_k = -H_k ∇ f_k where H_k is the current approximation H_k of the inverse of the Hessian of the objective function, and ∇ f_k is its gradient in θ^(k). The H_k is updated using the formula
        
    H_k+1 = H_k - H_k · y_k · y'_k · H_k/y'_k · H_k · y_k + s_k · s'_k/y'_k · s_k

where s_k = θ^(k+1) - θ^(k) and y_k = ∇ f_k+1 - ∇ f_k. If H_k is positive definite, H_k+1 will be positive definite as soon as the scalar product of s_k and y_k is positive: 

    s'_k · y_k > 0
    
This is the curvature condition. If the function is strongly convex it always holds after a linear search which successfully decreases the value of the objective function. In the general case, it may not hold, leading to the use of complex linear search methods to avoid this situation, such as the Nocedal line search, see <cit.> chapter 3. 

The QNEM starts by performing Baum-Welch algorithm steps, until the curvature condition (equation <ref>) is met; it then switches to BFGS algorithm. Our BFGS iterations use a simple backtracking search, until the Armijo condition is met. This condition ensures that the objective function decreased.

As long as the curvature condition holds, BFGS iterations are continued. If at some point it does not hold, the QNEM falls back to the Baum-Welch algorithm, until the curvature condition is met and the quasi-Newton iterations are performed again, until convergence.
The box constraint are simply dealt with by projecting the gradient on the box when θ is on the border.



    

 §.§ Stopping criterion


For the four algorithms, we choose to use the same stopping criterion to decide when the convergence point is reached. The criterion is the classic relative convergence tolerance (reltol), which is based on the likelihood of the current iteration and the previous one. At the k-th iteration, the algorithm stops if



    |ℓ^(k-1) - ℓ^(k)|/|ℓ^(k-1)| + reltol < reltol



where ℓ^(k) is the likelihood of the k-th iteration and reltol is a small constant, depending on the machine precision; in our experiment it is 1.49 · 10^-8. 
    
    

 §.§ Examples

        

  §.§.§ Umbrella example


We will first use a simple example based on a fictitious situation. Every day of the year, in a secret underground installation, a security guard observes whether the director comes in with an umbrella or not <cit.>. During 56 days, the security guard keep track of his observations on his note pad as follows: U if the director carries an umbrella; N otherwise. 
The guards aims to predict the daily weather conditions (hidden states) based on the umbrella observations: rainy (noted R) or dry (noted D) day. Here the hidden state is the weather. Since the security guard is in a underground installation, he can not observe directly the weather: he has to rely only on his observations of the umbrella status to determine the weather. The parameters to estimate are a the probability of weather state transition and b the probability of error in the umbrella status (not carrying an umbrella on rainy day or carrying an umbrella on dry day). The transition and emission matrix are presented respectively in Tables <ref> and <ref>.




To apply the quasi-Newton algorithm (L-BGFS-B) to this example, it was necessary to bound the values of a and b to the intervals [0.01, 0.99]. Using the interval [0, 1] would be preferable, but the likelihood is not always defined on the boundary, which leads to a failure of the algorithm. The algorithm can not accommodate open intervals such as (0,1).

        

  §.§.§ Old Faithful geyser example
 

Following an example presented in Zucchini <cit.>, we work with the Old Faithful geyser data <cit.>, available in the R package , considering the 272 duration of the eruptions. 
We make two uses of the data. First, we dichotomise at three minutes of eruption to divide the observations in two groups "Duration < 3min" and "Duration ≥ 3 min" (noted as "Dinf3" and "Dsup3" respectively) to work in a discrete observations case. This is similar to what is done in <cit.>. 
Second, we work with the continuous values, assuming they are drawn in Gaussian distributions whose parameters depend on the hidden states.
In the sequence of dichotomised values, it appears that each short eruption was followed by a long eruption, but there are sequences of long eruptions. We decided to assume that the hidden Markov chain has three states: "short", "long" and "steady long", noted respectively "S", "L" and "Sl". 
The parameters to estimate in the dichotomised model are a the probability to change from state "short" to "steady long", b the probability to stay in "steady long" state, c, d and e the probability to observe a "Duration ≥ 3 min" in respectively "short", "long" and "steady long" states. This emission matrix is presented in Table <ref>.
For the continuous model, the parameters to estimate are the same a and b and six other parameters μ_s, μ_l and μ_sl, the mean of the Gaussian distributions in which duration of eruptions are drawn, and σ_s, σ_l and σ_sl, their standard deviation. 
The transition matrix, common to both the dichotomised and continuous models is presented in Table <ref>



Again, for the quasi-Newton L-BFGS-B algorithm, we had to constraint the probabilities a, b, c, d, and e to the interval [0.01, 0.99], and to set a lower bound of 0.01 for  σ_s, σ_l, and σ_sl.
        
        

  §.§.§ HBD segments example

        
Consanguineous individuals are offspring of related parents. In the absence of pedigree information, consanguineous individuals can be identified through their genome which carries Homozygous-By-Descent (HBD) segments. HBD segments are regions of the genome were the two copies of the genome are identical and come from a common ancestor of the parents. These segments are characteristic of consanguineous individuals. Reconstructing these regions allow to determine whether an individual is consanguineous and his/her degree of consanguinity. The coefficient of consanguinity f can be defined as the expected proportion of HBD genome. 
HMM are used to compute for each variant of the genome, its probability to be in a HBD segment. These probabilities can be used to reconstruct the HBD segments, but this requires to first estimate the parameters of the HMM. The parameters to estimate are the coefficient of consanguinity f and another parameter a related to the mean length of HBD segments <cit.>.
cf. Figure <ref>



To test the algorithms, we will generate the genome of an individual from a given set of parameters θ = (f = 0.0625, a = 0.064), which corresponds to the expected values for an offsprings of first cousins. The data were created from an inhouse R package named  (available at <https://github.com/genostats/Mozza>, still in development). We then rounded the distances between positions to 0.1 centiMorgan (cM). For each of these new positions, we randomly selected only one observation, to create submaps of positions. This corresponds to a total of 1050 observations.


The transition and emission probabilities are presented in Tables <ref> and <ref>. We note d the distance between positions (fixed here at 0.1 cM), pA and pa the allele frequency of the reference and alternative allele respectively and ϵ the genotyping or mutation error. 



In this example, the probability f was constrained to the interval [0.01, 0.99] and the lower-bound was 0.01 for the parameter a when applying the quasi-Newton algorithm.

    

 §.§ Comparison of the algorithms

We compare the likelihood maximisation methods by comparing 

    
  * the number of iteration to convergence
    
  * the number of forwards step runs
    
  * the number of backwards step runs
    
  * the running time
    
  * the percentage of runs that failed to converge
    
  * the running time of forward step
    
  * the running time of backward step
    


We also studied the final estimations of the parameters and calculated the percentage of appearance of each corresponding likelihoods.

For all four examples, the tests were made on 1000 uniform draws within a wide possibility of values for the initial θ, some of them being unlikely.







§ RESULTS


    

 §.§ Umbrella example


The results are presented in Table <ref>. 



We can see that, overall, the QNEM algorithm achieves convergence with the lowest mean number of iterations, averaging 11.32 across all runs. However, the quasi-Newton algorithm has the lowest maximum number of iterations (37 iterations vs. 118-338 iterations for the other three algorithms). The number of iterations of the Baum-Welch and the SQUAREM algorithms are comparable and systematically higher than those of the QNEM and the quasi-Newton algorithms.

Apart from the total number of iterations, it is interesting to look at the mean number of forward and backward steps used by the four algorithms. The quasi-Newton algorithm does not perform any backward steps. For the quasi-Newton and Baum-Welch algorithms, the number of forward steps corresponds to the mean number of iterations on the 1000 draws of the initial starting point. The Baum-Welch algorithm compute the same number of forward and backward steps. For the SQUAREM algorithm, the mean number of forward and backward steps is slightly different, with 17.70 forward steps and 17.67 backward steps. Finally, the QNEM algorithm executes an average of 12.63 forward steps and 2.80 backward steps, where the backward steps correspond to the calls to the Baum-Welch algorithm. We note that a forward step takes about 1.5ms for the Baum-Welch and the SQUAREM algorithms and around 11.5ms for the other two, which is longer because of the additional computation of the derivatives. A backward step takes about 8ms for all algorithms. 
The total mean time spent in forward and backward steps for a run is almost the same for the quasi-Newton algorithm and the QNEM (∼ 18ms). It is higher for the SQUAREM, and even longer for the Baum-Welch algorithm (∼ 39ms and 61.5ms respectively). This trend is in accordance with the total mean execution time, where the Baum-Welch algorithm is the slowest one and the quasi-Newton and the QNEM algorithms are the fastest.






We also examined the final values obtained for θ. They can be divided into three subsets, corresponding to three different regions of the likelihood as shown in Figure <ref>. Two of these regions are points with a = 0.06 and b = 0.18 or b = 0.82, corresponding to the highest likelihood value (negative log-likelihood equal to 33.3). The third region is the segment b = 0.5 and a ∈ [0.5, 1] which corresponds to a crest line (negative log-likelihood equal to 33.8). The parameters and the corresponding proportion of runs converging to these parameters for each algorithm are presented in Table <ref>. The θ with the best likelihood is reached almost 90% of the time with the quasi-Newton algorithm but only 70% of the time for the other three algorithms. 



    

 §.§ Old Faithful geyser example


We will first present the results obtained with the dichotomised data, then with the continuous values.


        

  §.§.§ Dichotomised values

        
As previously done, we tested the algorithms on 1000 random initial sets of parameters. The results are presented in Table <ref>. 





The quasi-Newton and the QNEM algorithms had a lower number of iterations compared to the SQUAREM and the Baum-Welch algorithms. The latters showed a mean number of iterations almost 6 times higher. The lowest number of iterations was always reached by the QNEM. 

Moreover, the quickest algorithm was the QNEM as well, with a mean of 0.26 seconds by run, compared to the quasi-Newton, the Baum-Welch and the SQUAREM algorithms, which had respectively 0.33 seconds, 1.21 seconds and 0.71 seconds. The mean time spent in the forward and backward step is the lower for the QNEM as well, with 215ms compared to 280ms for the quasi-Newton algorithm, 1197ms for the Baum-Welch algorithm and 693ms for the SQUAREM.

We note that the QNEM did not converge for two runs. 




We then present the table of the convergence points in Table <ref>. We see that the optimum is obtained most of the time by all algorithms but it represents only ∼ 50% of the runs for all algorithms except the QNEM, for which it's closer to 40%. There are two other local optimum reached by the algorithms. Moreover, we observe that the QNEM converges to other local optima not reached by the other algorithms.

        

  §.§.§ Continuous values


For the continuous case, the summary of the results is presented in Table <ref>.




The quasi-Newton algorithm showed the biggest amount of iterations, with a maximum of 159 iterations. The QNEM keeps showing the less number of iterations globally. 

In term of running time, the lowest are reached by the EM-based algorithms: the Baum-Welch algorithm and the SQUAREM, with a mean running time near 0.23s. Therefore, the quasi-Newton was slow on this continuous case, with a mean running time of 1.47s. The hybrid QNEM showed an intermediate time of 0.66s. This difference of time can be explained by the time spent in the forward step, which is 1.5 millisecond for the forward used in the SQUAREM and the Baum-Welch (conditional probabilities forward) against 17.5 milliseconds for the forward used in the other two algorithms (conditional probabilities forward with derivatives) for a comparable mean number of forward steps.

The quasi-Newton algorithm did not converge for 68 sets of initial parameters and the QNEM for 5 sets of parameters.





In Table <ref>, the best likelihood is reached most of the time by the Baum-Welch, the SQUAREM and the QNEM algorithms but not by the quasi-Newton algorithm. The latter converges more often to local optima. Most of the time, the quasi-Newton converge to the point with the third best likelihood. This behavior can explain its slow convergence in this example.

    
    

 §.§ HBD segments example


For this last example, the results are presented in Table <ref>. 



The smallest number of iterations is always hit by the QNEM, with a maximum of 27 iterations. It is closely followed by the quasi-Newton algorithm. On the opposite, the Baum-Welch algorithm shows less good performances in term of number of iterations with a mean of 68.53 iterations. The SQUAREM algorithm, is a bit better than the Baum-Welch algorithm with a mean of 56.19 iterations. 

In term of running time, the QNEM and the quasi-Newton shows a similar mean running time of 0.32s, whereas the SQUAREM and the Baum-Welch algorithms are way longer with respectively 1.37s and 2.40s of running time. 
These running times are linked to the time spent in forward and backward steps. Even if the time spent in the forward step by the two EM-based algorithms (5.3ms) is lower than the time spent is the QNEM and quasi-Newton's forward step (around 20ms), the two latter need less forward steps (respectively a mean of 12.83 and 15.36 iterations) than the Baum-Welch and the SQUAREM algorithms (respectively a mean of 68.53 and 38.66 iterations).

For all algorithms, only one likelihood optimum is reach with all 1000 draws of the starting value theta θ. In fact, this likelihood admits only one optimum.





§ DISCUSSION


The Baum-Welch (or EM) algorithm is a natural choice when it comes to HMM parameter estimation, offering several advantages: its conceptual simplicity as well as the ease of implementation, and the fact it naturally accommodates the possible constraints on the parameters. However, while the first iterations may move quickly in the space of parameters, once it is close to the solution its convergence is usually slow. 

Another possibility is to use a direct maximisation of the likelihood, particularly via a quasi-Newton method. It is expected to be more efficient than the Baum-Welch algorithm. However, the computation of the likelihood, and of its derivatives, on the logarithmic scale, can be tedious. In the literature, it is usually done using a joint probability version of the forward algorithm (equations <ref> to <ref>), in which all quantities need to be computed in the logarithmic scale. We showed that it is possible to use the conditional probability version of the forward algorithm (equations <ref> to <ref>) to compute the log-likelihood in a simpler manner. 

A practical problem with the use of the quasi-Newton algorithm is the presence of bounds on the parameters. Typically, in our examples, many parameters take their values in [0,1], and others in [0, +∞), with possibly an infinite value of the log-likelihood on the border of the space of parameters. The off-the-shelf quasi-Newton method for this case is L-BFGS-B, which tends to favour in its first iterations points on the border of the space, leading to a failure of the algorithm. A common solution, but not totally satisfying, is to use constraints such as [0.01, 0.99] instead.

An appealing alternative is to use an EM acceleration, such as the SQUAREM algorithm. It is simple to implement from an existing Baum-Welch implementation. Bounds on the parameters are easily dealt with, including when the likelihood is infinite on the border.

We propose a new algorithm, QNEM, which alternates between Baum-Welch algorithm iterations and BFGS iterations. Our intuition was that the Baum-Welch algorithm can help to explore efficiently the space of parameters, and the BFGS can allow to achieve a fast convergence once the algorithm is near to a local optimum. The algorithm switches between the two kind of iterations using a convexity criterion. The box-constraints are simply managed with a projection of the gradient on the box when the parameter θ gets on the border.

We compared the four algorithms on four different examples. The umbrella example (cf. sections <ref> and <ref>) is a toy example with two hidden states and two possible observed states. As expected, the quasi-Newton was much faster than the Baum-Welch algorithm. The SQUAREM succeeds in accelerating the Baum-Welch algorithm, but is still much slower than the quasi-Newton. Our algorithm, QNEM, is equivalent to the quasi-Newton in term of computational time.

Our second example was based on the Old Faithful data, with three hidden states and two possible observed states (cf. sections <ref> and <ref>). In that case, the QNEM is the faster of all algorithms, followed by the quasi-Newton, then by the SQUAREM, and the Baum-Welch is the slowest.

In our third example, the Old Faithful data were used again, with the same model for the hidden Markov chain, but now the observed states are continuous (cf. sections <ref> and <ref>). Surprisingly, for this model, the quasi-Newton is by far the slowest method. The Baum-Welch algorithm is faster, and the SQUAREM is even slightly faster. The QNEM ranks third, with intermediate performances.

Our last example, which was our motivating example, addresses a concrete question in genetics: making inference about HBD segments on the genome. It has two hidden states, and the observations are discrete (cf. sections <ref> and <ref>). In that case, just as in the first example, the QNEM and the quasi-Newton were equivalent, and performed better than the SQUAREM and the Baum-Welch in term of computational time. Our results are in line with those of <cit.>, who also compared the L-BFGS-B algorithm to the EM algorithm in term of number of iterations and running time on a more complex HMM model aimed to identify HBD segments. In their work, the quasi-Newton method is several time faster than the EM algorithm.

In all examples except the HBD segments example, there are multiple local likelihood optimum, leading to multiple final estimations of the parameters. In a concrete situation, to deal with this issue, one would need to run the algorithms multiple times from different initialisation of the parameters, to retain the highest likelihood. In our examples, the quasi-Newton and the QNEM seems more prone to converge to a local optimum which is not a global optimum. It is worth to note that in our numeric experiments, the starting points are drawn in very large intervals, which may exacerbate this behaviour.

There was no algorithm uniformly more efficient. In particular, it was surprising to us to find that on our third example, the Baum-Welch and the SQUAREM are the more efficient ones. The presence of multiple optima in this example seems to play a role in this behaviour. The QNEM was the most efficient on our second example, and equivalent to the quasi-Newton on the first and fourth example. This new optimization algorithm for HMM parameter estimation seems to be an interesting alternative to the other methods, offering a good compromise in most situations.










§ CODE AVAILABILITY

Related to this paper, we created a R package available at <https://github.com/SidonieFoulon/steveHMM>.



§ ACKNOWLEDGMENT


This work is part of the Inserm cross-cutting program GOLD (GenOmics variability in health & Disease)



§ FUNDING SOURCES


Sidonie Foulon is funded for this work by a doctoral grant from the French Ministry of Research. 


plainnat



