Simultaneous model discovery and state estimation under high data corruption

1. Il Problema

Il documento affronta la sfida scientifica di scoprire le equazioni differenziali ordinarie (ODE) che governano un sistema fisico partendo da dati osservativi che sono sia rumorosi (contenenti errori di misura) sia incompleti (con valori mancanti o regioni non osservate).

I metodi esistenti presentano diverse limitazioni in questo contesto:

• Regressione sui residui (es. SINDy): Questi metodi sono estremamente sensibili al rumore perché richiedono la derivazione diretta dei dati.
• Stima dei parametri classica: Spesso ignora l'incertezza del modello, portando a problemi di ottimizzazione non convessi e instabili.
• Inferenza dello stato: In sistemi caotici, integrare un modello scoperto con condizioni iniziali rumorose porta a traiettorie che divergono rapidamente dallo stato reale.

L'obiettivo è quindi sviluppare un metodo che esegua simultaneamente la scoperta del modello (identificazione dei termini dell'equazione) e la stima dello stato (ricostruzione della soluzione vera), gestendo efficacemente alti livelli di corruzione dei dati.

2. Metodologia

Gli autori propongono una strategia di regressione sparsa ibrida basata su principi statistici e ottimizzazione di secondo ordine.

Formulazione del Problema

Il metodo si basa su una funzione di verosimiglianza (likelihood) che tiene conto delle incertezze sia nei dati ($\hat{u}$) sia nel modello ($N(u; \theta) = 0$).
La funzione obiettivo da minimizzare è definita come:
$$L_\lambda(\theta, u) = \|N(u; \theta)\|^2 + \lambda \|\hat{u} - u|_D\|^2$$
Dove:

• $u$ è il vettore di stato stimato.
• $\theta$ sono i parametri incogniti (coefficienti dei termini dell'equazione).
• $N(u; \theta)$ rappresenta l'operatore differenziale discretizzato.
• $\lambda$ è un parametro che bilancia l'incertezza relativa tra il modello e i dati.
• Il termine di errore sui dati è limitato al sottoinsieme $D$ dove le osservazioni sono disponibili.

Questa formulazione "ibrida" interpola tra la regressione ai minimi quadrati classici (quando $\lambda$ è piccolo) e la minimizzazione del residuo puro (quando $\lambda$ è grande), permettendo di inferire sia i parametri che lo stato latente.

Selezione Sparsa e Criteri di Informazione

Per garantire che il modello scoperto sia parsimonioso (evitando l'overfitting), viene imposta la sparsità sui parametri $\theta$.

• Penalità: Viene utilizzata un'approssimazione liscia della norma $L_0$ (conteggio dei termini non nulli) per rendere il problema ottimizzabile tramite gradienti.
• Algoritmo di Selezione: Invece di fissare soglie rigide, gli autori utilizzano i Criteri di Informazione Bayesiana (BIC). L'algoritmo procede in due fasi:
  1. Rifinitura in avanti (Forward Refinement): Rimuove iterativamente un blocco di termini con coefficienti piccoli finché il BIC non peggiora.
  2. Recupero all'indietro (Backward Recovery): Se la rimozione di un singolo termine peggiora il BIC, l'algoritmo torna indietro e aggiunge nuovamente termini per esplorare configurazioni vicine al minimo locale.

Ottimizzazione Numerica

La minimizzazione della funzione di perdita ad alta dimensionalità e non convessa viene risolta utilizzando una variante di secondo ordine del metodo di Levenberg-Marquardt.

• Derivazione Automatica (AD): I gradienti e gli Hessiani sono calcolati automaticamente.
• Sfruttamento della Sparsità: Sfruttando la struttura sparsa dell'Hessiano derivante dalla discretizzazione spaziale/temporale, viene utilizzata una tecnica di colorazione dei grafi per calcolare efficientemente i prodotti vettore-Hessiano, rendendo fattibile l'ottimizzazione di secondo ordine anche per sistemi complessi.

3. Contributi Chiave

1. Approccio Ibrido Statistico: Unificazione della scoperta del modello e dell'estimazione dello stato in un unico quadro di massima verosimiglianza, superando le limitazioni dei metodi basati solo sui residui.
2. Selezione del Modello basata su BIC: Sostituzione delle soglie arbitrarie con criteri statistici rigorosi (BIC) per la selezione dei termini, migliorando la robustezza nella scelta della struttura del modello.
3. Ottimizzazione di Secondo Ordine Efficiente: Dimostrazione pratica dell'uso di metodi di Newton/Levenberg-Marquardt con Hessiano esatto e derivazione automatica, superando i metodi di primo ordine (come L-BFGS) in termini di velocità di convergenza e precisione.
4. Gestione di Dati Incompleti e Rumorosi: Il metodo non richiede filtraggio diretto dei dati o interpolazioni complesse preliminari; l'inferenza dello stato avviene simultaneamente alla scoperta del modello, interpolando naturalmente tra i dati osservati e la soluzione liscia dell'ODE.
5. Generalità: Il metodo non è limitato a modelli lineari nei parametri, permettendo di scoprire equazioni con termini non lineari complessi (es. esponenziali con parametri ignoti).

4. Risultati Sperimentali

Il metodo è stato testato su diversi sistemi dinamici classici (Van der Pol, Lorenz, Lorenz-96, Oscillatore Colpitts) confrontandolo con l'algoritmo WSINDy (lo stato dell'arte per dati rumorosi).

• Robustezza al Rumore: In presenza di rumore fino al 50%, il metodo proposto ha mostrato tassi di recupero della struttura del modello (True Positive Rate - TPR) e errori di parametri significativamente inferiori rispetto a WSINDy.
• Dati Incompleti: Il metodo ha dimostrato capacità di ricostruire accuratamente le traiettorie dello stato anche quando fino al 30-50% dei dati mancava (sia in blocchi continui che casuali).
• Sistemi Caotici: Nel sistema di Lorenz, il metodo è riuscito a recuperare lo stato fondamentale con alta precisione, evitando la divergenza tipica dei metodi di integrazione diretta su modelli rumorosi.
• Termini Non Lineari: Nel caso dell'oscillatore Colpitts, il metodo ha identificato con successo un termine esponenziale con un parametro ignoto ($a$ in $e^{ax}$), dimostrando flessibilità oltre le librerie polinomiali standard.
• Analisi di Griglia: Gli esperimenti hanno mostrato che un raffinamento della griglia di discretizzazione del modello migliora l'accuratezza fino a un certo limite, oltre il quale il guadagno satura o peggiora a causa della sovrabbondanza di variabili di ottimizzazione.

5. Significato e Conclusioni

Questo lavoro rappresenta un avanzamento significativo nel campo della scoperta di modelli basata sui dati.

• Superamento delle limitazioni attuali: Risolve il compromesso tra la necessità di filtrare i dati (che può distorcere la dinamica) e la sensibilità al rumore dei metodi di regressione diretta.
• Flessibilità: Offre un quadro generale applicabile a sistemi non lineari nei parametri e a dati parzialmente osservati, scenari comuni nelle scienze applicate reali (biologia, scienza dei materiali, fisica).
• Efficienza Computazionale: L'uso intelligente della sparsità e della derivazione automatica rende l'approccio scalabile, rendendo pratici metodi di ottimizzazione di secondo ordine che erano precedentemente considerati troppo costosi per problemi di identificazione di sistemi.

In sintesi, gli autori dimostrano che un approccio statistico rigoroso, combinato con tecniche di ottimizzazione avanzate, può recuperare dinamiche fisiche complesse da dati di scarsa qualità, offrendo uno strumento potente per la modellazione scientifica quando le leggi fisiche sono sconosciute o parzialmente note.