Compactly-supported nonstationary kernels for computing exact Gaussian processes on big data

1. Problématique

Les Processus Gaussiens (GP) sont un cadre de référence pour la modélisation statistique, l'approximation de fonctions stochastiques et l'apprentissage automatique, offrant des mesures d'incertitude naturelles. Cependant, leur application aux données massives modernes (millions à milliards de points) se heurte à deux limitations majeures :

1. Hypothèse de stationnarité : Les implémentations traditionnelles utilisent des noyaux stationnaires (la covariance ne dépend que de la distance entre les points). Cela manque de flexibilité pour modéliser des phénomènes réels complexes où les dépendances varient spatialement ou temporellement (non-stationnarité).
2. Problème d'échelle (Scalabilité) : L'inférence exacte dans un GP nécessite l'inversion d'une matrice de covariance de taille $N \times N$, ce qui entraîne une complexité computationnelle de $O(N^3)$ et une consommation mémoire de $O(N^2)$. Cela rend les méthodes exactes inapplicables au-delà de quelques milliers de points.
3. Limites des approches existantes :
   • Les méthodes d'approximation (ex: Vecchia, processus prédictif, régression par noyaux) réduisent la complexité mais introduisent des biais, des choix subjectifs (ordre des points, sélection des voisins) et ne garantissent pas de bonnes propriétés théoriques pour les noyaux non stationnaires.
   • Les noyaux compacts (à support fini) existent pour créer des matrices creuses, mais ils sont généralement stationnaires et ne peuvent pas capturer des corrélations à longue distance non liées à la distance physique (ex: téléconnexions climatiques).

2. Méthodologie

Les auteurs proposent un cadre entièrement bayésien basé sur un noyau compactement supporté et non stationnaire capable de découvrir la structure de parcimonie (sparsity) des données.

A. Construction du Noyau

Le noyau proposé, noté $C_y(x, x')$, est défini comme le produit d'un noyau de base (core) et d'un noyau induisant la parcimonie :
$$C_y(x, x'; \theta) = C_{core}(x, x'; \theta_{core}) \times C_{sparse}(x, x'; \theta_{sparse})$$

• $C_{core}$ : Peut être stationnaire (ex: Matérn) ou non stationnaire. Il capture les structures de dépendance globales.
• $C_{sparse}$ : C'est l'innovation centrale. Il est construit à partir de sommes et produits de fonctions "bump" (fonctions en bosse lisses et à support compact).
  • La fonction bump $g(x)$ est non nulle uniquement si la distance à un centre $h$ est inférieure à un rayon $r$.
  • $C_{sparse}$ est une combinaison linéaire de produits de sommes de ces fonctions bump.
  • Cela permet de créer des corrélations non nulles entre des points éloignés (si ils partagent une fonction bump commune) et des corrélations exactement nulles au-delà d'un certain seuil, indépendamment de la distance euclidienne.

B. Modélisation Bayésienne et Régularisation

• Le modèle inclut une fonction de moyenne a priori, un bruit de mesure et le processus GP.
• Régularisation : Pour éviter le surapprentissage et gérer le grand nombre d'hyperparamètres (centres, rayons, amplitudes des fonctions bump), les auteurs utilisent des priors bayésiens.
  • Les amplitudes des fonctions bump sont modélisées comme des variables binaires (activées/désactivées) via une distribution de Bernoulli. Cela permet au modèle de "éteindre" automatiquement les fonctions bump inutiles, découvrant ainsi la structure de parcimonie optimale dans les données.
• Inférence : L'inférence est réalisée via des méthodes MCMC (Monte Carlo par Chaîne de Markov) pour estimer les hyperparamètres et les coefficients de la moyenne.

C. Optimisation Computationnelle (HPC)

Pour traiter des données massives (ex: 1 million de points) tout en restant exacts :

• Calcul de la matrice de covariance : Les calculs sont effectués par lots (batches) en parallèle sur des GPU. La matrice est stockée au format creux (seuls les éléments non nuls sont conservés).
• Résolution du système linéaire : Au lieu de la décomposition de Cholesky (trop coûteuse pour les matrices creuses non structurées), les auteurs utilisent l'algorithme MINRES (Minimum Residual), robuste aux valeurs propres mal conditionnées.
• Calcul du déterminant : Estimé via la méthode de Lanczos, permettant de calculer la vraisemblance marginale sans inversion complète de la matrice.

3. Résultats Clés

A. Données Synthétiques

• Détection de parcimonie et non-stationnarité : Sur des données synthétiques générées avec des structures de covariance complexes (creuses et non stationnaires), le noyau proposé surpasse les méthodes exactes (GP stationnaire) et les méthodes d'approximation (Vecchia, processus prédictif, régression par noyaux).
• Robustesse à la spécification : Même si le nombre de fonctions bump est surestimé, la régularisation bayésienne permet de désactiver les fonctions inutiles, produisant des résultats quasi identiques à ceux d'une configuration optimale.
• Performance : Le noyau atteint des scores d'erreur (RMSE) et de probabilité (CRPS) inférieurs, particulièrement lorsque le rapport signal/bruit est faible et que les données sont non stationnaires.

B. Application Réelle : Température Maximale Quotidienne

Les auteurs appliquent leur méthode à plus d'un million de mesures de température maximale quotidienne aux États-Unis (données GHCN-D, 2001-2005).

• Comparaison : Le modèle GP exact non stationnaire et creux est comparé à la méthode de référence L15 (moyenne pondérée par la distance inverse) et à d'autres approximations GP.
• Performance :
  • Réduction de l'erreur quadratique moyenne (RMSE) d'environ 12,6 % par rapport à la méthode L15.
  • Meilleure performance que les autres méthodes GP (Vecchia, processus prédictif) sur l'ensemble des 60 jours testés.
  • Capacité à capturer la variabilité topographique (ex: dans les montagnes du Nevada/Californie) là où les méthodes traditionnelles échouent.
• Interprétabilité : Le modèle révèle des structures de corrélation spatiale non stationnaires et des périodes d'autocorrélation temporelle variables selon les régions.

4. Contributions Principales

1. Noyau Hybride Innovant : Développement d'un noyau qui combine la flexibilité non stationnaire et la parcimonie exacte (support compact) via des fonctions bump, permettant de modéliser des corrélations à longue distance non liées à la distance physique.
2. GP Exact à Grande Échelle : Démonstration qu'il est possible d'effectuer une inférence exacte (sans approximations de vraisemblance) sur des jeux de données de plus d'un million de points en utilisant des ressources de calcul haute performance (HPC) et des algorithmes de résolution de systèmes creux.
3. Réduction de la Subjectivité : Contrairement aux méthodes d'approximation (comme Vecchia) qui nécessitent des choix subjectifs (ordre des points, nombre de voisins), la parcimonie est découverte directement par les données via l'apprentissage bayésien des paramètres du noyau.
4. Validation Scientifique : Application réussie dans les sciences de la Terre, surpassant les méthodes de pointe actuelles pour la production de données griddées (maillées) avec quantification d'incertitude.

5. Signification et Impact

Ce travail remet en question l'idée reçue selon laquelle les Processus Gaussiens ne peuvent pas être appliqués aux "Big Data". En démontrant que la densité de la matrice de covariance n'est pas une propriété fondamentale du GP mais le résultat d'un mauvais choix de noyau, les auteurs ouvrent la voie à l'utilisation de GP exacts pour des problèmes à grande échelle.

L'approche permet de concilier la rigueur statistique des GP (incertitude, inférence bayésienne) avec la puissance des méthodes d'apprentissage automatique modernes. Elle est particulièrement pertinente pour les sciences environnementales et climatiques où les données sont massives, irrégulièrement espacées et présentent des structures complexes non stationnaires. L'article suggère également que des échelles encore plus grandes (jusqu'à 100 millions de points) sont envisageables avec l'augmentation des capacités de calcul parallèle.