Fast and scalable inference in hidden Markov models with Gaussian fields

1. Problem Statement

Hidden Markov Models (HMMs) are widely used for analyzing time series data driven by unobserved, discrete states. However, empirical data in fields like ecology, medicine, and economics often contain complex latent structures such as spatial effects, temporal trends, or measurement perturbations that traditional parametric HMMs struggle to capture.

While Gaussian Fields (GFs) (or Gaussian Processes) offer a flexible framework for modeling such structured latent variation, integrating them into HMMs presents a significant computational challenge:

• The Integration Bottleneck: Frequentist inference requires evaluating the marginal likelihood, which involves integrating out high-dimensional latent variables.
• The Sparsity Conflict: Modern tools like Template Model Builder (TMB) and RTMB use the Laplace approximation combined with automatic differentiation to handle high-dimensional integrals efficiently. These tools rely on the Hessian matrix (second derivatives of the log-likelihood) being sparse to scale to large dimensions.
• The HMM Obstacle: Standard HMM likelihoods are computed using the forward algorithm, which creates global temporal dependencies. Consequently, the Hessian matrix with respect to the latent variables is dense, rendering standard sparse linear algebra techniques (and thus TMB/RTMB) computationally infeasible for high-dimensional GFs within an HMM framework.

2. Methodology

The paper proposes a unified framework that combines HMMs with Gaussian Fields using the Stochastic Partial Differential Equation (SPDE) approach and a modified forward algorithm.

A. Gaussian Fields via SPDE

Instead of using dense covariance matrices (which scale as $O(l^3)$), the authors utilize the Lindgren, Rue, and Lindström (2011) SPDE approach. This approximates Matérn Gaussian fields as Gaussian Markov Random Fields (GMRFs) on a triangulated mesh.

• Result: The precision matrix ($Q$) of the latent field is sparse, enabling efficient computation via sparse Cholesky factorization.

B. The Modified Forward Algorithm (Banded Forward Algorithm)

To resolve the density of the Hessian in the HMM component, the author modifies the standard forward algorithm:

1. Standard Issue: In a standard forward algorithm, the scaled forward variable $\phi_t$ depends on the entire history $x_1, \dots, x_{t-1}$. This creates a dense Hessian because any two time points are conditionally dependent.
2. The Solution: The time series is split into non-overlapping blocks of length $k$ (the bandwidth parameter).
   • Block 1: Computed exactly using the initial distribution.
   • Block 2: Computed exactly using the final state of Block 1.
   • Subsequent Blocks ($b \geq 3$): The recursion is initialized not with the exact history, but with an approximate scaled forward variable ($\tilde{\phi}$) derived by propagating a fixed prior vector $\rho$ (e.g., uniform distribution) through the immediately preceding block only.
3. Effect: This truncates the conditioning history. The likelihood contribution of a block depends only on the variables within that block and the previous block.
4. Theoretical Guarantee: Theorem 1 proves that the approximation error decays geometrically with the bandwidth $k$ due to the mixing properties of ergodic Markov chains.

C. Implementation

• Software: The method is implemented in the LaMa R package (for the forward algorithm) and integrated with RTMB (for automatic differentiation and Laplace approximation).
• SPDE: The fmesher package is used to construct the triangulated mesh and sparse precision matrices.
• Optimization: The marginal negative log-likelihood is minimized using gradient-based optimizers (e.g., nlminb), with the Laplace approximation providing the objective function and gradients.

3. Key Contributions

1. Sparse Hessian Construction: The primary contribution is the development of a banded forward algorithm that transforms the HMM log-likelihood into a form where the Hessian with respect to latent variables is banded (sparse). This bridges the gap between HMMs and sparse GMRF inference tools.
2. Scalable Inference: The approach enables the integration of high-dimensional Gaussian fields into HMMs using frequentist maximum likelihood estimation, a task previously limited by computational complexity.
3. Theoretical Justification: The paper provides a formal theorem demonstrating that the approximation error decreases exponentially with the bandwidth parameter $k$, ensuring the method is theoretically sound.
4. Practical Framework: The authors provide a complete, reproducible R-based workflow (LaMa + RTMB + fmesher) for fitting these complex models.

4. Results (Case Studies)

The authors validated the method through two distinct case studies:

Case Study 1: Stellar Flare Detection

• Context: Detecting stellar flares in photometric data (M dwarf TIC 031381302) where the signal is obscured by quasi-periodic oscillations.
• Model: A 3-state HMM (Quiet, Firing, Decaying) with an additive quasi-periodic trend modeled as a Gaussian Process.
• Performance:
  • Speed: The model fitted the full dataset ($T \approx 19,259$) in ~6 minutes on a single Apple M2 core.
  • Comparison: Previous Bayesian approaches (Hamiltonian Monte Carlo) required 1–4 hours per chunk and necessitated splitting the data. The proposed method handled the entire series without chunking.
  • Accuracy: Successfully detected 18 flare events while controlling for the quasi-periodic trend, outperforming simpler models that failed to distinguish low-energy flares from oscillations.

Case Study 2: Lion Movement Ecology

• Context: Modeling the movement of lions in the Central Kalahari to understand behavioral switching (Active vs. Resting) influenced by spatial location.
• Model: A 2-state HMM where transition probabilities depend on time of day (periodic functions) and a latent spatial Gaussian field.
• Performance:
  • Speed: Fitted $T \approx 53,453$ observations in ~6 minutes.
  • Insights: The model identified a specific region where lions are unlikely to transition to a resting state, suggesting subtle ecological drivers (e.g., prey distribution) not visible in satellite imagery. It also confirmed known diel patterns (nocturnal activity).

5. Significance

This paper represents a significant advancement in statistical modeling for time series and spatial data:

• Bridging Disciplines: It successfully merges the discrete state-switching capabilities of HMMs with the continuous, flexible modeling of Gaussian Fields.
• Computational Efficiency: By exploiting sparsity in the Hessian via the banded forward algorithm, it makes high-dimensional frequentist inference feasible where it was previously impossible or prohibitively slow.
• Scalability: The method scales linearly with the number of observations (due to the banded structure) and the number of spatial mesh nodes (due to the SPDE approach), allowing for the analysis of massive datasets in ecology, astronomy, and beyond.
• Future Potential: The framework opens the door for modeling complex phenomena involving non-stationary fields, non-Gaussian processes, and non-Euclidean manifolds within the HMM framework.