4.1 Redes neuronales - Machine Learning en las Geociencias

Una regresión lineal es una sola neurona totalmente conectada.

Figura: una unidad lógica de umbral (según A. Géron, Hands-On Machine Learning, 2019).

Una neurona artificial calcula una suma ponderada de sus entradas más un término de sesgo, y luego aplica una función de activación:

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

donde $y$ es la salida, $w_i$ son los pesos, $b$ es el sesgo de la neurona (uno por neurona, no uno por entrada) y $f$ es una función de activación (p. ej., sigmoide, ReLU, etc.).

La elección de $f$ da nombre a la unidad. Con una activación escalón que se dispara cuando la suma ponderada supera un umbral, la neurona es una unidad lógica de umbral (TLU, Threshold Logic Unit): produce una decisión dura positivo/negativo. Con la activación sigmoide, la misma neurona calcula una regresión logística y produce, en cambio, una probabilidad.

Un perceptrón es una sola capa de TLU, con cada TLU conectada a todas las entradas.

La regresión logística es simplemente un clasificador de red neuronal de una capa.

En este cuaderno construimos nuestra primera red neuronal en PyTorch y la entrenamos en un problema de clasificación de 4 clases. El flujo de trabajo tiene cinco pasos que se repiten en todo proyecto de aprendizaje profundo (deep learning):

Preparar el conjunto de datos.
Diseñar el modelo.
Elegir la función de pérdida.
Elegir el optimizador.
Escribir el bucle de entrenamiento.

🖥️ Diapositivas — Sesión 21 (mié 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. El conjunto de datos en PyTorch¶

Reutilizamos el conjunto de datos de fuentes sísmicas del capítulo 3: 62 características físicas (estadísticas de la forma de onda, de la envolvente y del espectro) extraídas de sismogramas de cuatro tipos de fuente registrados en el noroeste del Pacífico estadounidense: sismo (earthquake), explosión, ruido y evento superficial. En el capítulo 3 usted clasificó estos eventos con aprendizaje automático clásico (los bosques aleatorios alcanzaron aproximadamente 80-90 % de exactitud). Aquí entrenamos una red neuronal en el mismo problema, para que pueda comparar los dos enfoques directamente.

Los datos viven en Zenodo como cuatro archivos CSV, uno por clase, con 1 000 eventos cada uno. Los descargamos con pooch, que guarda los archivos en una caché local.

1.1 Cargar y etiquetar los datos¶

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 Limpiar la tabla de características¶

La tabla tiene 65 columnas. Tres no son características: Unnamed: 0 (un índice de fila remanente), serial_no (un identificador de evento) y source (el nombre de la clase, que se convierte en nuestra etiqueta). Mapeamos los cuatro nombres de clase a enteros 0-3, que es el formato de etiqueta que espera la pérdida de clasificación de PyTorch.

También revisamos si hay valores faltantes e infinitos. Una columna de características, Gamma, está vacía (toda NaN) en este conjunto de datos, así que eliminamos la columna en lugar de eliminar todas las filas o inventar valores por imputación. Después de la limpieza no quedan NaN ni inf, y conservamos 61 características numéricas.

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 Dividir y luego escalar¶

Dividimos los datos en tres partes: entrenamiento (60 %) para ajustar los pesos, validación (20 %) para monitorear la generalización durante el entrenamiento y prueba (20 %), reservada para la evaluación final. La división es estratificada, de modo que cada subconjunto conserva el mismo balance de clases.

Las redes neuronales entrenan mal cuando las características de entrada abarcan rangos muy distintos, así que estandarizamos cada característica a media cero y varianza unitaria. El escalador se ajusta solo con el conjunto de entrenamiento y luego se aplica a los conjuntos de validación y prueba. Ajustarlo con todos los datos filtraría información de los conjuntos reservados hacia el entrenamiento.

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 Crear una clase Dataset propia¶

PyTorch ofrece la clase Dataset, que puede usar para envolver sus arreglos. Debe implementar los métodos __len__ y __getitem__. La clase también convierte los arreglos de NumPy en tensors: float32 para las características, int64 (long) para las etiquetas de clase.

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 Crear un DataLoader¶

Una vez que tiene un conjunto de datos, un DataLoader se encarga de armar los lotes (batches) y de barajar. Barajamos el conjunto de entrenamiento en cada época para que los lotes sean distintos cada vez; los conjuntos de validación y prueba no necesitan barajarse.

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. Diseñar el modelo¶

Nuestra primera red neuronal tiene una sola capa oculta: un mapa lineal de las 61 características de entrada a 32 unidades ocultas, una activación ReLU y un mapa lineal de las unidades ocultas a 4 puntajes de salida, uno por clase.

Las salidas crudas se llaman logits. No aplicamos softmax dentro del modelo: la pérdida de entropía cruzada de PyTorch (sección 3) la aplica internamente, lo que es a la vez numéricamente más seguro y la práctica estándar.

Un modelo es una subclase de nn.Module. El método __init__ declara las capas; el método forward define cómo fluyen los datos a través de ellas.

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

Creamos un modelo de la clase NN1 y lo movemos al dispositivo de cómputo con .to(device). Todos los tensores que le pasemos deben vivir en el mismo dispositivo.

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. Función de pérdida¶

Para la clasificación binaria definimos la función logística o sigmoide,

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

como la función de probabilidad de un clasificador binario.

La función de pérdida logarítmica a minimizar es:

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

Queremos encontrar los pesos $w$ y los sesgos $b$ tales que $P(Y = y)$ sea máxima, lo que equivale a minimizar la pérdida $\mathcal{L}$ .

Para una clasificación multiclase con $K$ clases, la función sigmoide se extiende a una función softmax,

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

Una vez encontrados los valores de $w_k$ y $b_k$ para $k = 1, \cdots, K$ , podemos calcular los valores de $P(Y = k)$ . Luego buscamos el valor de $k$ para el cual $P(Y = k)$ es máximo y clasificamos la muestra como de clase $k$ . Aquí $K = 4$ .

De manera similar, definimos la pérdida con $K$ clases:

$\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)}$ ,

donde $y_k = 1$ si la etiqueta verdadera de la muestra $x$ es $k$ , y $y_k = 0$ en caso contrario.

Esto se llama la pérdida de entropía cruzada (cross-entropy loss). En PyTorch, nn.CrossEntropyLoss combina la softmax y la pérdida logarítmica en una sola función, y por eso el modelo produce logits crudos.

criterion = nn.CrossEntropyLoss()
criterion

CrossEntropyLoss()

4. Optimización¶

El descenso de gradiente usa el conjunto de datos completo para calcular el gradiente y dar el paso hacia el mínimo de la pérdida. En la práctica, el tamaño de lote es el conjunto de datos entero.

El descenso de gradiente por minilotes usa lotes de datos para calcular el gradiente, actualiza los pesos y pasa al siguiente lote. La pérdida se promedia sobre las muestras del lote antes de calcular el gradiente.

El descenso de gradiente estocástico (SGD) usa una muestra a la vez. En la práctica, el tamaño de lote es 1.

Usamos minilotes con el optimizador Adam, una variante del descenso de gradiente que adapta el tamaño del paso para cada parámetro a partir de promedios móviles del gradiente y de su cuadrado. Adam suele converger más rápido que el SGD simple y es la elección por defecto para un primer intento.

Los parámetros a elegir son:

el optimizador (SGD, Adam u otros),
la tasa de aprendizaje (la escala del paso que se da a lo largo del gradiente).

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

5. Entrenamiento¶

Los parámetros a elegir son:

tamaño de lote (batch size): el número de muestras usadas para estimar una pérdida promedio antes de cada actualización de pesos.
número de épocas: el número de pasadas sobre el conjunto de entrenamiento completo. Si el conjunto de entrenamiento se divide en lotes de tamaño batch_size, cada época recorre todos los lotes. Los modelos se entrenan típicamente durante decenas a miles de épocas.

Escribimos el bucle de entrenamiento como una función. En cada época:

recorre los lotes de entrenamiento: pasada hacia adelante, pérdida, pasada hacia atrás, actualización de pesos;
evalúa la pérdida y la exactitud en el conjunto de validación con los gradientes desactivados (torch.no_grad()), ya que ahí no ocurre aprendizaje.

Registramos la pérdida y la exactitud en ambos conjuntos para poder graficar después las curvas de aprendizaje.

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

Entrenamos durante 40 épocas, lo que toma unos segundos. Puede aumentar n_epochs o el tamaño de la capa oculta en su propia máquina para ver hasta dónde llega este modelo pequeño.

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%

Grafiquemos las curvas de aprendizaje: la pérdida en el eje izquierdo, la exactitud en el eje derecho, para los conjuntos de entrenamiento y de validación. Las líneas continuas son el conjunto de entrenamiento; las discontinuas, el de validación. Cuando las dos curvas se siguen mutuamente, el modelo generaliza; cuando la curva de validación se aplana o cambia de dirección mientras la de entrenamiento sigue mejorando, el modelo empieza a sobreajustarse.

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. Evaluar en el conjunto de prueba¶

El conjunto de validación guio el entrenamiento; el conjunto de prueba da la calificación final, sin sesgo. Reportamos la exactitud global y la matriz de confusión, que muestra qué clases se confunden entre sí.

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%

El ruido se separa limpiamente del resto. La mayor parte de la confusión ocurre entre las fuentes sísmicas, sobre todo entre sismos y explosiones, que producen formas de onda similares. Compare esta matriz con sus resultados del capítulo 3 sobre los mismos datos: una red de una capa oculta y un bosque aleatorio bien ajustado quedan en un rango similar, lo que es común en conjuntos de datos tabulares pequeños.

Resumen¶

Usted construyó un flujo de trabajo completo en PyTorch: un Dataset y un DataLoader para las características de fuentes sísmicas, una red de una capa oculta, la pérdida de entropía cruzada, el optimizador Adam y un bucle de entrenamiento con curvas de aprendizaje y una evaluación final en un conjunto de prueba reservado. En el cuaderno 4.2 apilamos más capas ocultas en un perceptrón multicapa y agregamos regularización para controlar el sobreajuste. En el cuaderno 4.5 ajustamos la arquitectura y los hiperparámetros de entrenamiento de manera sistemática en lugar de elegirlos a mano.