Scalable Gaussian Processes for Integrated and Overlapping Measurements Via Augmented State Space Models

Résumé technique : Processus gaussiens évolutifs pour des mesures intégrées et chevauchantes via des modèles d'espace d'état augmentés

Énoncé du problème
Les mesures astronomiques sont fréquemment obtenues par intégration de signaux sur des temps d'exposition finis (par exemple, l'accumulation de photoelectrons dans un CCD). Lorsque l'échelle de temps de la variabilité astrophysique sous-jacente est comparable à la durée d'exposition, les modèles standards de processus gaussiens (GP) instantanés deviennent des estimateurs biaisés. Une prise en compte correcte de cette intégration nécessite de modéliser la mesure comme une moyenne du processus latent sur l'intervalle d'exposition.

Bien que la régression GP soit un outil standard pour modéliser des phénomènes stochastiques (par exemple, les oscillations stellaires dans les données de vitesse radiale), la gestion de mesures intégrées pose des défis computationnels significatifs. Pour des expositions constantes et non chevauchantes, des noyaux spécialisés (comme ceux de celerite) peuvent maintenir une quasi-séparabilité, permettant des solutions en $O(N)$. Cependant, dans le cas général impliquant des temps d'exposition variables et, de manière cruciale, des expositions chevauchantes (fréquentes lors de la combinaison de données provenant de plusieurs instruments avec des cadences différentes), la matrice de covariance perd sa structure quasi-séparable. Cela force l'utilisation d'opérations matricielles denses, entraînant des coûts de mémoire en $O(N^2)$ et des coûts d'exécution en $O(N^3)$. Les méthodes existantes pour les noyaux intégrés (par exemple, Luhn et al. 2026) reposent sur la construction de la matrice de covariance dense complète, ce qui devient computationnellement prohibitif pour de grands ensembles de données ($N \sim 10^5$) typiques des relevés multi-instruments modernes.

Méthodologie
Les auteurs proposent d'étendre l'équivalence entre les processus gaussiens et les modèles d'espace d'état linéaires gaussiens (SSM) pour gérer des observations intégrées et chevauchantes tout en préservant l'évolutivité.

1. Formulation d'espace d'état augmenté :
   L'innovation centrale est l'augmentation du vecteur d'état du SSM standard. Pour un processus latent $x(t)$ régi par une équation différentielle stochastique linéaire (SDE), les auteurs introduisent un état intégral $z(t)$ tel que $dz/dt = x(t)$.

   • Le vecteur d'état augmenté devient $\tilde{x} = [x, z]^T$.
   • La dynamique du système est mise à jour pour inclure l'intégration de l'état latent.
   • Mécanisme de réinitialisation : Au début de chaque exposition $t_s$, l'état intégral $z$ est réinitialisé à zéro. À la fin de l'exposition $t_e$, l'observation $y$ est modélisée comme l'intégrale normalisée de l'état sur l'intervalle : $y = z(t_e)/\delta + \epsilon$, où $\delta$ est la durée d'exposition.

2. Gestion des chevauchements :
   Pour accommoder les expositions chevauchantes provenant de plusieurs instruments, le cadre attribue un état intégral distinct $z_i$ à chaque instrument $i$. Le vecteur d'état augmenté s'étend pour inclure tous les intégrales spécifiques aux instruments. Les matrices de réinitialisation sont appliquées de manière sélective : $z_i$ est réinitialisé uniquement au début d'une exposition appartenant à l'instrument $i$, permettant l'accumulation simultanée d'intégrales pour des instruments chevauchants. Cette construction encode naturellement les corrélations croisées requises pour les données chevauchantes sans nécessiter d'opérations matricielles denses.

3. Inférence évolutive :
   Les auteurs dérivent des algorithmes de filtre de Kalman et de lisseur Rauch–Tung–Striebel (RTS) modifiés pour ce système augmenté.

   • Complexité séquentielle : Les étapes de filtrage et de lissage opèrent sur la dimension de l'état augmenté ($d + N_{inst}$), maintenant une complexité computationnelle de $O(N)$ sur un CPU.
   • Complexité parallèle : En tirant parti de l'associativité des distributions de filtrage et de lissage, la méthode est reformulée comme un problème de somme préfixe totale (all-prefix-sums). Cela permet une implémentation parallèle sur GPU, réduisant le temps d'exécution à $O(N/T + \log T)$, où $T$ est le nombre de travailleurs parallèles.

4. Implémentation logicielle :
   Les auteurs introduisent smolgp, un package Python/JAX open-source. Il offre une compatibilité directe avec la bibliothèque existante tinygp, prenant en charge à la fois la modélisation GP standard et la modélisation GP consciente de l'exposition. Il implémente les solveurs séquentiels et parallèles décrits ci-dessus.

Résultats clés

• Équivalence numérique : Les auteurs valident l'approche SSM augmentée contre l'implémentation GP dense complète (Luhn et al. 2026) pour les scénarios à instrument unique et multi-instruments (chevauchants). Les résultats pour les moyennes conditionnelles, les variances et les log-vraisemblances concordent jusqu'à la précision machine (différences $\sim 10^{-14}$ à $10^{-15}$), confirmant l'équivalence théorique.
• Benchmarks de performance :
  • Mémoire : L'approche SSM est nettement plus économe en mémoire que la méthode GP dense. Alors que les GP denses nécessitent une mémoire en $O(N^2)$ (dépassant 100 Go pour $N \sim 10^5$), la méthode SSM s'étend linéairement, rendant les grands ensembles de données réalisables sur du matériel standard.
  • Temps d'exécution : Pour les données intégrées, où aucune alternative quasi-séparable (QSM) n'existe, la méthode SSM offre une accélération transformative. Sur un CPU, le solveur SSM séquentiel est le plus rapide pour la prédiction et compétitif pour le calcul de vraisemblance. Sur un GPU, le solveur SSM parallèle atteint les temps de conditionnement les plus rapides, rivalisant avec les méthodes QSM pour les vraisemblances.
  • Prédiction : La méthode prend en charge une prédiction efficace à des points de test arbitraires en utilisant un algorithme spécifique (Algorithme 1) qui est parallélisable de manière triviale et léger en mémoire par rapport aux multiplications matricielles denses.

Portée et affirmations
L'article affirme que ce travail permet des comparaisons croisées massives multi-instruments où le chevauchement des expositions est omniprésent, un scénario précédemment limité par le coût computationnel. En maintenant une complexité en $O(N)$ (ou $O(\log N)$ sur GPU), la méthode débloque le potentiel pour :

• L'analyse d'ensembles de données avec $N > 10^5$ observations (par exemple, la surveillance solaire quotidienne provenant de plusieurs instruments tels que KPF, NEID, EXPRES et HARPS-N).
• L'emploi de modèles plus complexes et d'ensembles de données de plus haute dimension qui étaient auparavant intraitables en raison du goulot d'étranglement en $O(N^3)$ des GP intégrés.
• L'extension de la modélisation GP évolutive à des noyaux qui manquent de quasi-séparabilité (par exemple, les noyaux quasipériodiques) lorsque des mesures intégrées sont impliquées.

Les auteurs soulignent que le cadre n'est pas limité aux séries temporelles mais s'applique à toutes les données binnées avec une coordonnée indépendante triable, et qu'il prend en charge naturellement les observations dérivées (par exemple, les modèles FF') grâce à l'inclusion des dérivées d'état dans la formulation SSM. Le travail positionne smolgp comme un outil évolutif et flexible pour l'analyse du domaine temporel astronomique de nouvelle génération.