4.1 Réseaux de neurones - Machine Learning en géosciences

Une régression linéaire est un neurone unique, entièrement connecté.

Figure : une unité logique à seuil (d’après A. Géron, Hands-On Machine Learning, 2019).

Un neurone artificiel unique calcule une somme pondérée de ses entrées plus un terme de biais, puis applique une fonction d’activation :

$y = f( \sum_{i=1}^3 w_i x_i + b)$ ,

où $y$ est la sortie, $w_i$ les poids, $b$ le biais du neurone (un par neurone, pas un par entrée), et $f$ une fonction d’activation (p. ex. sigmoïde, ReLU, etc.).

Le choix de $f$ donne son nom à l’unité. Avec une activation en échelon qui se déclenche quand la somme pondérée dépasse un seuil, le neurone est une unité logique à seuil (Threshold Logic Unit, TLU) : il produit une décision positive/négative tranchée. Avec l’activation sigmoïde, le même neurone calcule une régression logistique et produit une probabilité.

Un perceptron est une couche unique de TLU, chaque TLU étant connectée à toutes les entrées.

La régression logistique n’est rien d’autre qu’un classifieur à réseau de neurones à une couche.

Dans ce carnet, nous construisons notre premier réseau de neurones en PyTorch et l’entraînons sur un problème de classification à 4 classes. Le flux de travail comporte cinq étapes, qui se répètent dans chaque projet d’apprentissage profond (deep learning) :

Préparer le jeu de données.
Concevoir le modèle.
Choisir la fonction de perte.
Choisir l’optimiseur.
Écrire la boucle d’entraînement.

🖥️ Diapositives du cours — Séance 21 (mer. 18 nov.)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import torch
import torch.nn as nn
import torch.optim as optim

# Device-agnostic setup: use a GPU if one is available.
device = torch.device("cuda" if torch.cuda.is_available()
                      else "mps" if torch.backends.mps.is_available()
                      else "cpu")
print("Using device:", device)

torch.manual_seed(42)
np.random.seed(42)

Using device: mps

1. Le jeu de données en PyTorch¶

Nous réutilisons le jeu de données de sources sismiques du chapitre 3 : 62 caractéristiques physiques (statistiques de forme d’onde, d’enveloppe et de spectre) extraites de sismogrammes de quatre types de sources enregistrés dans le Nord-Ouest Pacifique des États-Unis : séisme (earthquake), explosion, bruit (noise) et événement de surface (surface event). Au chapitre 3, vous avez classé ces événements avec l’apprentissage automatique classique (les forêts aléatoires atteignaient environ 80 à 90 % d’exactitude). Ici, nous entraînons un réseau de neurones sur le même problème, pour que vous puissiez comparer directement les deux approches.

Les données résident sur Zenodo sous forme de quatre fichiers CSV, un par classe, de 1 000 événements chacun. Nous les téléchargeons avec pooch, qui met les fichiers en cache localement.

1.1 Charger et étiqueter les données¶

import pooch

base = "https://zenodo.org/api/records/14025693/files"
classes = ["earthquakes", "explosion", "noise", "surface"]

frames = []
for c in classes:
    fname = f"1000_{c}_physical_features.csv"
    path = pooch.retrieve(f"{base}/{fname}/content", known_hash=None, fname=fname)
    frames.append(pd.read_csv(path))

df = pd.concat(frames, ignore_index=True)
print(df.shape)
print(df["source"].value_counts())

(4000, 65)
source
earthquake       1000
explosion        1000
noise            1000
surface event    1000
Name: count, dtype: int64

1.2 Nettoyer la table de caractéristiques¶

La table compte 65 colonnes. Trois ne sont pas des caractéristiques : Unnamed: 0 (un index de lignes résiduel), serial_no (un identifiant d’événement) et source (le nom de la classe, qui devient notre étiquette). Nous convertissons les quatre noms de classes en entiers de 0 à 3, le format d’étiquettes qu’attend la fonction de perte de classification de PyTorch.

Nous vérifions aussi les valeurs manquantes et infinies. Une colonne de caractéristiques, Gamma, est vide (entièrement NaN) dans ce jeu de données : nous supprimons donc la colonne plutôt que de supprimer chaque ligne ou d’inventer des valeurs par imputation. Après nettoyage, il ne reste ni NaN ni inf, et nous conservons 61 caractéristiques numériques.

class_names = ["earthquake", "explosion", "noise", "surface event"]
labels = df["source"].map({name: i for i, name in enumerate(class_names)}).to_numpy()

features = df.drop(columns=["Unnamed: 0", "serial_no", "source"])

# Columns that are entirely NaN carry no information: drop them.
all_nan = features.columns[features.isna().all()]
print("Dropping all-NaN columns:", list(all_nan))
features = features.drop(columns=all_nan)

print("Remaining NaN:", int(features.isna().sum().sum()),
      "- inf:", int(np.isinf(features.to_numpy()).sum()))

X = features.to_numpy(dtype=np.float32)
n_features = X.shape[1]
print("Feature matrix:", X.shape, "- labels:", labels.shape)

Dropping all-NaN columns: ['Gamma']
Remaining NaN: 0 - inf: 0
Feature matrix: (4000, 61) - labels: (4000,)

1.3 Découper, puis normaliser¶

Nous découpons les données en trois : entraînement (60 %) pour ajuster les poids, validation (20 %) pour surveiller la généralisation pendant l’entraînement, et test (20 %) mis de côté pour l’évaluation finale. Le découpage est stratifié, pour que chaque sous-ensemble garde le même équilibre de classes.

Les réseaux de neurones s’entraînent mal quand les caractéristiques d’entrée couvrent des plages très différentes : nous standardisons donc chaque caractéristique à moyenne nulle et variance unité. La normalisation est ajustée sur l’ensemble d’entraînement seulement, puis appliquée aux ensembles de validation et de test. L’ajuster sur toutes les données ferait fuir de l’information des ensembles mis de côté vers l’entraînement.

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# 60% train, 20% validation, 20% test, stratified by class.
X_train, X_tmp, y_train, y_tmp = train_test_split(
    X, labels, test_size=0.4, random_state=42, stratify=labels)
X_val, X_test, y_val, y_test = train_test_split(
    X_tmp, y_tmp, test_size=0.5, random_state=42, stratify=y_tmp)

scaler = StandardScaler().fit(X_train)          # fit on train only
X_train = scaler.transform(X_train).astype(np.float32)
X_val = scaler.transform(X_val).astype(np.float32)
X_test = scaler.transform(X_test).astype(np.float32)

print("train:", X_train.shape, " val:", X_val.shape, " test:", X_test.shape)

train: (2400, 61)  val: (800, 61)  test: (800, 61)

1.4 Créer une classe Dataset personnalisée¶

PyTorch fournit la classe Dataset, que vous pouvez utiliser pour envelopper vos tableaux. Vous devez implémenter les méthodes __len__ et __getitem__. La classe convertit aussi les tableaux NumPy en tensors (tenseurs) : float32 pour les caractéristiques, int64 (long) pour les étiquettes de classes.

from torch.utils.data import Dataset

class SeismicFeatures(Dataset):
    def __init__(self, data, labels):
        self.data = torch.tensor(data, dtype=torch.float32)
        self.labels = torch.tensor(labels, dtype=torch.long)

    def __len__(self):
        return len(self.data)

    def __getitem__(self, index):
        return self.data[index], self.labels[index]

train_dataset = SeismicFeatures(X_train, y_train)
val_dataset = SeismicFeatures(X_val, y_val)
test_dataset = SeismicFeatures(X_test, y_test)

sample_x, sample_y = train_dataset[0]
print("one sample:", sample_x.shape, sample_x.dtype, "- label:", sample_y.item())

one sample: torch.Size([61]) torch.float32 - label: 1

1.5 Créer un DataLoader¶

Une fois le jeu de données en place, un DataLoader gère la constitution des lots (batches) et le mélange. Nous mélangeons l’ensemble d’entraînement à chaque époque, pour que les lots diffèrent d’une fois sur l’autre ; les ensembles de validation et de test n’ont pas besoin d’être mélangés.

from torch.utils.data import DataLoader

batch_size = 128
train_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True)
val_loader = DataLoader(val_dataset, batch_size=batch_size, shuffle=False)
test_loader = DataLoader(test_dataset, batch_size=batch_size, shuffle=False)

for batch_data, batch_labels in train_loader:
    print("batch data:", batch_data.shape)
    print("batch labels:", batch_labels.shape)
    break

batch data: torch.Size([128, 61])
batch labels: torch.Size([128])

2. Concevoir le modèle¶

Notre premier réseau de neurones a une seule couche cachée : une application linéaire des 61 caractéristiques d’entrée vers 32 unités cachées, une activation ReLU, et une application linéaire des unités cachées vers 4 scores de sortie, un par classe.

Les sorties brutes, non normalisées, s’appellent des logits. Nous n’appliquons pas de softmax dans le modèle : la perte d’entropie croisée de PyTorch (section 3) l’applique en interne, ce qui est à la fois plus sûr numériquement et la pratique standard.

Un modèle est une sous-classe de nn.Module. La méthode __init__ déclare les couches ; la méthode forward définit comment les données les traversent.

class NN1(nn.Module):

    def __init__(self, n_features, n_hidden, n_classes):
        # super() makes NN1 inherit the machinery of nn.Module
        super().__init__()
        self.hidden = nn.Linear(n_features, n_hidden)  # z1 = W1 x + b1
        self.activation = nn.ReLU()                    # a1 = max(0, z1)
        self.output = nn.Linear(n_hidden, n_classes)   # logits = W2 a1 + b2

    def forward(self, x):
        x = self.activation(self.hidden(x))
        return self.output(x)   # raw logits; softmax lives in the loss

Nous créons un modèle de la classe NN1 et le déplaçons vers l’appareil de calcul avec .to(device). Tous les tenseurs que nous lui fournissons doivent résider sur le même périphérique.

model = NN1(n_features=n_features, n_hidden=32, n_classes=4).to(device)
print(model)

n_params = sum(p.numel() for p in model.parameters())
print(f"trainable parameters: {n_params}")

NN1(
  (hidden): Linear(in_features=61, out_features=32, bias=True)
  (activation): ReLU()
  (output): Linear(in_features=32, out_features=4, bias=True)
)
trainable parameters: 2116

3. Fonction de perte¶

Pour la classification binaire, nous avons défini la fonction logistique ou sigmoïde,

$\hat{y} = \frac{1}{1 + e^{-z}}$ ,

comme fonction de probabilité d’un classifieur binaire.

La perte logarithmique (log loss) à minimiser est :

$\mathcal{L} = - (y \log \hat{y} + (1 - y) \log (1 - \hat{y}))$ avec $\hat{y} = \frac{1}{1 + e^{-z}}$ .

Nous cherchons les poids $w$ et les biais $b$ tels que $P(Y = y)$ soit maximale, ce qui équivaut à minimiser la perte $\mathcal{L}$ .

Pour une classification multi-classes à $K$ classes, la fonction sigmoïde s’étend en une fonction softmax,

$P(Y = k) = \frac{\exp(w_k x + b_k)}{\sum_{j = 1}^K \exp(w_j x + b_j)}$ .

Une fois trouvées les valeurs de $w_k$ et $b_k$ pour $k = 1, \cdots, K$ , nous pouvons calculer les valeurs de $P(Y = k)$ . Nous cherchons alors la valeur de $k$ pour laquelle $P(Y = k)$ est maximale et classons l’échantillon dans la classe $k$ . Ici $K = 4$ .

De même, nous définissons la perte à $K$ classes :

$\mathcal{L} = - \sum_{k = 1}^K y_k \log \frac{\exp(w_k x + b_k)}{\sum_{j = 1}^K \exp(w_j x + b_j)}$ ,

où $y_k = 1$ si la vraie étiquette de l’échantillon $x$ est $k$ , et $y_k = 0$ sinon.

C’est la perte d’entropie croisée (cross-entropy loss). En PyTorch, nn.CrossEntropyLoss combine le softmax et la perte logarithmique en une seule fonction, raison pour laquelle le modèle sort des logits bruts.

criterion = nn.CrossEntropyLoss()
criterion

CrossEntropyLoss()

4. Optimisation¶

La descente de gradient utilise le jeu de données entier pour calculer le gradient et avancer vers le minimum de la perte. De fait, la taille de lot est le jeu de données entier.

La descente de gradient par mini-lots (mini-batch) utilise des lots de données pour calculer le gradient, met à jour les poids, puis passe au lot suivant. La perte est moyennée sur les échantillons du lot avant le calcul du gradient.

La descente de gradient stochastique (SGD) utilise un échantillon à la fois. De fait, la taille de lot vaut 1.

Nous utilisons des mini-lots avec l’optimiseur Adam, une variante de la descente de gradient qui adapte la taille du pas pour chaque paramètre à partir de moyennes glissantes du gradient et de son carré. Adam converge généralement plus vite que la SGD simple et constitue le choix par défaut pour un premier essai.

Les paramètres à choisir sont :

l’optimiseur (SGD, Adam ou autres),
le taux d’apprentissage (l’échelle du pas effectué le long du gradient).

learning_rate = 1e-3
optimizer = optim.Adam(model.parameters(), lr=learning_rate)

5. Entraînement¶

Les paramètres à choisir sont :

la taille de lot (batch size) : le nombre d’échantillons utilisés pour estimer une perte moyenne avant chaque mise à jour des poids.
le nombre d’époques : le nombre de passages sur l’ensemble d’entraînement entier. Si l’ensemble d’entraînement est découpé en lots de taille batch_size, chaque époque parcourt tous les lots. Les modèles s’entraînent typiquement pendant des dizaines à des milliers d’époques.

Nous écrivons la boucle d’entraînement sous forme de fonction. À chaque époque, elle :

parcourt les lots d’entraînement : passe avant, perte, passe arrière, mise à jour des poids ;
évalue la perte et l’exactitude sur l’ensemble de validation avec les gradients désactivés (torch.no_grad()), puisqu’aucun apprentissage n’y a lieu.

Nous enregistrons la perte et l’exactitude sur les deux ensembles, pour tracer ensuite les courbes d’apprentissage.

def evaluate(model, loader, criterion):
    """Mean loss and accuracy (%) of the model over one DataLoader."""
    model.eval()
    total_loss, correct, total = 0.0, 0, 0
    with torch.no_grad():
        for inputs, targets in loader:
            inputs, targets = inputs.to(device), targets.to(device)
            outputs = model(inputs)
            total_loss += criterion(outputs, targets).item() * targets.size(0)
            predicted = outputs.argmax(dim=1)
            correct += (predicted == targets).sum().item()
            total += targets.size(0)
    return total_loss / total, 100 * correct / total


def train(model, n_epochs, train_loader, val_loader, criterion, optimizer):
    history = {"train_loss": [], "train_acc": [], "val_loss": [], "val_acc": []}

    for epoch in range(n_epochs):
        model.train()
        for inputs, targets in train_loader:
            inputs, targets = inputs.to(device), targets.to(device)
            optimizer.zero_grad()          # reset gradients
            outputs = model(inputs)        # forward pass
            loss = criterion(outputs, targets)
            loss.backward()                # backward pass
            optimizer.step()               # update the weights

        # Evaluate on both sets at the end of the epoch
        train_loss, train_acc = evaluate(model, train_loader, criterion)
        val_loss, val_acc = evaluate(model, val_loader, criterion)
        history["train_loss"].append(train_loss)
        history["train_acc"].append(train_acc)
        history["val_loss"].append(val_loss)
        history["val_acc"].append(val_acc)

        if (epoch + 1) % 5 == 0 or epoch == 0:
            print(f"[Epoch {epoch + 1:3d}] train loss: {train_loss:.3f} "
                  f"- val loss: {val_loss:.3f} - val accuracy: {val_acc:.1f}%")

    return history

Nous entraînons pendant 40 époques, ce qui prend quelques secondes. Vous pouvez augmenter n_epochs ou la taille de la couche cachée sur votre propre machine, pour voir jusqu’où ce petit modèle peut aller.

n_epochs = 40
history = train(model, n_epochs, train_loader, val_loader, criterion, optimizer)

[Epoch   1] train loss: 1.157 - val loss: 1.150 - val accuracy: 71.0%

[Epoch   5] train loss: 0.603 - val loss: 0.586 - val accuracy: 81.4%

[Epoch  10] train loss: 0.462 - val loss: 0.461 - val accuracy: 84.9%

[Epoch  15] train loss: 0.401 - val loss: 0.423 - val accuracy: 86.1%

[Epoch  20] train loss: 0.365 - val loss: 0.405 - val accuracy: 86.1%

[Epoch  25] train loss: 0.339 - val loss: 0.388 - val accuracy: 86.9%

[Epoch  30] train loss: 0.320 - val loss: 0.378 - val accuracy: 87.8%

[Epoch  35] train loss: 0.302 - val loss: 0.373 - val accuracy: 88.5%

[Epoch  40] train loss: 0.287 - val loss: 0.369 - val accuracy: 88.5%

Traçons les courbes d’apprentissage : la perte sur l’axe de gauche, l’exactitude sur l’axe de droite, pour l’ensemble d’entraînement et l’ensemble de validation. Les traits pleins correspondent à l’ensemble d’entraînement, les traits tiretés à l’ensemble de validation. Quand les deux courbes se suivent, le modèle généralise ; quand la courbe de validation s’aplatit ou s’inverse tandis que la courbe d’entraînement continue de s’améliorer, le modèle commence à surapprendre.

epochs = np.arange(1, n_epochs + 1)
fig, ax1 = plt.subplots(figsize=(7, 4))

loss_color = "tab:red"
ax1.set_xlabel("Epoch")
ax1.set_ylabel("Cross-entropy loss", color=loss_color)
ax1.plot(epochs, history["train_loss"], color=loss_color, label="train loss")
ax1.plot(epochs, history["val_loss"], color=loss_color, linestyle="--", label="validation loss")
ax1.tick_params(axis="y", labelcolor=loss_color)

ax2 = ax1.twinx()
acc_color = "tab:blue"
ax2.set_ylabel("Accuracy (%)", color=acc_color)
ax2.plot(epochs, history["train_acc"], color=acc_color, label="train accuracy")
ax2.plot(epochs, history["val_acc"], color=acc_color, linestyle="--", label="validation accuracy")
ax2.tick_params(axis="y", labelcolor=acc_color)

lines = ax1.get_lines() + ax2.get_lines()
ax1.legend(lines, [l.get_label() for l in lines], loc="center right", fontsize=9)
ax1.set_title("Learning curves")
fig.tight_layout()
plt.show()

6. Évaluer sur l’ensemble de test¶

L’ensemble de validation a guidé l’entraînement ; l’ensemble de test donne le score final, sans biais. Nous rapportons l’exactitude globale et la matrice de confusion, qui montre quelles classes se confondent.

from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay

test_loss, test_acc = evaluate(model, test_loader, criterion)
print(f"Test loss: {test_loss:.3f} - test accuracy: {test_acc:.1f}%")

# Predictions for the whole test set
model.eval()
with torch.no_grad():
    logits = model(torch.tensor(X_test, dtype=torch.float32).to(device))
    y_pred = logits.argmax(dim=1).cpu().numpy()

cm = confusion_matrix(y_test, y_pred)
disp = ConfusionMatrixDisplay(cm, display_labels=class_names)
fig, ax = plt.subplots(figsize=(5.5, 5))
disp.plot(ax=ax, cmap="Blues", colorbar=False, xticks_rotation=45)
ax.set_title("Test-set confusion matrix")
fig.tight_layout()
plt.show()

Test loss: 0.375 - test accuracy: 87.1%

Le bruit se sépare nettement du reste. L’essentiel de la confusion se produit entre les sources sismiques, surtout entre séismes et explosions, qui produisent des formes d’onde similaires. Comparez cette matrice avec vos résultats du chapitre 3 sur les mêmes données : un réseau à une couche cachée et une forêt aléatoire bien réglée aboutissent dans une plage similaire, ce qui est courant sur de petits jeux de données tabulaires.

Résumé¶

Vous avez construit un flux de travail PyTorch complet : un Dataset et un DataLoader pour les caractéristiques de sources sismiques, un réseau à une couche cachée, la perte d’entropie croisée, l’optimiseur Adam, et une boucle d’entraînement avec courbes d’apprentissage et évaluation finale sur un ensemble de test mis de côté. Dans le carnet 4.2, nous empilons davantage de couches cachées en un perceptron multicouche (MLP) et ajoutons de la régularisation pour contrôler le surapprentissage. Dans le carnet 4.5, nous réglons l’architecture et les hyperparamètres d’entraînement de façon systématique, au lieu de les choisir à la main.