Multivariate Time Series Forecasting with Dynamic Graph Neural ODEs

In plain words

Forecasting many related measurements at once, such as the electricity used by hundreds of customers or the speed of traffic at hundreds of road sensors, means tracking two things: how each measurement changes over time, and how the measurements affect one another. Most earlier models handle these with a stack of separate layers, one for time and one for links between measurements, and the links are often taken from a map someone drew by hand.

This paper replaces the stack with a continuous process. The model first turns each series into a set of learned features, then lets those features evolve smoothly over time, the way a quantity changes under a rule that is written as a rate of change. Two such processes run together: one mixes information along time, the other passes information between series over a link structure that the model learns during training. The authors argue that a smoothly changing state is more stable than one that jumps from layer to layer, and that the design needs fewer parameters.

For an implementer, the appeal is that there is no hand-made graph to supply and no special residual or skip wiring. The cost is that the model needs an ODE solver, with step sizes that have to be chosen carefully.

TL;DR

MTGODE forecasts several related series at once. It couples a temporal neural ODE, built from gated dilated convolutions, with a spatial neural ODE that runs over a graph learned from node embeddings §4.1 §4.2. It uses no predefined graph and no stacked discrete residual or skip layers §4.3. The paper reports the lowest RSE in all nine single-step settings on Electricity, Traffic and Solar-Energy, and competitive multi-step errors on Metr-La and Pems-Bay tab. 3 tab. 4. These are author-reported means over ten runs §5.2.

Problem

The task is to learn a spatio-temporal encoder that maps the last T observations of N series, each with D features, to a representation used for H-step forecasts §3 eq. 3. The authors name three limitations of prior work §1:

• Discrete architectures: alternating separately parameterised spatial and temporal blocks gives discontinuous latent trajectories, which they link to higher forecast error §1 fig. 1.
• High complexity: stacked modules need dedicated residual and skip connections and many redundant parameters, which raises compute and memory cost §1 §4.4.
• Reliance on graph priors: GNN forecasters usually assume a predefined static graph, which may be unavailable or a poor fit §1 §2.2.

Related methods each fix only part of this. STGODE uses a continuous graph ODE on a predefined graph but keeps a discrete temporal part §2.3. STG-NCDE is continuous and learns its graph, but is recursive and relies on interpolation preprocessing, which the authors say is costly for long inputs §2.3. MTGNN learns its graph but is discrete §1 §2.2.

Method

MTGODE has three parts: a pre-processing start convolution, a continuous spatio-temporal encoder, and a decoder §4 fig. 2. The encoder integrates an outer temporal ODE over horizon T_cta. The derivative of that ODE is computed by a gated dilated temporal convolution, which feeds an inner spatial ODE on a learned adjacency matrix. The inner ODE's intermediate states are combined by a learned linear map, and the result is zero-padded on the left back to the receptive field length R §4.4 eq. 14 eq. 15.

Steps, in forward order:

1. Start convolution. A 1×1 convolution maps the input series to latent channels, giving the initial temporal state H(0) §4.2 eq. 11.
2. Graph learning. Two embedding matrices per node, each passed through a learned projection and a tanh with saturation β. The adjacency is an antisymmetric score passed through tanh and ReLU, so an edge can point in one direction only §4.1 eq. 8.
3. Spatial ODE. The graph state obeys dH/dt = ( − I)H, where  is the normalised adjacency. It starts from the temporal convolution output and is integrated to time T_cgp with a black-box solver (Euler or Runge-Kutta) §4.1 §4.4 eq. 6.
4. Attentive transform. The output of the spatial ODE is a weighted sum of the states at selected times, including the initial state. The weights are learned matrices Φ_{t_i}, one per time §4.1 eq. 7.
5. Gated temporal dynamics. The temporal derivative is the element-wise product of a tanh filter and a sigmoid gate. Each is a bank of dilated 1-D convolutions with kernel widths 2, 3, 6 and 7 §4.2 eq. 12 eq. 13. The dilation at step t is floor(r^{t/Δt_cta}) eq. 13.
6. Padding. After each aggregation step, the state is zero-padded on the left to length R, where R = L(k−1)+1 when r = 1 and R = 1+(k−1)(r^L−1)/(r−1) otherwise, with L = T_cta/Δt_cta §4.2 eq. 10.
7. Output and loss. The state at the last time index is decoded by an end convolution. Training minimises the mean absolute error (MAE) of the forecast over valid time steps §4.4 eq. 16 alg. 1.

Complexity, as the authors state it: the graph-learning step costs O(Nd² + N²d); each graph propagation step costs O(ED′ + ND′²), with T_cgp/Δt_cgp steps; each temporal step costs O(NQ_l D′²/r), and the full temporal module costs O((T_cta/Δt_cta)·N·R·D′²/r) §4.4. The paper counts one spatial propagation per temporal step.

Forward pass (paper's structure, written in our own notation):

# MTGODE forward pass  [alg. 2] [eq. 14] [eq. 15]
# Input:  X (B, C_in, N, T); times T_cta, T_cgp; dilation r
# Output: forecasts (B, N, H)

def forecast(X):
    # Map each series to D' latent channels  [§4.2] [eq. 11]
    H = start_conv(X)                        # (B, D', N, T)
    # Learn a sparse directed adjacency from node embeddings  [eq. 8]
    A = learn_graph()                        # (N, N)
    # Outer ODE over the time axis, integrated to T_cta  [eq. 14]
    H = solve_ode(temporal_dynamics, H, T_cta)
    # Read the last time index and decode  [§4.4]
    return end_conv(H[..., -1])              # (B, N, H)

# Temporal derivative, one evaluation  [eq. 12] [eq. 13] [eq. 15]
# Input:  H (B, D', N, R), left-padded to length R
# Output: dH/dt (B, D', N, R)

def temporal_dynamics(H):
    # Gated multi-kernel filter: tanh(filter) * sigmoid(gate)  [eq. 12]
    G = tanh(filter_convs(H)) * sigmoid(gate_convs(H))
    # Inner ODE: spatial propagation on the learned graph  [eq. 6]
    S = solve_ode(spatial_dynamics, G, T_cgp)
    # Learned linear mix of the selected intermediate states  [eq. 7]
    S = attentive_mix(S)
    # Zero-pad on the left back to length R  [eq. 10]
    return left_zero_pad(S, R)

# Learned directed adjacency  [eq. 8]
# E1, E2: (N, d) node embeddings; G1, G2: (d, d) projections
# beta: saturation factor (not stated in the paper)

def learn_graph(E1, E2, G1, G2, beta):
    # Squash the two projections of the node embeddings  [eq. 8]
    M1 = tanh(beta * (E1 @ G1))              # (N, d)
    M2 = tanh(beta * (E2 @ G2))              # (N, d)
    # Antisymmetric score; ReLU keeps one direction  [eq. 8]
    score = M1 @ M2.T - M2 @ M1.T            # (N, N)
    A = relu(tanh(beta * score))
    # Sparsify to k neighbours per row  [code: layer.py#L204]
    return keep_top_k_per_row(A, k)

Hyperparameters

Values the paper states are cited to the paper. Values the code uses but the paper does not give are marked "not stated" and cited to the code. The paper gives no time-horizon T for the two solvers; the code uses 1.0.

Results

Author-reported, averaged over ten runs §5.2. Single-step forecasting uses RSE (lower is better) and CORR (higher is better) on Electricity, Traffic and Solar-Energy; the input length is 168 §5.2. Multi-step forecasting uses MAE, RMSE and MAPE on Metr-La and Pems-Bay with 12 input and 12 output steps §5.2. Baselines are in tab. 3 and tab. 4; the rows below show only the strongest ones.

Single-step tab. 3 (RSE ↓ / CORR ↑):

STG-NCDE runs out of memory on Traffic tab. 3.

Multi-step tab. 4 (strongest baselines; full table in tab. 4):

GMAN has a lower Pems-Bay 60-minute MAE than MTGODE (1.86 against 1.88) tab. 4.

Ablation tab. 5 (RSE / CORR; the horizon is not stated in the paper):

The paper's sensitivity studies (hidden size, batch size, the two integration times) are on Metr-La fig. 5. Depth and efficiency studies are on Solar-Energy and Metr-La fig. 4 fig. 6.

Paper vs code

Official code: GRAND-Lab/MTGODE at commit 9011a89c914823e6119c367d0ce7e533f2cc9652, MIT licence LICENSE. The repository links the arXiv paper and contains single-step, multi-step and model code. The paper's footnote points to the same repository §1. The repository ships no datasets README.md.

• Dynamic graph, one propagation vs two. The paper uses one spatial ODE on the learned graph A §4.1 eq. 6. The code runs two spatial ODEs, one on A and one on Aᵀ, and sums them model.py. The paper's complexity count §4.4 covers one.
• Graph learner. The paper's embedding projections have no bias eq. 8. The code's linear layers have biases layer.py layer.py. The paper writes eq. 8 element-wise over N×d matrices; the code multiplies them as matrices, which gives an N×N score layer.py. The card's pseudocode follows the code.
• Sparsity. The paper says the graph is "made sparse" §4.1 without a number. The code keeps the top k = 20 entries per row, after adding 0.01 times uniform noise to break ties layer.py run_multi_step.py.
• Propagation rate. Eq. 6 has unit rate, ( − I)H. The code computes 0.5·α·(ÂH − H) with  = D^(−1/2)(A + I)D^(−1/2) layer.py layer.py. With α = 2.0 (Metr-La) this equals eq. 6. With α = 1.0 (Electricity) the rate is halved electricity.sh. The paper never mentions α.
• Attentive transform. Eq. 7 has a learned matrix per time step and no bias. The code applies one 1×1 convolution with bias to the concatenation of all intermediate states, including the initial state layer.py layer.py. No softmax is used, so "attention" is a learned linear mix.
• Dilation. Eq. 13 gives dilation floor(r^{t/Δt}), a function of time. The code multiplies the dilation by r on every call of the temporal function model.py. For Euler this is once per step; for Runge-Kutta it is once per stage, so the two differ. The paper states no r for its experiments. Single-step defaults to r = 2 run_single_step.py and multi-step to r = 1 run_multi_step.py.
• Gate and dropout. The gate and filter are as in eq. 12–13 (tanh and sigmoid) model.py model.py. Dropout is applied after the gate, a placement the paper does not state model.py.
• Output normalisation. The code applies layer normalisation over channels, nodes and the length-1 time axis before the decoder model.py. Multi-step adds a per-node affine transform model.py; single-step does not run_single_step.py. The paper mentions neither.
• Decoder. The paper calls it an "end convolution" with no detail §4.4. The code uses two 1×1 convolutions with ReLU, widths C/2 and then end_channels, and a final 1×1 convolution model.py model.py.
• Loss. Multi-step trains on MAE in original units, after inverse scaling trainer.py trainer.py. Single-step trains on L1 over max-scaled values, summed rather than averaged run_single_step.py run_single_step.py. The paper gives one MAE loss eq. 16.
• Curriculum. Multi-step enables a curriculum on by default. The loss uses only the first m horizons, and m grows by one every 2500 iterations trainer.py run_multi_step.py run_multi_step.py. The paper does not mention it.
• Node permutation. The code reshuffles node order every 100 iterations run_multi_step.py run_single_step.py. The paper does not mention it.
• Optimiser. Adam with weight decay 1e-4 and gradient-norm clipping at 5 trainer.py trainer.py metr-la.sh. None of these is in the paper. The multi-step schedule is ×0.1 at epoch 100 metr-la.sh; the Electricity schedule is ×0.5 at epochs 20 and 40 electricity.sh. The paper gives only "decay" for multi-step.
• Model selection. The code keeps the checkpoint with the lowest validation loss and reports the test metrics of that checkpoint run_multi_step.py run_multi_step.py run_single_step.py. The paper does not say how the reported model was chosen.
• Input normalisation. Multi-step uses one global standard scaler, fitted on the first input channel of the training set util.py. Single-step divides each series by its maximum absolute value util.py. The paper states neither. Multi-step has two input channels (in_dim = 2) run_multi_step.py; the loader normalises only the first, and the paper does not say what the second is.
• Seeds and runs. The repository sets no random seed [grep: no manual_seed or np.random.seed in *.py]. Scripts run one seed by default run_multi_step.py run_single_step.py. The paper averages ten runs §5.2.
• Reproduced numbers. The README's single run gives Metr-La horizon-3 MAE 2.6491 README.md against 2.66 in tab. 4. It gives Electricity horizon-3 RSE 0.0728 and CORR 0.9431 README.md against 0.0736 and 0.9430 in tab. 3.
• Missing configurations. The repository has scripts for Electricity and Metr-La only. Solar-Energy's settings match the run_single_step.py defaults run_single_step.py. There is no Traffic script and no Pems-Bay script, so the Pems-Bay Runge-Kutta setting cannot be checked.
• Dependencies. The code pins torch 1.7.1 and torchdiffeq 0.2.3 requirements.txt. The paper does not state versions.

Pitfalls

• Dataset gaps in the paper. Section 5.1 says Pems-Bay has 325 sensors §5.1, but Table 2 lists 4,732 nodes tab. 2. The two disagree and the paper does not say which is right, so this card gives no node count for Pems-Bay.
• Integration times are missing. The paper gives only the step sizes Δt. The code uses T_cta = T_cgp = 1.0 in both scripts electricity.sh metr-la.sh. The number of solver steps is round(T/Δt) for Euler model.py layer.py.
• Receptive field can exceed the input. With r = 2 and five CTA steps (Electricity), the formula in the code gives R = 187, which is more than the 168-step input model.py. Inputs are then zero-padded on the left model.py. Multi-step with r = 1 and four steps gives R = 25 against 12 inputs.
• Solver choice changes the dilation schedule. Under Runge-Kutta, the dilation grows once per stage, not once per time unit model.py. The paper's Pems-Bay setting uses Runge-Kutta §5.2, and its script is not in the repository.
• Two spatial propagations cost twice the paper's count. The complexity analysis §4.4 counts one propagation per step. The code runs two model.py. Measured costs in fig. 6 are for the paper's implementation, not a count of the code.
• Attention is not softmax. The "attentive transformation" is a learned linear mix over time-stacked states layer.py. An implementer who expects attention weights that sum to one will get different behaviour.
• Layer normalisation couples nodes. Layer normalisation over channels and nodes together means one node's values affect the scale used for every other node model.py. This is not in the paper and affects the spatial claims.
• Loss scales differ. Single-step loss is a sum over batch and nodes; multi-step is a mean util.py. A gradient-norm threshold of 5 therefore acts on different scales in the two settings trainer.py.
• Zeros are missing values. The loss and metrics ignore target entries equal to zero util.py util.py. Real zero readings are excluded from MAE, RMSE and MAPE, and this choice is not in the paper.
• Model selection uses validation loss. The paper does not state this. Reproduced test numbers depend on it run_multi_step.py.
• Unseeded single runs. The scripts run one unseeded run by default, so expect run-to-run variation. Compare against the paper's ten-run means tab. 3 tab. 4 only after averaging several seeds.
• Ablation horizon is unknown. The full-model row of tab. 5 does not match tab. 3 at horizon 3: Electricity RSE 0.0727 against 0.0736, Traffic 0.4088 against 0.4127, Solar-Energy 0.1686 against 0.1693. The paper does not give the horizon used for the ablation.
• Graph-learner settings come from MTGNN. The paper says it uses "the settings suggested by" MTGNN for the graph learner §5.2 but does not restate them. This card gives only the code values, k = 20 and β = 3 run_multi_step.py run_multi_step.py.
• Numbering. Sections, tables and figures appear with Roman numerals in the printed paper (for example TABLE III). This card uses the LaTeX counters, Arabic numbers 1–6 for sections and tables and figures in print order.

Lineage

• Builds on: neural ODEs as a continuous generalisation of residual networks §2.3 eq. 1 eq. 2; simplified graph propagation (SGC) as the discrete form of the spatial step §4.1 eq. 4; the dilated gated temporal convolution and the multi-kernel filter of MTGNN and Graph WaveNet §4.2 §2.2.
• Compared against: VAR-MLP, GRU, LSTNet, TPA-LSTM, HyDCNN and STG-NCDE for single-step tab. 3; DCRNN, STGCN, Graph WaveNet, GMAN, MRA-BCGN, MTGNN, STGODE and STG-NCDE for multi-step tab. 4.
• Discrete variant: the paper's own discrete counterpart uses the discrete propagation and residual temporal stack fig. 4 tab. 5.
• Superseded by: not known from this wiki.