Paper deep dive
A Unified Uncertainty Representation for Graph Neural Networks via Doubly-Spectral Stochastic Expansion
Fred Xu, Thomas Markovich, Florence Regol, Yizhou Sun
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Reliable deployment of graph neural networks requires calibration, out-of-distribution (OOD) detection, and robustness to distribution shift, yet existing methods address these needs with separate models and objectives. We model uncertain node embeddings as random graph signals: graph Fourier filters capture structural variation, and a scalar orthogonal-polynomial chaos coordinate captures latent stochastic variation. The resulting doubly-spectral stochastic (DSS) expansion supplies task-matched readouts from one representation: the mean coefficient encodes class evidence for the energy-based OOD score, the higher-order coefficients encode structured logit variation, and quadrature averaging over the chaos coordinate defines the single predictive distribution used for prediction and calibration. A capacity theorem shows that, under a full-rank feature assumption, a restricted subfamily matches the chaos coefficients of any Gaussian-latent random graph signal, with exponentially decaying truncation error under a growth condition; the task-level claims are established empirically. DSS-GNN has two deployment modes: standalone, or as a residual branch beside a deterministic encoder (DSS-Hybrid). Standalone DSS-GNN achieves the lowest Brier score among the compared uncertainty-aware baselines on all 14 node classification benchmarks without post-hoc correction; DSS-Hybrid achieves the best AUROC on most node-OOD settings, competitive cross-graph OOD detection, and the strongest shifted accuracy on all 7 GOOD concept-shift benchmarks under standard empirical risk minimization (ERM). Cross-evaluating both modes on all three tasks shows that each remains effective on the other's tasks, with documented exceptions, and yields explicit deployment guidance.
Tags
Links
- Source: https://arxiv.org/abs/2609.35703v1
- Canonical: https://arxiv.org/abs/2609.35703v1
Trouble viewing inline? Open PDF directly →
Full Text
162,691 characters extracted from source content.
Expand or collapse full text
A Unified Uncertainty Representation for Graph Neural Networks via Doubly-Spectral Stochastic Expansion Fred Xu †thanks: Work done during an internship at Block, Inc. Affiliation: Department of Computer Science Affiliation: University of California, Los Angeles Email: fredxu@cs.ucla.edu Thomas Markovich Affiliation: Block, Inc. Email: tmarkovich@squareup.com Florence Regol Affiliation: Mila, Quebec AI Institute Email: florence.regol@mila.quebec Yizhou Sun Affiliation: Department of Computer Science Affiliation: University of California, Los Angeles Email: yzsun@cs.ucla.edu Abstract Reliable deployment of graph neural networks requires calibration, out-of-distribution (OOD) detection, and robustness to distribution shift, yet existing methods address these needs with separate models and objectives. We model uncertain node embeddings as random graph signals: graph Fourier filters capture structural variation, and a scalar orthogonal-polynomial chaos coordinate captures latent stochastic variation. The resulting doubly-spectral stochastic (DSS) expansion supplies task-matched readouts from one representation: the mean coefficient encodes class evidence for the energy-based OOD score, the higher-order coefficients encode structured logit variation, and quadrature averaging over the chaos coordinate defines the single predictive distribution used for prediction and calibration. A capacity theorem shows that, under a full-rank feature assumption, a restricted subfamily matches the chaos coefficients of any Gaussian-latent random graph signal, with exponentially decaying truncation error under a growth condition; the task-level claims are established empirically. DSS-GNN has two deployment modes: standalone, or as a residual branch beside a deterministic encoder (DSS-Hybrid). Standalone DSS-GNN achieves the lowest Brier score among the compared uncertainty-aware baselines on all 14 node classification benchmarks without post-hoc correction; DSS-Hybrid achieves the best AUROC on most node-OOD settings, competitive cross-graph OOD detection, and the strongest shifted accuracy on all 7 GOOD concept-shift benchmarks under standard empirical risk minimization (ERM). Cross-evaluating both modes on all three tasks shows that each remains effective on the other’s tasks, with documented exceptions, and yields explicit deployment guidance. 1 Introduction Graph neural networks (GNNs) perform well on graph-structured data, yet reliable deployment in drug discovery [31], fraud detection [7], or traffic forecasting [18] requires calibrated probabilities, out-of-distribution (OOD) detection, and robustness to distribution shift. On graphs the three needs are coupled, because uncertainty depends on neighborhood structure as well as node features, and they are instances of one question: when, and how confidently, can the model predict? Calibration measures the agreement between confidence and accuracy in distribution (ID), OOD detection requires low confidence outside the training support, and shift robustness requires that accuracy decrease slowly as the input distribution changes. A representation that explicitly encodes predictive uncertainty addresses all three, with calibrated probabilities, OOD evidence, and a regularized residual under shift obtained from one signal. Most existing methods target one requirement: sampling-based uncertainty methods add computational overhead [10, 17, 19]; evidential models such as Graph Posterior Network (GPN) [30] provide OOD evidence but not shift-robust classification, and intrinsic uncertainty models such as G-Δ [33] target calibration under shift; energy-based detectors [36, 9] provide separation scores rather than calibrated probabilities; and robust-training methods [27, 1, 38] do not model uncertainty explicitly. Spectral GNNs learn filters along the graph-frequency axis [6]; we add a second spectral axis for stochastic variation. DSS-GNN represents each node embedding as a polynomial chaos expansion in a scalar Gaussian latent [37] and propagates the joint graph-frequency/chaos-order coefficients through a learnable doubly-spectral layer (Figure 4, Appendix H). One deterministic forward pass gives task-matched readouts: the quadrature-averaged predictive distribution, whose argmax is the point predictor, is used for prediction and calibration, with the higher-order coefficients encoding structured logit variation; the zeroth-order mean logit gives an energy score for OOD detection that retains the common-logit shifts that softmax cancels; and, paired with a deterministic encoder (DSS-Hybrid, e.g., GCN [15]), the DSS branch acts as a regularized low-/high-pass residual under distribution shift. DSS-GNN is thus one uncertainty representation with two deployment modes; Section 4.5 evaluates both on all three tasks and states which to prefer where. Our contributions are: 1. A doubly-spectral stochastic representation for GNNs with a capacity theorem: under a full-rank feature assumption, a restricted DSS subfamily matches the chaos coefficients of any Gaussian-latent random graph signal, with exponentially decaying truncation error for targets with an entire extension of Gaussian growth (Theorem 2). The theorem (Section 3.3) establishes representational capacity, consistent with the empirical plateau at P=1P=1–22, but does not by itself guarantee calibration, OOD detection, or robustness, which are supported empirically. 2. A learnable doubly-spectral layer whose quadrature-averaged loss differs from the mean-logit loss by a curvature-weighted correction over the chaos covariance (Proposition 2), relating the chaos coefficients to calibration, bounding the gap between the quadrature predictive and the mean-logit softmax (Corollary 1), and identifying the mean-logit energy as the readout that retains common-logit shifts. 3. Empirically, the lowest Brier score among compared uncertainty-aware baselines on all 14 calibration benchmarks (standalone), the best AUROC on most node-OOD settings with competitive cross-graph OOD, and the strongest shifted accuracy on all 7 concept-shift benchmarks under ERM (hybrid), plus a cross-evaluation of both modes on all three tasks, component-wise ablations (Table 21), and explicit deployment guidance. 2 Background Spectral GNNs and uncertainty on graphs. A graph signal on =(V,E)G=(V,E) with N nodes is a vector x∈ℝNx ^N in the Hilbert space ℓ2(V) ^2(V); the Laplacian LGL_G is self-adjoint and its eigenvectors form the graph Fourier basis [29, 24]. Spectral GNNs define graph convolutions as polynomial filters of LGL_G; ChebNet [6] uses g(LG)=∑k=0KckTk(L~G)g(L_G)= _k=0^Kc_kT_k( L_G) with learned coefficients ckc_k, Chebyshev polynomials TkT_k, and the rescaled Laplacian L~G=2LG/λmax−I L_G=2L_G/ _ -I (λmax _ the largest eigenvalue), which keeps filtering localized and allows low-pass or high-pass responses under homophily or heterophily [2]. DSS-GNN adds a second spectral axis, chaos order, orthogonal to graph frequency, and embeds uncertainty directly into spectral filtering. Polynomial chaos and the doubly-spectral expansion. To organize stochastic variation spectrally, as the graph Fourier basis organizes structural variation, we use polynomial chaos expansion (PCE) [35, 4, 37], classical in uncertainty quantification [23] but, to our knowledge, not previously combined with graph spectral filtering inside a GNN (Appendix A.1 gives a self-contained introduction). PCE represents a square-integrable random variable Y(ω)Y(ω) as a series in the orthonormal Hermite polynomials Ψn(ω)n≥0\ _n(ω)\_n≥ 0 of a Gaussian latent ω∼(0,1)ω (0,1); truncating at order P gives Y(ω)≈∑n=0PYnΨn(ω)Y(ω)≈ _n=0^PY_n _n(ω), where Y0Y_0 is the mean and ∑n=1PYn2 _n=1^PY_n^2 (the chaos energy) the variance. Combined with the Laplacian eigenbasis uj\u_j\ on ℓ2(V) ^2(V), the Hermite basis on L2(ℝ,γ)L^2(R,γ) (γ the standard Gaussian measure) gives a tensor-product representation: a random graph signal determined by ω is in ℓ2(V)⊗L2(ℝ,γ) ^2(V) L^2(R,γ) and admits the doubly-spectral stochastic expansion Y(ω)=∑n=0∞∑j=1NujY^n,j⊤Ψn(ω),Y(ω)= _n=0^∞ _j=1^Nu_j\, Y_n,j \, _n(ω), (1) where Y^n,j∈ℝd Y_n,j ^d is the joint coefficient at graph frequency λj _j (the j-th Laplacian eigenvalue) and chaos order n, and d is the signal dimension; uj⊗Ψn\u_j _n\ is a complete orthonormal basis of this space (Remark 2), ω is scalar, and each chaos coefficient is a d-dimensional graph signal. DSS-GNN learns a truncation and a propagation rule on this tensor-product space. 3 Method DSS-GNN instantiates the doubly-spectral expansion (1) as a trainable network by learning a coupled propagation rule for the joint coefficients Y^n,j Y_n,j through parameterized graph filters and non-intrusive chaos projection (Figure 4, Appendix H). Intuition. A random graph signal assigns to each node a random variable with finite second moment. We model all randomness as driven by a single shared standard Gaussian factor ω: every node embedding is a deterministic function of ω, expanded in the Hermite polynomials Ψn(ω) _n(ω), which play along the stochastic axis exactly the role that the graph Fourier modes play along the structural axis. A DSS layer therefore stores, for every node, P+1P+1 coefficient vectors (one per chaos order) instead of a single hidden vector; a deterministic embedding is the special case in which only the order-0 coefficient is nonzero. Each layer decomposes the embeddings onto this basis, applies learned filters jointly per graph frequency and per chaos order (dual low/high-pass Chebyshev filters along the graph axis, per-order gates along the chaos axis), and reassembles. The core Gaussian is never materialized: because pointwise nonlinearities mix chaos orders, the layer evaluates the embedding field at S deterministic Gauss–Hermite nodes, applies the activation there, and projects back onto the Hermite basis. With S≥P+1S≥ P+1 the polynomial filter and gate terms are integrated exactly at the gate orders used in practice (Pg≤1P_g≤ 1, Appendix A), and the nonlinearity is handled by a controlled non-intrusive projection (Remark 3), so the network is deterministic end to end: no Monte-Carlo sampling and no sampling variance. The expansion itself requires only the shared Gaussian factor and finite second moments; the exponential truncation rate of Theorem 2 additionally requires targets with an entire (analytic) Gaussian extension. Readouts then select task-appropriate functionals of the same expansion: the quadrature-averaged predictive (8) for prediction and calibration, and the mean-logit energy (9) for OOD scores. 3.1 DSS layer Let =(V,E)G=(V,E) have N=|V|N=|V| nodes, features X∈ℝN×dinX ^N× d_ in, and Laplacian LGL_G. DSS-GNN stacks L layers on a truncated expansion at chaos order P; the input is lifted by Hn(0)=XWin(n)H_n^(0)=XW_in^(n) with Win(n)∈ℝdin×d0W_in^(n) ^d_ in× d_0, and at layer ℓ∈0,…,L−1 ∈\0,…,L-1\ the embeddings are H(ℓ)(ω)=∑n=0PHn(ℓ)Ψn(ω),Hn(ℓ)∈ℝN×dℓ,H^( )(ω)= _n=0^PH_n^( ) _n(ω), H_n^( ) ^N× d_ , (2) with dℓd_ the hidden width, so each node maintains P+1P+1 coefficient-valued hidden states. The layer updates them along both spectral axes. Along the graph-frequency axis, dual Chebyshev filters propagate each coefficient through localized low-pass and high-pass branches, Glp(LG)=∑k=0Klpck(lp)Tk(L~G),Ghp(LG)=∑k=0Khpck(hp)Tk(L~G),G_lp(L_G)= _k=0^K_lpc_k^(lp)T_k( L_G), G_hp(L_G)= _k=0^K_hpc_k^(hp)T_k( L_G), (3) with filter degrees Klp,KhpK_lp,K_hp and learned coefficients ck(lp),ck(hp)∈ℝc_k^(lp),c_k^(hp) . Along the chaos-order axis, per-order chaos gates αn(ω),βn(ω) _n(ω), _n(ω) couple the two branches; they are deterministic scalars or low-order chaos expansions of gate order PgP_g with learnable an,r,bn,r∈ℝa_n,r,b_n,r : αn(ω)=∑r=0Pgan,rΨr(ω),βn(ω)=∑r=0Pgbn,rΨr(ω). _n(ω)= _r=0^P_ga_n,r _r(ω), _n(ω)= _r=0^P_gb_n,r _r(ω). (4) Because pointwise nonlinearities mix chaos orders, the update has no closed form. A non-intrusive projection evaluates the field at S Gauss–Hermite nodes (ωs,μs)s=1S\( _s, _s)\_s=1^S (nodes ωs∈ℝ _s , weights μs>0 _s>0), applies the activation σ, and projects back: H~(ℓ+1)(ωs) H^( +1)( _s) =(∑n=0PΨn(ωs)(αn(ωs)Glp(LG)Hn(ℓ)+βn(ωs)Ghp(LG)Hn(ℓ)))Wℓ. = ( _n=0^P _n( _s) ( _n( _s)G_lp(L_G)H_n^( )+ _n( _s)G_hp(L_G)H_n^( ) ) )W_ . (5) Hm(ℓ+1)=∑s=1Sμsσ(H~(ℓ+1)(ωs))Ψm(ωs),m=0,…,P,H_m^( +1)= _s=1^S _s\,σ\! ( H^( +1)( _s) ) _m( _s), m=0,…,P, (6) with Wℓ∈ℝdℓ×dℓ+1W_ ^d_ × d_ +1 trainable (quadrature details in Appendix A.1; notation in Table 30, Appendix H). 3.2 Uncertainty readout After L layers a linear readout Wout∈ℝdL×CW_out ^d_L× C (C classes, labels yi∈1,…,Cy_i∈\1,…,C\) produces per-order logit coefficients Zi,n:=(Hn(L)Wout)i∈ℝCZ_i,n:=(H_n^(L)W_out)_i ^C for node i. The mean logit Zi,0Z_i,0 encodes the zeroth-order class evidence used by the energy score; the logit covariance and chaos energy are Σzi=∑n=1PZi,nZi,n⊤∈ℝC×C,ℰi=tr(Σzi)=∑n=1P‖Zi,n‖2. _z_i= _n=1^PZ_i,nZ_i,n ^C× C, _i=tr( _z_i)= _n=1^P\|Z_i,n\|^2. (7) The quadrature logits zi(s)=∑n=0PZi,nΨn(ωs)z_i^(s)= _n=0^PZ_i,n _n( _s) are deterministic evaluations of the learned expansion (no Gaussian noise is injected in training or inference) and define the readouts. The quadrature-averaged predictive p¯i=∑s=1Sμssoftmax(zi(s))∈ℝC p_i= _s=1^S _s\,softmax(z_i^(s)) ^C (8) is the model’s predictive distribution and the point prediction is its argmax, y^i=argmaxcp¯i,c y_i= _c p_i,c, so prediction and calibration are two functionals of one distribution. Its deviation from the mean-logit softmax is bounded by the chaos covariance, ‖p¯i−softmax(Zi,0)‖∞≤B2tr(Σzi)\| p_i-softmax(Z_i,0)\|_∞≤ B2tr( _z_i) for a universal constant B (Corollary 1), which training keeps small (Section 3.5); Section 4.2 measures the agreement. For OOD detection, the energy score [20] senergy(i)=−log∑c=1Cexp(Zi,0,c),s_energy(i)=- _c=1^C (Z_i,0,c), (9) with Zi,0,cZ_i,0,c the c-th mean logit, captures confidence in the mean prediction; predictive entropy and mutual information of the quadrature softmaxes are auxiliary summaries of the same readout. The OOD experiments use the energy score with graph propagation, following GNNSafe [36]. 3.3 Theoretical properties Four results connect the layer to the expansion (1); proofs are in Appendix A. Proposition 1 (Generalization). With P=0P=0 and Pg=0P_g=0, (2) has a single deterministic coefficient H0(ℓ)H_0^( ) with scalar gates, and the DSS layer reduces to a standard polynomial spectral GNN layer, recovering ChebNet, GCN, and related models as special cases. Theorem 1 (Doubly-spectral operator). With the pointwise nonlinearity omitted, gates of chaos degree at most PgP_g, states of degree at most P, and a Gauss–Hermite rule exact for all triple products ΨrΨnΨm _r _n _m (r≤Pgr≤ P_g; n,m≤Pn,m≤ P; S≥⌈(Pg+2P+1)/2⌉S≥ (P_g+2P+1)/2 suffices), the linearized DSS layer is a linear operator on the joint graph-frequency/chaos-order domain of (1), acting frequency by frequency: the Chebyshev filters set the graph-frequency response and the chaos gates the chaos-order mixing (Appendix A.4; the omitted nonlinearity’s projection error is decomposed in Proposition 3). Theorem 2 (Representational capacity and truncation under full-rank features). Let G have N nodes and features X∈ℝN×dinX ^N× d_ in with rank(X)=Nrank(X)=N (which requires din≥Nd_ in≥ N). With ReLU activation, a restricted subfamily of DSS-GNNs (one layer, lifted width 2d(P+1)2d(P+1), quadrature size S=P+1S=P+1, constant gates, identity low-pass filter) has the following properties. For any target Y∈L2(ℝ,γ,ℝN×d)Y∈ L^2(R,γ;R^N× d) and any chaos order P, some output Y~P Y_P matches the chaos coefficients of Y through order P; letting P and width vary, the resulting outputs are dense in L2(ℝ,γ,ℝN×d)L^2(R,γ;R^N× d). If Y in addition admits an entire Gaussian extension with controlled growth (Appendix A.6), the truncation error decays exponentially: ‖Y−Y~P‖L2(ℝ,γ,ℝN×d)≤CtrβP\|Y- Y_P\|_L^2(R,γ;\,R^N× d)≤ C_ tr\,β^P (10) for some Ctr>0C_ tr>0 and β∈(0,1)β∈(0,1) depending on Y. Theorems 1 and 2 are capacity statements: they establish the operator structure and the representational capacity of the layer and identify no inherent expressivity obstruction, but they do not by themselves guarantee calibration, OOD detection, or robustness, which are empirical claims supported task by task in Section 4. The full-rank case requires din≥Nd_ in≥ N; for rank(X)<Nrank(X)<N, the common case (Cora has 1,433 features for 2,708 nodes), Remark 5 shows that the same construction still matches every truncated target whose chaos-coefficient columns lie in range(X)range(X), under one positivity condition satisfiable by appending a constant feature column. Proposition 2 (Stochastic and mean-logit readouts). Assume the number of Gauss–Hermite quadrature nodes is at least one more than the chaos order, and that the per-node loss ℓ(⋅,yi) (·,y_i) satisfies the regularity conditions in Lemma 3. Then: 1. (Calibration channel.) The per-node quadrature-averaged loss ℓ¯i:=∑s=1Sμsℓ(zi(s),yi) _i:= _s=1^S _s\, (z_i^(s),y_i) admits the Taylor expansion ℓ¯i=ℓ(Zi,0,yi)+12tr(ℋiΣzi)+O(maxs‖zi(s)−Zi,0‖3), _i= (Z_i,0,y_i)+ 12tr(H_i\, _z_i)+O\! ( _s\|z_i^(s)-Z_i,0\|^3 ), where ℋi=∇z2ℓ(Zi,0,yi)H_i= _z^2 (Z_i,0,y_i) is the loss Hessian and Σzi _z_i is the chaos covariance (7). 2. (OOD channel.) Shifting every class logit by the same constant leaves the softmax output, and hence any softmax-derived summary such as predictive entropy or mutual information, unchanged; by contrast, the energy score (9) decreases by exactly that constant. Calibration: the curvature-weighted term 12tr(ℋiΣzi) 12tr(H_i _z_i) is the leading-order difference between the quadrature-averaged and the mean-logit objectives; it weights the chaos covariance by the loss curvature and, for cross-entropy, is at most 14ℰi 14E_i, the quantity that (11) regularizes; no post-hoc temperature scaling is used. The same expansion applied to the softmax readout bounds the gap between p¯i p_i and softmax(Zi,0)softmax(Z_i,0) by the chaos covariance (Corollary 1), which the objective keeps small; the calibration gains of P>0P>0 therefore arise through training rather than a readout-time correction, and the two trained readouts agree closely (Section 4.2). OOD detection: the energy score is linear in a common logit shift, while softmax-derived summaries cancel it, so an OOD input that uniformly lowers the class-evidence logits changes the energy but not the entropy; Section 4.3 confirms that detection is insensitive to the chaos order, as this readout predicts. 3.4 Hybrid deployment Distribution shift weakens the common-evidence signal and changes which graph frequencies are predictive; like OOD detection it requires a strong deterministic predictor on seen structure paired with a regularized residual that accounts for cross-distribution variability. DSS-Hybrid adds DSS-GNN as a residual branch beside a deterministic base encoder with one shared class readout, zi,0hyb=zibase+γresZi,0z_i,0^hyb=z_i^base+ _resZ_i,0, with zibase∈ℝCz_i^base ^C the base logit and γres∈ℝ _res a learned scale (lowercase for post-combination logits, uppercase Zi,nZ_i,n for the DSS coefficients); the quadrature logits combine the same way, zihyb,(s)=zibase+γres∑n=0PZi,nΨn(ωs)z_i^hyb,(s)=z_i^base+ _res _n=0^PZ_i,n _n( _s), so the readouts of Section 3.2 transfer directly. The chaos-energy regularizer keeps the residual from fitting in-distribution noise, so its dual filters retain capacity for shifts in graph-frequency response and its mean logit remains a common-evidence channel for the energy score. The base encoder, a 2-layer GCN for fair comparison with baselines, is warmed up before the DSS branch is activated. This design is an inductive bias rather than a distribution-free guarantee. 3.5 Training and downstream tasks Both modes are trained with the chaos-energy-regularized mean-logit loss ℒbase=1|tr|∑i∈trℓ(Zi,0,yi)+λreg1|tr|∑i∈trℰi,L_base= 1|V_tr| _i _tr (Z_i,0,\,y_i)\;+\; _reg 1|V_tr| _i _trE_i, (11) with trV_tr the training nodes, ℓ the per-node cross-entropy, and λreg≥0 _reg≥ 0; for DSS-Hybrid, Zi,0Z_i,0 is replaced by zi,0hybz_i,0^hyb and the regularizer applies to the DSS coefficients only. Only the mean logit is supervised, but the projection update (6) couples Zi,0Z_i,0 to all chaos orders, so the higher-order coefficients are trained implicitly while the regularizer constrains their variation. By Proposition 2, ℒbaseL_base agrees to first order with the quadrature-averaged loss ℒquadL_quad (35) of Appendix D; we use ℒbaseL_base throughout. The trace that Corollary 1 bounds is the regularized chaos energy ℰiE_i; where the selected λreg _reg is small or zero, the measured readout gap remains small (Table 27). 4 Experiments 4.1 Setup For calibration, standalone DSS-GNN is compared with TFE-GNN [8] and G-Δ [33] on 14 node classification benchmarks [15, 28, 25, 26], reporting accuracy and the Brier score [3] over 10 splits without post-hoc calibration. For OOD detection, DSS-Hybrid is evaluated on the GNNSafe benchmark [36] (nine node-OOD settings on Cora, Amazon-Photo, and Coauthor-CS, plus cross-graph Twitch and Arxiv; 3 seeds) against GNNSafe, GNNSafe++, and GPN [30] at the values reported by Wu et al. [36], and against Graph-EBM [9] (scored with its own diffusion), MC-dropout (M=20M=20) [10] and Deep Ensembles (M=5M=5) [17] on the same GCN backbone with the same score propagation. For distribution shift, DSS-Hybrid is trained with empirical risk minimization (ERM) on 7 concept-shift settings of GOOD [12] and compared with the 12 GOOD/TAR baselines [1, 27, 38] and G-Δ . Point predictions are argmaxp¯i p_i (Section 3.2); Table 1 reports the mean-logit readout, which Section 4.2 shows to agree with it. All experiments use ℒbaseL_base (11); in the OOD experiments DSS-Hybrid, like GNNSafe++, adds the GNNSafe++ energy-margin regularizer at per-dataset margins. Baselines, datasets, and configurations are in Appendices B–D; code is available at https://github.com/heraclixus/DSSGNN. 4.2 Node classification and calibration Table 1 shows that DSS-GNN achieves the lowest Brier score among the compared uncertainty-aware baselines on all 14 datasets, homophilous and heterophilous, without post-hoc calibration, and the best accuracy on 13. The Brier gains over the deterministic P=0P=0 arm are largest on Citeseer, Texas, and Minesweeper (Table 14), consistent with the curvature-weighted correction of Proposition 2. This is not an accuracy effect: standard GNNs match or exceed DSS-GNN in accuracy on Cora and Citeseer while remaining less calibrated (Appendix F), and at P=0P=0 DSS-GNN is a competitive dual-filter spectral GNN (Proposition 1), so the gains are attributable to the chaos expansion (Table 19). Figure 1 visualizes the correction ci=12tr(ℋiΣzi)c_i= 12tr(H_i _z_i) on Minesweeper: it is largest where stochastic variation aligns with high-curvature directions and concentrates on the hardest nodes. Prediction and calibration use one predictive distribution. Under the full protocol with both readouts evaluated at the same checkpoints, argmaxp¯i p_i and argmaxZi,0 Z_i,0 disagree on at most 0.21%0.21\% of test nodes (exactly zero on 10 of 14 datasets), with accuracy within 0.040.04 points and Brier within 0.0030.003 (Table 27, Appendix G); in hybrid mode the disagreement is 0.00000.0000 on all 14 datasets. The agreement is structural: for S≥P+1S≥ P+1 Corollary 1 bounds the readout gap by tr(Σzi)tr( _z_i), which the objective keeps small (Section 3.5), so the reported calibration is that of the distribution the predictions are drawn from. Chaos order. Sweeping P∈0,…,6,8P∈\0,…,6,8\ (Appendix E, Figure 3) improves Brier on 12 of 14 datasets at P>0P>0, and at least one of P=1P=1 or P=2P=2 is within 0.020.02 Brier of the best swept order on every dataset, consistent with the exponential truncation bound of Theorem 2. P is selected per dataset on validation Brier (Table 12); the larger selected orders are near-ties rather than requirements, and we recommend P=1–2P=1--2 as the default. Within the shared-factor family P=1P=1 is exactly a Gaussian embedding; the higher orders add non-Gaussian structure (Citeseer Brier 0.3460.346 to 0.3080.308 at P=2P=2; Tables 13–14). Figure 1: (a) Second-order correction cic_i (Proposition 2) on Minesweeper. Each node is positioned by loss curvature (tr(ℋi)tr(H_i)) and chaos energy (ℰiE_i); color shows correction magnitude. (b) Nodes binned by correction quartile: the correction concentrates on hard nodes (low accuracy, high Brier). Table 1: Node classification: test accuracy (%, ↑ ) and Brier score (↓ ). Mean ± standard deviation over 10 runs. Best result per dataset in bold. DSS-GNN uses the per-dataset configuration of Table 12 and the mean-logit readout. Accuracy (%) Brier Score Dataset DSS-GNN TFE-GNN G-Δ DSS-GNN TFE-GNN G-Δ Cora 86.63 ± 1.26 71.25 ± 1.90 84.35 ± 1.97 0.207 ± 0.019 0.449 ± 0.024 0.226 ± 0.023 Citeseer 80.20 ± 1.28 68.87 ± 1.87 70.85 ± 2.22 0.308 ± 0.010 0.463 ± 0.020 0.433 ± 0.026 PubMed 89.72 ± 0.31 82.84 ± 0.76 88.18 ± 0.66 0.156 ± 0.005 0.265 ± 0.011 0.179 ± 0.008 Texas 91.31 ± 3.44 84.92 ± 4.26 9.02 ± 1.83 0.240 ± 0.111 0.255 ± 0.083 0.779 ± 0.029 Cornell 85.11 ± 5.12 81.28 ± 5.45 21.06 ± 3.49 0.231 ± 0.079 0.290 ± 0.069 0.793 ± 0.023 Wisconsin 93.25 ± 3.39 76.50 ± 23.95 52.25 ± 17.80 0.104 ± 0.047 0.292 ± 0.186 0.571 ± 0.109 Chameleon 75.36 ± 1.19 73.41 ± 1.53 37.31 ± 2.36 0.339 ± 0.019 0.365 ± 0.016 0.765 ± 0.011 Squirrel 68.06 ± 2.21 71.56 ± 1.42 25.88 ± 1.91 0.394 ± 0.027 0.395 ± 0.017 0.795 ± 0.006 CS 96.17 ± 0.23 95.05 ± 0.24 94.37 ± 0.23 0.062 ± 0.004 0.082 ± 0.005 0.085 ± 0.004 Roman-Emp. 78.84 ± 0.61 70.35 ± 0.60 51.24 ± 1.07 0.300 ± 0.007 0.414 ± 0.006 0.641 ± 0.008 Amz-Rat. 49.72 ± 0.64 48.79 ± 0.43 44.43 ± 0.42 0.634 ± 0.005 0.637 ± 0.004 0.674 ± 0.003 Minesweeper 87.66 ± 1.10 82.42 ± 0.33 80.22 ± 0.22 0.169 ± 0.014 0.248 ± 0.003 0.288 ± 0.003 Tolokers 79.88 ± 0.49 78.18 ± 0.10 79.75 ± 0.55 0.268 ± 0.005 0.382 ± 0.004 0.289 ± 0.007 Questions 97.17 ± 0.04 97.07 ± 0.09 97.06 ± 0.05 0.053 ± 0.001 0.061 ± 0.001 0.054 ± 0.001 4.3 OOD detection Tables 2–3 compare DSS-Hybrid with GNNSafe++ [36], the strongest OOD baseline: it has the best AUROC on 6 of 9 node-OOD cells (Cora 94.32/97.60/94.1194.32/97.60/94.11 versus 90.62/95.56/92.7590.62/95.56/92.75 for structure/feature/label; Amazon-Photo feature/label 99.66/97.5299.66/97.52; Coauthor-CS label 98.0798.07) and is within 0.130.13 points of the best method on the three remaining cells, where both it and GNNSafe++ exceed 99.699.6. On Twitch it improves over GNNSafe++ on every metric (AUROC 95.7595.75 vs 95.3695.36, FPR95 22.4022.40 vs 33.5733.57); on Arxiv it gives the best FPR95 and ID accuracy while remaining within 1.411.41 AUROC points of GNNSafe++ (73.3673.36 vs 74.7774.77). DSS-Hybrid follows the GNNSafe++ protocol (same GCN backbone configuration, energy-margin objective with its per-dataset margins, and score propagation; Appendix D), so the two differ in the DSS residual branch. Against MC-dropout and Deep Ensembles, which receive the same propagation, DSS-Hybrid has the higher AUROC on all 11 settings in one forward pass versus their 20 passes and 5 models; the published GPN values are below GNNSafe++ throughout. Table 2: OOD detection (AUROC, %) on the nine node-OOD settings and the two cross-graph settings (AUPR, FPR95, and ID accuracy for the latter are in Table 3). GNNSafe/GNNSafe++ and GPN [30] values as reported by Wu et al. [36] on this benchmark (their Tables 1–2; their appendix lists 81.77/93.2481.77/93.24 for GPN on Coauthor-CS feature/label; GPN out-of-memory on Arxiv); Graph-EBM [9], MC-dropout (M=20M=20), and Deep Ensembles (M=5M=5) run in our pipeline on the same GCN backbone (3 seeds; the latter two with the same score propagation). Bold: best. Model Cora Amazon-Photo Coauthor-CS Cross-graph Structure Feature Label Structure Feature Label Structure Feature Label Twitch Arxiv GNNSafe 87.52 93.44 92.80 99.58 98.55 97.35 99.60 99.64 97.23 66.82 71.06 GNNSafe++ 90.62 95.56 92.75 99.82 99.64 97.51 99.99 99.97 97.89 95.36 74.77 Graph-EBM 61.14 72.42 92.69 75.22 86.49 97.24 72.75 89.28 97.91 44.43 52.80 MC-dropout 87.50 93.12 93.03 98.67 98.50 96.80 99.51 99.53 97.25 68.58 66.34 Deep Ensemble 87.87 93.63 93.75 98.58 98.43 97.35 98.18 98.45 94.97 72.25 67.72 GPN 77.47 85.88 90.34 97.17 87.91 92.72 34.67 72.56 83.65 51.73 OOM DSS-Hybrid 94.32 97.60 94.11 99.69 99.66 97.52 99.96 99.95 98.07 95.75 73.36 Table 3: Cross-graph OOD (AUROC, AUPR: ↑ ; FPR95, false positive rate at 95% true positive rate: ↓ ; ID: in-distribution; bold: best). GNNSafe/GNNSafe++ and GPN values as reported by Wu et al. [36] (OOM: out of memory on a 24 GB GPU); Graph-EBM [9], MC-dropout, and Deep Ensembles run in our pipeline (3 seeds). Twitch Arxiv Model AUROC AUPR FPR95 ID Acc. AUROC AUPR FPR95 ID Acc. GNNSafe 66.82 70.97 76.24 70.40 71.06 80.44 87.01 53.39 GNNSafe++ 95.36 97.12 33.57 70.18 74.77 83.21 77.43 53.50 Graph-EBM 44.43 57.91 94.84 63.04 52.80 65.88 97.40 53.45 MC-dropout 68.58 81.25 94.84 64.10 66.34 74.88 89.69 53.56 Deep Ensemble 72.25 83.19 94.29 66.94 67.72 75.98 87.32 47.58 GPN 51.73 66.36 95.51 68.09 OOM OOM OOM OOM DSS-Hybrid 95.75 97.51 22.40 70.52 73.36 79.11 76.67 53.99 Detection comes from the energy readout, with propagation, acting on the DSS mean logit (9), which Proposition 2 formalizes as the readout that retains common-evidence shifts. Under this identical readout the DSS representation separates ID from OOD more widely than the GCN encoder alone (the Cora columns of Table 2), and a controlled ablation that varies only the chaos order P∈0,…,3P∈\0,…,3\ at one fixed configuration leaves energy AUROC flat (Cora-structure 91.0/91.2/90.9/90.891.0/91.2/90.9/90.8; the level differs from Table 2 as the ablation is a separate set of runs; Appendix E), as the design predicts since the energy score reads the mean logit alone. Score propagation contributes up to +18.2+18.2 AUROC (Twitch) and +14.5+14.5 (Cora-structure) at identical checkpoints (Appendix E). Figure 2 shows the resulting wider energy-density gap between ID and OOD modes relative to the energy-based baselines. Figure 2: Energy density for ID (blue) vs. OOD (red). Top: Cora/structure. Bottom: Twitch (cross-graph). Lower energy indicates higher confidence. DSS-Hybrid separates the ID and OOD modes on both datasets. 4.4 Distribution-shifted classification The GOOD benchmark [12] tests generalization under distribution shift. Table 4 compares DSS-Hybrid, trained with standard ERM, against the GOOD/TAR baselines and G-Δ on seven concept-shift settings (covariate-shift splits in Appendix G.7): it has the best mean shifted accuracy in every case, by margins from 0.080.08 (GOOD-Cora/word) to 9.949.94 points (GOOD-WebKB/university). This is consistent with the inductive bias of Section 3.4: the regularized low-/high-pass residual can model the shift in which graph frequencies are predictive while the base predictor is preserved. Section 4.5 isolates the residual pathway as a unit (hybrid versus its identically trained GCN base), and Appendix E (Table 15) isolates the chaos order under shift: at the validation-selected standalone configurations, P=1P=1 improves shifted accuracy over P=0P=0 on GOOD-Cora/word (+1.7+1.7), GOOD-Cora/degree (+1.2+1.2), and GOOD-Arxiv/time (+0.6+0.6), so on those settings the uncertainty representation, not only the dual-filter structure, contributes. Table 4: Test accuracy (%) on the concept-shift GOOD node benchmark for DSS-Hybrid (standard ERM training) and 13 baselines; ERM through TAR as reported by Zheng et al. [38], G-Δ run in our pipeline. Best in bold, second-best underlined; OOM: out of memory. Model GOOD-CBAS GOOD-WebKB GOOD-Twitch GOOD-Cora GOOD-Cora GOOD-Arxiv GOOD-Arxiv color university language word degree time degree ERM 82.43 27.16 51.59 64.03 60.30 65.64 54.81 IRM 82.00 26.06 49.78 63.93 60.26 65.54 56.72 VREx 82.86 26.61 55.75 64.03 60.53 65.92 56.68 Coral 81.57 28.07 51.80 64.04 60.30 65.79 55.14 DANN 83.57 29.36 51.67 63.96 60.23 65.67 55.34 SRGNN 82.14 26.42 51.58 63.96 60.27 65.64 55.08 EERM 65.71 29.91 OOM 63.42 60.21 OOM OOM FLOOD 84.29 28.62 54.22 64.01 60.31 65.66 58.59 CIT 83.71 28.99 OOM 63.77 60.05 OOM OOM KL-DRO 81.14 29.54 51.87 64.03 60.52 65.51 54.70 GroupDRO 82.71 29.17 52.24 64.10 60.43 65.93 56.24 TAR 87.29 30.83 57.20 64.73 61.73 66.08 59.26 G-Δ 67.14 30.77 54.97 58.90 60.06 65.46 64.34 DSS-Hybrid 88.57 40.77 61.21 64.81 62.73 66.76 66.43 4.5 Standalone versus hybrid: cross-evaluation and deployment guidance Each mode is also evaluated on the other’s tasks under the same protocols. Table 5 summarizes the outcome, and Tables 24–26 in Appendix G give the per-dataset numbers with standard deviations and ID accuracies (protocols in Appendix G.1). Table 5: Both deployment modes on all three tasks. Each cell summarizes the full per-dataset comparison in Appendix G (Tables 24, 25, and 26); the standalone OOD and GOOD cells and the hybrid calibration cell are the cross-evaluations; the other cells restate Tables 1, 2–3, and 4. Task Standalone DSS-GNN DSS-Hybrid Calibration (14 datasets) Best Brier on all 14 and best accuracy on 13 (Table 1); stronger than the hybrid on 13 of 14 Never loses beyond noise to its own GCN base; Roman-Empire +27.7+27.7 accuracy, −0.33-0.33 Brier OOD detection (11 settings) Effective on 10 of 11: Cora 79.8/88.2/94.079.8/88.2/94.0, Photo 96.596.5–98.198.1, CS 94.394.3–97.897.8, Arxiv 73.773.7 (margin objective); the Twitch energy score does not separate Best published AUROC on 7 of 11 (Table 2); stronger than the standalone on 10 of 11 GOOD concept shift (7 settings) Above every published baseline on 5 of 7, at ERM level on the other 2 Strongest on all 7 (Table 4) Classification and calibration. The standalone is the stronger calibrated classifier on 13 of 14 datasets; Tolokers is the one exception (hybrid 81.281.2 accuracy / 0.2550.255 Brier versus 79.979.9 / 0.2680.268). Where the GCN encoder is strong the two are close (Cora 86.686.6 vs 84.384.3); on heterophilous graphs the DSS residual compensates for a weak encoder (Roman-Empire 78.878.8 vs 74.474.4); and where a plain GCN is unsuitable the standalone is higher (Texas 91.391.3 vs 45.645.6, Squirrel 68.168.1 vs 26.426.4). Against its own identically trained GCN base, the hybrid never loses beyond noise on either metric (worst accuracy delta −0.8-0.8 on Texas, inside one standard deviation of 7.77.7) and gains where the encoder is weak (Roman-Empire +27.7+27.7 accuracy / −0.33-0.33 Brier, Minesweeper +6.6+6.6 / −0.11-0.11): the residual never measurably reduces its encoder’s calibration or clean accuracy, and where the base is unsuitable (small heterophilous graphs) the hybrid matches it rather than improving on it. OOD detection. The hybrid is the stronger detector on 10 of 11 settings. Its GCN-plus-residual encoder is stronger in the sparse-label perturbation regime (Cora-structure 94.394.3 vs 79.879.8) and for cross-graph transfer (Twitch), while on Arxiv the standalone run under the GNNSafe++ margin objective is 0.40.4 higher with the best Arxiv ID accuracy (61.261.2 vs 54.054.0). The one standalone exception is Twitch: the classifier trains (ID accuracy 68.468.4 vs 70.570.5) but its energy score does not separate the cross-graph OOD inputs, so the hybrid (95.895.8) is preferred there. On the perturbation settings with tiny public splits the GNNSafe pipeline gives every backbone BatchNorm, and plain GCN in the same pipeline does not train without it (ID accuracy 75.775.7 with BatchNorm, 45.845.8 without, on Cora-structure); the standalone uses one BatchNorm1d per chaos channel (Appendix D) in all OOD and GOOD cells but OOD Twitch and GOOD-Arxiv/degree, and does not need it under the full-label calibration protocol (Appendix G.6). Distribution shift. On GOOD the standalone is above every baseline of Table 4 on 5 of 7 settings (WebKB 36.436.4 and Twitch 60.360.3 vs TAR 30.830.8 and 57.257.2) and at ERM’s level on the other two (CBAS 82.482.4, Arxiv/time 65.465.4), with ID accuracy between 63.563.5 and 96.296.2 (Appendix G); the hybrid is above every baseline on all 7. Deployment guidance. The two modes share one uncertainty representation and the same two readouts, (8) and (9), and differ in the encoder. Outside the exceptions above, each is the better mode on its deployment task and remains effective on the other’s. We recommend the standalone for calibration-critical single-graph use, in particular under heterophily, and the hybrid when a strong conventional encoder exists, for cross-graph transfer, and for OOD detection in sparse-label regimes; the standalone requires BatchNorm stabilization in sparse-label regimes, and the hybrid shares its encoder’s failures (Section 5). 5 Conclusion and Limitations DSS-GNN is one uncertainty representation, a doubly-spectral stochastic expansion over graph-frequency and chaos-order axes, with task-matched readouts and two deployment modes. Empirically, the standalone has the lowest Brier score among the compared uncertainty-aware baselines on all 14 calibration benchmarks, the hybrid the strongest AUROC on most node-OOD settings and the strongest shifted accuracy on all 7 concept-shift benchmarks under ERM, and both modes are effective on all three tasks with the exceptions noted in Section 4.5. Limitations. (i) Two deployment modes. The modes share one representation and the same readouts but are not interchangeable in every regime. The standalone model requires BatchNorm stabilization in sparse-label regimes, where every baseline backbone in the OOD pipeline also uses it, and uses it on 6 of 7 GOOD settings; under the full-label calibration protocol it is not needed and is not used (Appendix G.6). The hybrid is preferable when a strong conventional encoder exists and shares that encoder’s limits: on small heterophilous graphs the hybrid matches its base without improving on it. The standalone energy score does not separate cross-graph OOD inputs on Twitch, where only the hybrid is effective. (i) Scalar latent. All stochastic variation is due to one shared Gaussian factor, perfectly dependent across nodes, so the logit covariance across nodes has rank at most P. This suffices for the marginal readouts used here (per-node calibration and energy scores) and is what makes single-pass quadrature possible, but it cannot represent independent per-node noise, so joint cross-node uncertainty is not modeled; the sampling baselines (MC-dropout, Deep Ensembles) have higher Brier than the standalone and lower AUROC than the hybrid (Appendix F, Table 2). A vector latent of dimension q would cost SqS^q quadrature nodes on a full tensor grid, which sparse grids mitigate; we leave this to future work. (i) Compute. The chaos expansion adds a (P+1)×S(P+1)× S operation count per layer; the measured per-epoch overhead at P=2P=2 is 2.9×2.9× on Cora and 1.2×1.2× on Amazon-Ratings, and about 4%4\% per additional order on ogbn-arxiv, with memory linear in P+1P+1 (Table 20). Acknowledgments and Disclosure of Funding This work was supported in part by the National Science Foundation under grant 2531008. Part of this work was done while Fred Xu was an intern at Block, Inc., and we thank Block, Inc. for its support. References [1] M. Arjovsky, L. Bottou, I. Gulrajani, and D. Lopez-Paz (2019) Invariant risk minimization. arXiv preprint arXiv:1907.02893. External Links: Document, Link Cited by: 2nd item, Table 6, §1, §4.1. [2] D. Bo, X. Wang, C. Shi, and H. Shen (2021) Beyond low-frequency information in graph convolutional networks. Proceedings of the AAAI Conference on Artificial Intelligence 35 (5), p. 3950–3957. External Links: Document, Link Cited by: §2. [3] G. W. Brier (1950) Verification of forecasts expressed in terms of probability. Monthly Weather Review 78 (1), p. 1–3. External Links: Document, Link Cited by: §4.1. [4] R. H. Cameron and W. T. Martin (1947) The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals. Annals of Mathematics 48 (2), p. 385–392. External Links: Document, Link Cited by: §A.1, §2. [5] E. Chien, J. Peng, P. Li, and O. Milenkovic (2021) Adaptive universal generalized PageRank graph neural network. In International Conference on Learning Representations, External Links: Link Cited by: §A.3. [6] M. Defferrard, X. Bresson, and P. Vandergheynst (2016) Convolutional neural networks on graphs with fast localized spectral filtering. In Advances in Neural Information Processing Systems, External Links: Link Cited by: item 1, §A.3, §A.3, §1, §2. [7] Y. Dou, Z. Liu, L. Sun, Y. Deng, H. Peng, and P. S. Yu (2020) Enhancing graph neural network-based fraud detectors against camouflaged fraudsters. In Proceedings of the 29th ACM International Conference on Information & Knowledge Management, p. 315–324. External Links: Document, Link Cited by: §1. [8] R. Duan, M. Guang, J. Wang, C. Yan, H. Qi, W. Su, C. Tian, and H. Yang (2024) Unifying homophily and heterophily for spectral graph neural networks via triple filter ensembles. In Advances in Neural Information Processing Systems, External Links: Document, Link Cited by: 1st item, Table 6, §4.1. [9] D. Fuchsgruber, T. Wollschläger, and S. Günnemann (2024) Energy-based epistemic uncertainty for graph neural networks. In Advances in Neural Information Processing Systems, External Links: Document, Link Cited by: 3rd item, Table 6, §D.3, §1, §4.1, Table 2, Table 3. [10] Y. Gal and Z. Ghahramani (2016) Dropout as a Bayesian approximation: representing model uncertainty in deep learning. In Proceedings of the 33rd International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 48, p. 1050–1059. External Links: Link Cited by: 4th item, Table 6, §1, §4.1. [11] G. H. Golub and J. H. Welsch (1969) Calculation of Gauss quadrature rules. Mathematics of Computation 23 (106), p. 221–230. External Links: Document, Link Cited by: §A.1. [12] S. Gui, X. Li, L. Wang, and S. Ji (2022) GOOD: a graph out-of-distribution benchmark. In Advances in Neural Information Processing Systems, External Links: Link Cited by: §C.3, §4.1, §4.4. [13] M. He, Z. Wei, Z. Huang, and H. Xu (2021) BernNet: learning arbitrary graph spectral filters via Bernstein approximation. In Advances in Neural Information Processing Systems, External Links: Link Cited by: §A.3. [14] S. Janson (1997) Gaussian hilbert spaces. Cambridge University Press. External Links: Document, Link Cited by: §A.1, §A.1, §A.1, §A.1, §A.1, §A.2, §A.5. [15] T. N. Kipf and M. Welling (2017) Semi-supervised classification with graph convolutional networks. In International Conference on Learning Representations, External Links: Link Cited by: §A.3, Appendix F, §1, §4.1. [16] D. Krueger, E. Caballero, J. Jacobsen, A. Zhang, J. Binas, D. Zhang, R. Le Priol, and A. Courville (2021) Out-of-distribution generalization via risk extrapolation (REx). In Proceedings of the 38th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 139, p. 5815–5826. External Links: Link Cited by: 4th item. [17] B. Lakshminarayanan, A. Pritzel, and C. Blundell (2017) Simple and scalable predictive uncertainty estimation using deep ensembles. In Advances in Neural Information Processing Systems, External Links: Link Cited by: 4th item, Table 6, §1, §4.1. [18] Y. Li, R. Yu, C. Shahabi, and Y. Liu (2018) Diffusion convolutional recurrent neural network: data-driven traffic forecasting. In International Conference on Learning Representations, External Links: Link Cited by: §1. [19] X. Lin, W. Zhang, F. Shi, C. Zhou, L. Zou, X. Zhao, D. Yin, S. Pan, and Y. Cao (2024) Graph neural stochastic diffusion for estimating uncertainty in node classification. In Proceedings of the 41st International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 235, p. 30457–30478. External Links: Link Cited by: §1. [20] W. Liu, X. Wang, J. Owens, and Y. Li (2020) Energy-based out-of-distribution detection. In Advances in Neural Information Processing Systems, External Links: Link Cited by: §3.2. [21] E. Nelson (1973) The free Markoff field. Journal of Functional Analysis 12 (2), p. 211–227. External Links: Document, Link Cited by: §A.1. [22] D. Nualart (2006) The malliavin calculus and related topics. 2nd edition, Probability and Its Applications, Springer Berlin Heidelberg. External Links: Document, Link Cited by: §A.1. [23] S. Oladyshkin, T. Praditia, I. Kröker, F. Mohammadi, W. Nowak, and S. Otte (2023) The deep arbitrary polynomial chaos neural network or how Deep Artificial Neural Networks could benefit from data-driven homogeneous chaos theory. Neural Networks 166, p. 85–104. External Links: Document, Link Cited by: §2. [24] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, and P. Vandergheynst (2018) Graph signal processing: overview, challenges, and applications. Proceedings of the IEEE 106 (5), p. 808–828. External Links: Document, Link Cited by: §2. [25] H. Pei, B. Wei, K. C. Chang, Y. Lei, and B. Yang (2020) Geom-GCN: geometric graph convolutional networks. In International Conference on Learning Representations, External Links: Link Cited by: §4.1. [26] O. Platonov, D. Kuznedelev, M. Diskin, A. Babenko, and L. Prokhorenkova (2023) A critical look at the evaluation of GNNs under heterophily: are we really making progress?. In International Conference on Learning Representations, External Links: Link Cited by: Table 8, §D.1, §4.1. [27] S. Sagawa, P. W. Koh, T. B. Hashimoto, and P. Liang (2020) Distributionally robust neural networks. In International Conference on Learning Representations, External Links: Link Cited by: 3rd item, Table 6, §1, §4.1. [28] O. Shchur, M. Mumme, A. Bojchevski, and S. Günnemann (2018) Pitfalls of graph neural network evaluation. In Relational Representation Learning Workshop, NeurIPS 2018, External Links: Document, Link Cited by: §4.1. [29] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst (2013) The emerging field of signal processing on graphs: extending high-dimensional data analysis to networks and other irregular domains. IEEE Signal Processing Magazine 30 (3), p. 83–98. External Links: Document, Link Cited by: item 1, §A.3, §A.4, §A.6, §2. [30] M. Stadler, B. Charpentier, S. Geisler, D. Zügner, and S. Günnemann (2021) Graph Posterior Network: Bayesian predictive uncertainty for node classification. In Advances in Neural Information Processing Systems, External Links: Link Cited by: 5th item, Table 6, §1, §4.1, Table 2. [31] J. M. Stokes, K. Yang, K. Swanson, W. Jin, A. Cubillos-Ruiz, N. M. Donghia, C. R. MacNair, S. French, L. A. Carfrae, Z. Bloom-Ackermann, et al. (2020) A deep learning approach to antibiotic discovery. Cell 180 (4), p. 688–702.e13. External Links: Document, Link Cited by: §1. [32] G. Szegő (1975) Orthogonal polynomials. 4th edition, American Mathematical Society. Cited by: §A.1, §A.1, §A.1, §A.6. [33] P. Trivedi, M. Heimann, R. Anirudh, D. Koutra, and J. J. Thiagarajan (2024) Accurate and scalable estimation of epistemic uncertainty for graph neural networks. In International Conference on Learning Representations, External Links: Link Cited by: 2nd item, Table 6, §1, §4.1. [34] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio (2018) Graph attention networks. In International Conference on Learning Representations, External Links: Link Cited by: Appendix F. [35] N. Wiener (1938) The homogeneous chaos. American Journal of Mathematics 60 (4), p. 897–936. External Links: Document, Link Cited by: §A.1, §2. [36] Q. Wu, Y. Chen, C. Yang, and J. Yan (2023) Energy-based out-of-distribution detection for graph neural networks. In International Conference on Learning Representations, External Links: Link Cited by: 1st item, 2nd item, 5th item, Table 6, Table 6, §D.3, §G.1, §1, §3.2, §4.1, §4.3, Table 2, Table 3. [37] D. Xiu and G. E. Karniadakis (2002) The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing 24 (2), p. 619–644. External Links: Document, Link Cited by: §A.1, §A.1, §A.2, §A.2, §A.2, §A.3, §A.3, §1, §2, Remark 3. [38] W. Zheng, J. Liu, J. Li, J. Wu, P. Cui, and Y. Kong (2025) Topology-aware dynamic reweighting for distribution shifts on graph. In Proceedings of the 42nd International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 267, p. 78260–78275. External Links: Link Cited by: 5th item, §B.1, Table 6, §D.4, §1, §4.1, Table 4. Appendix: Table of Contents 1 Introduction 2 Background 3 Method 3.1 DSS layer 3.2 Uncertainty readout 3.3 Theoretical properties 3.4 Hybrid deployment 3.5 Training and downstream tasks 4 Experiments 4.1 Setup 4.2 Node classification and calibration 4.3 OOD detection 4.4 Distribution-shifted classification 4.5 Standalone versus hybrid: cross-evaluation and deployment guidance 5 Conclusion and Limitations References A Proofs and Derivations A.1 Background: Wiener chaos and Gauss–Hermite quadrature A.2 Supporting technical results A.3 Proof of Proposition (Generalization) A.4 Proof of Theorem (Doubly-spectral operator) A.5 Proof of Theorem (Representational capacity and truncation) A.6 Proof of the Truncation Bound in Theorem A.7 Proof of Proposition (Stochastic and mean-logit readouts) B Baseline Details B.1 Baseline methods C Dataset Statistics and Licenses C.1 Calibration benchmarks C.2 OOD detection benchmarks C.3 GOOD distribution-shift benchmarks D Implementation Details D.1 General training protocol D.2 Calibration (Table ) D.3 OOD detection (Section ) D.4 Distribution shift (Table ) D.5 Software, hardware, and compute budget E Ablation Study Details E.1 Polynomial order P across all datasets E.2 Chaos order on OOD detection and under distribution shift E.3 Energy-score propagation E.4 Chaos energy regularization λreg _reg E.5 Number of quadrature nodes S E.6 Contribution decomposition: backbone vs. chaos E.7 Computational overhead E.8 Summary of component evidence F Comparison with Standard GNN Baselines G Cross-Evaluation of the Two Deployment Modes G.1 Protocols G.2 Calibration: standalone, hybrid, and the hybrid’s GCN base G.3 OOD detection: standalone versus hybrid G.4 Distribution shift: standalone versus hybrid G.5 Readout coincidence G.6 BatchNorm under the calibration protocol G.7 Covariate shift on GOOD H Architecture Overview and Notation Appendix A Proofs and Derivations We begin with a self-contained introduction to Wiener chaos and Gauss–Hermite quadrature (Appendix A.1), then present supporting technical results (Appendix A.2), followed by individual proofs of each main-text theorem (Appendices A.3–A.7). A.1 Background: Wiener chaos and Gauss–Hermite quadrature This section collects the probability and approximation theory underlying the doubly-spectral expansion. The material is classical; we follow the presentation in Janson [14] and Xiu and Karniadakis [37]. This appendix is self-contained. Gaussian measure and Hermite polynomials. Let γ denote the standard Gaussian measure on ℝR, i.e., dγ(ω)=(2π)−1/2e−ω2/2dωdγ(ω)=(2π)^-1/2e^-ω^2/2\,dω. The Hilbert space L2(ℝ,γ)L^2(R,γ) consists of all square-integrable functions with respect to γ, equipped with the inner product ⟨f,g⟩=∫ℝf(ω)g(ω)γ(ω)=[f(ω)g(ω)] f,g = _Rf(ω)\,g(ω)\,dγ(ω)=E[f(ω)\,g(ω)] where ω∼(0,1)ω (0,1). The probabilists’ Hermite polynomials [32] are defined by the Rodrigues formula Hen(ω)=(−1)neω2/2dndωne−ω2/2,n=0,1,2,…He_n(ω)=(-1)^n\,e^ω^2/2\, d^ndω^n\,e^-ω^2/2, n=0,1,2,… (12) The first few are He0(ω)=1He_0(ω)=1, He1(ω)=ωHe_1(ω)=ω, He2(ω)=ω2−1He_2(ω)=ω^2-1, He3(ω)=ω3−3ωHe_3(ω)=ω^3-3ω. They satisfy the three-term recurrence Hen+1(ω)=ωHen(ω)−nHen−1(ω),He_n+1(ω)=ω\,He_n(ω)-n\,He_n-1(ω), (13) with He−1=0He_-1=0, and the orthogonality relation [Hem(ω)Hen(ω)]=n!δmnE[He_m(ω)\,He_n(ω)]=n!\, _mn. Normalizing gives the orthonormal Hermite basis Ψn(ω):=Hen(ω)n!,[ΨmΨn]=δmn, _n(ω):= He_n(ω) n!, [ _m _n]= _mn, (14) which forms a complete orthonormal system for L2(ℝ,γ)L^2(R,γ) [32, 14]. Every f∈L2(ℝ,γ)f∈ L^2(R,γ) can therefore be expanded as f(ω)=∑n=0∞fnΨn(ω)f(ω)= _n=0^∞f_n\, _n(ω) with fn=[fΨn]f_n=E[f _n] and Parseval’s identity ‖f‖L2(γ)2=∑nfn2\|f\|_L^2(γ)^2= _nf_n^2. Wiener–Itô chaos decomposition. The Hermite expansion generalizes to a fundamental structural result in Gaussian analysis. Let (Ω,ℱ,ℙ)( ,F,P) be a probability space carrying a standard Gaussian random variable ω, and define the associated Gaussian subspace Lω2:=f(ω):f∈L2(ℝ,γ)⊆L2(Ω,ℙ).L_ω^2:=\f(ω):f∈ L^2(R,γ)\ L^2( ,P). Define the n-th Wiener chaos nC_n as the closed linear span of Ψn(ω) _n(ω) in Lω2L_ω^2. Since we use a single Gaussian latent (q=1q=1), each nC_n is one-dimensional and spanned by Ψn(ω) _n(ω); with q>1q>1, the corresponding chaos subspace would be indexed by Hermite multi-indices. The Wiener–Itô chaos decomposition states [35, 4, 14]: Lω2=⨁n=0∞n,L_ω^2= _n=0^∞C_n, (15) that is, every square-integrable random variable measurable with respect to ω admits a unique, orthogonal expansion in the Hermite basis. Denoting the orthogonal projection onto nC_n by PnP_n, any f∈Lω2f∈ L_ω^2 satisfies f=∑n=0∞Pnf= _n=0^∞P_nf with ‖f‖2=∑n‖Pnf‖2\|f\|^2= _n\|P_nf\|^2. The zeroth chaos 0C_0 consists of constants (P0f=[f]P_0f=E[f]); the first chaos is the space of centered Gaussian random variables; higher chaoses capture progressively more complex nonlinear functionals of ω. Remark 1 (Analogy with graph spectral decomposition). The decomposition (15) is the stochastic counterpart of the graph spectral decomposition ℓ2(V)=⨁j=1Nspan(uj) ^2(V)= _j=1^Nspan(u_j). Both decompose a Hilbert space into eigenspaces of a self-adjoint operator: the graph Laplacian LGL_G for ℓ2(V) ^2(V), and the Ornstein–Uhlenbeck operator ℒL (defined below) for the Gaussian Hilbert space Lω2L_ω^2. Remark 2 (Completeness of the doubly-spectral basis). Take the signal domain to be the random graph signals generated by the shared latent: a random variable per node with finite second moment, measurable with respect to ω, i.e. an element of the tensor-product Hilbert space ℓ2(V)⊗L2(ℝ,γ) ^2(V) L^2(R,γ). The graph Fourier eigenvectors ujj=1N\u_j\_j=1^N form an orthonormal basis of ℓ2(V) ^2(V) and, by the Wiener–Itô decomposition (15), the normalized Hermite polynomials Ψnn≥0\ _n\_n≥ 0 form an orthonormal basis of L2(ℝ,γ)L^2(R,γ). The tensor product of orthonormal bases is an orthonormal basis of the tensor-product space, so uj⊗Ψnj≤N,n≥0\u_j _n\_j≤ N,\,n≥ 0 is a complete orthonormal basis of this domain, and the double expansion (1) is a bi-spectral decomposition in the strict sense: both factors are eigenbases of self-adjoint operators (Remark 1). The scope caveat is that the space is built over the scalar ω: signals with independent per-node latents are outside it (Section 5). The Ornstein–Uhlenbeck operator and hypercontractivity. The Ornstein–Uhlenbeck (O-U) operator ℒL is defined on its standard dense domain in L2(ℝ,γ)L^2(R,γ) by ℒf(ω)=f′(ω)−ωf′(ω),Lf(ω)=f (ω)-ω\,f (ω), (16) with eigenvalues and eigenfunctions ℒΨn=−nΨnL _n=-n\, _n [22]. This is directly analogous to the graph Laplacian having eigenpairs (uj,λj)(u_j, _j): −ℒ-L is the corresponding operator on L2(ℝ,γ)L^2(R,γ), with the chaos orders 0,1,2,…0,1,2,… as eigenvalues. The associated O-U semigroup Tt=etℒT_t=e^tL acts as TtΨn=e−ntΨnT_t _n=e^-nt _n, damping higher-order chaos components exponentially. A fundamental property of TtT_t is hypercontractivity: for 1<p≤q1<p≤ q and e2t≥(q−1)/(p−1)e^2t≥(q-1)/(p-1), ‖Ttf‖Lq(γ)≤‖f‖Lp(γ)\|T_tf\|_L^q(γ)≤\|f\|_L^p(γ) (17) for all f∈Lp(γ)f∈ L^p(γ) [21, 14]. A direct consequence for chaos components is the moment comparison inequality: for every g∈ng _n and q≥2q≥ 2, ‖g‖Lq(γ)≤(q−1)n/2‖g‖L2(γ).\|g\|_L^q(γ)≤(q-1)^n/2\,\|g\|_L^2(γ). (18) This follows by applying (17) with p=2p=2 and ρ=e−t=(q−1)−1/2ρ=e^-t=(q-1)^-1/2 to g∈ng _n, using Ttg=e−ntgT_tg=e^-ntg [14, Corollary 5.11]. Equation (18) bounds the higher moments of a fixed chaos component; the exponential decay of chaos coefficients with order n requires additional regularity and is established separately in Theorem 2 via a Cauchy-estimate argument. Gauss–Hermite quadrature. To evaluate integrals against γ numerically, we use Gauss–Hermite quadrature [11, 37]. An S-point rule (ωs,μs)s=1S\( _s, _s)\_s=1^S approximates ∫ℝf(ω)γ(ω)≈∑s=1Sμsf(ωs), _Rf(ω)\,dγ(ω)\;≈\; _s=1^S _s\,f( _s), (19) where the nodes ω1<⋯<ωS _1<·s< _S are the roots of HeS(ω)He_S(ω) and the weights μs>0 _s>0 sum to one. The main property is polynomial exactness: the approximation is exact whenever f is a polynomial of degree at most 2S−12S-1 [32]. In particular: • Orthonormality: ∑s=1SμsΨm(ωs)Ψn(ωs)=δmn _s=1^S _s\, _m( _s)\, _n( _s)= _mn for all m+n≤2S−1m+n≤ 2S-1. • Triple products: ∑s=1SμsΨr(ωs)Ψn(ωs)Ψm(ωs)=[ΨrΨnΨm] _s=1^S _s\, _r( _s) _n( _s) _m( _s)=E[ _r _n _m] whenever r+n+m≤2S−1r+n+m≤ 2S-1. For DSS-GNN with chaos order P and gate order PgP_g, the most demanding integrand is the triple product ΨrΨnΨm _r _n _m with r≤Pgr≤ P_g, n,m≤Pn,m≤ P, which has degree Pg+2P_g+2P. Exactness therefore requires S≥⌈(Pg+2P+1)/2⌉S≥ (P_g+2P+1)/2 quadrature nodes. For Pg≤1P_g≤ 1 the threshold is S≥P+1S≥ P+1; the experiments use S=4S=4, which is exact for P≤3P≤ 3. Remark 3 (Non-intrusive vs. intrusive projection). In the polynomial chaos literature [37], two strategies compute the chaos coefficients of a transformed variable σ(H(ω))σ(H(ω)): (i) intrusive methods analytically propagate the expansion through σ using triple-product coefficients (Lemma 1); (i) non-intrusive methods evaluate σ at quadrature nodes and project back (19). DSS-GNN uses the non-intrusive approach (6) because it applies unchanged to any pointwise activation σ without requiring closed-form triple-product expansions. Theorem 3 proves exact agreement for product terms α(ω)H(ω)α(ω)H(ω) when the quadrature rule integrates the relevant triple products exactly; for a general nonlinear activation σ, the non-intrusive projection remains an approximation whose error is decomposed in Proposition 3. A.2 Supporting technical results Lemma 1 (Products in chaos: triple-product contraction). Let α(ω)=∑r=0PgarΨr(ω)α(ω)= _r=0^P_ga_r _r(ω) and H(ω)=∑n=0PHnΨn(ω)H(ω)= _n=0^PH_n _n(ω). Then the order-m coefficient of the product α(ω)H(ω)α(ω)H(ω) is (αH)m=∑r=0Pg∑n=0ParHncrnm,crnm:=[ΨrΨnΨm].(α H)_m= _r=0^P_g _n=0^Pa_rH_n\,c_rnm, c_rnm:=E[ _r _n _m]. Proof. Expand the product in the Hermite basis: α(ω)H(ω)=(∑r=0PgarΨr(ω))(∑n=0PHnΨn(ω))=∑r,narHnΨr(ω)Ψn(ω).α(ω)H(ω)= ( _r=0^P_ga_r _r(ω) ) ( _n=0^PH_n _n(ω) )= _r,na_rH_n\, _r(ω) _n(ω). The order-m coefficient is obtained by projecting onto Ψm _m using the orthonormality [ΨiΨj]=δijE[ _i _j]= _ij: (αH)m=[α(ω)H(ω)Ψm(ω)]=∑r,narHn[ΨrΨnΨm]=∑r,narHncrnm.(α H)_m=E [α(ω)H(ω) _m(ω) ]= _r,na_rH_n\,E[ _r _n _m]= _r,na_rH_n\,c_rnm. The triple-product constants crnm=[ΨrΨnΨm]c_rnm=E[ _r _n _m] are computable in closed form from the three-term recurrence of the probabilists’ Hermite polynomials [37]. In particular, crnm=0c_rnm=0 unless |r−n|≤m≤r+n|r-n|≤ m≤ r+n and r+n+mr+n+m is even, which makes the coupling matrices sparse. ∎ Theorem 3 (Non-intrusive projection recovers product coefficients when quadrature is exact). Let (αH)^m:=∑s=1Sμsα(ωs)H(ωs)Ψm(ωs) (α H)_m:= _s=1^S _s\,α( _s)H( _s) _m( _s) be the discrete projection of the product using a Gauss–Hermite quadrature rule (ωs,μs)s=1S\( _s, _s)\_s=1^S. If the quadrature integrates all triple products Ψr(ω)Ψn(ω)Ψm(ω) _r(ω) _n(ω) _m(ω) exactly for r≤Pgr≤ P_g, n≤Pn≤ P, m≤Pm≤ P, then (αH)^m=(αH)m (α H)_m=(α H)_m for all m≤Pm≤ P. Proof. Substituting the expansions of α and H: (αH)^m (α H)_m =∑s=1Sμsα(ωs)H(ωs)Ψm(ωs) = _s=1^S _s\,α( _s)H( _s) _m( _s) =∑s=1Sμs∑r=0Pg∑n=0ParHnΨr(ωs)Ψn(ωs)Ψm(ωs) = _s=1^S _s _r=0^P_g _n=0^Pa_rH_n\, _r( _s) _n( _s) _m( _s) =∑r,narHn∑s=1SμsΨr(ωs)Ψn(ωs)Ψm(ωs)⏟QS[ΨrΨnΨm]. = _r,na_rH_n _s=1^S _s\, _r( _s) _n( _s) _m( _s)_ Q_S[ _r _n _m]. The integrand ΨrΨnΨm _r _n _m is a polynomial of degree r+n+m≤Pg+2Pr+n+m≤ P_g+2P under the Gaussian weight γ. An S-point Gauss–Hermite rule integrates polynomials of degree up to 2S−12S-1 exactly [37]. The condition 2S−1≥Pg+2P2S-1≥ P_g+2P, i.e. S≥⌈(Pg+2P+1)/2⌉S≥ (P_g+2P+1)/2 , therefore ensures QS[ΨrΨnΨm]=[ΨrΨnΨm]=crnmQ_S[ _r _n _m]=E[ _r _n _m]=c_rnm. By Lemma 1, (αH)^m=(αH)m (α H)_m=(α H)_m. ∎ Proposition 3 (Discrete projection error decomposition). Let ΠP _P be the exact L2(ℝ,γ)L^2(R,γ) projection onto the first P+1P+1 Wiener chaos orders, and let ΠP,SQ _P,SQ be the finite chaos expansion whose coefficients are induced by a quadrature rule: (ΠP,SQ)m:=∑s=1SμsQ(ωs)Ψm(ωs)( _P,SQ)_m:= _s=1^S _s\,Q( _s) _m( _s). Then for any Q∈L2(ℝ,γ)Q∈ L^2(R,γ), ‖Q−ΠP,SQ‖L2≤‖Q−ΠPQ‖L2⏟truncation error+‖ΠPQ−ΠP,SQ‖L2⏟quadrature error.\|Q- _P,SQ\|_L^2\;≤\; \|Q- _PQ\|_L^2_truncation error+ \| _PQ- _P,SQ\|_L^2_quadrature error. Moreover, Gauss–Hermite rules yield exact coefficients whenever Q(ω)Ψm(ω)Q(ω) _m(ω) is a polynomial of sufficiently low degree under the Gaussian weight. Proof. The triangle inequality in L2(ℝ,γ)L^2(R,γ) gives ‖Q−ΠP,SQ‖≤‖Q−ΠPQ‖+‖ΠPQ−ΠP,SQ‖\|Q- _P,SQ\|≤\|Q- _PQ\|+\| _PQ- _P,SQ\|. For the truncation term, ΠP _P is the orthogonal projection onto ⨁n=0Pn _n=0^PC_n, the direct sum of the first P+1P+1 Wiener chaos subspaces [14], so ‖Q−ΠPQ‖2=∑n>P‖Qn‖2\|Q- _PQ\|^2= _n>P\|Q_n\|^2 where Qn=[QΨn]Q_n=E[Q _n] are the exact chaos coefficients. For the quadrature term, the m-th coefficient of ΠP,SQ _P,SQ is ∑sμsQ(ωs)Ψm(ωs) _s _sQ( _s) _m( _s), which approximates ∫QΨmγ Q _m\,dγ. When QΨmQ _m is a polynomial of degree at most 2S−12S-1, Gauss–Hermite quadrature integrates it exactly [37], making the quadrature error zero. In the DSS layer context, the pre-activation H~(ℓ+1)(ω) H^( +1)(ω) before applying σ is a polynomial of degree at most Pg+P_g+P in ω; after applying σ it is no longer polynomial, and the quadrature error reflects the approximation quality of the S-point rule on this non-polynomial integrand. ∎ A.3 Proof of Proposition 1 (Generalization) See 1 Proof. When P=0P=0, the chaos expansion (2) contains a single term: H(ℓ)(ω)=H0(ℓ)Ψ0(ω)H^( )(ω)=H_0^( ) _0(ω). Since the zeroth normalized Hermite polynomial satisfies Ψ0≡1 _0≡ 1 [37], the representation is deterministic: H(ℓ)(ω)=H0(ℓ)H^( )(ω)=H_0^( ) for all ω∈Ωω∈ . Step 1 (Gate reduction). With P=0P=0, the mixing weights (4) are indexed by n=0n=0 only. Setting the gate order Pg=0P_g=0 (no stochastic gating when the representation itself is deterministic) yields scalar constants α0=a0,0∈ℝ _0=a_0,0 and β0=b0,0∈ℝ _0=b_0,0 . Step 2 (Pre-activation collapse). The quadrature-node pre-activation (5) reduces to H~(ℓ+1)(ωs) H^( +1)( _s) =Ψ0(ωs)(α0Glp(LG)H0(ℓ)+β0Ghp(LG)H0(ℓ))Wℓ = _0( _s) ( _0\,G_lp(L_G)H_0^( )+ _0\,G_hp(L_G)H_0^( ) )W_ =(α0Glp(LG)+β0Ghp(LG))H0(ℓ)Wℓ, = ( _0\,G_lp(L_G)+ _0\,G_hp(L_G) )H_0^( )W_ , which is independent of the quadrature node ωs _s. Step 3 (Projection collapse). The projection (6) for m=0m=0 becomes H0(ℓ+1)=∑s=1Sμsσ(H~(ℓ+1)(ωs))Ψ0(ωs)=σ(H~(ℓ+1))∑s=1Sμs=σ(H~(ℓ+1)),H_0^( +1)= _s=1^S _s\,σ\! ( H^( +1)( _s) ) _0( _s)=σ\! ( H^( +1) ) _s=1^S _s=σ\! ( H^( +1) ), where the second equality uses the fact that H~(ℓ+1) H^( +1) does not depend on ωs _s (Step 2), and the third uses the normalization ∑sμs=∫ℝγ=1 _s _s= _Rdγ=1 of the Gauss–Hermite weights [37]. Step 4 (Recovery of spectral GNN). Combining Steps 2 and 3, the full DSS layer update at P=0P=0 is H0(ℓ+1)=σ(g(LG)H0(ℓ)Wℓ),g(LG):=α0Glp(LG)+β0Ghp(LG).H_0^( +1)=σ\! (g(L_G)\,H_0^( )\,W_ ), g(L_G):= _0\,G_lp(L_G)+ _0\,G_hp(L_G). Since GlpG_lp and GhpG_hp are Chebyshev polynomials in LGL_G (3), their linear combination g(LG)g(L_G) is a polynomial in LGL_G of degree max(Klp,Khp) (K_lp,K_hp). This is exactly the form of a polynomial spectral graph convolution layer [6, 29]: the filter g acts in the graph-frequency domain as g(λj)=α0glp(λj)+β0ghp(λj)g( _j)= _0\,g_lp( _j)+ _0\,g_hp( _j), where glp(λ)=∑kck(lp)Tk(λ~)g_lp(λ)= _kc_k^(lp)T_k( λ) is the scalar frequency response of the low-pass branch and λ~=2λ/λmax−1 λ=2λ/ _ -1. Special cases. Setting β0=0 _0=0 and Klp=1K_lp=1 with the renormalization trick of Kipf and Welling [15] yields GCN as a single-layer DSS-GNN with P=0P=0. Allowing arbitrary Chebyshev degree with both filter branches active recovers ChebNet [6] in its scalar-coefficient form (filter weights Θk=ckWℓ _k=c_kW_ ). With generalized polynomial coefficient learning, GPR-GNN [5] and BernNet [13] are similarly obtained as P=0P=0 specializations. ∎ A.4 Proof of Theorem 1 (Doubly-spectral operator) Lemma 2 (Exact-projection conditions for the joint transfer). In Theorem 1, assume the pointwise nonlinearity is omitted, the gates have chaos degree at most PgP_g, the state has chaos degree at most P, and the Gauss–Hermite rule integrates all triple products ΨrΨnΨm _r _n _m with r≤Pgr≤ P_g and n,m≤Pn,m≤ P. A sufficient condition is S≥⌈(Pg+2P+1)/2⌉S≥ (P_g+2P+1)/2 . Under this condition, the stochastic coupling matrices in the joint-transfer formula are Mmn(ℓ,lp):=∑r=0Pgan,rcrnm,Mmn(ℓ,hp):=∑r=0Pgbn,rcrnm,crnm:=[ΨrΨnΨm].M_mn^( ,lp):= _r=0^P_ga_n,r\,c_rnm, M_mn^( ,hp):= _r=0^P_gb_n,r\,c_rnm, c_rnm:=E[ _r _n _m]. Proof. The product ΨrΨnΨm _r _n _m has degree at most Pg+2P_g+2P under the stated index ranges. An S-point Gauss–Hermite rule is exact for polynomials of degree at most 2S−12S-1, so 2S−1≥Pg+2P2S-1≥ P_g+2P suffices. The coupling matrix definitions are obtained by collecting the exactly projected triple-product coefficients in the linearized DSS update. ∎ See 1 Proof. Let ujj=1N\u_j\_j=1^N be the orthonormal eigenvectors of the graph Laplacian LGL_G with eigenvalues λjj=1N\ _j\_j=1^N, so that LG=∑j=1Nλjujuj⊤L_G= _j=1^N _j\,u_ju_j [29]. Because the Chebyshev filters (3) are polynomials in LGL_G, each eigenvector is simultaneously an eigenvector of both filters: Glp(LG)uj=glp(λj)uj,Ghp(LG)uj=ghp(λj)uj,G_lp(L_G)\,u_j=g_lp( _j)\,u_j, G_hp(L_G)\,u_j=g_hp( _j)\,u_j, (20) where glp(λ)=∑kck(lp)Tk(λ~)g_lp(λ)= _kc_k^(lp)T_k( λ) and ghp(λ)=∑kck(hp)Tk(λ~)g_hp(λ)= _kc_k^(hp)T_k( λ) are the scalar frequency responses. Step 1 (Graph Fourier transform of chaos coefficients). Define the joint spectral coefficient at graph frequency λj _j and chaos order n as H^n(ℓ)(λj):=uj⊤Hn(ℓ)∈ℝ1×dℓ H_n^( )( _j):=u_j H_n^( ) ^1× d_ . For proof convenience, write Unlp:=Glp(LG)Hn(ℓ)U_n^lp:=G_lp(L_G)H_n^( ) and Unhp:=Ghp(LG)Hn(ℓ)U_n^hp:=G_hp(L_G)H_n^( ). These spectrally filtered coefficients from (5) satisfy uj⊤Unlp=uj⊤Glp(LG)Hn(ℓ)=glp(λj)H^n(ℓ)(λj),u_j U_n^lp=u_j G_lp(L_G)\,H_n^( )=g_lp( _j)\; H_n^( )( _j), (21) and analogously uj⊤Unhp=ghp(λj)H^n(ℓ)(λj)u_j U_n^hp=g_hp( _j)\; H_n^( )( _j). Step 2 (Linearized update in the chaos basis). Omitting the pointwise nonlinearity σ (linearization), the DSS update (5)–(6) gives Hm(ℓ+1) H_m^( +1) =∑s=1SμsH~(ℓ+1)(ωs)Ψm(ωs) = _s=1^S _s\, H^( +1)( _s)\, _m( _s) =∑s=1Sμs∑n=0PΨn(ωs)(αn(ωs)Unlp+βn(ωs)Unhp)WℓΨm(ωs). = _s=1^S _s _n=0^P _n( _s) ( _n( _s)\,U_n^lp+ _n( _s)\,U_n^hp )W_ \; _m( _s). (22) Expanding the gates αn(ωs)=∑r=0Pgan,rΨr(ωs) _n( _s)= _r=0^P_ga_n,r _r( _s) from (4) and exchanging the order of summation: Hm(ℓ+1) H_m^( +1) =∑n=0P[(∑r=0Pgan,r∑s=1SμsΨr(ωs)Ψn(ωs)Ψm(ωs)⏟=crnm)Unlp = _n=0^P [ ( _r=0^P_ga_n,r _s=1^S _s\, _r( _s) _n( _s) _m( _s)_=\;c_rnm )U_n^lp +(∑r=0Pgbn,rcrnm)Unhp]Wℓ, \;+ ( _r=0^P_gb_n,r\,c_rnm )U_n^hp ]W_ , (23) where the identification ∑sμsΨr(ωs)Ψn(ωs)Ψm(ωs)=crnm:=[ΨrΨnΨm] _s _s\, _r( _s) _n( _s) _m( _s)=c_rnm:=E[ _r _n _m] follows from Theorem 3: the Gauss–Hermite rule integrates the triple product exactly when S is sufficiently large. Define the chaos coupling matrices M(ℓ,lp),M(ℓ,hp)∈ℝ(P+1)×(P+1)M^( ,lp),\,M^( ,hp) ^(P+1)×(P+1) with entries Mmn(ℓ,lp):=∑r=0Pgan,rcrnm,Mmn(ℓ,hp):=∑r=0Pgbn,rcrnm.M_mn^( ,lp):= _r=0^P_ga_n,r\,c_rnm, M_mn^( ,hp):= _r=0^P_gb_n,r\,c_rnm. (24) Then (23) simplifies to Hm(ℓ+1)=∑n=0P(Mmn(ℓ,lp)Unlp+Mmn(ℓ,hp)Unhp)Wℓ.H_m^( +1)= _n=0^P (M_mn^( ,lp)\,U_n^lp+M_mn^( ,hp)\,U_n^hp )\,W_ . (25) Step 3 (Joint spectral representation). Left-multiplying (25) by uj⊤u_j and applying (21): H^m(ℓ+1)(λj)=∑n=0P(Mmn(ℓ,lp)glp(λj)+Mmn(ℓ,hp)ghp(λj))H^n(ℓ)(λj)Wℓ. H_m^( +1)( _j)= _n=0^P (M_mn^( ,lp)\,g_lp( _j)+M_mn^( ,hp)\,g_hp( _j) )\; H_n^( )( _j)\;W_ . Stack the chaos orders into a matrix ^(ℓ)(λj):=[H^0(ℓ)(λj);…;H^P(ℓ)(λj)]∈ℝ(P+1)×dℓ H^( )( _j):= [ H_0^( )( _j);\;…\;;\; H_P^( )( _j) ] ^(P+1)× d_ . The equation above is the (m,⋅)(m,·) row of the matrix relation ^(ℓ+1)(λj)=ℋ(ℓ)(λj)^(ℓ)(λj)Wℓ, H^( +1)( _j)=H^( )( _j)\; H^( )( _j)\;W_ , (26) where the (P+1)×(P+1)(P+1)×(P+1) transfer matrix decomposes as ℋ(ℓ)(λ)=glp(λ)M(ℓ,lp)+ghp(λ)M(ℓ,hp).H^( )(λ)=g_lp(λ)\;M^( ,lp)+g_hp(λ)\;M^( ,hp). (27) Interpretation. Equation (26) shows that, in the linearized regime, the DSS layer implements a learned linear operator on the joint spectral domain of the doubly-spectral expansion (1). The decomposition (27) separates two independent mechanisms: 1. Graph-frequency response: the scalar functions glp(λ)g_lp(λ) and ghp(λ)g_hp(λ) shape how the layer responds to each graph eigenvalue λj _j, exactly as in standard spectral GNNs [6, 29]. 2. Chaos-order coupling: the matrices M(ℓ,lp)M^( ,lp) and M(ℓ,hp)M^( ,hp) control the coupling between chaos orders 0,…,P0,…,P, enabling the layer to redistribute energy across the chaos spectrum. Together, these two axes parameterize a bilinear family of operators on the tensor-product Hilbert space ℓ2(V)⊗L2(ℝ,γ) ^2(V) L^2(R,γ), thereby justifying the doubly-spectral designation. ∎ A.5 Proof of Theorem 2 (Representational capacity and truncation) See 2 Proof of the approximation claim. The proof proceeds in three stages: (i) truncation in the chaos-order axis, (i) a rank argument for the quadrature projection, and (i) a one-layer constructive realization using the final linear readout. Step 1 (Chaos-order truncation). By the Wiener–Itô chaos decomposition [14], any Y∈L2(ℝ,γ,ℝN×d)Y∈ L^2(R,γ;R^N× d) admits the expansion Y(ω)=∑n=0∞YnΨn(ω)Y(ω)= _n=0^∞Y_n _n(ω) with deterministic coefficients Yn=[YΨn]∈ℝN×dY_n=E[Y _n] ^N× d, and by Parseval’s identity ‖Y‖L22=∑n=0∞‖Yn‖F2<∞\|Y\|_L^2^2= _n=0^∞\|Y_n\|_F^2<∞. Fix a truncation order P for the construction. For the density claim, fix ε>0 >0 and take this P large enough that ∑n>P‖Yn‖F2<ε2 _n>P\|Y_n\|_F^2< ^2; for the explicit rate, keep P arbitrary and apply the tail estimate proved below. Step 2 (Quadrature projection rank). Take S=P+1S=P+1 and let (ωs,μs)s=1S\( _s, _s)\_s=1^S be the associated S-point Gauss–Hermite rule. Define the projection matrix Φ∈ℝ(P+1)×S ^(P+1)× S with entries Φns=μsΨn(ωs) _ns= _s\, _n( _s). The Hermite polynomials Ψ0,…,ΨP _0,…, _P are linearly independent on the S=P+1S=P+1 distinct quadrature nodes, and μs>0 _s>0 for all s, so Φ is invertible. Consequently, for any target coefficients Ynn=0P\Y_n\_n=0^P, there exist matrices Vss=1S\V_s\_s=1^S in ℝN×dR^N× d satisfying ∑s=1SμsVsΨn(ωs)=Yn,n=0,…,P. _s=1^S _s\,V_s\, _n( _s)=Y_n, n=0,…,P. (28) Step 3 (Constructive one-layer realization with ReLU). We construct a one-layer DSS-GNN (L=1L=1) whose post-readout output matches the truncated coefficients exactly. Let q:=2dq:=2d and define the nonnegative hidden targets Rs:=[Vs+,Vs−]∈ℝ+N×q,Vs=Vs+−Vs−,R_s:=[V_s^+,\,V_s^-] _+^N× q, V_s=V_s^+-V_s^-, where Vs+,Vs−V_s^+,V_s^- are the entrywise positive and negative parts of VsV_s. Choose the final readout Wout=[Id−Id]∈ℝq×d,W_out= bmatrixI_d\\[2.0pt] -I_d bmatrix ^q× d, so that RsWout=VsR_sW_out=V_s for every s. Next define the generalized Vandermonde matrix Θ∈ℝS×S,Θst:=Ψt−1(ωs),s,t=1,…,S. ^S× S, _st:= _t-1( _s), s,t=1,…,S. Because Ψ0,…,ΨS−1 _0,…, _S-1 are linearly independent and the quadrature nodes ω1,…,ωS _1,…, _S are distinct, Θ is invertible. Hence for the prescribed matrices RsR_s there exist unique matrices A1,…,AS∈ℝN×qA_1,…,A_S ^N× q such that ∑t=1SΨt−1(ωs)At=Rs,s=1,…,S. _t=1^S _t-1( _s)\,A_t=R_s, s=1,…,S. (29) Next, we realize these AtA_t with the DSS layer. Set the lifted input width to d0=Sqd_0=Sq and partition the lifted features into S blocks b1,…,bSb_1,…,b_S, each of width q. Choose the graph filters and gates as Glp(LG)=I(Klp=0,c0(lp)=1),Ghp(LG)=0,G_lp(L_G)=I (K_lp=0,\;c_0^(lp)=1), G_hp(L_G)=0, αn(ω)≡1(n=0,…,P),βn(ω)≡0, _n(ω)≡ 1\;\;(n=0,…,P), _n(ω)≡ 0, so Pg=0P_g=0 suffices. Since rank(X)=Nrank(X)=N, for each t=1,…,St=1,…,S there exists a matrix Bt∈ℝdin×qB_t ^d_ in× q with XBt=AtXB_t=A_t. Set Win(t−1)W_in^(t-1) to place BtB_t in block btb_t and zeros elsewhere. Then Ht−1(0)=XWin(t−1)H_t-1^(0)=XW_in^(t-1) is zero outside block btb_t and equals AtA_t on block btb_t. Finally, choose W0∈ℝd0×qW_0 ^d_0× q so that each row-block corresponding to btb_t equals the identity IqI_q. At quadrature node ωs _s, the quadrature-node pre-activation becomes H~(1)(ωs)=∑t=1SΨt−1(ωs)At=Rs H^(1)( _s)= _t=1^S _t-1( _s)\,A_t=R_s by (29). Since each RsR_s is entrywise nonnegative and σ=ReLUσ=ReLU, we obtain σ(H~(1)(ωs))=Rs,σ(H~(1)(ωs))Wout=Vs.σ\! ( H^(1)( _s) )=R_s, σ\! ( H^(1)( _s) )W_out=V_s. Projecting back to chaos coefficients gives, for every m=0,…,Pm=0,…,P, Hm(1)Wout=∑s=1Sμsσ(H~(1)(ωs))WoutΨm(ωs)=∑s=1SμsVsΨm(ωs)=Ym,H_m^(1)W_out= _s=1^S _s\,σ\! ( H^(1)( _s) )W_out _m( _s)= _s=1^S _s\,V_s\, _m( _s)=Y_m, where the last equality is exactly (28). Step 4 (Combining with the tail). The retained coefficients are matched exactly, so the order-P output Y~P Y_P lies in ⨁n=0Pn⊗ℝN×d _n=0^PC_n ^N× d and has no chaos coefficients for n>Pn>P. By orthogonality of the chaos subspaces and Parseval’s identity, it satisfies ‖Y−Y~P‖L22=∑n>P‖Yn‖F2.\|Y- Y_P\|_L^2^2= _n>P\|Y_n\|_F^2. Thus choosing P so that the right-hand side is below ε2 ^2 proves density in L2L^2. ∎ Remark 4 (Practical vs. theoretical capacity). The constructive part of Theorem 2 is an existence result for a restricted subfamily of DSS-GNNs: the proof requires rank(X)=Nrank(X)=N (i.e., din≥Nd_ in≥ N), a one-layer construction with lifted width d0=2dS=2d(P+1)d_0=2dS=2d(P+1) and output width d1=2d_1=2d, quadrature size S=P+1S=P+1, the ReLU activation, constant gates (Pg=0P_g=0), the identity low-pass filter, and a final linear readout that combines positive and negative hidden channels. The construction uses the identity low-pass filter and constant gates, so it bounds the capacity of the chaos axis alone; Theorem 1 describes the graph-frequency axis. The practical model differs from these conditions (dind_ in is typically smaller than N, P is small, and the filters are learned), so the theorem establishes representational capacity, while calibration, OOD detection, and robustness are established empirically. The explicit exponential bound adds that, for targets satisfying the entire-extension condition, a low truncation order should capture most of the gain, which is empirically consistent with the plateau at P=1–2P=1--2 observed in Figure 3(b). Remark 5 (Rank-deficient features). When rank(X)=r<Nrank(X)=r<N, the same construction still matches every truncated target whose chaos-coefficient columns lie in range(X)range(X), under one positivity condition required by the ReLU step: the nonnegative quadrature-node evaluations Rs=[Vs+,Vs−]R_s=[V_s^+,V_s^-] of Step 3 must themselves have columns in range(X)range(X), since At=XBtA_t=XB_t forces every AtA_t, and hence every RsR_s, into that range. The condition is satisfiable, for example, by appending a constant feature column to X: with Vs+:=Vs+c 11⊤V_s^+:=V_s+c\,11 and Vs−:=c 11⊤V_s^-:=c\,11 for c≥maxs‖Vs‖maxc≥ _s\|V_s\|_ , both parts are nonnegative and have columns in range([X,])⊇range(X)range([X,1]) (X), so every truncated target with columns in range([X,])range([X,1]) is matched. We state exactly this matching property; no claim is made for targets whose coefficients are outside that range. A.6 Proof of the Truncation Bound in Theorem 2 The analytic Gaussian-growth condition referenced in Theorem 2 is the following: each centered coordinate Y¯ik(ω):=Yik(ω)−[Yik] Y_ik(ω):=Y_ik(ω)-E[Y_ik] extends to an entire function satisfying |Y¯ik(z)|≤Aikexp(a|z|2)| Y_ik(z)|≤ A_ik (a|z|^2) for some a∈(0,14)a∈(0, 14) uniform in (i,k)(i,k), with ∑i,kAik2<∞ _i,kA_ik^2<∞. Under this condition, set α:=2a/(1−2a)<1α:=2a/(1-2a)<1; the proof yields the theorem for every β∈(α,1)β∈( α,1), with CtrC_ tr depending on a, β, and ∑i,kAik2 _i,kA_ik^2. Proof of the truncation bound. We prove the bound coordinate-wise and sum. Fix a centered coordinate function g:=Y¯ikg:= Y_ik with chaos coefficients cn:=[g(ω)Ψn(ω)]c_n:=E[g(ω) _n(ω)], so that g(ω)=∑n=1∞cnΨn(ω)g(ω)= _n=1^∞c_n _n(ω). Step 1 (Hermite generating function). The generating function of the normalized Hermite polynomials [32] gives ∑n=0∞tnn!Ψn(ω)=exp(tω−t2/2) _n=0^∞ t^n n! _n(ω)= (tω-t^2/2). Multiplying by g(ω)g(ω) and taking the Gaussian expectation yields ϕ(s):=∑n=1∞cnn!sn=[g(ω+s)]=∫ℝg(x+s)γ(x),φ(s)\;:=\; _n=1^∞ c_n n!\,s^n\;=\;E\! [g(ω+s) ]\;=\; _Rg(x+s)\,dγ(x), (30) which is verified by expanding g in the Hermite basis and using [Ψn(ω)exp(sω−s2/2)]=sn/n!E[ _n(ω) (sω-s^2/2)]=s^n/ n!. The series defining ϕφ converges for every complex s because |cn|≤‖g‖L2(γ)|c_n|≤\|g\|_L^2(γ), so ϕφ is entire; the right-hand side is entire by the growth bound below and dominated convergence, and the two sides agree for real s, hence for all s∈ℂs . Step 2 (Growth bound on ϕφ). The hypothesis |g(z)|≤Aexp(a|z|2)|g(z)|≤ A (a|z|^2) with a<14a< 14 gives |ϕ(s)| |φ(s)| ≤∫ℝAexp(a|x+s|2)γ(x)=Aea|s|2∫ℝexp(ax2+2axRe(s))γ(x). ≤ _RA\, \! (a|x+s|^2 )\,dγ(x)=A\,e^a|s|^2 _R \! (ax^2+2a\,xRe(s) )\,dγ(x). Completing the square in the exponent (ax2+2aRe(s)x−x2/2=−(12−a)x2+2aRe(s)xax^2+2aRe(s)\,x-x^2/2=-( 12-a)x^2+2aRe(s)\,x, which is a Gaussian integral convergent for a<12a< 12) yields |ϕ(s)|≤A1−2aexp(a1−2a|s|2).|φ(s)|\;≤\; A 1-2a\; \! ( a1-2a\,|s|^2 ). Step 3 (Cauchy estimate and Stirling bound). Since ϕφ is entire, Cauchy’s inequality on the circle |s|=R|s|=R gives |cn|n!≤R−nsup|s|=R|ϕ(s)|≤A1−2aR−nea′R2,a′:=a1−2a. |c_n| n!\;≤\;R^-n _|s|=R|φ(s)|\;≤\; A 1-2a\;R^-n\,e^a R^2, a := a1-2a. Setting R=n/(2a′)R= n/(2a ) (the minimizer of −nlnR+a′R2-n R+a R^2): |cn|≤A1−2an!(2a′n)n/2en/2.|c_n|\;≤\; A 1-2a\; n!\; ( 2a n )^n/2e^n/2. By Stirling’s inequality n!≤e1/2(2πn)1/4(n/e)n/2 n!≤ e^1/2(2π n)^1/4(n/e)^n/2, the product (n/e)n/2(2a′/n)n/2en/2(n/e)^n/2(2a /n)^n/2e^n/2 simplifies to (2a′)n/2(2a )^n/2, giving cn2≤A2e2πn1−2a(2a′)n=A2e2πn1−2aαn,α:=2a′=2a1−2a.c_n^2\;≤\; A^2\,e 2π n1-2a\;(2a )^n\;=\; A^2\,e 2π n1-2a\;α^n, α:=2a = 2a1-2a. (31) Since a<14a< 14, we have α<1α<1. Step 4 (Tail bound). For any η∈(α,1)η∈(α,1) the n n factor in (31) is absorbed into the geometric ratio: nαn≤Cηηn n\,α^n≤ C_η\,η^n, where Cη:=supn≥1n(α/η)n<∞C_η:= _n≥ 1 n\,(α/η)^n<∞ since α/η<1α/η<1. Substituting this into (31) and summing the resulting geometric series over n>Pn>P: ∑n=P+1∞cn2≤A2e2πCη1−2aηP+11−η. _n=P+1^∞c_n^2\;≤\; A^2\,e 2π\,C_η1-2a\; η^P+11-η. Summing over all coordinates (i,k)(i,k) and applying Parseval on the graph eigenbasis (∑j‖Y^n,j‖2=‖Yn‖ℓ2(V)2=∑i,kcn,ik2 _j\| Y_n,j\|^2=\|Y_n\|_ ^2(V)^2= _i,kc_n,ik^2) [29]: ∑n=P+1∞∑j=1N‖Y^n,j‖2≤CηP _n=P+1^∞ _j=1^N\| Y_n,j\|^2\;≤\;C\,η^P for a constant C depending on a, η, and ∑i,kAik2 _i,kA_ik^2. Now fix any β∈(α,1)β∈( α,1) and set η=β2η=β^2. Taking square roots gives (∑n>P∑j=1N‖Y^n,j‖2)1/2≤CtrβP,Ctr:=C. ( _n>P _j=1^N\| Y_n,j\|^2 )^1/2≤ C_ trβ^P, C_ tr:= C. Since the construction above matches all retained coefficients exactly, ‖Y−Y~P‖L2\|Y- Y_P\|_L^2 equals this omitted tail, which gives the bound stated in Theorem 2. Since β>αβ> α can be chosen arbitrarily close to α α, the effective decay rate is e−cPe^-cP with c=12ln(1/α)=12ln(1−2a2a)>0c= 12 (1/α)= 12 \! ( 1-2a2a )>0. ∎ A.7 Proof of Proposition 2 (Stochastic and mean-logit readouts) Lemma 3 (Readout regularity and quadrature conditions). Let zi(s)=∑n=0PZi,nΨn(ωs)z_i^(s)= _n=0^PZ_i,n _n( _s), where Zi,0Z_i,0 is the mean logit. Assume S≥P+1S≥ P+1, so the quadrature rule integrates ΨnΨm _n _m exactly for n,m≤Pn,m≤ P, and assume ℓ(⋅,yi) (·,y_i) is C3C^3 on a convex neighborhood of Zi,0Z_i,0 containing all quadrature logits. Then, with δi(s):=zi(s)−Zi,0 _i^(s):=z_i^(s)-Z_i,0, ∑sμsδi(s)=0,∑sμsδi(s)(δi(s))⊤=Σzi, _s _s\, _i^(s)=0, _s _s\, _i^(s)( _i^(s)) = _z_i, and Taylor expansion around Zi,0Z_i,0 has a remainder of order O(maxs‖δi(s)‖3)O( _s\| _i^(s)\|^3). For cross-entropy, ∇z2ℓ(z,y)=diag(p)−pp⊤ _z^2 (z,y)=diag(p)-p with p=softmax(z)p=softmax(z); for finite logits its nullspace is spanspan\1\. Proof. The first identity follows from ∑sμsΨn(ωs)=[Ψn]=0 _s _s _n( _s)=E[ _n]=0 for n≥1n≥ 1. The second follows from exact quadrature orthonormality, ∑sμsΨn(ωs)Ψm(ωs)=δnm _s _s _n( _s) _m( _s)= _nm, giving ∑sμsδi(s)(δi(s))⊤=∑n=1PZi,nZi,n⊤=Σzi _s _s _i^(s)( _i^(s)) = _n=1^PZ_i,nZ_i,n = _z_i. The C3C^3 condition gives the stated Taylor remainder. The cross-entropy Hessian is the covariance matrix of a categorical distribution with probabilities p, so its finite-logit nullspace consists exactly of common shifts. ∎ See 2 Proof. The quadrature logits at node ωs _s are zi(s)=∑n=0PZi,nΨn(ωs)z_i^(s)= _n=0^PZ_i,n _n( _s), and the quadrature-averaged loss is ℓ¯i=∑s=1Sμsℓ(zi(s),yi) _i= _s=1^S _s\, (z_i^(s),y_i). Step 1 (Taylor expansion at each quadrature node). Write zi(s)=Zi,0+δi(s)z_i^(s)=Z_i,0+ _i^(s) where δi(s)=∑n=1PZi,nΨn(ωs)∈ℝC _i^(s)= _n=1^PZ_i,n _n( _s) ^C is the stochastic perturbation. Expanding ℓ to second order around Zi,0Z_i,0, and using that ℓ(⋅,yi) (·,y_i) is C3C^3 in a neighborhood of Zi,0Z_i,0: ℓ(zi(s),yi)=ℓ(Zi,0,yi)+∇zℓ(Zi,0,yi)⊤δi(s)+12(δi(s))⊤∇z2ℓ(Zi,0,yi)δi(s)+O(‖δi(s)‖3). (z_i^(s),y_i)= (Z_i,0,y_i)+ _z (Z_i,0,y_i) _i^(s)+ 12( _i^(s)) _z^2 (Z_i,0,y_i)\, _i^(s)+O(\| _i^(s)\|^3). Step 2 (Quadrature averaging of the linear term). The first-order term vanishes under quadrature averaging: ∑s=1Sμs∇zℓ⊤δi(s)=∇zℓ⊤∑n=1PZi,n∑s=1SμsΨn(ωs)⏟=[Ψn]= 0=0, _s=1^S _s\, _z _i^(s)= _z _n=1^PZ_i,n _s=1^S _s\, _n( _s)_=\;E[ _n]\;=\;0=0, since [Ψn]=0E[ _n]=0 for n≥1n≥ 1 by the orthonormality of the Hermite basis (14), and the Gauss–Hermite rule integrates Ψn _n exactly for n≤2S−1n≤ 2S-1. Step 3 (Quadrature averaging of the quadratic term). For the second-order term, abbreviate ℋi:=∇z2ℓ(Zi,0,yi)H_i:= _z^2 (Z_i,0,y_i) and compute: ∑s=1Sμs(δi(s))⊤ℋiδi(s) _s=1^S _s\,( _i^(s)) H_i\, _i^(s) =∑s=1Sμs(∑n=1PZi,nΨn(ωs))⊤ℋi(∑m=1PZi,mΨm(ωs)) = _s=1^S _s ( _n=1^PZ_i,n _n( _s) )^\! H_i ( _m=1^PZ_i,m _m( _s) ) =∑n,m=1PZi,n⊤ℋiZi,m∑s=1SμsΨn(ωs)Ψm(ωs)⏟=δnm = _n,m=1^PZ_i,n H_i\,Z_i,m _s=1^S _s\, _n( _s) _m( _s)_=\; _nm =∑n=1PZi,n⊤ℋiZi,n=tr(ℋi∑n=1PZi,nZi,n⊤)=tr(ℋiΣzi), = _n=1^PZ_i,n H_i\,Z_i,n\;=\;tr\! (H_i _n=1^PZ_i,nZ_i,n )\;=\;tr(H_i\, _z_i), where the orthonormality ∑sμsΨn(ωs)Ψm(ωs)=δnm _s _s _n( _s) _m( _s)= _nm holds by exactness of the Gauss–Hermite rule for polynomials of degree n+m≤2P≤2S−1n+m≤ 2P≤ 2S-1, and Σzi=∑n=1PZi,nZi,n⊤ _z_i= _n=1^PZ_i,nZ_i,n is the logit covariance from (7). Step 4 (Combining). Substituting Steps 2 and 3 into the averaged Taylor expansion: ℓ¯i=∑sμsℓ(zi(s),yi)=ℓ(Zi,0,yi)+0+12tr(∇z2ℓ(Zi,0,yi)Σzi)+O(maxs‖δi(s)‖3). _i= _s _s\, (z_i^(s),y_i)= (Z_i,0,y_i)+0+ 12tr\! ( _z^2 (Z_i,0,y_i)\, _z_i )+O\! ( _s\| _i^(s)\|^3 ). (32) Step 5 (Positive semi-definiteness and nullspace characterization for cross-entropy). For the cross-entropy loss ℓ(z,y)=−logpy (z,y)=- p_y where p=softmax(z)p=softmax(z), the Hessian is ∇z2ℓ=diag(p)−pp⊤ _z^2 =diag(p)-p . This is the covariance matrix of a categorical distribution with probabilities p, hence positive semi-definite: for any v∈ℝCv ^C, v⊤(diag(p)−pp⊤)v=y^∼p[vy^2]−(y^∼p[vy^])2=Vary^∼p(vy^)≥0v (diag(p)-p )v=E_ y p[v_ y^2]-(E_ y p[v_ y])^2=Var_ y p(v_ y)≥ 0. Since ℋi⪰0H_i 0 and Σzi⪰0 _z_i 0, tr(ℋiΣzi)=tr(Σzi1/2ℋiΣzi1/2)≥0.tr(H_i _z_i)=tr\! ( _z_i^1/2H_i _z_i^1/2 )≥ 0. Moreover, tr(ℋiΣzi)=0tr(H_i _z_i)=0 if and only if Σzi1/2ℋiΣzi1/2=0 _z_i^1/2H_i _z_i^1/2=0, equivalently range(Σzi)⊆null(ℋi)range( _z_i) (H_i). For finite logits, softmax probabilities satisfy pc>0p_c>0 for every class, and the Hessian nullspace is null(ℋi)=spannull(H_i)=span\1\; thus purely common logit offsets do not contribute to the correction. The special case Σzi=0 _z_i=0 recovers the non-stochastic model, while in the fully confident limit p→eyp→ e_y we have ℋi→0H_i→ 0. Step 6 (Energy shift identity). For any scalar a, softmax(z+a)c=ezc+a∑r=1Cezr+a=ezc∑r=1Cezr=softmax(z)c.softmax(z+a1)_c= e^z_c+a _r=1^Ce^z_r+a= e^z_c _r=1^Ce^z_r=softmax(z)_c. Thus any uncertainty summary that is a function only of the softmax probabilities, such as predictive entropy, is invariant under common-logit shifts. By contrast, senergy(z+a)=−log∑c=1Cezc+a=−log(ea∑c=1Cezc)=senergy(z)−a.s_energy(z+a1)=- _c=1^Ce^z_c+a=- \! (e^a _c=1^Ce^z_c )=s_energy(z)-a. (33) For completeness, the analogous second-order expansion for the quadrature-averaged energy follows by the same calculation. Write zi(s)=Zi,0+δi(s)z_i^(s)=Z_i,0+ _i^(s) and A(z)=log∑cezcA(z)= _ce^z_c. The gradient and Hessian of A are ∇A(z)=p,∇2A(z)=diag(p)−pp⊤=:HA(z),∇ A(z)=p, ∇^2A(z)=diag(p)-p =:H_A(z), where p=softmax(z)p=softmax(z). Since senergy=−As_energy=-A, Taylor expansion around Zi,0Z_i,0 gives senergy(zi(s))=senergy(Zi,0)−pi⊤δi(s)−12(δi(s))⊤HA(Zi,0)δi(s)+O(‖δi(s)‖3).s_energy(z_i^(s))=s_energy(Z_i,0)-p_i _i^(s)- 12( _i^(s)) H_A(Z_i,0) _i^(s)+O(\| _i^(s)\|^3). The quadrature average of the linear term vanishes as in Step 2, and the quadratic term averages as in Step 3: ∑sμs(δi(s))⊤HA(Zi,0)δi(s)=tr(HA(Zi,0)Σzi), _s _s( _i^(s)) H_A(Z_i,0) _i^(s)=tr\! (H_A(Z_i,0) _z_i ), using exact quadrature orthonormality for pairwise products ΨnΨm _n _m with n,m≤Pn,m≤ P. Therefore ∑sμssenergy(zi(s))=senergy(Zi,0)−12tr(HA(Zi,0)Σzi)+O(maxs‖δi(s)‖3). _s _s\,s_energy(z_i^(s))=s_energy(Z_i,0)- 12tr\! (H_A(Z_i,0)\, _z_i )+O\! ( _s\| _i^(s)\|^3 ). (34) The Hessian HA(z)H_A(z) is the same categorical covariance matrix as the cross-entropy Hessian, so the stochastic energy correction also vanishes for common-logit stochastic variation. The non-stochastic common shift of the mean logit, however, is read out linearly by (33). ∎ Corollary 1 (Readout gap). Assume S≥P+1S≥ P+1 as in Lemma 3. Let f(z)=softmax(z)f(z)=softmax(z) and B:=supz∈ℝCmaxc‖∇z2fc(z)‖opB:= _z ^C _c\| _z^2f_c(z)\|_op, which is finite because every entry of ∇z2fc _z^2f_c is a polynomial in the probabilities f(z)∈[0,1]Cf(z)∈[0,1]^C; since ∇z2fc=fc[(ec−f)(ec−f)⊤−(diag(f)−ff⊤)] _z^2f_c=f_c [(e_c-f)(e_c-f) -(diag(f)-f ) ] is fcf_c times a difference of positive semidefinite matrices of operator norm at most 2(1−fc)22(1-f_c)^2 and 12 12, B≤12B≤ 12 for every C. Then for every node i, ‖p¯i−softmax(Zi,0)‖∞≤B2tr(Σzi),\| p_i-softmax(Z_i,0)\|_∞\;≤\; B2tr( _z_i), and argmaxcp¯i,c=argmaxcsoftmax(Zi,0)c _c p_i,c= _csoftmax(Z_i,0)_c whenever the gap between the two largest entries of softmax(Zi,0)softmax(Z_i,0) exceeds Btr(Σzi)Btr( _z_i). Proof. Write zi(s)=Zi,0+δi(s)z_i^(s)=Z_i,0+ _i^(s). Taylor’s theorem with integral remainder gives, for each class c, fc(zi(s))=fc(Zi,0)+∇fc(Zi,0)⊤δi(s)+∫01(1−t)(δi(s))⊤∇2fc(Zi,0+tδi(s))δi(s)dt,f_c(z_i^(s))=f_c(Z_i,0)+∇ f_c(Z_i,0) _i^(s)+ _0^1(1-t)\,( _i^(s)) ∇^2f_c (Z_i,0+t _i^(s) )\, _i^(s)\,dt, and the remainder is bounded in absolute value by B2‖δi(s)‖2 B2\| _i^(s)\|^2. Averaging with the quadrature weights, the linear term vanishes because ∑sμsδi(s)=0 _s _s _i^(s)=0 (Lemma 3), and ∑sμs‖δi(s)‖2=tr(Σzi) _s _s\| _i^(s)\|^2=tr( _z_i) by the same lemma, which gives the bound. Each component of p¯i p_i therefore lies within B2tr(Σzi) B2tr( _z_i) of the corresponding component of softmax(Zi,0)softmax(Z_i,0), so the two argmax rules agree whenever the top-two gap of softmax(Zi,0)softmax(Z_i,0) exceeds Btr(Σzi)Btr( _z_i). Since (11) penalizes ℰi=tr(Σzi)E_i=tr( _z_i) on training nodes, the bound is the structural reason the two readouts agree closely in practice (Section 4.2, Table 27). ∎ Appendix B Baseline Details B.1 Baseline methods Table 6 summarizes which tasks each method evaluates. G-Δ is the closest baseline in scope because it reports all three task families; we compare against it on calibration (Table 1) and distribution shift (Table 4). On OOD detection, the anchor-sampling mechanism in G-Δ causes out-of-memory errors in the GNNSafe OOD pipeline, preventing a direct comparison. Across the evaluated benchmarks, the DSS representation gives the strongest calibration and GOOD results among the compared methods, with the best AUROC on 7 of 11 OOD settings. The calibration and GOOD results use plain cross-entropy or ERM training without task-specific robust objectives; in the OOD detection experiments DSS-Hybrid is trained with the same energy-margin objective as GNNSafe++ (Appendix D). Table 6: Tasks evaluated by each method. ✓= evaluated in the original paper; (✓) = evaluated as a secondary metric or in the appendix; ✓∗ = evaluated in the original paper but OOM on our benchmark; ✗= not evaluated in the original paper or not designed for that task. Method Calibration OOD Detection Robust Classif. MC-dropout [10] (✓) ✗ ✗ Deep Ensembles [17] ✓ (✓) ✗ G-Δ [33] ✓ ✓∗ ✓ TFE-GNN [8] (✓) ✗ ✗ GNNSafe [36] ✗ ✓ ✗ GNNSafe++ [36] ✗ ✓ ✗ Graph-EBM [9] (✓) ✓ ✗ GPN [30] ✓ ✓ ✗ IRM [1] ✗ ✗ ✓ GroupDRO [27] ✗ ✗ ✓ TAR [38] ✗ ✗ ✓ DSS-GNN / DSS-Hybrid (ours) ✓ ✓ ✓ We compare against the following baselines across the three evaluation tasks. Calibration baselines. • TFE-GNN [8]: A spectral GNN that combines three complementary polynomial filters (low-pass, high-pass, and band-pass) via learned attention weights. It uses a propagate-first architecture where spectral filtering is applied once before an MLP backbone. • G-Δ [33]: An intrinsic uncertainty method that estimates predictive uncertainty through anchor distances in the learned embedding space, without requiring ensembles or Monte Carlo sampling. OOD detection baselines. • GNNSafe [36]: Uses the energy score E(x)=−Tlog∑cexp(zc/T)E(x)=-T _c (z_c/T) as an OOD detection signal, propagated over the graph adjacency to exploit neighborhood consistency. Does not require OOD exposure during training. • GNNSafe++ [36]: Extends GNNSafe with an energy-based regularization loss that uses OOD training data to separate in-distribution and OOD energy scores. Uses per-dataset margin hyperparameters (min,mout,λ)(m_in,m_out,λ). • Graph-EBM [9]: Models the joint distribution of node features and labels using an energy-based model with graph-aware components, including local and group energy terms propagated through graph diffusion. • MC-dropout [10] and Deep Ensembles [17]: Stochastic baselines on the same GCN backbone and training protocol, with M=20M=20 dropout passes at test time and M=5M=5 independently trained models respectively; both are scored with the same energy score and graph propagation as the energy-based methods (Appendix D). • GPN [30]: Graph Posterior Network, an evidential model that propagates per-node Dirichlet evidence over the graph; we report the values obtained on this benchmark by Wu et al. [36]. Distribution-shift baselines. Table 4 evaluates the GOOD/TAR baselines under concept shift and adds G-Δ . Results for ERM through TAR are taken from the TAR paper [38]. Representative baselines include: • ERM: Standard empirical risk minimization without any robustness modification. • IRM [1]: Invariant Risk Minimization, which penalizes features that are not invariant across training environments. • GroupDRO [27]: Distributionally robust optimization that upweights the worst-performing group during training. • VREx [16]: Variance Risk Extrapolation, which penalizes variance of risks across training environments. • TAR [38]: Topology-Aware Dynamic Reweighting, which uses graph-local transport to reweight training nodes based on loss and neighborhood structure. Appendix C Dataset Statistics and Licenses Tables 8–10 summarize the datasets used in our experiments. Table 7 lists the source repositories and licenses for all datasets. All datasets are publicly available and were used in accordance with their respective licenses. Table 7: Dataset sources and licenses. Dataset(s) Repository License Cora, Citeseer, PubMed https://github.com/kimiyoung/planetoid MIT Texas, Cornell, Wisconsin, https://github.com/bingzhewei/geom-gcn not stated Chameleon, Squirrel Coauthor-CS, Amazon-Photo https://github.com/shchur/gnn-benchmark MIT Roman-Empire, Amazon-Ratings, https://github.com/yandex-research/heterophilous-graphs MIT Minesweeper, Tolokers, Questions Twitch https://github.com/benedekrozemberczki/MUSAE GPL-3.0 ogbn-arxiv https://github.com/snap-stanford/ogb MIT GOOD benchmark https://github.com/divelab/GOOD GPL-3.0 C.1 Calibration benchmarks Table 8: Node classification datasets for calibration experiments (Section 4.2). Results use 10 splits per dataset (Appendix D); edges are counted in both directions. Cornell uses the corrected files of the geom-gcn repository, whose initial release duplicated the Texas features and labels. Datasets below the mid-rule are from the heterophilous suite of Platonov et al. [26]. Dataset Nodes Edges Classes Features Cora 2,708 10,556 7 1,433 Citeseer 3,327 9,228 6 3,703 PubMed 19,717 88,651 3 500 Texas 183 574 5 1,703 Cornell 183 557 5 1,703 Wisconsin 251 916 5 1,703 Chameleon 2,277 62,792 5 2,325 Squirrel 5,201 396,846 5 2,089 CS 18,333 163,788 15 6,805 Roman-Empire 22,662 65,854 18 300 Amazon-Ratings 24,492 186,100 5 300 Minesweeper 10,000 78,804 2 7 Tolokers 11,758 1,038,000 2 10 Questions 48,921 307,080 2 301 C.2 OOD detection benchmarks Table 9: OOD detection datasets (Section 4.3). For node-OOD (Cora, Amazon-Photo, Coauthor-CS), OOD nodes are generated via structure (stochastic block model), feature (random interpolation of node features), or label (held-out classes) perturbations. For cross-graph OOD, the ID and OOD graphs are from different domains. Statistics are for the in-distribution graph. Dataset OOD type Nodes Edges Classes Features Cora structure / feature / label 2,708 10,556 7 1,433 Amazon-Photo structure / feature / label 7,650 238,162 8 745 Coauthor-CS structure / feature / label 18,333 163,788 15 6,805 Twitch (DE→ /FR/RU) cross-graph 9,498 306,276 2 3,170 Arxiv (≤ 2015→ 2018–20) cross-graph (temporal) 169,343 2,315,598 40 128 C.3 GOOD distribution-shift benchmarks Table 10 reports the node-level concept-shift settings from GOOD [12] used in our distribution-shift experiments. Table 10: GOOD benchmark settings (Section 4.4). All use concept shift. The domain column indicates the covariate that defines the distribution shift. Node counts are for the full graph before splitting; GOOD-Cora is built on the full Cora citation graph (CoraFull), not on the 7-class Planetoid subset of Table 8. Setting Domain Nodes Edges Classes Features GOOD-Cora / degree node degree 19,793 126,842 70 8,710 GOOD-Cora / word word diversity 19,793 126,842 70 8,710 GOOD-Arxiv / degree node degree 169,343 2,315,598 40 128 GOOD-Arxiv / time publication year 169,343 2,315,598 40 128 GOOD-CBAS / color node color 700 3,962 4 4 GOOD-WebKB / university university 617 1,138 5 1,703 GOOD-Twitch / language user language 34,120 892,346 2 128 Appendix D Implementation Details D.1 General training protocol All DSS-GNN models are trained with the mean-logit cross-entropy loss ℒbaseL_base (11); the quadrature-averaged loss ℒquad=1|tr|∑i∈trℓ¯i+λreg1|tr|∑i∈trℰi,ℓ¯i=∑s=1Sμsℓ(zi(s),yi),L_quad= 1|V_tr| _i _tr _i\;+\; _reg 1|V_tr| _i _trE_i, _i= _s=1^S _s\, (z_i^(s),y_i), (35) supervises every quadrature node, agrees with ℒbaseL_base to first order by Proposition 2, and is used only in that analysis. The loss includes the chaos energy regularizer; λreg _reg is selected per dataset on validation Brier from 0,10−3,10−2,10−1\0,10^-3,10^-2,10^-1\ (Table 12; Table 17 reports its effect at P=2P=2), with 0.010.01 as the default. Training uses early stopping on validation loss with the patience listed per task below. All calibration results report mean ± standard deviation over 10 splits, identical for every method. For the nine datasets above the mid-rule of Table 8 the splits are random with fixed seeds: a class-balanced training set whose per-class quota totals 60% of the nodes when every class is large enough, 20% of the nodes for validation, and the remainder for testing. The five heterophilous-suite datasets use the first 10 public 50/25/25 splits of Platonov et al. [26]. D.2 Calibration (Table 1) Shared defaults. Chebyshev filter degrees Klp=Khp=4K_lp=K_hp=4, chaos order centered at P=2P=2 in the base grid, quadrature nodes S=4S=4, λreg=0.01 _reg=0.01, deterministic scalar gates (Pg=0P_g=0), loss type = mean, activation = ReLU, learning rate 0.010.01, weight decay 5×10−45×10^-4, propagation dropout 0.50.5, linear dropout 0.00.0, epochs 10001000, patience 200200. The best architecture variant per dataset was selected from a grid over propagation type (Chebyshev / TFE-adjacency / propagate-first), graph filter symmetry (symmetric / random-walk), optimizer (Adam / RMSprop), hidden dimension (6464 / 128128), number of layers (22 / 33), and dropout (0.50.5 / 0.60.6 / 0.80.8), evaluated on validation Brier score. Roman-Empire uses hidden 32 and Klp=Khp=2K_lp=K_hp=2, outside this grid, Minesweeper uses learning rate 0.0050.005, and Chameleon uses sum aggregation. The chaos order P is a regular hyperparameter selected per dataset on validation Brier over the sweep of Appendix E; Table 12 lists the selected order per dataset together with that single configuration’s accuracy and Brier. Per-dataset configurations. Table 11 lists the architecture variant and any hyperparameter overrides from the shared defaults for each dataset. Table 11: Per-dataset non-P DSS-GNN configuration for the calibration experiments (Table 1). “Prop-first” = propagate-once structure (spectral filtering followed by MLP); “TFE-adj” = TFE-style adjacency propagation; “Cheb” = Chebyshev filters on rescaled Laplacian. Unlisted hyperparameters use the shared defaults above; P, λreg _reg, and S per dataset are listed in Table 12. Dataset Architecture Hidden L Optimizer Overrides Cora Prop-first (rw) 64 2 Adam dropout 0.8 Citeseer Cheb 64 2 Adam none PubMed Cheb 128 2 Adam none Texas Prop-first 64 2 RMSprop none Cornell Cheb 128 2 Adam dropout 0.8 Wisconsin Cheb 64 2 Adam dropout 0.6 Chameleon Prop-first (sum) 64 2 Adam none Squirrel Prop-first (rw) 64 2 RMSprop none CS TFE-adj 64 2 Adam none Roman-Emp. TFE-adj 32 3 Adam hidden 32, K=2K=2 Amz-Rat. TFE-adj 64 3 Adam none Minesweeper Cheb 64 2 Adam lr 0.005 Tolokers Prop-first (rw) 64 2 RMSprop none Questions Prop-first 64 2 RMSprop none Selected chaos order per dataset. Table 12 lists, for each dataset, the chaos order and regularization strength selected on validation Brier, with the accuracy and Brier of that configuration over the same 10 splits, taken from the sweeps of Tables 13 and 14, and its validation Brier (recorded in a separate run of the same configuration on the same splits); Table 1 reports these configurations. DSS-GNN has the best mean Brier on all 14 datasets, including the closest comparisons (Texas 0.2400.240 vs 0.2550.255, Squirrel 0.3940.394 vs 0.3950.395, Amazon-Ratings 0.6340.634 vs 0.6370.637, all TFE-GNN; Questions 0.0530.053 vs 0.0540.054, G-Δ ). At least one of P=1P=1 or P=2P=2 is within 0.020.02 Brier of the best swept order on every dataset (Section 4.2), so the larger selected orders are near-ties rather than requirements. The quadrature size is S=4S=4 throughout: exact for the polynomial filter and gate terms at P≤3P≤ 3 (Appendix A.1), while at P≥4P≥ 4 the highest-order Hermite products are integrated approximately. The recommended default P=1P=1–22 is exact for these terms at S=4S=4, and the S sweep at P=2P=2 (Table 18) moves Brier by at most 0.0110.011 on 12 of 14 datasets. Table 12: Selected configuration per dataset for the calibration experiments (Table 1): chaos order P, chaos energy regularization λreg _reg, quadrature size S, and that single configuration’s test accuracy (%) and Brier score (mean ± SD over the same 10 splits). The last column reports the validation Brier of the selected configuration, recorded in a separate run on the same 10 splits (for Cornell, in the same run). Dataset P λreg _reg S Accuracy Brier Val. Brier Cora 1 0.01 4 86.63 ± 1.26 0.207 ± 0.019 0.203 Citeseer 2 0 4 80.20 ± 1.28 0.308 ± 0.010 0.315 PubMed 5 0.01 4 89.72 ± 0.31 0.156 ± 0.005 0.152 Texas 2 0.01 4 91.31 ± 3.44 0.240 ± 0.111 0.190 Cornell 1 0.01 4 85.11 ± 5.12 0.231 ± 0.079 0.166 Wisconsin 2 0 4 93.25 ± 3.39 0.104 ± 0.047 0.088 Chameleon 1 0.01 4 75.36 ± 1.19 0.339 ± 0.019 0.333 Squirrel 1 0.1 4 68.06 ± 2.21 0.394 ± 0.027 0.415 CS 1 0.01 4 96.17 ± 0.23 0.062 ± 0.004 0.060 Roman-Empire 1 0 4 78.84 ± 0.61 0.300 ± 0.007 0.304 Amazon-Ratings 8 0.01 4 49.72 ± 0.64 0.634 ± 0.005 0.637 Minesweeper 6 0.001 4 87.66 ± 1.10 0.169 ± 0.014 0.170 Tolokers 2 0.01 4 79.88 ± 0.49 0.268 ± 0.005 0.305 Questions 4 0.001 4 97.17 ± 0.04 0.053 ± 0.001 0.054 D.3 OOD detection (Section 4.3) DSS-Hybrid uses a 2-layer GCN base encoder (hidden 64, dropout 0.0) with a DSS residual branch (Klp=3K_lp=3, Khp=2K_hp=2, P=1P=1, S=4S=4, deterministic scalar gates). The base encoder is warmed up for 50 epochs before the DSS branch is activated; the residual scale γres _res is initialized to 0.10.1. All methods run in our pipeline use Adam with learning rate 0.010.01, weight decay 0.010.01, 200 epochs, and 3 seeds, following the GNNSafe protocol. The GNNSafe pipeline uses BatchNorm in its GCN backbones, and so does the DSS-Hybrid base; the standalone model of Appendix G adds per-chaos-channel BatchNorm in all OOD and GOOD cells but OOD Twitch and GOOD-Arxiv/degree. OOD detection scores are computed with the energy score (9) at temperature T=1T=1, propagated over the graph adjacency with K=2K=2 hops and mixing α=0.5α=0.5, following the GNNSafe protocol [36]. DSS-Hybrid is trained with the GNNSafe++ energy-margin regularizer, using the per-dataset margin hyperparameters (min,mout,λ)(m_in,m_out,λ) of Wu et al. [36], as GNNSafe++ is; Graph-EBM, MC-dropout, and Deep Ensembles are trained without it, as GNNSafe is. MC-dropout keeps dropout active at test time with BatchNorm frozen and aggregates M=20M=20 stochastic passes; Deep Ensembles aggregate M=5M=5 independently initialized backbones; both are scored with the same energy score and propagation as the energy-based methods, applied to the mean logits across passes or members. Standalone model in the OOD and GOOD pipelines. For the cross-evaluation of Appendix G, the standalone model adds one BatchNorm1d per chaos channel between hidden DSS layers where noted; Table 1 uses no BatchNorm. The per-cell configurations are listed in Appendix G. Graph-EBM reproduction. Graph-EBM is run with a port of the official GraphEBMWrapper [9] at the official default hyperparameters: the multi-scale corrected-conditional energy model with Mahalanobis logit correction, three energy scales, and label-propagation diffusion (diagonal untied covariance, γ=1γ=1, unit energy weights, α=0.5α=0.5 with k=10k=10 steps). The scorer is fitted on the same GCN encoder (with BatchNorm, as in the GNNSafe pipeline), trained with the same hyperparameters and seeds as the other baselines run in our pipeline; Graph-EBM’s own diffusion replaces the GNNSafe score propagation for this row. On an exported Cora/structure run (trained GCN logits and embeddings, both graphs, labels, and masks) the port and the official implementation agree at Spearman rank correlation 0.9940.994, and exactly (rank correlation 1.01.0) once the interval and covariance estimators of the official implementation are used. D.4 Distribution shift (Table 4) DSS-Hybrid uses a 2-layer GCN base encoder with a DSS residual branch (Klp=3K_lp=3, Khp=2K_hp=2, P=1P=1, S=4S=4, hidden 64, deterministic scalar gates, 3 layers total). Training uses Adam with learning rate 10−310^-3, weight decay 5×10−45×10^-4, dropout 0.50.5, 500 epochs, patience 200. The base encoder is warmed up for 25 epochs. Standard ERM training is used (no robust objective). Published baseline results (ERM through TAR) are taken from the TAR paper [38]; G-Δ is run in our pipeline. D.5 Software, hardware, and compute budget Experiments use PyTorch and PyTorch Geometric (the released environment pins PyTorch 2.5) on NVIDIA A100 (40 GB) and H100 GPUs; the ogbn-arxiv timing runs of Table 20 use one H100, and the Cornell experiments run on CPU. Mixed precision (FP16) is used for calibration experiments to reduce memory on larger graphs. On the A100 runs, per-epoch wall-clock times range from ∼35 35 ms (Cora, P=0P=0) to ∼650 650 ms (Amazon-Ratings, P=3P=3); see Table 20 for a detailed breakdown. Compute estimate. The reported experiments (calibration: 14 datasets × 10 runs; OOD: 5 datasets × 3 seeds; GOOD: 7 settings) total approximately 8 A100 GPU-hours. The ablation studies (chaos order P, regularization λreg _reg, quadrature nodes S, and contribution decomposition; Appendix E) add approximately 100 GPU-hours. Hyperparameter selection (grid search over architecture variant, hidden dimension, optimizer, and dropout) required an additional ∼25 25 GPU-hours. The total compute for the project, including preliminary experiments not reported in the paper, is approximately 200 A100 GPU-hours. Code is available at https://github.com/heraclixus/DSSGNN. Appendix E Ablation Study Details E.1 Polynomial order P across all datasets Figure 3 summarizes the chaos-order sweep referenced in Section 4.2, and Tables 13 and 14 report the per-dataset classification accuracy and Brier score for DSS-GNN across all 14 datasets for P∈0,1,2,3,4,5,6,8P∈\0,1,2,3,4,5,6,8\. All results are mean ± standard deviation over 10 random splits. For each dataset, the best-performing P is highlighted in bold. Figure 3: (a) Percentage Brier reduction of the best swept order P>0P>0 over P=0P=0 (Table 14). (b) Brier relative to P=0P=0 versus chaos order, for seven datasets. Table 13: Accuracy (%) of DSS-GNN for varying chaos order P across all datasets. Dataset P=0P=0 P=1P=1 P=2P=2 P=3P=3 P=4P=4 P=5P=5 P=6P=6 P=8P=8 Cora 85.12 ± 1.06 86.63 ± 1.26 86.33 ± 1.12 86.33 ± 0.82 86.63 ± 1.14 86.32 ± 1.13 86.61 ± 0.84 86.38 ± 0.96 Citeseer 65.39 ± 1.02 79.21 ± 0.79 80.20 ± 1.28 78.47 ± 0.80 71.50 ± 4.58 78.92 ± 1.13 76.77 ± 2.19 78.49 ± 1.08 PubMed 89.12 ± 0.57 89.33 ± 0.49 89.33 ± 0.66 89.38 ± 0.67 89.33 ± 0.54 89.72 ± 0.31 89.40 ± 0.57 89.39 ± 0.51 Texas 75.08 ± 23.70 88.36 ± 4.57 91.31 ± 3.44 89.18 ± 4.92 91.31 ± 3.30 87.54 ± 5.59 89.01 ± 2.82 87.37 ± 5.43 Cornell 83.62 ± 6.80 85.11 ± 5.12 85.74 ± 6.73 83.62 ± 5.39 84.26 ± 6.87 83.83 ± 7.38 83.62 ± 6.53 82.98 ± 7.61 Wisconsin 93.25 ± 3.46 92.87 ± 2.92 93.25 ± 3.39 88.50 ± 3.35 93.00 ± 3.46 90.75 ± 2.84 89.63 ± 4.82 90.50 ± 2.92 Chameleon 74.53 ± 2.32 75.36 ± 1.19 74.53 ± 1.74 74.51 ± 1.57 74.77 ± 1.86 74.92 ± 1.60 74.40 ± 1.65 73.63 ± 1.69 Squirrel 66.35 ± 2.56 68.06 ± 2.21 67.26 ± 2.13 66.99 ± 1.24 67.39 ± 2.44 67.23 ± 1.55 67.02 ± 1.50 67.71 ± 1.51 CS 95.94 ± 0.20 96.17 ± 0.23 96.01 ± 0.25 96.03 ± 0.19 96.13 ± 0.25 95.96 ± 0.29 96.02 ± 0.23 96.02 ± 0.17 Roman-Emp. 79.01 ± 0.27 78.84 ± 0.61 78.81 ± 0.61 78.98 ± 0.87 78.75 ± 0.32 78.77 ± 0.69 79.01 ± 0.27 78.97 ± 0.79 Amazon-Rat. 45.40 ± 0.76 48.13 ± 1.14 49.41 ± 0.71 49.65 ± 1.15 49.36 ± 0.70 49.05 ± 0.73 49.35 ± 0.83 49.72 ± 0.64 Minesweeper 86.14 ± 0.60 86.93 ± 0.51 87.37 ± 1.10 86.37 ± 1.06 86.53 ± 0.96 86.17 ± 0.45 87.66 ± 1.10 87.42 ± 0.93 Tolokers 79.33 ± 0.66 79.82 ± 0.45 79.88 ± 0.49 79.28 ± 0.79 79.62 ± 0.56 79.32 ± 0.66 79.21 ± 0.75 79.36 ± 0.79 Questions 97.15 ± 0.06 97.15 ± 0.06 97.13 ± 0.06 97.16 ± 0.06 97.17 ± 0.04 97.15 ± 0.06 97.15 ± 0.06 97.14 ± 0.09 Table 14: Brier Score (↓ ) of DSS-GNN for varying chaos order P across all datasets. Dataset P=0P=0 P=1P=1 P=2P=2 P=3P=3 P=4P=4 P=5P=5 P=6P=6 P=8P=8 Cora 0.228 ± 0.013 0.207 ± 0.019 0.212 ± 0.016 0.209 ± 0.016 0.208 ± 0.019 0.215 ± 0.015 0.211 ± 0.011 0.211 ± 0.015 Citeseer 0.427 ± 0.008 0.346 ± 0.008 0.308 ± 0.010 0.326 ± 0.012 0.422 ± 0.027 0.327 ± 0.012 0.370 ± 0.023 0.331 ± 0.013 PubMed 0.162 ± 0.007 0.158 ± 0.007 0.160 ± 0.008 0.159 ± 0.008 0.161 ± 0.007 0.156 ± 0.005 0.159 ± 0.008 0.159 ± 0.007 Texas 0.269 ± 0.082 0.213 ± 0.111 0.240 ± 0.111 0.233 ± 0.096 0.253 ± 0.096 0.253 ± 0.086 0.249 ± 0.114 0.242 ± 0.082 Cornell 0.241 ± 0.087 0.231 ± 0.079 0.228 ± 0.088 0.250 ± 0.090 0.263 ± 0.129 0.247 ± 0.107 0.252 ± 0.079 0.273 ± 0.111 Wisconsin 0.103 ± 0.048 0.107 ± 0.041 0.104 ± 0.047 0.161 ± 0.060 0.103 ± 0.048 0.129 ± 0.037 0.153 ± 0.052 0.135 ± 0.040 Chameleon 0.362 ± 0.019 0.339 ± 0.019 0.347 ± 0.022 0.350 ± 0.016 0.352 ± 0.016 0.351 ± 0.018 0.355 ± 0.013 0.360 ± 0.015 Squirrel 0.416 ± 0.028 0.394 ± 0.027 0.405 ± 0.024 0.411 ± 0.019 0.398 ± 0.028 0.400 ± 0.016 0.406 ± 0.014 0.393 ± 0.017 CS 0.065 ± 0.003 0.062 ± 0.004 0.063 ± 0.003 0.064 ± 0.003 0.063 ± 0.003 0.064 ± 0.004 0.063 ± 0.003 0.063 ± 0.003 Roman-Emp. 0.298 ± 0.005 0.300 ± 0.007 0.300 ± 0.008 0.298 ± 0.005 0.299 ± 0.005 0.300 ± 0.009 0.299 ± 0.007 0.298 ± 0.005 Amazon-Rat. 0.666 ± 0.003 0.642 ± 0.010 0.635 ± 0.006 0.633 ± 0.006 0.636 ± 0.006 0.637 ± 0.005 0.636 ± 0.004 0.634 ± 0.005 Minesweeper 0.190 ± 0.007 0.180 ± 0.006 0.170 ± 0.014 0.185 ± 0.011 0.181 ± 0.011 0.187 ± 0.005 0.169 ± 0.014 0.171 ± 0.012 Tolokers 0.278 ± 0.010 0.268 ± 0.004 0.268 ± 0.005 0.276 ± 0.012 0.280 ± 0.009 0.277 ± 0.010 0.278 ± 0.009 0.279 ± 0.009 Questions 0.055 ± 0.001 0.053 ± 0.001 0.054 ± 0.002 0.053 ± 0.001 0.053 ± 0.001 0.053 ± 0.001 0.053 ± 0.001 0.053 ± 0.002 On four datasets the best swept order exceeds two (PubMed, Squirrel, Amazon-Ratings, and Minesweeper); on every dataset P=1P=1 or P=2P=2 is within 0.02 Brier of the best swept value. E.2 Chaos order on OOD detection and under distribution shift OOD detection. Varying only the chaos order P∈0,1,2,3P∈\0,1,2,3\ in the DSS-Hybrid OOD configuration (GCN base, energy score with propagation, energy-margin training, every other hyperparameter fixed; 3 runs) leaves energy AUROC flat: Cora-structure 91.0/91.2/90.9/90.891.0/91.2/90.9/90.8 and Amazon-Photo/feature 99.4/99.5/99.5/99.599.4/99.5/99.5/99.5 for P=0,1,2,3P=0,1,2,3, Coauthor-CS likewise flat, and on ogbn-arxiv the AUROC moves by at most 0.80.8 across P=0P=0 to 33 (the runs of Table 20). The absolute levels differ from Table 2 because the ablation is a separate set of runs at one configuration per dataset. This is the behavior the design predicts: the energy score reads the mean logit alone (9), so P can be selected for calibration without affecting detection, and the detection margin comes from the DSS representation under a fixed readout (Section 4.3). Distribution shift. At the validation-selected standalone configurations of Appendix G, Table 15 varies only P (3 seeds). P=1P=1 is the argmax on both the OOD-validation and the test split on all three settings, so at these configurations the chaos order contributes +1.7+1.7, +1.2+1.2, and +0.6+0.6 points of shifted accuracy over the deterministic P=0P=0 arm. Table 15: Shifted test accuracy (%) of standalone DSS-GNN on GOOD concept shift as the chaos order varies at the validation-selected configuration (mean ± SD over 3 seeds). Setting P=0P=0 P=1P=1 P=2P=2 GOOD-Cora / word 63.24 ± 0.40 64.89 ± 0.21 64.01 ± 0.25 GOOD-Cora / degree 61.38 ± 0.81 62.60 ± 0.59 61.22 ± 0.07 GOOD-Arxiv / time 64.81 ± 0.27 65.43 ± 1.07 64.48 ± 0.27 E.3 Energy-score propagation Table 16 isolates the score-propagation step of the energy readout on DSS-Hybrid: the raw energy and the propagated energy are computed from the same trained checkpoints (one fixed configuration, 3 runs), so the difference is due to propagation alone. Propagation helps on 9 of 11 settings, by up to +18.2+18.2 AUROC (Twitch) and +14.5+14.5 (Cora-structure); the two exceptions are Cora/label (−6.7-6.7), where the fixed smoothing averages away part of the leave-out signal inside homophilous neighborhoods, and Amazon-Photo/structure (−0.4-0.4), at noise level. The same step transfers to the standalone model (Cora/label +0.9+0.9, Amazon-Photo/label +13.4+13.4 at identical checkpoints). We apply propagation uniformly for protocol parity with GNNSafe++. Table 16: Effect of energy-score propagation on DSS-Hybrid: change in AUROC (points) from the raw to the propagated energy score at identical checkpoints (3 runs, fixed configuration). Dataset Structure Feature Label Cora +14.5+14.5 +5.4+5.4 −6.7-6.7 Amazon-Photo −0.4-0.4 +1.7+1.7 +3.9+3.9 Coauthor-CS +1.9+1.9 +1.1+1.1 +1.3+1.3 Twitch (cross-graph) +18.2+18.2 Arxiv (cross-graph) +7.3+7.3 E.4 Chaos energy regularization λreg _reg Table 17 reports the Brier score as the chaos energy regularization strength λreg _reg varies from 0 to 0.1 with P=2P=2 fixed. Each column is an independent set of 10-split runs; at each dataset’s selected λreg _reg the cell agrees with the P=2P=2 column of Table 14 within one standard deviation. The regularizer of (11) penalizes large chaos coefficients; at λreg=0 _reg=0 the chaos expansion is still active (all P+1P+1 channels are learned), but without an explicit energy penalty. Moderate regularization (λreg∈10−3,10−2 _reg∈\10^-3,10^-2\) yields the best calibration on the majority of datasets. On Texas, λreg _reg matters most (0.3660.366 at 00 versus 0.2130.213 at 10−210^-2); on Citeseer, Wisconsin, and Roman-Empire λreg=0 _reg=0 is best and larger values degrade Brier, which is why λreg _reg is selected per dataset on validation Brier. Table 17: Effect of chaos energy regularization λreg _reg on Brier score (↓ ), with P=2P=2 fixed. Bold indicates the best λreg _reg per dataset. Dataset λreg=0 _reg=0 10−310^-3 10−210^-2 10−110^-1 Cora 0.219 0.219 0.207 0.220 Citeseer 0.308 0.310 0.318 0.499 PubMed 0.159 0.156 0.156 0.157 Texas 0.366 0.246 0.213 0.251 Cornell 0.230 0.237 0.228 0.342 Wisconsin 0.103 0.107 0.108 0.174 Chameleon 0.362 0.354 0.339 0.342 Squirrel 0.405 0.407 0.398 0.393 CS 0.064 0.066 0.062 0.083 Roman-Emp. 0.298 0.323 0.603 0.789 Amazon-Rat. 0.663 0.660 0.633 0.633 Minesweeper 0.186 0.169 0.190 0.228 Tolokers 0.281 0.275 0.268 0.280 Questions 0.054 0.053 0.054 0.055 E.5 Number of quadrature nodes S Table 18 reports the Brier score as the number of Gauss–Hermite quadrature nodes S varies from 2 to 8, with P=2P=2 fixed. Each column is an independent set of 10-split runs; the S=4S=4 column agrees with the P=2P=2 column of Table 14 within one standard deviation. For the deterministic gates used in these runs (Pg=0P_g=0), the theoretical minimum for exact chaos arithmetic at P=2P=2 is S≥⌈(Pg+2P+1)/2⌉=3S≥ (P_g+2P+1)/2 =3 (Appendix A.1). Brier varies by at most 0.0110.011 across S∈2,…,8S∈\2,…,8\ on 12 of 14 datasets. The exceptions are Wisconsin, where S=2S=2, below the theoretical minimum, gives a higher Brier score (0.1420.142 vs 0.1090.109 at S=3S=3), consistent with the exactness requirement of Theorem 3, and Cornell, where the spread across S (0.2180.218 to 0.2860.286) is within one standard deviation of its 10-split results. The chaos-order sweeps of Tables 13–14 use S=4S=4 throughout, exact for P≤3P≤ 3; at P≥4P≥ 4 the highest-order Hermite products are integrated approximately (Appendix D). Table 18: Brier score (↓ ) vs. number of quadrature nodes S, with P=2P=2 fixed. The theoretical minimum is S=3S=3. Dataset S=2S=2 S=3S=3 S=4S=4 S=6S=6 S=8S=8 Cora 0.216 0.207 0.217 0.209 0.212 Citeseer 0.310 0.308 0.312 0.313 0.311 PubMed 0.157 0.156 0.156 0.158 0.156 Texas 0.216 0.221 0.213 0.217 0.218 Cornell 0.235 0.286 0.228 0.218 0.242 Wisconsin 0.142 0.109 0.115 0.103 0.108 Chameleon 0.341 0.344 0.343 0.339 0.340 Squirrel 0.395 0.397 0.393 0.396 0.400 CS 0.062 0.062 0.062 0.063 0.063 Roman-Emp. 0.298 0.298 0.298 0.298 0.298 Amazon-Rat. 0.633 0.634 0.636 0.636 0.637 Minesweeper 0.170 0.169 0.170 0.171 0.170 Tolokers 0.269 0.269 0.273 0.270 0.268 Questions 0.055 0.053 0.054 0.055 0.055 E.6 Contribution decomposition: backbone vs. chaos Table 19 decomposes the calibration improvement into two factors: the graph-spectral backbone and the chaos expansion. The “ChebNet” column reports a vanilla single-branch Chebyshev filter (K=4K=4, hidden 64) as a diagnostic backbone baseline, not as part of the main uncertainty-aware baseline comparison in Table 1. The “P=0P=0” column uses the Table 1 architecture for each dataset (which may include propagate-first or dual branches) with the chaos expansion disabled. The “Best P>0P>0” column adds the chaos expansion on the same architecture. The P=0→P>0P=0→ P>0 comparison is strictly fair (same architecture, only P changes) and shows that the chaos expansion improves calibration on 12 of 14 datasets (3.63.6–27.9%27.9\% Brier reduction), indicating that the stochastic expansion is the primary source of the calibration gains over the P=0P=0 version in this controlled comparison. The ChebNet column gives useful context: on the small WebKB datasets Texas and Wisconsin, the single-filter ChebNet attains a lower absolute Brier score, while the controlled P=0P=0 versus P>0P>0 comparison still isolates the contribution of the chaos expansion. The ChebNet column is also the single-filter control for the filtering design: at P=0P=0 the per-dataset dual-filter architecture has lower or equal Brier than the single-branch ChebNet on 12 of 14 datasets (equal on Questions; Chameleon 0.3620.362 vs 0.7760.776, Squirrel 0.4160.416 vs 0.7980.798, Roman-Empire 0.2980.298 vs 0.3680.368), with Texas and Wisconsin as the exceptions noted above. Table 19: Contribution decomposition: Brier score (↓ ). ChebNet K=4K=4 is a vanilla single-filter diagnostic; P=0P=0 and Best P>0P>0 are the P=0P=0 column and the best P>0P>0 column of Table 14 (same per-dataset architecture). Boldface compares the controlled P=0P=0 versus Best P>0P>0 pair, and the rightmost column isolates the chaos contribution. Table 1 reports the selected order of Table 12, which can differ from the best swept order. Dataset ChebNet K=4K=4 P=0P=0 Best P>0P>0 Δ Chaos Cora 0.281 0.228 0.207 −9.2%-9.2\% Citeseer 0.470 0.427 0.308 −27.9%-27.9\% PubMed 0.171 0.162 0.156 −3.7%-3.7\% Texas 0.208 0.269 0.213 −20.8%-20.8\% Cornell 0.324 0.241 0.228 −5.4%-5.4\% Wisconsin 0.073 0.103 0.103 0% -0\% Chameleon 0.776 0.362 0.339 −6.4%-6.4\% Squirrel 0.798 0.416 0.393 −5.5%-5.5\% CS 0.069 0.065 0.062 −4.6%-4.6\% Roman-Emp. 0.368 0.298 0.298 0% -0\% Amazon-Rat. 0.676 0.666 0.633 −5.0%-5.0\% Minesweeper 0.208 0.190 0.169 −11.1%-11.1\% Tolokers 0.303 0.278 0.268 −3.6%-3.6\% Questions 0.055 0.055 0.053 −3.6%-3.6\% E.7 Computational overhead Asymptotic complexity. Although the doubly-spectral expansion is defined on the tensor product ℓ2(V)⊗L2(ℝ,γ) ^2(V) L^2(R,γ), truncating the chaos basis at order P retains only P+1P+1 coefficients along the stochastic axis, so the graph-filtering cost and the coefficient memory scale as (P+1)×(P+1)× the deterministic ChebNet baseline and the node-wise dense cost as S×S×. Concretely, one DSS layer takes ((P+1)(Klp+Khp)|E|d+SNd2+S(P+1)Nd)O\! ((P+1)(K_lp+K_hp)|E|d+SNd^2+S(P+1)Nd ) time and ((P+1+S)Nd)O((P+1+S)Nd) memory, where |E||E| is the number of edges, d the hidden width, and S the number of Gauss–Hermite quadrature nodes; the chaos extension is therefore additive, not combinatorial, in the stochastic direction, and no eigendecomposition is ever computed. With S at the exactness threshold (S=P+1S=P+1) the quadrature term is quadratic in P; this is the worst-case operation count, and the measured growth in P is much smaller (below). By Theorem 2 the truncation error decays exponentially in P for targets satisfying the growth condition of Appendix A.6, consistent with the empirical plateau at P=1P=1–22, so the constant factor over a deterministic backbone stays small. All chaos coefficients share the same Chebyshev filters and projection weights, so the parameter count per layer is (Klp+Khp+(P+1)(Pg+1)+dℓdℓ+1)O (K_lp+K_hp+(P+1)(P_g+1)+d_ d_ +1 ) (filter coefficients, gate coefficients an,r,bn,ra_n,r,b_n,r, and the linear map); the chaos order enters only through the gate coefficients and the input lift. By contrast, the MC-dropout and Deep Ensemble baselines require M forward passes or models at inference (M=20M=20 and M=5M=5 in our comparisons). Wall-clock measurements. Table 20 reports the wall-clock time per training epoch as the chaos order P increases, on three datasets spanning two orders of magnitude in size. The overhead grows linearly with P: each additional chaos order adds one extra graph filter evaluation and Gauss–Hermite projection per layer. At P=2P=2 the per-epoch cost ranges from 2.9×2.9× on Cora to 1.2×1.2× on Amazon-Ratings and 1.1×1.1× on ogbn-arxiv relative to the deterministic baseline (P=0P=0), within the asymptotic bound above and well below the M×M× overhead of MC-dropout or ensembles. The measured overhead decreases as graphs grow because the per-order arithmetic is batched into single kernels, and on larger graphs fixed per-epoch costs and better accelerator utilization account for a larger share of the epoch time. On ogbn-arxiv (169K nodes, one H100 GPU) each additional chaos order adds about 4%4\% per-epoch time, and peak GPU memory grows linearly in P+1P+1 (1.68/2.75/3.88/5.021.68/2.75/3.88/5.02 GiB for P=0,…,3P=0,…,3), consistent with the memory analysis above; peak memory on Cora is 75.975.9 MiB. Table 20: Wall-clock time per epoch (ms) vs. chaos order P (Cora and Amazon-Ratings on an A100; ogbn-arxiv on one H100, 3 runs, with peak GPU memory in MiB). Dataset P=0P=0 P=1P=1 P=2P=2 P=3P=3 Ratio P=2P=2/P=0P=0 Cora (2.7K nodes) 35 67 102 131 2.9×2.9× Amazon-Rat. (25K nodes) 464 526 576 646 1.2×1.2× ogbn-arxiv (169K nodes) 2002 2088 2174 2273 1.1×1.1× ogbn-arxiv peak memory (MiB) 1721 2813 3976 5141 2.3×2.3× E.8 Summary of component evidence Table 21 collects, for each architectural component, the controlled comparison that isolates it and the effect measured. The chaos expansion improves calibration (and, at the selected GOOD configurations, shifted accuracy); the spectral filtering design and the residual improve accuracy under heterophily and shift; the energy score with propagation provides OOD detection. Each task’s readout uses a different part of the same representation. Table 21: Component-wise evidence. Each row names one component, the comparison that isolates it, and the measured effect. Component Comparison Effect Chaos expansion (P) P=0P=0 vs best P>0P>0, same architecture (Tables 13, 14, 19) Calibration: Citeseer Brier −27.9%-27.9\% and accuracy +14.8+14.8 points (the endpoint also beats plain GCN, 0.3080.308 vs 0.3180.318 Brier, Tables 1/23); Texas −20.8%-20.8\%; Minesweeper −11.1%-11.1\% Chaos expansion (P) OOD ablation varying only P (Appendix E.2) Energy AUROC flat across P (Cora-structure 91.0/91.2/90.9/90.891.0/91.2/90.9/90.8 at one fixed configuration); detection reads the mean logit, so P can be selected for calibration without changing detection Chaos expansion (P) GOOD ablation at the selected configurations (Table 15) Shifted accuracy: P=1P=1 beats P=0P=0 on GOOD-Cora word (+1.7+1.7), degree (+1.2+1.2), and GOOD-Arxiv time (+0.6+0.6) Spectral backbone (chaos off) P=0P=0 vs GCN (Tables 13, 23) Heterophilous accuracy comes from the spectral backbone: Roman-Empire 79.079.0 at P=0P=0 vs GCN 51.351.3 Filtering design vs single filter Per-dataset architecture at P=0P=0 vs single-branch ChebNet (Table 19) Lower or equal Brier for the full design on 12 of 14 datasets (equal on Questions; Chameleon 0.3620.362 vs 0.7760.776, Squirrel 0.4160.416 vs 0.7980.798, Roman-Empire 0.2980.298 vs 0.3680.368); exceptions are Texas and Wisconsin DSS residual branch (as a unit) Hybrid vs identically trained GCN base, 14 datasets (Table 24) Accuracy: Roman-Empire +27.7+27.7, Minesweeper +6.6+6.6; no dataset loses beyond noise (worst −0.8-0.8) DSS representation under a fixed OOD readout DSS-Hybrid vs GNNSafe++ (same energy readout, propagation, objective, GCN encoder in both; Table 2) Cora 94.3/97.6/94.194.3/97.6/94.1 vs 90.6/95.6/92.890.6/95.6/92.8 Energy-score propagation Same-checkpoint on/off, all 11 settings (Table 16) Helps on 9 of 11, up to +18.2+18.2 AUROC (Twitch) and +14.5+14.5 (Cora-structure) Appendix F Comparison with Standard GNN Baselines Table 1 in the main text compares DSS-GNN against recent uncertainty-aware graph methods (TFE-GNN and G-Δ ). Here we additionally compare against GCN [15] and GAT [34], two widely used GNN architectures that do not incorporate any uncertainty mechanism. Both baselines use 2 layers, 64 hidden units, Adam optimizer (lr==0.01, weight decay==5×10-4), and the same 10 splits as DSS-GNN. Table 23 reports accuracy and Brier score. The main observation is that higher accuracy does not imply better calibration. On Cora, GCN and GAT both achieve higher accuracy than DSS-GNN (87.5%87.5\% and 87.8%87.8\% vs. 86.6%86.6\%) yet have worse Brier scores (0.2130.213 and 0.2220.222 vs. 0.2070.207). On Citeseer, GAT matches DSS-GNN in accuracy (80.6%80.6\% vs. 80.2%80.2\%) but is substantially worse in Brier (0.3530.353 vs. 0.3080.308). On many heterophilous benchmarks, GCN and GAT have much lower accuracy and higher Brier scores than DSS-GNN. The calibration gains of DSS-GNN are therefore not an accuracy effect. Stochastic baselines under the same protocol. Table 22 adds MC-dropout (M=20M=20 stochastic passes at test time, probabilities averaged) and Deep Ensembles (M=5M=5 independently trained models, probabilities averaged) on the GCN of Table 23, under the full 10-split protocol and the same code path; the plain-GCN control from that code path reproduces Table 23 (Cora 87.4787.47 vs 87.4587.45, Wisconsin 50.2550.25 vs 50.2550.25). DSS-GNN’s Brier is lower on all 14 datasets than either baseline, and ensembling improves the GCN’s Brier by at most 0.0080.008 on any dataset while MC-dropout sometimes worsens it (Cora): averaging a conventional GNN does not remove the calibration difference. On Cora and Questions the margins are small (0.2070.207 vs 0.2130.213 and 0.0530.053 vs 0.0550.055). The same multi-pass cost structure applies to sampling-based Bayesian GNNs and latent-variable GNNs; the chaos expansion is instead deterministic, with no sampling variance. Table 22: Brier score (↓ ) with test accuracy (%) in parentheses for DSS-GNN (Table 1) and stochastic GCN baselines under the same 10-split protocol (means over 10 splits). Dataset DSS-GNN MC-dropout Deep Ensemble Plain GCN Cora 0.207 (86.63) 0.2188 (87.36) 0.2134 (87.45) 0.2125 (87.47) Citeseer 0.308 (80.20) 0.3230 (80.03) 0.3170 (79.77) 0.3176 (80.00) PubMed 0.156 (89.72) 0.2176 (86.20) 0.2153 (86.16) 0.2155 (86.10) Texas 0.240 (91.31) 0.5862 (64.43) 0.5686 (62.79) 0.5758 (64.43) Cornell 0.231 (85.11) 0.6829 (52.77) 0.6826 (52.98) 0.6839 (52.34) Wisconsin 0.104 (93.25) 0.6370 (50.38) 0.6373 (51.25) 0.6435 (50.25) Chameleon 0.339 (75.36) 0.7378 (41.73) 0.7371 (42.71) 0.7374 (42.21) Squirrel 0.394 (68.06) 0.7850 (28.11) 0.7847 (28.81) 0.7849 (28.38) CS 0.062 (96.17) 0.0821 (94.55) 0.0809 (94.63) 0.0816 (94.58) Roman-Emp. 0.300 (78.84) 0.6554 (50.92) 0.6471 (51.49) 0.6477 (51.37) Amz-Rat. 0.634 (49.72) 0.6772 (44.08) 0.6769 (44.00) 0.6771 (44.07) Minesweeper 0.169 (87.66) 0.2874 (80.20) 0.2874 (80.20) 0.2874 (80.21) Tolokers 0.268 (79.88) 0.2970 (78.78) 0.2970 (78.78) 0.2971 (78.79) Questions 0.053 (97.17) 0.0548 (97.06) 0.0549 (97.07) 0.0549 (97.07) Table 23: Node classification: test accuracy (%, ↑ ) and Brier score (↓ ) for DSS-GNN vs. standard GNN baselines. Mean ± SD over 10 splits. Best per dataset in bold. Accuracy (%) Brier Score Dataset DSS-GNN GCN GAT DSS-GNN GCN GAT Cora 86.63 ± 1.26 87.45 ± 1.61 87.75 ± 1.71 0.207 ± 0.019 0.213 ± 0.014 0.222 ± 0.015 Citeseer 80.20 ± 1.28 79.84 ± 0.96 80.57 ± 1.12 0.308 ± 0.010 0.318 ± 0.007 0.353 ± 0.007 PubMed 89.72 ± 0.31 86.21 ± 0.31 85.72 ± 0.28 0.156 ± 0.005 0.215 ± 0.005 0.232 ± 0.004 Texas 91.31 ± 3.44 64.43 ± 6.80 71.80 ± 7.06 0.240 ± 0.111 0.576 ± 0.031 0.555 ± 0.038 Cornell 85.11 ± 5.12 52.34 ± 9.99 49.57 ± 12.63 0.231 ± 0.079 0.684 ± 0.080 0.631 ± 0.073 Wisconsin 93.25 ± 3.39 50.25 ± 5.30 56.62 ± 4.71 0.104 ± 0.047 0.643 ± 0.024 0.598 ± 0.021 Chameleon 75.36 ± 1.19 42.21 ± 1.78 43.83 ± 2.54 0.339 ± 0.019 0.737 ± 0.009 0.712 ± 0.015 Squirrel 68.06 ± 2.21 28.39 ± 1.17 27.57 ± 0.75 0.394 ± 0.027 0.785 ± 0.003 0.789 ± 0.002 CS 96.17 ± 0.23 94.61 ± 0.24 94.45 ± 0.18 0.062 ± 0.004 0.082 ± 0.004 0.085 ± 0.004 Roman-Emp. 78.84 ± 0.61 51.30 ± 0.31 42.85 ± 0.53 0.300 ± 0.007 0.648 ± 0.001 0.755 ± 0.003 Amz-Rat. 49.72 ± 0.64 43.97 ± 0.37 42.20 ± 0.43 0.634 ± 0.005 0.677 ± 0.002 0.690 ± 0.001 Minesweeper 87.66 ± 1.10 80.21 ± 0.19 80.10 ± 0.08 0.169 ± 0.014 0.287 ± 0.003 0.281 ± 0.003 Tolokers 79.88 ± 0.49 78.81 ± 0.47 78.19 ± 0.16 0.268 ± 0.005 0.297 ± 0.005 0.299 ± 0.002 Questions 97.17 ± 0.04 97.08 ± 0.02 97.02 ± 0.00 0.053 ± 0.001 0.055 ± 0.000 0.057 ± 0.000 Appendix G Cross-Evaluation of the Two Deployment Modes This appendix gives the full tables for Section 4.5: the hybrid under the calibration protocol, the standalone model under the OOD and GOOD protocols, the readout-coincidence measurement of Section 4.2, the effect of BatchNorm under the calibration protocol, and the covariate-shift counterpart of Table 4. Table 5 (Section 4.5) summarizes both deployment modes on all three tasks; the per-dataset comparisons summarized in its cells are Tables 24, 25, and 26 below. G.1 Protocols Hybrid under the calibration protocol. DSS-Hybrid (2-layer GCN base with the DSS residual branch of Appendix D) is trained under the identical 10-split protocol of Table 1, paired with its own GCN base trained identically with the residual removed and matched budgets, so the pair isolates the residual branch. Standalone under the OOD protocol. Standalone DSS-GNN is trained in the GNNSafe pipeline (200 epochs, 3 runs, energy score with propagation, pipeline defaults: hidden 64, 2 layers, weight decay 10−210^-2, learning rate 10−210^-2, P=1P=1, S=4S=4). Rows marked ++BN add the per-chaos-channel BatchNorm of Appendix D, with every other hyperparameter unchanged; the variant is selected on the validation split. Twitch uses learning rate 10−310^-3, hidden 128, and no BatchNorm, selected on validation ID accuracy. The Arxiv row is trained with the GNNSafe++ energy-margin objective at the Arxiv margin hyperparameters of Wu et al. [36] with mout=−1m_out=-1, selected on validation ID accuracy. All other rows use standard ID training. Standalone under the GOOD protocol. Standalone DSS-GNN is trained in the GOOD pipeline with L1L_1 row-normalized input features, 2 layers, Klp=3K_lp=3, Khp=2K_hp=2, 500 epochs, patience 200, 3 seeds; the epoch is selected by OOD-validation loss within each run and the configuration (BatchNorm on/off, learning rate, hidden width, P) on the GOOD OOD-validation split: CBAS and WebKB ++BN, learning rate 10−210^-2; Twitch ++BN, 10−310^-3; GOOD-Cora/word ++BN, 5×10−35× 10^-3, hidden 128, weight decay 10−310^-3; GOOD-Cora/degree ++BN, 3×10−33× 10^-3, hidden 128; GOOD-Arxiv/time ++BN, 10−210^-2, hidden 256; GOOD-Arxiv/degree no BatchNorm, 10−210^-2, hidden 64; P=1P=1 throughout. G.2 Calibration: standalone, hybrid, and the hybrid’s GCN base Table 24: Test accuracy (%) and Brier score, mean ± SD over the same 10 splits. The standalone column is Table 1; the hybrid and its GCN base are a matched pair (identical training, budgets, and splits; the only difference is the DSS residual). Standalone DSS-GNN DSS-Hybrid GCN base Dataset Accuracy Brier Accuracy Brier Accuracy Brier Cora 86.63 ± 1.26 0.207 ± 0.019 84.33 ± 0.97 0.241 ± 0.015 84.20 ± 1.02 0.243 ± 0.016 Citeseer 80.20 ± 1.28 0.308 ± 0.010 74.71 ± 1.05 0.375 ± 0.009 74.32 ± 1.13 0.378 ± 0.009 PubMed 89.72 ± 0.31 0.156 ± 0.005 88.27 ± 0.56 0.177 ± 0.006 87.93 ± 0.56 0.181 ± 0.006 Texas 91.31 ± 3.44 0.240 ± 0.111 45.57 ± 7.68 0.737 ± 0.034 46.39 ± 7.55 0.737 ± 0.034 Cornell 85.11 ± 5.12 0.231 ± 0.079 55.11 ± 14.03 0.656 ± 0.142 47.66 ± 12.01 0.727 ± 0.068 Wisconsin 93.25 ± 3.39 0.104 ± 0.047 48.12 ± 6.50 0.680 ± 0.029 48.50 ± 6.14 0.680 ± 0.028 Chameleon 75.36 ± 1.19 0.339 ± 0.019 33.50 ± 2.10 0.777 ± 0.004 34.00 ± 2.44 0.777 ± 0.004 Squirrel 68.06 ± 2.21 0.394 ± 0.027 26.35 ± 2.02 0.790 ± 0.003 26.23 ± 2.36 0.790 ± 0.003 CS 96.17 ± 0.23 0.062 ± 0.004 94.26 ± 0.24 0.087 ± 0.004 94.23 ± 0.21 0.088 ± 0.003 Roman-Empire 78.84 ± 0.61 0.300 ± 0.007 74.36 ± 0.61 0.361 ± 0.007 46.63 ± 0.26 0.690 ± 0.003 Amazon-Ratings 49.72 ± 0.64 0.634 ± 0.005 46.36 ± 0.57 0.661 ± 0.004 46.32 ± 0.48 0.661 ± 0.002 Minesweeper 87.66 ± 1.10 0.169 ± 0.014 86.91 ± 0.45 0.184 ± 0.006 80.31 ± 0.21 0.289 ± 0.003 Tolokers 79.88 ± 0.49 0.268 ± 0.005 81.19 ± 0.50 0.255 ± 0.006 79.46 ± 0.34 0.291 ± 0.006 Questions 97.17 ± 0.04 0.053 ± 0.001 97.09 ± 0.04 0.053 ± 0.001 97.08 ± 0.04 0.054 ± 0.001 The hybrid is at or above its own base on 11 of 14 datasets in accuracy; the three negative deltas (Texas −0.82-0.82, Chameleon −0.50-0.50, Wisconsin −0.38-0.38) are far inside one standard deviation, and every positive Brier delta (hybrid worse) is at most 0.00080.0008, also far inside one standard deviation. The gains are structural and large where the GCN inductive bias fails: Roman-Empire +27.7+27.7 accuracy and −0.33-0.33 Brier, Minesweeper +6.6+6.6 and −0.11-0.11, Tolokers +1.7+1.7 and −0.04-0.04. On the small heterophilous graphs where the GCN backbone itself fails (Texas, Cornell, Wisconsin, Chameleon, Squirrel), the residual matches the base rather than improving on it, and the standalone model is the deployment choice. In hybrid mode the argmax of the quadrature predictive and of the mean logit agree on every test node of every dataset (disagreement 0.00000.0000 on all 14). G.3 OOD detection: standalone versus hybrid Table 25: AUROC (%) and ID test accuracy (%) on the 11 OOD settings. Standalone: 3 runs, mean ± SD; hybrid: Tables 2–3 (AUROC only for the nine node-OOD settings∗). The last column gives the standalone training objective. Setting Standalone AUROC Standalone ID acc. Hybrid AUROC Hybrid ID acc. Objective (standalone) Cora / structure 79.79 ± 1.46 (++BN) 69.50 94.32 n/r∗ standard (no OOD exposure) Cora / feature 88.15 ± 1.02 (++BN) 69.17 97.60 n/r∗ standard (no OOD exposure) Cora / label 94.03 ± 1.09 (++BN) 89.24 94.11 n/r∗ standard (no OOD exposure) Amazon-Photo / structure 96.53 ± 1.41 (++BN) 91.59 99.69 n/r∗ standard (no OOD exposure) Amazon-Photo / feature 98.10 ± 0.11 (++BN) 92.00 99.66 n/r∗ standard (no OOD exposure) Amazon-Photo / label 96.53 ± 0.10 (++BN) 94.81 97.52 n/r∗ standard (no OOD exposure) Coauthor-CS / structure 94.33 ± 0.29 (++BN) 92.64 99.96 n/r∗ standard (no OOD exposure) Coauthor-CS / feature 97.79 ± 0.16 (++BN) 92.47 99.95 n/r∗ standard (no OOD exposure) Coauthor-CS / label 95.87 ± 1.83 (++BN) 94.00 98.07 n/r∗ standard (no OOD exposure) Twitch (cross-graph) 51.85 ± 4.83 (no separation) 68.39 ± 0.20 95.75 70.52 standard (no OOD exposure) Arxiv (cross-graph) 73.73 ± 0.24 (++BN) 61.16 ± 0.58 73.36 53.99 GNNSafe++ margin objective ∗Table 2 reports AUROC only. Hybrid ID accuracy on these nine settings (3 runs under the same protocol) is Cora 81.6/81.1/90.681.6/81.1/90.6, Amazon-Photo 92.4/91.0/95.492.4/91.0/95.4, Coauthor-CS 93.0/93.0/95.593.0/93.0/95.5 (structure/feature/label). The hybrid is the stronger detector on 10 of 11 settings; the standalone trains on all 11 (ID accuracies in Table 25) and is competitive on the label-shift settings (Cora/label 94.0394.03 vs 94.1194.11) and on Arxiv, where under the GNNSafe++ margin objective it reaches 73.7373.73 AUROC, above the hybrid’s 73.3673.36, with the best Arxiv ID accuracy of any method in Table 3 (61.1661.16). The one standalone exception is Twitch: the classifier itself trains (ID accuracy 68.3968.39 vs the hybrid’s 70.5270.52), but its energy score does not separate the cross-graph OOD inputs, so the hybrid (95.7595.75) is preferred there. The per-chaos-channel BatchNorm of Appendix D is used on all nine perturbation settings, where the sparse public splits require it for every backbone in the pipeline (Section 4.5). G.4 Distribution shift: standalone versus hybrid Table 26: Shifted (OOD) test accuracy (%) on the GOOD concept-shift settings. Standalone: mean ± SD over 3 seeds, configuration selected on the GOOD OOD-validation split (Appendix G.1); hybrid and best baseline from Table 4. Setting Standalone Hybrid Best baseline GOOD-CBAS / color 82.38 ± 13.17 (++BN) 88.57 TAR 87.29 GOOD-WebKB / university 36.41 ± 1.60 (++BN) 40.77 TAR 30.83 GOOD-Twitch / language 60.28 ± 0.51 (++BN) 61.21 TAR 57.20 GOOD-Cora / word 64.89 ± 0.21 (++BN) 64.81 TAR 64.73 GOOD-Cora / degree 62.60 ± 0.59 (++BN) 62.73 TAR 61.73 GOOD-Arxiv / time 65.43 ± 1.07 (++BN) 66.76 TAR 66.08 GOOD-Arxiv / degree 64.62 ± 0.15 66.43 G-Δ 64.34 The standalone model is above every baseline of Table 4 on 5 of 7 settings (WebKB, Twitch, Cora word, Cora degree, Arxiv degree) and at ERM’s level on the remaining two (ERM: 82.4382.43 on CBAS and 65.6465.64 on Arxiv/time), while the hybrid is above every baseline on all 7; on GOOD-Cora/word the two modes agree within noise (64.8164.81 vs 64.89±0.2164.89± 0.21). GOOD-CBAS is a small synthetic graph with the largest seed variance of the seven settings (standalone SD 13.213.2 over 3 seeds). ID accuracies of the standalone runs are between 63.563.5 and 96.296.2 across the seven settings (CBAS 96.296.2, WebKB 76.076.0, Twitch 63.563.5, Cora word 68.868.8, Cora degree 68.968.9, Arxiv time 74.174.1, Arxiv degree 71.571.5). G.5 Readout coincidence Table 27 compares the two readouts of Section 3.2 at identical checkpoints. Table 27: Quadrature readout versus mean-logit readout at identical checkpoints under the full 10-split protocol of Table 1: absolute Brier difference and the fraction of test nodes on which argmaxp¯i p_i and argmaxZi,0 Z_i,0 disagree (mean over 10 splits). Accuracy differs by at most 0.040.04 points on every dataset. In hybrid mode the disagreement is 0.00000.0000 on all 14 datasets. Dataset |Δ| || Disagree (%) Dataset |Δ| || Disagree (%) Cora 0.0002 0.03 Roman-Emp. 0.0000 0.00 Citeseer 0.0000 0.00 Amz-Rat. 0.0000 0.00 PubMed 0.0000 0.02 Minesweeper 0.0003 0.21 Texas 0.0021 0.00 Tolokers 0.0000 0.00 Cornell 0.0000 0.00 Questions 0.0000 0.00 Wisconsin 0.0000 0.00 CS 0.0001 0.00 Chameleon 0.0001 0.00 Squirrel 0.0007 0.08 G.6 BatchNorm under the calibration protocol Table 28 re-runs the Table 1 protocol with and without the per-chaos-channel BatchNorm on four datasets (10 splits each). The BatchNorm-free controls reproduce Table 1 within a fraction of one standard deviation (the Citeseer run uses the selected P=2P=2 configuration of Table 12, with λreg=0 _reg=0, and reproduces that row as well). With BatchNorm, Cora and Roman-Empire are unchanged within noise and Wisconsin is within one standard deviation, while Citeseer degrades far beyond noise (−13.2-13.2 accuracy). BatchNorm is therefore a stabilization for the sparse-label OOD pipeline, where every baseline backbone also uses it, and is not used under the full-label calibration protocol. Table 28: Test accuracy (%) / Brier, mean ± SD over 10 splits, under the protocol of Table 1 at one calibration configuration per dataset, without and with per-chaos-channel BatchNorm. Dataset BatchNorm off (control) Table 1 BatchNorm on Cora 86.44 ± 1.56 / 0.216 ± 0.022 86.63 ± 1.26 / 0.207 ± 0.019 85.78 ± 1.76 / 0.220 ± 0.017 Citeseer 80.27 ± 1.18 / 0.307 ± 0.010 80.20 ± 1.28 / 0.308 ± 0.010 67.03 ± 1.85 / 0.480 ± 0.024 Wisconsin 93.62 ± 1.72 / 0.109 ± 0.029 93.25 ± 3.39 / 0.104 ± 0.047 91.75 ± 4.75 / 0.134 ± 0.061 Roman-Empire 78.82 ± 1.08 / 0.300 ± 0.012 78.84 ± 0.61 / 0.300 ± 0.007 78.60 ± 0.64 / 0.300 ± 0.009 G.7 Covariate shift on GOOD Table 29 is the covariate-shift counterpart of Table 4: the same seven GOOD settings under their covariate splits, comparing ERM-trained DSS-Hybrid (GCN encoder, configuration selected by validation loss) with the same-protocol ERM GCN of the GNNSafe pipeline (3 seeds). DSS-Hybrid is above the ERM control on 6 of 7 settings; the one loss, GOOD-Cora/degree, is within 1.81.8 points. The two Arxiv ERM baselines did not converge within the shared 120-epoch budget (ID accuracy 45 to 50, versus 72 to 77 for the DSS runs on the same splits), so those two margins are not informative on their own. With the seven concept-shift settings of Table 4 this covers both shift families that GOOD defines. Table 29: Shifted test accuracy (%) under GOOD covariate shift: ERM-trained DSS-Hybrid (GCN encoder, validation-loss-selected configuration) versus the same-protocol ERM GCN. ∗ERM baseline not converged within the shared budget (see text). Setting (covariate split) ERM GCN DSS-Hybrid (ERM) GOOD-CBAS / color 49.05 53.33 GOOD-WebKB / university 13.49 29.10 GOOD-Twitch / language 42.78 52.64 GOOD-Cora / word 61.84 64.18 GOOD-Cora / degree 56.66 54.94 GOOD-Arxiv / time 44.80∗ 70.39 GOOD-Arxiv / degree 25.68∗ 57.60 Appendix H Architecture Overview and Notation Figure 4 gives a schematic of the DSS-GNN architecture of Section 3: the two spectral axes that index the chaos coefficients, one DSS layer, the two readouts and the chaos energy, and the two deployment modes. Table 30 summarizes the notation of the DSS layer, Eqs. (2)–(6). Figure 4: Overview of DSS-GNN. Input features X on graph G are lifted to chaos coefficients Hn(0)=XWin(n)H_n^(0)=XW_in^(n), which are indexed by two spectral axes: the graph Fourier axis (Laplacian eigenpairs (λj,uj)( _j,u_j), filtered by dual Chebyshev filters Glp,GhpG_lp,G_hp) and the Wiener chaos axis (Hermite basis Ψn _n of a scalar Gaussian latent ω). Each DSS layer propagates every chaos channel Hn(ℓ)H_n^( ) through both filters, couples them with per-order chaos gates αn,βn _n, _n, and applies a shared linear map and ReLU at Gauss–Hermite nodes before projecting back (Eqs. (3)–(6)). The final coefficients Zn=Hn(L)WoutZ_n=H_n^(L)W_out are the inputs of the two task readouts, the mean-logit energy score for OOD detection (9) and the quadrature-averaged predictive p¯ p for class prediction and calibration (8), together with the chaos energy ℰE, which is used as regularizer and uncertainty scalar (7) (bottom left). Bottom right: the two deployment modes, standalone DSS-GNN and DSS-Hybrid, which adds the DSS branch residually (scale γres _res) to a warmed-up GCN base (Section 3.4); both use the same readouts. The node index i is suppressed for readability. Table 30: Notation for the DSS layer, Eqs. (2)–(6). Symbol Meaning ω the single shared standard Gaussian latent; never sampled, integrated out by quadrature Ψn _n normalized Hermite polynomial of order n; orthonormal basis of L2(ℝ,γ)L^2(R,γ) P chaos truncation order; n,m=0,…,Pn,m=0,…,P index input/output chaos orders Hn(ℓ)∈ℝN×dℓH_n^( ) ^N× d_ order-n chaos coefficient of the layer-ℓ embeddings (2); the embedding field is H(ℓ)(ω)=∑nHn(ℓ)Ψn(ω)H^( )(ω)= _nH_n^( ) _n(ω) L~G L_G, TkT_k rescaled Laplacian 2LG/λmax−I2L_G/ _ -I and Chebyshev polynomial of degree k (Section 2) Glp,GhpG_lp,G_hp learned low/high-pass filter polynomials of degrees Klp,KhpK_lp,K_hp with coefficients ck(lp),ck(hp)c_k^(lp),c_k^(hp) (3) αn(ω),βn(ω) _n(ω), _n(ω) per-order chaos gates coupling the two filter branches; deterministic scalars, or order-PgP_g expansions with coefficients an,r,bn,ra_n,r,b_n,r (4) ωs,μs,S _s, _s,S Gauss–Hermite quadrature nodes and weights, s=1,…,Ss=1,…,S; exact for the filter/gate polynomials when S≥P+1S≥ P+1 (for Pg≤1P_g≤ 1) H~(ℓ+1)(ωs) H^( +1)( _s) pre-activation field evaluated at quadrature node ωs _s (5) WℓW_ , σ trainable linear map dℓ→dℓ+1d_ → d_ +1 and pointwise activation; Eq. (6) projects σ(H~(ℓ+1)(ωs))σ( H^( +1)( _s)) back onto the Ψm _m basis