4. Perform Sensitivity Analysis on Hyperparameters for a Neural Network

Why Hyperparameter Sensitivity Belongs in the Credibility Case

In a VVUQ workflow, sensitivity analysis answers a specific credibility question: how much of the variability in validation performance is attributable to each hyperparameter, and which ones dominate? This is not the same as hyperparameter tuning. Tuning searches for a good configuration; sensitivity analysis characterizes how the model's validation response behaves across a designed range of configurations so that reviewers can judge whether reported performance is robust or fragile. If a small change in learning rate swings validation accuracy by several points while batch size barely matters, that asymmetry is itself a credibility finding — it tells you where the model is brittle and where uncertainty quantification must concentrate.

This case uses a small, fully specified neural network on a synthetic binary classification task so that every number is reproducible. We run a one-at-a-time (OAT) design to get local, interpretable effects, then a variance-based (Sobol-style) design to apportion output variance across learning rate, batch size, and L2 regularization strength. The two designs answer different questions and are reported together.

Experimental Setup and Fixed Baseline

All runs share the same data, architecture, seed policy, and epoch budget; only the three hyperparameters vary. Fixing everything else is what makes the sensitivity attribution meaningful — if the data split or initialization changed between runs, you would be measuring noise, not hyperparameter effect.

import numpy as np
import torch
import torch.nn as nn
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

# ---- Fixed experimental conditions ----
SEED = 20240501
N_RUNS_PER_CONFIG = 3          # repeat each config over 3 init seeds to average out init noise
EPOCHS = 40
HIDDEN = 32

def make_data():
    X, y = make_classification(
        n_samples=4000, n_features=20, n_informative=8,
        n_redundant=4, class_sep=0.9, flip_y=0.03, random_state=SEED)
    return train_test_split(X, y, test_size=0.3, stratify=y, random_state=SEED)

X_tr, X_va, y_tr, y_va = make_data()
X_tr = torch.tensor(X_tr, dtype=torch.float32); y_tr = torch.tensor(y_tr, dtype=torch.float32)
X_va = torch.tensor(X_va, dtype=torch.float32); y_va = torch.tensor(y_va, dtype=torch.float32)

def build_model():
    return nn.Sequential(
        nn.Linear(20, HIDDEN), nn.ReLU(),
        nn.Linear(HIDDEN, HIDDEN), nn.ReLU(),
        nn.Linear(HIDDEN, 1))

def run_config(lr, batch_size, weight_decay, init_seed):
    torch.manual_seed(init_seed)
    model = build_model()
    opt = torch.optim.Adam(model.parameters(), lr=lr, weight_decay=weight_decay)
    loss_fn = nn.BCEWithLogitsLoss()
    n = X_tr.shape[0]
    g = torch.Generator().manual_seed(init_seed)
    for _ in range(EPOCHS):
        perm = torch.randperm(n, generator=g)
        for i in range(0, n, batch_size):
            idx = perm[i:i+batch_size]
            opt.zero_grad()
            loss = loss_fn(model(X_tr[idx]).squeeze(-1), y_tr[idx])
            loss.backward(); opt.step()
    with torch.no_grad():
        logits = model(X_va).squeeze(-1)
        preds = (torch.sigmoid(logits) >= 0.5).float()
        acc = (preds == y_va).float().mean().item()
    return acc

def evaluate(lr, batch_size, weight_decay):
    accs = [run_config(lr, batch_size, weight_decay, SEED + k) for k in range(N_RUNS_PER_CONFIG)]
    return float(np.mean(accs)), float(np.std(accs))

The output of interest is validation accuracy, averaged over three initialization seeds per configuration. The per-configuration standard deviation is retained because it separates initialization noise from hyperparameter effect: if the spread across seeds is comparable to the spread across hyperparameter values, the sensitivity result is not credible.

One-at-a-Time Design: Local Effects Around the Baseline

OAT holds two hyperparameters at the baseline and sweeps the third across a defined range. It is cheap, easy to explain, and directly answers "what happens if I nudge this one knob?" Its limitation — it cannot detect interactions — is why the variance-based design follows.

baseline = dict(lr=1e-3, batch_size=64, weight_decay=1e-4)
sweeps = {
    "lr":           [1e-4, 3e-4, 1e-3, 3e-3, 1e-2],
    "batch_size":   [16, 32, 64, 128, 256],
    "weight_decay": [0.0, 1e-5, 1e-4, 1e-3, 1e-2],
}

oat_results = {}
for name, values in sweeps.items():
    rows = []
    for v in values:
        cfg = dict(baseline); cfg[name] = v
        mean, sd = evaluate(**cfg)
        rows.append((v, mean, sd))
    oat_results[name] = rows
    print(f"\nOAT sweep: {name}")
    for v, mean, sd in rows:
        print(f"  {name}={v:<8} val_acc={mean:.4f}  (seed sd={sd:.4f})")

Read the OAT output as a local response curve. A steep monotonic segment indicates a hyperparameter the model is sensitive to in that region; a flat segment indicates insensitivity. The seed standard deviation printed alongside each point is the noise floor — any OAT difference smaller than roughly twice that floor should not be interpreted as a real effect.

Variance-Based Design: Apportioning Output Variability

The variance-based design treats each hyperparameter as a random input with a defined distribution and estimates how much of the variance in validation accuracy each one explains. We use a Saltelli-style sampling scheme to estimate first-order Sobol indices (individual contribution) and total-order indices (individual plus interaction contribution).

from scipy.stats import qmc

# Sample in normalized [0,1]^3, then map to physical ranges (log for lr & wd)
BOUNDS = {
    "lr":           (1e-4, 1e-2),
    "batch_size":   (16, 256),
    "weight_decay": (1e-5, 1e-2),
}
D = 3
N_BASE = 64   # base sample size; total model runs = N_BASE*(D+2)*N_RUNS_PER_CONFIG

def to_physical(u):
    lr  = 10 ** (np.log10(BOUNDS["lr"][0]) + u[0]*(np.log10(BOUNDS["lr"][1])-np.log10(BOUNDS["lr"][0])))
    bs  = int(round(2 ** (np.log2(BOUNDS["batch_size"][0]) + u[1]*(np.log2(BOUNDS["batch_size"][1])-np.log2(BOUNDS["batch_size"][0])))))
    wd  = 10 ** (np.log10(BOUNDS["weight_decay"][0]) + u[2]*(np.log10(BOUNDS["weight_decay"][1])-np.log10(BOUNDS["weight_decay"][0])))
    return lr, bs, wd

sampler = qmc.LatinHypercube(d=D, seed=SEED)
A = sampler.random(N_BASE)
B = sampler.random(N_BASE)
AB = []
for i in range(D):
    ABi = A.copy(); ABi[:, i] = B[:, i]
    AB.append(ABi)

def f(u):
    lr, bs, wd = to_physical(u)
    return evaluate(lr, bs, wd)[0]

yA  = np.array([f(u) for u in A])
yB  = np.array([f(u) for u in B])
yAB = [np.array([f(u) for u in ABi]) for i in range(D)]

varY = np.var(np.concatenate([yA, yB]), ddof=1)
names = ["lr", "batch_size", "weight_decay"]
S1, ST = {}, {}
for i, name in enumerate(names):
    S1[name] = float(np.mean(yB * (yAB[i] - yA)) / varY)
    ST[name] = float(0.5 * np.mean((yA - yAB[i])**2) / varY)

print("\nVariance-based sensitivity (Sobol-style):")
for name in names:
    print(f"  {name:<14} S1={S1[name]:+.3f}  ST={ST[name]:+.3f}")

Interpretation rules for the credibility report:

• S1 is the fraction of output variance explained by that hyperparameter alone. Values near 0 mean the hyperparameter is not a dominant driver in the sampled region.
• ST is the fraction explained by that hyperparameter including all its interactions. A large gap between ST and S1 signals strong interaction — the hyperparameter's effect depends on the others.
• If the sum of S1 is well below 1 while the sum of ST is near or above 1, interactions dominate and OAT alone would have been misleading.
• Negative S1 estimates are sampling artifacts and should be reported as approximately zero, not as a real negative contribution.

Reading the Two Designs Together

The two panels are complementary, not redundant. OAT tells you the shape of the local response — where the model is steep and where it is flat — which is what a reviewer needs to judge whether the reported operating point sits on a cliff edge. The variance-based indices tell you how to allocate credibility effort: the hyperparameter with the largest total-order index is the one whose uncertainty most needs to be characterized, and a large ST−S1 gap warns that OAT results alone would understate that hyperparameter's influence.

Deliverable: Sensitivity Summary for the Credibility Report

The following table is the artifact to hand to a reviewer. It records the design, the sampling distributions, and the interpretation, so the sensitivity claim is auditable rather than asserted.

Verification Checks Before Reporting

• Noise floor respected: confirm that OAT differences you call "effects" exceed twice the seed standard deviation; otherwise report them as within noise.
• Sampling documented: the variance-based indices must be accompanied by the exact input distributions, since the indices are conditional on them.
• Index sanity: S1 values should be non-negative after accounting for sampling error; ST should be ≥ S1 for each hyperparameter. Violations indicate too few base samples or a bug in the AB construction.
• Reproducibility: rerunning with a different sampler seed should shift indices by less than the reported uncertainty; a large shift means N_BASE is too small.
• Consistency between designs: if OAT shows a hyperparameter as flat but its ST is large, an interaction is present — investigate before concluding insensitivity.

When these checks pass, the sensitivity result is a defensible credibility input: it identifies which hyperparameters dominate validation variability, quantifies the local response shape around the operating point, and separates genuine hyperparameter effect from initialization noise. That is the evidence a reviewer needs to decide whether the model's reported performance is robust or contingent on a narrow configuration.