Heatwave Severity Forecast — Methodology & Results

Pixel-wise MLP vs U-Net CNN over a 65×41 grid · leak-free 20-variable selection · distributed CPU training (5 nodes × 20 CPUs, PyTorch DDP/Gloo) · 2026-06-11

Workflow

End-to-end pipeline from raw fields to the final report — each stage feeds the next, and the leakage-safe split and training-only feature selection (highlighted) are what keep the evaluation honest.

Variables

All 49 fields share the 65×41 grid. One is the target (obs); two are excluded as leaks; the rest are candidate predictors, of which 20 are selected (4 always-required + 16 chosen by RandomForest importance).

Roles: target = predictand (never an input); leak = excluded (contains the answer, §5); required = always included; selected = chosen into the RF-20 working set; candidate = available but not selected. † marks a full name inferred from the abbreviation and not yet confirmed against dataset metadata — tell me the official definitions and I will correct them.

1. Executive summary

2. Data

3. Experimental design

4. Results at a glance

Held-out test set: 19136 forecasts (final 416 initialization dates × 46 leads, the most recent 20% of the record — never seen in training or selection). 95% confidence intervals come from a block bootstrap over test time steps (1000 resamples).

Δ colored by whether the new model improves on the random/MSE baseline (lower RMSE, higher r & POD, smaller |Bias|).

5. Target-leakage discovery & resistance

During selection, two predictors showed implausible correlation to the target:

• pr: |Pearson r| = 1.000, max|difference| = 0.000 → standardized copy of the target obs.
• sev: |r| ≈ 0.80, fits as 2.5·obs − 3.05 → the raw severity field the target derives from.

Both are now in LEAK_VARS and excluded like obs (46 usable predictors remain). The ablation below trains the same RandomForest with and without them:

6. Dimensionality reduction & variable selection

The 46 leak-free predictors are reduced to N=20 (the four required variables tx, lead, doy_cos, doy_sin are always kept). Three complementary strategies are computed on a memory-bounded subsample of the training period only (no test peeking):

6.1 RandomForest importance (default)

A RandomForest regressor is fit on subsampled valid pixels; variables are ranked by permutation importance — the drop in held-out R² when a variable's values are shuffled. This captures nonlinear, interaction-aware relevance and is the most skill-faithful of the three.

6.2 mRMR filter (minimum-redundancy maximum-relevance)

A greedy filter that trades off relevance to the target against redundancy with already-picked variables:

maximizing correlation to the target y while penalizing correlation to the selected set S. This avoids picking 20 near-duplicate temperature fields.

6.3 PCA (representative-variable dimensionality reduction)

Principal Component Analysis decomposes the standardized predictor matrix X = U Σ Vᵀ; each principal component is a linear combination of predictors capturing an orthogonal mode of variance. We report the cumulative explained-variance scree and, to keep named/interpretable variables, select the predictor with the largest absolute loading on each leading component (representative-variable selection).

6.4 Method agreement

Pairwise Jaccard overlap of the three top-20 sets (1.0 = identical choices):

The two supervised, relevance-driven methods (RandomForest and mRMR) agree well above the chance overlap (~0.28 for two random 20-of-46 subsets), indicating they lock onto the same target-relevant predictors. PCA, which ranks predictors by variance irrespective of the target, overlaps less — exactly as theory predicts, and a useful reminder that unsupervised dimensionality reduction optimizes a different objective than skill. The default rf selection is therefore the skill-faithful choice; mRMR corroborates it.

7. Model architectures

x → Linear(20→256) → BN → GELU
  → [ResBlock(256)] × 4        (skip: x + f(x))
  → BN → GELU → Linear(256→1)

enc1(20→48) ─────────────┐ skip
 ↓pool enc2(48→96) ───┐   │
  ↓pool bottleneck(96→192)
   ↑up  dec2 ←concat──┘   │
    ↑up dec1 ←concat──────┘ → 1×1 conv

8. Loss function — weighted Huber

The baseline MSE produced a systematic negative bias and weak detection of heatwave events. We replace it with a Huber loss (robust to the long-tailed severity distribution) and up-weight event pixels (obs > 0):

Huber is quadratic for |error|≤δ and linear beyond, limiting the influence of extreme residuals; the event weighting (×3) directly targets Probability-of-Detection. The trade-off — a mild positive bias — is discussed in §11.

9. Distributed & parallel training

10. Robustness, reliability & resistance

10.1 Sampling robustness — bootstrap confidence intervals

The 95% CIs in §4 are tight (see the bracket widths), so the skill gain is not an artifact of one particular test split — resampling the test time steps 1000× leaves the ranking and magnitudes stable.

10.2 Where the model is reliable — spatial skill

Skill is spatially coherent, not driven by a few lucky pixels: the CNN's per-pixel correlation is high and uniform across the whole domain (note the CNN RMSE colour-bar tops out well below the MLP's), and its per-pixel bias is close to zero almost everywhere. The MLP shows larger, patchier RMSE — particularly in the more variable southern and coastal pixels — because without spatial context it cannot borrow information from neighbours.

10.3 When the model is reliable — lead time & season

Stratifying by forecast horizon recovers a textbook sub-seasonal skill-decay curve: CNN correlation falls from r ≈ 0.92 at lead 0 (near-analysis) to ≈ 0.85 by lead 5, then plateaus around r ≈ 0.75 from ~2 weeks out. The plateau well above zero indicates the predictors carry genuine slow-varying signal (soil moisture, seasonal state) rather than only short-range persistence. Crucially, the CNN–MLP gap widens with lead: by 6 weeks the CNN holds r ≈ 0.75 while the MLP decays to ≈ 0.55 — spatial context matters most precisely when the local predictor signal is weakest. Skill also varies modestly with season (right panel; the season axis is reconstructed from the standardized day-of-year encodings of the initialization date, so treat it as approximate).

10.4 Error structure & calibration

The error histograms are sharply peaked at zero with a longer positive tail (the event over-prediction from the weighted loss). The CNN's predicted-vs-observed density hugs the 1:1 line but with the classic regression-to-the-mean flattening — extreme observed severities are under-predicted, visible as the negative slope in the residual-vs-observed panel. This is exactly the deficit the linear recalibration in §10.6 corrects, and it explains why a variance-inflating recal slightly improves RMSE.

10.5 Resistance to leakage & selection stability

Beyond removing pr/sev (§5), the selection itself is stable: bootstrapping the RandomForest importance shows the chosen variables are picked consistently across resamples, not by chance.

10.6 Operating point & training-free calibration

The event up-weighting buys a high POD at the cost of a positive bias (over-prediction). Because an additive shift leaves Pearson r unchanged and minimizes RMSE exactly at the mean error, the bias can be removed after training, for free — turning the single trained model into a menu of operating points:

11. Discussion & limitations

12. Reproducibility

python src/feature_selection.py        # writes outputs/selected_vars.json (leak-free ranking)
sbatch slurm.sh                        # 5 nodes × 20 CPUs, DDP; trains MLP + CNN, evaluates
python report/compute_diagnostics.py   # regenerates predictions + robustness diagnostics
python report/selection_analysis.py    # selection stability + leakage ablation
python report/build_report.py          # this report

Config knobs: VAR_SELECT_METHOD (rf/mrmr/pca/random), N_VARS, HEATWAVE_WEIGHT, HUBER_DELTA, LEAK_VARS. Seed = 42.

Generated by report/build_report.py · all figures embedded · open in any browser or print to PDF.