libhmm: A Modern C++20 Library for Hidden Markov Models with Correct MLE Emission M-Steps

Technical Summary: libhmm – A Modern C++20 Library for Hidden Markov Models

Problem Statement
The paper identifies two primary gaps in the existing Hidden Markov Model (HMM) software landscape:

1. Lack of Embeddable C++ Infrastructure: Dominant HMM tools (e.g., moveHMM, fHMM, hmmlearn) rely on R or Python runtimes, making them unsuitable for embedding in production C++ systems without significant overhead. Existing C/C++ libraries (e.g., GHMM, HMMLib) are either abandoned, lack modern C++ idioms, or are restricted to discrete emissions.
2. Incorrect Emission M-Steps: Many existing packages utilize Method-of-Moments (MOM) approximations for the emission distribution M-step in the Baum-Welch algorithm when closed-form Maximum Likelihood Estimators (MLE) are unavailable. The authors argue this is a computational convenience that compromises statistical correctness, particularly for distributions like Gamma, Beta, Weibull, Negative Binomial, and Student-t, where MOM can converge to different fixed points than exact MLE.

Methodology
The paper describes libhmm, a zero-dependency C++20 library designed for parameter estimation, sequence decoding, and model selection.

• Inference Algorithms: The library implements full log-space forward-backward and Viterbi algorithms. Unlike scaled algorithms that require explicit normalization constants (which can underflow), libhmm uses the numerically stable log-sum-exp operation throughout, ensuring unconditional stability for sequences of arbitrary length.
• Emission Distributions: The library supports sixteen continuous and discrete emission distributions.
  • Closed-form MLE: Implemented for Gaussian, Log-normal, Exponential, Poisson, Rayleigh, Uniform, Discrete, and Von Mises (mean direction).
  • Newton-Raphson MLE: Implemented for Gamma, Beta, Weibull, Negative Binomial, and Von Mises (concentration parameter). These use MOM estimates as seeds and typically converge in 2–4 iterations.
  • ECME Algorithm: For the location-scale Student-t distribution, the library implements the Expectation/Conditional Maximization Either (ECME) algorithm (Liu and Rubin, 1994). This addresses the known failure of MOM kurtosis estimation to converge to the true MLE for near-Gaussian tails.
• Architecture:
  • Zero Dependencies: The library compiles with only the C++20 standard library, avoiding external dependencies like Eigen or GSL.
  • Polymorphism vs. Templates: The design uses runtime polymorphism (virtual methods) for emission distributions to allow mixing distribution types within a single model and simplify serialization/Python bindings, accepting a minor virtual call overhead.
  • SIMD Acceleration: Compile-time dispatch provides AVX-512, AVX2, SSE2, and ARM NEON support for the forward-backward inner loop, with scalar fallbacks.
  • Layered Design: The code is organized into eight layers with strictly one-way dependency flow to ensure maintainability and auditable dependency graphs.

Key Contributions

1. Correct MLE M-Steps: The derivation and implementation of correct maximum likelihood estimators for sixteen distributions, specifically replacing MOM approximations with Newton-Raphson and ECME methods.
2. Log-Space Stability: A complete implementation of forward-backward and Viterbi decoding in log-space to prevent numerical underflow without the bookkeeping complexity of rescaling.
3. Modern C++20 Implementation: A zero-dependency, embeddable library that supports modern build systems (CMake) and idioms, filling the void left by abandoned legacy C/C++ libraries.
4. Python Bindings: A companion package, pylibhmm, wrapping the C++ API via nanobind to provide correct M-step capabilities to Python users.

Results and Benchmarks
The paper compares libhmm against established R reference packages and other C/C++ libraries across five real-data benchmarks and synthetic data:

• Throughput (Synthetic): On discrete models, libhmm is slower than HMMLib (approx. 2.8×) and GHMM (approx. 4.9×) on raw forward-backward throughput. The authors attribute this to libhmm's use of log-space (vs. scaled) and runtime polymorphism. However, the authors note this gap narrows significantly for continuous distributions where emission evaluation dominates the computation.
• Statistical Accuracy (Real Data):
  • DAX Financial Regimes: On a 3-state Student-t HMM for DAX returns, the ECME implementation achieved a log-likelihood of 17,487.2, surpassing the fHMM reference (17,485.7) and the previous MOM-based implementation (17,334.9) by 152 nats. This confirms ECME finds a superior solution compared to MOM.
  • Elk Movement: Parameters fitted to GPS tracks matched moveHMM reference values within <2% difference.
  • Wind Direction: A Von Mises HMM correctly handled circular data boundaries (0°/360°), whereas a Normal approximation (used by HiddenMarkov) resulted in 100% misclassification for observations near the boundary.
• Performance (Real Data): libhmm demonstrated significant speedups over R packages due to compiled C++ execution and lack of interpreter overhead (e.g., 23× faster on Elk movement data, 5–12× faster on DAX regimes).

Significance and Claims
The paper positions libhmm as a distinct tool in the HMM ecosystem: it is the only C++ library combining broad distribution support (16 types), statistically correct M-step implementations (ECME/Newton-Raphson), and modern, zero-dependency build conventions.

The authors explicitly acknowledge architectural tradeoffs: the library sacrifices raw throughput on discrete models (where GHMM/HMMLib excel) to gain API usability, distribution breadth, and numerical correctness. The primary significance lies in providing a "correct" and "embeddable" alternative for researchers and engineers integrating HMMs into production C++ pipelines, particularly for continuous emission distributions where standard MOM approximations are statistically suboptimal. The library is released under the MIT license with Python bindings available via pylibhmm.