Sequence Modeling with Recurrent Neural Networks (RNN, LSTM, GRU) - 2026

Course: Machine Learning for Geologists & Geocomputation (Spatial Ecology 2026)
Instructor: Raffaele Scaringi
Dataset: NOAA Global Historical Climatology - Seattle Daily Meteorological Observations
Format: Single Self-Contained Lecture & Practical Notebook

Course Progression Recap

  • Class 1-3: Perceptrons, Multi-Layer Perceptrons (MLPs), and Optimizers (Gradient Descent, Momentum, Adam).

  • Class 4: Deep Feedforward Networks on Tabular Earth Data.

  • Class 5: Convolutional Neural Networks (CNNs) for Spatial Satellite Imagery (Sentinel-2 EuroSAT).

  • Class 6 (Today): Sequence Modeling & Recurrent Architectures for dynamic temporal and depth-dependent geological processes.

Pedagogical Learning Objectives

  1. The Sequence Modeling Paradigm: Understand why 2D feedforward networks fail on temporal and spatial sequences, and how recurrent connections maintain an internal hidden memory state.

  2. Real-World Environmental Data Ingestion: Load real-world daily meteorological observations directly from a public URL into Pandas without external credentials.

  3. Exploratory Data Analysis (EDA): Inspect atmospheric dynamics, seasonal thermal cycles, and the heavy-tailed intermittent nature of precipitation.

  4. Strict Temporal Validation (Zero Data Leakage): Why standard randomized cross-validation (shuffle=True) is fatal for geoscientific time series, and how to execute proper chronological block splitting.

  5. Geoscientific Feature Preprocessing: Apply physical transforms (e.g., logarithmic scaling for zero-inflated precipitation) and feature normalization fitted strictly on historical observations.

  6. PyTorch 3D Sequence Windowing: Master the PyTorch 3D recurrent tensor contract: (Batch_Size, Sequence_Length, Num_Features).

  7. Vanilla RNN & The Vanishing Gradient Pathology: Build a Vanilla RNN from scratch with torch.nn.RNN and mathematically/empirically observe gradient decay during Backpropagation Through Time (BPTT).

  8. Defensive Engineering with Gradient Clipping: Prevent exploding gradients using torch.nn.utils.clip_grad_norm_.

  9. Long Short-Term Memory (LSTM): Deconstruct the LSTM cell architecture: Forget, Input, and Output gates protecting the linear Cell State Highway (Constant Error Carousel).

  10. Gated Recurrent Unit (GRU): Compare GRU’s merged hidden/cell state and evaluate its parameter efficiency vs. LSTM.

  11. Geoscientific Model Evaluation: Benchmark models against the Persistence Baseline using hydrometeorological metrics: RMSE, MAE, \(R^2\), and the Nash-Sutcliffe Efficiency (:math:`NSE`).

  12. Physical Interpretability: Probe the internal cell state channels of the trained LSTM to reveal how deep recurrent layers store seasonal energy vs. short-term synoptic fronts.

  13. Bidirectional Recurrent Networks (BiLSTM) for Spatial Sequences: Understand why spatial geological profiles (borehole well-logging, stratigraphic sequences) require non-causal bidirectional context.

[1]:

# 1. Imports, Determinism & Device Configuration
import os
import random
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import Dataset, DataLoader

# Configure Matplotlib styling for high-contrast geoscientific figures
plt.rcParams['figure.dpi'] = 120
plt.rcParams['font.size'] = 10
plt.rcParams['axes.grid'] = True
plt.rcParams['grid.alpha'] = 0.3
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False

# Ensure complete reproducibility
SEED = 42
random.seed(SEED)
np.random.seed(SEED)
torch.manual_seed(SEED)
if torch.cuda.is_available():
    torch.cuda.manual_seed_all(SEED)

device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(f"[INFO] PyTorch Version: {torch.__version__}")
print(f"[INFO] Computing Device: {device}")

[INFO] PyTorch Version: 2.12.0+cu130
[INFO] Computing Device: cuda

1. Real-World Environmental Data Acquisition

In geoscience and environmental computing, synthetic toy data often masks the true complexities of natural processes:

  • Intermittency: Rainfall consists of prolonged zeros punctuated by intense non-linear bursts.

  • Thermal Inertia & Atmospheric Memory: Surface temperatures do not vary randomly; soil thermal mass and synoptic pressure systems store memory across days and weeks.

  • Seasonal Periodicity: Strong annual solar forcing superimposed with daily weather noise.

The Dataset: NOAA Seattle Daily Weather Record (2012-2015)

We acquire a real-world multi-variate meteorological record directly from a public permanent CDN:

  • Source: National Oceanic and Atmospheric Administration (NOAA) Global Historical Climatology Network.

  • Span: 4 continuous calendar years (2012-01-01 to 2015-12-31 = 1,461 daily records).

  • Observed Variables:

    • date: Daily observation timestamp.

    • precipitation: Total daily rainfall (mm).

    • temp_max: Maximum daily air temperature (deg C).

    • temp_min: Minimum daily air temperature (deg C).

    • wind: Mean daily wind speed (m/s).

    • weather: Qualitative weather summary (drizzle, rain, sun, snow, fog).

[2]:

# 2. Ingest Dataset Directly from Public Source URL
DATA_URL = "https://raw.githubusercontent.com/vega/vega-datasets/master/data/seattle-weather.csv"

print(f"[INFO] Fetching real-world meteorological data from: {DATA_URL}")
df_raw = pd.read_csv(DATA_URL)
df_raw['date'] = pd.to_datetime(df_raw['date'])
df_raw = df_raw.sort_values('date').reset_index(drop=True)

print(f"[INFO] Successfully loaded {len(df_raw)} daily observations.")
print(f"[INFO] Date range: {df_raw['date'].min().strftime('%Y-%m-%d')} to {df_raw['date'].max().strftime('%Y-%m-%d')}")
print("[INFO] Missing values across all columns:")
print(df_raw.isnull().sum())

# Inspect first 5 records
display(df_raw.head())

[INFO] Fetching real-world meteorological data from: https://raw.githubusercontent.com/vega/vega-datasets/master/data/seattle-weather.csv
[INFO] Successfully loaded 1461 daily observations.
[INFO] Date range: 2012-01-01 to 2015-12-31
[INFO] Missing values across all columns:
date             0
precipitation    0
temp_max         0
temp_min         0
wind             0
weather          0
dtype: int64
date precipitation temp_max temp_min wind weather
0 2012-01-01 0.0 12.8 5.0 4.7 drizzle
1 2012-01-02 10.9 10.6 2.8 4.5 rain
2 2012-01-03 0.8 11.7 7.2 2.3 rain
3 2012-01-04 20.3 12.2 5.6 4.7 rain
4 2012-01-05 1.3 8.9 2.8 6.1 rain

2. Exploratory Data Analysis & Geological/Atmospheric Dynamics

Before training any neural network, an Earth scientist must visualize the raw data to understand physical dynamics:

  1. Seasonal Thermal Envelope: The smooth sinusoidal oscillation of solar radiation creates the temperature envelope (\(T_{\max}\) and \(T_{\min}\)).

  2. Episodic Precipitation: Precipitation is highly skewed and sparse (zero-inflated). Traditional neural networks with mean squared error struggle with extreme positive values unless properly scaled.

  3. Atmospheric Inertia & Cross-Correlation: Cold/warm fronts persist over multi-day windows. Today’s temperature is strongly correlated with yesterday’s, but non-linear shifts occur when precipitation fronts pass through.

[3]:

# 3. Exploratory Data Analysis: Time Series Dynamics and Distribution Skewness
fig = plt.figure(figsize=(14, 9))
gs = fig.add_gridspec(3, 2, height_ratios=[1.2, 1.0, 1.0])

# Subplot 1: Temperature Envelope over 4 Years
ax1 = fig.add_subplot(gs[0, :])
ax1.plot(df_raw['date'], df_raw['temp_max'], color='#d95f02', lw=1.2, label='Max Temperature ($T_{max}$)')
ax1.plot(df_raw['date'], df_raw['temp_min'], color='#1b9e77', lw=1.0, alpha=0.8, label='Min Temperature ($T_{min}$)')
ax1.fill_between(df_raw['date'], df_raw['temp_min'], df_raw['temp_max'], color='#d95f02', alpha=0.15, label='Diurnal Thermal Range')
ax1.set_title("A. Four-Year Daily Thermal Envelope (Seattle NOAA Station)", fontweight='bold', loc='left')
ax1.set_ylabel("Temperature (deg C)")
ax1.legend(loc='upper right', framealpha=0.9)

# Subplot 2: Daily Precipitation Hyetograph
ax2 = fig.add_subplot(gs[1, :])
ax2.bar(df_raw['date'], df_raw['precipitation'], color='#386cb0', width=1.0, alpha=0.85, label='Precipitation (mm/day)')
ax2.set_title("B. Daily Precipitation Hyetograph (Episodic Storms vs. Dry Summer Intervals)", fontweight='bold', loc='left')
ax2.set_ylabel("Precipitation (mm)")
ax2.legend(loc='upper right', framealpha=0.9)

# Subplot 3: Raw Precipitation Distribution (Heavy-Tailed Skewness)
ax3 = fig.add_subplot(gs[2, 0])
ax3.hist(df_raw['precipitation'], bins=35, color='#386cb0', edgecolor='black', alpha=0.7)
ax3.set_yscale('log')
ax3.set_title("C. Raw Precipitation Distribution (Log-Scale Counts)", fontweight='bold', loc='left')
ax3.set_xlabel("Precipitation (mm)")
ax3.set_ylabel("Count (log)")
zero_pct = ((df_raw['precipitation'] == 0).sum() / len(df_raw)) * 100
ax3.annotate(f"Zero-Rain Days: {(df_raw['precipitation'] == 0).sum()} / {len(df_raw)} ({zero_pct:.1f}%)",
             xy=(15, 50), fontsize=9, bbox=dict(boxstyle="round,pad=0.3", fc="yellow", alpha=0.3))

# Subplot 4: Log-Transformed Precipitation: log(1 + P)
ax4 = fig.add_subplot(gs[2, 1])
log_precip = np.log1p(df_raw['precipitation'])
ax4.hist(log_precip, bins=35, color='#7570b3', edgecolor='black', alpha=0.7)
ax4.set_title("D. Transformed Precipitation: $\\ln(1 + P)$", fontweight='bold', loc='left')
ax4.set_xlabel("$\\ln(1 + \\text{Precipitation})$")
ax4.set_ylabel("Count")

plt.tight_layout()
plt.show()
../_images/CASESTUDY_RNN_LSTM_Class6_2026_5_0.png

3. The Golden Rule of Earth Time Series: Zero Data Leakage

In standard tabular machine learning (e.g., classifying rock specimens or predicting soil properties), random shuffling (train_test_split(shuffle=True)) is standard practice.

[WARNING] Critical Pitfall: Random Shuffling on Time Series

In environmental and Earth observation time series, random shuffling is completely invalid.
If day \(t\) is in the training set and day \(t+1\) is in the test set, the model simply interpolates between neighboring days due to strong temporal autocorrelation. This results in unrealistically low validation loss that catastrophically collapses in real-world deployment.

Chronological Block Split

We enforce a strict temporal boundary:

  • Training Set: 2012-01-01 to 2014-12-31 (3 full calendar years = 1,096 days, ~75% of data).

  • Test Set: 2015-01-01 to 2015-12-31 (1 full calendar year = 365 days, ~25% of data).

  • The model must forecast the unseen future of 2015 using exclusively patterns learned from the past (2012-2014).

[4]:

# 4. Strict Chronological Block Splitting
SPLIT_DATE = pd.to_datetime('2015-01-01')

train_mask = df_raw['date'] < SPLIT_DATE
test_mask = df_raw['date'] >= SPLIT_DATE

df_train = df_raw[train_mask].copy().reset_index(drop=True)
df_test = df_raw[test_mask].copy().reset_index(drop=True)

print(f"[INFO] Training Period: {df_train['date'].min().strftime('%Y-%m-%d')} to {df_train['date'].max().strftime('%Y-%m-%d')} ({len(df_train)} days)")
print(f"[INFO] Test Period:     {df_test['date'].min().strftime('%Y-%m-%d')} to {df_test['date'].max().strftime('%Y-%m-%d')} ({len(df_test)} days)")

# Visualize the temporal split
plt.figure(figsize=(13, 3.5))
plt.plot(df_train['date'], df_train['temp_max'], color='#2b83ba', label='Training Set (2012-2014, 1096 days)', lw=1.2)
plt.plot(df_test['date'], df_test['temp_max'], color='#d7191c', label='Test Set (2015, 365 days)', lw=1.2)
plt.axvline(SPLIT_DATE, color='black', linestyle='--', lw=1.5, label='Temporal Boundary (No Leakage)')
plt.title("Chronological Train / Test Partitioning on Maximum Temperature ($T_{max}$)", fontweight='bold')
plt.ylabel("Temperature (deg C)")
plt.legend(loc='upper right', framealpha=0.9)
plt.tight_layout()
plt.show()

[INFO] Training Period: 2012-01-01 to 2014-12-31 (1096 days)
[INFO] Test Period:     2015-01-01 to 2015-12-31 (365 days)
../_images/CASESTUDY_RNN_LSTM_Class6_2026_7_1.png

4. Geoscientific Feature Preprocessing & Domain Scaling

Recurrent neural networks compute affine transformations followed by non-linear activations (\(\tanh\) and \(\sigma\)). Unscaled variables with wildly different physical units (e.g., wind in m/s vs temperature in deg C vs precipitation in mm) cause gradient oscillations.

Transformations Applied:

  1. Precipitation: As demonstrated in EDA, precipitation is zero-inflated and positively skewed. We compress its dynamic range using the natural logarithmic transform:

    \[P_{\log} = \ln(1 + P)\]

  2. Z-Score Normalization:

    \[z = \frac{x - \mu_{\text{train}}}{\sigma_{\text{train}}}\]

    Crucial Rule: The scaler parameters (mean \(\mu\) and standard deviation \(\sigma\)) are calculated strictly on the training partition and then applied to the test partition. Calculating \(\mu\) on the entire dataset would constitute subtle future data leakage!

[5]:

# 5. Geoscientific Feature Engineering & Leak-Free Normalization
FEATURE_COLS = ['precip_log', 'temp_max', 'temp_min', 'wind']
TARGET_COL = 'temp_max'

# Apply log1p transform to precipitation on the full dataframe
df_raw['precip_log'] = np.log1p(df_raw['precipitation'])

# Extract raw numpy feature matrices for train and test
train_raw_feats = df_raw.loc[train_mask, FEATURE_COLS].values
test_raw_feats = df_raw.loc[test_mask, FEATURE_COLS].values

# Fit StandardScaler STRICTLY on the training set
scaler = StandardScaler()
train_scaled_feats = scaler.fit_transform(train_raw_feats)
test_scaled_feats = scaler.transform(test_raw_feats)

# Store target mean and std for exact inverse transformation during evaluation
target_idx = FEATURE_COLS.index(TARGET_COL)
TARGET_MEAN = scaler.mean_[target_idx]
TARGET_STD = scaler.scale_[target_idx]

print("[INFO] Feature Scaler Parameters (Fitted on Training Set Only):")
for col_name, m, s in zip(FEATURE_COLS, scaler.mean_, scaler.scale_):
    print(f"  - {col_name:12s}: Mean = {m:7.3f}, Std = {s:7.3f}")

print()
print(f"[INFO] Target Variable '{TARGET_COL}' has unscaling equation: y_celsius = (y_scaled * {TARGET_STD:.3f}) + {TARGET_MEAN:.3f}")

[INFO] Feature Scaler Parameters (Fitted on Training Set Only):
  - precip_log  : Mean =   0.712, Std =   1.014
  - temp_max    : Mean =  16.110, Std =   7.330
  - temp_min    : Mean =   8.035, Std =   5.073
  - wind        : Mean =   3.268, Std =   1.471

[INFO] Target Variable 'temp_max' has unscaling equation: y_celsius = (y_scaled * 7.330) + 16.110

5. 3D Sequence Tensor Windowing in PyTorch

In feedforward networks, tabular data has shape (Batch_Size, Num_Features).
In convolutional networks, spatial images have shape (Batch_Size, Channels, Height, Width).
In recurrent networks (RNN, LSTM, GRU), time series tensors require a 3D contract:
\[\mathbf{X} \in \mathbb{R}^{\text{Batch Size} \times \text{Sequence Length} \times \text{Number of Features}}\]
                      Sequence Length (L = 30 days)
               [ t-29,  t-28,  t-27,  ...,   t   ]
Feature 0 (P)   [ 0.1,   0.0,   1.2,  ...,  0.0  ]  --> Predict Target:
Feature 1 (Tmax)[ 12.3, 13.1,  11.8,  ...,  14.2 ]      Tmax at day t+1
Feature 2 (Tmin)[  5.1,  6.0,   4.5,  ...,   7.1 ]
Feature 3 (Wind)[  3.2,  2.8,   4.1,  ...,   3.0 ]

Warm-up History for the Test Year

If we test on 2015, we need 30 days of prior history to predict day 1 (2015-01-01). We concatenate the final \(L=30\) days of 2014 as the warm-up sequence for the test set. This ensures we evaluate all 365 days of 2015 without throwing away the first month!

[6]:

# 6. PyTorch Dataset with Sliding Windows
SEQ_LEN = 30   # 30-day lookback window
BATCH_SIZE = 32

class WeatherSequenceDataset(Dataset):
    # Custom PyTorch Dataset for sliding-window sequence modeling
    def __init__(self, feature_array, target_idx, seq_len=30):
        self.X, self.y = [], []
        n_samples = len(feature_array) - seq_len
        for i in range(n_samples):
            # Sequence from day i to day i + seq_len
            self.X.append(feature_array[i : i + seq_len])
            # Target is the temp_max at day i + seq_len (next day)
            self.y.append(feature_array[i + seq_len, target_idx])

        self.X = torch.tensor(np.array(self.X), dtype=torch.float32)
        self.y = torch.tensor(np.array(self.y), dtype=torch.float32).unsqueeze(1)

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

    def __getitem__(self, idx):
        return self.X[idx], self.y[idx]

# Create Training Dataset
train_dataset = WeatherSequenceDataset(train_scaled_feats, target_idx=target_idx, seq_len=SEQ_LEN)
train_loader = DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True, drop_last=False)

# Prepare Test Set with 30-day warm-up from late 2014
test_warmup_feats = np.vstack([train_scaled_feats[-SEQ_LEN:], test_scaled_feats])
test_dataset = WeatherSequenceDataset(test_warmup_feats, target_idx=target_idx, seq_len=SEQ_LEN)
test_loader = DataLoader(test_dataset, batch_size=BATCH_SIZE, shuffle=False)

print(f"[INFO] Train Dataset: {len(train_dataset)} sliding windows of shape (L={SEQ_LEN}, D={len(FEATURE_COLS)})")
print(f"[INFO] Test Dataset:  {len(test_dataset)} sliding windows (Exactly 365 daily forecasts for year 2015)")

# Verify sample batch dimensions
sample_X, sample_y = next(iter(train_loader))
print(f"[INFO] Batch Input Tensor Shape:  {sample_X.shape} -> (Batch, Seq_Len, Features)")
print(f"[INFO] Batch Target Tensor Shape: {sample_y.shape} -> (Batch, 1)")

[INFO] Train Dataset: 1066 sliding windows of shape (L=30, D=4)
[INFO] Test Dataset:  365 sliding windows (Exactly 365 daily forecasts for year 2015)
[INFO] Batch Input Tensor Shape:  torch.Size([32, 30, 4]) -> (Batch, Seq_Len, Features)
[INFO] Batch Target Tensor Shape: torch.Size([32, 1]) -> (Batch, 1)

6. Part 1: Vanilla RNN & The Vanishing Gradient Pathology

A standard Vanilla Recurrent Neural Network processes input vectors \(\mathbf{x}_t\) sequentially using a recurrent recurrence relation:

\[\mathbf{h}_t = \tanh\left(\mathbf{W}_{ih} \mathbf{x}_t + \mathbf{W}_{hh} \mathbf{h}_{t-1} + \mathbf{b}_h\right)\]
\[\hat{y}_t = \mathbf{W}_{hy} \mathbf{h}_t + b_y\]

The Mathematical Pathology: Backpropagation Through Time (BPTT)

To compute the gradient of the loss \(\mathcal{L}_T\) at time step \(T\) with respect to the initial hidden state \(\mathbf{h}_1\), we apply the multivariable chain rule across all intermediate time steps:

\[\frac{\partial \mathcal{L}_T}{\partial \mathbf{h}_1} = \frac{\partial \mathcal{L}_T}{\partial \mathbf{h}_T} \prod_{k=2}^T \frac{\partial \mathbf{h}_k}{\partial \mathbf{h}_{k-1}} = \frac{\partial \mathcal{L}_T}{\partial \mathbf{h}_T} \prod_{k=2}^T \text{diag}\left(1 - \mathbf{h}_k^2\right) \mathbf{W}_{hh}^T\]

Notice the product of Jacobians:

  1. Vanishing Gradients: The derivative of \(\tanh(z)\) is \(1 - \tanh^2(z) \in (0, 1]\). Multiplying numbers smaller than 1 across \(T=30\) steps results in exponential decay: \(0.5^{30} \approx 9.3 \times 10^{-10}\). The network forgets early sequence signals and cannot learn long-term seasonal trends!

  2. Exploding Gradients: If the eigenvalues of \(\mathbf{W}_{hh}\) exceed 1, the product explodes toward infinity: \(\|\mathbf{W}_{hh}\|^T \to \infty\), causing numerical overflow (NaN loss) and weight divergence.

[7]:

# 7. Vanilla RNN Model Definition with Gradient Norm Logging
class VanillaRNNModel(nn.Module):
    def __init__(self, input_size, hidden_size=64, num_layers=2):
        super(VanillaRNNModel, self).__init__()
        self.hidden_size = hidden_size
        self.num_layers = num_layers
        # nn.RNN with batch_first=True accepts (Batch, Seq_Len, Features)
        self.rnn = nn.RNN(input_size, hidden_size, num_layers=num_layers,
                          nonlinearity='tanh', batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        # rnn_out: (Batch, Seq_Len, Hidden)
        # h_n: (Num_Layers, Batch, Hidden)
        rnn_out, _ = self.rnn(x)
        # We take the representation at the final time step t=L
        last_hidden = rnn_out[:, -1, :]
        out = self.fc(last_hidden)
        return out

# Instantiate Vanilla RNN
rnn_model = VanillaRNNModel(input_size=len(FEATURE_COLS), hidden_size=64, num_layers=2).to(device)
criterion = nn.MSELoss()
optimizer_rnn = optim.Adam(rnn_model.parameters(), lr=0.003)

# Training loop logging gradient norms of recurrent weights
NUM_EPOCHS = 40
rnn_train_losses = []
rnn_grad_norms = []

print("[INFO] Training Vanilla RNN over 40 Epochs...")
for epoch in range(1, NUM_EPOCHS + 1):
    rnn_model.train()
    running_loss = 0.0
    epoch_grad_norms = []

    for X_batch, y_batch in train_loader:
        X_batch, y_batch = X_batch.to(device), y_batch.to(device)
        optimizer_rnn.zero_grad()

        preds = rnn_model(X_batch)
        loss = criterion(preds, y_batch)
        loss.backward()

        # Calculate total gradient norm across all model parameters
        total_norm = 0.0
        for p in rnn_model.parameters():
            if p.grad is not None:
                total_norm += p.grad.data.norm(2).item() ** 2
        total_norm = math.sqrt(total_norm)
        epoch_grad_norms.append(total_norm)

        optimizer_rnn.step()
        running_loss += loss.item() * len(X_batch)

    avg_loss = running_loss / len(train_dataset)
    rnn_train_losses.append(avg_loss)
    rnn_grad_norms.append(np.mean(epoch_grad_norms))

    if epoch % 10 == 0 or epoch == 1:
        print(f"  Epoch {epoch:2d}/{NUM_EPOCHS} | Train MSE Loss: {avg_loss:.4f} | Avg Grad Norm: {rnn_grad_norms[-1]:.4f}")

[INFO] Training Vanilla RNN over 40 Epochs...
  Epoch  1/40 | Train MSE Loss: 0.2686 | Avg Grad Norm: 1.1827
  Epoch 10/40 | Train MSE Loss: 0.1435 | Avg Grad Norm: 0.5475
  Epoch 20/40 | Train MSE Loss: 0.1411 | Avg Grad Norm: 0.6151
  Epoch 30/40 | Train MSE Loss: 0.1343 | Avg Grad Norm: 0.6482
  Epoch 40/40 | Train MSE Loss: 0.1194 | Avg Grad Norm: 0.6441

7. Defensive Engineering: Gradient Norm Clipping

When training recurrent models on multi-variate environmental sequences, sharp sudden events (e.g. violent rainstorms or temperature drops) can induce massive gradient spikes:

\[\mathbf{g} = \nabla_\theta \mathcal{L}\]

If \(\|\mathbf{g}\| > \gamma\), where \(\gamma\) is a maximum threshold (commonly \(\gamma = 1.0\)), gradient clipping rescales the gradient vector back to length \(\gamma\) while preserving its directional trajectory:

\[\mathbf{g} \leftarrow \mathbf{g} \cdot \frac{\gamma}{\max(\gamma, \|\mathbf{g}\|)}\]

In PyTorch, this is implemented cleanly via:

torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)

This single line is essential defensive engineering for sequence modeling, preventing NaN gradient crashes.

[8]:

# 8. Visualizing Vanilla RNN Training Loss & Gradient Norm Dynamics
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))

ax1.plot(range(1, NUM_EPOCHS + 1), rnn_train_losses, marker='o', color='#e41a1c', lw=1.5)
ax1.set_title("Vanilla RNN Training Loss (MSE)", fontweight='bold')
ax1.set_xlabel("Epoch")
ax1.set_ylabel("Loss")

ax2.plot(range(1, NUM_EPOCHS + 1), rnn_grad_norms, marker='s', color='#377eb8', lw=1.5)
ax2.axhline(1.0, color='gray', linestyle='--', label='Standard Clipping Threshold ($\\gamma = 1.0$)')
ax2.set_title("Vanilla RNN Gradient Norm $\\|g\\|_2$ across Epochs", fontweight='bold')
ax2.set_xlabel("Epoch")
ax2.set_ylabel("Gradient L2 Norm")
ax2.legend(loc='upper right')

plt.tight_layout()
plt.show()
../_images/CASESTUDY_RNN_LSTM_Class6_2026_15_0.png

8. Part 2: Long Short-Term Memory (LSTM) Architecture

To solve the vanishing gradient problem, Hochreiter & Schmidhuber (1997) introduced the Long Short-Term Memory (LSTM) network.

The Geoscientific Exception: Spatial & Stratigraphic Sequences

In geology, sequences are frequently spatial rather than temporal:

  • Borehole Geophysical Well Logs: Gamma ray, sonic velocity, and electrical resistivity recorded as a function of depth \(z\).

  • Stratigraphic Core Profiles: Layer facies and grain size measured from bottom to top.

  • Topographic Elevation Transects: Soil moisture transects along a hill slope.

In a borehole log at depth \(z\):

  • The rock layer at depth \(z\) is physically bounded and influenced both by the shallower overburden rocks above (:math:`z - Delta z`) AND the underlying basement formation below (:math:`z + Delta z`)!

  • Modeling spatial depth sequences with a unidirectional RNN is suboptimal because it ignores the forward geological context.

PyTorch Bidirectional LSTM Contract

\[\vec{\mathbf{h}}_z = \text{LSTM}_{\text{forward}}(\mathbf{x}_z, \vec{\mathbf{h}}_{z-1})\]
\[\overleftarrow{\mathbf{h}}_z = \text{LSTM}_{\text{backward}}(\mathbf{x}_z, \overleftarrow{\mathbf{h}}_{z+1})\]
\[\mathbf{h}_z^{\text{Bi}} = \left[ \vec{\mathbf{h}}_z \,\Vert\, \overleftarrow{\mathbf{h}}_z \right] \in \mathbb{R}^{2 \times \text{Hidden}}\]
[14]:

# 14. Bidirectional LSTM for Spatial Well-Logging / Depth Profiles
class BoreholeBiLSTM(nn.Module):
    # Bidirectional LSTM for geological depth sequences (e.g. well-log lithology classification)
    def __init__(self, n_logs, hidden_size=32, n_classes=4):
        super(BoreholeBiLSTM, self).__init__()
        # Setting bidirectional=True creates forward and backward recurrent units
        self.bilstm = nn.LSTM(input_size=n_logs, hidden_size=hidden_size,
                              num_layers=1, batch_first=True, bidirectional=True)
        # Note: hidden representation has dimension 2 * hidden_size!
        self.classifier = nn.Linear(hidden_size * 2, n_classes)

    def forward(self, depth_profiles):
        # depth_profiles: (Batch, Depth_Steps, N_Logs)
        lstm_out, _ = self.bilstm(depth_profiles)
        # lstm_out shape: (Batch, Depth_Steps, 2 * Hidden)
        logits = self.classifier(lstm_out)
        return logits

# Simulate a borehole depth tensor: 4 geophysical logs over 100 meters depth (0.5m spacing = 200 depth steps)
BATCH_WELLS = 4
DEPTH_STEPS = 200
N_GEOPHYSICAL_LOGS = 4  # e.g., Gamma Ray, Resistivity, Neutron Porosity, Density

synthetic_well_logs = torch.randn(BATCH_WELLS, DEPTH_STEPS, N_GEOPHYSICAL_LOGS)
bilstm_net = BoreholeBiLSTM(n_logs=N_GEOPHYSICAL_LOGS, hidden_size=32, n_classes=4)

with torch.no_grad():
    facies_logits = bilstm_net(synthetic_well_logs)

print(f"[INFO] Input Borehole Tensor:  {synthetic_well_logs.shape} -> (Wells, Depth_Steps, Geophysical_Logs)")
print(f"[INFO] Output Lithology Logits: {facies_logits.shape} -> (Wells, Depth_Steps, Lithology_Classes)")
print("[INFO] BiLSTM successfully combined both upward and downward spatial geological context!")

[INFO] Input Borehole Tensor:  torch.Size([4, 200, 4]) -> (Wells, Depth_Steps, Geophysical_Logs)
[INFO] Output Lithology Logits: torch.Size([4, 200, 4]) -> (Wells, Depth_Steps, Lithology_Classes)
[INFO] BiLSTM successfully combined both upward and downward spatial geological context!

14. Hands-on Student Exercises & Research References

[EXERCISE] Practical Challenges for PhD & Post-Doc Researchers

  1. Challenge 1: Multi-Step Autoregressive Rollout

    • Objective: Modify the training loop or inference script to predict a 7-day temperature trajectory rather than a single day (\(t+1, t+2, \dots, t+7\)). Compare an autoregressive iterative rollout (feeding predictions back as inputs) against a direct multi-head output (Linear(hidden_size, 7)).

  2. Challenge 2: Extreme Precipitation Event Classification

    • Objective: Convert the problem into an episodic event detection classifier. Predict whether tomorrow will have heavy rain (\(P > 10\text{ mm}\)) using nn.BCEWithLogitsLoss. Implement class weighting to handle the severe positive class imbalance.

  3. Challenge 3: Memory Length Ablation (:math:`L in {3, 7, 14, 30, 60}`)

    • Objective: Train GeoLSTM with varying sequence history lengths. At what lookback window does test performance plateau? Does a 60-day window improve over a 14-day window for temperature prediction?

Key Takeaways for Geocomputation

  1. Never Shuffle Time Series: Always enforce chronological block splits to prevent optimistic autocorrelation leakage.

  2. The Constant Error Carousel: LSTMs solve vanishing gradients through the additive linear cell state \(\mathbf{c}_t\), allowing gradients to flow backward through dozens of steps.

  3. Gradient Norm Clipping: Essential defensive engineering (clip_grad_norm_) for stabilizing recurrent backpropagation against extreme natural events.

  4. GRU vs. LSTM: GRU achieves comparable predictive skill with 25% fewer parameters, making it ideal for compact geological datasets.

  5. Bidirectional Recurrence: While time series are strictly causal (\(t-1 \to t\)), spatial profiles (boreholes, seismic traces, stratigraphic sections) require bidirectional context.

Landmark Literature & Reading

  • Hochreiter, S., & Schmidhuber, J. (1997). Long Short-Term Memory. Neural Computation, 9(8), 1735-1780.

  • Cho, K., et al. (2014). Learning Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation. EMNLP.

  • Kratzert, F., et al. (2018). Rainfall-runoff modelling using Long Short-Term Memory (LSTM) networks. Hydrology and Earth System Sciences (HESS), 22(11), 6005-6022.

  • Kratzert, F., et al. (2019). Towards learning universal, regional, and local hydrological behaviors via machine learning applied to large-sample datasets. HESS, 23(12), 5083-5110.