Fast and scalable inference in hidden Markov models with Gaussian fields

1. Il Problema

I Modelli a Markov Nascosti (HMM) sono strumenti potenti per l'analisi di serie temporali guidate da stati discreti non osservati, ampiamente utilizzati in ecologia, medicina ed economia. Tuttavia, la crescente complessità dei dati empirici richiede l'integrazione di strutture latenti aggiuntive, come effetti spaziali, trend temporali o perturbazioni di misurazione.

Le Gaussian Fields (GF) o Processi Gaussiani (GP) offrono un framework flessibile per modellare queste strutture latenti complesse. Sebbene l'approccio SPDE (Stochastic Partial Differential Equation) di Lindgren et al. (2011) abbia reso i GP computazionalmente gestibili trasformandoli in Campi Aleatori Markoviani Gaussiani (GMRF) con matrici di precisione sparse, l'integrazione di questi campi negli HMM presenta una sfida fondamentale:

• Gli HMM richiedono la marginalizzazione sugli stati latenti discreti per calcolare la verosimiglianza.
• I metodi di inferenza frequentista moderni (come quelli basati su TMB/RTMB) utilizzano l'approssimazione di Laplace e la differenziazione automatica. Questi metodi richiedono il calcolo della matrice Hessiana della log-verosimiglianza rispetto alle variabili latenti.
• Mentre la parte GMRF produce una Hessiana sparsa, la struttura classica degli HMM genera una Hessiana densa (piena) a causa delle dipendenze temporali globali introdotte dall'algoritmo forward. Questo rende il calcolo proibitivo in dimensioni elevate, vanificando i vantaggi della sparsità degli SPDE.

2. Metodologia

L'autore propone una soluzione che combina tre componenti chiave per rendere l'inferenza veloce e scalabile:

A. Approssimazione di Laplace e SPDE

Il framework utilizza il pacchetto R RTMB (Template Model Builder) per eseguire l'approssimazione di Laplace. Le Gaussian Fields sono approssimate tramite l'approccio SPDE su una mesh triangolata, risultando in una densità GMRF con una matrice di precisione sparsa ($Q_{\tau, \kappa}$).

B. L'Algoritmo Forward a Banda (Banded Forward Algorithm)

Il contributo metodologico principale è una modifica all'algoritmo forward classico per il calcolo della log-verosimiglianza degli HMM.

• Il problema: Nell'algoritmo forward standard, la variabile forward $\phi_t$ dipende da tutta la storia passata ($x_1, \dots, x_{t-1}$), creando dipendenze globali che rendono la Hessiana densa.
• La soluzione: L'algoritmo viene modificato dividendo la serie temporale in blocchi non sovrapposti di lunghezza $k$ (parametro di banda).
  1. Il primo blocco è calcolato esattamente.
  2. Per i blocchi successivi ($b \geq 3$), la dipendenza dalla storia precedente viene "troncata". Invece di usare la variabile forward esatta (che dipende da tutti i dati precedenti), si utilizza una variabile forward approssimata ($\tilde{\phi}$) inizializzata con una distribuzione fissa (es. uniforme) e propagata solo attraverso il blocco precedente.
• Risultato: Questa troncatura locale garantisce che le derivate incrociate della log-verosimiglianza siano nulle per punti temporali sufficientemente distanti ($|t - t'| > k$). Di conseguenza, la matrice Hessiana risultante diventa a banda (banded) e quindi sparsa, permettendo l'uso efficiente di fattorizzazioni di Cholesky sparse.

C. Teorema di Convergenza

L'autore dimostra (Teorema 1) che l'errore di approssimazione decade geometricamente all'aumentare della banda $k$. Se la catena di Markov è uniformemente ergodica, l'errore è limitato da $C \rho^k$, garantendo che con una banda sufficiente, l'approssimazione sia numericamente indistinguibile dal calcolo esatto.

3. Contributi Chiave

1. Integrazione Scalabile: Ha risolto il collo di bottiglia computazionale che impediva l'uso di campi Gaussiani ad alta dimensionalità negli HMM tramite metodi frequentisti.
2. Algoritmo Ibrido: Ha sviluppato un algoritmo forward modificato che preserva la sparsità della matrice Hessiana, rendendo compatibile la struttura degli HMM con gli strumenti di ottimizzazione basati su TMB/RTMB.
3. Implementazione Pratica: L'approccio è implementato nel pacchetto R LaMa, che facilita l'uso di questi modelli complessi senza richiedere la scrittura di codice C++ personalizzato.

4. Risultati e Casi di Studio

L'efficacia del metodo è stata dimostrata attraverso simulazioni e due studi di caso reali:

• Caso 1: Rilevamento di Brillamenti Stellari (Stellar Flares)

  • Dati: Serie temporali di luminosità di una nana M (TIC 031381302) con 19.259 osservazioni.
  • Modello: Un HMM a 3 stati (quieto, in fiamme, in decadimento) con un trend quasi-periodico modellato come un Processo Gaussiano (SPDE).
  • Risultato: Il modello ha rilevato 18 eventi di brillamento in circa 6 minuti su un singolo core Apple M2.
  • Confronto: Questo rappresenta un miglioramento drastico rispetto ai 1-4 ore per blocco richiesti dai metodi Bayesiani (HMC) precedenti, permettendo l'analisi dell'intera serie temporale senza suddivisione in chunk.

• Caso 2: Movimento dei Leoni nel Kalahari Centrale

  • Dati: 53.453 posizioni GPS orarie di leoni.
  • Modello: HMM a 2 stati (riposo/attivo) dove le probabilità di transizione dipendono da covariate temporali (cicli diurni) e da un campo spaziale latente (GMRF).
  • Risultato: Il modello è stato adattato in circa 6 minuti. Ha rivelato una struttura spaziale non ovvia: una regione centrale dove i leoni hanno una probabilità molto bassa di passare allo stato di riposo, suggerendo driver ecologici sottili (es. distribuzione delle prede) non visibili nelle immagini satellitari.

5. Significato e Implicazioni

Questo lavoro è significativo perché:

• Democratizza l'uso dei GP negli HMM: Permette ai ricercatori di incorporare strutture latenti complesse (spaziali, temporali, non stazionarie) in modelli a stati nascosti mantenendo tempi di calcolo accettabili.
• Alternativa Scalabile ai Metodi Bayesiani: Offre un'alternativa frequentista efficiente ai metodi MCMC (come HMC), che spesso diventano impraticabili con grandi dataset o modelli ad alta dimensionalità.
• Flessibilità Applicativa: Il framework è applicabile a una vasta gamma di discipline, dall'ecologia del movimento all'astrofisica, permettendo di modellare fenomeni dove il "rumore" o le tendenze di fondo non sono semplici errori indipendenti ma strutture correlate complesse.

In sintesi, Fischer ha colmato il divario tra la flessibilità dei campi Gaussiani e l'efficienza computazionale degli HMM, fornendo uno strumento robusto per l'analisi di dati temporali e spaziali complessi.