Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?

Published: 10 décembre 2025
DMCA.com Protection Status

Published: 10 décembre 2025

DMCA.com Protection Status

Introduction

Les Processus Gaussiens sont de puissants modèles probabilistes qui représentent des distributions sur des fonctions.
Ils excellent en régression et en classification, surtout lorsque la quantification de l’incertitude et la performance avec peu de données sont importantes.

Table des matières

Leurs principaux inconvénients sont le coût computationnel et la dépendance à un bon choix de noyau (kernel).

Qu’est-ce qu’un Processus Gaussien ? (GP)

Un Processus Gaussien (GP) est un modèle probabiliste qui définit une distribution sur des fonctions.
Au lieu de supposer une forme fonctionnelle spécifique (comme un polynôme ou un réseau de neurones), un GP suppose que :

Tout ensemble fini de valeurs de la fonction suit une distribution gaussienne multivariée.

Un GP est entièrement défini par :

  • une fonction de moyenne ( m(x) )
  • une fonction de covariance ou kernel ( k(x, x') )

\begin{equation}
[
f(x) \sim \mathrm{GP}(m(x), k(x,x'))
]
\end{equation}

Le kernel encode la lissité, la similarité, la périodicité, le bruit, ainsi que d’autres hypothèses structurelles sur la fonction.

Régression par Processus Gaussien (GPR)

La GPR modélise une relation :

\begin{equation}
y = f(x) + \epsilon, \quad \epsilon \sim \mathrm{N}(0, \sigma^2)
\end{equation}

À partir des données d’entraînement ((X, y)), le GP infère une distribution a posteriori sur les fonctions.
Cela fournit :

  • une moyenne prédictive (meilleure estimation)
  • une variance prédictive (incertitude)

Ainsi, la GPR n’est pas seulement une méthode de prédiction : c’est un cadre complet de quantification de l’incertitude.

Pourquoi la GPR est efficace pour la régression ?

  • Fonctionne bien avec des jeux de données petits ou moyens
  • Fournit de l’incertitude (intervalles de confiance)
  • Équilibre automatiquement biais et complexité via le kernel

Classification par Processus Gaussien (GPC)

Pour la classification, les sorties ne sont pas gaussiennes. Les GPs utilisent donc une fonction de lien (logistique ou probit) :

\begin{equation}
[
p(y=1 \mid f) = \sigma(f)
]
\end{equation}

Le GP définit toujours un prior sur ( f(x) ), mais la vraisemblance est non gaussienne.
L’inférence nécessite des approximations :

  • Approximation de Laplace
  • Propagation d’attente (Expectation Propagation)
  • Inférence variationnelle

La GPC fournit des prédictions probabilistes, permettant d’évaluer la confiance dans la classification et des probabilités bien calibrées.

Avantages des Processus Gaussiens

  1. Quantification d’incertitude intégrée
    Les GPs retournent une distribution complète (moyenne + variance) et pas seulement une estimation ponctuelle.

  2. Excellente performance avec peu ou du bruit dans les données
    Idéal lorsque les données sont coûteuses ou limitées.

  3. Grande flexibilité grâce aux kernels
    Permet de modéliser :

  4. fonctions lisses

  5. motifs périodiques
  6. corrélations spatiales

  7. structures complexes via somme/produit de kernels

  8. Peu d’hyperparamètres
    L’entraînement ne concerne généralement que les paramètres du kernel, souvent moins nombreux que dans les réseaux de neurones.

  9. Non-paramétrique
    La complexité du modèle croît avec les données, sans avoir à spécifier une forme fonctionnelle fixe.

Limitations des Processus Gaussiens

  1. Mauvaise scalabilité avec de grandes données
    Coût d’entraînement :

\begin{equation}
[
O(n^3) \quad \text{en temps}, \qquad O(n^2) \quad \text{en mémoire}
]
\end{equation}

Peu adaptés aux grandes bases de données sauf si l’on utilise des approximations (sparse GPs, points induits…).

  1. Choix du kernel complexe
    La performance dépend fortement du type de kernel et de son optimisation.

  2. La GPC nécessite des approximations
    La classification est plus coûteuse computationnellement à cause des vraisemblances non gaussiennes.

  3. Vitesse de prédiction limitée
    La prédiction dépend du nombre de points d’entraînement ; les grandes bases ralentissent l’inférence.

## Régression par Processus Gaussien (GPR)

### Bruit homoscédastique – Exemple 1

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
import torch
 import torch.nn as nn
 import numpy as np
 import matplotlib.pyplot as plt

 # Utiliser la double précision pour la stabilité numérique avec l'algèbre linéaire
 dtype = torch.double
 torch.manual_seed(0)
 np.random.seed(0)

 # 1) Génération de données synthétiques
 def generate_data(n_train=25, noise_std=0.2):
     X = np.linspace(-3, 3, n_train)[:, None]
     y = np.sin(X).ravel() + noise_std * np.random.randn(n_train)
     return X, y

 X_train_np, y_train_np = generate_data(n_train=25, noise_std=0.2)
 X_test_np = np.linspace(-5, 5, 400)[:, None]

 X_train = torch.tensor(X_train_np, dtype=dtype)
 y_train = torch.tensor(y_train_np[:, None], dtype=dtype)  # vecteur colonne
 X_test = torch.tensor(X_test_np, dtype=dtype)

 # 2) Kernel RBF implémenté en PyTorch (paramétré en log-space)
 class RBFKernel(nn.Module):
     def __init__(self, init_lengthscale=1.0, init_variance=1.0):
         super().__init__()
         self.log_lengthscale = nn.Parameter(torch.tensor(np.log(init_lengthscale), dtype=dtype))
         self.log_variance = nn.Parameter(torch.tensor(np.log(init_variance), dtype=dtype))

     def forward(self, x1, x2):
         # Matrice des distances au carré
         x1_sq = (x1**2).sum(dim=1, keepdim=True)
         x2_sq = (x2**2).sum(dim=1, keepdim=True)
         sqdist = x1_sq + x2_sq.T - 2.0 * (x1 @ x2.T)
         lengthscale = torch.exp(self.log_lengthscale)
         variance = torch.exp(self.log_variance)
         K = variance * torch.exp(-0.5 * sqdist / (lengthscale**2))
         return K, lengthscale, variance

 # 3) Classe Régression par Processus Gaussien (NLML et prédiction)
 class GaussianProcessRegression(nn.Module):
     def __init__(self, kernel: RBFKernel, noise_init=0.2):
         super().__init__()
         self.kernel = kernel
         self.log_noise = nn.Parameter(torch.tensor(np.log(noise_init), dtype=dtype))

     def negative_log_marginal_likelihood(self, X, y, jitter=1e-6):
         K, _, _ = self.kernel(X, X)
         noise = torch.exp(self.log_noise)
         n = X.shape[0]
         Kt = K + (noise**2 + jitter) * torch.eye(n, dtype=dtype)

         # Décomposition de Cholesky
         L = torch.linalg.cholesky(Kt)
         # Résolution pour alpha = K^{-1} y via systèmes triangulaires
         v = torch.linalg.solve_triangular(L, y, upper=False)
         alpha = torch.linalg.solve_triangular(L.T, v, upper=True)

         data_fit = 0.5 * (y.T @ alpha)
         log_det = torch.sum(torch.log(torch.diag(L)))
         constant = 0.5 * n * torch.log(torch.tensor(2.0 * np.pi, dtype=dtype))
         nlml = data_fit + log_det + constant
         return nlml.squeeze()

     def predict(self, X_train, y_train, X_test, jitter=1e-6):
         K_tt, _, _ = self.kernel(X_train, X_train)
         noise = torch.exp(self.log_noise)
         n = X_train.shape[0]
         Kt = K_tt + (noise**2 + jitter) * torch.eye(n, dtype=dtype)
         L = torch.linalg.cholesky(Kt)

         K_s, _, _ = self.kernel(X_train, X_test)
         K_ss, _, _ = self.kernel(X_test, X_test)

         # moyenne prédictive
         v = torch.linalg.solve_triangular(L, y_train, upper=False)
         alpha = torch.linalg.solve_triangular(L.T, v, upper=True)
         pred_mean = K_s.T @ alpha

         # covariance prédictive via complément de Schur
         w = torch.linalg.solve_triangular(L, K_s, upper=False)
         z = torch.linalg.solve_triangular(L.T, w, upper=True)
         pred_cov = K_ss - K_s.T @ z
         pred_var = torch.clamp(torch.diag(pred_cov), min=1e-12).unsqueeze(1)
         return pred_mean.detach(), pred_var.detach()

 # 4) Instanciation et entraînement
 kernel = RBFKernel(init_lengthscale=1.0, init_variance=1.0)
 gpr = GaussianProcessRegression(kernel=kernel, noise_init=0.2).double()

 optimizer = torch.optim.Adam(gpr.parameters(), lr=0.05)
 n_iters = 300  # réduire pour des exécutions rapides ; augmenter pour meilleure convergence
 loss_history = []

 for i in range(n_iters):
     optimizer.zero_grad()
     nlml = gpr.negative_log_marginal_likelihood(X_train, y_train)
     nlml.backward()
     optimizer.step()
     loss_history.append(nlml.item())
     if (i+1) % 50 == 0 or i==0:
         with torch.no_grad():
             ls = torch.exp(gpr.kernel.log_lengthscale).item()
             var = torch.exp(gpr.kernel.log_variance).item()
             noise = torch.exp(gpr.log_noise).item()
         print(f"Itération {i+1}/{n_iters} - NLML: {nlml.item():.4f}  lengthscale: {ls:.4f}  variance: {var:.4f}  noise: {noise:.4f}")

 # 5) Prédiction et visualisation
 pred_mean, pred_var = gpr.predict(X_train, y_train, X_test)
 pred_mean = pred_mean.numpy().ravel()
 pred_std = np.sqrt(pred_var.numpy().ravel())

 plt.figure(figsize=(8,5))
 plt.scatter(X_train_np.ravel(), y_train_np.ravel(), marker='o', label='Données d’entraînement')
 plt.plot(X_test_np.ravel(), pred_mean, label='Moyenne prédictive')
 upper = pred_mean + 1.96 * pred_std
 lower = pred_mean - 1.96 * pred_std
 plt.fill_between(X_test_np.ravel(), lower, upper, alpha=0.3)
 plt.title("Régression par Processus Gaussien (PyTorch)")
 plt.xlabel("x")
 plt.ylabel("y")
 plt.legend()
 plt.show()

 with torch.no_grad():
     final_ls = torch.exp(gpr.kernel.log_lengthscale).item()
     final_var = torch.exp(gpr.kernel.log_variance).item()
     final_noise = torch.exp(gpr.log_noise).item()
 print(f"Hyperparamètres finaux:\n lengthscale = {final_ls:.4f}\n variance = {final_var:.4f}\n noise = {final_noise:.4f}")

!

Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?
Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?

### Explication du code de Régression par Processus Gaussien

Ce code implémente la régression par processus gaussien (GPR) à partir de zéro en PyTorch, incluant :

  • définition du kernel
  • optimisation de la vraisemblance marginale
  • prédiction
  • visualisation

1. Imports et configuration

On importe :

  • PyTorch pour la différentiation automatique et l’algèbre linéaire
  • NumPy pour générer les données synthétiques
  • Matplotlib pour tracer les graphiques

On utilise :

1
2
3
dtype = torch.double
 torch.manual_seed(0)
 np.random.seed(0)

La double précision est recommandée pour les Processus Gaussiens à cause de la décomposition de Cholesky et des inversions de matrices. Les graines (seed) assurent la reproductibilité des résultats.

2. Génération de données synthétiques

1
2
3
4
def generate_data(n_train=25, noise_std=0.2):
     X = np.linspace(-3, 3, n_train)[:, None]
     y = np.sin(X).ravel() + noise_std * np.random.randn(n_train)
     return X, y
  • Les entrées (X) sont des points de –3 à 3.
  • Les cibles (y) suivent une sinusoïde bruitée.
  • Les données de test couvrent un intervalle plus large pour la visualisation.

Les tableaux sont convertis en tenseurs PyTorch en double précision.

3. Classe Kernel RBF

Le kernel RBF (Radial Basis Function) est une fonction de covariance classique :

\begin{equation}
[
k(x,x') = \sigma^2 \exp\left( -\frac{|x-x'|^2}{2\ell^2} \right)
]
\end{equation}

  • Les paramètres log_lengthscale et log_variance sont appris en log-space pour la stabilité numérique.
  • La méthode forward calcule les distances au carré entre tous les points et construit la matrice du kernel (K).

4. Classe Régression par Processus Gaussien

Cette classe gère :

(a) Log-vraisemblance marginale négative (NLML)

\begin{equation}
[
\log p(y|X,\theta) = -\frac12 y^\top K^{-1} y - \frac12 \log |K| - \frac{n}{2}\log(2\pi)
]
\end{equation}

Implémentée avec :

  • calcul de la matrice du kernel
  • ajout de la variance du bruit
  • décomposition de Cholesky
  • résolution de systèmes triangulaires pour calculer (\alpha = K^{-1} y)

La NLML est utilisée comme fonction de perte pour l’entraînement.

(b) Fonction de prédiction

Pour de nouveaux points (X_*) :

\begin{equation}
\mu_* = K_{*}^{T} K^{-1} y
\end{equation}

\begin{equation}
\Sigma_{*} = K_{**} - K_{*}^{T} K^{-1} K_*
\end{equation}

Le code calcule :

  • la moyenne prédictive (pred_mean)
  • la variance prédictive (pred_var)

Les tenseurs sont détachés pour ne pas suivre les gradients.

5. Instanciation et entraînement du modèle

1
2
3
kernel = RBFKernel(...)
 gpr = GaussianProcessRegression(...)
 optimizer = torch.optim.Adam(gpr.parameters(), lr=0.05)

Le modèle apprend :

  • le lengthscale
  • la variance
  • la variance du bruit

en minimisant la NLML.

La boucle d’entraînement :

  • calcule la NLML
  • rétropropagation des gradients
  • mise à jour des paramètres via Adam
  • affichage de l’avancement tous les 50 itérations

6. Prédiction et visualisation

Après entraînement :

1
2
pred_mean, pred_var = gpr.predict(...)
 pred_std = np.sqrt(pred_var...)

Le graphique affiche :

  • les points de données d’entraînement
  • la moyenne a posteriori du GP (courbe lisse)
  • l’intervalle de confiance 95 % (zone ombrée)

Comment changer le bruit d’observation

Le bruit d’observation est entièrement contrôlé par cette ligne à l’intérieur de la classe GaussianProcessRegression :

1
self.log_noise = nn.Parameter(torch.tensor(np.log(noise_init), dtype=dtype))

et utilisé ici pendant l’entraînement :

1
2
noise = torch.exp(self.log_noise)
 Kt = K + (noise**2 + jitter) * torch.eye(n, dtype=dtype)

Cela signifie que le modèle suppose une seule valeur globale de bruit, qui est apprise pendant l’entraînement.

OPTION 1 — Changer la valeur initiale du bruit

Modifier :

1
gpr = GaussianProcessRegression(kernel=kernel, noise_init=0.2)

Par exemple :

1
gpr = GaussianProcessRegression(kernel=kernel, noise_init=0.8)

Mais rappelez-vous : le GP apprendra toujours un seul niveau global de bruit, quelle que soit la valeur initiale.

OPTION 2 — Geler le bruit et le fixer manuellement

Si vous voulez un niveau de bruit fixe (non appris), faites :

1
2
gpr.log_noise.requires_grad_(False)
 gpr.log_noise.data = torch.log(torch.tensor(0.5, dtype=dtype))

Cela force le GP à utiliser :

\begin{equation}
[
\sigma_n = 0.5
]
\end{equation}

pour toutes les observations.

Bruit homoscédastique – Exemple 2

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
import torch
import numpy as np
import matplotlib.pyplot as plt
torch.set_default_dtype(torch.double)

# 1) Données
X = torch.linspace(-3, 3, 25).unsqueeze(1)
y = torch.sin(X) + 0.2 * torch.randn_like(X)

# incertitude des mesures (exemple)
y_err = 0.5 * torch.ones_like(y)

X_test = torch.linspace(-5, 5, 400).unsqueeze(1)

# 2) Paramètres (log-space)
log_length = torch.tensor(0.0, requires_grad=True)
log_var    = torch.tensor(0.0, requires_grad=True)
log_noise  = torch.tensor(-1.0, requires_grad=True)

# 3) Kernel RBF
def rbf(x1, x2, log_l, log_v):
    l = torch.exp(log_l)
    v = torch.exp(log_v)
    sqdist = (x1**2).sum(1, keepdim=True) + (x2**2).sum(1) - 2*x1@x2.T
    return v * torch.exp(-0.5 * sqdist / l**2)

# 4) Log-vraisemblance marginale négative
def nlml():
    K = rbf(X, X, log_length, log_var)
    noise = torch.exp(log_noise)
    K = K + noise**2 * torch.eye(len(X))
    L = torch.linalg.cholesky(K)
    alpha = torch.cholesky_solve(y, L)
    return (
        0.5 * (y.T @ alpha) +
        torch.sum(torch.log(torch.diag(L))) +
        0.5 * len(X) * torch.log(torch.tensor(2*np.pi))
    )

# 5) Entraînement des hyperparamètres
opt = torch.optim.Adam([log_length, log_var, log_noise], lr=0.05)
for i in range(300):
    opt.zero_grad()
    loss = nlml()
    loss.backward()
    opt.step()
    if i % 50 == 0: print(i, loss.item())

# 6) Distribution prédictive
def predict(X_new):
    K = rbf(X, X, log_length, log_var)
    noise = torch.exp(log_noise)
    K = K + torch.diag(noise**2 + y_err.squeeze()**2)

    L = torch.linalg.cholesky(K)
    Ks = rbf(X, X_new, log_length, log_var)
    Kss = rbf(X_new, X_new, log_length, log_var)

    alpha = torch.cholesky_solve(y, L)
    mean = Ks.T @ alpha

    v = torch.linalg.solve(L, Ks)
    cov = Kss - v.T @ v
    std = torch.diag(cov).clamp(min=1e-9).sqrt()

    return mean.squeeze().detach(), std.detach()

mu, std = predict(X_test)

# 7) Visualisation
plt.figure(figsize=(8,5))

# Données d'entraînement AVEC barres d'erreur
plt.errorbar(
    X.squeeze().numpy(),
    y.squeeze().numpy(),
    yerr=y_err.squeeze().numpy(),
    fmt='o',
    color='black',
    ecolor='gray',
    elinewidth=1.5,
    capsize=3,
    label="Données d'entraînement ± erreur"
)

# Moyenne prédictive
plt.plot(X_test.numpy(), mu.numpy(), label="Moyenne prédictive")

# Région prédictive ombrée
plt.fill_between(
    X_test.squeeze().numpy(),
    (mu - 2*std).numpy(),
    (mu + 2*std).numpy(),
    alpha=0.3,
    label="GP ±2σ"
)

plt.legend()
plt.title("Régression par Processus Gaussien avec barres d'erreur")
plt.show()

Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?
Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?

Bruit hétéroscédastique – Exemple 3

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt

# Utiliser la double précision pour la stabilité numérique
dtype = torch.double
torch.manual_seed(0)
np.random.seed(0)

# -------------------------------------------------------
# 1) Générer des données synthétiques + erreurs hétéroscédastiques
# -------------------------------------------------------
def generate_data(n_train=25, noise_std=0.2):
    X = np.linspace(-3, 3, n_train)[:, None]
    y = np.sin(X).ravel() + noise_std * np.random.randn(n_train)
    return X, y

X_train_np, y_train_np = generate_data(n_train=25, noise_std=0.2)
X_test_np = np.linspace(-5, 5, 400)[:, None]

# Exemple : incertitude hétéroscédastique
# On augmente l'erreur vers les bords pour démonstration
obs_error_std_np = 0.1 + 0.3 * np.abs(X_train_np.ravel()) / 3.0

X_train = torch.tensor(X_train_np, dtype=dtype)
y_train = torch.tensor(y_train_np[:, None], dtype=dtype)
obs_error_std = torch.tensor(obs_error_std_np[:, None], dtype=dtype)

X_test = torch.tensor(X_test_np, dtype=dtype)

# -------------------------------------------------------
# 2) Kernel RBF
# -------------------------------------------------------
class RBFKernel(nn.Module):
    def __init__(self, init_lengthscale=1.0, init_variance=1.0):
        super().__init__()
        self.log_lengthscale = nn.Parameter(torch.tensor(np.log(init_lengthscale), dtype=dtype))
        self.log_variance = nn.Parameter(torch.tensor(np.log(init_variance), dtype=dtype))

    def forward(self, x1, x2):
        x1_sq = (x1**2).sum(dim=1, keepdim=True)
        x2_sq = (x2**2).sum(dim=1, keepdim=True)
        sqdist = x1_sq + x2_sq.T - 2.0 * (x1 @ x2.T)
        lengthscale = torch.exp(self.log_lengthscale)
        variance = torch.exp(self.log_variance)
        K = variance * torch.exp(-0.5 * sqdist / (lengthscale**2))
        return K, lengthscale, variance

# -------------------------------------------------------
# 3) Régression par Processus Gaussien avec erreurs par observation
# -------------------------------------------------------
class GaussianProcessRegression(nn.Module):
    def __init__(self, kernel: RBFKernel):
        super().__init__()
        self.kernel = kernel

    def negative_log_marginal_likelihood(self, X, y, obs_error_std, jitter=1e-6):
        K, _, _ = self.kernel(X, X)
        n = X.shape[0]

        # --- NOUVEAU : utiliser les variances par observation ---
        noise_diag = obs_error_std.squeeze()**2
        Kt = K + torch.diag(noise_diag + jitter)

        # Décomposition de Cholesky
        L = torch.linalg.cholesky(Kt)

        v = torch.linalg.solve_triangular(L, y, upper=False)
        alpha = torch.linalg.solve_triangular(L.T, v, upper=True)

        data_fit = 0.5 * (y.T @ alpha)
        log_det = torch.sum(torch.log(torch.diag(L)))
        constant = 0.5 * n * torch.log(torch.tensor(2.0 * np.pi, dtype=dtype))
        nlml = data_fit + log_det + constant
        return nlml.squeeze()

    def predict(self, X_train, y_train, X_test, obs_error_std, jitter=1e-6):
        K_tt, _, _ = self.kernel(X_train, X_train)

        noise_diag = obs_error_std.squeeze()**2
        Kt = K_tt + torch.diag(noise_diag + jitter)

        L = torch.linalg.cholesky(Kt)

        K_s, _, _ = self.kernel(X_train, X_test)
        K_ss, _, _ = self.kernel(X_test, X_test)

        v = torch.linalg.solve_triangular(L, y_train, upper=False)
        alpha = torch.linalg.solve_triangular(L.T, v, upper=True)
        pred_mean = K_s.T @ alpha

        w = torch.linalg.solve_triangular(L, K_s, upper=False)
        z = torch.linalg.solve_triangular(L.T, w, upper=True)
        pred_cov = K_ss - K_s.T @ z

        pred_var = torch.clamp(torch.diag(pred_cov), min=1e-12).unsqueeze(1)
        return pred_mean.detach(), pred_var.detach()

# -------------------------------------------------------
# 4) Instanciation et entraînement
# -------------------------------------------------------
kernel = RBFKernel(init_lengthscale=1.0, init_variance=1.0)
gpr = GaussianProcessRegression(kernel=kernel).double()

optimizer = torch.optim.Adam(gpr.parameters(), lr=0.05)
n_iters = 300
loss_history = []

for i in range(n_iters):
    optimizer.zero_grad()
    nlml = gpr.negative_log_marginal_likelihood(X_train, y_train, obs_error_std)
    nlml.backward()
    optimizer.step()
    loss_history.append(nlml.item())
    if (i+1) % 50 == 0:
        with torch.no_grad():
            ls = torch.exp(gpr.kernel.log_lengthscale).item()
            var = torch.exp(gpr.kernel.log_variance).item()
        print(f"Iter {i+1}/{n_iters} | NLML={nlml:.4f} | ls={ls:.4f} | var={var:.4f}")

# -------------------------------------------------------
# 5) Prédiction et visualisation
# -------------------------------------------------------
pred_mean, pred_var = gpr.predict(X_train, y_train, X_test, obs_error_std)
pred_mean = pred_mean.numpy().ravel()
pred_std = np.sqrt(pred_var.numpy().ravel())

plt.figure(figsize=(10,5))
plt.scatter(X_train_np.ravel(), y_train_np.ravel(), c='k', label='Train', s=40)
plt.errorbar(
    X_train_np.ravel(), 
    y_train_np.ravel(),
    yerr=obs_error_std_np,
    fmt='none',
    ecolor='gray',
    alpha=0.5,
    label='Erreurs d’observation'
)
plt.plot(X_test_np.ravel(), pred_mean, label="Moyenne")
plt.fill_between(
    X_test_np.ravel(),
    pred_mean - 1.96*pred_std,
    pred_mean + 1.96*pred_std,
    alpha=0.3,
    label="IC 95%"
)
plt.legend()
plt.title("Régression par Processus Gaussien avec erreurs par observation")
plt.show()

Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?
Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?

Processus Gaussiens avec GPyTorch

Qu’est-ce que GPyTorch ?

GPyTorch est une bibliothèque moderne et scalable pour les processus gaussiens (GP), construite sur PyTorch.
Elle fournit une accélération GPU, la différentiation automatique, et des méthodes efficaces d’algèbre linéaire, permettant d’entraîner des GP sur des ensembles de données beaucoup plus grands que les implémentations classiques.

Elle a été développée chez Uber AI Labs et supporte :

  • Les processus gaussiens exacts (pour petits/moyens ensembles de données)
  • Les GP approximatifs / variationnels (pour grands ensembles de données)
  • Les kernels additifs, les GP multitâches, les kernels de mélange spectral
  • Une intégration naturelle avec les modèles PyTorch

Comment installer GPyTorch

Installation avec pip :

1
pip install gpytorch

Ou avec conda :

1
conda install -c conda-forge gpytorch

Vous devez avoir PyTorch installé au préalable :

1
pip install torch

Pourquoi GPyTorch est intéressant

Différentiation automatique
Les paramètres du kernel, le bruit et les paramètres du modèle sont optimisés automatiquement grâce à l’autograd de PyTorch.

Accélération GPU
Les décompositions de Cholesky, la résolution de matrices et les opérations sur les kernels s’exécutent sur GPU.

Kernels intégrés
RBF, Matern, Spectral Mixture, RQ, Polynômes, etc.

Support pour GP à grande échelle
Via les kernels SKI, les GP variationnels et l’interpolation de kernels structurée.

API propre et modulaire
Les modèles sont des sous-classes de gpytorch.models.ExactGP.

Support du bruit hétéroscédastique
Via likelihood(noise) ou des modèles de bruit personnalisés.

Implémentation complète de GPyTorch (avec bruit hétéroscédastique)

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
import torch
import gpytorch
import numpy as np
import matplotlib.pyplot as plt

dtype = torch.double
torch.manual_seed(0)
np.random.seed(0)

# -------------------------------------------------------
# 1) Génération des données
# -------------------------------------------------------
def generate_data(n_train=25, noise_std=0.2):
    X = np.linspace(-3, 3, n_train)[:, None]
    y = np.sin(X).ravel() + noise_std * np.random.randn(n_train)
    return X, y

X_train_np, y_train_np = generate_data()
X_test_np = np.linspace(-5, 5, 400)[:, None]

# Exemple de bruit hétéroscédastique
obs_error_std_np = 0.1 + 0.3 * np.abs(X_train_np.ravel()) / 3.0

X_train = torch.tensor(X_train_np, dtype=dtype)
y_train = torch.tensor(y_train_np, dtype=dtype)
obs_noise = torch.tensor(obs_error_std_np**2, dtype=dtype)

X_test = torch.tensor(X_test_np, dtype=dtype)

# -------------------------------------------------------
# 2) Définition du modèle GP (Exact GP)
# -------------------------------------------------------
class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel()
        )

    def forward(self, x):
        mean_x = self.mean_module(x)
        cov_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, cov_x)

# -------------------------------------------------------
# 3) Likelihood avec bruit hétéroscédastique
# -------------------------------------------------------
likelihood = gpytorch.likelihoods.FixedNoiseGaussianLikelihood(
    noise=obs_noise,
    learn_additional_noise=False  # empêche l'apprentissage d'un terme homoscédastique supplémentaire
)

model = ExactGPModel(X_train, y_train, likelihood).double()

model.train()
likelihood.train()

optimizer = torch.optim.Adam(model.parameters(), lr=0.1)
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

# -------------------------------------------------------
# 4) Entraînement des hyperparamètres du GP
# -------------------------------------------------------
training_iter = 300
for i in range(training_iter):
    optimizer.zero_grad()
    output = model(X_train)
    loss = -mll(output, y_train)
    loss.backward()
    if (i + 1) % 50 == 0:
        ls = model.covar_module.base_kernel.lengthscale.item()
        var = model.covar_module.outputscale.item()
        print(f"Iter {i+1}/{training_iter} | Loss={loss.item():.4f} | ls={ls:.3f} | var={var:.3f}")
    optimizer.step()

# -------------------------------------------------------
# 5) Évaluation
# -------------------------------------------------------
model.eval()
likelihood.eval()

with torch.no_grad():
    preds = likelihood(model(X_test))
    mean = preds.mean.numpy()
    lower, upper = preds.confidence_region()

# -------------------------------------------------------
# 6) Visualisation
# -------------------------------------------------------
plt.figure(figsize=(10, 5))
plt.scatter(X_train_np, y_train_np, c='k', s=40, label="Train")
plt.errorbar(
    X_train_np, y_train_np,
    yerr=obs_error_std_np,
    fmt='none',
    ecolor='gray',
    alpha=0.6,
    label="Bruit observé"
)

plt.plot(X_test_np, mean, label='Moyenne prédite')
plt.fill_between(
    X_test_np.ravel(),
    lower.numpy(),
    upper.numpy(),
    alpha=0.3,
    label="Intervalle de confiance 95%"
)
plt.title("GPyTorch — Régression GP avec bruit hétéroscédastique")
plt.legend()
plt.show()

Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?
Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?

Un autre exemple

Nous revisitons maintenant l’exemple de l’article précédent où nous avions implémenté un modèle de GP manuellement sans PyTorch.
Cela permet de comparer cette approche avec une implémentation propre utilisant GPyTorch.

Exemple 1D

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import torch
import gpytorch
import matplotlib.pyplot as plt

# 1) Données synthétiques
X = torch.tensor([1., 3., 5., 6., 7., 8.]).unsqueeze(1)
Y = X.squeeze() * torch.sin(X.squeeze())

# Nouveau : bruit différent par observation
obs_noise_std = torch.tensor([0.2, 0.1, 0.5, 0.3, 0.15, 0.4])
obs_noise_var = obs_noise_std**2

# 2) Modèle GP avec FixedNoiseGaussianLikelihood
class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel()
        )

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

likelihood = gpytorch.likelihoods.FixedNoiseGaussianLikelihood(
    noise=obs_noise_var
)
model = ExactGPModel(X, Y, likelihood)

# 3) Entraînement
model.train()
likelihood.train()
optimizer = torch.optim.Adam(model.parameters(), lr=0.1)
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

for i in range(120):
    optimizer.zero_grad()
    output = model(X)
    loss = -mll(output, Y)
    loss.backward()
    optimizer.step()

# 4) Prédiction
model.eval()
likelihood.eval()
Xtest = torch.linspace(0, 10, 200).unsqueeze(1)

with torch.no_grad(), gpytorch.settings.fast_pred_var():
    pred = likelihood(model(Xtest))

pred_mean = pred.mean.numpy()
pred_std = pred.stddev.numpy()

# 5) Visualisation
plt.figure(figsize=(9,5))
plt.errorbar(
    X.squeeze().numpy(),
    Y.numpy(),
    yerr=obs_noise_std.numpy(),
    fmt='o',
    capsize=4,
    label="Données d'entraînement (bruit hétéroscédastique)"
)
plt.plot(Xtest.numpy(), pred_mean, "r--", label="Moyenne GP")
plt.fill_between(
    Xtest.squeeze().numpy(),
    pred_mean - 1.96*pred_std,
    pred_mean + 1.96*pred_std,
    alpha=0.3,
    label="Intervalle 95%"
)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Régression GP (avec bruit par point)")
plt.legend()
plt.grid(True)
plt.show()

Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?
Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?

Exemple 2D

Modèle de régression GP (kernel RBF 2D)

\begin{equation}
k(x, x') = \sigma_f^2 \exp\left(-\frac{(x_1-x_1')^2}{2\ell_1^2}-\frac{(x_2-x_2')^2}{2\ell_2^2}\right)
\end{equation}

  • Régression 2D avec GP exact (RBF avec ARD)
  • Classification binaire avec GP variationnel sparse/approx.

Updated Code

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
import torch
import gpytorch
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

torch.manual_seed(0)

# -----------------------------
# Data (same as our original example)
# -----------------------------
X_np = np.array([
    [-8.0,-8.0], [-6.0,-3.0], [-7.0,2.0], [-4.0,4.0],
    [ 2.0, 3.0], [ 5.0, 7.0], [ 1.0,-1.0], [ 3.0,-4.0], [ 7.0,-7.0]
], dtype=np.float32)

# labels originally -1 / +1 -> convert to 0/1 for Bernoulli likelihood
Y_class_np = np.array([-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0], dtype=np.float32)
Y_class_01_np = (Y_class_np > 0).astype(np.float32)  # 0 or 1

# We'll also use Y_reg as same numeric labels for regression (float)
Y_reg_np = Y_class_np.copy()

# Torch tensors (float)
X = torch.from_numpy(X_np)
Y_reg = torch.from_numpy(Y_reg_np)
Y_class = torch.from_numpy(Y_class_01_np)  # 0/1

# -----------------------------
# Regression: Exact GP (RBF with ARD)
# -----------------------------
class ExactGPRegression(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        # RBF with ARD -> separate lengthscale per input dim
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel(ard_num_dims=train_x.shape[1])
        )

    def forward(self, x):
        mean = self.mean_module(x)
        cov = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean, cov)

# Setup likelihood & model
likelihood_reg = gpytorch.likelihoods.GaussianLikelihood()
model_reg = ExactGPRegression(X, Y_reg, likelihood_reg)

# Train regression model
model_reg.train()
likelihood_reg.train()
optimizer = torch.optim.Adam(model_reg.parameters(), lr=0.1)
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood_reg, model_reg)

n_iter_reg = 200
for i in range(n_iter_reg):
    optimizer.zero_grad()
    output = model_reg(X)
    loss = -mll(output, Y_reg)
    loss.backward()
    optimizer.step()
    if (i+1) % 50 == 0:
        print(f"[Reg] Iter {i+1}/{n_iter_reg}  NLML: {loss.item():.4f}")

# -----------------------------
# Regression: predictions on grid
# -----------------------------
model_reg.eval()
likelihood_reg.eval()

# grid settings
grid_x1 = np.arange(-10.0, 10.0, 0.12, dtype=np.float32)
grid_x2 = np.arange(-10.0, 10.0, 0.12, dtype=np.float32)
X1g, X2g = np.meshgrid(grid_x1, grid_x2, indexing='xy')
X_grid = torch.from_numpy(np.column_stack([X1g.ravel(), X2g.ravel()]).astype(np.float32))

with torch.no_grad(), gpytorch.settings.fast_pred_var():
    pred_reg = likelihood_reg(model_reg(X_grid))
    mean_grid = pred_reg.mean.reshape(X1g.shape).numpy()
    var_grid = pred_reg.variance.reshape(X1g.shape).numpy()

# Plot regression mean
plt.figure(figsize=(6,5))
plt.imshow(mean_grid, origin='lower', extent=[-10,10,-10,10], cmap=cm.jet)
plt.colorbar(label='Predictive mean')
plt.scatter(X_np[:,0], X_np[:,1], c=Y_reg_np, edgecolors='k', cmap=cm.coolwarm, s=60)
plt.title("GP Regression — Predictive Mean")
plt.xlabel("x1"); plt.ylabel("x2")
plt.grid(True, linestyle='--', alpha=0.4)
plt.show()

# Plot regression variance
plt.figure(figsize=(6,5))
plt.imshow(var_grid, origin='lower', extent=[-10,10,-10,10], cmap=cm.jet)
plt.colorbar(label='Predictive variance')
plt.scatter(X_np[:,0], X_np[:,1], c='white', edgecolors='k', s=60)
plt.title("GP Regression — Predictive Variance")
plt.xlabel("x1"); plt.ylabel("x2")
plt.grid(True, linestyle='--', alpha=0.4)
plt.show()

# -----------------------------
# Classification: Variational GP (Sparse/Approximate)
# -----------------------------
# Convert dataset for classifier: X (same), Y_class (0/1)
# Use inducing points (we'll use training inputs as inducing for this tiny example)
inducing_points = X.clone()

class VariationalGPClassifier(gpytorch.models.ApproximateGP):
    def __init__(self, inducing_points):
        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(inducing_points.size(0))
        variational_strategy = gpytorch.variational.VariationalStrategy(
            self, inducing_points, variational_distribution, learn_inducing_locations=True
        )
        super().__init__(variational_strategy)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel(ard_num_dims=inducing_points.shape[1])
        )

    def forward(self, x):
        mean = self.mean_module(x)
        cov = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean, cov)

# Likelihood and model
likelihood_clf = gpytorch.likelihoods.BernoulliLikelihood()
model_clf = VariationalGPClassifier(inducing_points)

# Put into training mode
model_clf.train()
likelihood_clf.train()

# Optimizer and ELBO
optimizer_clf = torch.optim.Adam(model_clf.parameters(), lr=0.05)
num_data = X.shape[0]
elbo = gpytorch.mlls.VariationalELBO(likelihood_clf, model_clf, num_data)

# Training loop
n_iter_clf = 500
for i in range(n_iter_clf):
    optimizer_clf.zero_grad()
    output = model_clf(X)
    loss = -elbo(output, Y_class)
    loss.backward()
    optimizer_clf.step()
    if (i+1) % 100 == 0:
        print(f"[Clf] Iter {i+1}/{n_iter_clf}  ELBO loss: {loss.item():.4f}")

# -----------------------------
# Classification: predictions on grid
# -----------------------------
model_clf.eval()
likelihood_clf.eval()

with torch.no_grad(), gpytorch.settings.fast_pred_var():
    preds_clf = likelihood_clf(model_clf(X_grid))
    probs = preds_clf.mean.reshape(X1g.shape).numpy()  # P(Y=1)

# Plot classification probability map
plt.figure(figsize=(6,5))
plt.imshow(probs, origin='lower', extent=[-10,10,-10,10], cmap=cm.jet, vmin=0, vmax=1)
plt.colorbar(label='P(Y=1)')
# plot training points (use marker styles like original)
for i, y in enumerate(Y_class_np):
    if y == -1.0:
        plt.scatter(X_np[i,0], X_np[i,1], marker='o', color='#1f77b4', edgecolors='k', s=80)
    else:
        plt.scatter(X_np[i,0], X_np[i,1], marker='x', color='#ff7f0e', s=80, linewidths=2)
plt.title("GP Classification — Predicted Probability P(Y=1)")
plt.xlabel("x1"); plt.ylabel("x2")
plt.grid(True, linestyle='--', alpha=0.4)
plt.show()

# Contour of decision surface
plt.figure(figsize=(6,5))
CS = plt.contour(X1g, X2g, probs, levels=np.linspace(0,1,11), cmap=cm.jet)
plt.clabel(CS, inline=True, fontsize=8)
plt.title("GP Classification — Probability Contours")
plt.xlabel("x1"); plt.ylabel("x2")
plt.grid(True, linestyle='--', alpha=0.4)
plt.scatter(X_np[:,0], X_np[:,1], c='k', s=30)
plt.show()

Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ? Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?
Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ? Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?
Comment Implémenter un Gaussian Process Simple en Python avec PyTorch ?

Références

Liens Site
PyTorch Documentation Documentation officielle PyTorch
Cholesky PyTorch Décomposition Cholesky
Optimisateurs PyTorch Optimiseurs PyTorch
Scikit-Learn GP Implémentation de référence des GP
GPyTorch Bibliothèque officielle GPyTorch
GPML Textbook Livre GPML (Rasmussen & Williams)
Wikipedia - Gaussian Process Vue d’ensemble des processus gaussiens
Creative Commons License

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Author profile-image

Benjamin

Développeur à mes heures perdues, je m’intéresse à la télédétection terrestre et à la détection des incendies par satellite. Passionné par Python, le machine learning et la science ouverte.

Skills
  • Remote sensing
  • Python
  • Machine Learning
  • Tensorflow 2.0
Image

of