RNNs and LSTMs — Sequence Modelling

Reading this sentence word by word — your understanding of each new word depends on everything you have read before. "The bank was steep" versus "The bank was closed" — the word "bank" means something different based on prior context. You carry a mental model forward as you read. That mental model is the hidden state.

An RNN does exactly this — it maintains a hidden state vector that gets updated at every word (or time step). The hidden state at the end of the sequence encodes the full context. The problem: RNNs forget things from 20+ steps ago. LSTMs fix this with explicit memory management using gates.

import numpy as np
import torch
import torch.nn as nn

# ── RNN cell from scratch ─────────────────────────────────────────────
class RNNCellScratch:
    def __init__(self, input_size, hidden_size):
        scale = np.sqrt(1.0 / hidden_size)
        self.Wx = np.random.randn(hidden_size, input_size)  * scale
        self.Wh = np.random.randn(hidden_size, hidden_size) * scale
        self.b  = np.zeros(hidden_size)

    def forward(self, x, h_prev):
        # hₜ = tanh(Wx @ xₜ + Wh @ hₜ₋₁ + b)
        return np.tanh(self.Wx @ x + self.Wh @ h_prev + self.b)

# ── Process a Amazon session sequence ───────────────────────────────
np.random.seed(42)
input_size  = 8    # embedding dimension per action
hidden_size = 16

cell = RNNCellScratch(input_size, hidden_size)

# Simulate: search → product_page → add_to_cart → remove → add_to_cart
session = [np.random.randn(input_size) for _ in range(5)]

h = np.zeros(hidden_size)   # initial hidden state h₀
print("RNN hidden state norm at each step:")
for t, x in enumerate(session):
    h = cell.forward(x, h)
    print(f"  t={t+1}: ||h|| = {np.linalg.norm(h):.4f}  "
          f"range=[{h.min():.3f}, {h.max():.3f}]")

print(f"
Final hidden state (summarises full session): {h[:4].round(4)} ...")

# ── Vanishing gradient demonstration ─────────────────────────────────
print("
Vanishing gradient across time steps (tanh chain):")
grad = 1.0
for t in range(20):
    # tanh derivative ≤ 1, typically ~0.5 in saturated regions
    tanh_deriv = 0.7   # typical value
    grad *= tanh_deriv
    if t % 4 == 3:
        print(f"  After {t+1:2d} steps: gradient = {grad:.6f}")

import numpy as np

# ── LSTM cell from scratch — all four gates explicit ──────────────────
class LSTMCellScratch:
    def __init__(self, input_size, hidden_size):
        self.hs = hidden_size
        scale   = np.sqrt(1.0 / hidden_size)
        # Concatenate [h_prev, x] as input — so weight matrix is (hs, hs+is)
        concat_size = hidden_size + input_size
        self.Wf = np.random.randn(hidden_size, concat_size) * scale
        self.Wi = np.random.randn(hidden_size, concat_size) * scale
        self.Wg = np.random.randn(hidden_size, concat_size) * scale
        self.Wo = np.random.randn(hidden_size, concat_size) * scale
        self.bf = np.zeros(hidden_size)
        self.bi = np.zeros(hidden_size)
        self.bg = np.zeros(hidden_size)
        self.bo = np.zeros(hidden_size)

    def sigmoid(self, z): return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

    def forward(self, x, h_prev, C_prev):
        # Concatenate previous hidden state and current input
        combined = np.concatenate([h_prev, x])

        # ── Three gates ────────────────────────────────────────────────
        f = self.sigmoid(self.Wf @ combined + self.bf)  # forget gate
        i = self.sigmoid(self.Wi @ combined + self.bi)  # input gate
        g = np.tanh(   self.Wg @ combined + self.bg)   # candidate values
        o = self.sigmoid(self.Wo @ combined + self.bo)  # output gate

        # ── Cell state update (additive — key to gradient flow) ────────
        C = f * C_prev + i * g

        # ── New hidden state ───────────────────────────────────────────
        h = o * np.tanh(C)

        return h, C, {'f': f, 'i': i, 'g': g, 'o': o}

# ── Process Amazon session with LSTM ────────────────────────────────
np.random.seed(42)
input_size  = 8
hidden_size = 16

lstm = LSTMCellScratch(input_size, hidden_size)

session_actions = ['home', 'search', 'product', 'add_cart', 'remove', 'add_cart']
session = [np.random.randn(input_size) for _ in session_actions]

h = np.zeros(hidden_size)
C = np.zeros(hidden_size)   # ← cell state: LSTMs have TWO states

print("LSTM states at each step:")
print(f"{'Step':<12} {'Action':<12} {'||h||':>8} {'||C||':>8} {'forget':>8} {'input':>8}")
print("─" * 58)

for t, (x, action) in enumerate(zip(session, session_actions)):
    h, C, gates = lstm.forward(x, h, C)
    print(f"  t={t+1:<9} {action:<12} {np.linalg.norm(h):>8.4f} "
          f"{np.linalg.norm(C):>8.4f} {gates['f'].mean():>8.4f} {gates['i'].mean():>8.4f}")

# ── Compare RNN vs LSTM gradient flow ─────────────────────────────────
print("
Gradient flow comparison (100 time steps):")
rnn_grad  = 1.0
lstm_grad = 1.0
for t in range(100):
    rnn_grad  *= 0.7       # tanh derivative, typical value
    lstm_grad *= 0.95      # additive cell state — much slower decay
    if t % 19 == 19:
        print(f"  Step {t+1:3d}: RNN={rnn_grad:.8f}  LSTM={lstm_grad:.6f}")

import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings('ignore')

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

# ── Shape exploration first — the most important LSTM topic ───────────
lstm = nn.LSTM(
    input_size=8,
    hidden_size=32,
    num_layers=2,
    batch_first=True,    # (batch, seq, features) — matches DataLoader output
    dropout=0.2,
    bidirectional=False,
)

batch, seq_len, input_size = 16, 10, 8
x = torch.randn(batch, seq_len, input_size)   # (batch, seq, input)

output, (h_n, c_n) = lstm(x)

print("nn.LSTM shape guide (batch=16, seq=10, input=8, hidden=32, layers=2):")
print(f"  x shape:      {tuple(x.shape)}")
print(f"  output shape: {tuple(output.shape)}  ← hidden state at EVERY step")
print(f"  h_n shape:    {tuple(h_n.shape)}     ← final hidden (layers, batch, hidden)")
print(f"  c_n shape:    {tuple(c_n.shape)}     ← final cell   (layers, batch, hidden)")
print()
print(f"  output[:, -1, :] shape: {tuple(output[:, -1, :].shape)}  ← last step only")
print(f"  h_n[-1] shape:          {tuple(h_n[-1].shape)}            ← last layer final h")
print("  Both give the same values for single-direction LSTM")

# ── Bidirectional LSTM ─────────────────────────────────────────────────
bilstm = nn.LSTM(input_size=8, hidden_size=32, batch_first=True, bidirectional=True)
out_bi, (h_bi, c_bi) = bilstm(x)
print(f"
Bidirectional LSTM:")
print(f"  output shape: {tuple(out_bi.shape)}  ← hidden_size × 2 = 64")
print(f"  h_n shape:    {tuple(h_bi.shape)}    ← 2 directions × 1 layer")

import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np
from torch.utils.data import Dataset, DataLoader
from sklearn.metrics import roc_auc_score
import copy, warnings
warnings.filterwarnings('ignore')

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

# ── Simulate Amazon browsing sessions ───────────────────────────────
# Each session: variable-length sequence of actions
# Each action: 16-dim feature vector (action type, time on page, etc.)
# Label: 1 = purchased within 30 min, 0 = did not

N_SESSIONS  = 2000
INPUT_SIZE  = 16
MAX_SEQ_LEN = 20

class AmazonSessionDataset(Dataset):
    def __init__(self, n=N_SESSIONS):
        np.random.seed(42)
        self.sessions = []
        self.labels   = []
        for _ in range(n):
            # Buyers have slightly different patterns
            will_buy = np.random.random() < 0.3   # 30% purchase rate
            seq_len  = np.random.randint(3, MAX_SEQ_LEN + 1)
            # Buyers: more time on product pages, more add-to-cart events
            signal   = 0.4 if will_buy else 0.0
            session  = np.random.randn(seq_len, INPUT_SIZE).astype(np.float32)
            session[:, 0] += signal    # feature 0: add-to-cart frequency
            session[:, 1] += signal    # feature 1: time on product page
            self.sessions.append(torch.FloatTensor(session))
            self.labels.append(int(will_buy))

    def __len__(self): return len(self.sessions)
    def __getitem__(self, i): return self.sessions[i], self.labels[i]

def collate_fn(batch):
    """Pad variable-length sequences to the same length within a batch."""
    sequences, labels = zip(*batch)
    lengths = torch.tensor([len(s) for s in sequences])
    # Pad to max length in this batch
    padded = nn.utils.rnn.pad_sequence(sequences, batch_first=True)
    return padded, torch.tensor(labels, dtype=torch.long), lengths

dataset    = AmazonSessionDataset()
train_ds, val_ds = torch.utils.data.random_split(dataset, [1600, 400])
train_ld   = DataLoader(train_ds, batch_size=64, shuffle=True,  collate_fn=collate_fn)
val_ld     = DataLoader(val_ds,   batch_size=64, shuffle=False, collate_fn=collate_fn)

# ── LSTM classifier ───────────────────────────────────────────────────
class SessionLSTM(nn.Module):
    def __init__(self, input_size=INPUT_SIZE, hidden_size=64,
                 num_layers=2, dropout=0.3):
        super().__init__()
        self.lstm = nn.LSTM(
            input_size, hidden_size,
            num_layers=num_layers,
            batch_first=True,
            dropout=dropout if num_layers > 1 else 0,
        )
        self.classifier = nn.Sequential(
            nn.Linear(hidden_size, 32),
            nn.ReLU(),
            nn.Dropout(0.3),
            nn.Linear(32, 1),   # binary: buy or not
        )

    def forward(self, x, lengths):
        # Pack padded sequence — LSTM ignores padding positions
        packed = nn.utils.rnn.pack_padded_sequence(
            x, lengths.cpu(), batch_first=True, enforce_sorted=False,
        )
        _, (h_n, _) = self.lstm(packed)
        # h_n: (num_layers, batch, hidden) — take last layer
        last_h = h_n[-1]          # (batch, hidden_size)
        return self.classifier(last_h).squeeze(1)

model     = SessionLSTM()
criterion = nn.BCEWithLogitsLoss()
optimizer = optim.AdamW(model.parameters(), lr=1e-3, weight_decay=0.01)
scheduler = optim.lr_scheduler.ReduceLROnPlateau(optimizer, patience=5, factor=0.5)

# ── Training loop ─────────────────────────────────────────────────────
best_auc, best_wts, patience_count = 0.0, None, 0

print("Training SessionLSTM on Amazon purchase intent:")
print(f"{'Epoch':>6} {'Train loss':>12} {'Val AUC':>10}")
print("─" * 32)

for epoch in range(1, 31):
    model.train()
    total_loss = 0
    for Xb, yb, lengths in train_ld:
        optimizer.zero_grad()
        logits = model(Xb, lengths)
        loss   = criterion(logits, yb.float())
        loss.backward()
        nn.utils.clip_grad_norm_(model.parameters(), 1.0)  # gradient clipping
        optimizer.step()
        total_loss += loss.item()

    model.eval()
    all_probs, all_labels = [], []
    with torch.no_grad():
        for Xb, yb, lengths in val_ld:
            probs = torch.sigmoid(model(Xb, lengths))
            all_probs.extend(probs.numpy())
            all_labels.extend(yb.numpy())
    val_auc = roc_auc_score(all_labels, all_probs)
    scheduler.step(1 - val_auc)

    if epoch % 5 == 0:
        print(f"  {epoch:>4}  {total_loss/len(train_ld):>12.4f}  {val_auc:>10.4f}")

    if val_auc > best_auc:
        best_auc, best_wts = val_auc, copy.deepcopy(model.state_dict())
        patience_count = 0
    else:
        patience_count += 1
        if patience_count >= 10:
            print(f"  Early stop at epoch {epoch}")
            break

print(f"
Best val AUC: {best_auc:.4f}")

import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np
from sklearn.metrics import mean_absolute_error
import warnings
warnings.filterwarnings('ignore')

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

# ── Simulate Instacart hourly demand data ─────────────────────────────────
# 1 SKU at 1 dark store, hourly demand for 60 days = 1440 hours
hours = np.arange(1440)
# Demand = base + daily pattern + weekly pattern + noise
demand = (
    50
    + 20 * np.sin(2 * np.pi * hours / 24)           # daily cycle
    + 10 * np.sin(2 * np.pi * hours / (24 * 7))     # weekly cycle
    + np.abs(np.random.normal(0, 8, 1440))           # noise
).clip(0, None)

# ── Create sliding window sequences ───────────────────────────────────
SEQ_LEN = 24   # use last 24 hours to predict next hour

X, y = [], []
for i in range(len(demand) - SEQ_LEN):
    X.append(demand[i:i+SEQ_LEN])
    y.append(demand[i+SEQ_LEN])

X = np.array(X, dtype=np.float32)
y = np.array(y, dtype=np.float32)

# Normalise
X_mean, X_std = X.mean(), X.std()
X = (X - X_mean) / X_std
y = (y - X_mean) / X_std

# Train/test split (time-ordered — no shuffle)
split  = int(len(X) * 0.8)
X_tr, X_te = X[:split], X[split:]
y_tr, y_te = y[:split], y[split:]

X_tr_t = torch.FloatTensor(X_tr).unsqueeze(-1)   # (N, seq, 1)
X_te_t = torch.FloatTensor(X_te).unsqueeze(-1)
y_tr_t = torch.FloatTensor(y_tr).unsqueeze(-1)
y_te_t = torch.FloatTensor(y_te).unsqueeze(-1)

loader = torch.utils.data.DataLoader(
    torch.utils.data.TensorDataset(X_tr_t, y_tr_t),
    batch_size=64, shuffle=True,
)

# ── LSTM forecaster ────────────────────────────────────────────────────
class DemandLSTM(nn.Module):
    def __init__(self):
        super().__init__()
        self.lstm = nn.LSTM(
            input_size=1,
            hidden_size=64,
            num_layers=2,
            batch_first=True,
            dropout=0.2,
        )
        self.head = nn.Linear(64, 1)

    def forward(self, x):
        out, _ = self.lstm(x)
        # Use only the last time step's hidden state for prediction
        return self.head(out[:, -1, :])

model     = DemandLSTM()
criterion = nn.MSELoss()
optimizer = optim.AdamW(model.parameters(), lr=1e-3, weight_decay=0.01)

print("Training Instacart demand forecaster (LSTM):")
for epoch in range(1, 51):
    model.train()
    for Xb, yb in loader:
        optimizer.zero_grad()
        criterion(model(Xb), yb).backward()
        optimizer.step()

    if epoch % 10 == 0:
        model.eval()
        with torch.no_grad():
            pred_norm = model(X_te_t).numpy()
        pred = pred_norm * X_std + X_mean
        true = y_te_t.numpy() * X_std + X_mean
        mae  = mean_absolute_error(true, pred)
        print(f"  Epoch {epoch:2d}: MAE = {mae:.2f} units/hour")

# ── Baseline comparison ────────────────────────────────────────────────
# Naive baseline: predict last observed value
naive_pred = X_te[:, -1] * X_std + X_mean
naive_mae  = mean_absolute_error(y_te * X_std + X_mean, naive_pred)
print(f"
Naive baseline MAE: {naive_mae:.2f} units/hour")

import numpy as np

# ── Why a live anomaly detector reaches for an LSTM, not a Transformer ──
# A sensor feed emits a reading every second and never stops. The model
# must score each new reading using history, forever, without unbounded
# memory growth. An LSTM's hidden and cell state do exactly this.

class StreamingLSTMScorer:
    def __init__(self, input_size, hidden_size):
        scale = np.sqrt(1.0 / hidden_size)
        concat = hidden_size + input_size
        self.Wf = np.random.randn(hidden_size, concat) * scale
        self.Wi = np.random.randn(hidden_size, concat) * scale
        self.Wg = np.random.randn(hidden_size, concat) * scale
        self.Wo = np.random.randn(hidden_size, concat) * scale
        self.h  = np.zeros(hidden_size)
        self.C  = np.zeros(hidden_size)     # fixed size — never grows

    def sigmoid(self, z): return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

    def step(self, x_reading):
        combined = np.concatenate([self.h, x_reading])
        f = self.sigmoid(self.Wf @ combined)
        i = self.sigmoid(self.Wi @ combined)
        g = np.tanh(self.Wg @ combined)
        o = self.sigmoid(self.Wo @ combined)
        self.C = f * self.C + i * g
        self.h = o * np.tanh(self.C)
        return np.linalg.norm(self.h)        # anomaly score for this reading

np.random.seed(42)
scorer = StreamingLSTMScorer(input_size=4, hidden_size=16)

# Simulate a sensor stream running for a very long time — memory used
# to hold state (self.h, self.C) never grows, no matter how many
# readings have already been processed.
for reading_number in [1, 100, 10_000, 1_000_000]:
    reading = np.random.randn(4)
    score = scorer.step(reading)
    state_size = scorer.h.nbytes + scorer.C.nbytes
    print(f"  reading #{reading_number:>9,}: score={score:.4f}  "
          f"state size held in memory = {state_size} bytes (constant)")

print("
A Transformer scoring the same stream would need to grow its")
print("attended context (or manage a KV cache) as the stream gets")
print("longer — an LSTM's state size here never changes.")

LSTMs dramatically mitigate vanishing gradients through the additive cell-state update (C = f×C_prev + i×g), which lets gradients flow backward through many time steps without being repeatedly multiplied through a squashing non-linearity — but "mitigate" is not "eliminate." Extremely long sequences can still see gradient decay, since the forget gate itself is a sigmoid that can push toward zero and cut the memory pathway, and exploding gradients remain a live risk through the recurrent weight matrices regardless of gating. That is exactly why gradient clipping is still considered mandatory, standard practice for LSTM training, not an optional safety net — the gates reduce the severity of the problem, they don't make gradient management unnecessary.

Each gate outputs a sigmoid value between 0 and 1 for every dimension of the hidden state, which does make them "gate-like" in spirit, but they are continuously valued and jointly computed from the same shared inputs, not independent binary decisions. A single dimension might have a forget value of 0.3 (mostly, not fully, erase this memory slot) combined with an input value of 0.6 (partially write new information into the same slot) in the same time step — the gates blend old and new information smoothly rather than switching cleanly between "keep everything" and "erase everything." Thinking of them as binary flags misses exactly the mechanism that makes gating powerful: fine-grained, differentiable control over what's remembered.

A GRU merges the LSTM's separate hidden state and cell state into one, and uses two gates (reset and update) instead of three, giving it noticeably fewer parameters and typically faster training and inference per step. Empirically the two often perform comparably on many mid-sized sequence tasks, which is where the "doesn't matter" intuition comes from — but "often comparable" is not "always equivalent." LSTMs' explicit, separate cell-state memory pathway tends to have an edge on tasks with longer-range dependencies, while GRUs' smaller footprint is genuinely preferable when compute- or latency-constrained. Treat the choice as an empirical question for your specific task, not an interchangeable default.

Transformers dominate for tasks where you can afford to process the whole sequence at once and attention's quadratic cost in sequence length is acceptable — but LSTMs remain the better engineering choice in specific, common situations. For streaming or online inference, an LSTM's hidden/cell state gives constant memory and compute per new token, while a Transformer needs to manage a growing context window or KV cache. For very long sequences — tens of thousands of steps, common in genomics or long-form sensor telemetry — full self-attention's quadratic cost can be computationally infeasible where an LSTM's linear cost is not. "Transformers won" for language modelling at scale; that is not a blanket replacement for recurrent architectures in every setting.

A bidirectional LSTM runs one LSTM forward and a second backward over the sequence and concatenates their hidden states — genuinely useful when the full sequence is available upfront and you want each position informed by both what came before and after it. But that requirement makes it structurally unusable for real-time, causal, or streaming prediction: you cannot run the "backward" direction on future tokens that haven't happened yet. For online tasks like next-token prediction, live user-session scoring, or step-by-step time-series forecasting, only a unidirectional LSTM is a valid architecture at all, regardless of any accuracy benefit bidirectionality might offer on paper.

A vanilla RNN's hidden state is recomputed every step through a tanh non-linearity — backpropagating through time means repeatedly multiplying by tanh's derivative (bounded well under 1 in saturated regions) at every step, so the gradient shrinks exponentially with sequence length. The LSTM's key structural change is the cell state, updated additively — C = f×C_prev + i×g — rather than being pushed back through a squashing function at every step. Because addition, unlike repeated multiplication by sub-1 values, doesn't inherently shrink a gradient, gradients can flow backward through many time steps along the cell-state pathway largely intact, as long as the forget gate stays close to 1 for the relevant span. The gates decide when to use that additive pathway versus when to actually forget — but the additive update itself is what solves the underlying math problem.

With batch_first=True, input is (batch, seq_len, input_size) — matching what a standard DataLoader naturally produces. The call returns output, (h_n, c_n): output has shape (batch, seq_len, hidden_size) and contains the hidden state at every single time step. h_n and c_n, by contrast, are (num_layers, batch, hidden_size) — batch, not seq_len, is not first here even with batch_first=True, because these represent only the final states, indexed by layer rather than by time step. The two most common bugs are forgetting that h_n's first dimension is num_layers, so h_n[-1] gives the last layer's final hidden state, and conflating output[:, -1, :] with h_n[-1] — equal for a single-direction, single-layer LSTM, but they diverge the moment you add layers or bidirectionality.

Three cases come up in practice. First, streaming or real-time inference where you need to process a new data point as it arrives with constant per-step latency and memory — an LSTM's hidden/cell state naturally supports this, while a Transformer needs a growing context window or KV cache. Second, very long sequences — tens of thousands of time steps, common in genomics or extended time series — where full self-attention's quadratic cost becomes computationally prohibitive and an LSTM's linear cost is the pragmatic choice. Third, resource-constrained deployment, where a small LSTM is simply cheaper to run than even a modest Transformer. For most NLP and any task with abundant compute and the full sequence available upfront, Transformers remain the default.

A feedforward network's gradient flows through depth — a fixed, architecturally-bounded number of layers. A recurrent network's gradient flows through both depth and time: processing a sequence of length 100 is, from backprop's perspective, equivalent to a 100-layer-deep unrolled network, and that unrolled depth is determined by your data, not a fixed architectural choice. This makes exploding gradients — where the same recurrent weight matrix, applied repeatedly, amplifies the gradient each step — a much more common and less predictable failure mode for RNNs/LSTMs than for a network with a handful of fixed layers. Gradient clipping rescales the gradient vector when its norm exceeds a threshold, cheap insurance applied by default in essentially every serious RNN/LSTM training loop, independent of whether NaNs have actually been observed yet.

For a single-layer, unidirectional LSTM, h_n[-1] and output[:, -1, :] contain identical values — both are the hidden state at the final time step. They diverge the moment you add multiple layers or bidirectionality. With multiple layers, output only ever reflects the top layer's hidden state at each time step, while h_n contains the final hidden state of every layer stacked — so h_n[-1] is what matches output[:, -1, :], and grabbing the wrong index silently gives you an intermediate layer's representation instead. With bidirectionality it's worse: h_n's last two entries are the final forward-direction and final backward-direction states, and the backward direction's "final" state actually corresponds to the first time step of the sequence — using output[:, -1, :] directly in that setting mixes a real forward-final-step signal with the backward pass's very first computed state, a substantive and easy-to-miss correctness bug.