TWIG: A Time-Causal Wavelet Operator for Autoregressive Forecasting on Irregular Graphs

Subashree Venkatasubramanian; David A. Barajas-Solano; Chuyang Liu; Daniel M. Tartakovsky; Dipankar Dwivedi

cs.LGSep 18, 2026

TWIG: A Time-Causal Wavelet Operator for Autoregressive Forecasting on Irregular Graphs

Authors: Subashree Venkatasubramanian, David A. Barajas-Solano, Chuyang Liu, Daniel M. Tartakovsky, Dipankar Dwivedi

Organizations: Lawrence Berkeley National Laboratory, Berkeley, CA, USA · Stanford University, Stanford, CA, USA · Pacific Northwest National Laboratory, Richland, WA, USA

Abstract

We introduce TWIG (Time-Causal Wavelet Operator for Irregular Graphs), a graph-native neural operator for autoregressive surrogate modeling on static irregular graphs. TWIG transforms each node history into causal multiscale temporal features that separate recent variation from progressively slower memory components, then propagates these features through graph-wavelet operator blocks with gated pointwise channel mixing. The architecture is causal by construction and designed for closed-loop forecasting, where predictions are recursively reused as future inputs. We evaluate TWIG on three irregular-domain forecasting problems spanning regional diffusion, three-dimensional subsurface hydrology, and aerodynamic flow, with graphs ranging from 400 to 5,233 nodes and model capacities from approximately 70k to 10M parameters. TWIG achieves the lowest aggregate rollout errors on the subsurface-hydrology and regional-diffusion benchmarks and ranks second on the 10M-parameter aerodynamic-flow benchmark, behind the GPS Transformer. Across all three settings, TWIG consistently outperforms the corresponding non-time-causal Graph WNO baseline. These results demonstrate that TWIG provides an effective and scalable approach to stable autoregressive forecasting of dynamical fields on irregular graphs.

Figures & tables

Panel of Figure 1 — Figure 1: Architecture of TWIG (Time-Causal Wavelet Operator for Irregular Graphs). a , Overall TWIG architecture and block-autoregressive forecasting procedure. b , Structure of the GraphWNO block used for spatial propagation and pointwise channel mixing.

Figure 2: Autoregressive forecasting across three irregular-domain benchmarks. a , PFLOTRAN test RMSE over a 60-step rollout; the inset enlarges the final five forecast steps to distinguish the leading models. b , SI-diffusion test RMSE over a 100-step rollout. c , Airfoil test RMSE over a 180-step rollout. Curves denote the mean across three independently trained runs and shaded regions denote one standard deviation. Selected competitive baselines and the non-time-causal Graph WNO control are shown for visual clarity; complete benchmark comparisons are reported in the Appendix.

Panel of Figure 3 — Figure 3: Spatial distribution of final-timestep forecast error across the three benchmarks. a , Nodewise PFLOTRAN pressure RMSE at forecast step 60. b , Nodewise absolute infected-fraction error for SI diffusion at forecast step 100. c , Nodewise normalized error for the Airfoil $x$ - and $y$ -velocity components at forecast step 180. Each field is computed over the complete held-out test set and averaged over the three independently trained model seeds. All models within a panel use a common color scale, allowing the spatial extent and concentration of forecast error to be compared directly.

Figure 4: Time-causal multiscale representation and sensitivity to temporal and spatial resolution. a , Causal memory traces for a representative PFLOTRAN salinity history using $K=5$ , showing progressively slower responses across temporal scales. b , Corresponding residual features $u^{(\tau)}-m_{k}^{(\tau)}$ , which encode deviations from memories with different temporal spans. c , PFLOTRAN mean rollout and final-timestep RMSE for Graph WNO and TWIG with $K=3,5,7$ . d , Airfoil mean rollout RMSE as a function of the number $S$ of graph-wavelet spatial scales for TWIG and Graph WNO. Error bars denote one standard deviation where multiple independently trained runs are available.

Figure 5: Robustness of PFLOTRAN forecasts to perturbations of the initial history. The initial ten-state context is perturbed by the indicated percentage of each channel’s training-set standard deviation. The left panel reports mean RMSE over the $60$ -step rollout and the right panel reports final-timestep RMSE. Models are trained only on clean data. Curves report means across three independently trained runs.

Setting	PFLOTRAN	SI diffusion	Airfoil CFD
Dataset
Graph nodes	$2{,}500$	400	$5{,}233$
Scenarios	100	25	$1{,}200$
States per scenario	70	100	200
Train / validation / test	80 / 10 / 10	18 / 3 / 4	$1{,}000$ / 100 / 100
Forecasting and model scale

Table 1: Dataset, forecasting, and training settings for the three benchmarks. All datasets use fixed graph topology within each benchmark, scenario-level data partitions, and training-set normalization.

Method	Parameters (M)	Mean rollout RMSE	Final-timestep RMSE
TWIG ( $K=5,S=4$ )	1.001	$\mathbf{0.0136\pm 0.0001}$	$\mathbf{0.0173\pm 0.0004}$
GPS Transformer	1.015	$0.0140\pm 0.0006$	$0.0177\pm 0.0014$
GAT	1.004	$0.0143\pm 0.0002$	$0.0182\pm 0.0006$
MeshGraphNet	0.976	$0.0145\pm 0.0001$	$0.0177\pm 0.0005$
Graph WNO ( $S=4$ )	1.024	$0.0145\pm 0.0005$	$0.0185\pm 0.0013$
GATv2	0.995	$0.0146\pm 0.0001$	$0.0185\pm 0.0007$

Table 2: Complete PFLOTRAN autoregressive forecasting results. Mean rollout RMSE is averaged over forecast steps 1–60. Parameter counts are reported in millions.

Method	Parameters (M)	Mean rollout RMSE	Final-timestep RMSE
TWIG ( $K=5,S=4$ )	0.07023	$\mathbf{0.0576\pm 0.0049}$	$\mathbf{0.0470\pm 0.0143}$
Graph FNO	0.07086	$0.0602\pm 0.0123$	$0.0541\pm 0.0050$
GPS Transformer	0.06791	$0.0606\pm 0.0111$	$0.0701\pm 0.0111$
Graph WNO ( $S=4$ )	0.06420	$0.0617\pm 0.0089$	$0.0612\pm 0.0163$
GATv2	0.06945	$0.0664\pm 0.0057$	$0.0574\pm 0.0058$
GAT	0.06751	$0.0673\pm 0.0027$	$0.0581\pm 0.0076$

Table 3: Complete SI-diffusion autoregressive forecasting results. Mean rollout RMSE is averaged over forecast steps 1–100. Parameter counts are reported in millions.

Method	Parameters (M)	Mean rollout RMSE	Final-timestep RMSE
GPS Transformer	9.85	$\mathbf{0.1038\pm 0.0080}$	$\mathbf{0.1856\pm 0.0299}$
TWIG ( $K=5,S=12$ )	9.88	$0.1563\pm 0.0035$	$0.2980\pm 0.0152$
MeshGraphNet	10.04	$0.1909\pm 0.0104$	$0.3303\pm 0.0107$
Graph WNO ( $S=4$ )	10.03	$0.2011\pm 0.0042$	$0.3276\pm 0.0054$
Graph FNO	10.00	$0.2026\pm 0.0066$	$0.3363\pm 0.0185$
GAT	9.99	$0.2755\pm 0.0150$	$0.5438\pm 0.0566$

Table 4: Complete Airfoil autoregressive forecasting results. Mean rollout RMSE is averaged over forecast steps 1–180. The main table reports the best tested graph-wavelet spatial-scale configuration for each graph-wavelet architecture:

S=12

for TWIG and

S=4

for Graph WNO.

Model	SI diffusion	PFLOTRAN	Airfoil
TWIG	$0.07023\;(K=5,S=4)$	$1.001\;(K=5,S=4)$	$9.88\;(K=5,S=12)$
Graph WNO	$0.06420\;(S=4)$	$1.024\;(S=4)$	$10.03\;(S=4)$
Graph FNO	$0.07086$	$1.048$	$10.00$
MeshGraphNet	$0.06901$	$0.976$	$10.04$
GAT	$0.06751$	$1.004$	$9.99$
GATv2	$0.06945$	$0.995$	$9.97$

Table 5: Capacity-matched model configurations across the three benchmarks. Parameter counts are reported in millions. Parenthetical values denote the number of time-causal temporal bands

K

and graph-wavelet spatial scales

S

, where applicable.

Experiment	Dataset	Tested settings
Temporal resolution	PFLOTRAN	$K\in\{3,5,7\}$
Spatial resolution	Airfoil	$S\in\{4,8,12\}$
Initial-history noise	PFLOTRAN	$0,1,2.5,5,10\%$

Table 6: Configurations used for the multiscale and robustness experiments.

Neural Surrogate Modeling

Neural surrogate modeling uses learned models to approximate costly simulations, enabling faster prediction, design, control, and inference.

See all 412 papers

Physics-Refined Spatiotemporal Forecasting on Open-Boundary Hydrologic Graphs
Oct 1, 2026Haoyang Jiang, Zhengui Wang, Shenghan Gao +3Spatiotemporal Forecasting Physics-Informed GNNs
William & Mary, Williamsburg, VA, USA · Virginia Institute of Marine Science, William & Mary, Gloucester Point, VA, USA · Florida Atlantic University, Boca Raton, FL, USA
GraphToolbox: A Configurable Python Framework for Graph Neural Network Forecasting
Sep 21, 2026Eloi Campagne, Yvenn Amara-Ouali, Yannig Goude +1Time Series Forecasting Electricity Load Forecasting
CB · LMO, CELESTE · EDF R\&D +1
Towards Scalable One-Step Generative Modeling for Autoregressive Dynamical System Forecasting
May 7, 2026Tianyue Yang, Xiao XuePDE Surrogate Modeling One-Step Generative Modeling
The Center for Computational Science University College London

Oct 1, 2026J/K move · Enter open · S save

Physics-Refined Spatiotemporal Forecasting on Open-Boundary Hydrologic Graphs

Haoyang Jiang, Zhengui Wang, Shenghan Gao, Y. Joseph Zhang+2

William & Mary, Williamsburg, VA, USA · Virginia Institute of Marine Science, William & Mary, Gloucester Point, VA, USA · Florida Atlantic University, Boca Raton, FL, USA

Spatiotemporal Forecasting Physics-Informed GNNs Spatiotemporal GNNs Hydrology

Open paper AskPDF arXiv

Spatiotemporal forecasting on hydrologic graphs is especially prone to instability in open-boundary systems, where the forecast domain exchanges fluxes with an unobserved exterior. In such systems, boundary nodes receive external forcing, e.g., upstream inflows in rivers or tidal signals in coastal regions, that is typically unavailable at prediction time. The absence of this information can compound errors as forecasts unfold in an autoregressive fashion, leading to inferior long-horizon performance. This paper dissects this instability issue by exploring two questions. 1) What boundary forcing enters the forecast domain when information beyond the boundary is missing? 2) How should this forcing propagate through the domain without incurring error amplification under autoregressive rollout? To address both, we propose a new computing framework comprising two key components. First, to compensate for the boundary forcing, our framework learns ghost node proxies from the boundary and interior nodes, striving to approximate unobserved external inputs. Second, to control error accumulation from these learned proxies, we leverage two physics refiners. In particular, one refiner enforces local consistency by aligning ghost proxies with their two-hop neighbors (i.e., boundary nodes and their immediate interiors). The other refiner enhances global stability by correcting the model forecasts through a physics-guided graph neural operator, reducing long-horizon numerical drift. Two real-world hydrologic graphs are employed for empirical evaluation. Comparative results show that our proposal enjoys higher prediction accuracy and long-horizon stability over both learning-based and physics-informed model competitors.

Abstract

Neural Surrogate Modeling

Explore similar work

Physics-Refined Spatiotemporal Forecasting on Open-Boundary Hydrologic Graphs

GraphToolbox: A Configurable Python Framework for Graph Neural Network Forecasting

Towards Scalable One-Step Generative Modeling for Autoregressive Dynamical System Forecasting