Scalable Uncertainty Quantification for Extreme Weather Forecasting via Empirical Neural Tangent Kernels

# Scalable Uncertainty Quantification for Extreme Weather Forecasting via Empirical Neural Tangent Kernels
Source: [https://arxiv.org/html/2606.02886](https://arxiv.org/html/2606.02886)
Jose Marie Antonio Miñoza,Rex Gregor LayloDepartment of EducationCenter for AI ResearchMakatiPhilippines[ecair\.rlaylo@deped\.gov\.ph](https://arxiv.org/html/2606.02886v1/mailto:[email protected])andSebastian C\. IbañezDepartment of EducationCenter for AI ResearchMakatiPhilippines[ecair\.sibanez@deped\.gov\.ph](https://arxiv.org/html/2606.02886v1/mailto:[email protected])

\(2026\)

###### Abstract\.

Deep learning weather models now match numerical weather prediction accuracy while running orders of magnitude faster, but produce deterministic forecasts without uncertainty estimates, a critical gap for high\-stakes decisions during extreme weather events\. This paper proposes Neural Tangent Kernel\-based uncertainty quantification \(NTK\-UQ\) using last\-layer empirical features\. Theoretical analysis predicts that UQ quality is architecture\-dependent through two mechanisms\. First, a variance collapse mechanism explains when UQ fails: when the eigenvalue truncation rank approaches the effective rank of the feature space, the GP correction term consumes nearly all prior variance, destroying discrimination between tropical cyclones and routine conditions; architectures with concentrated spectra \(spectral operators\) require aggressive truncation \(k≤10k\\leq 10\), while attention\-based models tolerate full\-rank computation\. Second, decomposition performance depends on the non\-Gaussian, heavy\-tailed structure of extreme weather: Independent Component Analysis exploits higher\-order statistics \(kurtosis, negentropy\) to isolate heavy\-tailed extreme\-event features, achieving higher discrimination than singular value decomposition, which captures only second\-order variance\. A data\-driven selection rule chooses ICA or SVD from the feature eigenspectrum concentration ratio, correctly prescribing the superior decomposition for all four evaluated architectures\. Compared to split conformal prediction \(the natural post\-hoc baseline\), NTK\-UQ achieves 31–37% sharper prediction intervals at 90% coverage, and uniquely produces*adaptive*intervals that scale with extreme event severity, which conformal prediction cannot achieve by construction\. The framework requires no retraining; inference\-time uncertainty requires only a single matrix\-vector product per sample\.

uncertainty quantification, neural tangent kernel, Gaussian processes, deep learning, calibration, weather forecasting

††journalyear:2026††copyright:othergov††conference:Proceedings of the 32nd ACM SIGKDD Conference on Knowledge Discovery and Data Mining V\.2; August 09–13, 2026; Jeju Island, Republic of Korea††booktitle:Proceedings of the 32nd ACM SIGKDD Conference on Knowledge Discovery and Data Mining V\.2 \(KDD ’26\), August 09–13, 2026, Jeju Island, Republic of Korea††doi:10\.1145/3770855\.3818106††isbn:979\-8\-4007\-2259\-2/2026/08††ccs:Computing methodologies Uncertainty quantification††ccs:Computing methodologies Gaussian processes††ccs:Computing methodologies Neural networks††ccs:Computing methodologies Spectral methods††ccs:Applied computing Earth and atmospheric sciences![Refer to caption](https://arxiv.org/html/2606.02886v1/figures/pipeline.png)Figure 1\.Overview of the NTK\-UQ pipeline for extreme weather forecasting\. Atmospheric variables from extreme weather events are processed through four foundation AI weather models \(FourCastNetV2, Aurora, AIFS, Pangu\-Weather\) to extract last\-layer features\. These features construct the empirical Neural Tangent Kernel matrix, which is decomposed via SVD or ICA \(shown:U​Σ​V⊤U\\Sigma V^\{\\top\}decomposition\) to obtain rank\-kkapproximation\. At inference, the GP posterior variance formula yields calibrated prediction intervals that quantify epistemic uncertainty per variable\.Diagram showing the NTK\-UQ pipeline: feature extraction from four AI weather models, empirical NTK kernel construction, SVD/ICA decomposition, and GP posterior uncertainty estimation at inference\.## 1\.Introduction

Extreme weather events cause an estimated US$143 billion per year in climate\-attributable damages\(Newman and Noy,[2023](https://arxiv.org/html/2606.02886#bib.bib40)\), with the EM\-DAT database recording 399 disasters in 2023 alone, affecting 93\.1 million people\(Delforgeet al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib44)\)\. Accurate forecasting of these events is essential, yet the value of a forecast depends not only on its accuracy but on knowing*how much to trust it*, a question that requires calibrated uncertainty estimates\.

Deep learning has transformed weather forecasting\. Models such as FourCastNetV2\(Pathaket al\.,[2022](https://arxiv.org/html/2606.02886#bib.bib4)\), Pangu\-Weather\(Biet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib6)\), GraphCast\(Lamet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib7)\), and Aurora\(Bodnaret al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib8)\)now match or exceed the accuracy of traditional numerical weather prediction \(NWP\) systems while running orders of magnitude faster, generating 10\-day global forecasts in seconds rather than hours\. However, these models produce deterministic point forecasts without calibrated uncertainty estimates\. Uncertainty quantification is essential for risk\-sensitive applications: decision\-makers require not only point predictions but probabilistic intervals that correlate with actual forecast errors\. For extreme events where forecast errors carry the highest consequences, the absence of reliable uncertainty estimates limits model utility\.

Existing approaches to uncertainty quantification \(UQ\) for neural networks face significant limitations when applied to large\-scale weather models\. Deep ensembles\(Lakshminarayananet al\.,[2017](https://arxiv.org/html/2606.02886#bib.bib13)\)require training multiple copies of billion\-parameter models from scratch, which is computationally prohibitive for foundation weather models\. Monte Carlo dropout\(Gal and Ghahramani,[2016](https://arxiv.org/html/2606.02886#bib.bib17)\)can produce miscalibrated uncertainties\(Ovadiaet al\.,[2019](https://arxiv.org/html/2606.02886#bib.bib26)\)and requires architectural modifications incompatible with pre\-trained checkpoints\. Bayesian neural networks\(Blundellet al\.,[2015](https://arxiv.org/html/2606.02886#bib.bib18)\)add substantial memory and compute overhead, scaling poorly to operational\-size models\. Conformal prediction\(Angelopoulos and Bates,[2021](https://arxiv.org/html/2606.02886#bib.bib19)\)provides distribution\-free coverage guarantees but in its standard form produces uniform interval widths that do not correlate with actual forecast errors\.

This paper proposes last\-layer Neural Tangent Kernel \(NTK\) based uncertainty quantification for AI weather models\. The key insight is that a weather model’s last\-layer featuresϕ​\(x\)\\phi\(x\)—learned from decades of ERA5 reanalysis—encode physically meaningful atmospheric structure\. Under the last\-layer NTK–GP correspondence, the feature kernelK​\(x,x′\)=ϕ​\(x\)⊤​ϕ​\(x′\)K\(x,x^\{\\prime\}\)=\\phi\(x\)^\{\\top\}\\phi\(x^\{\\prime\}\)acts as an*ERA5\-informed similarity measure*: a test input receives high uncertainty when its atmospheric state is unusual relative to both the model’s learned feature manifold and the calibration distribution\. This two\-level epistemic signal is inaccessible to purely statistical baselines such as conformal prediction\. Critically, UQ quality is architecture\-dependent and decomposition\-dependent: a data\-driven selection rule determines whether Independent Component Analysis or Singular Value Decomposition is appropriate from the feature eigenspectrum, correctly prescribing the superior method without exhaustive comparison\.

Throughout this paper, the term*NTK uncertainty*refers to the posterior variance obtained by treating the frozen model’s last\-layer features as an*empirical*Neural Tangent Kernel and applying Gaussian Process posterior theory\. This usage differs from the full infinite\-width NTK formulation and should be interpreted as a finite\-width, post\-hoc kernel approximation induced by the learned feature representations\. The theoretical results in this paper \(Propositions[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)and[2](https://arxiv.org/html/2606.02886#S3.Thmtheorem2), Theorem[1](https://arxiv.org/html/2606.02886#S5.Thmtheorem1)\) are proved directly under this empirical kernel without invoking the infinite\-width limit; prior work\(Heet al\.,[2020](https://arxiv.org/html/2606.02886#bib.bib12); Huanget al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib32)\)shows that finite\-width networks behave approximately as kernel machines, and post\-hoc calibration corrects for residual approximation error\.

NTK\-UQ has several properties that make it suitable for studying UQ across large\-scale weather models\. First, the method requires no model retraining or architectural changes; it works with any pre\-trained checkpoint as a purely post\-hoc procedure\. Second, after one\-time offline calibration, inference\-time UQ requires only a matrix\-vector product, adding minimal overhead to the forward pass\. Third, uncertainties are computed per output variable, enabling variable\-level uncertainty estimates\.

Theoretical analysis predicts that UQ quality depends on both neural architecture \(through eigenspectrum concentration\) and decomposition method \(through higher\-order statistics exploitation\)\. The framework is evaluated on four architecturally diverse AI weather models: FourCastNetV2 \(SFNO\), Pangu\-Weather \(Swin Transformer\), Aurora \(Perceiver\), and AIFS \(GNN\-Transformer\), using ERA5 reanalysis\(Hersbachet al\.,[2020](https://arxiv.org/html/2606.02886#bib.bib10)\)as ground truth\. Evaluation focuses on extreme weather events from the EM\-DAT International Disaster Database, including tropical cyclones, floods, droughts, and extreme temperature events\. Experiments span forecast lead times from 6 to 120 hours\. Results validate these predictions: uncertainty discrimination quality follows architecture\-dependent patterns, with Independent Component Analysis achieving adaptive intervals that scale with extreme event severity, while singular value decomposition produces more uniform intervals that fail to distinguish tropical cyclone forecasts from routine conditions\.

#### Contributions\.

This paper makes five contributions: \(1\)Variance Collapse Characterization: formal analysis of how eigenvalue spectrum concentration causes UQ failure, with diagnostic criterionRk=Ck/P<0\.9R\_\{k\}=C\_\{k\}/P<0\.9for maintaining discrimination \(Proposition[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)\), linking neural architecture \(SFNO vs Transformer\) to effective rank and optimal truncation strategy; \(2\)Non\-Gaussian Discrimination Theory: explanation of why Independent Component Analysis outperforms singular value decomposition for extreme weather through higher\-order statistics exploitation \(Proposition[2](https://arxiv.org/html/2606.02886#S3.Thmtheorem2)\), providing theoretical justification for decomposition method selection based on feature distribution properties; \(3\)Architecture\-UQ Interaction Framework: systematic characterization of how neural architecture determines NTK eigenspectrum properties, which govern UQ quality, enabling predictive diagnosis without exhaustive experimentation; \(4\)Decomposition Selection Rule: Algorithm[1](https://arxiv.org/html/2606.02886#alg1)provides a data\-driven recipe that selects ICA or SVD from the feature eigenspectrum concentration ratio, correctly prescribing the superior method for all four evaluated architectures, validated against split conformal prediction with 31–37% sharper intervals in 81% of valid comparisons; and \(5\)Empirical Validation: evaluation across four foundation weather models \(FourCastNetV2, Pangu\-Weather, Aurora, AIFS\) on 100 extreme weather events from EM\-DAT confirms theoretical predictions and demonstrates that NTK\-UQ produces adaptive intervals \(CV\>0\>0\) that conformal prediction cannot achieve by construction\.

## 2\.Related Work

AI weather foundation models\(Pathaket al\.,[2022](https://arxiv.org/html/2606.02886#bib.bib4); Biet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib6); Lamet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib7); Bodnaret al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib8)\)now match numerical weather prediction accuracy\. A subset produce probabilistic forecasts natively, but each at a cost: ECMWF’s operational ensemble \(ENS\) requires 51\-member perturbation runs at deployment; GenCast\(Priceet al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib1)\)trains a diffusion model from scratch; and SEEDS\(Liet al\.,[2024](https://arxiv.org/html/2606.02886#bib.bib2)\)requires a pre\-existing ensemble to emulate\. All are tied to specific architectures\. By contrast, the large majority of AI weather checkpoints—including FourCastNetV2, Pangu\-Weather, Aurora, and AIFS—are deterministic and lack native uncertainty estimates\. NTK\-UQ targets this majority: it applies post\-hoc to any pre\-trained deterministic checkpoint without retraining, enabling*checkpoint reusability*across the rapidly growing ecosystem of foundation weather models\.

Existing post\-hoc UQ methods face significant barriers for billion\-parameter models: deep ensembles\(Lakshminarayananet al\.,[2017](https://arxiv.org/html/2606.02886#bib.bib13)\)require training multiple copies \(prohibitively expensive\), Bayesian methods\(Blundellet al\.,[2015](https://arxiv.org/html/2606.02886#bib.bib18); Gal and Ghahramani,[2016](https://arxiv.org/html/2606.02886#bib.bib17)\)need architectural modifications and yield poorly calibrated uncertainties\(Ovadiaet al\.,[2019](https://arxiv.org/html/2606.02886#bib.bib26)\), and conformal prediction\(Angelopoulos and Bates,[2021](https://arxiv.org/html/2606.02886#bib.bib19)\)provides coverage guarantees but lacks per\-sample discrimination\.

The Neural Tangent Kernel\(Jacotet al\.,[2018](https://arxiv.org/html/2606.02886#bib.bib11)\)shows that infinitely wide networks behave as Gaussian Processes, enabling closed\-form uncertainty quantification\. For tractability, the last\-layer empirical NTK uses the feature kernelK​\(x,x′\)=ϕ​\(x\)⊤​ϕ​\(x′\)K\(x,x^\{\\prime\}\)=\\phi\(x\)^\{\\top\}\\phi\(x^\{\\prime\}\)rather than the full gradient\-based NTK\. This coincides with last\-layer Laplace approximation\(MacKay,[1992](https://arxiv.org/html/2606.02886#bib.bib23); Daxbergeret al\.,[2021](https://arxiv.org/html/2606.02886#bib.bib24)\)for linear output heads\. Recent work\(Heet al\.,[2020](https://arxiv.org/html/2606.02886#bib.bib12); Huanget al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib32)\)demonstrates that NTK\-based GP posteriors capture epistemic uncertainty even in finite\-width networks\. UnlikeΔ\\Delta\-UQ\(Thiagarajanet al\.,[2022](https://arxiv.org/html/2606.02886#bib.bib33)\), which requires retraining with anchor perturbation, NTK\-UQ operates entirely post\-hoc on pre\-trained models\. Detailed comparisons are provided in Appendix[A](https://arxiv.org/html/2606.02886#A1)\.

## 3\.Method

NTK\-UQ is a framework for post\-hoc uncertainty quantification in pre\-trained neural weather models\. The method consists of three phases: \(1\) last\-layer feature extraction, \(2\) offline GP posterior construction via kernel decomposition, and \(3\) post\-hoc scaling to achieve target coverage\.

### 3\.1\.Problem Setup

Letfθ:𝒳→𝒴f\_\{\\theta\}:\\mathcal\{X\}\\to\\mathcal\{Y\}be a pre\-trained weather model that maps atmospheric statesx∈𝒳⊂ℝC×H×Wx\\in\\mathcal\{X\}\\subset\\mathbb\{R\}^\{C\\times H\\times W\}to predictionsy∈𝒴⊂ℝC′×H×Wy\\in\\mathcal\{Y\}\\subset\\mathbb\{R\}^\{C^\{\\prime\}\\times H\\times W\}, whereCCandC′C^\{\\prime\}are input and output channels, andH×WH\\times Wis the spatial grid\. Given a calibration dataset𝒟cal=\{\(xi,yi∗\)\}i=1N\\mathcal\{D\}\_\{\\text\{cal\}\}=\\\{\(x\_\{i\},y\_\{i\}^\{\*\}\)\\\}\_\{i=1\}^\{N\}with ground truthyi∗y\_\{i\}^\{\*\}\(used to construct the GP posterior and determine post\-hoc scaling\), the goal is to estimate predictive uncertaintyσ2​\(x\)\\sigma^\{2\}\(x\)such that prediction intervals achieve a target coverage level \(e\.g\., 90% of ground truth values fall within the 90% prediction interval\)\.

### 3\.2\.Gaussian Process Interpretation

Under the last\-layer NTK–GP correspondence, a neural network with last\-layer featuresϕ​\(x\)\\phi\(x\)can be viewed as a Gaussian Process:

\(1\)f​\(x\)∼𝒢​𝒫​\(0,K​\(x,x′\)\),f\(x\)\\sim\\mathcal\{GP\}\(0,K\(x,x^\{\\prime\}\)\),whereK​\(x,x′\)=ϕ​\(x\)⊤​ϕ​\(x′\)K\(x,x^\{\\prime\}\)=\\phi\(x\)^\{\\top\}\\phi\(x^\{\\prime\}\)is the last\-layer empirical NTK \(the feature kernel\)\. Given calibration data, the GP predictive variance at a new pointx∗x\_\{\*\}is:

\(2\)σ2​\(x∗\)=K​\(x∗,x∗\)\+σn2−𝐤∗⊤​\(K\+σn2​I\)−1​𝐤∗,\\sigma^\{2\}\(x\_\{\*\}\)=K\(x\_\{\*\},x\_\{\*\}\)\+\\sigma\_\{n\}^\{2\}\-\\mathbf\{k\}\_\{\*\}^\{\\top\}\(K\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\mathbf\{k\}\_\{\*\},whereK​\(x∗,x∗\)=‖ϕ​\(x∗\)‖2K\(x\_\{\*\},x\_\{\*\}\)=\\\|\\phi\(x\_\{\*\}\)\\\|^\{2\}is the prior variance,σn2\\sigma\_\{n\}^\{2\}is the observation noise variance,𝐤∗=\[K​\(x∗,x1\),…,K​\(x∗,xN\)\]⊤\\mathbf\{k\}\_\{\*\}=\[K\(x\_\{\*\},x\_\{1\}\),\\ldots,K\(x\_\{\*\},x\_\{N\}\)\]^\{\\top\}is the kernel vector to calibration points, andKi​j=K​\(xi,xj\)K\_\{ij\}=K\(x\_\{i\},x\_\{j\}\)is the kernel matrix\. The termσn2\\sigma\_\{n\}^\{2\}in the predictive variance accounts for irreducible noise in the observations and is estimated from the eigenvalue spectrum \(Section[D\.2](https://arxiv.org/html/2606.02886#A4.SS2)\)\.

#### Interpretation\.

The GP posterior variance has a natural interpretation: the kernelKKencodes a prior over the model’s function space shaped by the calibration data geometry\(Leeet al\.,[2018](https://arxiv.org/html/2606.02886#bib.bib36)\)\. For a test inputx∗x\_\{\*\}, the posterior varianceσ2​\(x∗\)\\sigma^\{2\}\(x\_\{\*\}\)quantifies similarity to the calibration distribution in feature space\. Whenx∗x\_\{\*\}is dissimilar to calibration inputs, the correction term𝐤∗⊤​\(K\+σn2​I\)−1​𝐤∗\\mathbf\{k\}\_\{\*\}^\{\\top\}\(K\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\mathbf\{k\}\_\{\*\}is small, and the posterior variance remains close to the prior, yielding high epistemic uncertainty\(Kendall and Gal,[2017](https://arxiv.org/html/2606.02886#bib.bib35)\)\. Conversely, inputs similar to the calibration set receive large corrections, yielding low uncertainty\.

Crucially, the feature mapϕ\\phiis not hand\-crafted but learned from ERA5 reanalysis data spanning decades of global atmospheric observations\. The kernelK​\(x,x′\)=ϕ​\(x\)⊤​ϕ​\(x′\)K\(x,x^\{\\prime\}\)=\\phi\(x\)^\{\\top\}\\phi\(x^\{\\prime\}\)is therefore an*ERA5\-informed similarity measure*: it encodes physically meaningful atmospheric structure—the geometry of realizable weather states as the model learned it from training\. A test input receives high uncertainty when it is unusual relative to both \(1\) the ERA5\-learned feature manifold \(encoded in the frozen weightsϕ\\phi\) and \(2\) the calibration distribution \(encoded in the GP posterior fromnnsamples\)\. This two\-level epistemic signal is inaccessible to purely statistical baselines such as conformal prediction, which operate in prediction\-error space without access to the model’s learned atmospheric representation\.

### 3\.3\.Last\-Layer Feature Extraction

Modern neural weather models decompose asfθ=gψ∘ϕωf\_\{\\theta\}=g\_\{\\psi\}\\circ\\phi\_\{\\omega\}, whereϕω:𝒳→ℝd\\phi\_\{\\omega\}:\\mathcal\{X\}\\to\\mathbb\{R\}^\{d\}extracts features andgψ:ℝd→𝒴g\_\{\\psi\}:\\mathbb\{R\}^\{d\}\\to\\mathcal\{Y\}is the final prediction head\. Last\-layer features are extracted by registering forward hooks during inference\. For spatial feature maps, multi\-statistic aggregation computes six statistics per channel \(mean, standard deviation, minimum, maximum, 25th and 75th percentiles\), yielding a fixed\-dimensional feature vector regardless of spatial resolution\. Architecture\-specific extraction details are provided in Appendix[B\.1](https://arxiv.org/html/2606.02886#A2.SS1)\.

### 3\.4\.Kernel Decomposition

Direct inversion of the kernel matrixKKisO​\(N3\)O\(N^\{3\}\), prohibitive for large calibration sets\. Before decomposition, features are centered by subtracting the calibration mean:ϕ¯=1N​∑iϕ​\(xi\)\\bar\{\\phi\}=\\frac\{1\}\{N\}\\sum\_\{i\}\\phi\(x\_\{i\}\)andϕ~​\(x\)=ϕ​\(x\)−ϕ¯\\tilde\{\\phi\}\(x\)=\\phi\(x\)\-\\bar\{\\phi\}\. This removes the dominant mean direction from the spectrum, ensuring the decomposition captures directions of*variation*rather than the shared mean signal\. The choice of decomposition method significantly affects UQ quality; this work compares Singular Value Decomposition \(SVD\) and Independent Component Analysis \(ICA\) to characterize these effects\.

#### SVD Decomposition\.

The standard approach uses singular value decomposition on the centered feature matrixΦ~∈ℝN×d\\tilde\{\\Phi\}\\in\\mathbb\{R\}^\{N\\times d\}:

\(3\)Φ~=U​S​V⊤\.\\tilde\{\\Phi\}=USV^\{\\top\}\.This decomposition yields the centered kernel eigenstructure directly, sinceK~=Φ~​Φ~⊤=U​S2​U⊤\\tilde\{K\}=\\tilde\{\\Phi\}\\tilde\{\\Phi\}^\{\\top\}=US^\{2\}U^\{\\top\}, meaning the eigenvalues areλj=sj2\\lambda\_\{j\}=s\_\{j\}^\{2\}\(squared singular values\) and the eigenvectors are the columns ofUU\. SVD finds orthogonal directions ofmaximum variancein the feature space\. Retaining only the top\-kkcomponents \(wherek≪Nk\\ll N\) captures the dominant directions of variation\.

#### ICA Decomposition\.

An alternative approach uses Independent Component Analysis \(ICA\)\(Hyvärinen and Oja,[2000](https://arxiv.org/html/2606.02886#bib.bib45)\)to decompose features intostatistically independentcomponents rather than orthogonal directions of maximum variance\. ICA assumes that the observed featuresϕ~​\(x\)\\tilde\{\\phi\}\(x\)are linear mixtures of independent source signals:ϕ~​\(x\)=A​s​\(x\)\\tilde\{\\phi\}\(x\)=As\(x\)wheres​\(x\)s\(x\)are the independent components andAAis the mixing matrix\. The FastICA algorithm\(Hyvärinen,[1999](https://arxiv.org/html/2606.02886#bib.bib46)\)recovers the unmixing matrixW=A−1W=A^\{\-1\}by maximizing non\-Gaussianity of the sources, yielding componentss​\(x\)=W​ϕ~​\(x\)s\(x\)=W\\tilde\{\\phi\}\(x\)\. Unlike SVD, which prioritizes variance, ICA exploits higher\-order statistics \(kurtosis, skewness\) to separate sources\.

For extreme weather events, ICA offers a critical advantage: while SVD’s maximum\-variance criterion biases the decomposition toward typical weather patterns \(high\-frequency, high\-variance modes\), ICA’s independence criterion can isolate rare extreme event signatures that occur as independent factors in the joint distribution, even when they contribute low marginal variance\. Empirical results show that ICA outperforms SVD for uncertainty quantification in extreme events for three of four architectures \(AIFS, Aurora, FourCastNetV2\); SVD achieves coverage for Pangu\-Weather while ICA fails \(Section[5](https://arxiv.org/html/2606.02886#S5)\)\.

The predictive variance formula becomes:

\(4\)σraw2​\(x∗\)=‖ϕ~​\(x∗\)‖2\+σn2⏟prior \+ noise−∑j=1kλj⋅\(ϕ~​\(x∗\)⊤​vj\)2λj\+σn2⏟GP correction,\\sigma^\{2\}\_\{\\text\{raw\}\}\(x\_\{\*\}\)=\\underbrace\{\\\|\\tilde\{\\phi\}\(x\_\{\*\}\)\\\|^\{2\}\+\\sigma\_\{n\}^\{2\}\}\_\{\\text\{prior \+ noise\}\}\-\\underbrace\{\\sum\_\{j=1\}^\{k\}\\frac\{\\lambda\_\{j\}\\cdot\(\\tilde\{\\phi\}\(x\_\{\*\}\)^\{\\top\}v\_\{j\}\)^\{2\}\}\{\\lambda\_\{j\}\+\\sigma\_\{n\}^\{2\}\}\}\_\{\\text\{GP correction\}\},whereϕ~​\(x∗\)=ϕ​\(x∗\)−ϕ¯\\tilde\{\\phi\}\(x\_\{\*\}\)=\\phi\(x\_\{\*\}\)\-\\bar\{\\phi\}is the centered test feature,vjv\_\{j\}are the right singular vectors ofΦ~\\tilde\{\\Phi\},λj=sj2\\lambda\_\{j\}=s\_\{j\}^\{2\}are the eigenvalues, andσn2\\sigma\_\{n\}^\{2\}is the noise variance\. Projections onto high\-variance directions receive large corrections \(low uncertainty\); dissimilar inputs receive small corrections \(high uncertainty\)\.

#### Noise Variance Estimation\.​​

The noise parameterσn2\\sigma\_\{n\}^\{2\}is estimated as the mean of the residual eigenvalues\{λk\+1,…,λd\}\\\{\\lambda\_\{k\+1\},\\ldots,\\lambda\_\{d\}\\\}\. When the top\-kkcomponents exhaust all variance this estimate approaches zero, causing the correction to fully cancel the prior and destroying discrimination\. The method falls back to the mean of the retained eigenvalues as a regularization nugget\(Cressie,[1993](https://arxiv.org/html/2606.02886#bib.bib15)\), preserving posterior variation even when the feature space is low\-rank\.

###### Proposition 1 \(Variance Collapse\)\.

LetΦ~=U​S​V⊤\\tilde\{\\Phi\}=USV^\{\\top\}be the SVD of the centered calibration features and define the correction\-to\-prior ratioRk=Ck/PR\_\{k\}=C\_\{k\}/PwhereCk=∑j=1kλj​cj2/\(λj\+σn2\)C\_\{k\}=\\sum\_\{j=1\}^\{k\}\\lambda\_\{j\}c\_\{j\}^\{2\}/\(\\lambda\_\{j\}\+\\sigma\_\{n\}^\{2\}\)andP=‖ϕ~​\(x∗\)‖2P=\\\|\\tilde\{\\phi\}\(x\_\{\*\}\)\\\|^\{2\}\. When the noise regularizerσn2\>0\\sigma\_\{n\}^\{2\}\>0, each shrinkage weightwj=λj/\(λj\+σn2\)<1w\_\{j\}=\\lambda\_\{j\}/\(\\lambda\_\{j\}\+\\sigma\_\{n\}^\{2\}\)<1, soRk<1R\_\{k\}<1andσ2​\(x∗\)\>0\\sigma^\{2\}\(x\_\{\*\}\)\>0for all rankskk– no collapse occurs\. Whenσn2=0\\sigma\_\{n\}^\{2\}=0, askkapproaches the true rankrr,Rk→1R\_\{k\}\\to 1andσ2​\(x∗\)→0\\sigma^\{2\}\(x\_\{\*\}\)\\to 0, destroying uncertainty discrimination\. The actionable diagnostic is to maintainRk<0\.9R\_\{k\}<0\.9before deployment; this threshold is empirically validated in Table[7](https://arxiv.org/html/2606.02886#A5.T7)\(collapse atk=100k=100whereRk=0\.92R\_\{k\}=0\.92\) rather than derived from the proof\. Proof in Appendix[C](https://arxiv.org/html/2606.02886#A3)\.

###### Proposition 2 \(Non\-Gaussian Discrimination\)\.

When feature distributions exhibit joint non\-Gaussianity \(higher\-order cumulantsκi1,…,im≠0\\kappa\_\{i\_\{1\},\\ldots,i\_\{m\}\}\\neq 0form≥3m\\geq 3\), SVD captures only second\-order structure \(the covariance matrix\), discarding tail behavior and higher\-order dependencies, whereas ICA exploits kurtosis and negentropy to isolate statistically independent sources\. For extreme weather events with heavy\-tailed marginals, ICA components aligned with extreme directions achieve higher kurtosis, producing adaptive uncertainty estimates that SVD cannot recover\. Full formalization and proof in Appendix[B\.6](https://arxiv.org/html/2606.02886#A2.SS6); main\-text propositions are accessible summaries with complete proofs in the appendix\.

#### Decomposition Selection\.

Which method to use depends on eigenspectrum concentration\. Letλ1≥⋯≥λd\\lambda\_\{1\}\\geq\\cdots\\geq\\lambda\_\{d\}be the eigenvalues of the centered feature covariance\. Algorithm[1](https://arxiv.org/html/2606.02886#alg1)gives a data\-driven selection rule validated empirically in Section[5\.2](https://arxiv.org/html/2606.02886#S5.SS2)\.

Algorithm 1ICA/SVD Decomposition Selection1:Input:Centered feature matrix

Φ~∈ℝn×d\\tilde\{\\Phi\}\\in\\mathbb\{R\}^\{n\\times d\}, calibration set

2:Compute eigenvalues

λ1≥⋯≥λd\\lambda\_\{1\}\\geq\\cdots\\geq\\lambda\_\{d\}of

Φ~⊤​Φ~\\tilde\{\\Phi\}^\{\\top\}\\tilde\{\\Phi\}
3:Compute concentration ratio

ρ=λ1/∑j=1dλj\\rho=\\lambda\_\{1\}/\\sum\_\{j=1\}^\{d\}\\lambda\_\{j\}
4:if

ρ\>0\.8\\rho\>0\.8then⊳\\trianglerightConcentrated spectrum \(e\.g\., SFNO\)

5:UseSVDwith

k≤10k\\leq 10; select

kkby

Rk=Ck/P<0\.9R\_\{k\}=C\_\{k\}/P<0\.9\(Prop\. 1\)

6:elseif

ρ<0\.5\\rho<0\.5then⊳\\trianglerightDistributed spectrum \(e\.g\., GNN\-Transformer, Perceiver\)

7:UseICA; select

kkby CRPS on a held\-out validation split, subject to coverage

≥85%\\geq 85\\%
8:else⊳\\trianglerightIntermediate: compare both by CRPS

9:Run both on held\-out validation split; use method with lower CRPS and coverage

≥85%\\geq 85\\%
10:endif

11:Output:Decomposition method and rank

kk

### 3\.5\.Post\-hoc Calibration Scaling

Raw NTK uncertainties capturerelativeuncertainty ordering across samples but not theabsolutemagnitude: empirical coverage is typically well below the target level \(e\.g\., 50% instead of 90%\)\. A scaling factorα\\alphais learned per lead time via binary search such thatσcal=α⋅σraw\\sigma\_\{\\text\{cal\}\}=\\alpha\\cdot\\sigma\_\{\\text\{raw\}\}achieves target 90% coverage\. This is equivalent to temperature scaling\(Guoet al\.,[2017](https://arxiv.org/html/2606.02886#bib.bib21)\)applied to the GP variance, withα\\alphaplaying the role of the temperature parameter\. Per\-variable calibration learns separate scalesαv\\alpha\_\{v\}for each meteorological variable, accommodating their different error characteristics\. The binary search algorithm is detailed in Appendix[B\.3](https://arxiv.org/html/2606.02886#A2.SS3)\.

### 3\.6\.Autoregressive Feature Extraction

Rather than running separate forward passes for each forecast horizon, features are extracted at multiple checkpoints during a single autoregressive rollout, reducing computational cost by a factor of\|𝒯\|\|\\mathcal\{T\}\|\(the number of target horizons\)\. Implementation details are provided in Appendix[B\.4](https://arxiv.org/html/2606.02886#A2.SS4)\.

## 4\.Experimental Setup

### 4\.1\.Models

Experiments evaluate NTK\-UQ on four production AI weather models representing diverse architectural approaches\. FourCastNetV2\(Pathaket al\.,[2022](https://arxiv.org/html/2606.02886#bib.bib4)\)uses NVIDIA’s Spherical Fourier Neural Operator \(SFNO\)\(Bonevet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib5)\)architecture with 73 input channels at 0\.25° resolution\. Pangu\-Weather\(Biet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib6)\)employs Huawei’s 3D Swin Transformer with separate 6\-hour and 24\-hour prediction models in ONNX format, using 69 input channels\. Aurora\(Bodnaret al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib8)\)is Microsoft’s foundation model combining a 3D Swin Transformer backbone with Perceiver\-based encoders and decoders, fine\-tuned for 0\.25° ERA5 data with 69 input channels\. AIFS\(Langet al\.,[2024](https://arxiv.org/html/2606.02886#bib.bib3)\)is ECMWF’s operational model combining graph neural network encoding on an icosahedral mesh with transformer\-based processing, using 69 input channels\. These models were selected to demonstrate that NTK\-UQ generalizes across fundamentally different neural architectures \(Fourier operators, vision transformers, perceiver networks, and GNN\-transformer hybrids\)\.

### 4\.2\.Data

Experiments use ERA5 reanalysis\(Raspet al\.,[2024](https://arxiv.org/html/2606.02886#bib.bib9)\)at 0\.25° resolution, following standard practice in AI weather model evaluation\. Evaluation focuses on extreme weather events from 2021, ensuring out\-of\-distribution temporal evaluation \(all four models were trained on data ending before 2021\)\. The dataset comprises initialization dates from the EM\-DAT International Disaster Database\(Delforgeet al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib44)\), constituting a near\-complete census of verified high\-impact events in 2021 \(not a random sample\): 136 flood events, 63 storms \(tropical cyclones Tauktae, Ida, Rai, Elsa\), 5 droughts, and 2 extreme temperature events \(June 2021 Pacific Northwest heat wave\) across 82 countries, yieldingn=100n=100distinct initialization dates after deduplication\. Initialization dates are selected 3 days before event onset to capture the development phase where forecast uncertainty is most critical\. Features are extracted at lead timesτ∈\{6,12,24,48,72,120\}\\tau\\in\\\{6,12,24,48,72,120\\\}hours during autoregressive rollouts\. Detailed dataset construction and training data overlap analysis are provided in Appendix[B\.2](https://arxiv.org/html/2606.02886#A2.SS2)\.

### 4\.3\.Evaluation Metrics

Uncertainty quantification quality is evaluated using a principled three\-tier framework with Sharpness as the primary optimization target, Coverage as a hard constraint, and CRPS as the overall score\.

#### Sharpness \(Primary Metric\)\.

Sharpness measures the tightness of prediction intervals, computed as the mean uncertainty width:

\(5\)Sharpness=1\|𝒱\|​∑i∈𝒱σi\\text\{Sharpness\}=\\frac\{1\}\{\|\\mathcal\{V\}\|\}\\sum\_\{i\\in\\mathcal\{V\}\}\\sigma\_\{i\}Lower sharpness is better—narrower intervals provide more informative forecasts\. Sharpness directly quantifies the primary goal of UQ: to minimize uncertainty while maintaining reliability\. However, sharpness alone is insufficient; intervals can be arbitrarily narrow \(sharp\) but miscalibrated\. This motivates the coverage constraint\.

#### Coverage \(Constraint\)\.

Coverage measures the fraction of ground truth values falling within thepp% prediction interval:

\(6\)Coverage​\(p\)=1\|𝒱\|​∑i∈𝒱𝟏​\[\|yi∗−y^i\|≤zp⋅σi\]\\text\{Coverage\}\(p\)=\\frac\{1\}\{\|\\mathcal\{V\}\|\}\\sum\_\{i\\in\\mathcal\{V\}\}\\mathbf\{1\}\\left\[\|y\_\{i\}^\{\*\}\-\\hat\{y\}\_\{i\}\|\\leq z\_\{p\}\\cdot\\sigma\_\{i\}\\right\]wherezpz\_\{p\}is the corresponding normal quantile \(e\.g\.,z0\.95≈1\.645z\_\{0\.95\}\\approx 1\.645for two\-sided 90% intervals\)\. Well\-calibrated UQ satisfies Coverage\(90%\)∈\[0\.85,0\.95\]\\in\[0\.85,0\.95\]; values below 0\.85 indicate overconfidence \(intervals too narrow\), while values above 0\.95 indicate underconfidence \(intervals too wide\)\. Coverage is treated as a hard constraint rather than an optimization target: methods must achieve the target coverage to be considered valid, but among valid methods, the sharpest \(tightest\) intervals are preferred\.

#### CRPS \(Overall Score\)\.

The CRPS\(Gneiting and Raftery,[2007](https://arxiv.org/html/2606.02886#bib.bib27)\)is a proper scoring rule that jointly evaluates sharpness and calibration:

\(7\)CRPS=𝔼​\[∫−∞\+∞\(F​\(y\)−𝟏​\[y≥y∗\]\)2​𝑑y\]\.\\text\{CRPS\}=\\mathbb\{E\}\\\!\\left\[\\textstyle\\int\_\{\-\\infty\}^\{\+\\infty\}\\bigl\(F\(y\)\-\\mathbf\{1\}\[y\\geq y^\{\*\}\]\\bigr\)^\{2\}dy\\right\]\.For Gaussian predictive distributions𝒩​\(y^,σ2\)\\mathcal\{N\}\(\\hat\{y\},\\sigma^\{2\}\), CRPS has a closed form\. Lower CRPS indicates better overall probabilistic forecast quality\. CRPS rewards both accuracy \(low bias\) and sharpness \(low variance\) while penalizing miscalibration\.

#### Error\-Uncertainty Correlation \(Diagnostic\)\.

Spearman rank correlation between absolute errors and predicted uncertainties\(Tranet al\.,[2020](https://arxiv.org/html/2606.02886#bib.bib16)\)provides a diagnostic measure of discrimination:

\(8\)ρs=1−6​∑i=1Ndi2N​\(N2−1\),\\rho\_\{s\}=1\-\\frac\{6\\sum\_\{i=1\}^\{N\}d\_\{i\}^\{2\}\}\{N\(N^\{2\}\-1\)\},wheredi=rank​\(\|ei\|\)−rank​\(σi\)d\_\{i\}=\\text\{rank\}\(\|e\_\{i\}\|\)\-\\text\{rank\}\(\\sigma\_\{i\}\)is the difference between the rank of the absolute error\|ei\|=\|yi∗−y^i\|\|e\_\{i\}\|=\|y\_\{i\}^\{\*\}\-\\hat\{y\}\_\{i\}\|and the rank of the predicted uncertaintyσi\\sigma\_\{i\}\. Higherρs\\rho\_\{s\}indicates that uncertainty estimates meaningfully rank extreme event difficulty: intense tropical cyclones and atmospheric rivers should receive higher uncertainty than typical synoptic conditions\. Values above 0\.3 are generally considered adequate\. This work reportsρs\\rho\_\{s\}as supplementary evidence of UQ quality but does not use it as a primary evaluation criterion, as it can be high even for poorly calibrated intervals\.

#### Discrimination via Uncertainty Variation\.

The coefficient of variation \(CV\) of predicted uncertainties measures the method’s capacity to distinguish forecast difficulty:

\(9\)CV=1N​∑i=1N\(σi−σ¯\)2σ¯,\\text\{CV\}=\\frac\{\\sqrt\{\\tfrac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\(\\sigma\_\{i\}\-\\bar\{\\sigma\}\)^\{2\}\}\}\{\\bar\{\\sigma\}\},whereσ¯=1N​∑i=1Nσi\\bar\{\\sigma\}=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\sigma\_\{i\}and\{σi\}i=1N\\\{\\sigma\_\{i\}\\\}\_\{i=1\}^\{N\}are the GP posterior standard deviations\. Higher CV indicates that the method produces heterogeneous rather than uniform intervals\. CV\>0\.3\>0\.3indicates substantial per\-sample variation, while CV<0\.1<0\.1indicates nearly uniform intervals\. Note that CV measures variation but not directionality: whether high\-uncertainty samples correspond to genuinely difficult forecasts is verified separately by Spearmanρs\\rho\_\{s\}\(Table[6](https://arxiv.org/html/2606.02886#A5.T6)\)\.

## 5\.Results

Experiments validate the theoretical predictions \(Propositions[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)and[2](https://arxiv.org/html/2606.02886#S3.Thmtheorem2)\) on four AI weather models using disaster\-precursor dates from 2021\. For each model and lead time, the GP posterior is constructed from extracted features, and post\-hoc calibration scales uncertainties to achieve target coverage\. Results are reported for six lead times: 6, 12, 24, 48, 72, and 120 hours, across 17 meteorological variables \(6 surface \+ 11 pressure\-level\)\.

### 5\.1\.Calibration Quality

Table[1](https://arxiv.org/html/2606.02886#S5.T1)presents 90% prediction interval coverage for 2 m temperature across models and lead times\. Post\-hoc calibration achieves the target 89–91% coverage for all four models across all forecast horizons\. However, achieving target coverage is necessary but not sufficient: Table[5](https://arxiv.org/html/2606.02886#A5.T5)shows that discrimination quality depends critically on the decomposition method\. ICA produces adaptive intervals \(higher coefficient of variation\) that scale with extreme event severity, while SVD produces more uniform intervals that fail to distinguish tropical cyclones from routine weather\. The discrimination quality is captured by the Spearman correlation \(Table[6](https://arxiv.org/html/2606.02886#A5.T6)\)\.

Table 1\.Coverage at 90% prediction interval for 2 m temperature \(t2m\) by model and lead time\. Post\-hoc scaling achieves near\-target coverage for all models\.All models use per\-variable post\-hoc scaling \(Section[3\.5](https://arxiv.org/html/2606.02886#S3.SS5)\)\. Coverage is achieved with method\-dependent discrimination quality: ICA produces adaptive intervals \(CV=0\.07=0\.07–1\.81\), while SVD produces more uniform intervals \(CV=0\.01=0\.01–0\.49\)\. See Table[5](https://arxiv.org/html/2606.02886#A5.T5)for details\.

###### Theorem 1 \(Post\-Hoc Coverage Bound\)\.

Letc^n\\hat\{c\}\_\{n\}be the empirical coverage onnni\.i\.d\. calibration samples\. For anyδ∈\(0,1\)\\delta\\in\(0,1\), with probability at least1−δ1\-\\delta:

\(10\)ctrue≥c^n−ln⁡\(1/δ\)2​n\.c\_\{\\mathrm\{true\}\}\\geq\\hat\{c\}\_\{n\}\-\\sqrt\{\\frac\{\\ln\(1/\\delta\)\}\{2n\}\}\.

\(Proof via one\-sided Hoeffding inequality applied to Bernoulli coverage indicators; see Appendix[C](https://arxiv.org/html/2606.02886#A3)\.\)

Withn=100n=100held\-out evaluation samples achievingc^n=0\.90\\hat\{c\}\_\{n\}=0\.90empirical coverage \(post\-hoc scaleα\\alphais fixed from the calibration set; coverage is evaluated on independent data\), true coverage exceeds0\.7780\.778with 95% confidence \(δ=0\.05\\delta=0\.05\)\. The 85% floor used to filter valid comparisons is a practical threshold: it excludes clearly miscalibrated configurations \(e\.g\., Aurora under SVD at 58\.1% coverage\) while providing a margin above the 77\.8% Hoeffding worst\-case floor\.

### 5\.2\.Decomposition Method Comparison

Proposition[2](https://arxiv.org/html/2606.02886#S3.Thmtheorem2)predicts that ICA exploits non\-Gaussian, heavy\-tailed structure in extreme weather events to achieve higher discrimination than SVD\. Table[2](https://arxiv.org/html/2606.02886#S5.T2)validates this prediction empirically by comparing the two kernel decomposition methods across all four models, with each method evaluated at its optimal rankk∗k^\{\*\}selected via coverage\-constrained CRPS minimization on held\-out data\. Results confirm that optimal method selection depends on model architecture and feature distribution properties\.

Table 2\.ICA vs SVD decomposition comparison at optimal rankk∗k^\{\*\}per method\. Coverage must satisfy 85–95% constraint; sharpness \(meanσ\\sigma\) should be minimized subject to coverage; CRPS provides overall score\. Bold indicates method satisfying coverage constraint\.Optimal rankk∗k^\{\*\}selected per method via coverage\-constrained CRPS minimization\. Coverage is the primary constraint \(must be 85–95%\)\. ICA achieves proper coverage for 3/4 models \(AIFS, Aurora, FCNv2\) with lower CRPS\. SVD fails coverage for Aurora \(58\.1%\)\. For FCNv2, both methods satisfy coverage, but ICA achieves lower CRPS \(61\.5 vs 64\.5\) at higher sharpness\. Pangu exhibits numerical instabilities \(SVDσ\>35,000\\sigma\>35\{,\}000\) but SVD satisfies coverage while ICA fails\.†Pangu\-ICA \(68\.4% coverage\) does not satisfy the 85% constraint and is excluded from valid comparisons\.

The coverage constraint \(85–95%\) serves as the primary filter: methods failing this constraint produce unreliable prediction intervals regardless of sharpness or CRPS\. Among methods satisfying coverage, sharpness quantifies interval tightness \(lower is better\), while CRPS provides an aggregate score combining calibration and sharpness\.

At optimal ranks,ICA satisfies coverage for three models\(AIFS, Aurora, FourCastNetV2\) whileSVD satisfies coverage for only two\(AIFS, Pangu\)\. ForAurora\(PerceiverIO\), ICA achieves target coverage \(90\.1%,k∗=50k^\{\*\}=50\) and lower CRPS, while SVD severely underfits \(58\.1% coverage atk∗=2k^\{\*\}=2\), indicating intervals too narrow to capture forecast errors\. ForFourCastNetV2\(SFNO\), both methods satisfy coverage \(89\.5%\), but ICA achieves lower CRPS \(61\.5 vs 64\.5 atk∗=3k^\{\*\}=3vsk∗=1k^\{\*\}=1\)\. ForAIFS\(GNN\-Transformer\), both methods satisfy coverage, but ICA achieves lower CRPS \(129\.8 vs 133\.5 atk∗=7k^\{\*\}=7vs 50\)\. Only forPangu\-Weather\(Swin Transformer\) does SVD outperform, achieving 91\.1% coverage while ICA fails \(68\.4% atk∗=40k^\{\*\}=40\), though both exhibit numerical instabilities \(extremely largeσ\\sigmavalues\)\.

Figure[2](https://arxiv.org/html/2606.02886#S5.F2)illustrates CRPS evolution across forecast horizons \(6h to 120h\) for all four models using both decomposition methods, evaluated on five major extreme weather events from the EM\-DAT international disaster database representing operational scenarios where accurate probabilistic forecasts are most critical\.

FourCastNetV2 with ICA achieves the lowest CRPS \(20–150\) across all horizons, indicating both sharp and well\-calibrated intervals for extreme weather events\. Aurora with ICA shows moderate CRPS \(300–1,400\), while AIFS with ICA achieves CRPS in the 100–250 range\. Pangu\-Weather exhibits elevated CRPS values \(SVD: 10,000–40,000; ICA: 15,000–80,000\) in its 69\-dimensional feature space, though SVD maintains target coverage \(Table[2](https://arxiv.org/html/2606.02886#S5.T2)\)\. The consistent separation between ICA and SVD curves demonstrates that decomposition method selection impacts not only coverage calibration but also the overall probabilistic forecast quality measured by CRPS\.

![Refer to caption](https://arxiv.org/html/2606.02886v1/x1.png)Figure 2\.CRPS vs lead time for all four models using ICA \(left\) and SVD \(right\) decomposition\. Lower CRPS indicates better probabilistic forecast quality\. Evaluated on five EM\-DAT extreme weather events \(Tropical Cyclone Tauktae, Tropical Cyclone Ida, Pacific Northwest heat wave, Central European floods, Typhoon Rai\)\. FourCastNetV2 with ICA achieves the lowest CRPS \(20–150\) across all horizons\. Note: Pangu\-Weather plotted on log scale due to numerical instabilities\.Figures[3](https://arxiv.org/html/2606.02886#S5.F3)and[4](https://arxiv.org/html/2606.02886#S5.F4)show sharpness evolution \(meanσ\\sigmawith median and IQR bands\) for 2\-meter temperature and mean sea level pressure\. ICA achieves lower CRPS than SVD for AIFS, FourCastNetV2, and Aurora, indicating better overall probabilistic quality despite similar or wider meanσ\\sigmafor AIFS and FourCastNetV2 \(where SVD achieves lower meanσ\\sigmabut higher CRPS\)\. The wider IQR for ICA indicates adaptive intervals that scale with event severity \(tropical cyclones like Typhoon Rai receiveσ\>500\\sigma\>500, routine conditions receiveσ<100\\sigma<100\), while SVD’s narrower IQR indicates more uniform intervals\. This adaptive behavior validates Proposition[2](https://arxiv.org/html/2606.02886#S3.Thmtheorem2): ICA exploits higher\-order statistics in non\-Gaussian extreme weather features to discriminate event difficulty, while SVD captures only second\-order variance structure\.

![Refer to caption](https://arxiv.org/html/2606.02886v1/x2.png)Figure 3\.Sharpness \(mean uncertaintyσ\\sigmawith median and IQR bands\) vs lead time for 2\-meter temperature\. Each row shows one model; columns compare SVD \(left\) vs ICA \(right\) decomposition\. ICA achieves lower CRPS than SVD for most models \(Table[2](https://arxiv.org/html/2606.02886#S5.T2)\); SVD achieves lower meanσ\\sigmafor AIFS and FourCastNetV2 but higher CRPS\. Wider IQR for ICA indicates adaptive intervals that scale with extreme event severity\.![Refer to caption](https://arxiv.org/html/2606.02886v1/x3.png)Figure 4\.Sharpness \(mean uncertaintyσ\\sigmawith median and IQR bands\) vs lead time for mean sea level pressure\. Layout identical to Figure[3](https://arxiv.org/html/2606.02886#S5.F3)\. ICA’s wider IQR reflects event\-specific adaptation: tropical cyclones like Typhoon Rai receive wider intervals \(σ\>1000\\sigma\>1000Pa\) while routine mid\-latitude conditions receive narrower intervals \(σ<500\\sigma<500Pa\)\.ICA produces substantially higher uncertainty variation than SVD \(CV=0\.27=0\.27–1\.811\.81vs\.0\.010\.01–0\.490\.49\), with AIFS maintaining the strongest directional discrimination \(ρs=0\.25\\rho\_\{s\}=0\.25–0\.330\.33across lead times\)\. Full CV and Spearmanρs\\rho\_\{s\}results are in Appendix Tables[5](https://arxiv.org/html/2606.02886#A5.T5)and[6](https://arxiv.org/html/2606.02886#A5.T6)\.

### 5\.3\.Comparison with Conformal Prediction

To situate NTK\-UQ relative to an established post\-hoc baseline, we compare against split conformal prediction \(80/20 split of the samen=100n=100calibration samples,q0\.90q\_\{0\.90\}nonconformity score equal to the per\-variable empirical RMSE quantile\)\. Conformal prediction provides distribution\-free coverage guarantees but produces*uniform*prediction intervals – a singleq^\\hat\{q\}per variable applied to all inputs regardless of event severity\. Table[3](https://arxiv.org/html/2606.02886#S5.T3)reports mean uncertainty widthσ\\sigmaand observed coverage for key variables and lead times, using per\-variable post\-hoc calibration for all three methods\.

Table 3\.Sharpness comparison: NTK\-UQ \(ICA and SVD\) vs\. conformal prediction\. Each cell shows meanσ\\sigma\(coverage%\)\.Bold= sharpest method with coverage≥85%\\geq 85\\%\. ‘–’ = unavailable\. Pangu excluded \(numerical instabilities\); Aurora in model\-normalized units\.All three methods achieve≈\\approx90% empirical coverage\. Across the full evaluation spanning 17 meteorological variables and six lead times, NTK\-UQ achieves lowerσ\\sigmathan conformal prediction in81% of valid comparisons\(230/284, coverage≥85%\\geq 85\\%\), using the better\-performing NTK\-UQ variant \(ICA or SVD\) per comparison as determined by Algorithm[1](https://arxiv.org/html/2606.02886#alg1)\. Table[3](https://arxiv.org/html/2606.02886#S5.T3)shows representative cases: SVD is 31–37% sharper than conformal for AIFS and FourCastNetV2 on t2m and msl\. For AIFS and Aurora, ICA achieves better CRPS than SVD despite similar meanσ\\sigma, consistent with Proposition[2](https://arxiv.org/html/2606.02886#S3.Thmtheorem2): ICA’s exploitation of higher\-order statistics produces better\-calibrated intervals under non\-Gaussian heavy\-tailed features, even when raw interval widths are comparable\.

The key distinction is*adaptive sharpness*: while conformal prediction assigns a single interval width per variable \(CV=0=0by construction\), NTK\-UQ produces heterogeneous intervals \(CV=0\.07=0\.07–1\.81, Appendix Table[5](https://arxiv.org/html/2606.02886#A5.T5)\)\. This variation is a necessary condition for distinguishing a tropical cyclone from routine conditions; it is not sufficient alone – whether the variation is correctly directed \(highσ\\sigmafor difficult forecasts, lowσ\\sigmafor easy ones\) is measured by Spearmanρs\\rho\_\{s\}\. AIFS achieves meaningful directional discrimination \(ρs=0\.25\\rho\_\{s\}=0\.25–0\.330\.33across lead times\); other models show weaker but positive correlation\. Conformal, with CV=0=0by construction, cannot achieve positiveρs\\rho\_\{s\}regardless of sample size\.

Variance collapse \(Proposition[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)\) is empirically validated for FourCastNetV2 in Appendix Table[7](https://arxiv.org/html/2606.02886#A5.T7): atk=100k=100,Rk=Ck/P=0\.92R\_\{k\}=C\_\{k\}/P=0\.92, leaving only 8% residual variance;k≤10k\\leq 10maintainsRk<15%R\_\{k\}<15\\%\.

## 6\.Discussion

![Refer to caption](https://arxiv.org/html/2606.02886v1/figures/fig_aifs_typhoon_odette.png)Figure 5\.AIFS spatial uncertainty for Typhoon Odette at t\+12h \(2021\-12\-16, Philippines region\)\. \(a\) ERA5 ground truth shows mean sea level pressure\. \(b\) AIFS forecast\. \(c\) Absolute forecast error concentrates near the cyclone track, with maximum error exceeding 5200 Pa\. \(d\) NTK\-UQ uncertainty map \(σ\\sigma\) showing spatial variation in epistemic uncertainty\. AIFS with ICA atk=20k=20exhibits spatially\-varying uncertainty patterns that correlate with forecast error magnitude, demonstrating that ICA’s exploitation of non\-Gaussian structure enables finer\-grained discrimination between high\-error \(typhoon core\) and low\-error \(surrounding regions\) areas\.Four\-panel map of Philippines showing AIFS forecast, ground truth, error, and spatial NTK uncertainty for Typhoon Odette at t\+12h\. Uncertainty pattern shows spatial variation correlated with forecast errors\.Empirical validation across four foundation weather models confirms the theoretical predictions: NTK\-UQ achieves calibrated coverage \(89\-91%\), with discrimination quality following the architecture\-dependent and decomposition\-dependent patterns predicted by Propositions[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)and[2](https://arxiv.org/html/2606.02886#S3.Thmtheorem2)\. The following subsections analyze these findings and their implications\.

### 6\.1\.Architecture\-Dependent Behavior

A key finding is the strong dependence of NTK\-UQ behavior on both neural network architecture and decomposition method\. All four models achieve target 90% coverage after post\-hoc scaling, but the ability to discriminate extreme event difficulty varies substantially\. This behavior is explained by the spectral characterization in Section[C\.1](https://arxiv.org/html/2606.02886#A3.SS1): architectures with concentrated eigenvalue spectra \(e\.g\., SFNO’s global Fourier basis\) yield low effective rank, and when the truncation rankkkapproaches the effective rank, the correction term consumes nearly all prior variance \(Proposition[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)\), collapsing extreme event discrimination\.

Rank selection must match the architecture: spectral models needk≤10k\\leq 10; attention\-based models tolerate full rank\. MaintainRk<0\.9R\_\{k\}<0\.9as a pre\-deployment check for uncertainty discrimination\.

#### Implications for Model Design\.

These findings have significant implications for designing the next generation of weather AI architectures\. If post\-hoc uncertainty quantification is a deployment requirement, architectural choices should favor inductive biases that produce desirable eigenspectrum properties\. The variance collapse analysis \(Proposition[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)\) provides a predictive diagnostic: architectures with concentrated eigenspectra \(e\.g\., SFNO’s global Fourier basis\) will require aggressive rank truncation for NTK\-UQ, while architectures with distributed spectra \(attention\-based models, graph networks\) tolerate full\-rank computation and exhibit more robust UQ behavior\. This suggests that uncertainty\-aware architectural design should consider not only forecast accuracy but also the spectral properties of learned feature representations, selecting inductive biases that enable efficient post\-hoc UQ without retraining\.

#### Architecture\-Dependent Spatial Uncertainty Patterns\.

Figure[5](https://arxiv.org/html/2606.02886#S6.F5)demonstrates architecture\-dependent spatial uncertainty structure\. AIFS \(GNN\-Transformer\) with ICA decomposition exhibits spatially\-varying uncertainty that aligns with forecast error concentrations, while models with global pooling \(e\.g\., Pangu\-Weather’sd=69d=69features\) produce scalar uncertainty per variable\. This difference stems from feature representation: AIFS’s graph\-based architecture preserves local spatial structure in itsd=1024d=1024dimensional feature space, enabling ICA to isolate spatially\-coherent independent components that correlate with regional forecast difficulty\. For operational typhoon forecasting, spatially\-varying uncertainty enables targeted warnings for high\-risk regions \(landfall zones, population centers\) rather than uniform domain\-wide alerts\.

### 6\.2\.Limitations and Practical Considerations

The NTK\-UQ framework relies on three approximations: the last\-layer NTK–GP correspondence holds rigorously only at infinite width, the last\-layer restriction ignores early\-layer uncertainty contributions, and unlike conformal prediction it lacks distribution\-free coverage guarantees\. However, empirical coverage consistently achieves 90% on held\-out data, and Theorem[1](https://arxiv.org/html/2606.02886#S5.Thmtheorem1)provides a worst\-case floor of 77\.8% with 95% confidence \(n=100n=100\)\. The evaluation is also bounded in scope: a single out\-of\-distribution year \(2021,n=100n=100deduplicated events\) with temporally autocorrelated precursors\. Cross\-year and cross\-resolution validation remains future work; the post\-hoc design makes such recalibration straightforward, requiring only a re\-run of the offline stage on expanded data\.

The framework’s post\-hoc nature and negligible inference overhead \(a single matrix\-vector product per sample\) make it applicable in resource\-constrained settings where ensemble methods are infeasible\. For operational extreme weather warning systems, ICA is preferred over SVD despite slightly higher computational cost, as it produces adaptive intervals that distinguish tropical cyclone forecasts from routine conditions—critical for disaster preparedness\.

#### Broader Impact\.

Calibrated, spatially\-adaptive uncertainty is directly actionable for extreme\-weather early warning: per\-grid\-point intervals support targeted alerts for high\-risk regions \(landfall zones, population centers\) rather than uniform domain\-wide warnings, and the post\-hoc, model\-agnostic design lets any deployed deterministic checkpoint gain uncertainty estimates without retraining, lowering the barrier to trustworthy forecasting for under\-resourced meteorological agencies\. Because the intervals are calibrated to historical ERA5 reanalysis rather than direct observations, they should be validated against local station data before operational deployment, particularly in regions with sparse observational coverage\.

## 7\.Conclusion

This paper presents a systematic study of last\-layer NTK\-based uncertainty quantification across four foundation weather models, comparing SVD and ICA decomposition methods\. The framework requires no retraining, adds minimal inference overhead, and achieves calibrated prediction intervals when properly matched to model architecture\.

Two key findings emerge\. First, no universal decomposition method succeeds across all architectures: ICA achieves proper coverage \(89–91%\) for three models \(AIFS, FourCastNetV2, Aurora\) by exploiting non\-Gaussian structure, while SVD achieves coverage for only two models with severe underfitting for Aurora \(58% coverage\)\. Second, eigenvalue concentration determines discrimination capacity\. The variance collapse proposition shows that when correction terms consume\>90%\>90\\%of prior variance, the ability to distinguish tropical cyclones from routine weather fails\. ICA consistently produces higher coefficient of variation \(0\.07–1\.81\) than SVD \(0\.01–0\.49\), yielding adaptive intervals that scale with extreme event severity\.

These findings provide actionable guidance: practitioners should validate decomposition methods on held\-out extreme events before deployment, prioritizing coverage constraints over sharpness optimization alone\. The theoretical characterizations enable predictive diagnosis of UQ quality from architectural properties and data statistics\.

###### Acknowledgements\.

This research was supported by the Department of Education \(DepEd\), Philippines, under Department Order No\. 013, s\. 2025, which established the Education Center for AI Research \(ECAIR\), implemented through SEAMEO INNOTECH\. Code, calibration matrices, and the EM\-DAT date list are publicly released\.

## References

- A\. N\. Angelopoulos and S\. Bates \(2021\)A gentle introduction to conformal prediction and distribution\-free uncertainty quantification\.arXiv preprint arXiv:2107\.07511\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2107.07511)Cited by:[§1](https://arxiv.org/html/2606.02886#S1.p3.1),[§2](https://arxiv.org/html/2606.02886#S2.p2.1)\.
- K\. Bi, L\. Xie, H\. Zhang, X\. Chen, X\. Gu, and Q\. Tian \(2023\)Accurate medium\-range global weather forecasting with 3d neural networks\.Nature619\(7970\),pp\. 533–538\.External Links:[Document](https://dx.doi.org/10.1038/s41586-023-06185-3)Cited by:[§A\.1](https://arxiv.org/html/2606.02886#A1.SS1.p1.1),[§B\.2](https://arxiv.org/html/2606.02886#A2.SS2.SSS0.Px2.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p2.1),[§2](https://arxiv.org/html/2606.02886#S2.p1.1),[§4\.1](https://arxiv.org/html/2606.02886#S4.SS1.p1.1)\.
- C\. Blundell, J\. Cornebise, K\. Kavukcuoglu, and D\. Wierstra \(2015\)Weight uncertainty in neural networks\.InInternational Conference on Machine Learning,Lille, France,pp\. 1613–1622\.Cited by:[§A\.2](https://arxiv.org/html/2606.02886#A1.SS2.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p3.1),[§2](https://arxiv.org/html/2606.02886#S2.p2.1)\.
- C\. Bodnar, W\. P\. Bruinsma, A\. Lucic, M\. Stanley, A\. Allen, J\. Brandstetter, P\. Garvan, M\. Riechert, J\. A\. Weyn, H\. Dong, J\. K\. Gupta, K\. Thambiratnam, A\. T\. Archibald, C\. Wu, E\. Heider, M\. Welling, R\. E\. Turner, and P\. Perdikaris \(2025\)A foundation model for the Earth system\.Nature641,pp\. 1180–1187\.External Links:[Document](https://dx.doi.org/10.1038/s41586-025-09005-y)Cited by:[§A\.1](https://arxiv.org/html/2606.02886#A1.SS1.p1.1),[§B\.2](https://arxiv.org/html/2606.02886#A2.SS2.SSS0.Px2.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p2.1),[§2](https://arxiv.org/html/2606.02886#S2.p1.1),[§4\.1](https://arxiv.org/html/2606.02886#S4.SS1.p1.1)\.
- B\. Bonev, T\. Kurth, C\. Hundt, J\. Pathak, M\. Baust, K\. Kashinath, and A\. Anandkumar \(2023\)Spherical fourier neural operators: learning stable dynamics on the sphere\.InInternational Conference on Machine Learning,Honolulu, HI, USA,pp\. 2806–2823\.Cited by:[§4\.1](https://arxiv.org/html/2606.02886#S4.SS1.p1.1)\.
- J\. Cardoso \(1999\)High\-order contrasts for independent component analysis\.Neural Computation11\(1\),pp\. 157–192\.External Links:[Document](https://dx.doi.org/10.1162/089976699300016863)Cited by:[§B\.6](https://arxiv.org/html/2606.02886#A2.SS6.SSS0.Px1.p1.2),[§C\.2](https://arxiv.org/html/2606.02886#A3.SS2.p2.1)\.
- N\. A\. C\. Cressie \(1993\)Statistics for spatial data\.Revised edition,John Wiley & Sons\.External Links:[Document](https://dx.doi.org/10.1002/9781119115151)Cited by:[§3\.4](https://arxiv.org/html/2606.02886#S3.SS4.SSS0.Px3.p1.3)\.
- E\. Daxberger, A\. Kristiadi, A\. Immer, R\. Eschenhagen, M\. Bauer, and P\. Hennig \(2021\)Laplace redux – effortless bayesian deep learning\.InAdvances in Neural Information Processing Systems,Vol\.34,Red Hook, NY, USA,pp\. 20089–20103\.Cited by:[§2](https://arxiv.org/html/2606.02886#S2.p3.2)\.
- D\. Delforge, V\. Wathelet, R\. Below, C\. Lanfredi Sofia, M\. Tonnelier, J\. A\. F\. van Loenhout, and N\. Speybroeck \(2025\)EM\-DAT: the emergency events database\.International Journal of Disaster Risk Reduction124,pp\. 105509\.External Links:[Document](https://dx.doi.org/10.1016/j.ijdrr.2025.105509)Cited by:[§B\.2](https://arxiv.org/html/2606.02886#A2.SS2.SSS0.Px3.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p1.1),[§4\.2](https://arxiv.org/html/2606.02886#S4.SS2.p1.2)\.
- Y\. Gal and Z\. Ghahramani \(2016\)Dropout as a bayesian approximation: representing model uncertainty in deep learning\.InInternational Conference on Machine Learning,New York, NY, USA,pp\. 1050–1059\.Cited by:[§A\.2](https://arxiv.org/html/2606.02886#A1.SS2.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p3.1),[§2](https://arxiv.org/html/2606.02886#S2.p2.1)\.
- T\. Gneiting and A\. E\. Raftery \(2007\)Strictly proper scoring rules, prediction, and estimation\.Journal of the American Statistical Association102\(477\),pp\. 359–378\.External Links:[Document](https://dx.doi.org/10.1198/016214506000001437)Cited by:[§4\.3](https://arxiv.org/html/2606.02886#S4.SS3.SSS0.Px3.p1.2)\.
- A\. Graves \(2011\)Practical variational inference for neural networks\.InAdvances in Neural Information Processing Systems,Vol\.24,Red Hook, NY, USA,pp\. 2348–2356\.Cited by:[§A\.2](https://arxiv.org/html/2606.02886#A1.SS2.p1.1)\.
- C\. Guo, G\. Pleiss, Y\. Sun, and K\. Q\. Weinberger \(2017\)On calibration of modern neural networks\.InInternational Conference on Machine Learning,Sydney, Australia,pp\. 1321–1330\.Cited by:[§3\.5](https://arxiv.org/html/2606.02886#S3.SS5.p1.4)\.
- B\. He, B\. Lakshminarayanan, and Y\. W\. Teh \(2020\)Bayesian deep ensembles via the neural tangent kernel\.InAdvances in Neural Information Processing Systems,Vol\.33,Red Hook, NY, USA,pp\. 1010–1022\.Cited by:[Appendix C](https://arxiv.org/html/2606.02886#A3.SS0.SSS0.Px2.p1.11),[§1](https://arxiv.org/html/2606.02886#S1.p5.1),[§2](https://arxiv.org/html/2606.02886#S2.p3.2)\.
- H\. Hersbach, B\. Bell, P\. Berrisford, S\. Hirahara, A\. Horányi, J\. Muñoz\-Sabater, J\. Nicolas, C\. Peubey, R\. Radu, D\. Schepers,et al\.\(2020\)The era5 global reanalysis\.Quarterly Journal of the Royal Meteorological Society146\(730\),pp\. 1999–2049\.External Links:[Document](https://dx.doi.org/10.1002/qj.3803)Cited by:[§1](https://arxiv.org/html/2606.02886#S1.p7.1)\.
- Z\. Huang, H\. Lam, and H\. Zhang \(2023\)Efficient uncertainty quantification and reduction for over\-parameterized neural networks\.InAdvances in Neural Information Processing Systems,Vol\.36,Red Hook, NY, USA,pp\. 64428–64467\.Cited by:[Appendix C](https://arxiv.org/html/2606.02886#A3.SS0.SSS0.Px2.p1.11),[§1](https://arxiv.org/html/2606.02886#S1.p5.1),[§2](https://arxiv.org/html/2606.02886#S2.p3.2)\.
- A\. Hyvärinen and E\. Oja \(2000\)Independent component analysis: algorithms and applications\.Neural Networks13\(4\-5\),pp\. 411–430\.External Links:[Document](https://dx.doi.org/10.1016/S0893-6080%2800%2900026-5)Cited by:[§B\.6](https://arxiv.org/html/2606.02886#A2.SS6.SSS0.Px1.p1.2),[§B\.6](https://arxiv.org/html/2606.02886#A2.SS6.SSS0.Px2.p1.3),[§C\.2](https://arxiv.org/html/2606.02886#A3.SS2.p2.1),[§3\.4](https://arxiv.org/html/2606.02886#S3.SS4.SSS0.Px2.p1.6)\.
- A\. Hyvärinen \(1999\)Fast and robust fixed\-point algorithms for independent component analysis\.IEEE Transactions on Neural Networks10\(3\),pp\. 626–634\.External Links:[Document](https://dx.doi.org/10.1109/72.761722)Cited by:[§B\.6](https://arxiv.org/html/2606.02886#A2.SS6.SSS0.Px1.p2.7),[§3\.4](https://arxiv.org/html/2606.02886#S3.SS4.SSS0.Px2.p1.6)\.
- A\. Jacot, F\. Gabriel, and C\. Hongler \(2018\)Neural tangent kernel: convergence and generalization in neural networks\.InAdvances in Neural Information Processing Systems,Vol\.31,Red Hook, NY, USA,pp\. 8571–8580\.Cited by:[Appendix C](https://arxiv.org/html/2606.02886#A3.SS0.SSS0.Px2.p1.1),[§2](https://arxiv.org/html/2606.02886#S2.p3.2)\.
- A\. Kendall and Y\. Gal \(2017\)What uncertainties do we need in Bayesian deep learning for computer vision?\.InAdvances in Neural Information Processing Systems,Vol\.30,Red Hook, NY, USA,pp\. 5580–5590\.Cited by:[§3\.2](https://arxiv.org/html/2606.02886#S3.SS2.SSS0.Px1.p1.5)\.
- B\. Lakshminarayanan, A\. Pritzel, and C\. Blundell \(2017\)Simple and scalable predictive uncertainty estimation using deep ensembles\.InAdvances in Neural Information Processing Systems,Vol\.30,Red Hook, NY, USA,pp\. 6402–6413\.Cited by:[§A\.2](https://arxiv.org/html/2606.02886#A1.SS2.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p3.1),[§2](https://arxiv.org/html/2606.02886#S2.p2.1)\.
- R\. Lam, A\. Sanchez\-Gonzalez, M\. Willson, P\. Wirnsberger, M\. Fortunato, F\. Alet, S\. Ravuri, T\. Ewalds, Z\. Eaton\-Rosen, W\. Hu,et al\.\(2023\)Learning skillful medium\-range global weather forecasting\.Science382\(6677\),pp\. 1416–1421\.External Links:[Document](https://dx.doi.org/10.1126/science.adi2336)Cited by:[§A\.1](https://arxiv.org/html/2606.02886#A1.SS1.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p2.1),[§2](https://arxiv.org/html/2606.02886#S2.p1.1)\.
- S\. Lang, M\. Alexe, M\. Chantry, J\. Dramsch, F\. Pinault, B\. Raoult, M\. Clare, C\. Lessig, M\. Maier\-Gerber,et al\.\(2024\)AIFS – ecmwf’s data\-driven forecasting system\.arXiv preprint arXiv:2406\.01465\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2406.01465)Cited by:[§B\.2](https://arxiv.org/html/2606.02886#A2.SS2.SSS0.Px2.p1.1),[§4\.1](https://arxiv.org/html/2606.02886#S4.SS1.p1.1)\.
- J\. Lee, Y\. Bahri, R\. Novak, S\. S\. Schoenholz, J\. Pennington, and J\. Sohl\-Dickstein \(2018\)Deep neural networks as Gaussian processes\.InInternational Conference on Learning Representations,Cited by:[§3\.2](https://arxiv.org/html/2606.02886#S3.SS2.SSS0.Px1.p1.5)\.
- L\. Li, R\. Carver, I\. Lopez\-Gomez, F\. Sha, and J\. Anderson \(2024\)Generative emulation of weather forecast ensembles with diffusion models\.Science Advances10\(13\),pp\. eadk4489\.External Links:[Document](https://dx.doi.org/10.1126/sciadv.adk4489)Cited by:[§2](https://arxiv.org/html/2606.02886#S2.p1.1)\.
- D\. J\. MacKay \(1992\)A practical bayesian framework for backpropagation networks\.Neural Computation4\(3\),pp\. 448–472\.Cited by:[§2](https://arxiv.org/html/2606.02886#S2.p3.2)\.
- K\. V\. Mardia \(1970\)Measures of multivariate skewness and kurtosis with applications\.Biometrika57\(3\),pp\. 519–530\.External Links:[Document](https://dx.doi.org/10.1093/biomet/57.3.519)Cited by:[2nd item](https://arxiv.org/html/2606.02886#A3.I1.i2.p1.1)\.
- R\. Newman and I\. Noy \(2023\)The global costs of extreme weather that are attributable to climate change\.Nature Communications14\(1\),pp\. 6103\.External Links:[Document](https://dx.doi.org/10.1038/s41467-023-41888-1)Cited by:[§1](https://arxiv.org/html/2606.02886#S1.p1.1)\.
- Y\. Ovadia, E\. Fertig, J\. Ren, Z\. Nado, D\. Sculley, S\. Nowozin, J\. V\. Dillon, B\. Lakshminarayanan, and J\. Snoek \(2019\)Can you trust your model’s uncertainty? evaluating predictive uncertainty under dataset shift\.InAdvances in Neural Information Processing Systems,Vol\.32,Red Hook, NY, USA,pp\. 13969–13980\.Cited by:[§A\.2](https://arxiv.org/html/2606.02886#A1.SS2.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p3.1),[§2](https://arxiv.org/html/2606.02886#S2.p2.1)\.
- J\. Pathak, S\. Subramanian, P\. Harrington, S\. Raja, A\. Chattopadhyay, M\. Mardani, T\. Kurth, D\. Hall, Z\. Li, K\. Azizzadenesheli, P\. Hassanzadeh, K\. Kashinath, and A\. Anandkumar \(2022\)FourCastNet: a global data\-driven high\-resolution weather model using adaptive fourier neural operators\.arXiv preprint arXiv:2202\.11214\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2202.11214)Cited by:[§A\.1](https://arxiv.org/html/2606.02886#A1.SS1.p1.1),[§B\.2](https://arxiv.org/html/2606.02886#A2.SS2.SSS0.Px2.p1.1),[§1](https://arxiv.org/html/2606.02886#S1.p2.1),[§2](https://arxiv.org/html/2606.02886#S2.p1.1),[§4\.1](https://arxiv.org/html/2606.02886#S4.SS1.p1.1)\.
- I\. Price, A\. Sanchez\-Gonzalez, F\. Alet, T\. R\. Andersson, A\. El\-Kadi, D\. Masters, T\. Ewalds, J\. Stott, S\. Mohamed, P\. Battaglia, R\. Lam, and M\. Willson \(2025\)GenCast: diffusion\-based ensemble weather forecasting at scale\.Nature637,pp\. 84–90\.External Links:[Document](https://dx.doi.org/10.1038/s41586-024-08252-9)Cited by:[§2](https://arxiv.org/html/2606.02886#S2.p1.1)\.
- S\. Rasp, S\. Hoyer, A\. Merose, I\. Langmore, P\. Battaglia, T\. Russell, A\. Sanchez\-Gonzalez, V\. Yang, R\. Carver, S\. Agrawal,et al\.\(2024\)WeatherBench 2: a benchmark for the next generation of data\-driven global weather models\.Journal of Advances in Modeling Earth Systems16\(6\),pp\. e2023MS004019\.External Links:[Document](https://dx.doi.org/10.1029/2023MS004019)Cited by:[§B\.2](https://arxiv.org/html/2606.02886#A2.SS2.SSS0.Px1.p1.1),[§4\.2](https://arxiv.org/html/2606.02886#S4.SS2.p1.2)\.
- J\. J\. Thiagarajan, R\. Anirudh, V\. Narayanaswamy, and P\. Bremer \(2022\)Single model uncertainty estimation via stochastic data centering\.InAdvances in Neural Information Processing Systems,Vol\.35,Red Hook, NY, USA,pp\. 25967–25981\.Cited by:[§2](https://arxiv.org/html/2606.02886#S2.p3.2)\.
- K\. Tran, W\. Neiswanger, J\. Yoon, Q\. Zhang, E\. Xing, and Z\. W\. Ulissi \(2020\)Methods for comparing uncertainty quantifications for material property predictions\.Machine Learning: Science and Technology1\(2\),pp\. 025006\.External Links:[Document](https://dx.doi.org/10.1088/2632-2153/ab7e1a)Cited by:[§4\.3](https://arxiv.org/html/2606.02886#S4.SS3.SSS0.Px4.p1.6)\.
- B\. Yu \(1994\)Rates of convergence for empirical processes of stationary mixing sequences\.Annals of Probability22\(1\),pp\. 94–116\.External Links:[Document](https://dx.doi.org/10.1214/aop/1176988849)Cited by:[§C\.2](https://arxiv.org/html/2606.02886#A3.SS2.SSS0.Px4.p1.1)\.

## Appendix AExtended Related Work

### A\.1\.AI Weather Foundation Models

Deep learning weather models have achieved competitive accuracy with numerical weather prediction while offering orders\-of\-magnitude speedups\. FourCastNet\(Pathaket al\.,[2022](https://arxiv.org/html/2606.02886#bib.bib4)\)pioneered the use of vision transformers for global weather forecasting at 0\.25° resolution, demonstrating that data\-driven models can match ECMWF’s Integrated Forecast System on many variables\. Pangu\-Weather\(Biet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib6)\)introduced a 3D Swin Transformer with hierarchical patch merging, achieving state\-of\-the\-art scores on multiple benchmarks\. GraphCast\(Lamet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib7)\)employed graph neural networks on an icosahedral mesh representation, winning multiple WeatherBench2 categories\. Aurora\(Bodnaret al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib8)\)extended this paradigm with a Perceiver architecture handling heterogeneous input sources and variable atmospheric conditions\. These models share common design principles: pretraining on decades of ERA5 reanalysis data, autoregressive rollout for multi\-step forecasting, and deterministic outputs\. The absence of uncertainty estimates limits their applicability for high\-stakes decision\-making in extreme weather scenarios\.

### A\.2\.Comparison with Existing UQ Methods

Deep ensembles\(Lakshminarayananet al\.,[2017](https://arxiv.org/html/2606.02886#bib.bib13)\)train multiple independent networks with different random initializations, using prediction variance as uncertainty\. This captures both epistemic uncertainty \(model disagreement\) and aleatoric uncertainty \(inherent stochasticity\)\. However, for billion\-parameter weather models with week\-long training times, even modest ensembles \(M=5M=5\) require impractical computational resources\. Bayesian neural networks\(Blundellet al\.,[2015](https://arxiv.org/html/2606.02886#bib.bib18); Graves,[2011](https://arxiv.org/html/2606.02886#bib.bib25)\)place distributions over weights but face similar scalability challenges\. Monte Carlo dropout\(Gal and Ghahramani,[2016](https://arxiv.org/html/2606.02886#bib.bib17)\)approximates Bayesian inference via stochastic forward passes but requires dropout layers \(incompatible with many pretrained architectures\) and produces poorly calibrated uncertainties on complex regression tasks\(Ovadiaet al\.,[2019](https://arxiv.org/html/2606.02886#bib.bib26)\)\. Temperature scaling and Platt scaling adjust output distributions but do not provide error\-correlated uncertainties\.

## Appendix BImplementation Details

### B\.1\.Feature Extraction

For each model, forward hooks register on the last layer to capture activations during inference\. The hook captures the output tensor before global average pooling is applied, then aggregates to produce a fixed\-size feature vector using multi\-statistic pooling \(mean, std, min, max, q25, q75 per channel\)\.

#### Per\-Architecture Details\.

Feature extraction differs by architecture\. For FourCastNetV2, a hook captures the last Spherical Fourier Neural Operator \(SFNO\) block output with shape\(B,256,H,W\)\(B,256,H,W\), where 256 is the channel dimension\. Multi\-statistic pooling across spatial dimensions yieldsd=256×6=1536d=256\\times 6=1536features\. For Pangu\-Weather, the ONNX model’s 69\-channel prediction tensor is pooled via global averaging tod=69d=69dimensions\. For Aurora, a hook on the Perceiver decoder captures the latent representation with shape\(B,2,65\)\(B,2,65\)for two time steps; global averaging yieldsd=65d=65features\. For AIFS, global\-average pooling over the final graph neural network layer yieldsd=1024d=1024features\.

### B\.2\.Dataset Details

#### ERA5 Reanalysis\.

All experiments use the WeatherBench2\(Raspet al\.,[2024](https://arxiv.org/html/2606.02886#bib.bib9)\)ERA5 dataset at 0\.25° resolution \(721×1440721\\times 1440grid points\), spanning 1959–2021 with 6\-hourly temporal resolution\. ERA5 is a reanalysis product, not direct observations: it is produced by assimilating historical observations into a numerical weather model, yielding a physically consistent but model\-dependent gridded dataset\. Reanalysis errors are generally small for well\-observed variables \(temperature, geopotential\) but may be larger for quantities with sparse observational coverage \(humidity, polar regions\)\. Following standard practice in AI weather model evaluation\(Raspet al\.,[2024](https://arxiv.org/html/2606.02886#bib.bib9)\), ERA5 is treated as ground truth throughout\.

#### Training Data Overlap\.

The four evaluated models were trained on overlapping subsets of ERA5: FourCastNetV2 on 1979–2015\(Pathaket al\.,[2022](https://arxiv.org/html/2606.02886#bib.bib4)\), Pangu\-Weather on 1979–2017\(Biet al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib6)\), Aurora on 1979–2020 \(including fine\-tuning\)\(Bodnaret al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib8)\), and AIFS on 1979–2020\(Langet al\.,[2024](https://arxiv.org/html/2606.02886#bib.bib3)\)\. The year 2021 falls outside all models’ training and fine\-tuning periods, ensuring that the evaluation dates represent genuinely unseen data for every model\.

#### Detailed Date Selection Methodology\.

The extreme weather events dataset comprises initialization dates from 2021, the only year in WeatherBench2 falling outside all models’ training periods and thus representing out\-of\-distribution temporal evaluation\. Dates are drawn from the EM\-DAT International Disaster Database\(Delforgeet al\.,[2025](https://arxiv.org/html/2606.02886#bib.bib44)\)and stratified across hazard types to ensure diverse scenarios\. Multiple disasters often occur simultaneously across different regions, allowing a single initialization date to capture several concurrent extreme events\.

For each disaster event, the initialization date is selected using a 3\-day lookback period from the reported landfall or onset date\. Since ERA5 provides data at 6\-hour intervals \(00:00, 06:00, 12:00, 18:00 UTC\), this corresponds to 12 timesteps or 72 hours prior to the event peak\. This lookback ensures that forecast initialization occurs during the event’s development phase rather than after landfall, capturing the operational forecasting scenario where prediction uncertainty is most critical for disaster preparedness\.

#### Event\-Type Distribution\.

Table[4](https://arxiv.org/html/2606.02886#A2.T4)reports the distribution of the 2021 EM\-DAT events used for calibration, by disaster type\. At the time of data filtering, the meteorologically relevant subset comprised 206 verified disaster events across 82 countries: 136 floods, 63 storms \(including tropical cyclones Tauktae, Ida, Rai, and Elsa\), 5 droughts, and 2 extreme\-temperature events \(the June 2021 Pacific Northwest heat wave\)\. These deduplicate to 100 distinct initialization dates\. Floods and storms dominate, consistent with their global frequency among hydro\-meteorological hazards\.

Table 4\.Distribution of 2021 EM\-DAT events used for calibration, by disaster type, as captured at the time of data filtering\. The 206 events deduplicate to 100 distinct initialization dates across 82 countries\. Counts are a snapshot of a continuously updated database \(see notes\)\.Counts reflect the EM\-DAT snapshot at filtering time\. EM\-DAT is a living database: records are revised and historical events are added retrospectively, so a later download yields different totals \(e\.g\., a subsequent snapshot reported 362 events for 2021\)\. The post\-hoc design allows recalibration on any updated or future snapshot without retraining\.

Two properties of this dataset warrant emphasis\. First, the 206 events deduplicate to only 100 distinct initialization dates\. Because multiple extreme events frequently co\-occur on the same calendar dates across different regions, several concurrent disasters can map to a single initialization date\. A given date in the calibration set may therefore represent simultaneous hydro\-meteorological hazards in distinct parts of the globe, and the event\-type counts in Table[4](https://arxiv.org/html/2606.02886#A2.T4)sum to more than the number of unique dates by construction\. Second, the counts are a snapshot\. EM\-DAT is continuously curated, with disaster records revised and late\-reported 2021 events added over time, so the distribution reported here reflects the database state at filtering time and a future re\-download would yield somewhat different totals and shares\. Because NTK\-UQ is purely post\-hoc, recalibrating on an updated or expanded event set requires no model retraining, only a re\-run of the offline calibration stage\.

### B\.3\.Binary Search Calibration Algorithm

The scaling factorα\\alphais found via binary search on the calibration set\. Given errors\{ei\}\\\{e\_\{i\}\\\}and raw uncertainties\{σi\}\\\{\\sigma\_\{i\}\\\}, the algorithm findsα\\alphasuch that:

\(11\)1N​∑i=1N𝟏​\[\|ei\|<z⋅α⋅σi\]≈0\.90,\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\mathbf\{1\}\\left\[\|e\_\{i\}\|<z\\cdot\\alpha\\cdot\\sigma\_\{i\}\\right\]\\approx 0\.90,wherez≈1\.645z\\approx 1\.645is the standard normal quantile at probability 0\.95 \(for a two\-sided 90% prediction interval\)\. The binary search maintains bounds\[αlow,αhigh\]\[\\alpha\_\{\\text\{low\}\},\\alpha\_\{\\text\{high\}\}\]initialized to\[0\.1,100\]\[0\.1,100\]and iteratively refines the interval until the empirical coverage converges to the target within toleranceϵ=0\.01\\epsilon=0\.01\. Typically, convergence occurs within 10\-15 iterations\.

### B\.4\.Autoregressive Feature Extraction Implementation

AI weather models generate forecasts autoregressively:y^t\+τ=fθ\(τ/Δ​t\)​\(xt\)\\hat\{y\}\_\{t\+\\tau\}=f\_\{\\theta\}^\{\(\\tau/\\Delta t\)\}\(x\_\{t\}\)whereΔ​t\\Delta tis the model’s native time step and the superscript denotes iterated application\. For a set of target horizons𝒯=\{6,12,24,48,72,120\}\\mathcal\{T\}=\\\{6,12,24,48,72,120\\\}hours, features at each horizonτ∈𝒯\\tau\\in\\mathcal\{T\}are collected as:

\(12\)ϕτ=ϕ​\(fθ\(τ/Δ​t\)​\(xt\)\),\\phi\_\{\\tau\}=\\phi\\left\(f\_\{\\theta\}^\{\(\\tau/\\Delta t\)\}\(x\_\{t\}\)\\right\),whereϕ​\(⋅\)\\phi\(\\cdot\)denotes last\-layer feature extraction\. This yields a collection\{ϕτ\}τ∈𝒯\\\{\\phi\_\{\\tau\}\\\}\_\{\\tau\\in\\mathcal\{T\}\}from a single rollout, reducing computational cost by a factor of\|𝒯\|=6\|\\mathcal\{T\}\|=6compared to independent forward passes to each horizon\.

### B\.5\.Computational Requirements

Calibration costs depend on model complexity and GPU hardware\. With multi\-lead\-time rollout extracting 6 lead times per sample, calibration time scales linearly with sample count\. FourCastNetV2 processes approximately 7 samples per hour on an L4 GPU \(24GB\)\. Pangu\-Weather and Aurora require an A100 GPU \(40GB\) due to higher memory requirements, processing approximately 5 samples per hour\.

Inference overhead is minimal: computing uncertainty for a single sample requires one matrix\-vector productϕ⋅Vk\\phi\\cdot V\_\{k\}\(sized×kd\\times k, withk≤50k\\leq 50,d≤1536d\\leq 1536\) followed by element\-wise operations—negligible relative to the model forward pass\. Storage requirements are approximately 50MB per model \(right singular vectorsVkV\_\{k\}and eigenvaluesΛk\\Lambda\_\{k\}for each calibrated lead time\)\.

### B\.6\.ICA Theory: Non\-Gaussian Discrimination

###### Proposition 1 \(Non\-Gaussian Discrimination \(Formal\)\)\.

Letϕ​\(x\)∈ℝd\\phi\(x\)\\in\\mathbb\{R\}^\{d\}denote last\-layer features withϕ∼P\\phi\\sim P\. SVD decomposes centered featuresΦ~\\tilde\{\\Phi\}by maximizing varianceVar​\(vj⊤​ϕ~\)\\mathrm\{Var\}\(v\_\{j\}^\{\\top\}\\tilde\{\\phi\}\)subject to orthogonality, yielding principal components ordered by decreasingλj=Var​\(vj⊤​ϕ~\)\\lambda\_\{j\}=\\mathrm\{Var\}\(v\_\{j\}^\{\\top\}\\tilde\{\\phi\}\)\. ICA instead maximizes statistical independence, findingS=W​Φ~whiteS=W\\tilde\{\\Phi\}\_\{\\mathrm\{white\}\}such that∑j=1kI​\(Sj;S−j\)\\sum\_\{j=1\}^\{k\}I\(S\_\{j\};S\_\{\-j\}\)is minimized\. When feature distributions exhibit joint non\-Gaussianity \(higher\-order cumulantsκi1,…,im≠0\\kappa\_\{i\_\{1\},\\ldots,i\_\{m\}\}\\neq 0form≥3m\\geq 3\), SVD captures only second\-order structure \(the covarianceΣ=𝔼​\[ϕ~​ϕ~⊤\]\\Sigma=\\mathbb\{E\}\[\\tilde\{\\phi\}\\tilde\{\\phi\}^\{\\top\}\]\), discarding higher\-order information, whereas ICA exploits kurtosis and negentropy via contrast functionsJ​\(s\)=𝔼​\[G​\(s\)\]J\(s\)=\\mathbb\{E\}\[G\(s\)\]\(e\.g\.,G​\(s\)=s4G\(s\)=s^\{4\}\)\. For extreme events with heavy\-tailed marginals \(𝔼​\[s4\]≫3​σ4\\mathbb\{E\}\[s^\{4\}\]\\gg 3\\sigma^\{4\}\), ICA’s leading components aligned with extreme directions achieve higher kurtosis than PCA’s leading \(high\-variance\) components, which are dominated by typical synoptic modes of large variance but low kurtosis\. This enables discrimination that SVD’s top\-kkcomponents cannot achieve when extreme event signals carry low variance relative to background variability\.

#### Detailed Justification\.

SVD decomposes features via eigendecomposition of the covariance matrixΣ\\Sigma, which captures only pairwise correlations \(second\-order statistics\)\. For Gaussian data,Σ\\Sigmafully characterizes the distribution\. For non\-Gaussian data, higher\-order moments \(third, fourth, etc\.\) carry essential information about the distribution’s shape, tail behavior, and dependencies\(Hyvärinen and Oja,[2000](https://arxiv.org/html/2606.02886#bib.bib45); Cardoso,[1999](https://arxiv.org/html/2606.02886#bib.bib47)\)\. By ignoring these higher\-order statistics, SVD produces principal components that maximize variance but fail to isolate directions of non\-Gaussian extreme events\.

ICA, by contrast, finds a linear transformationWWsuch that transformed featuresS=W​ΦwhiteS=W\\Phi\_\{\\mathrm\{white\}\}are maximally independent\. Independence is a*stronger*condition than decorrelation: while decorrelation \(enforced by PCA whitening\) ensures𝔼​\[Si​Sj\]=0\\mathbb\{E\}\[S\_\{i\}S\_\{j\}\]=0, independence requires𝔼​\[g​\(Si\)​h​\(Sj\)\]=𝔼​\[g​\(Si\)\]​𝔼​\[h​\(Sj\)\]\\mathbb\{E\}\[g\(S\_\{i\}\)h\(S\_\{j\}\)\]=\\mathbb\{E\}\[g\(S\_\{i\}\)\]\\mathbb\{E\}\[h\(S\_\{j\}\)\]for all functionsg,hg,h\. FastICA\(Hyvärinen,[1999](https://arxiv.org/html/2606.02886#bib.bib46)\)achieves this by maximizing negentropyJ​\(s\)=H​\(sGauss\)−H​\(s\)J\(s\)=H\(s\_\{\\mathrm\{Gauss\}\}\)\-H\(s\)whereHHis differential entropy\. This criterion is sensitive to non\-Gaussianity: sources with high kurtosis \(heavy tails\) yield high negentropy, guiding ICA to isolate extreme event directions\.

#### ICA Identifiability\.

A fundamental result from ICA theory\(Hyvärinen and Oja,[2000](https://arxiv.org/html/2606.02886#bib.bib45)\)states that if observed features arise from a linear mixtureϕ=A​s\\phi=Aswheressare independent sources with*at most one Gaussian component*, then the mixing matrixAAis identifiable up to permutation and scaling\. Critically,ICA fails if all sources are Gaussian, because any rotation of jointly Gaussian variables remains Gaussian with the same likelihood\. The identifiability theorem thus requires non\-Gaussian sources, which in this context correspond to physical drivers of extreme events: vorticity \(heavy\-tailed during tropical cyclones\), moisture advection \(bimodal during atmospheric rivers\), and diabatic heating \(positively skewed during convective extremes\)\. By assumption, neural weather models learn to encode these non\-Gaussian physical processes in their feature representations\. ICA recovers these sources by unmixing the learned features, enabling adaptive uncertainty that scales with the magnitude of extreme event drivers\.

## Appendix CTheoretical Guarantees

This section establishes the theoretical foundations of NTK\-UQ\. Full proofs are provided in the supplementary material\.

#### Architecture\-Agnostic Last\-Layer Kernel\.

The NTK\-UQ framework applies to any neural network architecture that admits a decompositionfθ=gψ∘ϕωf\_\{\\theta\}=g\_\{\\psi\}\\circ\\phi\_\{\\omega\}whereϕω:𝒳→ℝd\\phi\_\{\\omega\}:\\mathcal\{X\}\\to\\mathbb\{R\}^\{d\}is a feature extractor andgψ:ℝd→𝒴g\_\{\\psi\}:\\mathbb\{R\}^\{d\}\\to\\mathcal\{Y\}is the prediction head\. Given such a decomposition, the last\-layer kernelK​\(x,x′\)=ϕω​\(x\)⊤​ϕω​\(x′\)K\(x,x^\{\\prime\}\)=\\phi\_\{\\omega\}\(x\)^\{\\top\}\\phi\_\{\\omega\}\(x^\{\\prime\}\)is a valid positive semi\-definite kernel by construction, since for any set of points\{xi\}i=1n\\\{x\_\{i\}\\\}\_\{i=1\}^\{n\}, the Gram matrixKi​j=ϕ​\(xi\)⊤​ϕ​\(xj\)=Φ​Φ⊤K\_\{ij\}=\\phi\(x\_\{i\}\)^\{\\top\}\\phi\(x\_\{j\}\)=\\Phi\\Phi^\{\\top\}whereΦ∈ℝn×d\\Phi\\in\\mathbb\{R\}^\{n\\times d\}is the feature matrix\. This holds regardless of the internal structure ofϕω\\phi\_\{\\omega\}, whether it is an SFNO \(FourCastNetV2\), Swin Transformer \(Pangu\-Weather\), or Perceiver \(Aurora\)\.

#### Predictive Variance\.

Under the NTK\-GP correspondence\(Jacotet al\.,[2018](https://arxiv.org/html/2606.02886#bib.bib11)\), the GP posterior variance at test pointx∗x\_\{\*\}is:

\(13\)σ2​\(x∗\)=‖ϕ~​\(x∗\)‖2\+σn2−∑j=1kλj​cj2λj\+σn2\\begin\{split\}\\sigma^\{2\}\(x\_\{\*\}\)&=\\\|\\tilde\{\\phi\}\(x\_\{\*\}\)\\\|^\{2\}\+\\sigma\_\{n\}^\{2\}\\\\ &\\quad\-\\sum\\nolimits\_\{j=1\}^\{k\}\\frac\{\\lambda\_\{j\}c\_\{j\}^\{2\}\}\{\\lambda\_\{j\}\+\\sigma\_\{n\}^\{2\}\}\\end\{split\}wherecj=ϕ~​\(x∗\)⊤​vjc\_\{j\}=\\tilde\{\\phi\}\(x\_\{\*\}\)^\{\\top\}v\_\{j\}are the projection coefficients onto the right singular vectorsVVofΦ~=U​S​V⊤\\tilde\{\\Phi\}=USV^\{\\top\},λj=sj2\\lambda\_\{j\}=s\_\{j\}^\{2\}are the squared singular values, andσn2\\sigma\_\{n\}^\{2\}is estimated from the eigenvalue tail \(Section[D\.2](https://arxiv.org/html/2606.02886#A4.SS2)\)\. For centered data, SVD ofΦ~\\tilde\{\\Phi\}and PCA eigendecomposition ofΣ=Φ~⊤​Φ~\\Sigma=\\tilde\{\\Phi\}^\{\\top\}\\tilde\{\\Phi\}are equivalent; SVD is used throughout for numerical stability\. The\+σn2\+\\sigma\_\{n\}^\{2\}term ensures the predictive variance never drops below the noise floor\. The derivation follows from the Woodbury identity applied toσ2=k​\(x∗,x∗\)\+σn2−k∗⊤​\(K\+σn2​I\)−1​k∗\\sigma^\{2\}=k\(x\_\{\*\},x\_\{\*\}\)\+\\sigma\_\{n\}^\{2\}\-k\_\{\*\}^\{\\top\}\(K\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}k\_\{\*\}withK=Φ~​Φ~⊤K=\\tilde\{\\Phi\}\\tilde\{\\Phi\}^\{\\top\}; application to finite\-width networks is justified by He et al\.\(Heet al\.,[2020](https://arxiv.org/html/2606.02886#bib.bib12)\)and Huang et al\.\(Huanget al\.,[2023](https://arxiv.org/html/2606.02886#bib.bib32)\)\.

#### Approximations\.

The framework relies on: \(1\) the infinite\-width NTK approximation, which holds approximately for wide networks and is corrected by post\-hoc calibration; and \(2\) last\-layer\-only features, which may underestimate total epistemic uncertainty\. Empirical validation on held\-out data confirms these approximations are acceptable in practice\.

#### Rank Selection and Variance Collapse \(Proof of Proposition[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)\)\.

The correction isCk=∑j=1kwj​cj2C\_\{k\}=\\sum\_\{j=1\}^\{k\}w\_\{j\}c\_\{j\}^\{2\}wherewj=λj/\(λj\+σn2\)w\_\{j\}=\\lambda\_\{j\}/\(\\lambda\_\{j\}\+\\sigma\_\{n\}^\{2\}\)\.

###### Proof\.

Whenσn2\>0\\sigma\_\{n\}^\{2\}\>0: sincewj<1w\_\{j\}<1, we haveCk<∑j=1kcj2=‖ϕ~​\(x∗\)‖2−‖ϕ~⟂‖2C\_\{k\}<\\sum\_\{j=1\}^\{k\}c\_\{j\}^\{2\}=\\\|\\tilde\{\\phi\}\(x\_\{\*\}\)\\\|^\{2\}\-\\\|\\tilde\{\\phi\}\_\{\\perp\}\\\|^\{2\}\(Pythagorean identity\), soσ2​\(x∗\)\>‖ϕ~⟂‖2\+σn2\>0\\sigma^\{2\}\(x\_\{\*\}\)\>\\\|\\tilde\{\\phi\}\_\{\\perp\}\\\|^\{2\}\+\\sigma\_\{n\}^\{2\}\>0\. No collapse occurs\. Whenσn2=0\\sigma\_\{n\}^\{2\}=0:wj=1w\_\{j\}=1for allλj\>0\\lambda\_\{j\}\>0, soCr=∑j=1rcj2=‖ϕ~​\(x∗\)‖2−‖ϕ~⟂‖2C\_\{r\}=\\sum\_\{j=1\}^\{r\}c\_\{j\}^\{2\}=\\\|\\tilde\{\\phi\}\(x\_\{\*\}\)\\\|^\{2\}\-\\\|\\tilde\{\\phi\}\_\{\\perp\}\\\|^\{2\}, givingσ2​\(x∗\)=‖ϕ~⟂‖2\\sigma^\{2\}\(x\_\{\*\}\)=\\\|\\tilde\{\\phi\}\_\{\\perp\}\\\|^\{2\}\. Whenn≥dn\\geq d, the orthogonal residualϕ~⟂=0\\tilde\{\\phi\}\_\{\\perp\}=0andσ2​\(x∗\)→0\\sigma^\{2\}\(x\_\{\*\}\)\\to 0\. ∎

A useful rank selection heuristic is∑j=1kλj/∑j=1dλj≥0\.99\\sum\_\{j=1\}^\{k\}\\lambda\_\{j\}/\\sum\_\{j=1\}^\{d\}\\lambda\_\{j\}\\geq 0\.99\. For concentrated spectra \(SFNO\), this requiresk=2k=2–1010; for distributed spectra \(ViT, Perceiver\),kkmay equal the full rank\.

### C\.1\.Architecture\-Dependent Spectral Structure

The architecture\-dependent UQ behavior reduces to a single mechanism: the decay rate of the feature covariance spectrum controls the effective rank, and hence \(via Proposition[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)\) the truncation budget before variance collapse\. We first establish this mechanism unconditionally, then show how each architecture’s representation geometry induces the relevant decay\.

###### Lemma 1 \(Spectral Decay Bounds Effective Rank\)\.

LetΣ=𝔼​\[ϕ~​\(x\)​ϕ~​\(x\)⊤\]\\Sigma=\\mathbb\{E\}\[\\tilde\{\\phi\}\(x\)\\tilde\{\\phi\}\(x\)^\{\\top\}\]be the centered feature covariance with eigenvaluesλ1≥λ2≥⋯≥0\\lambda\_\{1\}\\geq\\lambda\_\{2\}\\geq\\cdots\\geq 0and total energyT=∑jλjT=\\sum\_\{j\}\\lambda\_\{j\}, and letrα=min⁡\{k:∑j=1kλj≥α​T\}r\_\{\\alpha\}=\\min\\\{k:\\sum\_\{j=1\}^\{k\}\\lambda\_\{j\}\\geq\\alpha T\\\}be the effective rank at thresholdα∈\(0,1\)\\alpha\\in\(0,1\)\. If the eigenvalues decay polynomially,λk≤C​k−β\\lambda\_\{k\}\\leq C\\,k^\{\-\\beta\}withβ\>1\\beta\>1, then

\(14\)rα≤⌈\(C\(β−1\)​\(1−α\)​T\)1/\(β−1\)⌉=O​\(\(1−α\)−1/\(β−1\)\),r\_\{\\alpha\}\\;\\leq\\;\\left\\lceil\\left\(\\frac\{C\}\{\(\\beta\-1\)\(1\-\\alpha\)\\,T\}\\right\)^\{1/\(\\beta\-1\)\}\\right\\rceil=O\\\!\\big\(\(1\-\\alpha\)^\{\-1/\(\\beta\-1\)\}\\big\),*independent of the ambient feature dimensiondd\.*

###### Proof\.

Sincej−βj^\{\-\\beta\}is decreasing, the tail energy satisfies

\(15\)∑j\>kλj≤C​∑j\>kj−β≤C​∫k∞t−β​𝑑t=C​k−\(β−1\)β−1\.\\sum\_\{j\>k\}\\lambda\_\{j\}\\;\\leq\\;C\\sum\_\{j\>k\}j^\{\-\\beta\}\\;\\leq\\;C\\int\_\{k\}^\{\\infty\}t^\{\-\\beta\}\\,dt\\;=\\;\\frac\{C\\,k^\{\-\(\\beta\-1\)\}\}\{\\beta\-1\}\.The thresholdrαr\_\{\\alpha\}is attained once the retained fraction reachesα\\alpha, i\.e\. once the tail∑j\>kλj≤\(1−α\)​T\\sum\_\{j\>k\}\\lambda\_\{j\}\\leq\(1\-\\alpha\)T\. It therefore suffices thatC​k−\(β−1\)/\(β−1\)≤\(1−α\)​TCk^\{\-\(\\beta\-1\)\}/\(\\beta\-1\)\\leq\(1\-\\alpha\)T, which rearranges tok≥\(C/\[\(β−1\)​\(1−α\)​T\]\)1/\(β−1\)k\\geq\\big\(C/\[\(\\beta\-1\)\(1\-\\alpha\)T\]\\big\)^\{1/\(\\beta\-1\)\}, giving \([14](https://arxiv.org/html/2606.02886#A3.E14)\)\. The bound depends onβ,C,T\\beta,C,Tbut not ondd\. ∎

Lemma[1](https://arxiv.org/html/2606.02886#A3.Thmtheorem1)makes the qualitative claim precise:*faster decay \(largerβ\\beta\) yields smaller effective rank, with no dependence on the ambient dimension\.*The architecture enters only through the decay exponentβ\\beta, which we now characterize per family under one explicit, empirically checkable hypothesis each\.

###### Proposition 2 \(Architecture\-Dependent Spectral Structure\)\.

Letϕω:𝒳→ℝd\\phi\_\{\\omega\}:\\mathcal\{X\}\\to\\mathbb\{R\}^\{d\}be a last\-layer feature extractor andΣ\\Sigmaits centered feature covariance with eigenvaluesλ1≥⋯≥λd≥0\\lambda\_\{1\}\\geq\\cdots\\geq\\lambda\_\{d\}\\geq 0and effective rankrαr\_\{\\alpha\}as in Lemma[1](https://arxiv.org/html/2606.02886#A3.Thmtheorem1)\. For spectral operators such as the SFNO, suppose the covariance eigenvalues decay polynomially,λk≤C​k−β\\lambda\_\{k\}\\leq C\\,k^\{\-\\beta\}withβ\>1\\beta\>1\(the*spectral\-decay hypothesis*\); thenrα=O​\(\(1−α\)−1/\(β−1\)\)r\_\{\\alpha\}=O\\\!\\big\(\(1\-\\alpha\)^\{\-1/\(\\beta\-1\)\}\\big\), bounded independent of the spatial resolution and ambient dimensiondd\. For attention\-based architectures, suppose features are input\-dependent convex combinationsϕ​\(x\)=∑i=1mai​\(x\)​vi\\phi\(x\)=\\sum\_\{i=1\}^\{m\}a\_\{i\}\(x\)\\,v\_\{i\}of value vectors\{vi\}\\\{v\_\{i\}\\\}, with value matrixV=\[v1,…,vm\]V=\[v\_\{1\},\\dots,v\_\{m\}\]of full column rankdeffd\_\{\\mathrm\{eff\}\}and attention\-weight covarianceCovx​\[a​\(x\)\]\\mathrm\{Cov\}\_\{x\}\[a\(x\)\]nonsingular on the range ofV⊤V^\{\\top\}\(the*non\-degeneracy hypothesis*, i\.e\. no token or latent collapse\); thenΣ\\Sigmahas rankdeffd\_\{\\mathrm\{eff\}\}, exhibits no polynomial decay, andr0\.99=Θ​\(deff\)r\_\{0\.99\}=\\Theta\(d\_\{\\mathrm\{eff\}\}\), wheredeffd\_\{\\mathrm\{eff\}\}isdheadd\_\{\\mathrm\{head\}\},dlatentd\_\{\\mathrm\{latent\}\}, ordhiddend\_\{\\mathrm\{hidden\}\}for ViT, Perceiver, and GNN\-Transformer architectures respectively\.

###### Proof\.

The spectral\-decay hypothesis is exactly the premise of Lemma[1](https://arxiv.org/html/2606.02886#A3.Thmtheorem1)with exponentβ\>1\\beta\>1, so the bound \([14](https://arxiv.org/html/2606.02886#A3.E14)\) holds and its right\-hand side depends only onβ,C,T\\beta,C,T, not ondd\. For the attention case, centeringϕ\\phiand writinga​\(x\)a\(x\)for the centered attention weights givesϕ~​\(x\)=V​a​\(x\)\\tilde\{\\phi\}\(x\)=Va\(x\), hence

\(16\)Σ=𝔼​\[V​a​\(x\)​a​\(x\)⊤​V⊤\]=V​Covx​\[a​\(x\)\]​V⊤\.\\Sigma=\\mathbb\{E\}\[Va\(x\)a\(x\)^\{\\top\}V^\{\\top\}\]=V\\,\\mathrm\{Cov\}\_\{x\}\[a\(x\)\]\\,V^\{\\top\}\.By Sylvester’s rank inequality,rank​\(Σ\)\\mathrm\{rank\}\(\\Sigma\)equalsdeffd\_\{\\mathrm\{eff\}\}whenVVhas full column rankdeffd\_\{\\mathrm\{eff\}\}andCovx​\[a\]\\mathrm\{Cov\}\_\{x\}\[a\]is nonsingular onrange​\(V⊤\)\\mathrm\{range\}\(V^\{\\top\}\)\. A full\-rank covariance has no vanishing eigenvalues to induce polynomial decay; the energy is spread across alldeffd\_\{\\mathrm\{eff\}\}directions, sorα=Θ​\(deff\)r\_\{\\alpha\}=\\Theta\(d\_\{\\mathrm\{eff\}\}\)\. ∎

Each hypothesis is physically grounded and*checkable from the data*, so the empirics verify the*premises*of Proposition[2](https://arxiv.org/html/2606.02886#A3.Thmtheorem2)while its conclusions follow rigorously from Lemma[1](https://arxiv.org/html/2606.02886#A3.Thmtheorem1)\. The spectral\-decay hypothesis is the spectral signature of the finite Sobolev energy of atmospheric fields under SFNO’s band\-limited spherical convolutions, and is confirmed in Table[7](https://arxiv.org/html/2606.02886#A5.T7)\(FourCastNetV2’s rapidCk/PC\_\{k\}/Pgrowth reflects steep eigenvalue decay\)\. The non\-degeneracy hypothesis is confirmed by the distributed centered spectra of AIFS and Aurora \(λ1≈27\\lambda\_\{1\}\\approx 27–29%29\\%\)\.

This explains why FourCastNetV2 \(SFNO\) requires strict rank truncation \(k≤10k\\leq 10\) to avoid collapse \(Proposition[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)\), while Aurora and AIFS tolerate higher or full\-rank computation\. Pangu\-Weather’s Swin Transformer nominally satisfies the non\-degeneracy hypothesis, but global\-average pooling tod=69d=69dimensions collapses spatial structure and produces a concentrated effective spectrum \(λ1≈99\.6%\\lambda\_\{1\}\\approx 99\.6\\%\), placing it in regime \(i\) despite its attention\-based architecture—an instructive boundary case where the pooling operator, not the backbone, determines the spectral regime\.

### C\.2\.ICA vs SVD: Higher\-Order Statistics for Extreme Events

The empirical superiority of ICA over SVD \(Table[2](https://arxiv.org/html/2606.02886#S5.T2)\) is explained by the*non\-Gaussian structure*of extreme weather events in feature space\. Extreme weather events exhibit heavy\-tailed distributions with positive excess kurtosis, skewness, and multimodal characteristics\. These properties propagate to the learned feature representations when neural networks encode forecast difficulty\.

The key distinction \(Proposition[2](https://arxiv.org/html/2606.02886#S3.Thmtheorem2)\): SVD maximizes variance \(second\-order\), while ICA maximizes statistical independence \(all orders\)\. For Gaussian data, these are equivalent\. For non\-Gaussian extreme weather events, ICA’s higher\-order criterion isolates physical drivers \(vorticity, moisture advection, diabatic heating\) that govern event severity, while SVD’s variance criterion biases toward typical high\-frequency patterns\(Hyvärinen and Oja,[2000](https://arxiv.org/html/2606.02886#bib.bib45); Cardoso,[1999](https://arxiv.org/html/2606.02886#bib.bib47)\)\.

Empirical verification shows joint non\-Gaussianity in the feature datasets:

- •Marginal Gaussianity: 60–75% of individual features pass Shapiro\-Wilk normality tests \(p\>0\.05p\>0\.05\)\.
- •Joint Non\-Gaussianity: 54–77% of feature pairs exhibit significant multivariate third\-order moments\(Mardia,[1970](https://arxiv.org/html/2606.02886#bib.bib48)\), indicating non\-Gaussian dependencies\.
- •Excess Kurtosis: Surface variables \(t2m, msl, winds\) show kurtosis 5–15 during extreme weather events vs\. 3 for Gaussian\.

This joint non\-Gaussianity validates the ICA assumption: features arise from*non\-Gaussian independent sources*\(e\.g\., vorticity, moisture advection, diabatic heating\) mixed by the neural network’s forward pass\. ICA unmixes these sources, isolating extreme event drivers and enabling adaptive uncertainty intervals\.

#### Practical Implication\.

For operational extreme weather forecasting, ICA’s exploitation of higher\-order statistics produces uncertainty estimates that scale with event severity \(Table[5](https://arxiv.org/html/2606.02886#A5.T5), CV=0\.07=0\.07–1\.81\), while SVD’s variance\-only criterion yields more uniform intervals \(CV=0\.01=0\.01–0\.49\)\. This≈5×\{\\approx\}5\{\\times\}improvement in coefficient of variation for AIFS and FourCastNetV2 \(where both ICA and SVD achieve valid coverage\) translates directly to better discrimination of tropical cyclones, atmospheric rivers, and heat waves from routine synoptic conditions, as required for effective early warning systems\.

#### Remark: Non\-Gaussianity vs\. GP Assumption\.

The GP posterior formula \(Eq\. 4\) assumes Gaussian process priorsf∼𝒢​𝒫​\(0,K\)f\\sim\\mathcal\{GP\}\(0,K\)over functions, not Gaussianity of the feature distributionP​\(ϕ\)P\(\\phi\)\. The kernelK​\(x,x′\)=ϕ​\(x\)⊤​ϕ​\(x′\)K\(x,x^\{\\prime\}\)=\\phi\(x\)^\{\\top\}\\phi\(x^\{\\prime\}\)is a valid positive semi\-definite kernel regardless of whether features exhibit non\-Gaussian marginals or heavy tails\. ICA and SVD differ in*how they decompose*the feature matrixΦ\\Phi\(independence vs\. variance maximization\), not in the validity of the GP variance formula itself\. Both methods produce componentscjc\_\{j\}that are plugged into the same GP posterior variance estimator; the non\-Gaussianity affects component selection, not the variance computation\. Furthermore, post\-hoc calibration \(Section[3\.5](https://arxiv.org/html/2606.02886#S3.SS5)\) empirically corrects for any misspecification of the GP prior, ensuring target coverage even when the Gaussian process assumption is violated\. The key advantage of ICA is that by exploiting non\-Gaussian structure during decomposition, it identifies components aligned with extreme event physics, leading to more informative uncertainty estimates after calibration\.

#### Proof of Theorem[1](https://arxiv.org/html/2606.02886#S5.Thmtheorem1)\.

Coverage indicatorsZi=𝟏​\[\|ei\|<z⋅σi\]Z\_\{i\}=\\mathbf\{1\}\[\|e\_\{i\}\|<z\\cdot\\sigma\_\{i\}\]are i\.i\.d\. Bernoulli with𝔼​\[Zi\]=ctrue\\mathbb\{E\}\[Z\_\{i\}\]=c\_\{\\mathrm\{true\}\}\. By the one\-sided Hoeffding inequality:

\(17\)ℙ​\(c^n−ctrue≥t\)≤exp⁡\(−2​n​t2\)\.\\mathbb\{P\}\\left\(\\hat\{c\}\_\{n\}\-c\_\{\\mathrm\{true\}\}\\geq t\\right\)\\leq\\exp\(\-2nt^\{2\}\)\.Settingexp⁡\(−2​n​t2\)=δ\\exp\(\-2nt^\{2\}\)=\\deltaand solving yieldst=ln⁡\(1/δ\)/\(2​n\)t=\\sqrt\{\\ln\(1/\\delta\)/\(2n\)\}, completing the proof\.□\\square

#### Remark on I\.I\.D\. Assumption\.

Weather data exhibits temporal autocorrelation, violating the i\.i\.d\. assumption\. When consecutive samples are positively correlated, the effective sample sizeneff<nn\_\{\\mathrm\{eff\}\}<n, and the bound becomes conservative \(wider\)\. To mitigate this, calibration samples should be temporally spaced \(e\.g\., one sample per week\) or the bound adjusted using techniques for dependent data\(Yu,[1994](https://arxiv.org/html/2606.02886#bib.bib39)\)\. The bound remains valid as an upper bound on coverage deviation even under weak dependence, though it may not be tight\.

## Appendix DProofs of Theoretical Results

### D\.1\.Proof: Last\-Layer Kernel is Positive Semi\-Definite

###### Proof\.

Letϕ:𝒳→ℝd\\phi:\\mathcal\{X\}\\to\\mathbb\{R\}^\{d\}be any feature extractor \(regardless of internal architecture\)\. Define the kernelK​\(x,x′\)=ϕ​\(x\)⊤​ϕ​\(x′\)K\(x,x^\{\\prime\}\)=\\phi\(x\)^\{\\top\}\\phi\(x^\{\\prime\}\)\. For any finite set of points\{x1,…,xn\}⊂𝒳\\\{x\_\{1\},\\ldots,x\_\{n\}\\\}\\subset\\mathcal\{X\}and any vector𝐜∈ℝn\\mathbf\{c\}\\in\\mathbb\{R\}^\{n\}:

\(18\)𝐜⊤​K​𝐜\\displaystyle\\mathbf\{c\}^\{\\top\}K\\mathbf\{c\}=∑i,jci​cj​K​\(xi,xj\)=∑i,jci​cj​ϕ​\(xi\)⊤​ϕ​\(xj\)\\displaystyle=\\sum\_\{i,j\}c\_\{i\}c\_\{j\}K\(x\_\{i\},x\_\{j\}\)=\\sum\_\{i,j\}c\_\{i\}c\_\{j\}\\phi\(x\_\{i\}\)^\{\\top\}\\phi\(x\_\{j\}\)\(19\)=\(∑ici​ϕ​\(xi\)\)⊤​\(∑jcj​ϕ​\(xj\)\)=‖∑ici​ϕ​\(xi\)‖2≥0\\displaystyle=\\left\(\\sum\_\{i\}c\_\{i\}\\phi\(x\_\{i\}\)\\right\)^\{\\top\}\\left\(\\sum\_\{j\}c\_\{j\}\\phi\(x\_\{j\}\)\\right\)=\\left\\\|\\sum\_\{i\}c\_\{i\}\\phi\(x\_\{i\}\)\\right\\\|^\{2\}\\geq 0Since𝐜⊤​K​𝐜≥0\\mathbf\{c\}^\{\\top\}K\\mathbf\{c\}\\geq 0for all𝐜\\mathbf\{c\}, the Gram matrixKKis positive semi\-definite\. This holds for any feature mapϕ\\phi, independent of architecture\. ∎

### D\.2\.Derivation: SVD\-Based Predictive Variance

The GP predictive variance with centered kernel matrixK~=Φ~​Φ~⊤\\tilde\{K\}=\\tilde\{\\Phi\}\\tilde\{\\Phi\}^\{\\top\}whereΦ~∈ℝn×d\\tilde\{\\Phi\}\\in\\mathbb\{R\}^\{n\\times d\}is the centered feature matrix:

\(20\)σ2​\(x∗\)=K~​\(x∗,x∗\)\+σn2−𝐤~∗⊤​\(K~\+σn2​I\)−1​𝐤~∗\\sigma^\{2\}\(x\_\{\*\}\)=\\tilde\{K\}\(x\_\{\*\},x\_\{\*\}\)\+\\sigma\_\{n\}^\{2\}\-\\tilde\{\\mathbf\{k\}\}\_\{\*\}^\{\\top\}\(\\tilde\{K\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\tilde\{\\mathbf\{k\}\}\_\{\*\}where𝐤~∗=Φ~​ϕ~∗\\tilde\{\\mathbf\{k\}\}\_\{\*\}=\\tilde\{\\Phi\}\\tilde\{\\phi\}\_\{\*\}\. Applying the push\-through identityΦ~⊤​\(Φ~​Φ~⊤\+σn2​I\)−1=\(Φ~⊤​Φ~\+σn2​I\)−1​Φ~⊤\\tilde\{\\Phi\}^\{\\top\}\(\\tilde\{\\Phi\}\\tilde\{\\Phi\}^\{\\top\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}=\(\\tilde\{\\Phi\}^\{\\top\}\\tilde\{\\Phi\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\tilde\{\\Phi\}^\{\\top\}, the correction term becomes:

𝐤~∗⊤​\(K~\+σn2​I\)−1​𝐤~∗\\displaystyle\\tilde\{\\mathbf\{k\}\}\_\{\*\}^\{\\top\}\(\\tilde\{K\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\tilde\{\\mathbf\{k\}\}\_\{\*\}=ϕ~∗⊤​Φ~⊤​\(Φ~​Φ~⊤\+σn2​I\)−1​Φ~​ϕ~∗\\displaystyle=\\tilde\{\\phi\}\_\{\*\}^\{\\top\}\\tilde\{\\Phi\}^\{\\top\}\(\\tilde\{\\Phi\}\\tilde\{\\Phi\}^\{\\top\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\tilde\{\\Phi\}\\tilde\{\\phi\}\_\{\*\}\(21\)=ϕ~∗⊤​\(Φ~⊤​Φ~\+σn2​I\)−1​Φ~⊤​Φ~​ϕ~∗\\displaystyle=\\tilde\{\\phi\}\_\{\*\}^\{\\top\}\(\\tilde\{\\Phi\}^\{\\top\}\\tilde\{\\Phi\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\tilde\{\\Phi\}^\{\\top\}\\tilde\{\\Phi\}\\,\\tilde\{\\phi\}\_\{\*\}LetΣ~=Φ~⊤​Φ~∈ℝd×d\\tilde\{\\Sigma\}=\\tilde\{\\Phi\}^\{\\top\}\\tilde\{\\Phi\}\\in\\mathbb\{R\}^\{d\\times d\}with eigendecompositionΣ~=V​Λ​V⊤\\tilde\{\\Sigma\}=V\\Lambda V^\{\\top\}, whereV∈ℝd×dV\\in\\mathbb\{R\}^\{d\\times d\}is unitary andΛ=diag​\(λ1,…,λd\)\\Lambda=\\mathrm\{diag\}\(\\lambda\_\{1\},\\ldots,\\lambda\_\{d\}\)\. SinceVVis square unitary,\(Σ~\+σn2​I\)−1=V​\(Λ\+σn2​I\)−1​V⊤\(\\tilde\{\\Sigma\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}=V\(\\Lambda\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}V^\{\\top\}, giving:

\(22\)ϕ~∗⊤​\(Σ~\+σn2​I\)−1​Σ~​ϕ~∗=ϕ~∗⊤​V​diag​\(λjλj\+σn2\)​V⊤​ϕ~∗=∑j=1kλj​cj2λj\+σn2\\displaystyle\\tilde\{\\phi\}\_\{\*\}^\{\\top\}\(\\tilde\{\\Sigma\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\tilde\{\\Sigma\}\\,\\tilde\{\\phi\}\_\{\*\}=\\tilde\{\\phi\}\_\{\*\}^\{\\top\}V\\,\\mathrm\{diag\}\\\!\\left\(\\tfrac\{\\lambda\_\{j\}\}\{\\lambda\_\{j\}\+\\sigma\_\{n\}^\{2\}\}\\right\)V^\{\\top\}\\tilde\{\\phi\}\_\{\*\}=\\sum\_\{j=1\}^\{k\}\\frac\{\\lambda\_\{j\}\\,c\_\{j\}^\{2\}\}\{\\lambda\_\{j\}\+\\sigma\_\{n\}^\{2\}\}wherecj=ϕ~∗⊤​vjc\_\{j\}=\\tilde\{\\phi\}\_\{\*\}^\{\\top\}v\_\{j\}are the projection coefficients onto the eigenvectorsvjv\_\{j\}ofΣ~\\tilde\{\\Sigma\}, which coincide with the right singular vectors ofΦ~\\tilde\{\\Phi\}\. Substituting back yields the predictive variance formula used throughout:

\(23\)σraw2​\(x∗\)=‖ϕ~∗‖2\+σn2−∑j=1kλj⋅cj2λj\+σn2\\sigma^\{2\}\_\{\\text\{raw\}\}\(x\_\{\*\}\)=\\\|\\tilde\{\\phi\}\_\{\*\}\\\|^\{2\}\+\\sigma\_\{n\}^\{2\}\-\\sum\_\{j=1\}^\{k\}\\frac\{\\lambda\_\{j\}\\cdot c\_\{j\}^\{2\}\}\{\\lambda\_\{j\}\+\\sigma\_\{n\}^\{2\}\}whereσn2\\sigma\_\{n\}^\{2\}is estimated from the eigenvalue tail as described in Section[D\.2](https://arxiv.org/html/2606.02886#A4.SS2)\. Asσn2→0\\sigma\_\{n\}^\{2\}\\to 0, the correction approaches‖Vk⊤​ϕ~∗‖2\\\|V\_\{k\}^\{\\top\}\\tilde\{\\phi\}\_\{\*\}\\\|^\{2\}and the predictive variance reduces to the residual norm‖ϕ~∗‖2−‖Vk⊤​ϕ~∗‖2\\\|\\tilde\{\\phi\}\_\{\*\}\\\|^\{2\}\-\\\|V\_\{k\}^\{\\top\}\\tilde\{\\phi\}\_\{\*\}\\\|^\{2\}\.

#### Remark: Spearman Correlation Invariance under Post\-Hoc Scaling\.

The Spearman rank correlationρs\\rho\_\{s\}between absolute errors\{\|ei\|\}\\\{\|e\_\{i\}\|\\\}and uncertainties\{σi\}\\\{\\sigma\_\{i\}\\\}is invariant under any positive scalar multiple ofσi\\sigma\_\{i\}: sincerank​\(α​σi\)=rank​\(σi\)\\mathrm\{rank\}\(\\alpha\\sigma\_\{i\}\)=\\mathrm\{rank\}\(\\sigma\_\{i\}\)for allα\>0\\alpha\>0, applying the post\-hoc scaleα\\alpha\(Section[3\.5](https://arxiv.org/html/2606.02886#S3.SS5)\) does not alterρs\\rho\_\{s\}\. Consequently,ρs\\rho\_\{s\}measures the discrimination quality of the raw NTK kernel independently of the calibration scale chosen to achieve target coverage\.

## Appendix EAdditional Results

### E\.1\.Discrimination Metrics: CV and Spearmanρs\\rho\_\{s\}

Table[5](https://arxiv.org/html/2606.02886#A5.T5)reports the coefficient of variation \(CV\) of calibrated uncertaintyσ\\sigmaat t\+6h for each decomposition method at optimal rankk∗k^\{\*\}\. ICA produces substantially higher CV than SVD across all models, indicating adaptive intervals that scale with event severity\. Table[6](https://arxiv.org/html/2606.02886#A5.T6)reports Spearmanρs\\rho\_\{s\}between absolute errors and uncertainties for 850 hPa temperature\. AIFS maintains the strongest directional discrimination \(ρs=0\.25\\rho\_\{s\}=0\.25–0\.330\.33\); Pangu\-Weather achievesρs=0\.56\\rho\_\{s\}=0\.56at 6h but degrades to zero at 12–72h, reflecting instability of the single\-component configuration\.

Table 5\.Coefficient of variation \(CV\) of calibrated uncertaintyσ\\sigmaat t\+6h\. Higher CV indicates adaptive intervals that distinguish extreme events from routine conditions\.ICA consistently produces higher CV than SVD\. SVD shows CV<0\.1<0\.1\(nearly uniform intervals\)\.†Pangu\-ICA \(68\.4% coverage\) excluded from valid comparisons\.

Table 6\.Spearman correlation \(ρs\\rho\_\{s\}\) between errors and uncertainties for 850 hPa temperature \(t\_850\) at each model’s optimal rankkk\. Bold values indicateρs\>0\.3\\rho\_\{s\}\>0\.3\.k∗k^\{\*\}: optimal rank byρs\\rho\_\{s\}\. Pangu zero correlations at 12–72h reflect 69\-dim feature instability; high CV does not guarantee highρs\\rho\_\{s\}\.

### E\.2\.Variance Collapse Empirical Validation

Table[7](https://arxiv.org/html/2606.02886#A5.T7)quantifies the correction\-to\-prior ratioRk=Ck/PR\_\{k\}=C\_\{k\}/Pfor FourCastNetV2 at increasing truncation rankkk, empirically validating Proposition[1](https://arxiv.org/html/2606.02886#S3.Thmtheorem1)\. Atk=100k=100, the correction term consumes 92% of the prior variance, leaving only 8% residual—insufficient for meaningful uncertainty discrimination\. The calibration algorithm compensates by scaling uncertainties byα=10,000\\alpha=10\{,\}000, but this uniform scaling cannot recover discriminative power\. Maintainingk≤10k\\leq 10keepsCk/P<15%C\_\{k\}/P<15\\%, preserving discrimination while achieving target coverage\. TheRk<0\.9R\_\{k\}<0\.9threshold is empirically validated here \(collapse occurs atk=100k=100whereRk=0\.92R\_\{k\}=0\.92\)\.

Table 7\.Correction ratioCk/PC\_\{k\}/Pfor FourCastNetV2 at increasing rankkk\.