MillGNN: Multi-Scale Lead-Lag Forecasting

• The paper introduces MillGNN, a graph neural network that models multivariate time series by explicitly capturing fine- and coarse-scale lead-lag dependencies via FFT-based cross-correlation.
• It employs a structured pipeline combining statistical initialization, decay-aware attention, and hierarchical message passing to efficiently propagate messages across multiple grouping scales.
• Empirical evaluations on diverse benchmarks demonstrate that MillGNN outperforms state-of-the-art models in both long-term and short-term forecasting tasks.

Searching arXiv for MillGNN and closely related forecasting papers. MillGNN is a graph neural network-based method for multivariate time series forecasting that learns multiple grouping scale lead-lag dependencies. It is designed for settings in which hierarchical propagation produces delayed effects across variables and groups of variables, and it combines statistically grounded lead-lag detection with structured multi-scale message passing. Formally, given historical observations $\boldsymbol{X}\in\mathbb{R}^{N\times L}$ for $N$ variables over input length $L$, MillGNN learns a forecasting function $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$ with $\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$ for prediction horizon $H$ (Wu et al., 4 Sep 2025).

1. Problem setting and conceptual basis

MillGNN addresses multivariate time series forecasting under explicit lead-lag structure. In this setting, a directed dependency means that segments of a leading variate affect a lagging variate after a time delay $\tau$; for a pair $(i,j)$, the effect from $x_i(t)$ to $x_j(t+\tau)$ means that $N$0 leads $N$1 at lag $N$2. The method is motivated by the observation that existing methods often overlook lead-lag dependencies at multiple grouping scales, and therefore fail to capture hierarchical lead-lag effects in complex systems (Wu et al., 4 Sep 2025).

The central modeling assumption is that complex systems can exhibit both fine-scale and coarse-scale propagation. At fine scales, individual variate interactions are sensitive and dynamic; at coarse scales, fluctuations are smoothed and lead-lag relations become more stable. MillGNN therefore models both variate-wise and group-wise lead-lag dependencies. This suggests a division of labor across scales: fine scales capture localized and rapidly changing dependencies, whereas coarser scales summarize longer-range and more stable propagation patterns.

The underlying graph view is time-resolved and causal. Nodes are defined as $N$3, where time steps are nodes, and a directed edge $N$4 encodes influence from $N$5 to $N$6 with lag $N$7, subject to causality. The weighted adjacency $N$8 is upper triangular per pair $N$9 because

$L$0

This causal structure is retained after scale-specific aggregation and patchification (Wu et al., 4 Sep 2025).

2. Multi-scale grouping and scale-specific lead-lag graph learning

A defining feature of MillGNN is its multi-scale grouping mechanism. Variates are grouped into $L$1 scales based on a binarized Dynamic Time Warping similarity graph. If $L$2 denotes the initial binary similarity graph from the training set, then the assignment matrix at scale $L$3 is

$L$4

The similarity structure and signals at scale $L$5 are then

$L$6

Patchification is scale-dependent:

$L$7

Larger patch lengths $L$8 are used at coarser scales to capture longer delays (Wu et al., 4 Sep 2025).

The scale-specific lead-lag graph learning module, SiLL-GL, begins with FFT-based cross-correlation. Using Wiener–Khinchin, MillGNN computes

$L$9

where $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$0 is FFT, $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$1 is inverse FFT, $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$2 is element-wise product, and the overline denotes complex conjugation. In practice, cross-correlations are computed on patch-averaged signals per scale, and the Top-$$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$3 lag values are selected for each candidate pair (Wu et al., 4 Sep 2025).

Candidate pairs at scale $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$4 are restricted to groups that share the same parent at scale $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$5. For each such pair $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$6, the lag set is

$$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$7

The initial binary adjacency $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$8 over patches $$\hat{\boldsymbol{X}}=f(\boldsymbol{X};\boldsymbol{\Theta})$$9 is defined by

$\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$0

Because of causality, $\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$1 is upper triangular (Wu et al., 4 Sep 2025).

MillGNN then augments this statistically initialized graph with a decay-aware attention mechanism. Group embeddings for excitatory and inhibitory characteristics are introduced as $\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$2, together with time embeddings $\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$3. Broadcasting addition yields

$\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$4

with $\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$5. Rate matrices are computed through SoftMax over patch interactions:

$\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$6

Real-time input attention is derived from patchwise features:

$\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$7

A lag tensor $\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$8 stores absolute patch differences. The dynamic decaying weight tensor is then

$\hat{\boldsymbol{X}}\in\mathbb{R}^{N\times H}$9

and the scale-specific lead-lag adjacency becomes

$H$0

The intended interpretation is explicit: larger $H$1 strengthen decay, while larger $H$2 and $H$3 counteract decay, modeling periodic revival (Wu et al., 4 Sep 2025).

3. Hierarchical lead-lag message passing

MillGNN’s second core module is hierarchical lead-lag message passing, HiLL-MP. At each scale, patchified signals are embedded as

$H$4

Hierarchical aggregation is then performed back to the variate level. The aggregate message is

$H$5

Here, $H$6 is an activation function; $H$7 are linear maps; $H$8 maps scale-$H$9 groups to variates; $\tau$0 is the number of variates under each scale-$\tau$1 group; and $\tau$2 is a duplication operator that efficiently broadcasts group messages to assigned variates (Wu et al., 4 Sep 2025).

This structure separates two propagation modes. Intra-scale propagation acts among similar variates or groups within the same scale. Inter-scale propagation moves information progressively from coarser groups down to finer ones through duplication along the learned hierarchy. The paper characterizes this design as balancing comprehensiveness and efficiency, since it transmits lead-lag effects across multiple grouping scales without exhaustive pairwise multiplication (Wu et al., 4 Sep 2025).

After aggregation, variate-level and group-level features are updated. For the variate level,

$\tau$3

For groups,

$\tau$4

Multiple HiLL-MP layers are stacked with skip connections, producing the final variate-level output $\tau$5 (Wu et al., 4 Sep 2025).

The forecasting head follows PatchTST. MillGNN flattens patches and hidden channels and linearly projects to the forecast horizon:

$\tau$6

Training uses Mean Squared Error. Evaluation uses MSE and MAE for long-term tasks and MAE/MAPE/RMSE for short-term traffic and air-quality tasks (Wu et al., 4 Sep 2025).

4. Computational properties and algorithmic pipeline

MillGNN is designed to preserve tractability despite explicit lag modeling. FFT-based cross-correlation reduces per-pair complexity at scale $\tau$7 from $\tau$8 to $\tau$9 after patch pooling, and cross-correlation is computed only for candidate pairs. The total message passing complexity per forward pass is stated as $(i,j)$0, based on the facts that total nodes satisfy $(i,j)$1, edges per scale are linear in $(i,j)$2 given bounded lag window and group size $(i,j)$3, graph message passing scales as $(i,j)$4, and inter-scale duplication is $(i,j)$5 per group-to-variates broadcast. Adjacency tensors per scale require $(i,j)$6 space; practical memory is controlled by candidate-pair filtering, patch sizes $(i,j)$7, and small $(i,j)$8, typically $(i,j)$9–$x_i(t)$0 in practice (Wu et al., 4 Sep 2025).

The training and inference workflow has four main stages. First, data preprocessing normalizes series, computes DTW similarities on the training set, binarizes them to form $x_i(t)$1, clusters to build the assignment matrices $x_i(t)$2, aggregates signals per scale, and patchifies with scale-specific patch length $x_i(t)$3. Second, SiLL-GL identifies candidate group pairs, computes FFT-based cross-correlation, selects Top-$x_i(t)$4 lags, builds $x_i(t)$5, computes the excitatory and inhibitory embeddings and real-time attention, forms the lag tensor, and constructs $x_i(t)$6. Third, HiLL-MP computes $x_i(t)$7, aggregates intra-scale and inter-scale messages, updates variate-level and group-level features, stacks layers with skip connections, and retains $x_i(t)$8. Fourth, the forecast head computes $x_i(t)$9 and optimizes MSE with Adam, using early stopping via validation; inference repeats the same stages and outputs $x_j(t+\tau)$0 (Wu et al., 4 Sep 2025).

Implementation details reported for the experiments include PyTorch, a single NVIDIA A100 80GB GPU, Adam optimization, initial learning rates in $x_j(t+\tau)$1, batch size $x_j(t+\tau)$2, $x_j(t+\tau)$3 epochs, and early stopping tolerance $x_j(t+\tau)$4. Hyperparameters are searched with NNI using a Bayesian tuner with $x_j(t+\tau)$5 trials and a curve fitting assessor for early termination. Search spaces are dataset dependent: for small-$x_j(t+\tau)$6 datasets, scales in $x_j(t+\tau)$7 and groups at scale $x_j(t+\tau)$8 in $x_j(t+\tau)$9; for large-$N$00 datasets, scales in $N$01, groups at scale $N$02 in $N$03, and groups at scale $N$04 in $N$05. Patch length satisfies $N$06 and Top-$N$07 lags satisfy $N$08 (Wu et al., 4 Sep 2025).

5. Empirical performance and ablation results

MillGNN is evaluated on 11 datasets. The long-term forecasting benchmarks are ETTh1, ETTh2, ETTm1, ETTm2, Weather, Electricity, Exchange, and Traffic, with input length $N$09 unless noted and horizons $N$10. The short-term benchmarks are PEMS04, PEMS08, and China-AQI, with $N$11 for PEMS and $N$12 for China-AQI. The baselines comprise 16 state-of-the-art methods: Autoformer, FEDformer, PatchTST, Stationary, DLinear, SCINet, TimesNet, TimeMixer, MTGNN, Crossformer, CrossGNN, MSGNet, iTransformer, FilterTS, LightMTS/LIFT, and VCformer (Wu et al., 4 Sep 2025).

Across the long-term benchmarks, MillGNN reports the following average results: ETTm1, MSE $N$13 and MAE $N$14; ETTm2, MSE $N$15 and MAE $N$16; ETTh1, MSE $N$17 and MAE $N$18; ETTh2, MSE $N$19, described as close to best while TimeMixer attains $N$20; Electricity, MSE $N$21 and MAE $N$22; Exchange, MSE $N$23 and MAE $N$24; Traffic, MSE $N$25 and MAE $N$26; and Weather, MSE $N$27 with MAE approximately $N$28 near best. On ETTm1 at horizon $N$29, MillGNN attains MSE $N$30 and MAE $N$31; at $N$32, MSE $N$33 and MAE $N$34. The reported “Count” line assigns MillGNN $N$35 best scores for both MSE and MAE (Wu et al., 4 Sep 2025).

On short-term forecasting, MillGNN reports the best performance across all three metrics on all three datasets: PEMS04 with MAE $N$36, MAPE $N$37, RMSE $N$38; PEMS08 with MAE $N$39, MAPE $N$40, RMSE $N$41; and China-AQI with MAE $N$42, MAPE $N$43, RMSE $N$44 (Wu et al., 4 Sep 2025).

The ablation results isolate the contribution of the major components. Removing multi-scale modeling yields ETTm1 MSE $N$45 instead of $N$46, Electricity MSE $N$47 instead of $N$48, and China-AQI MAE $N$49 instead of $N$50. Removing the initial graphs based on cross-correlation yields ETTm1 $N$51, Electricity $N$52, and China-AQI $N$53. Removing decay-aware weights yields ETTm1 $N$54, Electricity $N$55, and China-AQI $N$56. Replacing HiLL-MP with a vanilla GNN yields ETTm1 $N$57, Electricity $N$58, and China-AQI $N$59. The full design consistently yields the lowest errors, and the decay-aware weights are reported as particularly impactful on China-AQI (Wu et al., 4 Sep 2025).

Robustness analyses further qualify the performance claims. On ETTh1 and ETTh2, MillGNN excels when train-test lag distributions are stable, as on ETTh1, but performance degrades when lag distributions drift, as observed on ETTh2. Hyperparameter analysis on China-AQI identifies an optimal number of scales around $N$60, Top-$N$61 lags around $N$62, and decay embedding dimension around $N$63; too small $N$64 or $N$65 underfit, while too large $N$66 or $N$67 risk overfitting. Among inter-scale strategies, duplication is reported as both accurate and efficient relative to learnable assignment or bi-directional passing, the latter being described as overmixing and slowing training (Wu et al., 4 Sep 2025).

6. Interpretability, limitations, and relation to adjacent literature

MillGNN is presented as combining statistical interpretability with learned dynamic weighting. Interpretability begins with the initial edges, which arise from FFT-based cross-correlation, a well-understood statistic for time-lagged co-movements. It continues through the decay-aware weights, whose explicit lag-dependent form

$N$68

encodes decay over longer lags, revival under strong attention, and suppression under high inhibitory rates. Case studies on China-AQI report learned variate-to-group assignments, including cities grouped by similar AQI dynamics, and lead-lag edges matching observed delays; one example states that wind direction over Shanghai explains lagged AQI peaks in Ningbo, and that MillGNN detects these lags purely from AQI without exogenous wind inputs while outperforming the listed state-of-the-art methods (Wu et al., 4 Sep 2025).

The reported limitations are correspondingly specific. MillGNN depends on the presence of meaningful lead-lag structure; if lead-lag relations are weak or dominated by synchronous interactions, gains may be limited. Performance can degrade under regime shifts and lag distribution drift. Multi-scale graphs and FFT-based cross-correlation add overhead relative to purely channel-independent models, although HiLL-MP is intended to mitigate this. Grouping quality also matters: spectral clustering is described as accurate but costly, whereas K-means offers better efficiency with a modest accuracy trade-off (Wu et al., 4 Sep 2025).

In relation to other forecasting architectures, MillGNN is distinguished from graph-based models such as MTGNN, CrossGNN, and MSGNet by explicitly modeling lead-lag edges across patches and scales, using decay-aware dynamic weights informed by real-time inputs and lag magnitude, and employing structured hierarchical message passing through coarse-to-fine duplication. Relative to Transformer-based patching models such as PatchTST, iTransformer, and Crossformer, it retains patching but adds multi-scale lead-lag graphs with explicit lag-aware decay and hierarchical propagation. Relative to decomposition-based approaches such as DLinear, Autoformer, FEDformer, and TimesNet, it is characterized by explicit lag modeling rather than predominantly synchronous trend, seasonal, or frequency decomposition. Relative to single-scale lead-lag methods such as LightMTS/LIFT and VCformer, it advances by discovering and using lead-lag dependencies at multiple grouping scales and by learning dynamic decay-aware weights (Wu et al., 4 Sep 2025).

The name “MillGNN” also has a narrower terminological boundary than might be assumed from secondary discussions around graph neural networks for mixed-integer linear programs. The paper "On Representing Mixed-Integer Linear Programs by Graph Neural Networks" studies MILP graph encodings, WL-bounded expressivity, foldable versus unfoldable instances, and random-feature symmetry breaking, but it does not introduce a single named model called “MillGNN” (Chen et al., 2022). A plausible implication is that the established use of “MillGNN” in the arXiv literature refers to the multivariate time series forecasting model rather than to the MILP-representation theory framework.

7. Practical use cases and methodological significance

The practical guidance reported for MillGNN centers on systems with hierarchical propagation and evident delays, including transportation networks, power grids, and environmental diffusion, especially when datasets contain many variables and known group structures such as spatial regions or functional subsystems. Recommended starting values from the hyperparameter study are $N$69 and $N$70, with larger $N$71 suggested for large $N$72 when clear group hierarchies are present. Patch lengths $N$73 should reflect sampling rates and expected lag ranges, with longer patches at coarser scales for longer delays. Spectral clustering is preferred when compute permits, while K-means is suggested for efficiency. Domain priors, such as geography or network topology, may be incorporated into $N$74 to guide grouping. Efficiency recommendations include restricting candidate pairs to siblings under the same parent at scale $N$75, keeping $N$76 and $N$77 modest, and avoiding bi-directional inter-scale passing (Wu et al., 4 Sep 2025).

Methodologically, MillGNN occupies a specific position in the forecasting literature. Its distinctive claim is not simply that graph neural networks improve multivariate time series prediction, but that lead-lag structure should be represented as a causal, multi-scale, lag-aware graph whose edges are initialized by cross-correlation and refined by dynamic, decay-aware weights. The empirical results indicate that this combination is especially effective on datasets with spatially structured propagation, including Traffic, PEMS, and China-AQI. This suggests that the method is particularly well matched to domains in which delayed transmission across related subsystems is an organizing feature of the data-generating process rather than a minor deviation from synchronous dependence (Wu et al., 4 Sep 2025).