Refonte du krigeage et des classes associées

On se propose de simplifier les classes liées à l’estimation paramétrique de processus Gaussiens (GeneralLinearModelAlgorithm abbrégé en GLMAlgo) et le conditionnement de tels processus (KrigingAlgorithm) sur les points suivants:

  • La calibration des modèles de covariance n’est plus déclenchée par l’algorithme de krigeage: on doit faire un appel explicite à GLMAlgo. C’est en particulier pertinent en vue d’introduire un krigeage séquentiel.
  • On ajoute un constructeur de KrigingAlgorithm prenant un GLMResult en argument pour pouvoir extraire la covariance calibrée ainsi que les sous-produits de la calibration (matrice de covariance + facteur de Cholesky)
  • On sépare le cas où les observations sont traitées séquentiellement comme des obsrevations de dimension 1 du cas où elles sont traitées globalement, résultant en deux algorithmes KrigingAlgorithm et KrigingAlgorithmND
  • On externalise l’adaptation du modèle de covariance aux dimensions du problème, qui correspondent à (explications issues du code C++):
    // There are 4 cases:
    // 1) Both the input dimension and the dimension of the model match the dimensions of the problem, in which case the model is used as-is
    // 2) The input dimension of the model is 1 and different from the input dimension of the problem, and the dimension of both the model and the problem are 1. The actual model is a product of the given model.
    covarianceModel_ = ProductCovarianceModel(ProductCovarianceModel::CovarianceModelCollection(inputDimension, covarianceModel));
    // 3) The input dimension of the model and the problem match, but the dimension of the model is 1, different from the dimension of the problem. The actual model is a tensorization of the given model.
    covarianceModel_ = TensorizedCovarianceModel(TensorizedCovarianceModel::CovarianceModelCollection(dimension, covarianceModel));
    // 4) The input dimension of the model is 1 and different from the input dimension of the problem, and the dimension of the model is 1 and different from the dimension of the problem. The actual model is a tensorization of products of the given model.
    covarianceModel_ = TensorizedCovarianceModel(TensorizedCovarianceModel::CovarianceModelCollection(dimension, ProductCovarianceModel(ProductCovarianceModel::CovarianceModelCollection(inputDimension, covarianceModel))));
    // The other situations are invalid.
  • Ne plus utiliser de collections de Basis: on fournit une Basis de la bonne dimension
  • Ne plus tester le caractère centré de l’échantillon des observations (@sofianehaddad on faisait ça pour quelle raison?)
  • On se propose de changer les noms GeneralLinearModelAlgorithm->GaussianProcessFitter, KrigingAlgorithm->GaussianProcessRegression. Plus généralement, on souhaite remplacer les suffixes Algorithm et Factory par des termes plus parlant (par exemple les DistributionFactory->DistributionFitter)

Pour se débarasser de GLM::setCovModel

// 1) Both the input dimension and the dimension of the model match the dimensions of the problem, in which case the model is used as-is
  // 2) The input dimension of the model is 1 and different from the input dimension of the problem, and the dimension of both the model and the problem are 1. The actual model is a product of the given model.
  // 3) The input dimension of the model and the problem match, but the dimension of the model is 1, different from the dimension of the problem. The actual model is a tensorization of the given model.
  // 4) The input dimension of the model is 1 and different from the input dimension of the problem, and the dimension of the model is 1 and different from the dimension of the problem. The actual model is a tensorization of products of the given model.
  // The other situations are invalid.

  const Bool sameDimension = dimension == covarianceModel.getOutputDimension();
  const Bool unitModelDimension = covarianceModel.getOutputDimension() == 1;
  const Bool sameSpatialDimension = inputDimension == covarianceModel.getInputDimension();
  const Bool unitModelSpatialDimension = covarianceModel.getInputDimension() == 1;
  // Case 1
  if (sameSpatialDimension && sameDimension)
    covarianceModel_ = covarianceModel;
  // Case 2
  else if (unitModelSpatialDimension && sameDimension && unitModelDimension)
    covarianceModel_ = ProductCovarianceModel(ProductCovarianceModel::CovarianceModelCollection(inputDimension, covarianceModel));
  // Case 3
  else if (sameSpatialDimension && unitModelDimension)
    covarianceModel_ = TensorizedCovarianceModel(TensorizedCovarianceModel::CovarianceModelCollection(dimension, covarianceModel));
  // Case 4
  else if (unitModelSpatialDimension && unitModelDimension)
    covarianceModel_ = TensorizedCovarianceModel(TensorizedCovarianceModel::CovarianceModelCollection(dimension, ProductCovarianceModel(ProductCovarianceModel::CovarianceModelCollection(inputDimension, covarianceModel))));
cov = AbsoluteExponential(...)
cov1 = cov
cov2 = ProductCovarianceModel([cov] * inputDimension)
cov3 = TensorizedCovarianceModel([cov] * dimension)
cov4 = TensorizedCovarianceModel([ProductCovarianceModel([cov] * inputDimension)] * dimension)

Quelle sera la portée de GaussianProcessFitter? Du moins comment sera utilisée cette classe? Comment s’articule le krigeage séquentiel autour des classes existantes? Est il possible de partager des uc?
@regislebrun on teste le caractère centré des échantillons car on autorise (aujourd’hui) l’utilisateur à ne pas fournir de Basis. Dans cette situation on veut être sur que l’échantillon est “centré” plutôt que de mettre par défaut une base constante

le GPFitter produit une classe result qui est passée au constructeur de KrigingAlgo pour réutiliser le modèle de covariance calibré et éventuellement des résultats de calcul intermédiaire.
Le KrigingAlgo ne fait plus la calibration; le modèle passé en entrée est considéré comme calibré.

@sofianehaddad En fait il n’y a aucune raison de demander à ce que les données soient centrées lorsqu’on fait du krigeage dans le cas où la fonction moyenne est nulle. En effet, meme si le processus Gaussien sous-jacent est centré, la moyenne empirique d’un échantillon d’un tel processus ne sera certainement pas nulle. C’est la même chose que pour un échantillon de loi normale centrée…

Pour compléter cette discussion je pense qu’il faudrait garder à l’esprit comment doivent /peuvent s’imbriquer de nouveaux développements/fonctionnalités. Je pense notamment à :

  • Leave-One-Out. Aujourd’hui on se base essentiellement sur la maximisation de vraisemblance. Comment implémenter un estimateur LOO pour l’estimation des paramètres dans cette nouvelle architecture?
  • Comment faire en sorte de considérer un modèle de krigeage prenant en compte l’estimation du bruit blanc (kernel + Dirac) ? Je ne pense pas que cela doit être spécifiquement un problème de “paramètrage” de noyau.
  • Aussi si on veut ajouter des estimations (approximations) pour le cas de grandes données (outre H-Mat) comment ajouter cet estimateur pour qu’il mutualise au maximum avec l’existant ?

En somme garder à l’esprit les potentiels futurs travaux (outre ceux déjà identifiés)

RAPPEL! On crée deux nouvelles classes:

  • GaussianProcessFitter est la classe en charge de l’estimation paramétrique d’un processus Gaussien dont la tendance est décomposée sur une base fonctionnelle donnée et dont le modèle de covariance est un modèle paramétrique. L’estimation des coefficients de la tendance et des paramètres actifs du modèle de covariance est faite mar maximisation de la vraisemblance. La tendance est accessible dans la structure résultat GaussianProcessFitterResult à la fois via l’accès à ses coefficients dans la base et comme fonction. En particulier, on a toutes les informations pour instancier un GaussianProcess
  • GaussianProcessRegression est la classe rendant interpolant un processus Gaussien donné par sa fonction tendance, son modèle de covariance et les échantillons d’entrée/sortie en lesquels se fait l’intrpolation. Aucun paramètre n’est estimé, et la tendance n’a pas à être décrite comme combinaison de fonctions de base. Un constructeur particulier consomme un GaussianProcessFitterResult afin d’éviter la discrétisation et la factorisation du modèle de covariance.
  • Pour construire un méta-modèle en traitant les sorties une à une il convient d’effectuer une boucle du type:
allMarginalMetaModels = list()
for index in range(outputSample.getDimension()):
  algoFit = ot.GaussianProcessFitter(inputSample, outputSample.getMarginal(index), basis, covarianceModel)
  algoFit.run()
  algoRegression = ot.GaussianProcessRegression(algoFit.getResult())
  algoRegression.run() 
  allMarginalMetaModels.append(algoRegression.getResult().getMetaModel())

metaModel = ot.AggregatedFunction(allMarginalMetaModels)