Causal discovery from nonlinear multivariate time series is challenging in high-dimensional systems, where the number of candidate directed relations grows quadratically with the number of variables. Existing methods are often limited by target-wise model fitting, repeated conditional testing, or expensive graph-search procedures. We propose SCDM, a shared-parameter meta-learning framework for scalable temporal causal discovery. SCDM treats each target variable as a task while learning a shared temporal predictor and a shared causal-strength matrix across tasks. Each meta-episode updates only a subset of target-variable tasks, allowing the model to accumulate structural evidence across the full system without exhaustive target-wise optimization. After training, a post-training graph readout converts the learned predictor into continuous directed causal scores for threshold-independent evaluation. We also provide a finite-error recovery principle showing that, under fully observed delayed-system assumptions and a positive risk-gap condition, thresholding SCDM scores recovers the graph when statistical, meta-optimization, and readout errors are below the population margin. Experiments on controlled synthetic, neuroimaging-inspired, high-dimensional, and real-world temporal benchmarks demonstrate competitive causal recovery and improved scalability in large temporal systems.
Figures & tables
Figure 1: Overview of SCDM. Each target variable is treated as a task; sampled tasks undergo temporary support adaptation, and their query losses update the shared temporal predictor and causal-strength matrix.
Method
NC8
Simulated EEG
Simulated fMRI
AUROC ↑
AUPRC ↑
AUROC ↑
AUPRC ↑
AUROC ↑
AUPRC ↑
UnCLe
0.947 ± 0.014
0.756 ± 0.080
0.571 ± 0.008
0.155 ± 0.008
0.661 ± 0.012
0.268 ± 0.006
PCMCI
0.909 ± 0.003
0.594 ± 0.022
0.743 ± 0.012
0.390 ± 0.018
0.569 ± 0.012
0.160 ± 0.010
TCDF
0.740 ± 0.019
0.415 ± 0.020
0.685 ± 0.005
0.224 ± 0.007
0.585 ± 0.019
0.178 ± 0.020
GVAR
0.750 ± 0.013
0.413 ± 0.041
0.750 ± 0.004
0.344 ± 0.019
0.619 ± 0.017
0.219 ± 0.019
CUTS+
0.795 ± 0.025
0.354 ± 0.035
0.589 ± 0.013
0.160 ± 0.007
0.532 ± 0.009
0.141 ± 0.009
Table 1: Causal discovery performance on low- and medium-dimensional temporal benchmarks with known ground-truth graphs. We report AUROC and AUPRC as threshold-independent metrics. Results are shown as mean ± 95% confidence interval. Best results are shown in bold and second-best results are underlined .
Figure 2: Qualitative causal recovery on the atlas-grounded simulated fMRI benchmark. (a) Schaefer-100 ROI layout with Yeo-network labels. (b) Network-sorted ground-truth directed graph. (c) Continuous causal scores recovered by SCDM, shown in the same source–target ordering as (b). (d) Representative high-confidence recovered interactions mapped back to ROI space.
Method
200 s
400 s
600 s
800 s
UnCLe (lag 1)
0.573 / 2.05 ×
0.556 / 1.85 ×
0.556 / 1.85 ×
0.541 / 1.53 ×
CUTS+ ( L=3 )
0.571 / 1.31 ×
0.563 / 1.27 ×
0.577 / 1.33 ×
0.586 / 1.36 ×
SCDM
0.566 / 1.96 ×
0.566 / 1.96 ×
0.605 / 2.47 ×
0.645 / 4.01 ×
Table 2: Fixed-budget results on voxel-level fMRI ( d=1000 ; mean over three seeds). Each entry is AUROC / AUPRC enrichment over random (random AUPRC =0.0071 ). At each budget we evaluate the latest checkpoint finished within that budget, so equal values at consecutive budgets mean that no new checkpoint finished in between. SCDM uses the structural readout ∣Wj,i∣ . Best value per budget in bold .
Figure 3: Spatial-plausibility diagnostic on LargeST-GLA. (a) Locations of the sensors used in this experiment; this panel shows geography only. (b) Representative high-scoring directed relations from SCDM. Each arrow points from the inferred source sensor to the inferred target sensor; only relations shorter than 10 km are drawn to reduce clutter, and arrow curvature is for readability only. (c) Distribution of the geographic distance between the two sensors of each top-scoring relation. The median distance falls from about 26 km for random sensor pairs to 9 km after 200 training steps and 2.9 km after 1,000 steps, so high-scoring relations become more local as training proceeds.
Method
Outcome
Time to outcome (s)
Peak GPU (GB)
SCDM
Completed
56.96±4.23
6.9
UnCLe
Out of memory
349.25±13.30
15.6
CUTS+
Time limit
>4800
72.7
JRNGC
Time limit
>4800
14.3
cMLP
Time limit
>4800
38.1
Table 3: Feasibility on LargeST-GLA-3000 (3,000 sensors, 8,997,000 candidate pairs). One A100 80 GB GPU per run, 80-minute limit, three seeds. A run counts as completed only if it writes directed scores for all candidate pairs. Time to outcome is reported as mean ± standard deviation (SD) across three seeds.
Task batch B
AUROC ↑
AUPRC ↑
F1 ↑
SHD ↓
Train (s) ↓
Speedup ↑
10
0.819±0.007
0.656±0.017
0.541±0.030
841.2±30.9
226.9±2.5
9.0×
32
0.867±0.006
0.728±0.019
0.675±0.025
770.6±22.2
668.4±7.2
3.1×
64
0.822±0.007
0.680±0.012
0.588±0.015
785.2±14.9
1301.0±24.8
1.6×
100
0.747±0.008
0.552±0.015
0.392±0.033
1020.8±29.9
2049.9±13.3
1.0×
Table 4: Task-batch-size ablation on the 100-variable simulated fMRI benchmark. All variants use the full method with the structural readout and differ only in the number of target tasks per meta-episode. Results are mean ± 95% CI over five seeds; runtime reports training time only.
Appendix figures & tables32 assets
Supplementary material from the paper’s appendix.
Appendix
Symbol
Definition
T
Number of time points
d
Number of observed variables, which is also the number of target tasks
p
Observed-history order: the input at time t is xt−1,…,xt−p
ℓ
Lag index, ℓ=1,…,p
loss(⋅,⋅)
Pointwise prediction loss (squared error in all experiments)
θ
Shared temporal predictor parameters
Appendix
Table 5: Notation used throughout the paper.
Method
Temporal setting
NC8
Sim. EEG
Sim. fMRI
How chosen
SCDM
observed-history order p
16
5
3
True maximum lag
UnCLe
latent AR lag / kernel
1 / 8
3 / 5
3 / 5
NC8: UnCLe benchmark configuration; others: cover the true lag
cMLP
observed lag order
16
3
3
NC8: UnCLe benchmark configuration; others: lag 3
GVAR
observed lag order
16
3
3
As for cMLP
TCDF
kernel / dilation / max delay
16 / – / –
4 / 4 / 5
4 / 4 / 5
Receptive field covers the true lag
PCMCI
maximum lag τmax
16
5
3
True maximum lag
Appendix
Table 6: Temporal settings of each method. These quantities are not directly comparable: SCDM’s p and the observed lag orders set how many past observations enter the model, UnCLe’s lag acts on latent representations after a convolutional encoder, and convolution kernels set a receptive field. The true maximum lag is 16 for NC8, 5 for simulated EEG, and 3 for simulated fMRI.
Dataset
d
p
B
Lsup
Lqry
λW
N
Split
NC8
8
16
8
64
64
0.05
1,000
70/15/15
Simulated fMRI
100
3
32
64
64
0.05
5,000
70/15/15
Simulated EEG
100
5
32
64
64
0.05
3,000
70/15/15
Voxel-level fMRI
1,000
3
32
32
32
0.03
500
70/15/15
LargeST-GLA-1500
1,500
3
64
64
64
10−4
200
60/20/20
LargeST-GLA-3000
3,000
3
64
64
64
10−4
1,000
60/20/20
Appendix
Table 7: Dataset-specific hyperparameters for SCDM. d : number of observed variables; p : observed-history order; B : target tasks per outer update; Lsup / Lqry : support/query window length; λW : sparsity coefficient (same in support and query losses); N : number of outer updates; Split: train/validation/test fraction.
Dataset
d
q
GT graph
Primary role
NC8
8
0
Yes
Controlled causal recovery
Masked NC8
6
2
Yes
Controlled latent diagnostic
Simulated fMRI
100
0
Yes
Neuroimaging causal recovery
Latent-view simulated fMRI
96
4
Yes
Latent-view observed-block recovery
Simulated EEG
100
0
Yes
Source-space causal recovery
Latent-view simulated EEG
96
4
Yes
Latent-view observed-block recovery
Appendix
Table 8: Summary of datasets used in the experiments. The symbol d denotes the number of observed variables provided to the learner, and q denotes withheld or latent variables when applicable. For LargeST, chronological train/validation/test splits are used and the complete directed causal ground truth is unavailable.
Dataset group
Split protocol
NC8 and masked NC8
Training, validation, and test sequences or replicas are separated before model selection and final reporting.
Simulated fMRI / EEG and latent-view variants
200/50/50 train/validation/test samples are used for controlled recovery experiments.
Voxel-level fMRI
Train, validation, and test intervals are separated; ground-truth graph evaluation is performed on held-out test scores.
LargeST-GLA
Chronological train/validation/test intervals are used to avoid future-to-past leakage.
Appendix
Table 9: Dataset split protocol. Synthetic datasets use train, validation, and test partitions for model fitting, threshold selection, and final reporting. LargeST uses chronological splits to avoid leakage from future intervals.
Dataset group
Causal recovery metrics
Auxiliary use
NC8, simulated fMRI, simulated EEG
AUROC, AUPRC, F1, SHD against known directed graphs
Runtime and controlled ablations
Masked NC8 and latent-view neuroimaging benchmarks
Block-level recovery diagnostics under controlled masking
Latent-extension analysis
Voxel-level fMRI
AUROC/AUPRC against known high-dimensional synthetic graph
Runtime–performance tradeoff and scalability analysis
LargeST-GLA
No causal AUROC/AUPRC because complete directed ground truth is unavailable
Real-world scalability, forecasting, and qualitative graph visualization
Appendix
Table 10: Evaluation role of each dataset group. Controlled synthetic datasets support causal recovery metrics because directed ground-truth graphs are available. LargeST is used for real-world scalability, forecasting, and qualitative structural analysis because complete directed causal ground truth is unavailable.
Figure 4: Overview of the simulated fMRI benchmark. The benchmark uses Schaefer-100 cortical ROIs grouped by Yeo-7 functional networks. The visualization summarizes the atlas-based ROI layout, the network-sorted directed ground-truth graph, representative BOLD-like time series, and the structured lagged effective-connectivity pattern used for evaluation. This figure illustrates the construction of the controlled TVB-inspired, atlas-grounded fMRI simulation rather than empirical subject-level fMRI connectivity.
Figure 5: Overview of the simulated EEG benchmark. The benchmark is constructed as a Schaefer-100 source-space EEG-like simulation with Yeo-network-modular directed source interactions. The visualization summarizes the source ROI organization, the network-sorted source-level ground-truth graph, representative oscillatory source signals, and the synthetic 64-channel source-to-sensor projection used for qualitative visualization. The projection is a simplified source-mixing model and is not intended to represent a subject-specific BEM/FEM EEG forward solution.
Figure 6: Runtime–accuracy tradeoff on the high-dimensional voxel-level fMRI benchmark under fixed wall-clock budgets. Each point reports the mean over three random seeds, without confidence intervals. Left: AUROC computed from continuous causal scores. Right: AUPRC enrichment over random, defined as AUPRC normalized by the random-ranking baseline. UnCLe is slightly ahead at 200 seconds; from 400 seconds on, SCDM has higher AUROC and AUPRC enrichment.
Graph prior
MAE ↓
RMSE ↓
Δ MAE vs. Self-only ↑
Self-only AR
20.60
32.35
—
Random
19.77 ± 0.06
31.11 ± 0.16
+0.83
Distance-kNN
20.37
31.98
+0.23
Road-network
20.68
32.91
− 0.08
UnCLe
19.90
31.38
+0.70
Permuted SCDM
19.92 ± 0.33
31.44 ± 0.77
+0.68
Appendix
Table 11: Supplementary one-step graph-conditioned forecasting on LargeST-GLA-1500 active sensors. All variants use the same graph-masked RidgeVAR forecaster with k=5 , lag p=3 , and fixed ridge regularization α=1.0 . Active sensors are selected using train-only standard deviation σtrain≥10 ( d=1489 ).
Dataset
Readout
AUROC ↑
AUPRC ↑
Simulated fMRI
Learned ∣W∣ score
0.867±0.006
0.728±0.019
Simulated fMRI
Source-isolated perturbation readout
0.832±0.003
0.693±0.010
Simulated EEG
Learned ∣W∣ score
0.713±0.053
0.426±0.086
Simulated EEG
Source-isolated perturbation readout
0.717±0.037
0.476±0.067
Appendix
Table 12: Readout ablation on the simulated fMRI and simulated EEG benchmarks. Both readouts are applied to SCDM outputs. The learned ∣W∣ score evaluates whether the trained causal-strength matrix itself contains useful structural information, while the source-isolated perturbation readout evaluates an alternative post-training graph-extraction procedure. Results are reported as mean ± 95% confidence interval.
Driver a
Mean ∣corr(X1,X2)∣
Condition number
AUPRC ↑
Reversals ↓
0.0
0.037
8.45
0.741
0/5
0.2
0.252
9.91
0.960
0/5
0.4
0.572
12.76
0.955
1/5
0.6
0.750
17.09
0.947
1/5
0.8
0.841
21.81
0.965
1/5
1.0
0.892
26.69
0.915
2/5
Appendix
Table 13: NC8-Collinear stress test (five seeds per level). Reversals count runs in which the non-edge X2→X3 scores above the true edge X1→X3 .
Variant
Task subsample
Shared W
AUROC ↑
AUPRC ↑
F1 ↑
SHD ↓
Full SCDM
✓
✓
0.867±0.006
0.728±0.019
0.675
770.6
w/o task subsampling
–
✓
0.747±0.008
0.552±0.015
0.392
1020.8
w/o shared W (MLP-only)
✓
–
0.732±0.012
0.390±0.022
0.401
1486
Appendix
Table 14: Component ablation on the Schaefer100 TVB-inspired simulated fMRI benchmark. All variants use the same training budget and evaluation protocol. AUROC and AUPRC are reported as mean ± 95% confidence interval; F1 and SHD are reported as means.
Training variant
Inner steps in training
Post-adaptation query MSE ↓
Joint MTL (sampled targets)
–
4.305
Episodic, no inner adaptation
K=0
5.487
Full SCDM
K=3
3.782
Appendix
Table 15: Held-out-target adaptation. All variants share architecture, data, initialization, training budget, and test-time adaptation; only the training objective differs.
Aggregate and latent-to-observed recovery using structural scores
Evaluation target
AUROC ↑
AUPRC ↑
F1 ↑
SHD ↓
Full augmented graph
0.926±0.008
0.634±0.026
0.598±0.057
4.87±0.49
Latent → observed (ZX)
0.606±0.040
0.680±0.028
0.648±0.059
3.80±0.44
Observed-block score diagnostic
Evaluation target
Score readout
AUROC ↑
Observed → observed (XX)
Structural score ∣W∣
0.943±0.007
Appendix
Table 16: Controlled latent recovery on masked NC8. Upper block: full augmented graph and latent-to-observed (ZX) recovery with structural scores. Lower block: observed-to-observed (XX) recovery with the structural score and with the source-isolated readout. Mean ± 95% CI over five graph replicas and three seeds.
Method
Simulated EEG
Simulated fMRI
XX AUROC ↑
XX AUPRC ↑
XX AUROC ↑
XX AUPRC ↑
UnCLe
0.566 ± 0.005
0.146 ± 0.002
0.626 ± 0.033
0.228 ± 0.033
PCMCI
0.745 ± 0.011
0.388 ± 0.019
0.566 ± 0.015
0.156 ± 0.012
TCDF
0.681 ± 0.006
0.213 ± 0.012
0.590 ± 0.009
0.179 ± 0.011
GVAR
0.747 ± 0.003
0.339 ± 0.016
0.619 ± 0.014
0.214 ± 0.014
CUTS+
0.589 ± 0.001
0.154 ± 0.000
0.527 ± 0.012
0.137 ± 0.006
Appendix
Table 17: Observed-block causal recovery under controlled latent-view masking on simulated EEG and fMRI. Results are reported as mean ± 95% confidence interval. Best results are shown in bold and second-best results are underlined .
Statistic
NC8
Masked NC8
Total variables in full system
8
8
Observed variables provided to learner
8
6
Withheld / latent variables
0
2
Time points
2000
2000
Ground-truth graph available
Yes
Yes
Maximum lag
16
16
Appendix
Table 18: Statistics of the NC8 benchmark and its controlled latent-view variant.
Component
Setting
Parcellation
Schaefer-100 cortical atlas
Number of ROIs
100
Coordinate space
MNI152
Functional labels
Yeo-7 canonical networks
Hemisphere labels
Left/right hemisphere labels retained
ROI metadata
ROI name, hemisphere, Yeo network, and MNI centroid
Appendix
Table 19: Atlas and anatomical metadata used in the simulated fMRI benchmark.
Denser directed interactions within the same functional network
Cross-network structure
Sparser directed interactions across functional networks
Lag structure
Directed effects occur over lags 1–5
Appendix
Table 26: Ground-truth source graph construction in the simulated EEG benchmark.
Parameter
Value
Dynamics type
Stable lagged source dynamics with oscillatory components
Lag order
1–5
Time points per sample
2000
Frequency structure
Theta, alpha, and beta-band source activity
Broadband component
1/f -like background activity
Nonstationarity
Slow amplitude modulation and ROI-level variability
Appendix
Table 27: Source dynamics and spectral diagnostics for the simulated EEG benchmark.
Component
Setting
Number of synthetic sensors
64
Sensor layout
Synthetic 10–10-like layout
Leadfield / mixing matrix
leadfield_synthetic.npy , shape 64×100
Sensor positions
sensor_positions.npy , shape 64×3
Full sensor signal
EEG.npy , generated from all 100 source ROIs
Observed-only sensor signal
EEG_obs96.npy , generated from 96 observed source ROIs
Appendix
Table 28: Synthetic source-to-sensor projection used in the simulated EEG benchmark.
Statistic
Full-view
Latent-view
Total source ROIs in full system
100
100
Observed variables provided to learner
100
96
Withheld / latent variables
0
4
Time points per sample
2000
2000
Ground-truth graph available
Yes
Yes
Graph type
Directed lagged source EC
Directed lagged source EC
Appendix
Table 29: Full-view and latent-view versions of the simulated EEG benchmark.
Index
ROI name
5
7Networks_LH_Vis_6
58
7Networks_RH_SomMot_1
86
7Networks_RH_Cont_Cing_1
41
7Networks_LH_Default_PFC_1
Appendix
Table 30: Withheld ROIs in the latent-view simulated EEG setting.
Block
Description
A_xx.npy
Directed edges among observed source ROIs, X→X
A_xz.npy
Directed edges from observed source ROIs to hidden source ROIs, X→Z
A_zx.npy
Directed edges from hidden source ROIs to observed source ROIs, Z→X
A_zz.npy
Directed edges among hidden source ROIs, Z→Z
Appendix
Table 31: Ground-truth graph blocks in the latent-view simulated EEG setting.
Method
AUROC / AUPRC
Interpretation
Random baseline
0.500 / 0.1149
Expected from graph density
Lagged correlation
0.6292 / 0.1826
Simple temporal association has partial signal
Ridge-VAR
0.6983 / 0.2668
Linear VAR recovers part of the lagged structure
Ours, 3000 steps
0.6933 / 0.4598
Stronger sparse-edge ranking in AUPRC
Appendix
Table 32: Recoverability diagnostics for the simulated EEG benchmark. The “Ours” row is a one-seed smoke diagnostic on the full 100-source non-latent setting and is not a replacement for the final multi-seed experimental table.
Causal discovery with instantaneous effects in multivariate time series is challenging, as the instantaneous structure must be acyclic. Prior methods enforce this by either separating instantaneous and lagged estimation into multi-stage pipelines or imposing algebraic acyclicity constraints via complex augmented Lagrangian optimization, both of which incur high computational cost. In this work, we propose a different approach: we learn a differentiable permutation of variables using the Gumbel--Sinkhorn operator and triangularize the instantaneous coefficient matrix of a Structural Vector Autoregressive (SVAR) model in the learned order. This converts acyclicity from a hard constraint into a parameterization and keeps it valid throughout optimization. In doing so, our method enables unified, continuous optimization with gradient-based learning, leading to improved efficiency in time--series causal discovery. Across three real-world benchmarks, our method achieves the best overall performance compared with 12 baselines in both discovery accuracy and efficiency. On the large-scale benchmark, it further demonstrates strong scalability, achieving more than a 6x speedup over competing methods.
Causal discovery for both cross-sectional and temporal data has traditionally followed a dataset-specific paradigm, where a new model is fitted for each individual dataset. Such an approach limits the potential of multi-dataset pretraining. The concept of large causal models (LCMs) envisions a class of pre-trained neural architectures specifically designed for temporal causal discovery. Prior approaches are constrained to small variable counts, degrade with larger inputs, and rely heavily on synthetic data, limiting generalization. We propose a principled framework for LCMs, combining diverse synthetic generators with realistic time-series datasets, allowing learning at scale. Extensive experiments on synthetic, semi-synthetic and realistic benchmarks show that LCMs scale effectively to higher variable counts and deeper architectures while maintaining strong performance. Trained models achieve competitive or superior accuracy compared to classical and neural baselines, particularly in out-of-distribution settings, while enabling fast, single-pass inference. Results demonstrate LCMs as a promising foundation-model paradigm for temporal causal discovery. Experiments and model weights are available at https://github.com/kougioulis/LCM/.
Nikolaos Kougioulis, Nikolaos Gkorgkolis, MingXue Wang +4
Institute of Applied & Computational Mathematics, FORTH · Computer Science Department, University of Crete · Huawei Ireland Research Centre, Dublin, Ireland
Scalable nonlinear causal discovery requires methods that combine flexible mechanism estimators with efficient search over large graph spaces. Several algorithmic families have been proposed to address this challenge, yet their accuracy-runtime trade-offs remain poorly understood. We empirically compare the four major approaches: differentiable structure learning, amortized structure learning, score-matching, and combinatorial search. Our results reveal complementary bottlenecks: differentiable and amortized methods scale well but exhibit an accuracy gap, score-matching methods can be accurate in low dimensions but degrade quickly for increasing feature sizes, and combinatorial methods remain accurate but are slowed by repeated and redundant local scoring. Motivated by this bottleneck, we develop SPADE, a spline-based score-evaluation scheme that compiles sufficient statistics once and reuses them throughout combinatorial search. Under bounded indegree, its Gaussian variant reduces algorithmic complexity from O(nd^3) to O(nd^2+d^3). Empirically, SPADE shifts the observed scalability-accuracy frontier by orders of magnitude: it solves 100-variable problems with 160K samples in seconds and 1600-variable problems with 2.5K samples in minutes, while retaining high structural accuracy across synthetic and real-world benchmarks. These results reveal a substantial shift in the practical scale of combinatorial search and highlight the importance of evaluating scalable causal-discovery methods along the full accuracy-runtime frontier.