Machine Learning in Engineering Science | NCKU + LAiMM | Fall 2026
Teaching date: Monday 19 October 2026. Guided notebook time: 25 minutes within a 60-minute module.
- Define a sequence, a forecast origin, and a prediction horizon.
- Trace how the four LSTM gates update cell and hidden state.
- Compare a trained LSTM with persistence and ridge baselines without future leakage.
1. In Colab choose Runtime > Change runtime type > CPU. In Jupyter choose a Python 3 kernel.
2. Run cells from top to bottom with Shift+Enter. Read the Markdown before each code cell.
3. The setup installs only missing free packages. A network connection is needed for first installation; the data itself is generated locally.
4. If a variable is undefined, restart the kernel and run all earlier cells. Do not skip a setup cell.
5. Save a personal copy before changing parameters. Every core cell is complete: there are no hidden files, tokens, paid services, or unfinished code blocks.
Data: synthetic teaching data, not engineering measurements. Record your actual output; a model is not required to outperform a baseline.
Practice policy: in-class exercise and optional Homework extension, not a new weekly graded assignment.
Let $u_t=[h_{t-1};x_t]$ concatenate the previous hidden vector and current input. A standard LSTM computes
$$f_t=\sigma(W_fu_t+b_f),\quad i_t=\sigma(W_iu_t+b_i),$$
$$g_t=\tanh(W_gu_t+b_g),\quad o_t=\sigma(W_ou_t+b_o),$$
$$c_t=f_t\odot c_{t-1}+i_t\odot g_t,\quad h_t=o_t\odot\tanh(c_t).$$
Sigmoid gives coefficients between 0 and 1; $\odot$ multiplies matching coordinates. Forget $f$ retains old memory, input $i$ scales a signed write $g$, and output $o$ controls exposed state. Weights are shared across time. Each gate has one coefficient per hidden dimension.
The additive path can improve gradient flow but does not guarantee infinite memory. With gates held fixed, the direct retained contribution includes a product of forget coefficients; full gradients also include gate dependencies. Gradient clipping limits explosions but does not restore vanished gradients.
Run setup once. The random seed fixes the synthetic data; recording package versions makes differences between runs traceable.
import sys, subprocess, importlib.util, platform, json, random, hashlib
from pathlib import Path
required = {"numpy": "numpy", "matplotlib": "matplotlib", "sklearn": "scikit-learn"}
missing = [package for module, package in required.items() if importlib.util.find_spec(module) is None]
if missing:
subprocess.check_call([sys.executable, "-m", "pip", "install", *missing])
import numpy as np
import matplotlib.pyplot as plt
SEED = 2026
random.seed(SEED)
np.random.seed(SEED)
rng = np.random.default_rng(SEED)
ARTIFACTS = Path("artifacts")
ARTIFACTS.mkdir(exist_ok=True)
print("Python:", platform.python_version(), "NumPy:", np.__version__, "Seed:", SEED)
Checkpoint. You should see the Python version, NumPy version, and seed 2026. If installation fails, check internet access and rerun this cell.
No GPU is required. We use one CPU thread and deterministic operations so this small example is inexpensive to rerun.
if importlib.util.find_spec("torch") is None:
subprocess.check_call([sys.executable, "-m", "pip", "install", "torch"])
import torch
from torch import nn
import copy, time
torch.set_num_threads(1)
torch.manual_seed(SEED)
torch.use_deterministic_algorithms(True)
device = torch.device("cpu")
print("PyTorch:", torch.__version__, "Device:", device)
Checkpoint. The printed device must be cpu. Small numerical differences can still occur across package versions or hardware.
Imagine a sensor sampled at equal time intervals. We predict the next sample from the preceding 24 samples. The signal below is synthetic: a slow trend, two oscillations, and noise. Index t represents a sample, not a calibrated second.
N, WINDOW = 1600, 24
t = np.arange(N)
signal = (0.0004*t + np.sin(2*np.pi*t/40) + 0.35*np.sin(2*np.pi*t/13)
+ rng.normal(0, 0.12, N)).astype(np.float32)
TRAIN_END, VAL_END = 960, 1280
plt.figure(figsize=(11,3))
plt.plot(t, signal, lw=1, label="synthetic sensor")
plt.axvline(TRAIN_END,color="tab:orange",label="train/validation boundary")
plt.axvline(VAL_END,color="tab:red",label="validation/test boundary")
plt.xlabel("Sample index"); plt.ylabel("Synthetic sensor units"); plt.legend(); plt.show()
print("Dataset SHA256:", hashlib.sha256(signal.tobytes()).hexdigest())
Checkpoint. Training occurs on earlier samples. Validation and test occur later. Do not shuffle the raw timeline before splitting.
A row is a complete history window, not a single sensor value. With length L=24 and target index t, input is x[t-24:t] and the target is x[t]. Thus every input index is strictly smaller than its target index. A test window may use previously observed validation/test sensor history: this is rolling one-step forecasting, not a forecast of the entire future from one fixed origin.
target_index = np.arange(WINDOW, N)
X_raw = np.stack([signal[i-WINDOW:i] for i in target_index])
y_raw = signal[target_index]
train = target_index < TRAIN_END
valid = (target_index >= TRAIN_END) & (target_index < VAL_END)
test = target_index >= VAL_END
mu, sigma = signal[:TRAIN_END].mean(), signal[:TRAIN_END].std()
X = ((X_raw-mu)/sigma).astype(np.float32)[..., None]
y = ((y_raw-mu)/sigma).astype(np.float32)
assert X.shape == (N-WINDOW, WINDOW, 1)
assert target_index[train].max() < target_index[valid].min()
assert target_index[valid].max() < target_index[test].min()
print("X shape (batch,time,features):", X.shape)
print("Counts:", train.sum(), valid.sum(), test.sum())
print("First window uses indices", target_index[0]-WINDOW, "to", target_index[0]-1,
"and predicts", target_index[0])
Checkpoint. The shape is (1576, 24, 1), with 936 training, 320 validation, and 320 test targets. Mean and standard deviation use only raw training samples.
Persistence predicts that the next value equals the last observation. Ridge regression can exploit the full window without a recurrent state. Choose ridge regularization on validation only, then evaluate the selected baseline on test once.
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_absolute_error, mean_squared_error
ridge_candidates = []
for alpha in [0.01, 0.1, 1.0, 10.0]:
candidate = Ridge(alpha=alpha).fit(X[train,:,0], y[train])
val_mae = mean_absolute_error(y[valid], candidate.predict(X[valid,:,0]))
ridge_candidates.append((val_mae, alpha, candidate))
_, best_alpha, ridge = min(ridge_candidates, key=lambda item:item[0])
predictions = {"Persistence": X[test,-1,0]*sigma+mu,
"Ridge": ridge.predict(X[test,:,0])*sigma+mu}
print("Validation-selected ridge alpha:", best_alpha)
Checkpoint. Persistence has no fitted parameters. Both baselines predict exactly the same target indices as the LSTM.
The forget gate f retains old cell memory; input gate i controls writing; candidate g supplies new content; output gate o exposes the updated cell. All products are elementwise. Here the gate values are chosen, not trained. For c_previous=0.8, f=0.9, i=0.2, g=-0.5, and o=0.7, compute c_new=f*c_previous+i*g before running.
c_previous, f, i, g, o = 0.8, 0.9, 0.2, -0.5, 0.7
c_new = f*c_previous + i*g
h_new = o*np.tanh(c_new)
print("New cell:", c_new, "New hidden:", h_new)
fig, ax = plt.subplots(figsize=(8,3))
ax.bar(["retained old memory", "new written content", "updated cell", "exposed hidden"],
[f*c_previous, i*g, c_new, h_new], color=["#1e8798","#e4002b","#12233f","#f3b23a"])
ax.axhline(0,color="gray"); ax.set_ylabel("Illustrative scalar state"); plt.show()
assert np.isclose(c_new, 0.62)
Checkpoint. The new cell is 0.62 and hidden state is about 0.386. The cell and hidden state are different quantities; an LSTM does not simply copy the last input.
With batch_first=True, input and output are (B,T,F) and (B,T,H). Hidden state has shape (layers,B,H), regardless of batch_first. A linear head turns the final hidden output into one forecast. Model weights are shared across all 24 time steps.
class SensorLSTM(nn.Module):
def __init__(self, hidden=16):
super().__init__()
self.rnn = nn.LSTM(input_size=1, hidden_size=hidden, batch_first=True)
self.head = nn.Linear(hidden, 1)
def forward(self, x):
states, (h_n, c_n) = self.rnn(x)
return self.head(states[:, -1, :]).squeeze(-1)
torch.manual_seed(SEED)
model = SensorLSTM().to(device)
X_t = torch.from_numpy(X)
y_t = torch.from_numpy(y)
print("Parameters:", sum(p.numel() for p in model.parameters()))
print("Output shape:", tuple(model(X_t[:4]).shape))
assert model(X_t[:4]).shape == (4,)
Checkpoint. There must be one prediction per window. A bidirectional model is unnecessary here; every input window already contains only known past values.
Each epoch performs gradient updates on training batches. Validation is evaluated with model.eval() and no gradients. We retain the state with the lowest validation loss. The test set is not consulted during model selection.
optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
loss_fn = nn.MSELoss()
train_ids = np.flatnonzero(train)
history, best_val, best_state = [], float("inf"), None
start = time.perf_counter()
for epoch in range(24):
model.train()
order = rng.permutation(train_ids)
epoch_losses = []
for batch in np.array_split(order, 12):
optimizer.zero_grad()
loss = loss_fn(model(X_t[batch]), y_t[batch])
loss.backward()
nn.utils.clip_grad_norm_(model.parameters(), 1.0)
optimizer.step()
epoch_losses.append(loss.item())
model.eval()
with torch.no_grad():
val_loss = loss_fn(model(X_t[valid]), y_t[valid]).item()
history.append((float(np.mean(epoch_losses)), val_loss))
if val_loss < best_val:
best_val, best_state = val_loss, copy.deepcopy(model.state_dict())
model.load_state_dict(best_state)
print("Training seconds:", round(time.perf_counter()-start,2), "Best validation MSE:", best_val)
plt.plot(np.array(history), label=None)
plt.legend(["training MSE", "validation MSE"]); plt.xlabel("Epoch"); plt.ylabel("Standardized MSE"); plt.show()
Checkpoint. Your actual loss curve is the evidence. If validation worsens while training improves, discuss overfitting rather than increasing epochs blindly.
Report errors in original synthetic sensor units. MAE is the average absolute error; RMSE penalizes large residuals more strongly. A lower score on this toy signal is not evidence of reliability on a real sensor.
model.eval()
with torch.no_grad():
predictions["LSTM"] = model(X_t[test]).numpy()*sigma+mu
metrics = {}
for name, prediction in predictions.items():
mae = mean_absolute_error(y_raw[test], prediction)
rmse = np.sqrt(mean_squared_error(y_raw[test], prediction))
metrics[name] = {"MAE": float(mae), "RMSE": float(rmse)}
print(f"{name:12s} MAE={mae:.4f} RMSE={rmse:.4f}")
fig, axes = plt.subplots(2,1,figsize=(10,5),sharex=True)
axes[0].plot(target_index[test], y_raw[test], label="observed",color="black")
for name, prediction in predictions.items(): axes[0].plot(target_index[test],prediction,label=name,alpha=.7)
axes[0].legend(); axes[0].set_ylabel("Sensor units")
axes[1].plot(target_index[test], predictions["LSTM"]-y_raw[test])
axes[1].axhline(0,color="black"); axes[1].set_ylabel("LSTM residual"); axes[1].set_xlabel("Target sample index")
plt.tight_layout(); plt.show()
Checkpoint. Explain the comparison even if persistence or ridge wins. Inspect whether errors cluster at peaks; one aggregate score hides this structure.
PyTorch stores gate chunks in i,f,g,o order. Reconstruct one cell from learned weights and compare it to nn.LSTM. This connects the equation to an actual implementation, not an illustrative diagram.
x_first = X_t[test][:1,:1,:]
H = model.rnn.hidden_size
h0 = torch.zeros(1,H); c0 = torch.zeros(1,H)
with torch.no_grad():
pre = (x_first[:,0,:] @ model.rnn.weight_ih_l0.T + model.rnn.bias_ih_l0
+ h0 @ model.rnn.weight_hh_l0.T + model.rnn.bias_hh_l0)
pi, pf, pg, po = pre.chunk(4,dim=-1)
gate_i, gate_f, candidate, gate_o = pi.sigmoid(), pf.sigmoid(), pg.tanh(), po.sigmoid()
cell = gate_f*c0 + gate_i*candidate
hidden = gate_o*cell.tanh()
reference, (_, reference_cell) = model.rnn(x_first)
assert torch.allclose(hidden, reference[:,0,:], atol=1e-6)
assert torch.allclose(cell, reference_cell[0], atol=1e-6)
print("Manual learned-cell update matches PyTorch.")
plt.imshow(torch.stack([gate_i[0],gate_f[0],candidate[0],gate_o[0]]).numpy(),aspect="auto",cmap="coolwarm",vmin=-1,vmax=1)
plt.yticks(range(4),["input i","forget f","candidate g","output o"])
plt.xlabel("Hidden dimension"); plt.colorbar(); plt.show()
Checkpoint. The gates are vectors with one value per hidden dimension, not a single universal decision. Their values do not prove causal importance of input samples.
Save data, boundaries, normalization, versions, model state, and measured results. These files stay in your notebook runtime until downloaded. Colab files disappear after the session ends.
record = {"seed":SEED,"window":WINDOW,"train_end":TRAIN_END,"validation_end":VAL_END,
"forecast_protocol":"rolling one-step, observed history available",
"data_sha256":hashlib.sha256(signal.tobytes()).hexdigest(),
"numpy":np.__version__,"torch":torch.__version__,"python":platform.python_version(),
"normalization":{"mean":float(mu),"std":float(sigma)},"metrics":metrics}
(ARTIFACTS/"w07-h1-evidence.json").write_text(json.dumps(record,indent=2))
np.savez(ARTIFACTS/"w07-h1-data.npz",signal=signal,target_index=target_index)
torch.save(model.state_dict(),ARTIFACTS/"w07-h1-lstm.pt")
print(json.dumps(record,indent=2))
Checkpoint. Download the artifacts folder or zip it. A classmate should be able to rebuild the dataset and obtain a comparable baseline from the recorded configuration.
1. Predict how WINDOW=12 versus 48 changes the information available. Rerun from Step 3 with one change only; keep the original result record.
2. Explain why fitting normalization on all 1600 samples leaks future information even though no labels are directly passed to a model.
3. What does the forget gate do when it is 0, 0.5, or 1? Does f near 1 guarantee unlimited memory? Explain the additive path and remaining optimization limits.
4. Why might a linear model beat an LSTM on this signal? State model capacity, data size, and signal structure.
Stop when: you have a held-out comparison table, a residual plot, an evidence JSON, and a two-sentence limitation.
Bridge to H2: an encoder that passes only its final state to a decoder compresses the entire history. Can we retrieve relevant earlier states directly for each output?
- Hochreiter & Schmidhuber (1997), Long Short-Term Memory