Vecchia approximated Bayesian heteroskedastic Gaussian processes

1. Problématique

Les simulations informatiques modernes sont de plus en plus stochastiques (aléatoires), produisant un bruit dépendant des entrées (hétéroscédasticité). Dans de tels cas, les processus gaussiens (GP) standards, qui supposent un bruit constant (homoscédatique), sont inadéquats. Les processus gaussiens hétéroscédastiques (hetGPs) sont la solution idéale car ils modélisent à la fois la moyenne et la variance latente non constante.

Cependant, les implémentations existantes de hetGP souffrent de deux limitations majeures :

1. Incapacité à gérer les grandes campagnes de simulation : Les méthodes actuelles nécessitent la décomposition d'une matrice de covariance dense de taille $N \times N$ (où $N$ est le nombre total de simulations), entraînant un coût computationnel en $O(N^3)$. Cela limite l'applicabilité à quelques milliers de points.
2. Estimation par points (Point-estimates) : La plupart des méthodes (comme hetGP de Binois et al., 2018) utilisent une estimation par maximum de vraisemblance (MLE) ou EM pour les variances latentes, au lieu d'une intégration bayésienne complète. Cela sous-estime l'incertitude (UQ) et limite la fidélité inférentielle.

L'article se concentre sur un cas d'étude concret : la prévision de la température des lacs via le modèle NOAA-GLM, générant plus de 9 millions de simulations avec des niveaux de bruit variables selon la profondeur et le jour de l'année.

2. Méthodologie

Les auteurs proposent une nouvelle approche, nommée bhetGP, qui combine quatre ingrédients clés pour rendre l'inférence bayésienne hétéroscédastique à grande échelle possible :

A. Modélisation Hiérarchique Couplée

Le modèle repose sur deux processus gaussiens couplés :

• Un GP pour la moyenne des sorties.
• Un GP pour le logarithme des variances latentes ($\log \Lambda$).
  Contrairement aux approches antérieures qui découplent ces processus pour des raisons de calcul, bhetGP les maintient liés pour une inférence conjointe.

B. Inférence par Échantillonnage Elliptique (ESS)

Au lieu d'utiliser des méthodes d'optimisation (MLE/EM) ou des chaînes de Markov (MCMC) traditionnelles avec des propositions de marche aléatoire (qui convergent mal en haute dimension), les auteurs utilisent l'Échantillonnage Elliptique (Elliptical Slice Sampling - ESS).

• L'ESS est particulièrement efficace pour les priors gaussiens multivariés.
• Il permet un échantillonnage conjoint de toutes les variances latentes sans rejet, évitant les chaînes « collantes » (sticky chains) et assurant une meilleure exploration de l'espace des paramètres.

C. Statistiques Suffisantes et Identité de Woodbury

Pour gérer les simulations avec réplications (plusieurs exécutions pour une même entrée unique), l'article utilise l'approche de Binois et al. (2018).

• Au lieu de travailler sur les $N$ observations brutes, le modèle travaille sur les $n$ entrées uniques ($n \ll N$).
• L'identité de Woodbury est utilisée pour reformuler la vraisemblance et les équations de krigeage, réduisant la complexité de $O(N^3)$ à $O(n^3)$.
• Cela permet d'utiliser des statistiques suffisantes (moyennes et variances empiriques des réplications) plutôt que les données brutes.

D. Approximation de Vecchia

Pour surmonter le goulot d'étranglement lorsque le nombre d'entrées uniques $n$ est très grand (ex: $n > 300\,000$), les auteurs intègrent l'approximation de Vecchia.

• Cette méthode factorise la densité de probabilité conjointe en un produit de conditionnelles, en ne considérant qu'un petit nombre de voisins ($m$) pour chaque point.
• Elle induit une structure creuse (sparse) dans la décomposition de Cholesky de la matrice de précision, réduisant la complexité à $O(nm^3)$.
• Innovation clé : L'approximation de Vecchia est appliquée directement sur la vraisemblance de Woodbury (basée sur les $n$ entrées uniques), et non sur les $N$ données brutes. Cela permet de capturer correctement la structure de réplication tout en restant computationnellement efficace.

3. Contributions Clés

1. Première intégration bayésienne complète pour hetGP à grande échelle : L'article est le premier à combiner l'ESS, l'identité de Woodbury et l'approximation de Vecchia pour fournir une inférence postérieure complète (et non par points) sur des jeux de données massifs.
2. Algorithme hybride (Woodbury + Vecchia) : Les auteurs montrent comment insérer l'approximation de Vecchia dans le cadre de vraisemblance de Woodbury, ce qui est crucial pour gérer les réplications sans perdre en efficacité.
3. Implémentation Open Source : Le package R bhetGP est disponible sur CRAN, permettant la reproduction des résultats et l'application à d'autres problèmes.
4. Gestion de l'incertitude (UQ) : En passant d'une estimation par point à une intégration bayésienne, le modèle fournit une quantification de l'incertitude beaucoup plus fiable, essentielle pour la prise de décision.

4. Résultats

Les performances de bhetGP ont été évaluées sur deux types de données :

• Données de benchmark (Assemble-to-Order) :

  • Sur des problèmes de simulation stochastique classiques, bhetGP surpasse les méthodes hetGP standards (MLE) en termes de précision (RMSE) et de score de vraisemblance (Proper Score).
  • Contrairement aux méthodes classiques, bhetGP peut utiliser des ensembles de données d'entraînement massifs (jusqu'à 11 000 points uniques) pour améliorer la précision, là où les méthodes classiques échouent.

• Données NOAA-GLM (Température des lacs) :

  • Sur un corpus de 9 millions de simulations ($N \approx 9M$, $n \approx 290\,000$), bhetGP est capable de s'ajuster et de prédire, là où les méthodes hetGP standards sont inapplicables.
  • Comparé à l'approche précédente de Holthuijzen et al. (2024) qui utilisait une approximation de Vecchia mais découplait moyenne et variance (approche SK-like), bhetGP offre :
    • Une précision supérieure (RMSE plus faible).
    • Des intervalles de prédiction plus étroits tout en maintenant une couverture nominale (90%).
    • Une meilleure calibration de l'incertitude, évitant les sous-estimations fréquentes des méthodes par point.

5. Signification et Impact

Cet article représente une avancée significative dans le domaine de la modélisation par substituts (surrogates) pour les simulations stochastiques massives.

• Passage à l'échelle (Scalability) : Il démontre qu'il est possible de réaliser une inférence bayésienne complète (et non approximative par optimisation) sur des millions de points de données, ce qui était auparavant impossible.
• Fiabilité de l'incertitude : En fournissant une distribution postérieure complète pour les variances latentes, le modèle évite les pièges des estimations ponctuelles, offrant des intervalles de confiance plus robustes pour la gestion des risques (ex: qualité de l'eau, prévisions écologiques).
• Futur : Les auteurs prévoient d'utiliser bhetGP pour la calibration de modèles (ajustement des simulations aux données de terrain) et l'apprentissage actif, des tâches où une quantification précise de l'incertitude est critique.

En résumé, bhetGP comble le fossé entre la rigueur théorique de l'inférence bayésienne hétéroscédastique et les contraintes pratiques des simulations informatiques modernes à très grande échelle.