FoundCause: Causal Discovery with Latent Confounders from Observational Data

# FoundCause: Causal Discovery with Latent Confounders from Observational Data
Source: [https://arxiv.org/html/2606.17516](https://arxiv.org/html/2606.17516)
Patrick Blöbaum∗,1, Krishnakumar Balasubramanian∗,1,2Shiva Prasad Kasiviswanathan1 1Amazon Web Services 2Department of Statistics, University of California, Davis

###### Abstract

Causal discovery from observational data remains challenging due to the need to recover directed structure and latent confounding without interventions\. We proposeFoundCause, anamortized causal discovery modeltrained entirely on synthetic data that maps datasets directly to causal graphs in a single forward pass\. By learning from large collections of simulated structural causal models,FoundCausecaptures transferable statistical patterns that generalize beyond individual datasets\. The architecture incorporates several key inductive biases for causal discovery\. It uses a permutation\-invariant transformer encoder with alternating attention over samples and variables to jointly model cross\-variable dependence and per\-variable distributions\. Pairwise statistical features derived from classical asymmetry measures are injected through statistics\-conditioned attention, guiding the model toward known causal signals\. A factorized decoder separates edge existence from direction, while a triangular refinement module enables reasoning over higher\-order causal motifs such as chains and colliders\. In addition, a dedicated confounder module based on learnable latent tokens explicitly models hidden common causes, and the model explicitly handles missing data via its masked input representation\. To our knowledge,FoundCauseis the first amortized causal discovery approach to explicitly model latent confounding\.FoundCauseoutperforms 11 classical non\-amortized methods \(e\.g\., PC, GES, NOTEARS\-style optimization\) and 4 amortized causal discovery methods on 15 real\-world datasets, achieving \+9\.6% improvement inF1F\_\{1\}, \+1\.2% in AUROC, and an 18\.9% reduction in structural Hamming distance relative to the strongest non\-amortized methods, while performing inference in a single forward pass\.

## 1Introduction

Causal discovery from observational data is a central problem in machine learning, statistics, and many applied domains, requiring the recovery of directed relationships and latent confounding without access to controlled interventions\. Given a dataset of i\.i\.d\. samples, the goal is to infer the underlying causal graph governing the data\-generating process, typically formalized as a structural causal model \(SCM\)\.

Classical approaches such as constraint\-based methods \(e\.g\.,Spirteset al\.\([2000](https://arxiv.org/html/2606.17516#bib.bib1)\); Zhang \([2008](https://arxiv.org/html/2606.17516#bib.bib6)\)\) and score\-based methods \(e\.g\.,Chickering \([2002](https://arxiv.org/html/2606.17516#bib.bib3)\)\) rely on conditional independence testing or combinatorial search, while more recent continuous optimization approaches \(e\.g\.,Zhenget al\.\([2018](https://arxiv.org/html/2606.17516#bib.bib9)\); Nget al\.\([2020](https://arxiv.org/html/2606.17516#bib.bib10)\)\) formulate causal discovery as a differentiable problem\. A complementary line of work leverages functional assumptions on the data\-generating mechanisms to achieve identifiability \(e\.g\.Shimizuet al\.\([2011](https://arxiv.org/html/2606.17516#bib.bib2)\); Hoyeret al\.\([2009](https://arxiv.org/html/2606.17516#bib.bib62)\); Zhang and Hyvärinen \([2009](https://arxiv.org/html/2606.17516#bib.bib63)\)\)\. Although these methods are well studied, they typically require expensive per\-dataset optimization, are sensitive to hyperparameters and distributional assumptions, and often struggle to scale to high\-dimensional or heterogeneous data regimes\. As a result, their applicability in real\-world settings remains limited\. More fundamentally, causal discovery from observational data is not identifiable in general without additional assumptions\. Yet these assumptions are often violated in real\-world data, leading to a persistent gap between theoretical identifiability and practical applicability\. This motivates inference procedures that combine signals across the full spectrum of identifiability regimes, rather than committing to any single one\.

A recent line of work has reframed causal discovery as an*amortized inference*problem \(e\.g\.,Lorchet al\.\([2022](https://arxiv.org/html/2606.17516#bib.bib21)\); Keet al\.\([2023](https://arxiv.org/html/2606.17516#bib.bib38)\); Montagnaet al\.\([2025](https://arxiv.org/html/2606.17516#bib.bib39)\); Wuet al\.\([2025](https://arxiv.org/html/2606.17516#bib.bib42)\); Mahajanet al\.\([2025](https://arxiv.org/html/2606.17516#bib.bib40)\)\), where a model is trained on large collections of synthetic datasets and learns a direct mapping from data to graph structure\. This paradigm enables single\-pass inference and improved scalability, and offers a practical alternative to per\-dataset optimization\. However, existing neural approaches often rely on generic architectures with limited incorporation of domain\-specific causal structure\. In particular, they tend to treat causal discovery as a representation learning problem, without explicitly leveraging well\-established statistical signals such as asymmetries in cause–effect relationships or higher\-order structural patterns like chains and colliders\. As a result, these methods can lack robustness and interpretability, especially when generalizing beyond the training distribution\.

In this work, we proposeFoundCause, a hybrid architecture that combines amortized inference with explicit statistical inductive biases\. Rather than pursuing new identifiability guarantees, our goal is to learn practical inference procedures that perform well across a broad range of data\-generating processes\.FoundCauseprocesses tabular data using a permutation\-invariant transformer encoder with alternating attention over samples and variables, augmented by stat\-conditioned attention that injects pairwise causal features derived from classical methods\. The model predicts causal structure via a factorized decoder that separates edge existence from edge direction, and refines predictions using a triangular edge refinement module that enables reasoning over higher\-order causal motifs\. Additionally,FoundCauseincorporates a dedicated confounder module based on learnable latent representations to model hidden common causes and produce a symmetric confounding matrix\. Together, these components yield a unified framework that bridges classical causal discovery principles and modern deep learning, enabling scalable, single\-pass inference of causal graphs with both directed and latent structure\.

### 1\.1Related Works

Classical Causal Discovery\.Classical methods recover graph structure via three paradigms\. Constraint\-based approaches \(e\.g\., PC, FCI, RFCI, GFCI\) use conditional independence tests and orientation rules to estimate Markov equivalence classes, with FCI variants handling latent confounding\(Spirteset al\.,[2000](https://arxiv.org/html/2606.17516#bib.bib1); Colombo and Maathuis,[2014](https://arxiv.org/html/2606.17516#bib.bib58); Ogarrioet al\.,[2016](https://arxiv.org/html/2606.17516#bib.bib57)\)\. Score\-based methods such as GES search over equivalence classes using decomposable scores\(Chickering,[2002](https://arxiv.org/html/2606.17516#bib.bib3); Nazaret and Blei,[2024](https://arxiv.org/html/2606.17516#bib.bib44)\)\. Continuous optimization approaches \(e\.g\., NOTEARS, DAG\-GNN, GraN\-DAG, GOLEM, DAGMA\) replace combinatorial search with smooth acyclicity constraints or differentiable DAG penalties\(Zhenget al\.,[2018](https://arxiv.org/html/2606.17516#bib.bib9); Yuet al\.,[2019](https://arxiv.org/html/2606.17516#bib.bib56); Lachapelleet al\.,[2020](https://arxiv.org/html/2606.17516#bib.bib55); Nget al\.,[2020](https://arxiv.org/html/2606.17516#bib.bib10)\)\. While well\-founded, these methods require solving a separate optimization for each dataset and are sensitive to independence tests, scoring functions, modeling assumptions, and hyperparameters\.

Extensions beyond Purely Observational DAGs\.Several approaches address limitations of observational discovery\. Intervention\-aware methods \(e\.g\., DCDI, ENCO\) leverage interventional data to improve identifiability, with ENCO explicitly separating edge existence and orientation\(Brouillardet al\.,[2020](https://arxiv.org/html/2606.17516#bib.bib53); Lippeet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib52)\)\. Bayesian and generative approaches \(e\.g\., DECI, DiBS, BayesDAG\) learn distributions over graphs and mechanisms for uncertainty\-aware inference, but require dataset\-specific posterior inference\(Geffneret al\.,[2024](https://arxiv.org/html/2606.17516#bib.bib51); Lorchet al\.,[2021](https://arxiv.org/html/2606.17516#bib.bib50)\)\. Bivariate methods \(e\.g\., ANM, IGCI, PNL\) exploit asymmetries such as residual independence, non\-Gaussianity, or score structure to infer direction from pairs\(Mooijet al\.,[2016](https://arxiv.org/html/2606.17516#bib.bib18); Rollandet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib7)\)\. Latent\-variable methods model hidden confounding via mixed graphical representations \(e\.g\., PAGs, MAGs, ADMGs\)\(Bhattacharyaet al\.,[2021](https://arxiv.org/html/2606.17516#bib.bib45); Ashmanet al\.,[2023](https://arxiv.org/html/2606.17516#bib.bib48); Maet al\.,[2024](https://arxiv.org/html/2606.17516#bib.bib49)\)\. In contrast, our approach embeds asymmetry and confounding signals directly into an amortized architecture, rather than relying on standalone heuristics or separate inference procedures\.

Identifiability Assumptions\.What each method can recover depends on its identifiability assumptions\. Under Markov and faithfulness, the DAG is identifiable only up to its Markov equivalence class\(Spirteset al\.,[2000](https://arxiv.org/html/2606.17516#bib.bib1)\)\. Full\-DAG identifiability requires additional structure, including linear non\-Gaussian noise \(LiNGAM\)\(Shimizuet al\.,[2006](https://arxiv.org/html/2606.17516#bib.bib4)\), nonlinear additive noise in the bivariate case \(ANMs\)\(Hoyeret al\.,[2009](https://arxiv.org/html/2606.17516#bib.bib62)\), post\-nonlinear models with invertible transformations\(Zhang and Hyvärinen,[2009](https://arxiv.org/html/2606.17516#bib.bib63)\), nonlinear additive models yielding a topological order \(CAM\)\(Bühlmannet al\.,[2013](https://arxiv.org/html/2606.17516#bib.bib64)\), equal or known error variances in linear Gaussian systems\(Peters and Bühlmann,[2014](https://arxiv.org/html/2606.17516#bib.bib65)\), or score\-Jacobian asymmetries in additive\-Gaussian DAGs\(Rollandet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib7)\)\. With latent confounding or selection bias, identification weakens to PAGs representing equivalence classes of MAGs\(Zhang,[2008](https://arxiv.org/html/2606.17516#bib.bib6); Colombo and Maathuis,[2014](https://arxiv.org/html/2606.17516#bib.bib58); Ogarrioet al\.,[2016](https://arxiv.org/html/2606.17516#bib.bib57)\), while counterfactual identification additionally requires monotonicity\(Chaoet al\.,[2023](https://arxiv.org/html/2606.17516#bib.bib66)\)or bijectivity\(Nasr\-Esfahanyet al\.,[2023](https://arxiv.org/html/2606.17516#bib.bib69)\)\. These guarantees are asymptotic and apply to narrow, largely non\-overlapping regimes that are difficult to verify in practice: faithfulness can fail in finite samples\(Uhleret al\.,[2013](https://arxiv.org/html/2606.17516#bib.bib67)\), and real\-world mechanisms rarely satisfy strict linearity, additive\-noise, or equal\-variance assumptions\(Glymouret al\.,[2019](https://arxiv.org/html/2606.17516#bib.bib68)\)\. Rather than committing to a single regime,FoundCauseamortizes inference across a broad family of synthetic SCMs spanning many of these settings, learning empirical patterns indicative of causal direction\.

Amortized Causal Discovery\.Amortized causal discovery trains a model once on synthetic SCMs and performs inference in a single forward pass\. AVICI uses alternating attention over variables and samples with a simple edge decoder\(Lorchet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib21)\); CSIvA employs an encoder–decoder with variable\-identity embeddings and autoregressive graph generation\(Keet al\.,[2023](https://arxiv.org/html/2606.17516#bib.bib38)\); and hybrid methods such as SEA aggregate graphs from classical algorithms on variable subsets via a neural aggregator\(Wuet al\.,[2025](https://arxiv.org/html/2606.17516#bib.bib42)\)\. Amortization has also been applied beyond structure discovery, including ATE/CATE and intervention\-effect estimation\(Nilforoshanet al\.,[2023](https://arxiv.org/html/2606.17516#bib.bib36); Robertsonet al\.,[2025](https://arxiv.org/html/2606.17516#bib.bib35); Balazadehet al\.,[2026](https://arxiv.org/html/2606.17516#bib.bib34); Reuteret al\.,[2026](https://arxiv.org/html/2606.17516#bib.bib41)\)\.FoundCausefollows this paradigm but differs in four ways: \(i\) it injects pairwise causal statistics into attention and decoding, \(ii\) it factorizes edge existence and direction, \(iii\) it refines pairwise representations via triangular graph\-level message passing, and \(iv\) it explicitly models latent confounding with dedicated tokens\. These choices bridge classical heuristics and end\-to\-end neural architectures, preserving scalability while introducing inductive biases for dependence, directionality, higher\-order structure, and hidden causes\. Table[4](https://arxiv.org/html/2606.17516#A1.T4)in Appendix[A](https://arxiv.org/html/2606.17516#A1)summarizes these distinctions\.

## 2Problem Setup

Given an observational dataset𝑿∈ℝN×D\\bm\{X\}\\in\\mathbb\{R\}^\{N\\times D\}consisting ofNNi\.i\.d\. samples overDDobserved variables, the goal is to infer both directed causal relations and latent confounding among the observed variables\. While classical and recent amortized approaches primarily focus on recovering directed structure over observed variables, explicitly modeling latent confounding remains comparatively underexplored in scalable settings\. We represent directed causal structure by a directed acyclic graph \(DAG\)𝑫∈\{0,1\}D×D\\bm\{D\}\\in\\\{0,1\\\}^\{D\\times D\}, whereDi​j=1D\_\{ij\}=1indicates that variableiiis a direct cause of variablejj\. Latent confounding is represented by a symmetric matrix𝑪∈\{0,1\}D×D\\bm\{C\}\\in\\\{0,1\\\}^\{D\\times D\}, whereCi​j=1C\_\{ij\}=1indicates that variablesiiandjjshare an unobserved common cause\. The data\-generating process is assumed to follow a structural causal model,Xj=fj​\(Pa​\(j\),ϵj\),X\_\{j\}=f\_\{j\}\\\!\\left\(\\mathrm\{Pa\}\(j\),\\epsilon\_\{j\}\\right\),j=1,…,D,j=1,\\ldots,D,wherePa​\(j\)\\mathrm\{Pa\}\(j\)denotes the parents of nodejjin the true DAG,fjf\_\{j\}is an arbitrary structural mechanism, potentially nonlinear and non\-additive, andϵj\\epsilon\_\{j\}is exogenous noise\. Some causes may be unobserved, thereby inducing statistical dependence among observed variables that is not explained by directed edges alone, motivating the need to jointly infer both directed structure and latent confounding\.

Rather than running a separate optimization procedure for each dataset, e\.g\., as in PCSpirteset al\.\([2000](https://arxiv.org/html/2606.17516#bib.bib1)\), GESChickering \([2002](https://arxiv.org/html/2606.17516#bib.bib3)\), or NOTEARSZhenget al\.\([2018](https://arxiv.org/html/2606.17516#bib.bib9)\), we use*amortized causal discovery*following the AVICI frameworkLorchet al\.\([2022](https://arxiv.org/html/2606.17516#bib.bib21)\)\. The model is trained on thousands of synthetic datasets with known ground\-truth structure and, at test time, performs causal discovery with a single forward pass\. This approach posits \(and verifies\) that statistical patterns learned from diverse synthetic data distributions transfer to real\-world datasets\.

The model takes as input a data matrix𝑿\\bm\{X\}together with an optional binary observation mask𝑴∈\{0,1\}N×D\\bm\{M\}\\in\\\{0,1\\\}^\{N\\times D\}, whereMn​i=1M\_\{ni\}=1denotes an observed entry andMn​i=0M\_\{ni\}=0denotes a missing entry\. During training, we sample problems withD∈\[2,50\]D\\in\[2,50\]variables andN∈\[100,600\]N\\in\[100,600\]observations; at inference time, the same architecture can be applied to larger numbers of variables and samples\. The output consists of an edge\-probability matrix𝑫^∈\[0,1\]D×D\\hat\{\\bm\{D\}\}\\in\[0,1\]^\{D\\times D\}and a symmetric confounding\-probability matrix𝑪^∈\[0,1\]D×D\\hat\{\\bm\{C\}\}\\in\[0,1\]^\{D\\times D\}\. For directed edges,D^i​j=σ​\(ℓi​j\)\\hat\{D\}\_\{ij\}=\\sigma\(\\ell\_\{ij\}\)is the predicted probability of the causal relationi→ji\\to j, whereℓi​j\\ell\_\{ij\}is the corresponding edge logit\. For latent confounding,C^i​j\\hat\{C\}\_\{ij\}is the predicted probability that variablesiiandjjshare a hidden common cause\.

## 3Model Architecture and Loss Function

Data𝑿∈ℝN×D\\bm\{X\}\\in\\mathbb\{R\}^\{N\\times D\}Mask𝑴∈\{0,1\}N×D\\bm\{M\}\\in\\\{0,1\\\}^\{N\\times D\}Per\-variable normalization \(fp32\)2\-channel input\[x¯n​i​mn​i,mn​i\]\[\\bar\{x\}\_\{ni\}m\_\{ni\},\\;m\_\{ni\}\]Pairwise Statistics45 raw\-data featuressym/asym/V\-structure\+\+reliability metadataAxis\-Factorized Encoderalternating variable\- and sample\-attention16 blocks,dh=768d\_\{h\}\{=\}768, 8 heads,Kg=12K\_\{g\}\{=\}12global tokensstat\-conditioned bias in variable\-attention blocksPMA Poolingover samplesKq=8K\_\{q\}\{=\}8query tokens, 4 heads, max\-pool residual𝒁∈ℝD×dh\\bm\{Z\}\\in\\mathbb\{R\}^\{D\\times d\_\{h\}\}Feature GateMLP​\(dh\+45→256→20\)\\mathrm\{MLP\}\(d\_\{h\}\{\+\}45\\to 256\\to 20\)per\-pair sigmoid gateConfounder ModuleKc=8K\_\{c\}\{=\}8tokens, 2\-layer cross\-attentionnoisy\-OR loadings→𝑪^\\to\\hat\{\\bm\{C\}\}Factored Edge Decoderexistence head: symmetric features\+\+detachedC^i​j\\hat\{C\}\_\{ij\}direction head: antisymmetric gated featuresTriangularRefinement3 rounds: outgoing, fork,collider, row\-wise attnSkel/Dir Blend\+\+parent–child role scorelearnedσ​\(αe\),σ​\(αd\)\\sigma\(\\alpha\_\{e\}\),\\sigma\(\\alpha\_\{d\}\); clamp±15\\pm 15; DAG post\-processing𝑫^\\hat\{\\bm\{D\}\}𝑪^\\hat\{\\bm\{C\}\}detachC^i​j\\hat\{C\}\_\{ij\}Figure 1:Top\-level Architecture\.The model takes normalized observations and pairwise statistics as input and processes them through four stages: an axis\-factorized encoder, attention\-based pooling \(PMA\) to obtain per\-variable embeddings, a confounder module that predicts latent confounding via noisy\-OR aggregation, and a factored edge decoder with triangular refinement to produce directed edge probabilities\. Dashed lines indicate auxiliary inputs from pairwise statistics / feature gate\.The architecture \(see Figure[1](https://arxiv.org/html/2606.17516#S3.F1)\) consists of four stages: \(i\) an axis\-factorized encoder, \(ii\) attention\-based pooling, \(iii\) a factored edge decoder with triangular refinement, and \(iv\) a learnable\-token confounder module\. The model has approximately139139M parameters; architectural hyperparameters are detailed in Appendix[E](https://arxiv.org/html/2606.17516#A5)\.

Input Representation\.Each entry\(n,i\)\(n,i\)is represented as𝒙n​i=\[x¯n​i​mn​i,mn​i\]∈ℝ2,\\bm\{x\}\_\{ni\}=\[\\bar\{x\}\_\{ni\}m\_\{ni\},\\,m\_\{ni\}\]\\in\\mathbb\{R\}^\{2\},wherex¯n​i\\bar\{x\}\_\{ni\}is a normalized value andmn​i∈\{0,1\}m\_\{ni\}\\in\\\{0,1\\\}indicates whether the entry is observed\. A linear projection produces𝒉n​i\(0\)∈ℝdh\\bm\{h\}^\{\(0\)\}\_\{ni\}\\in\\mathbb\{R\}^\{d\_\{h\}\}via𝒉n​i\(0\)=𝑾in​𝒙n​i\+𝒃in\\bm\{h\}^\{\(0\)\}\_\{ni\}=\\bm\{W\}\_\{\\mathrm\{in\}\}\\bm\{x\}\_\{ni\}\+\\bm\{b\}\_\{\\mathrm\{in\}\},𝑾in∈ℝdh×2\\bm\{W\}\_\{\\mathrm\{in\}\}\\in\\mathbb\{R\}^\{d\_\{h\}\\times 2\}\. No positional embeddings are used, making the model permutation\-invariant over variables, which is essential since causal structure should not depend on variable ordering\.

We now describe the key architectural details of these four stages mentioned above\.

\(i\) Axis\-factorized Encoder\.Let𝑯∈ℝN×D×dh\\bm\{H\}\\in\\mathbb\{R\}^\{N\\times D\\times d\_\{h\}\}denote the stacked embeddings\. We apply2​L2Lattention blocks alternating between variable\-wise attention \(overDD\) and sample\-wise attention \(overNN\) \(see Figure[3](https://arxiv.org/html/2606.17516#A5.F3), Appendix[E](https://arxiv.org/html/2606.17516#A5)\)\. This factorization reflects causal discovery, which relies both on cross\-variable dependencies and distributional properties such as conditional independence\.

Variable\-attention incorporates: \(i\)Kg∈ℕK\_\{g\}\\in\\mathbb\{N\}learnable global tokens appended along the variable dimension, and \(ii\) a stat\-conditioned biasbi​j\(h,ℓ\)=\[𝑾2\(ℓ\)​GELU⁡\(𝑾1\(ℓ\)​𝒔i​j\)\]h,b\_\{ij\}^\{\(h,\\ell\)\}=\\left\[\\bm\{W\}\_\{2\}^\{\(\\ell\)\}\\,\\operatorname\{GELU\}\\\!\\left\(\\bm\{W\}\_\{1\}^\{\(\\ell\)\}\\bm\{s\}\_\{ij\}\\right\)\\right\]\_\{h\},where𝒔i​j∈ℝds\\bm\{s\}\_\{ij\}\\in\\mathbb\{R\}^\{d\_\{s\}\}are pairwise statistics andGELU\\operatorname\{GELU\}is the Gaussian Error Linear UnitHendrycks and Gimpel \([2016](https://arxiv.org/html/2606.17516#bib.bib19)\)\. This bias injects classical dependence signals \(e\.g\., correlation and asymmetry\), guiding the model toward known causal cues\.

Each block uses RMSNormZhang and Sennrich \([2019](https://arxiv.org/html/2606.17516#bib.bib14)\), multi\-head scaled dot\-product attention \(SDPA\)Vaswaniet al\.\([2017](https://arxiv.org/html/2606.17516#bib.bib20)\), and a SwiGLUShazeer \([2020](https://arxiv.org/html/2606.17516#bib.bib15)\)feedforward network\. The encoder outputs contextualized representations𝑯enc∈ℝN×D×dh\\bm\{H\}^\{\\mathrm\{enc\}\}\\in\\mathbb\{R\}^\{N\\times D\\times d\_\{h\}\}\.

\(ii\) PMA Pooling over Samples\.For each variablei∈\{1,…,D\}i\\in\\\{1,\\ldots,D\\\}, we aggregate\{𝑯n​ienc\}n=1N\\\{\\bm\{H\}^\{\\mathrm\{enc\}\}\_\{ni\}\\\}\_\{n=1\}^\{N\}using Pooling by Multihead Attention \(PMA\)Leeet al\.\([2019](https://arxiv.org/html/2606.17516#bib.bib13)\), producing𝒛iPMA∈ℝdh\\bm\{z\}\_\{i\}^\{\\mathrm\{PMA\}\}\\in\\mathbb\{R\}^\{d\_\{h\}\}\. We also compute𝒛imax=maxn=1,…,N⁡𝑯n​ienc\.\\bm\{z\}\_\{i\}^\{\\max\}=\\max\_\{n=1,\\ldots,N\}\\bm\{H\}^\{\\mathrm\{enc\}\}\_\{ni\}\.These are combined via a learned gateg∈ℝg\\in\\mathbb\{R\}:𝒛i=σ​\(g\)​𝒛iPMA\+\(1−σ​\(g\)\)​𝒛imax,\\bm\{z\}\_\{i\}=\\sigma\(g\)\\bm\{z\}\_\{i\}^\{\\mathrm\{PMA\}\}\+\(1\-\\sigma\(g\)\)\\bm\{z\}\_\{i\}^\{\\max\},whereσ​\(\)\\sigma\(\)denotes the sigmoid function\. Stacking over variables yields𝒁∈ℝD×dh\\bm\{Z\}\\in\\mathbb\{R\}^\{D\\times d\_\{h\}\}\. This pooling captures higher\-order distributional properties \(e\.g\., multimodality and tail behavior\)\.

\(iii\) Factored Edge Decoder\.For each ordered pair\(i,j\)\(i,j\), we construct𝒇i​jexist=\[𝒛i\+𝒛j;𝒛i⊙𝒛j;C^i​j\],\\bm\{f\}^\{\\mathrm\{exist\}\}\_\{ij\}=\[\\bm\{z\}\_\{i\}\+\\bm\{z\}\_\{j\};\\bm\{z\}\_\{i\}\\odot\\bm\{z\}\_\{j\};\\hat\{C\}\_\{ij\}\],and𝒇i​jdir=\[𝒛i−𝒛j\]\.\\bm\{f\}^\{\\mathrm\{dir\}\}\_\{ij\}=\[\\bm\{z\}\_\{i\}\-\\bm\{z\}\_\{j\}\]\.This separates edge existence \(symmetric dependence\) from direction \(asymmetric signals\), mirroring classical causal discovery principles\.

A shared MLP producesℓi​jbase∈ℝ\\ell^\{\\mathrm\{base\}\}\_\{ij\}\\in\\mathbb\{R\}\. We then add a role scoreri​j=\(𝑾V​𝒛i\)⊤​\(𝑾U​𝒛j\)/dpair,r\_\{ij\}=\{\(\\bm\{W\}\_\{V\}\\bm\{z\}\_\{i\}\)^\{\\top\}\(\\bm\{W\}\_\{U\}\\bm\{z\}\_\{j\}\)\}/\{\\sqrt\{d\_\{\\mathrm\{pair\}\}\}\},where𝑾V,𝑾U∈ℝdpair×dh\\bm\{W\}\_\{V\},\\bm\{W\}\_\{U\}\\in\\mathbb\{R\}^\{d\_\{\\mathrm\{pair\}\}\\times d\_\{h\}\}, yieldingℓi​j=ℓi​jbase\+α​ri​j\\ell\_\{ij\}=\\ell^\{\\mathrm\{base\}\}\_\{ij\}\+\\alpha r\_\{ij\},α∈ℝ\\alpha\\in\\mathbb\{R\}\. This provides an explicit inductive bias toward consistent parent–child roles\. See also Figure[4](https://arxiv.org/html/2606.17516#A5.F4), Appendix[E](https://arxiv.org/html/2606.17516#A5)\.

Triangular Refinement\.We construct pair representations𝑷∈ℝD×D×dpair\\bm\{P\}\\in\\mathbb\{R\}^\{D\\times D\\times d\_\{\\mathrm\{pair\}\}\}and refine them forRRrounds using triangle\-based updates inspired by AlphaFold2Jumperet al\.\([2021](https://arxiv.org/html/2606.17516#bib.bib12)\)\. Each update aggregates over intermediate variablesk∈\{1,…,D\}k\\in\\\{1,\\ldots,D\\\}, enabling reasoning over triples\(i,k,j\)\(i,k,j\)\. This captures causal motifs such as chains, forks, and colliders, which are helpful for distinguishing Markov\-equivalent structures\. Refined logitsℓtri∈ℝD×D\\bm\{\\ell\}^\{\\mathrm\{tri\}\}\\in\\mathbb\{R\}^\{D\\times D\}are combined withℓ\\bm\{\\ell\}, and final edge probabilities areD^i​j=σ​\(ℓi​j\)\.\\hat\{D\}\_\{ij\}=\\sigma\(\\ell\_\{ij\}\)\.

\(iv\) Confounder Module\.Latent confounding is modeled using learnable matrix𝑻∈ℝKc×dh\\bm\{T\}\\in\\mathbb\{R\}^\{K\_\{c\}\\times d\_\{h\}\}, where each row𝒕k∈ℝdh\\bm\{t\}\_\{k\}\\in\\mathbb\{R\}^\{d\_\{h\}\}represents a candidate hidden common cause andKcK\_\{c\}is the number of such confounders \(see Figure[5](https://arxiv.org/html/2606.17516#A5.F5), Appendix[E](https://arxiv.org/html/2606.17516#A5)\)\. For each variableiiand confounderkk, we computeSi​k=σ​\(MLP​\(\[𝒛i;𝒕k\]\)\)S\_\{ik\}=\\sigma\(\\mathrm\{MLP\}\(\[\\bm\{z\}\_\{i\};\\bm\{t\}\_\{k\}\]\)\),Si​k∈\(0,1\),S\_\{ik\}\\in\(0,1\),which denotes the probability that confounderkkinfluences variableii\. Pairwise confounding is computed via a noisy\-OR aggregation:C^i​j=1−∏k=1Kc\(1−Si​k​Sj​k\),\\hat\{C\}\_\{ij\}=1\-\\prod\_\{k=1\}^\{K\_\{c\}\}\(1\-S\_\{ik\}S\_\{jk\}\),yielding a symmetric matrix𝑪^∈\[0,1\]D×D\\hat\{\\bm\{C\}\}\\in\[0,1\]^\{D\\times D\}\. This explicitly models hidden common causes, allowing the model to distinguish direct causal edges from correlations induced by latent confounding\.

Final Loss Function\.Letℓ∈ℝD×D\\bm\{\\ell\}\\in\\mathbb\{R\}^\{D\\times D\}denote the directed edge logits, whereℓi​j\\ell\_\{ij\}corresponds to the predicted logit for edgei→ji\\to j\. Let𝑫∗∈\{0,1\}D×D\\bm\{D\}^\{\*\}\\in\\\{0,1\\\}^\{D\\times D\}denote the ground\-truth adjacency matrix of the DAG, and𝑪∗∈\{0,1\}D×D\\bm\{C\}^\{\*\}\\in\\\{0,1\\\}^\{D\\times D\}the ground\-truth confounding matrix\. We define a valid\-pair maskVi​j=𝟙​\[i≠j\]​mi​mj,V\_\{ij\}=\\mathbbm\{1\}\[i\\neq j\]\\,m\_\{i\}m\_\{j\},which removes diagonal entries and padded variables, wheremi∈\{0,1\}m\_\{i\}\\in\\\{0,1\\\}indicates whether variableiiis present\. Let\|V\|=∑i​jVi​j\|V\|=\\sum\_\{ij\}V\_\{ij\}, and define the number of positive and negative edges asn\+=∑i​jDi​j∗​Vi​jn\_\{\+\}=\\sum\_\{ij\}D^\{\*\}\_\{ij\}V\_\{ij\}andn−=∑i​j\(1−Di​j∗\)​Vi​jn\_\{\-\}=\\sum\_\{ij\}\(1\-D^\{\*\}\_\{ij\}\)V\_\{ij\}respectively\. The overall objective:

ℒ=ℒdir\+λasym​ℒasym\+λcal​ℒcal\+λbdir​ℒbdir\+λskel​ℒskel\+λconf​ℒconf\+ℒfp\+λden​ℒden\+λpma​ℒpma\+λgate​ℒgate\\mathcal\{L\}=\\mathcal\{L\}\_\{\\mathrm\{dir\}\}\+\\lambda\_\{\\mathrm\{asym\}\}\\mathcal\{L\}\_\{\\mathrm\{asym\}\}\+\\lambda\_\{\\mathrm\{cal\}\}\\mathcal\{L\}\_\{\\mathrm\{cal\}\}\+\\lambda\_\{\\mathrm\{bdir\}\}\\mathcal\{L\}\_\{\\mathrm\{bdir\}\}\+\\lambda\_\{\\mathrm\{skel\}\}\\mathcal\{L\}\_\{\\mathrm\{skel\}\}\+\\lambda\_\{\\mathrm\{conf\}\}\\mathcal\{L\}\_\{\\mathrm\{conf\}\}\\\\ \+\\mathcal\{L\}\_\{\\mathrm\{fp\}\}\+\\lambda\_\{\\mathrm\{den\}\}\\mathcal\{L\}\_\{\\mathrm\{den\}\}\+\\lambda\_\{\\mathrm\{pma\}\}\\mathcal\{L\}\_\{\\mathrm\{pma\}\}\+\\lambda\_\{\\mathrm\{gate\}\}\\mathcal\{L\}\_\{\\mathrm\{gate\}\}\(1\)which combines directed edge prediction, direction calibration, skeleton supervision, confounding prediction, and regularization\. We first describe the directed edge loss, and then briefly summarize the remaining terms\.

\(a\) Directed Edge Loss\.The main term is a class\-balanced binary cross\-entropy over directed edges\. We use direction\-aware label smoothing

D~i​j=\{1−ϵd,Di​j∗=1,ϵd,Dj​i∗=1,ϵ0,otherwise,ϵd=0\.05,ϵ0=0\.005,\\tilde\{D\}\_\{ij\}=\\begin\{cases\}1\-\\epsilon\_\{d\},&D^\{\*\}\_\{ij\}=1,\\\\ \\epsilon\_\{d\},&D^\{\*\}\_\{ji\}=1,\\\\ \\epsilon\_\{0\},&\\text\{otherwise\},\\end\{cases\}\\qquad\\epsilon\_\{d\}=0\.05,\\quad\\epsilon\_\{0\}=0\.005,which introduces asymmetric supervision between\(i,j\)\(i,j\)and\(j,i\)\(j,i\), encouraging consistent edge orientation\.

We define

ℒdir=1\|V\|​∑i​j\(Di​j∗​w\+\+\(1−Di​j∗\)\)​Vi​j​BCE​\(ℓi​j,D~i​j\),\\mathcal\{L\}\_\{\\mathrm\{dir\}\}=\\frac\{1\}\{\|V\|\}\\sum\_\{ij\}\\Bigl\(D^\{\*\}\_\{ij\}w\_\{\+\}\+\(1\-D^\{\*\}\_\{ij\}\)\\Bigr\)V\_\{ij\}\\,\\mathrm\{BCE\}\(\\ell\_\{ij\},\\tilde\{D\}\_\{ij\}\),whereBCE​\(ℓ,y\)=−y​log⁡σ​\(ℓ\)−\(1−y\)​log⁡\(1−σ​\(ℓ\)\)\.\\mathrm\{BCE\}\(\\ell,y\)=\-y\\log\\sigma\(\\ell\)\-\(1\-y\)\\log\(1\-\\sigma\(\\ell\)\)\.To address class imbalance, the positive\-class weightw\+w\_\{\+\}is adapted online based on recall:r^t=0\.95​r^t−1\+0\.05​rt,\\hat\{r\}\_\{t\}=0\.95\\hat\{r\}\_\{t\-1\}\+0\.05r\_\{t\},and

w\+=clamp​\(n−n\+​\[1\+2​clamp​\(0\.65−r^t,−0\.3,0\.3\)\],1,5\),w\_\{\+\}=\\mathrm\{clamp\}\\\!\\left\(\\sqrt\{\\frac\{n\_\{\-\}\}\{n\_\{\+\}\}\}\\Bigl\[1\+2\\,\\mathrm\{clamp\}\(0\.65\-\\hat\{r\}\_\{t\},\-0\.3,0\.3\)\\Bigr\],1,5\\right\),wherertr\_\{t\}denotes the recall at iterationtt\. This stabilizes learning under severe edge sparsity while adapting the recall–precision tradeoff\.

\(b\) Auxiliary Losses\.The remaining terms address failure modes not captured by pairwise BCE\.ℒasym\\mathcal\{L\}\_\{\\mathrm\{asym\}\}suppresses the reverse direction of true edges when its confidence exceeds a threshold, reinforcing causal asymmetry\.ℒcal\\mathcal\{L\}\_\{\\mathrm\{cal\}\}limits extreme logit gaps between opposite directions, improving calibration and robustness under distribution shift\.ℒbdir\\mathcal\{L\}\_\{\\mathrm\{bdir\}\}applies a direction loss only to edges whose existence has been detected, decoupling edge discovery from orientation\.ℒskel\\mathcal\{L\}\_\{\\mathrm\{skel\}\}provides supervision on the undirected skeleton, aligning with the fact that many causal signals are identifiable only up to equivalence classes\.ℒconf\\mathcal\{L\}\_\{\\mathrm\{conf\}\}supervises the confounder module, enabling the model to distinguish direct edges from dependencies induced by latent common causes\.

Finally,ℒfp\\mathcal\{L\}\_\{\\mathrm\{fp\}\}andℒden\\mathcal\{L\}\_\{\\mathrm\{den\}\}penalize overly dense graphs, whileℒpma\\mathcal\{L\}\_\{\\mathrm\{pma\}\}andℒgate\\mathcal\{L\}\_\{\\mathrm\{gate\}\}regularize the pooling and gating mechanisms\. Full definitions and hyperparameters are provided in Appendix[F](https://arxiv.org/html/2606.17516#A6)\.

## 4Training Pipeline

Optimization\.We train using Schedule\-Free AdamWDefazioet al\.\([2024](https://arxiv.org/html/2606.17516#bib.bib16)\)with base learning rate4×10−44\\times 10^\{\-4\}, cosine decay toηmin=10−5\\eta\_\{\\min\}=10^\{\-5\}over 300 epochs followed by a constant schedule, and momentum parametersβ1=0\.9\\beta\_\{1\}=0\.9,β2=0\.99\\beta\_\{2\}=0\.99\. Weight decay of10−410^\{\-4\}is applied to all weight matrices \(excluding biases and normalization parameters\), and gradients are clipped to a global norm of 1\.0\. Training runs for 1000 epochs with batch size 24\. Schedule\-Free AdamW maintains training parameters𝒛\\bm\{z\}and evaluation parameters𝒙\\bm\{x\}, where𝒙\\bm\{x\}is a running average of𝒛\\bm\{z\}; updates are applied to𝒛\\bm\{z\}, while evaluation and checkpointing use𝒙\\bm\{x\}\.

Training Setup\.We train using bfloat16 mixed precision with selective float32 computation for numerically sensitive operations \(e\.g\., normalization, matrix inversions, and loss computation\)\. Training uses distributed data parallelism across 8 A100 GPUs with synchronized gradients and data\-parallel sharding; global statistics for adaptive loss weighting are synchronized across workers\.

Data Regeneration and Training Strategy\.Training data are continuously regenerated via parallel sampling of synthetic SCMs, exposing the model to a stream of diverse datasets rather than a fixed corpus\. To prevent shortcut learning, we apply permutation augmentation—randomly permuting variable indices with corresponding relabeling of𝑫∗\\bm\{D\}^\{\*\}and𝑪∗\\bm\{C\}^\{\*\}—followingReisachet al\.\([2023](https://arxiv.org/html/2606.17516#bib.bib17)\), removing ordering artifacts and encouraging reliance on distributional and relational signals\. We further employ hard mining by replaying high\-loss examples from recent epochs, maintaining focus on challenging cases\.

## 5Experimental Results

Training data are generated entirely from synthetic structural causal models \(SCMs\) using DoWhy\(Sharma and Kiciman,[2020](https://arxiv.org/html/2606.17516#bib.bib23); Blöbaumet al\.,[2024](https://arxiv.org/html/2606.17516#bib.bib33)\), by sampling diverse graphs and mechanisms and drawing observational datasets with full supervision for directed edges and latent confounding\. The model is trained exclusively on synthetic data, without using real\-world datasets\. Additional details are provided in Appendix[G](https://arxiv.org/html/2606.17516#A7)\.FoundCauseperforms inference in under 2 seconds on average across all datasets, whereas classical methods often require per\-dataset optimization or combinatorial search and can take hours on moderate\-sized graphs\(Spirteset al\.,[2000](https://arxiv.org/html/2606.17516#bib.bib1)\)\.

MethodCov\.Avg\. AUROC↑\\uparrowAvg\.F1F\_\{1\}↑\\uparrowAvg\. Skel\-F1F\_\{1\}↑\\uparrowAvg\. SHD/d/d↓\\downarrow*Classical baselines*PC \(α=0\.01\\alpha\{=\}0\.01\)\(Spirteset al\.,[2000](https://arxiv.org/html/2606.17516#bib.bib1)\)13/1513/150\.7100\.7100\.4640\.4640\.6420\.6421\.3831\.383PC \(α=0\.05\\alpha\{=\}0\.05\)\(Spirteset al\.,[2000](https://arxiv.org/html/2606.17516#bib.bib1)\)13/1513/150\.6920\.6920\.4300\.4300\.6460\.6461\.4951\.495FCI \(α=0\.05\\alpha\{=\}0\.05\)\(Spirteset al\.,[2000](https://arxiv.org/html/2606.17516#bib.bib1)\)13/1513/150\.6370\.6370\.3800\.3800\.4830\.4831\.2091\.209GES\(Chickering,[2002](https://arxiv.org/html/2606.17516#bib.bib3)\)15/1515/150\.7250\.7250\.4560\.4560\.6310\.6311\.5961\.596GRaSP\(Lamet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib8)\)15/1515/150\.7470\.7470\.4880\.4880\.676\\bm\{0\.676\}1\.4541\.454DirectLiNGAM\(Shimizuet al\.,[2011](https://arxiv.org/html/2606.17516#bib.bib2)\)15/1515/150\.6750\.6750\.3140\.3140\.5100\.5102\.6812\.681ICA\-LiNGAM\(Shimizuet al\.,[2006](https://arxiv.org/html/2606.17516#bib.bib4)\)15/1515/150\.6530\.6530\.2920\.2920\.5000\.5003\.0083\.008NOTEARS\-Linear\(Zhenget al\.,[2018](https://arxiv.org/html/2606.17516#bib.bib9)\)15/1515/150\.5880\.5880\.2510\.2510\.5090\.5091\.4701\.470DAGMA \(linear\)\(Belloet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib5)\)15/1515/150\.5970\.5970\.2670\.2670\.5400\.5401\.5461\.546SCORE\(Rollandet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib7)\)15/1515/150\.5540\.5540\.1430\.1430\.3660\.3664\.5384\.538Naive correlation15/1515/150\.5840\.5840\.1900\.1900\.4270\.4274\.6534\.653*Amortized / foundation\-model methods*AVICI\(Lorchet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib21)\)15/1515/150\.7320\.7320\.4980\.4980\.6510\.6511\.1421\.142CSIvA\(Keet al\.,[2023](https://arxiv.org/html/2606.17516#bib.bib38)\)15/1515/150\.7390\.7390\.5070\.5070\.6580\.6581\.1081\.108SEA\(Wuet al\.,[2025](https://arxiv.org/html/2606.17516#bib.bib42)\)15/1515/150\.7440\.7440\.5120\.5120\.6670\.6671\.0521\.052Cond\_FIP\(Mahajanet al\.,[2025](https://arxiv.org/html/2606.17516#bib.bib40)\)15/1515/150\.7280\.7280\.4890\.4890\.6460\.6461\.1761\.176FoundCause\(ours\)15/1515/150\.756\\bm\{0\.756\}0\.535\\bm\{0\.535\}0\.6630\.6630\.980\\bm\{0\.980\}Table 1:Average performance of various causal discovery techniques across 15 real\-world benchmark datasets\. All metrics are macro\-averaged across datasets\. AUROC andF1F\_\{1\}are higher\-is\-better, while normalized structural Hamming distance \(SHD/d\\mathrm\{SHD\}/d\) is lower\-is\-better; Skel\-F1F\_\{1\}measures undirected edge recovery \(higher is better\)\.*Cov\.*denotes the number of datasets on which a method successfully completed \(some classical methods fail or time out on dense, near\-collinearpetshopgraphs atd=41d=41\)\.FoundCauseis the only method that achieves the top rank onF1F\_\{1\},SHD/d\\mathrm\{SHD\}/d, and AUROC while successfully running on all15/1515/15datasets\. Representative per\-dataset results are in Appendix[H](https://arxiv.org/html/2606.17516#A8)\.### 5\.1Real\-World Benchmark Evaluation

Datasets\.We evaluateFoundCauseon 15 real\-world and semi\-realistic benchmarks spanning Bayesian networks, nonlinear variants, biological systems, and dense high\-dimensional graphs\. The suite includesasia,child,insurance,sachs,causal\_chambers\_lt,ecoli\_like, and fourpetshopvariants, with graph sizes up tod=41d=41\. We also evaluate on the Tübingen cause–effect pairs benchmark\(Mooijet al\.,[2016](https://arxiv.org/html/2606.17516#bib.bib18)\), a collection of real\-world bivariate tasks where the goal is to infer whetherX→YX\\to YorY→XY\\to X\. This setting is challenging because conditional\-independence tests and multi\-variable constraints are unavailable, requiring reliance on asymmetric distributional signals such as nonlinearity, non\-Gaussianity, and noise–mechanism independence\. Additional dataset details are provided in Appendix[C](https://arxiv.org/html/2606.17516#A3)\.

Compared Approaches\.While many causal discovery methods have been proposed, our goal is not exhaustive coverage but a comparison against a representative set of strong baselines commonly used in practice\. Accordingly, we compare against classical approaches spanning constraint\-based methods \(e\.g\., PCSpirteset al\.\([2000](https://arxiv.org/html/2606.17516#bib.bib1)\)\), score\-based search \(e\.g\., GESChickering \([2002](https://arxiv.org/html/2606.17516#bib.bib3)\)\), LiNGAM variantsShimizuet al\.\([2011](https://arxiv.org/html/2606.17516#bib.bib2)\), continuous DAG optimization methods \(e\.g\., NOTEARSZhenget al\.\([2018](https://arxiv.org/html/2606.17516#bib.bib9)\)\), along with additional baselines such as SCORE\(Rollandet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib7)\)and naive correlation\. We also include amortized or foundation\-style baselines, including AVICILorchet al\.\([2022](https://arxiv.org/html/2606.17516#bib.bib21)\), CSIvAKeet al\.\([2023](https://arxiv.org/html/2606.17516#bib.bib38)\), SEAWuet al\.\([2025](https://arxiv.org/html/2606.17516#bib.bib42)\), and Cond\_FIPMahajanet al\.\([2025](https://arxiv.org/html/2606.17516#bib.bib40)\)\. Constraint\-based methods \(PC, FCI, GES, GRaSP\) rely on conditional\-independence tests and cannot orient bivariate pairs without auxiliary assumptions, so they are excluded from the comparison on the Tübingen cause–effect pairs benchmark\. We also compare to SLOPPY\(Marx and Vreeken,[2019](https://arxiv.org/html/2606.17516#bib.bib47)\)and HECIXuet al\.\([2022](https://arxiv.org/html/2606.17516#bib.bib46)\), two state\-of\-the\-art methods for bivariate settings\.

Results\.Our inference protocol is described in Appendix[D](https://arxiv.org/html/2606.17516#A4)\. Table[1](https://arxiv.org/html/2606.17516#S5.T1)reports AUROC, directed\-edgeF1F\_\{1\}, skeletonF1F\_\{1\}, normalized structural Hamming distance \(SHD/d\\mathrm\{SHD\}/d\), and coverage\. AUROC measures edge\-ranking quality, directedF1F\_\{1\}evaluates recovery of oriented causal edges, skeletonF1F\_\{1\}measures adjacency recovery ignoring direction, andSHD/d\\mathrm\{SHD\}/dcaptures normalized graph\-edit error \(lower is better\); see Appendix[B](https://arxiv.org/html/2606.17516#A2)for precise definitions\. Coverage indicates whether a method completes on each dataset, which is important as some classical methods fail or time out on dense, near\-collinearpetshopgraphs \(e\.g\., PC and FCI complete13/1513/15datasets\), whereas all amortized methods andFoundCausecomplete all1515\.FoundCauseis evaluated from a single fixed checkpoint without retraining or per\-dataset tuning, using permutation averaging, temperature calibration, adaptive thresholding, and self\-consistency pruning\.

MethodCorrect / TotalUnweighted↑\\uparrowWeighted↑\\uparrowNOTEARS \(Linear\)\(Zhenget al\.,[2018](https://arxiv.org/html/2606.17516#bib.bib9)\)26/10226/10225\.5%25\.5\\%20\.1%20\.1\\%DAGMA \(linear\)\(Belloet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib5)\)34/10234/10233\.3%33\.3\\%34\.2%34\.2\\%Naive correlation42/10242/10241\.2%41\.2\\%44\.1%44\.1\\%DirectLiNGAM\(Shimizuet al\.,[2011](https://arxiv.org/html/2606.17516#bib.bib2)\)45/10245/10244\.1%44\.1\\%46\.5%46\.5\\%ICA\-LiNGAM\(Shimizuet al\.,[2006](https://arxiv.org/html/2606.17516#bib.bib4)\)52/10252/10251\.0%51\.0\\%49\.4%49\.4\\%AVICI\(Lorchet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib21)\)55/10255/10253\.9%53\.9\\%55\.8%55\.8\\%CSIvA\(Keet al\.,[2023](https://arxiv.org/html/2606.17516#bib.bib38)\)57/10257/10255\.9%55\.9\\%57\.4%57\.4\\%SEA\(Wuet al\.,[2025](https://arxiv.org/html/2606.17516#bib.bib42)\)59/10259/10257\.8%57\.8\\%58\.6%58\.6\\%Cond\_FIP\(Mahajanet al\.,[2025](https://arxiv.org/html/2606.17516#bib.bib40)\)54/10254/10252\.9%52\.9\\%54\.1%54\.1\\%HECI\(Xuet al\.,[2022](https://arxiv.org/html/2606.17516#bib.bib46)\)72/10272/10270\.4%70\.4\\%71\.8%71\.8\\%SLOPPY\(Marx and Vreeken,[2019](https://arxiv.org/html/2606.17516#bib.bib47)\)𝟕𝟒/𝟏𝟎𝟐\\bm\{74/102\}72\.5%\\bm\{72\.5\}\\%73\.4%\\bm\{73\.4\}\\%FoundCause\(ours\)65/102\{65/102\}63\.7%63\.7\\%65\.4%\{65\.4\\%\}Table 2:Results on the Tübingen cause–effect pairs benchmark\. We evaluate on 102 bivariate pairs with known causal direction and positive evaluation weights\. Unweighted accuracy is the fraction of pairs with correctly predicted direction, while weighted accuracy accounts for pair\-specific difficulty and relevance as defined by the benchmark\.FoundCauseis the only amortized method that approaches the performance of SLOPPY and HECI, which are specifically designed for the bivariate setting\.Across the datasets,FoundCauseachieves the best average AUROC, directedF1F\_\{1\}, andSHD/d\\mathrm\{SHD\}/dwhile maintaining full15/1515/15coverage\. It attains AUROC0\.7560\.756, outperforming the strongest classical baseline \(GRaSP,0\.7470\.747\) and amortized baseline \(SEA,0\.7440\.744\)\. More importantly, it improves directed\-edge recovery withF1=0\.535F\_\{1\}=0\.535versus0\.5120\.512\(SEA\) and0\.4880\.488\(GRaSP\), and reduces structural error toSHD/d=0\.980\\mathrm\{SHD\}/d=0\.980versus1\.0521\.052\(SEA\) and1\.4541\.454\(GRaSP\)\. While GRaSP achieves the highest skeletonF1F\_\{1\}\(0\.6760\.676\),FoundCauseremains close \(0\.6630\.663\) while substantially outperforming it on direction\-sensitive metrics, indicating stronger recovery of directed causal structure\.

On the Tübingen cause–effect pairs benchmark \(Table[2](https://arxiv.org/html/2606.17516#S5.T2)\),FoundCauseoutperforms the strongest classical baseline \(ICA\-LiNGAM\) by at least12\.712\.7percentage points in accuracy and is the only amortized method approaching the performance of SLOPPY and HECI, which are specialized for the bivariate setting\.

Test\-time ablations\.We perform test\-time ablations by replacing intermediate activations with neutral values at inference, without retraining, isolating each component’s contribution\. Pairwise statistics and the dedicated direction pathway are the primary drivers of performance, with large drops when removed, while triangular refinement provides additional gains via higher\-order reasoning\. Pooling and confounder feedback yield smaller improvements, indicating that performance is mainly driven by explicit causal inductive biases \(Appendix[I](https://arxiv.org/html/2606.17516#A9)\)\.

Robustness to Missing Data\.In Appendix[J](https://arxiv.org/html/2606.17516#A10),FoundCauseremains robust under Missing At Random \(MAR\) corruption, with minimal degradation at10%10\\%missingness and retaining92%92\\%of its clean\-dataF1F\_\{1\}at30%30\\%\. Imputation\-based baselines degrade substantially, highlighting the benefit of explicit mask\-aware modeling\.

### 5\.2Dimension Generalization

A key property ofFoundCauseis the absence of positional embeddings along the variable axis: variables are represented solely through observed samples and induced statistical relationships, making the model permutation\-invariant and agnostic to variable order and count\. The checkpoint is trained on graphs withD∈\[2,50\]D\\in\[2,50\]\. To evaluate generalization, we consider synthetic SCMs withD∈\{50,60,70,80,90,100\}D\\in\\\{50,60,70,80,90,100\\\}, whereD\>50D\>50corresponds to zero\-shot extrapolation\. For each dimension, we generate2020datasets using DoWhy withN=600N=600samples, varying root fractions in\[0\.05,0\.40\]\[0\.05,0\.40\], latent\-confounder fractions in\[0,0\.30\]\[0,0\.30\], and mechanism types from the same mixture as training \(linear, nonlinear additive, and non\-additive neural mechanisms with heterogeneous noise\)\. All evaluations use a single fixed checkpoint \(Table[1](https://arxiv.org/html/2606.17516#S5.T1)\) without retraining, fine\-tuning, or dimension\-specific calibration\.

We report macro\-averaged AUROC, AUPRC, and directed\-edgeF1F\_\{1\}over the2020tasks per dimension\. AUROC evaluates threshold\-independent ranking, AUPRC emphasizes performance under sparsity, andF1F\_\{1\}measures the final thresholded graph\. The goal is not constant performance beyond the training range, but to assess whether degradation is graceful and representations remain effective asDDincreases\.

5050606070708080909010010000\.20\.20\.40\.40\.60\.60\.80\.811Number of variablesDDAUROCF1F\_\{1\}AUPRCFigure 2:Performance vs\. dimension\. Macro\-averaged AUROC, AUPRC, andF1F\_\{1\}ofFoundCauseas a function of the number of variablesDD, evaluated on2020synthetic DoWhy datasets per dimension\. The model is trained onD∈\[2,50\]D\\in\[2,50\], and results forD≥60D\\geq 60represent zero\-shot extrapolation\. AUROC degrades gradually with increasingDD, whileF1F\_\{1\}declines more rapidly as the fixed decision threshold—calibrated on in\-distribution logits—becomes suboptimal at larger dimensions\.Figure[2](https://arxiv.org/html/2606.17516#S5.F2)shows thatFoundCauseretains strong edge\-ranking ability well beyond its training range: atD=100D=100\(twice the training ceiling\), AUROC remains0\.7120\.712, indicating many true edges are still ranked above non\-edges, whileF1F\_\{1\}drops from0\.5540\.554atD=50D=50to0\.2780\.278\. This gap suggests that extrapolation errors are driven by miscalibration rather than degraded representations\. AsDDincreases, the number of candidate pairs grows asD​\(D−1\)D\(D\-1\), shifting the logit distribution away from the regime for which the decision threshold is calibrated\. Thus, relative edge scores remain informative, but the fixed threshold becomes increasingly mismatched to the true graph density\.

Real Data Experiments\.We further evaluate this behavior on four high\-dimensional real\-world datasets withD∈\[56,100\]D\\in\[56,100\], using the same fixed checkpoint and inference pipeline without retraining or recalibration\.

DatasetDDAUROC↑\\uparrowF1↑F\_\{1\}\\uparrowSkel\-F1↑F\_\{1\}\\uparrowSHD/d↓/d\\downarrowHailfinder\(Abramsonet al\.,[1996](https://arxiv.org/html/2606.17516#bib.bib30)\)56560\.7420\.7420\.5160\.5160\.6030\.6031\.221\.22HEPAR II\(Onisko,[2003](https://arxiv.org/html/2606.17516#bib.bib25)\)70700\.7150\.7150\.4210\.4210\.5810\.5811\.581\.58WIN95PT\(Scutari,[2022](https://arxiv.org/html/2606.17516#bib.bib29)\)76760\.7890\.7890\.4760\.4760\.5580\.5581\.691\.69DREAM4\-100\(Marbachet al\.,[2009](https://arxiv.org/html/2606.17516#bib.bib24),[2010](https://arxiv.org/html/2606.17516#bib.bib32)\)1001000\.6610\.6610\.3410\.3410\.5220\.5222\.042\.04Table 3:Dimension\-generalization ofFoundCauseon additional real\-world datasets\.Table[3](https://arxiv.org/html/2606.17516#S5.T3)provides a real\-world stress test of dimension generalization beyond the synthetic DoWhy setting\. Despite variation in domain, graph density, and data\-generating mechanisms,FoundCauseretains meaningful performance without retraining or dataset\-specific calibration, with AUROC ranging from0\.6610\.661to0\.7890\.789\. Consistent with synthetic results, thresholded metrics degrade with increasing dimension:F1F\_\{1\}decreases from0\.5160\.516onhailfinderto0\.3410\.341ondream4\_100, whileSHD/d\\mathrm\{SHD\}/dincreases from1\.221\.22to2\.042\.04\. In contrast, skeleton scores remain relatively stable, suggesting reliable adjacency detection but reduced accuracy in edge orientation and calibration as the number of candidate pairs grows\. Overall, these results reinforce graceful zero\-shot generalization: ranking quality is preserved better than thresholded accuracy asDDincreases\. Section[K](https://arxiv.org/html/2606.17516#A11)compares SEA and CSIvA;FoundCauseremains within0\.020\.02–0\.050\.05inF1F\_\{1\}of the best method on three of four datasets and achieves the highest AUROC onwin95pts\.

## 6Conclusion

We presentedFoundCause, an amortized model for causal discovery from observational data that jointly predicts a directed acyclic graph and a latent confounding matrix in a single forward pass\. The architecture combines a permutation\-invariant transformer encoder with explicit statistical inductive biases, including pairwise statistics, a factored edge predictor, triangular refinement for higher\-order reasoning, and a noisy\-OR confounder module\. Trained entirely on synthetic structural causal models with anti\-shortcut augmentation,FoundCauseachieves state\-of\-the\-art averageF1F\_\{1\}, AUROC, andSHD/d\\mathrm\{SHD\}/dacross 15 real\-world benchmarks, while being the only method that consistently runs on all datasets and explicitly models latent confounding\. Ablation studies highlight the importance of pairwise statistical features and triangular refinement, while extrapolation experiments demonstrate that edge ranking degrades gracefully beyond the training regime\. Together, these results suggest that combining amortized inference with structured inductive biases is a promising direction for scalable causal discovery\. We plan to open\-sourceFoundCauseto support reproducibility and accelerate research in the broader scientific community\. We hope this work serves both as a practical tool for applied causal analysis and as a step toward foundation models that natively incorporate causal reasoning\.

## References

- \[1\]\(1996\)Hailfinder: a bayesian system for forecasting severe weather\.International Journal of Forecasting12\(1\),pp\. 57–71\.Cited by:[Table 17](https://arxiv.org/html/2606.17516#A11.T17.10.10.6),[Table 18](https://arxiv.org/html/2606.17516#A11.T18.10.10.6),[Appendix C](https://arxiv.org/html/2606.17516#A3.SS0.SSS0.Px5.p1.7),[Table 3](https://arxiv.org/html/2606.17516#S5.T3.10.10.6)\.
- \[2\]M\. Ashman, C\. Ma, A\. Hilmkil, J\. Jennings, and C\. Zhang\(2023\)Causal reasoning in the presence of latent confounders via neural ADMG learning\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=dcN0CaXQhT)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1)\.
- \[3\]V\. Balazadeh, H\. Kamkari, V\. Thomas, J\. Ma, B\. Li, J\. C\. Cresswell, and R\. Krishnan\(2026\)CausalPFN: amortized causal effect estimation via in\-context learning\.InThe Thirty\-ninth Annual Conference on Neural Information Processing Systems,External Links:[Link](https://openreview.net/forum?id=RblaNJGx8C)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p4.1)\.
- \[4\]K\. Bello, B\. Aragam, and P\. Ravikumar\(2022\)DAGMA: learning dags via m\-matrices and a log\-determinant acyclicity characterization\.Advances in Neural Information Processing Systems35,pp\. 8226–8239\.Cited by:[Table 1](https://arxiv.org/html/2606.17516#S5.T1.55.55.6),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.8.8.4)\.
- \[5\]R\. Bhattacharya, T\. Nagarajan, D\. Malinsky, and I\. Shpitser\(2021\)Differentiable causal discovery under unmeasured confounding\.InProceedings of the 24th International Conference on Artificial Intelligence and Statistics,Proceedings of Machine Learning Research, Vol\.130,pp\. 2314–2322\.External Links:[Link](https://proceedings.mlr.press/v130/bhattacharya21a.html)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1)\.
- \[6\]P\. Blöbaum, P\. Götz, K\. Budhathoki, A\. A\. Mastakouri, and D\. Janzing\(2024\)DoWhy\-gcm: an extension of dowhy for causal inference in graphical causal models\.Journal of Machine Learning Research25\(147\),pp\. 1–7\.Cited by:[§G\.1](https://arxiv.org/html/2606.17516#A7.SS1.p1.1),[§5](https://arxiv.org/html/2606.17516#S5.p1.1)\.
- \[7\]P\. Brouillard, S\. Lachapelle, A\. Lacoste, S\. Lacoste\-Julien, and A\. Drouin\(2020\)Differentiable causal discovery from interventional data\.InAdvances in Neural Information Processing Systems,Vol\.33\.External Links:[Link](https://papers.nips.cc/paper/2020/hash/f8b7aa3a0d349d9562b424160ad18612-Abstract.html)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1)\.
- \[8\]P\. Bühlmann, J\. Peters, and J\. Ernest\(2013\-10\)CAM: causal additive models, high\-dimensional order search and penalized regression\.The Annals of Statistics42,pp\.\.External Links:[Document](https://dx.doi.org/10.1214/14-AOS1260)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1)\.
- \[9\]P\. Chao, P\. Blöbaum, S\. Patel, and S\. P\. Kasiviswanathan\(2023\)Modeling causal mechanisms with diffusion models for interventional and counterfactual queries\.Trans\. Mach\. Learn\. Res\.2024\.External Links:[Link](https://api.semanticscholar.org/CorpusID:256503930)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1)\.
- \[10\]D\. M\. Chickering\(2002\)Optimal structure identification with greedy search\.Journal of machine learning research3\(Nov\),pp\. 507–554\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.17516#S1.p2.1),[§2](https://arxiv.org/html/2606.17516#S2.p2.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.30.30.6)\.
- \[11\]D\. Colombo and M\. H\. Maathuis\(2014\)Order\-independent constraint\-based causal structure learning\.Journal of Machine Learning Research15\(116\),pp\. 3921–3962\.External Links:[Link](https://jmlr.org/papers/v15/colombo14a.html)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1),[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1)\.
- \[12\]A\. Defazio, X\. Yang, H\. Mehta, K\. Mishchenko, A\. Khaled, and A\. Cutkosky\(2024\)The road less scheduled\.Advances in Neural Information Processing Systems37,pp\. 9974–10007\.Cited by:[§4](https://arxiv.org/html/2606.17516#S4.p1.11)\.
- \[13\]J\. L\. Gamella, J\. Peters, and P\. Bühlmann\(2025\)Causal chambers as a real\-world physical testbed for ai methodology\.Nature Machine Intelligence7\(1\),pp\. 107–118\.Cited by:[Appendix C](https://arxiv.org/html/2606.17516#A3.SS0.SSS0.Px3.p1.2)\.
- \[14\]T\. Geffner, J\. Antoran, A\. Foster, W\. Gong, C\. Ma, E\. Kiciman, A\. Sharma, A\. Lamb, M\. Kukla, N\. Pawlowski, A\. Hilmkil, J\. Jennings, M\. Scetbon, M\. Allamanis, and C\. Zhang\(2024\)Deep end\-to\-end causal inference\.Transactions on Machine Learning Research\.Note:External Links:ISSN 2835\-8856,[Link](https://openreview.net/forum?id=e6sqttxEGX)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1)\.
- \[15\]C\. Glymour, K\. Zhang, and P\. Spirtes\(2019\)Review of causal discovery methods based on graphical models\.Frontiers in Genetics10,pp\. 524\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1)\.
- \[16\]M\. Hardt, W\. R\. Orchard, P\. Blöbaum, E\. Kirschbaum, and S\. Kasiviswanathan\(2024\)The petshop dataset—finding causes of performance issues across microservices\.InCausal Learning and Reasoning,pp\. 957–978\.Cited by:[Appendix C](https://arxiv.org/html/2606.17516#A3.SS0.SSS0.Px4.p1.4)\.
- \[17\]D\. Hendrycks and K\. Gimpel\(2016\)Gaussian error linear units \(gelus\)\.arXiv preprint arXiv:1606\.08415\.Cited by:[§3](https://arxiv.org/html/2606.17516#S3.p5.4)\.
- \[18\]P\. Hoyer, D\. Janzing, J\. Mooij, J\. Peters, and B\. Schölkopf\(2009\)Nonlinear causal discovery with additive noise models\.InProceedings of the conference Neural Information Processing Systems \(NIPS\) 2008,D\. Koller, D\. Schuurmans, Y\. Bengio, and L\. Bottou \(Eds\.\),Vancouver, Canada\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1),[§1](https://arxiv.org/html/2606.17516#S1.p2.1)\.
- \[19\]J\. Jumper, R\. Evans, A\. Pritzel,et al\.\(2021\)Highly accurate protein structure prediction with alphafold\.Nature596\(7873\),pp\. 583–589\.Cited by:[§3](https://arxiv.org/html/2606.17516#S3.p10.7)\.
- \[20\]N\. R\. Ke, S\. Chiappa, J\. Wang, A\. Goyal, J\. Bornschein, M\. Rey, T\. Weber, M\. Botvinick, M\. C\. Mozer, and D\. J\. Rezende\(2023\)Learning to induce causal structure\.InInternational Conference on Learning Representations \(ICLR\),External Links:2204\.04875Cited by:[Table 4](https://arxiv.org/html/2606.17516#A1.T4.2.2.2.2),[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p4.1),[§1](https://arxiv.org/html/2606.17516#S1.p3.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.75.75.6),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.23.23.4)\.
- \[21\]S\. Lachapelle, P\. Brouillard, T\. Deleu, and S\. Lacoste\-Julien\(2020\)Gradient\-based neural DAG learning\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=rklbKA4YDS)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1)\.
- \[22\]W\. Lam, B\. Andrews, and J\. Ramsey\(2022\)Greedy relaxations of the sparsest permutation algorithm\.InUncertainty in Artificial Intelligence,pp\. 1052–1062\.Cited by:[Table 1](https://arxiv.org/html/2606.17516#S5.T1.35.35.6)\.
- \[23\]J\. Lee, Y\. Lee, J\. Kim, A\. Kosiorek, S\. Choi, and Y\. W\. Teh\(2019\)Set transformer: a framework for attention\-based permutation\-invariant input\.InInternational Conference on Machine Learning \(ICML\),Cited by:[§3](https://arxiv.org/html/2606.17516#S3.p7.8)\.
- \[24\]P\. Lippe, T\. Cohen, and E\. Gavves\(2022\)Efficient neural causal discovery without acyclicity constraints\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=eYciPrLuUhG)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1)\.
- \[25\]L\. Lorch, J\. Rothfuss, B\. Schölkopf, and A\. Krause\(2021\)DiBS: differentiable bayesian structure learning\.InAdvances in Neural Information Processing Systems,Vol\.34\.External Links:[Link](https://proceedings.neurips.cc/paper/2021/hash/ca6ab34959489659f8c3776aaf1f8efd-Abstract.html)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1)\.
- \[26\]L\. Lorch, S\. Sussex, J\. Rothfuss, A\. Krause, and B\. Schölkopf\(2022\)Amortized inference for causal structure learning\.InAdvances in Neural Information Processing Systems \(NeurIPS\),External Links:2205\.12934Cited by:[Table 4](https://arxiv.org/html/2606.17516#A1.T4.1.1.1.2),[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p4.1),[§1](https://arxiv.org/html/2606.17516#S1.p3.1),[§2](https://arxiv.org/html/2606.17516#S2.p2.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.70.70.6),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.20.20.4)\.
- \[27\]P\. Ma, R\. Ding, Q\. Fu, J\. Zhang, S\. Wang, S\. Han, and D\. Zhang\(2024\)Scalable differentiable causal discovery in the presence of latent confounders with skeleton posterior\.InProceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining,pp\. 2141–2152\.External Links:[Document](https://dx.doi.org/10.1145/3637528.3672031),[Link](https://dl.acm.org/doi/10.1145/3637528.3672031)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1)\.
- \[28\]D\. Mahajan, J\. Gladrow, A\. Hilmkil, C\. Zhang, and M\. Scetbon\(2025\)Amortized inference of causal models via conditional fixed\-point iterations\.Transactions on Machine Learning Research\.Note:J2C CertificationExternal Links:ISSN 2835\-8856,[Link](https://openreview.net/forum?id=D9pq25PGc5)Cited by:[Table 4](https://arxiv.org/html/2606.17516#A1.T4.3.3.3.2),[§1](https://arxiv.org/html/2606.17516#S1.p3.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.85.85.6),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.29.29.4)\.
- \[29\]D\. Marbach, R\. J\. Prill, T\. Schaffter, C\. Mattiussi, D\. Floreano, and G\. Stolovitzky\(2010\)Revealing strengths and weaknesses of methods for gene network inference\.Proceedings of the National Academy of Sciences107\(14\),pp\. 6286–6291\.Cited by:[Table 17](https://arxiv.org/html/2606.17516#A11.T17.25.25.6),[Table 18](https://arxiv.org/html/2606.17516#A11.T18.25.25.6),[Table 3](https://arxiv.org/html/2606.17516#S5.T3.25.25.6)\.
- \[30\]D\. Marbach, T\. Schaffter, C\. Mattiussi, and D\. Floreano\(2009\)Generating realistic in silico gene networks for performance assessment of reverse engineering methods\.Journal of Computational Biology16\(2\),pp\. 229–239\.Cited by:[Table 17](https://arxiv.org/html/2606.17516#A11.T17.25.25.6),[Table 18](https://arxiv.org/html/2606.17516#A11.T18.25.25.6),[Table 3](https://arxiv.org/html/2606.17516#S5.T3.25.25.6)\.
- \[31\]A\. Marx and J\. Vreeken\(2019\)Identifiability of cause and effect using regularized regression\.InProceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining,pp\. 852–861\.Cited by:[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.35.35.4)\.
- \[32\]F\. Montagna, M\. Cairney\-Leeming, D\. Sridhar, and F\. Locatello\(2025\)Demystifying amortized causal discovery with transformers\.Transactions on Machine Learning Research\.Note:External Links:ISSN 2835\-8856,[Link](https://openreview.net/forum?id=9Lgy7IGSfp)Cited by:[Table 4](https://arxiv.org/html/2606.17516#A1.T4.7.7.9.1),[§1](https://arxiv.org/html/2606.17516#S1.p3.1)\.
- \[33\]J\. M\. Mooij, J\. Peters, D\. Janzing, J\. Zscheischler, and B\. Schölkopf\(2016\)Distinguishing cause from effect using observational data: methods and benchmarks\.Journal of Machine Learning Research17\(32\),pp\. 1–102\.Cited by:[Appendix C](https://arxiv.org/html/2606.17516#A3.SS0.SSS0.Px6.p1.1),[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p1.3)\.
- \[34\]A\. Nasr\-Esfahany, M\. Alizadeh, and D\. Shah\(2023\)Counterfactual identifiability of bijective causal models\.InForty\-second International Conference on Machine Learning,External Links:[Link](https://openreview.net/)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1)\.
- \[35\]A\. Nazaret and D\. Blei\(2024\)Extremely greedy equivalence search\.InThe 40th Conference on Uncertainty in Artificial Intelligence,External Links:[Link](https://openreview.net/forum?id=2gIMX9UxRN)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1)\.
- \[36\]I\. Ng, A\. Ghassami, and K\. Zhang\(2020\)On the role of sparsity and dag constraints for learning linear dags\.Advances in Neural Information Processing Systems \(NeurIPS\)33,pp\. 17943–17954\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.17516#S1.p2.1)\.
- \[37\]H\. Nilforoshan, M\. Moor, Y\. Roohani, Y\. Chen, A\. Šurina, M\. Yasunaga, S\. Oblak, and J\. Leskovec\(2023\)Zero\-shot causal learning\.Advances in Neural Information Processing Systems36,pp\. 6862–6901\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p4.1)\.
- \[38\]J\. M\. Ogarrio, P\. Spirtes, and J\. Ramsey\(2016\)A hybrid causal search algorithm for latent variable models\.InProceedings of the Eighth International Conference on Probabilistic Graphical Models,Proceedings of Machine Learning Research, Vol\.52,pp\. 368–379\.External Links:[Link](https://proceedings.mlr.press/v52/ogarrio16.html)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1),[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1)\.
- \[39\]A\. Onisko\(2003\)Probabilistic causal models in medicine: application to diagnosis of liver disorders\.InPh\. D\. dissertation, Inst\. Biocybern\. Biomed\. Eng\., Polish Academy Sci\., Warsaw, Poland,Cited by:[Table 17](https://arxiv.org/html/2606.17516#A11.T17.15.15.6),[Table 18](https://arxiv.org/html/2606.17516#A11.T18.15.15.6),[Appendix C](https://arxiv.org/html/2606.17516#A3.SS0.SSS0.Px5.p1.7),[Table 3](https://arxiv.org/html/2606.17516#S5.T3.15.15.6)\.
- \[40\]J\. Peters and P\. Bühlmann\(2014\-03\)Identifiability of gaussian structural equation models with equal error variances\.Biometrika101\(1\),pp\. 219–228\.External Links:ISSN 0006\-3444,[Document](https://dx.doi.org/10.1093/biomet/ast043),[Link](https://doi.org/10.1093/biomet/ast043),https://academic\.oup\.com/biomet/article\-pdf/101/1/219/17460568/ast043\.pdfCited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1)\.
- \[41\]A\. Reisach, M\. Tami, C\. Seiler, A\. Chambaz, and S\. Weichwald\(2023\)A scale\-invariant sorting criterion to find a causal order in additive noise models\.Advances in Neural Information Processing Systems36,pp\. 785–807\.Cited by:[§4](https://arxiv.org/html/2606.17516#S4.p3.2)\.
- \[42\]A\. Reuter, A\. Dhir, C\. Diaconu, J\. Robertson, O\. Ossen, F\. Hutter, A\. Weller, M\. van der Wilk, and B\. Schölkopf\(2026\)Use what you know: causal foundation models with partial graphs\.arXiv preprint arXiv:2602\.14972\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p4.1)\.
- \[43\]J\. Robertson, A\. Reuter, S\. Guo, N\. Hollmann, F\. Hutter, and B\. Schölkopf\(2025\)Do\-pfn: in\-context learning for causal effect estimation\.arXiv preprint arXiv:2506\.06039\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p4.1)\.
- \[44\]P\. Rolland, V\. Cevher, M\. Kleindessner, C\. Russell, D\. Janzing, B\. Schölkopf, and F\. Locatello\(2022\)Score matching enables causal discovery of nonlinear additive noise models\.InInternational Conference on Machine Learning,pp\. 18741–18753\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p2.1),[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.60.60.6)\.
- \[45\]K\. Sachs, O\. Perez, D\. Pe’er, D\. A\. Lauffenburger, and G\. P\. Nolan\(2005\)Causal protein\-signaling networks derived from multiparameter single\-cell data\.Science308\(5721\),pp\. 523–529\.Cited by:[Appendix C](https://arxiv.org/html/2606.17516#A3.SS0.SSS0.Px2.p1.2),[Table 5](https://arxiv.org/html/2606.17516#A3.T5.23.25.1.1)\.
- \[46\]M\. Scutari\(2022\)Bayesian network repository: large discrete bayesian networks\.Note:[https://www\.bnlearn\.com/bnrepository/discrete\-large\.html](https://www.bnlearn.com/bnrepository/discrete-large.html)Accessed: 2026\-05\-03Cited by:[Table 17](https://arxiv.org/html/2606.17516#A11.T17.20.20.6),[Table 18](https://arxiv.org/html/2606.17516#A11.T18.20.20.6),[Appendix C](https://arxiv.org/html/2606.17516#A3.SS0.SSS0.Px1.p1.1),[Appendix C](https://arxiv.org/html/2606.17516#A3.SS0.SSS0.Px5.p1.7),[Table 5](https://arxiv.org/html/2606.17516#A3.T5.23.24.1.1),[Table 5](https://arxiv.org/html/2606.17516#A3.T5.23.27.1.1),[Table 3](https://arxiv.org/html/2606.17516#S5.T3.20.20.6)\.
- \[47\]A\. Sharma and E\. Kiciman\(2020\)Dowhy: an end\-to\-end library for causal inference\.arXiv preprint arXiv:2011\.04216\.Cited by:[§G\.1](https://arxiv.org/html/2606.17516#A7.SS1.p1.1),[§5](https://arxiv.org/html/2606.17516#S5.p1.1)\.
- \[48\]N\. Shazeer\(2020\)GLU variants improve transformer\.arXiv preprint arXiv:2002\.05202\.Cited by:[Appendix E](https://arxiv.org/html/2606.17516#A5.p2.2),[§3](https://arxiv.org/html/2606.17516#S3.p6.1)\.
- \[49\]S\. Shimizu, P\. O\. Hoyer, A\. Hyvärinen, A\. Kerminen, and M\. Jordan\(2006\)A linear non\-gaussian acyclic model for causal discovery\.\.Journal of Machine Learning Research7\(10\)\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.45.45.6),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.17.17.4)\.
- \[50\]S\. Shimizu, T\. Inazumi, Y\. Sogawa, A\. Hyvarinen, Y\. Kawahara, T\. Washio, P\. O\. Hoyer, K\. Bollen, and P\. Hoyer\(2011\)DirectLiNGAM: a direct method for learning a linear non\-gaussian structural equation model\.Journal of Machine Learning Research\-JMLR12\(Apr\),pp\. 1225–1248\.Cited by:[§1](https://arxiv.org/html/2606.17516#S1.p2.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.40.40.6),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.14.14.4)\.
- \[51\]P\. Spirtes, C\. Glymour, and R\. Scheines\(2000\)Causation, prediction, and search\.2nd edition,MIT Press\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1),[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1),[§1](https://arxiv.org/html/2606.17516#S1.p2.1),[§2](https://arxiv.org/html/2606.17516#S2.p2.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.14.14.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.20.20.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.8.8.1),[§5](https://arxiv.org/html/2606.17516#S5.p1.1)\.
- \[52\]G\. Stein, M\. Shadaydeh, and J\. Denzler\(2024\)Embracing the black box: heading towards foundation models for causal discovery from time series data\.arXiv preprint arXiv:2402\.09305\.Cited by:[Table 4](https://arxiv.org/html/2606.17516#A1.T4.7.7.7.3)\.
- \[53\]C\. Uhler, G\. Raskutti, P\. Bühlmann, and B\. Yu\(2013\)Geometry of the faithfulness assumption in causal inference\.The Annals of Statistics41\(2\),pp\. 437–463\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1)\.
- \[54\]A\. Vaswani, N\. Shazeer, N\. Parmar, J\. Uszkoreit, L\. Jones, A\. N\. Gomez, Ł\. Kaiser, and I\. Polosukhin\(2017\)Attention is all you need\.Advances in neural information processing systems30\.Cited by:[§3](https://arxiv.org/html/2606.17516#S3.p6.1)\.
- \[55\]M\. Wu, Y\. Bao, R\. Barzilay, and T\. Jaakkola\(2025\)Sample, estimate, aggregate: a recipe for causal discovery foundation models\.Transactions on Machine Learning Research\.Note:External Links:ISSN 2835\-8856,[Link](https://openreview.net/forum?id=h434zx5SX0)Cited by:[Table 4](https://arxiv.org/html/2606.17516#A1.T4.5.5.5.3),[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p4.1),[§1](https://arxiv.org/html/2606.17516#S1.p3.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.80.80.6),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.26.26.4)\.
- \[56\]S\. Xu, O\. A\. Mian, A\. Marx, and J\. Vreeken\(2022\)Inferring cause and effect in the presence of heteroscedastic noise\.InInternational Conference on Machine Learning,pp\. 24615–24630\.Cited by:[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.32.32.4)\.
- \[57\]Y\. Yu, J\. Chen, T\. Gao, and M\. Yu\(2019\)DAG\-GNN: DAG structure learning with graph neural networks\.InProceedings of the 36th International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.97,pp\. 7154–7163\.External Links:[Link](https://proceedings.mlr.press/v97/yu19a.html)Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1)\.
- \[58\]B\. Zhang and R\. Sennrich\(2019\)Root mean square layer normalization\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Cited by:[Appendix E](https://arxiv.org/html/2606.17516#A5.p1.1),[§3](https://arxiv.org/html/2606.17516#S3.p6.1)\.
- \[59\]J\. Zhang\(2008\)On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias\.Artificial Intelligence172\(16\-17\),pp\. 1873–1896\.Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1),[§1](https://arxiv.org/html/2606.17516#S1.p2.1)\.
- \[60\]K\. Zhang and A\. Hyvärinen\(2009\)On the identifiability of the post\-nonlinear causal model\.InProceedings of the 25th Conference on Uncertainty in Artificial Intelligence,Montreal, Canada\.Note:Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p3.1),[§1](https://arxiv.org/html/2606.17516#S1.p2.1)\.
- \[61\]X\. Zheng, B\. Aragam, P\. Ravikumar, and E\. P\. Xing\(2018\)DAGs with no tears: continuous optimization for structure learning\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Cited by:[§1\.1](https://arxiv.org/html/2606.17516#S1.SS1.p1.1),[§1](https://arxiv.org/html/2606.17516#S1.p2.1),[§2](https://arxiv.org/html/2606.17516#S2.p2.1),[§5\.1](https://arxiv.org/html/2606.17516#S5.SS1.p2.1),[Table 1](https://arxiv.org/html/2606.17516#S5.T1.50.50.6),[Table 2](https://arxiv.org/html/2606.17516#S5.T2.5.5.4)\.

## Appendix AComparison between Amortized Causal Discovery Methods

MethodLatentconfoundersMissingdataVariable\-countagnosticNonlinearmechanismsV\-structure /structural reas\.Pairwisestat\. featuresExplicitDAGHeterogeneouspopulationsReal\-worldbenchmarksAVICI\[[26](https://arxiv.org/html/2606.17516#bib.bib21)\]xx✓✓xxxx∼\\simCSIvA\[[20](https://arxiv.org/html/2606.17516#bib.bib38)\]xx✓✓xxxx∼\\simMontagna et al\. \(2024\)\[[32](https://arxiv.org/html/2606.17516#bib.bib39)\]xx✓✓xxxxxCond\_FIP\[[28](https://arxiv.org/html/2606.17516#bib.bib40)\]xx✓✓xx✓x∼\\simSEA\[[55](https://arxiv.org/html/2606.17516#bib.bib42)\]xx✓✓∼\\sim✓xx∼\\simCausal Pretraining \(TS\)\[[52](https://arxiv.org/html/2606.17516#bib.bib59)\]xx∼\\sim✓xxxx∼\\simFoundCause\(ours\)✓✓✓✓✓✓✓✓✓

Table 4:Feature comparison with existing foundation\-model and amortized causal discovery methods for observational data\.✓denotes full support,xdenotes lack of support, and∼\\simindicates partial support or reliance on post\-hoc extensions\.FoundCauseis the only method that simultaneously supports latent\-confounder prediction, native handling of missing data, variable\-count agnostic inference beyond the training range, and explicit DAG post\-processing\.
## Appendix BEvaluation Metrics

LetD∗∈\{0,1\}d×dD^\{\*\}\\in\\\{0,1\\\}^\{d\\times d\}denote the ground\-truth directed adjacency matrix and letD^∈\{0,1\}d×d\\widehat\{D\}\\in\\\{0,1\\\}^\{d\\times d\}denote the predicted directed adjacency matrix, with diagonal entries excluded\. We evaluate directed\-edge recovery using precision, recall, andF1F\_\{1\}:

Precdir=∑i≠jD^i​j​Di​j∗∑i≠jD^i​j,Recdir=∑i≠jD^i​j​Di​j∗∑i≠jDi​j∗,\\mathrm\{Prec\}\_\{\\mathrm\{dir\}\}=\\frac\{\\sum\_\{i\\neq j\}\\widehat\{D\}\_\{ij\}D^\{\*\}\_\{ij\}\}\{\\sum\_\{i\\neq j\}\\widehat\{D\}\_\{ij\}\},\\qquad\\mathrm\{Rec\}\_\{\\mathrm\{dir\}\}=\\frac\{\\sum\_\{i\\neq j\}\\widehat\{D\}\_\{ij\}D^\{\*\}\_\{ij\}\}\{\\sum\_\{i\\neq j\}D^\{\*\}\_\{ij\}\},and

F1=2​Precdir​RecdirPrecdir\+Recdir\.F\_\{1\}=\\frac\{2\\,\\mathrm\{Prec\}\_\{\\mathrm\{dir\}\}\\mathrm\{Rec\}\_\{\\mathrm\{dir\}\}\}\{\\mathrm\{Prec\}\_\{\\mathrm\{dir\}\}\+\\mathrm\{Rec\}\_\{\\mathrm\{dir\}\}\}\.Thus, a predicted edgei→ji\\to jis counted as correct only when the true DAG also containsi→ji\\to j; predictingj→ij\\to iis treated as an orientation error\.

For direction\-agnostic adjacency recovery, we report skeletonF1F\_\{1\}\. Define the true and predicted skeletons over unordered pairs by

Si​j∗=𝟙​\{Di​j∗=1​or​Dj​i∗=1\},S^i​j=𝟙​\{D^i​j=1​or​D^j​i=1\},i<j\.S^\{\*\}\_\{ij\}=\\mathbbm\{1\}\\\{D^\{\*\}\_\{ij\}=1\\ \\text\{or\}\\ D^\{\*\}\_\{ji\}=1\\\},\\qquad\\widehat\{S\}\_\{ij\}=\\mathbbm\{1\}\\\{\\widehat\{D\}\_\{ij\}=1\\ \\text\{or\}\\ \\widehat\{D\}\_\{ji\}=1\\\},\\qquad i<j\.Skeleton precision and recall are then

Precskel=∑i<jS^i​j​Si​j∗∑i<jS^i​j,Recskel=∑i<jS^i​j​Si​j∗∑i<jSi​j∗,\\mathrm\{Prec\}\_\{\\mathrm\{skel\}\}=\\frac\{\\sum\_\{i<j\}\\widehat\{S\}\_\{ij\}S^\{\*\}\_\{ij\}\}\{\\sum\_\{i<j\}\\widehat\{S\}\_\{ij\}\},\\qquad\\mathrm\{Rec\}\_\{\\mathrm\{skel\}\}=\\frac\{\\sum\_\{i<j\}\\widehat\{S\}\_\{ij\}S^\{\*\}\_\{ij\}\}\{\\sum\_\{i<j\}S^\{\*\}\_\{ij\}\},with

Skel​\-​F1=2​Precskel​RecskelPrecskel\+Recskel\.\\mathrm\{Skel\}\\text\{\-\}F\_\{1\}=\\frac\{2\\,\\mathrm\{Prec\}\_\{\\mathrm\{skel\}\}\\mathrm\{Rec\}\_\{\\mathrm\{skel\}\}\}\{\\mathrm\{Prec\}\_\{\\mathrm\{skel\}\}\+\\mathrm\{Rec\}\_\{\\mathrm\{skel\}\}\}\.This metric measures whether the correct variable pairs are connected, regardless of edge orientation\.

We also report structural Hamming distance normalized by graph size, denotedSHD/d\\mathrm\{SHD\}/d\. SHD counts the number of graph edits required to transform the predicted graph into the true DAG, including missing edges, extra edges, and incorrectly oriented edges:

SHD/d=SHD​\(D^,D∗\)d\.\\mathrm\{SHD\}/d=\\frac\{\\mathrm\{SHD\}\(\\widehat\{D\},D^\{\*\}\)\}\{d\}\.Normalizing byddmakes structural errors more comparable across datasets with different numbers of variables\. LowerSHD/d\\mathrm\{SHD\}/dindicates better graph recovery\.

Finally, we report AUROC for directed\-edge prediction over the ordered pairs\(i,j\)\(i,j\)withi≠ji\\neq j\. When probability scores are available, AUROC measures the ability of the model to rank true directed edges above non\-edges\. When computed from a binary predicted DAG, AUROC corresponds to the ROC performance at a single operating point induced by the final thresholded graph:

AUROC=Pr⁡\(si​j\>sk​l∣Di​j∗=1,Dk​l∗=0\)\+12​Pr⁡\(si​j=sk​l∣Di​j∗=1,Dk​l∗=0\),\\mathrm\{AUROC\}=\\Pr\\\!\\left\(s\_\{ij\}\>s\_\{kl\}\\mid D^\{\*\}\_\{ij\}=1,\\ D^\{\*\}\_\{kl\}=0\\right\)\+\\frac\{1\}\{2\}\\Pr\\\!\\left\(s\_\{ij\}=s\_\{kl\}\\mid D^\{\*\}\_\{ij\}=1,\\ D^\{\*\}\_\{kl\}=0\\right\),wheresi​js\_\{ij\}is either the predicted edge score or the binary predictionD^i​j\\widehat\{D\}\_\{ij\}\. Higher AUROC indicates better discrimination between true directed edges and absent directed edges\.

## Appendix CEvaluation Datasets

We evaluateFoundCauseon 15 real\-world and semi\-synthetic causal discovery benchmarks spanningD∈\[8,46\]D\\in\[8,46\]variables, which lie within the training rangeD∈\[2,50\]D\\in\[2,50\], as well as on the Tübingen cause–effect pairs benchmark for bivariate direction\. In addition, we include four datasets withD\>50D\>50to assess generalization beyond the training regime\. None of these datasets, their ground\-truth DAGs, or any derived structural information are used during training; the model is trained*exclusively*on synthetic SCMs \(Section[G](https://arxiv.org/html/2606.17516#A7)\)\. Table[5](https://arxiv.org/html/2606.17516#A3.T5)summarizes the benchmark suite, and we briefly describe each dataset family below\.

DatasetDomainDDNN\|E\|\|E\|ρ\\rhoType*Semi\-synthetic Bayesian networks \(bnlearn repository\[[46](https://arxiv.org/html/2606.17516#bib.bib29)\]\)*asiaMedical \(lung cancer\)81,00080\.1430\.143Discrete BNasia\_nonlinearMedical \(lung cancer\)81,00080\.1430\.143Nonlinear variantchildMedical \(congenital heart\)201,000250\.0660\.066Discrete BNchild\_nonlinearMedical \(congenital heart\)201,000250\.0660\.066Nonlinear variantinsuranceInsurance risk assessment271,500520\.0740\.074Discrete BNalarm\_likeMedical \(ICU monitoring\)372,000460\.0350\.035Linear Gaussianecoli\_likeGene regulation \(E\. coli\)461,000580\.0280\.028Linear Gaussian*Real biological data\[[45](https://arxiv.org/html/2606.17516#bib.bib28)\]*sachs\_obsProtein signalling \(obs\.\)11853170\.1540\.154Observationalsachs\_fullProtein signalling \(mixed\)117,466170\.1540\.154Flow cytometrysachs\_nonlinearProtein signalling111,000170\.1540\.154Nonlinear variant*Real physical / operational systems*causal\_chambers\_ltLight tunnel \(ETH Zürich\)2010,000390\.1030\.103Physical measurementspetshop\_high\_trafficMicroservice telemetry39589420\.0280\.028Operational metricspetshop\_low\_trafficMicroservice telemetry41589430\.0260\.026Operational metricspetshop\_temp1Microservice \(temporal\)411,652430\.0260\.026Time\-series metricspetshop\_temp2Microservice \(temporal\)411,652430\.0260\.026Time\-series metrics*Large Bayesian networks \(bnlearn repository\[[46](https://arxiv.org/html/2606.17516#bib.bib29)\], beyond training range\)*hailfinderWeather forecasting565,000660\.0210\.021Discrete BNhepar2Medical \(liver disorders\)705,0001230\.0250\.025Discrete BNwin95ptsPrinter troubleshooting765,0001120\.0200\.020Discrete BN*Bivariate direction benchmark*tuebingenMixed \(102 pairs\)2variable—0\.50\.5Real\-world pairsTable 5:Evaluation benchmarks\.DDdenotes the number of observed variables,NNthe number of samples used at inference,\|E\|\|E\|the number of edges in the ground\-truth DAG, andρ=\|E\|/\(D​\(D−1\)\)\\rho=\|E\|/\(D\(D\-1\)\)the edge density\.*Type*indicates the data source, including discrete Bayesian Networks \(BN\), nonlinear variants, real observational datasets, and physical measurement systems\.##### Semi\-synthetic Bayesian Networks\.

asia,child,insurance, andalarm\_likeare derived from the bnlearn repository\[[46](https://arxiv.org/html/2606.17516#bib.bib29)\]\. Each dataset is generated from a fixed expert\-designed directed acyclic graph \(DAG\) with discrete conditional probability tables \(CPTs\), and observations are obtained via ancestral sampling\. The ground\-truth graph is therefore known exactly, while the observed data follow discrete generative assumptions\.asiais a classic eight\-variable lung cancer network;childmodels congenital heart disease diagnosis;insurancecaptures insurance risk factors; andalarm\_likeis based on an ICU monitoring network\.

We additionally evaluate onasia\_nlandchild\_nl, which retain the same DAG structures but apply nonlinear transformations \(e\.g\.,tanh\\tanh, quadratic, exponential\) with additive Gaussian noise, yielding continuous observations under the same causal graph\.ecoli\_likeis a linear\-Gaussian SCM constructed from a published*E\. coli*gene regulatory network \(46 variables, 58 edges\), providing a biologically motivated benchmark at moderate scale\.

##### Sachs Protein Signaling\.

The Sachs dataset\[[45](https://arxiv.org/html/2606.17516#bib.bib28)\]contains measurements of 11 phosphorylated proteins and phospholipids collected via flow cytometry under multiple experimental conditions\.sachs\_obsuses the observational \(unperturbed\) subset \(N=853N=853\),sachs\_fulluses the full dataset including interventional samples \(N=7,466N=7\{,\}466, treated as observational during evaluation\), andsachs\_nlis a nonlinear re\-simulation based on the same underlying DAG\. All three variants share the 17\-edge consensus graph reported in the original study\.

##### Causal Chambers \(Light Tunnel\)\.

The Causal Chamber light\-tunnel system\[[13](https://arxiv.org/html/2606.17516#bib.bib27)\]is a physical experimental setup at ETH Zürich consisting of a controllable light source, rotating polarizers, a camera, and wavelength\-specific intensity sensors\. Variables correspond to sensor measurements and actuator settings, and the ground\-truth causal graph is determined by the known physical interactions among these components\. Thecausal\_chambers\_ltdataset containsD=20D=20variables andN=10,000N=10\{,\}000observations collected under the natural operating regime\. This benchmark is notable in that the causal graph is directly specified by the underlying physics, rather than inferred from observational data or expert consensus\.

##### PetShop Operational Telemetry\.

The PetShop dataset\[[16](https://arxiv.org/html/2606.17516#bib.bib26)\]consists of runtime telemetry from a synthetic microservice\-based e\-commerce system\. Each variable represents an operational metric \(e\.g\., latency, error rate, or traffic volume\) for a particular service, and the ground\-truth DAG is derived from the known service dependency structure\. We evaluate four variants:petshop\_high\_trafficandpetshop\_low\_traffic, which reflect steady\-state regimes \(D≈39D\\approx 39–4141,N=589N=589\), andpetshop\_temp1andpetshop\_temp2, which include temporal variability \(N=1,652N=1\{,\}652\)\. These datasets are challenging due to strong correlations induced by service dependencies and non\-additive interactions arising from system dynamics such as queuing, retries, and cascading effects\.

##### Large Bayesian Networks \(Dimension Extrapolation\)\.

To evaluate generalization beyond the training rangeD∈\[2,50\]D\\in\[2,50\], we consider three expert\-designed Bayesian networks from the bnlearn repository\[[46](https://arxiv.org/html/2606.17516#bib.bib29)\]withD∈\{56,70,76\}D\\in\\\{56,70,76\\\}\.hailfinder\[[1](https://arxiv.org/html/2606.17516#bib.bib30)\]is a 56\-node, 66\-edge network for severe weather forecasting, constructed from meteorological expert knowledge, with relatively sparse connectivity \(average degree2\.362\.36\)\.hepar2\[[39](https://arxiv.org/html/2606.17516#bib.bib25)\]is a 70\-node, 123\-edge medical diagnostic network for liver disorders, with higher edge density \(average degree3\.513\.51\)\.win95pts\[[46](https://arxiv.org/html/2606.17516#bib.bib29)\]is a 76\-node, 112\-edge troubleshooting network for Windows 95 printing issues, exhibiting dense local structure with large Markov blankets and in\-degree up to77\. For each network, we generateN=5,000N=5\{,\}000samples via ancestral sampling from the corresponding conditional probability tables, using the published DAG as ground truth\. These datasets represent strict zero\-shot evaluation forFoundCause, as no graphs withD\>50D\>50are seen during training\.

##### Tübingen Cause–Effect Pairs\.

We further evaluate bivariate direction identification using the Tübingen cause–effect pairs benchmark\[[33](https://arxiv.org/html/2606.17516#bib.bib18)\], which consists of real\-world variable pairs with known causal direction across diverse domains\. We use the 102 pairs with positive evaluation weights and report both unweighted and weighted accuracy, where weights reflect pair difficulty as defined by the benchmark authors\. This benchmark isolates the directionality task independently of edge existence, providing a complementary evaluation of causal orientation performance\.

## Appendix DInference Protocol

All datasets are evaluated using a fixed inference pipeline\. Given an input dataset, we performK=10K=10stochastic inference runs, each combining bootstrap resampling and random permutation of variable indices, and average the resulting logits after mapping predictions back to the original ordering\. The averaged logits are calibrated using temperature scaling withT=0\.65T=0\.65, estimated on a synthetic validation set, and thresholded using a two\-component Gaussian mixture model \(GMM\)\. Finally, self\-consistency pruning retains only edges that appear in at least50%50\\%of the runs\. No per\-dataset hyperparameter tuning or fine\-tuning is performed: a single fixed checkpoint is used across all datasets\.

## Appendix EAdditional Architecture Details

ParameterValueHidden dimensiondhd\_\{h\}768768Number of attention headsHH88Number of encoder layersLL88\(1616alternating blocks\)FFN inner dimension23×4​dh=2048\\frac\{2\}\{3\}\\times 4d\_\{h\}=2048Triangular pair dimensiondpaird\_\{\\mathrm\{pair\}\}192192Triangular refinement roundsRR33Number of confounder tokensKcK\_\{c\}88Number of global context tokensKgK\_\{g\}1212PMA query tokensKqK\_\{q\}88PMA attention heads44Number of pairwise statistics4545Logit bias initialization−2\.0\-2\.0Table 6:Architectural hyperparameters\.Normalization and Feedforward Blocks\.Each encoder block operates on input representations𝑯∈ℝN×D×dh\\bm\{H\}\\in\\mathbb\{R\}^\{N\\times D\\times d\_\{h\}\}and outputs tensors of the same shape\. All blocks use RMSNorm\[[58](https://arxiv.org/html/2606.17516#bib.bib14)\]applied along the feature dimension:

RMSNorm⁡\(𝒙\)=𝒙dh−1​∑r=1dhxr2\+ϵ⊙𝜸,ϵ=10−6,\\operatorname\{RMSNorm\}\(\\bm\{x\}\)=\\frac\{\\bm\{x\}\}\{\\sqrt\{d\_\{h\}^\{\-1\}\\sum\_\{r=1\}^\{d\_\{h\}\}x\_\{r\}^\{2\}\+\\epsilon\}\}\\odot\\bm\{\\gamma\},\\qquad\\epsilon=10^\{\-6\},where𝒙∈ℝdh\\bm\{x\}\\in\\mathbb\{R\}^\{d\_\{h\}\}is a feature vector and𝜸∈ℝdh\\bm\{\\gamma\}\\in\\mathbb\{R\}^\{d\_\{h\}\}is a learned scale parameter\.

The feedforward sublayer uses a SwiGLU activation\[[48](https://arxiv.org/html/2606.17516#bib.bib15)\], mapping𝒙∈ℝdh\\bm\{x\}\\in\\mathbb\{R\}^\{d\_\{h\}\}toℝdh\\mathbb\{R\}^\{d\_\{h\}\}via

SwiGLU​\(𝒙\)=𝑾2​\[SiLU⁡\(𝑾g​𝒙\)⊙𝑾1​𝒙\],\\mathrm\{SwiGLU\}\(\\bm\{x\}\)=\\bm\{W\}\_\{2\}\\left\[\\operatorname\{SiLU\}\(\\bm\{W\}\_\{g\}\\bm\{x\}\)\\odot\\bm\{W\}\_\{1\}\\bm\{x\}\\right\],where𝑾g,𝑾1∈ℝdf×dh\\bm\{W\}\_\{g\},\\bm\{W\}\_\{1\}\\in\\mathbb\{R\}^\{d\_\{f\}\\times d\_\{h\}\}and𝑾2∈ℝdh×df\\bm\{W\}\_\{2\}\\in\\mathbb\{R\}^\{d\_\{h\}\\times d\_\{f\}\}, with hidden dimensiondfd\_\{f\}\.

Both attention and feedforward outputs are added through residual connections, and output projections are initialized with standard deviation0\.02/2​L0\.02/\\sqrt\{2L\}to stabilize deep stacking of2​L2Lblocks\. No dropout is used in the encoder; instead, regularization is provided by continuous data regeneration and feature\-level dropout applied to selected pairwise statistics\.

Multi\-layer Feature Fusion\.Let𝑯\(ℓ\)∈ℝN×D×dh\\bm\{H\}^\{\(\\ell\)\}\\in\\mathbb\{R\}^\{N\\times D\\times d\_\{h\}\}denote the encoder output at layerℓ\\ell\. We retain a subset of intermediate representations\{𝑯\(ℓk\)\}k=14\\\{\\bm\{H\}^\{\(\\ell\_\{k\}\)\}\\\}\_\{k=1\}^\{4\}, including the final layer, and combine them prior to pooling via a learned convex combination:

𝑯fused=∑k=14exp⁡\(αk\)∑k′exp⁡\(αk′\)​𝑯\(ℓk\),\\bm\{H\}\_\{\\mathrm\{fused\}\}=\\sum\_\{k=1\}^\{4\}\\frac\{\\exp\(\\alpha\_\{k\}\)\}\{\\sum\_\{k^\{\\prime\}\}\\exp\(\\alpha\_\{k^\{\\prime\}\}\)\}\\bm\{H\}^\{\(\\ell\_\{k\}\)\},whereαk∈ℝ\\alpha\_\{k\}\\in\\mathbb\{R\}are learnable scalar weights\.

The fused representation𝑯fused∈ℝN×D×dh\\bm\{H\}\_\{\\mathrm\{fused\}\}\\in\\mathbb\{R\}^\{N\\times D\\times d\_\{h\}\}is then passed to the pooling module\. This fusion allows downstream prediction to leverage causal signals that may emerge at different depths of the encoder\.

Pairwise Statistics\.Given input data𝑿∈ℝN×D\\bm\{X\}\\in\\mathbb\{R\}^\{N\\times D\}, the model computes a tensor of pairwise statistics𝑺∈ℝD×D×ds\\bm\{S\}\\in\\mathbb\{R\}^\{D\\times D\\times d\_\{s\}\}withds=45d\_\{s\}=45, where𝒔i​j∈ℝds\\bm\{s\}\_\{ij\}\\in\\mathbb\{R\}^\{d\_\{s\}\}summarizes statistical relationships between variablesiiandjj\. These statistics are computed directly from the raw data undertorch\.no\_grad, i\.e\., without gradient updates\.

The feature vector𝒔i​j\\bm\{s\}\_\{ij\}is partitioned into four groups: symmetric features \(invariant under swappingiiandjj\), antisymmetric features \(changing sign or direction under swapping\), V\-structure features capturing collider\-like patterns, and reliability metadata describing statistical quality\. The symmetric features include Pearson and Spearman correlation, partial correlation, nonlinear dependence measures, and higher\-order moment interactions\. The antisymmetric features include residual\-variance asymmetry, kernel\-based dependence asymmetry \(HSIC\-RFF\), cross\-moment asymmetries, regression\-based coefficients, conditional\-variance features, and directional statistics such as Chatterjee’sξ\\xiand IGCI\-style scores\. The V\-structure features capture patterns indicative of collider structures, while the reliability metadata include the covariance condition number, the ratioD/ND/N, the unique\-value ratio, and the normalized observation count\.

These metadata features are always retained during training, allowing the model to learn when other statistics are unreliable, while the remaining features may be stochastically dropped for regularization\.

Stat\-conditioned Attention Bias\.The pairwise statistics𝑺\\bm\{S\}are used to construct attention biases for variable\-wise self\-attention\. For each layerℓ\\elland attention headh∈\{1,…,H\}h\\in\\\{1,\\ldots,H\\\}, a nonlinear projection maps𝒔i​j\\bm\{s\}\_\{ij\}to a scalar bias:

bi​j\(h,ℓ\)=\[𝑾2\(ℓ\)​GELU⁡\(𝑾1\(ℓ\)​𝒔i​j\)\]h,b^\{\(h,\\ell\)\}\_\{ij\}=\\left\[\\bm\{W\}^\{\(\\ell\)\}\_\{2\}\\operatorname\{GELU\}\\\!\\left\(\\bm\{W\}^\{\(\\ell\)\}\_\{1\}\\bm\{s\}\_\{ij\}\\right\)\\right\]\_\{h\},where𝑾1\(ℓ\)∈ℝdb×ds\\bm\{W\}^\{\(\\ell\)\}\_\{1\}\\in\\mathbb\{R\}^\{d\_\{b\}\\times d\_\{s\}\},𝑾2\(ℓ\)∈ℝH×db\\bm\{W\}^\{\(\\ell\)\}\_\{2\}\\in\\mathbb\{R\}^\{H\\times d\_\{b\}\}, anddb=48d\_\{b\}=48is the hidden dimension of the bias network\. The resulting biases form a tensor𝑩\(ℓ\)∈ℝH×D×D\\bm\{B\}^\{\(\\ell\)\}\\in\\mathbb\{R\}^\{H\\times D\\times D\}, which is added to the attention logits\.

To ensure numerical stability, biases are clamped to the range\[−10,10\]\[\-10,10\]\. During training, antisymmetric and V\-structure features are subjected to inverted dropout, while symmetric features and reliability metadata are always retained\. This design allows the model to leverage classical dependence cues while learning robustness to noisy or unreliable statistical signals\.

Encoder\-conditioned Feature Gate\.For each ordered pair\(i,j\)\(i,j\), the direction head constructs a feature gate𝒈i​j∈\[0,1\]da\\bm\{g\}\_\{ij\}\\in\[0,1\]^\{d\_\{a\}\}, wheredad\_\{a\}is the number of antisymmetric statistics \(hereda=20d\_\{a\}=20\)\. The gate is computed from the variable embeddings𝒛i,𝒛j∈ℝdh\\bm\{z\}\_\{i\},\\bm\{z\}\_\{j\}\\in\\mathbb\{R\}^\{d\_\{h\}\}and pairwise statistics𝒔i​j∈ℝds\\bm\{s\}\_\{ij\}\\in\\mathbb\{R\}^\{d\_\{s\}\}as

𝒈i​j=σ​\(𝑾2​GELU⁡\(𝑾1​\[𝒛i−𝒛j;𝒔i​j\]\)\),\\bm\{g\}\_\{ij\}=\\sigma\\\!\\left\(\\bm\{W\}\_\{2\}\\operatorname\{GELU\}\\\!\\left\(\\bm\{W\}\_\{1\}\[\\bm\{z\}\_\{i\}\-\\bm\{z\}\_\{j\};\\bm\{s\}\_\{ij\}\]\\right\)\\right\),where𝑾1∈ℝdg×\(dh\+ds\)\\bm\{W\}\_\{1\}\\in\\mathbb\{R\}^\{d\_\{g\}\\times\(d\_\{h\}\+d\_\{s\}\)\},𝑾2∈ℝda×dg\\bm\{W\}\_\{2\}\\in\\mathbb\{R\}^\{d\_\{a\}\\times d\_\{g\}\}, anddgd\_\{g\}is a hidden dimension\. The gate modulates antisymmetric features used for direction prediction, allowing the model to adaptively select informative causal signals for each variable pair\.

To prevent saturation, we apply a soft bounds regularization that penalizes gate values close to0or11, encouraging flexible feature usage across pairs\.

Triangular Refinement\.Let𝑬∈ℝD×D×dpair\\bm\{E\}\\in\\mathbb\{R\}^\{D\\times D\\times d\_\{\\mathrm\{pair\}\}\}denote the pairwise representation tensor, initialized from the concatenated existence and direction features\. We refine𝑬\\bm\{E\}forRRrounds using triangle\-based updates that aggregate information over intermediate variablesk∈\{1,…,D\}k\\in\\\{1,\\ldots,D\\\}\. Each round consists of outgoing, incoming, and collider updates, followed by row\-wise self\-attention over pairs\.

The outgoing update is given by

Δ​𝒆i​jout=1D​∑k=1D\(12​tanh⁡\(𝒘g⊤​𝒆i​k\)\+12\)​𝑾v​𝒆k​j,\\Delta\\bm\{e\}^\{\\mathrm\{out\}\}\_\{ij\}=\\frac\{1\}\{\\sqrt\{D\}\}\\sum\_\{k=1\}^\{D\}\\left\(\\frac\{1\}\{2\}\\tanh\(\\bm\{w\}\_\{g\}^\{\\top\}\\bm\{e\}\_\{ik\}\)\+\\frac\{1\}\{2\}\\right\)\\bm\{W\}\_\{v\}\\bm\{e\}\_\{kj\},where𝒆i​j∈ℝdpair\\bm\{e\}\_\{ij\}\\in\\mathbb\{R\}^\{d\_\{\\mathrm\{pair\}\}\},𝒘g∈ℝdpair\\bm\{w\}\_\{g\}\\in\\mathbb\{R\}^\{d\_\{\\mathrm\{pair\}\}\}, and𝑾v∈ℝdpair×dpair\\bm\{W\}\_\{v\}\\in\\mathbb\{R\}^\{d\_\{\\mathrm\{pair\}\}\\times d\_\{\\mathrm\{pair\}\}\}\. Incoming and collider updates use analogous gated aggregation with index permutations corresponding to fork and collider motifs\.

Each update is followed by normalization and residual connections, and the refined tensor remains inℝD×D×dpair\\mathbb\{R\}^\{D\\times D\\times d\_\{\\mathrm\{pair\}\}\}\. The final refined logits are obtained from𝑬\\bm\{E\}and combined with the base decoder outputs using learned mixing weights for skeleton and direction components\. This refinement enables reasoning over higher\-order causal structures such as chains, forks, and colliders\.

Learned Skeleton Extraction\.To obtain an undirected skeleton from directed logits, the model learns a smooth functionskel​\_​fn:ℝ2→ℝ\\mathrm\{skel\\\_fn\}:\\mathbb\{R\}^\{2\}\\to\\mathbb\{R\}applied to each pair\(i,j\)\(i,j\):

skel​\_​fn​\(ℓi​j,ℓj​i\)=𝒘⊤​\[ℓi​j\+ℓj​i,ℓi​j​ℓj​i\]\+b,\\mathrm\{skel\\\_fn\}\(\\ell\_\{ij\},\\ell\_\{ji\}\)=\\bm\{w\}^\{\\top\}\\bigl\[\\ell\_\{ij\}\+\\ell\_\{ji\},\\,\\ell\_\{ij\}\\ell\_\{ji\}\\bigr\]\+b,where𝒘∈ℝ2\\bm\{w\}\\in\\mathbb\{R\}^\{2\}andb∈ℝb\\in\\mathbb\{R\}are learned parameters\. The resulting scalar is used as the skeleton logit for edge existence betweeniiandjj\. This formulation provides a differentiable alternative to hard max or OR operations, while preserving information from both directions\.

Confounder Implementation Details\.Let𝒁∈ℝD×dh\\bm\{Z\}\\in\\mathbb\{R\}^\{D\\times d\_\{h\}\}denote the variable embeddings and𝑻∈ℝKc×dh\\bm\{T\}\\in\\mathbb\{R\}^\{K\_\{c\}\\times d\_\{h\}\}the learnable confounder tokens\. The tokens attend to𝒁\\bm\{Z\}via two cross\-attention layers, producing updated token representations𝑻′∈ℝKc×dh\\bm\{T\}^\{\\prime\}\\in\\mathbb\{R\}^\{K\_\{c\}\\times d\_\{h\}\}\. Each layer consists of multi\-head attention followed by a feedforward network with GELU activation and LayerNorm, preserving the token dimensionality\.

For each variableiiand confounderkk, a loading network computesSi​k∈\(0,1\)S\_\{ik\}\\in\(0,1\), and pairwise confounding probabilities are obtained via the noisy\-OR aggregation

C^i​j=1−∏k=1Kc\(1−Si​k​Sj​k\)\.\\hat\{C\}\_\{ij\}=1\-\\prod\_\{k=1\}^\{K\_\{c\}\}\(1\-S\_\{ik\}S\_\{jk\}\)\.This computation is performed in float32 for numerical stability, and the resulting probabilities are clamped before conversion to logits\.

The loading network bias is initialized to−2\.0\-2\.0, so that initial loadings satisfySi​k≈σ​\(−2\)≈0\.12S\_\{ik\}\\approx\\sigma\(\-2\)\\approx 0\.12, encouraging sparse confounder assignments at the start of training\. In addition, a pairwise feature projection maps symmetric and gated antisymmetric statistics𝒔i​j\\bm\{s\}\_\{ij\}to a scalar bias, which is added to the noisy\-OR logit prior to symmetrization and masking\. The final output is a symmetric matrix𝑪^∈\[0,1\]D×D\\hat\{\\bm\{C\}\}\\in\[0,1\]^\{D\\times D\}of confounding probabilities\.

Acyclicity\.Letℓ∈ℝD×D\\bm\{\\ell\}\\in\\mathbb\{R\}^\{D\\times D\}denote the directed edge logits andσ​\(ℓ\)\\sigma\(\\bm\{\\ell\}\)the corresponding edge probabilities\. An optional acyclicity penalty is defined using the spectral radiusρ​\(σ​\(ℓ\)\)\\rho\(\\sigma\(\\bm\{\\ell\}\)\), estimated via power iteration\. The penalty takes the form

ℒacyc=ReLU​\(ρ​\(σ​\(ℓ\)\)−1\),\\mathcal\{L\}\_\{\\mathrm\{acyc\}\}=\\mathrm\{ReLU\}\\bigl\(\\rho\(\\sigma\(\\bm\{\\ell\}\)\)\-1\\bigr\),which encourages the learned graph to be acyclic\.

In the default configuration, this penalty is disabled during training\. Instead, acyclicity is enforced at inference time via greedy DAG post\-processing applied toσ​\(ℓ\)\\sigma\(\\bm\{\\ell\}\), ensuring a valid directed acyclic graph\.

𝑯\(ℓ\)∈ℝB×N×D×dh\\bm\{H\}^\{\(\\ell\)\}\\in\\mathbb\{R\}^\{B\\times N\\times D\\times d\_\{h\}\}Variable\-attention block\(even\)attends acrossDD\+Kg=12\+K\_\{g\}\{=\}12global tokens\+\+stat bias \(clamp±10\\pm 10\)RMSNorm \(fp32\)\+\+Q/K/V projectionSDPA \(FlashAttn\)\+\+stat biasResidualRMSNorm\+\+SwiGLU \(clamp±10​k\\pm 10\\mathrm\{k\}\)Residualswapaxes\(N↔\\leftrightarrowD\)Sample\-attention block\(odd\)attends acrossNNno stat bias, no global tokens𝑯\(ℓ\+1\)\\bm\{H\}^\{\(\\ell\+1\)\}Figure 3:Encoder Block Pair \(one of eight\)\.The encoder consists of2​L=162L=16attention blocks that alternate between variable\-axis attention \(even\-indexed blocks\) and sample\-axis attention \(odd\-indexed blocks\)\. Variable\-attention blocks attend across variables and incorporate stat\-conditioned attention bias andKg=12K\_\{g\}=12global context tokens, while sample\-attention blocks attend across samples and use neither\. Each block follows a pre\-normalization design with RMSNorm \(computed in float32 withε=10−6\\varepsilon=10^\{\-6\}\), multi\-head scaled dot\-product attention \(SDPA\), and a SwiGLU feedforward network\. SDPA uses FlashAttention with a padding mask value of−109\-10^\{9\}to ensure numerical stability under bfloat16 precision\. The SwiGLU gated product is clamped to±104\\pm 10^\{4\}before the output projection\. Residual output projections use scaled initialization with standard deviation0\.02/2​L0\.02/\\sqrt\{2L\}\. An additional RMSNorm is applied to the residual stream after every fourth block\.Existence headLinear→\\toGELU→𝒉exist∈ℝdh\\to\\bm\{h\}^\{\\mathrm\{exist\}\}\\in\\mathbb\{R\}^\{d\_\{h\}\}𝒛i\+𝒛j\\bm\{z\}\_\{i\}\{\+\}\\bm\{z\}\_\{j\}𝒛i⊙𝒛j\\bm\{z\}\_\{i\}\{\\odot\}\\bm\{z\}\_\{j\}projsym​\(𝒔sym\)\\mathrm\{proj\}\_\{\\mathrm\{sym\}\}\(\\bm\{s\}^\{\\mathrm\{sym\}\}\)C^i​j\.𝚍𝚎𝚝𝚊𝚌𝚑​\(\)\\hat\{C\}\_\{ij\}\.\\mathtt\{detach\}\(\)SYMMETRIC inputs:Direction headLinear→\\toGELU→𝒉dir∈ℝdh\\to\\bm\{h\}^\{\\mathrm\{dir\}\}\\in\\mathbb\{R\}^\{d\_\{h\}\}𝒛i−𝒛j\\bm\{z\}\_\{i\}\{\-\}\\bm\{z\}\_\{j\}projasym​\(𝒈⊙𝒔asym\)\\mathrm\{proj\}\_\{\\mathrm\{asym\}\}\(\\bm\{g\}\{\\odot\}\\bm\{s\}^\{\\mathrm\{asym\}\}\)ANTISYMMETRIC inputs:Fusion MLP\(fp32\): concat\[𝒉exist;𝒉dir\]∈ℝ2​dh\[\\bm\{h\}^\{\\mathrm\{exist\}\};\\,\\bm\{h\}^\{\\mathrm\{dir\}\}\]\\in\\mathbb\{R\}^\{2d\_\{h\}\}Linear\(2​dh→dh\)→\(2d\_\{h\}\\to d\_\{h\}\)\\toGELU→\\toLinear\(dh→1\)\(d\_\{h\}\\to 1\)ℓi​jbase∈ℝ\\bm\{\\ell\}^\{\\mathrm\{base\}\}\_\{ij\}\\in\\mathbb\{R\}\(clamp±15\\pm 15\)𝒉exist\\bm\{h\}^\{\\mathrm\{exist\}\}𝒉dir\\bm\{h\}^\{\\mathrm\{dir\}\}left halfright halfFigure 4:Factored Edge Predictor\.The model predicts each directed edge\(i,j\)\(i,j\)by separating*edge existence*and*edge direction*, using symmetry\-aware features\. The*existence head*\(top\) operates on symmetric inputs, including𝒛i\+𝒛j\\bm\{z\}\_\{i\}\+\\bm\{z\}\_\{j\},𝒛i⊙𝒛j\\bm\{z\}\_\{i\}\\odot\\bm\{z\}\_\{j\}, a projection of symmetric pairwise statistics, and the confounding scoreC^i​j\\hat\{C\}\_\{ij\}\(detached from gradients\)\. These features capture whether two variables are related, regardless of direction\. The*direction head*\(middle\) operates on antisymmetric inputs, including𝒛i−𝒛j\\bm\{z\}\_\{i\}\-\\bm\{z\}\_\{j\}and a projection of encoder\-gated asymmetric pairwise statistics\. These features capture directional asymmetry between variables\. Each head produces adhd\_\{h\}\-dimensional hidden representation \(computed in float32 for numerical stability\)\. These are concatenated into a2​dh2d\_\{h\}\-dimensional vector and passed to a fusion MLP \(bottom\), which combines existence and directional evidence to produce a single logitℓi​jbase\\ell^\{\\mathrm\{base\}\}\_\{ij\}for the ordered pair\(i,j\)\(i,j\)\. The logit is clamped to±15\\pm 15before triangular refinement\.𝒁∈ℝB×D×dh\\bm\{Z\}\\in\\mathbb\{R\}^\{B\\times D\\times d\_\{h\}\}Learnable tokens𝑻∈ℝKc×dh\\bm\{T\}\\in\\mathbb\{R\}^\{K\_\{c\}\\times d\_\{h\}\},Kc=8K\_\{c\}\{=\}8\[𝒔\(9\)sym;\(𝒈⊙𝒔asym\)\(20\)\]\[\\,\\bm\{s\}^\{\\mathrm\{sym\}\}\_\{\(9\)\}\\;;\\;\(\\bm\{g\}\\\!\\odot\\\!\\bm\{s\}^\{\\mathrm\{asym\}\}\)\_\{\(20\)\}\\,\]Cross\-attn layer 1\(fp32\)Q=Q\{=\}tokens,K=V=𝒁K\{=\}V\{=\}\\bm\{Z\}\+\+GELU FFN \(2×2\\timeswidth\)Cross\-attn layer 2\(fp32\)clamp tokens±5000\\pm 5000between layersLoading networkSi​k=σ​\(LoadNet​\(\[𝒛i;𝒕k\]\)\)S\_\{ik\}=\\sigma\(\\mathrm\{LoadNet\}\(\[\\bm\{z\}\_\{i\};\\,\\bm\{t\}\_\{k\}\]\)\)conf\_stat\_projMLP​\(29→48→1\)\\mathrm\{MLP\}\(29\\to 48\\to 1\), bias init−1\-1Noisy\-OR\(fp32\)C^i​j=1−∏k=1Kc\(1−Si​k​Sj​k\)\\hat\{C\}\_\{ij\}=1\-\\prod\_\{k=1\}^\{K\_\{c\}\}\(1\-S\_\{ik\}S\_\{jk\}\)clamp\[10−3,1−10−3\]\[10^\{\-3\},\\,1\{\-\}10^\{\-3\}\]\+\+𝑪^∈\[0,1\]D×D\\hat\{\\bm\{C\}\}\\in\[0,1\]^\{D\\times D\}Figure 5:Confounder Module\.The model represents latent confounders usingKc=8K\_\{c\}=8learnable tokens, which cross\-attend to the variable embeddings𝒁∈ℝD×dh\\bm\{Z\}\\in\\mathbb\{R\}^\{D\\times d\_\{h\}\}over two layers\. Each layer consists of multi\-head attention followed by a GELU feedforward network \(with2×2\\timeshidden width\) and LayerNorm, with computations performed in float32 for numerical stability\. Activations are clamped to±5000\\pm 5000between layers\. A loading network maps each variable–token pair to a probabilitySi​k∈\[0,1\]S\_\{ik\}\\in\[0,1\], indicating whether confounderkkinfluences variableii\. Pairwise confounding probabilities are then computed via noisy\-OR aggregation, with computation in float32 and outputs clamped away from0and11\. In parallel, a pairwise feature projection maps symmetric and gated antisymmetric statistics to a scalar logit bias, which is added to the noisy\-OR logits\. The resulting confounding matrix𝑪^∈\[0,1\]D×D\\hat\{\\bm\{C\}\}\\in\[0,1\]^\{D\\times D\}is used both as an output and as input to the edge predictor, whereC^i​j\\hat\{C\}\_\{ij\}is detached to prevent edge\-prediction gradients from suppressing confounder learning\.
## Appendix FLoss Function Details

As introduced in \([1](https://arxiv.org/html/2606.17516#S3.E1)\), the training objective combines a primary directed\-edge binary cross\-entropy loss with a set of auxiliary terms that separately regularize directionality, skeleton recovery, latent confounding, calibration, sparsity, and representation diversity\. This decomposition reflects the fact that causal discovery contains several coupled but distinct subproblems: detecting whether two variables are adjacent, orienting the edge correctly, avoiding overconfident reverse predictions, accounting for hidden common causes, and controlling graph density\. In particular,ℒasym\\mathcal\{L\}\_\{\\mathrm\{asym\}\}andℒbdir\\mathcal\{L\}\_\{\\mathrm\{bdir\}\}provide explicit supervision for edge orientation,ℒskel\\mathcal\{L\}\_\{\\mathrm\{skel\}\}gives a direction\-agnostic signal for adjacency recovery, andℒconf\\mathcal\{L\}\_\{\\mathrm\{conf\}\}trains the confounder module to identify latent common causes\. The calibration, false\-positive, density, PMA\-diversity, and gate\-regularization terms stabilize training and improve transfer by discouraging pathological solutions such as saturated logits, overly dense graphs, collapsed pooling queries, or feature gates that over\-specialize to synthetic mechanisms\.

TermDefinitionWeight / settingℒasym\\mathcal\{L\}\_\{\\mathrm\{asym\}\}1n\+​∑\(i,j\):Di​j∗=1ReLU​\(σ​\(ℓj​i\)−δ\)​Vi​j\\displaystyle\\frac\{1\}\{n\_\{\+\}\}\\sum\_\{\(i,j\):D^\{\*\}\_\{ij\}=1\}\\mathrm\{ReLU\}\\\!\\left\(\\sigma\(\\ell\_\{ji\}\)\-\\delta\\right\)V\_\{ij\}λasym=0\.08\\lambda\_\{\\mathrm\{asym\}\}=0\.08,δ=0\.15\\delta=0\.15ℒcal\\mathcal\{L\}\_\{\\mathrm\{cal\}\}1n\+​∑\(i,j\):Di​j∗=1\[ReLU​\(\|ℓi​j−ℓj​i\|−Δmax\)\]2\\displaystyle\\frac\{1\}\{n\_\{\+\}\}\\sum\_\{\(i,j\):D^\{\*\}\_\{ij\}=1\}\\left\[\\mathrm\{ReLU\}\\\!\\left\(\|\\ell\_\{ij\}\-\\ell\_\{ji\}\|\-\\Delta\_\{\\max\}\\right\)\\right\]^\{2\}λcal=0\.1\\lambda\_\{\\mathrm\{cal\}\}=0\.1,Δmax=4\.0\\Delta\_\{\\max\}=4\.0ℒbdir\\mathcal\{L\}\_\{\\mathrm\{bdir\}\}1\|𝒜\|​∑\(i,j\)∈𝒜BCE​\(ℓi​j−ℓj​i,1\)\\displaystyle\\frac\{1\}\{\|\\mathcal\{A\}\|\}\\sum\_\{\(i,j\)\\in\\mathcal\{A\}\}\\mathrm\{BCE\}\(\\ell\_\{ij\}\-\\ell\_\{ji\},1\), where𝒜=\{\(i,j\):Di​j∗=1,ℓi​j\>−3​or​ℓj​i\>−3\}\\mathcal\{A\}=\\\{\(i,j\):D^\{\*\}\_\{ij\}=1,\\ \\ell\_\{ij\}\>\-3\\ \\text\{or\}\\ \\ell\_\{ji\}\>\-3\\\}λbdir=0\.25\\lambda\_\{\\mathrm\{bdir\}\}=0\.25ℒskel\\mathcal\{L\}\_\{\\mathrm\{skel\}\}1\|V△\|​∑i<jwi​jskel​Vi​j​BCE​\(skel​\_​fn​\(ℓi​j,ℓj​i\),skel​\_​fn​\(Di​j∗,Dj​i∗\)\)\\displaystyle\\frac\{1\}\{\|V\_\{\\triangle\}\|\}\\sum\_\{i<j\}w^\{\\mathrm\{skel\}\}\_\{ij\}V\_\{ij\}\\mathrm\{BCE\}\\\!\\left\(\\mathrm\{skel\\\_fn\}\(\\ell\_\{ij\},\\ell\_\{ji\}\),\\mathrm\{skel\\\_fn\}\(D^\{\*\}\_\{ij\},D^\{\*\}\_\{ji\}\)\\right\)λskel=0\.5\\lambda\_\{\\mathrm\{skel\}\}=0\.5; positive weight capped at2\.02\.0ℒconf\\mathcal\{L\}\_\{\\mathrm\{conf\}\}1B​∑b=1B1\|V\(b\)\|​∑i​jwi​j\(b\)​Vi​j\(b\)​BCE​\(C^i​j\(b\),Ci​j∗\)\\displaystyle\\frac\{1\}\{B\}\\sum\_\{b=1\}^\{B\}\\frac\{1\}\{\|V^\{\(b\)\}\|\}\\sum\_\{ij\}w^\{\(b\)\}\_\{ij\}V^\{\(b\)\}\_\{ij\}\\mathrm\{BCE\}\\\!\\left\(\\hat\{C\}^\{\(b\)\}\_\{ij\},C^\{\*\}\_\{ij\}\\right\)λconf=0\.2\\lambda\_\{\\mathrm\{conf\}\}=0\.2; per\-graph positive weight capped at5\.05\.0ℒfp\\mathcal\{L\}\_\{\\mathrm\{fp\}\}λfpeff\|V\|​∑i​jReLU​\(σ​\(ℓi​j\)​\(1−Di​j∗\)​Vi​j−0\.15\)\\displaystyle\\frac\{\\lambda\_\{\\mathrm\{fp\}\}^\{\\mathrm\{eff\}\}\}\{\|V\|\}\\sum\_\{ij\}\\mathrm\{ReLU\}\\\!\\left\(\\sigma\(\\ell\_\{ij\}\)\(1\-D^\{\*\}\_\{ij\}\)V\_\{ij\}\-0\.15\\right\), whereλfpeff=λfp​\[1\+3​clamp​\(p∗−p^t,−0\.3,0\.3\)\]\\lambda\_\{\\mathrm\{fp\}\}^\{\\mathrm\{eff\}\}=\\lambda\_\{\\mathrm\{fp\}\}\\left\[1\+3\\,\\mathrm\{clamp\}\(p^\{\*\}\-\\hat\{p\}\_\{t\},\-0\.3,0\.3\)\\right\]λfp=0\.15\\lambda\_\{\\mathrm\{fp\}\}=0\.15,p∗=0\.50p^\{\*\}=0\.50ℒden\\mathcal\{L\}\_\{\\mathrm\{den\}\}1B​∑b=1Bsoftplus​\(∑i​jσ​\(ℓi​j\(b\)\)​Vi​j\(b\)−deff\(b\)​e¯\)2​deff\(b\)\\displaystyle\\frac\{1\}\{B\}\\sum\_\{b=1\}^\{B\}\\frac\{\\mathrm\{softplus\}\\\!\\left\(\\sum\_\{ij\}\\sigma\(\\ell^\{\(b\)\}\_\{ij\}\)V^\{\(b\)\}\_\{ij\}\-d\_\{\\mathrm\{eff\}\}^\{\(b\)\}\\bar\{e\}\\right\)\}\{2d\_\{\\mathrm\{eff\}\}^\{\(b\)\}\}λden=0\.05\\lambda\_\{\\mathrm\{den\}\}=0\.05,e¯=max⁡\(2\.0,e^ema\)\\bar\{e\}=\\max\(2\.0,\\hat\{e\}\_\{\\mathrm\{ema\}\}\)ℒpma\\mathcal\{L\}\_\{\\mathrm\{pma\}\}1Kq​\(Kq−1\)/2​∑a<bReLU​\(cos⁡\(𝒒a,𝒒b\)−0\.5\)\\displaystyle\\frac\{1\}\{K\_\{q\}\(K\_\{q\}\-1\)/2\}\\sum\_\{a<b\}\\mathrm\{ReLU\}\\\!\\left\(\\cos\(\\bm\{q\}\_\{a\},\\bm\{q\}\_\{b\}\)\-0\.5\\right\)λpma=0\.01\\lambda\_\{\\mathrm\{pma\}\}=0\.01ℒgate\\mathcal\{L\}\_\{\\mathrm\{gate\}\}Soft penalty on gate values outside\[0\.1,0\.9\]\[0\.1,0\.9\]λgate=min⁡\(0\.5,max⁡\(0,0\.8−t^ema\)\)\\lambda\_\{\\mathrm\{gate\}\}=\\min\(0\.5,\\max\(0,0\.8\-\\hat\{t\}\_\{\\mathrm\{ema\}\}\)\)Table 7:Auxiliary loss terms\. All sums are masked by valid pairs unless stated otherwise\.The auxiliary losses in Table[7](https://arxiv.org/html/2606.17516#A6.T7)are weighted to provide targeted gradients without dominating the primary edge\-prediction objective\. Directional losses use dead zones or activation gates so that they focus on ambiguous or incorrectly oriented edges rather than indefinitely pushing already\-correct logits to extremes\. The skeleton and confounding losses use capped positive\-class weights to address class imbalance while preventing rare positive labels from overwhelming shared encoder representations\. The false\-positive and density penalties control graph sparsity at local and global levels, respectively, while the calibration penalty limits excessive logit gaps that can harm transfer under distribution shift\. Finally,ℒpma\\mathcal\{L\}\_\{\\mathrm\{pma\}\}encourages the attention\-pooling queries to specialize to distinct sample\-distribution patterns, andℒgate\\mathcal\{L\}\_\{\\mathrm\{gate\}\}keeps the statistic gate from saturating, preserving adaptive use of pairwise causal features across datasets\.

## Appendix GSynthetic Training Dataset Creation

FoundCauseis trained*exclusively on synthetic*structural causal models \(SCMs\); no real\-world dataset is ever used for training\. All 15 benchmarks reported in Section[5](https://arxiv.org/html/2606.17516#S5)and the Tübingen cause–effect pairs are strictly held out\. The training distribution is constructed to span the breadth of generative assumptions encountered in practice, namely, graph topology, functional mechanism, noise family, dimensionality, sample size, missingness, confounding, so that the learned amortized inference transfers zero\-shot to real data\.

##### Continuous Regeneration\.

Rather than keeping a fixed training set, the data pipeline regenerates∼\\sim12,000 fresh tasks \(10,000 multivariate plus 2,000 bivariate\) at the start of every epoch using a pool of CPU workers\. Each epoch therefore exposes the model to previously unseen SCMs; the model never sees the same dataset twice\. This continuous regeneration acts as the sole regulariser \(no dropout is used anywhere in the architecture\)\.

##### Task Definition\.

A training task is a tuple\(𝑿,𝑴,𝑫∗,𝑪∗,𝒗\)\(\\bm\{X\},\\,\\bm\{M\},\\,\\bm\{D\}^\{\*\},\\,\\bm\{C\}^\{\*\},\\,\\bm\{v\}\)where

- •𝑿∈ℝN×Dmax\\bm\{X\}\\in\\mathbb\{R\}^\{N\\times D\_\{\\max\}\}is the observational data matrix \(padded toDmax=65D\_\{\\max\}\{=\}65\);
- •𝑴∈\{0,1\}N×Dmax\\bm\{M\}\\in\\\{0,1\\\}^\{N\\times D\_\{\\max\}\}is the observation mask \(observed versus missing\);
- •𝑫∗∈\{0,1\}Dmax×Dmax\\bm\{D\}^\{\*\}\\in\\\{0,1\\\}^\{D\_\{\\max\}\\times D\_\{\\max\}\}is the ground\-truth DAG among*observed*variables \(entries outside the observed sub\-block are zero\);
- •𝑪∗∈\{0,1\}Dmax×Dmax\\bm\{C\}^\{\*\}\\in\\\{0,1\\\}^\{D\_\{\\max\}\\times D\_\{\\max\}\}is the ground\-truth confounding matrix \(symmetric\) indicating pairs sharing a hidden common cause through hidden\-only paths;
- •𝒗∈\{0,1\}Dmax\\bm\{v\}\\in\\\{0,1\\\}^\{D\_\{\\max\}\}is the variable\-validity mask\.

Per task, the number of observed variables is sampled uniformly from\[2,50\]\[2,50\], the number of samples from\[100,600\]\[100,600\], and the latent confounder fraction from\[0,0\.30\]\[0,\\,0\.30\]of the observed count\.

To improve robustness in low\-edge\-density regimes,8%8\\%of synthetic training tasks are generated as ultra\-sparse SCMs, with graph density sampled in the range\[0\.01,0\.04\]\[0\.01,0\.04\]\. In addition,5%5\\%of synthetic tasks are generated as discrete Bayesian networks: conditional probability tables are sampled from Dirichlet distributions, and observations are produced by ancestral sampling from the resulting directed graphical model\.

##### Two Generation Paths\.

Half of the tasks are produced by theDoWhypath described in Section[G\.1](https://arxiv.org/html/2606.17516#A7.SS1); the other half use a custom diverse\-graph generator supporting scale\-free \(Barabási–Albert\), small\-world \(Watts–Strogatz\), stochastic block model, geometric, and star/hub topologies, together with 18\+ nonlinear mechanism families \(GP via random Fourier features, tanh and sigmoid MLPs, Michaelis–Menten, Hill, piecewise\-linear, polynomial, step, interaction, competitive inhibition, sinusoidal, exponential, logarithmic, log\-space ANM, threshold\-binary, categorical\) and noise distributions including Gaussian, Laplace, uniform, Gaussian mixtures, Cauchy, and Student\-tt\. A separate post\-processing pipeline injects real\-world artefacts: MCAR / MNAR / MAR missingness \(15% of tasks\), measurement noise, zero\-inflation, discretization, rounding, heteroscedasticity, and log\-uniform scale randomization\. These anti\-shortcut perturbations are applied symmetrically to causes and effects so the model cannot recover direction from marginal\-distribution signatures\.

##### Anti\-shortcut Augmentation\.

At batch assembly time, variables are randomly permuted on 100% of training batches \(with consistent relabelling of𝑫∗\\bm\{D\}^\{\*\}and𝑪∗\\bm\{C\}^\{\*\}\), eliminating the variance\-ordering shortcut documented by Reisach et al\. \(2023\)\. Additionally, 5% of multivariate tasks are replaced with subpopulation\-mixture tasks \(samples drawn from two independent SCMs, labels set to the union of both DAGs\) to teach robustness against heterogeneous populations\. A 12% “correlation fog” mode replaces base tasks with near\-deterministic \(noise\_std∈\[0\.001,0\.05\]\\text\{noise\\\_std\}\\in\[0\.001,0\.05\]\) dense linear SCMs, forcing the model to distinguish direct edges from strong indirect correlation\.

##### Noise Models\.

To encourage robustness across heterogeneous data\-generating processes, synthetic SCMs are generated with a mixture of additive noise distributions\. At full training difficulty, Gaussian noise receives weight30%30\\%, while the remaining probability mass is split equally across Laplace, uniform, mixture, and Cauchy noise\. This design exposes the model to light\-tailed, heavy\-tailed, bounded, multimodal, and pathological noise regimes\. In particular, the inclusion of Laplace and Cauchy noise tests robustness to heavy\-tailed perturbations, uniform noise tests bounded\-support mechanisms, and mixture noise introduces multimodality that cannot be summarized by low\-order moments alone\. Cauchy noise is clipped at50×50\\timesthe nominal standard\-deviation scale used for sampling to avoid numerical instabilities while retaining its extreme\-tail behavior\.

Noise modelWeightPropertiesGaussian30%30\\%Standard additive noiseLaplace17\.5%17\.5\\%Heavy\-tailed noiseUniform17\.5%17\.5\\%Bounded\-support noiseMixture17\.5%17\.5\\%Multimodal noise with22–44componentsCauchy17\.5%17\.5\\%Undefined variance; clipped at50×50\\timesnominal scaleTable 8:Noise models used for synthetic SCM generation\.Weights correspond to the full\-difficulty training distribution\.
##### Post\-processing and Special Task Types\.

After mechanism generation, each synthetic dataset is passed through a stochastic post\-processing pipeline designed to mimic common artifacts in real observational data\. The pipeline injects autocorrelation with probability5%5\\%and batch effects with probability3%3\\%, both using per\-variable random splits; adds independent noise with probability5%5\\%, measurement noise with probability7\.5%7\.5\\%, and outliers with probability5%5\\%at magnitudes between33and88standard deviations, skipping these continuous perturbations for discrete variables; applies discretization with probability5%5\\%using55–100100bins; introduces bounded or positive transforms with probability4%4\\%; applies monotone post\-transforms such as log, square\-root, power, sigmoid, and softplus with probability3%3\\%; adds zero\-inflation or censoring with probability4%4\\%; applies rounding noise in the range0–15%15\\%; randomizes variable scales with probability10%10\\%; and introduces extreme scale heterogeneity with probability5%5\\%using multiplicative factors in10\[−3,3\]10^\{\[\-3,3\]\}\. To model latent dependence and near\-collinearity, the generator also introduces correlated exogenous noise with probability12%12\\%, updating the confounding matrix accordingly, and grouped collinearity with probability5%5\\%, where one source variable is copied into33–66noisy variants, again updating the confounding matrix\. Additional augmentations include Tübingen\-like bivariate transformations for35%35\\%ofd=2d=2tasks and log\-transform augmentation with probability6%6\\%to create positive, skewed, Sachs\-like data\. Discrete variables are automatically detected and protected from incompatible continuous\-noise perturbations\. Finally, the synthetic mixture includes several special task types: bivariate tasks comprise15%15\\%of training data, with10%10\\%easy and90%90\\%full\-range difficulty; correlation\-fog tasks comprise12%12\\%and use all\-linear mechanisms with very low noise in\[0\.001,0\.05\]\[0\.001,0\.05\]; subpopulation\-mixing tasks comprise5%5\\%and combine two SCMs using the union of their edge sets; and intervention\-mixture tasks comprise8%8\\%and pool observations across multiple experimental conditions\.

### G\.1TheDoWhyGenerator Path

Half of each training batch is generated using the random structural causal model \(SCM\) generator provided by DoWhy\[[47](https://arxiv.org/html/2606.17516#bib.bib23),[6](https://arxiv.org/html/2606.17516#bib.bib33)\]\. Data are sampled directly from these SCMs using the library’s reference implementation, ensuring consistency with the underlying data\-generation procedures\.

##### Graph Topology\.

DoWhy constructs a DAG by sequential random attachment, yielding graphs approximately following an Erdős–Rényi distribution with edge density drawn from\[0\.02,0\.40\]\[0\.02,0\.40\]per task\. The root fraction is sampled from\[0\.05,0\.40\]\[0\.05,0\.40\]\. A scale\-aware cap restricts expected in\-degree to≤4\\leq 4at largeDDviamin⁡\(ρ,4/\(D−1\)\)\\min\(\\rho,\\,4/\(D\-1\)\), keeping direct edges detectable atD=50D\{=\}50\.

##### Mechanism Types\.

For each non\-root nodeYY, the data\-generating process assigns a structural mechanism based on a fixed probability of selecting a linear model \(0\.300\.30\), with the remaining probability split evenly across two nonlinear families\. Let𝑿Pa​\(Y\)\\bm\{X\}\_\{\\mathrm\{Pa\}\(Y\)\}denote the parent variables ofYY\.

- •Linear\(30%\): Y=𝒘⊤​𝑿Pa​\(Y\)\+ε,Y=\\bm\{w\}^\{\\top\}\\bm\{X\}\_\{\\mathrm\{Pa\}\(Y\)\}\+\\varepsilon,where𝒘\\bm\{w\}is sampled from a Gaussian distribution andε\\varepsilonis independent noise\.
- •Nonlinear additive MLP\(35%\): Y=fNN​\(𝑿Pa​\(Y\)\)\+ε,Y=f\_\{\\mathrm\{NN\}\}\(\\bm\{X\}\_\{\\mathrm\{Pa\}\(Y\)\}\)\+\\varepsilon,wherefNNf\_\{\\mathrm\{NN\}\}is a multilayer perceptron \(2–4 layers, 4–64 hidden units per layer\) with spectral normalization applied to the weights and a scaling factor to induce strong nonlinearity\.
- •Non\-additive noise MLP\(35%\): Y=fNN​\(\[𝑿Pa​\(Y\);ε\]\),Y=f\_\{\\mathrm\{NN\}\}\(\[\\bm\{X\}\_\{\\mathrm\{Pa\}\(Y\)\};\\varepsilon\]\),where the noiseε\\varepsilonis concatenated with the parent variables before being processed by the network, resulting in non\-additive noise effects\.

To stabilize the range of generated values, we perform an initial calibration step by drawing1,0001\{,\}000samples from the SCM and recording per\-variable output ranges\. Subsequently, outputs of neural mechanisms are rescaled to lie within\[−2,2\]\[\-2,2\], preventing value explosion as signals propagate through the graph\.

##### Noise Distributions\.

DoWhy supports Gaussian, Laplace, uniform, and Student\-ttadditive noise \(weighted0\.30/0\.30/0\.20/0\.200\.30/0\.30/0\.20/0\.20respectively\) with log\-uniform scale drawn from\[0\.01,0\.25\]\[0\.01,\\,0\.25\]\. Root variables additionally sample from log\-normal, exponential,β\\beta, andχ2\\chi^\{2\}distributions with65%65\\%of root marginals replaced by a22–66component Gaussian mixture \(componentstd=0\.4\\mathrm\{std\}\{=\}0\.4, means in\[−3,3\]\[\-3,3\]\) to ensure multimodal root behaviour\.

##### Hidden Confounders\.

To induce latent\-confounding structure, DoWhy generates an extended SCM withD\+KHD\+K\_\{H\}variables whereKHK\_\{H\}equals the requested latent count; a weighted selection favouring nodes with≥2\\geq 2children is then applied to hide variables\. The ground\-truth DAG𝑫∗\\bm\{D\}^\{\*\}over observed variables is computed as the transitive closure through hidden\-only paths \(so an observediicauses observedjjwhenever there is a directed pathi→⋯→ji\\to\\cdots\\to jpassing only through hidden nodes\), and𝑪i​j∗=1\\bm\{C\}^\{\*\}\_\{ij\}\{=\}1whenever observediiandjjshare a hidden common ancestor reachable only through hidden nodes\. Both𝑫∗\\bm\{D\}^\{\*\}and𝑪∗\\bm\{C\}^\{\*\}are therefore correctly labelled under the marginalised\-observed process, not the full latent process\.

##### Sanitization\.

Deep nonlinear chains at largeDDoccasionally produce NaN or Inf values; any such entries are replaced with zero and marked as missing in𝑴\\bm\{M\}\. A post\-sanitization guard enforcesmax⁡\(10,0\.10​N\)\\max\(10,\\,0\.10N\)observed samples per variable to prevent degenerate columns that would cause all\-masked softmaxes downstream\.

##### Configuration\.

The full DoWhy configuration is summarised in Table[9](https://arxiv.org/html/2606.17516#A7.T9)\.

ParameterValue / range\|𝒱obs\|\|\\mathcal\{V\}\_\{\\mathrm\{obs\}\}\|\(observed variables\)uniform\[2,50\]\[2,50\]\|𝒱hidden\|/\|𝒱obs\|\|\\mathcal\{V\}\_\{\\mathrm\{hidden\}\}\|/\|\\mathcal\{V\}\_\{\\mathrm\{obs\}\}\|\(latent frac\.\)\[0\.00,0\.30\]\[0\.00,0\.30\]Root fraction\[0\.05,0\.40\]\[0\.05,0\.40\]Edge density\[0\.02,0\.40\]\[0\.02,0\.40\], capped at4/\(D−1\)4/\(D\-1\)Number of samplesNNuniform\[100,600\]\[100,600\]Mechanism: linear / tanh\-NN / non\-additive30%/35%/35%30\\%/35\\%/35\\%NN hidden layers\[2,4\]\[2,4\]NN hidden units/layer\[4,64\]\[4,64\]NN weight scale \(spectral\-normed\)5\.05\.0NN output normalization\[−2,2\]\[\-2,2\]Noise std \(log\-uniform\)\[0\.01,0\.25\]\[0\.01,0\.25\]Noise: Gaussian / Laplace / uniform /tt30%/30%/20%/20%30\\%/30\\%/20\\%/20\\%Root mixture fraction65%65\\%\(22–66Gaussian components\)Prob\. missing data injection15%15\\%of tasks \(2–12% of cells\)Table 9:DoWhygenerator configuration used for 50% of training tasks\.

## Appendix HPer\-Dataset Results

Table 10:Per\-dataset results for GRaSP\.GRaSP is the strongest classical baseline overall, with average directed\-edgeF1=0\.488F\_\{1\}=0\.488across the 15 benchmark datasets\.DatasetF1↑F\_\{1\}\\uparrowSHD↓\\downarrowFoundOrientedExtraalarm\_like0\.8220\.822191946/4646/4644/4644/4644asia0\.8420\.842338/88/88/88/80asia\_nonlinear0\.7140\.714445/85/85/55/50causal\_chambers0\.3890\.389444422/3922/3914/2214/2299child0\.7540\.754151525/2525/2523/2523/2577child\_nonlinear0\.2380\.238323212/2512/255/125/1233ecoli\_like0\.8990\.899131358/5858/5858/5858/580insurance0\.3550\.355808047/5247/5222/4722/471919petshop\_high0\.3930\.393686827/4227/4222/2722/273434petshop\_low0\.3550\.355696922/4322/4319/2219/223939petshop\_temp10\.2960\.296767625/4325/4316/2516/254040petshop\_temp20\.3390\.339747424/4324/4319/2419/244343sachs\_full0\.3640\.36421218/178/176/86/866sachs\_nonlinear0\.2220\.22221215/175/173/53/544sachs\_obs0\.3450\.34519197/177/175/75/70Average0\.4880\.488––––Table 11:Per\-dataset results for GES\.GES completes all 15 benchmark datasets and obtains average directed\-edgeF1=0\.456F\_\{1\}=0\.456\.DatasetF1↑F\_\{1\}\\uparrowSHD↓\\downarrowFoundOrientedExtraalarm\_like0\.8070\.807212145/4645/4644/4544/4555asia0\.8420\.842338/88/88/88/80asia\_nonlinear0\.7140\.714445/85/85/55/50causal\_chambers0\.5000\.500343417/3917/3917/1717/171111child0\.4000\.400424221/2521/2514/2114/211818child\_nonlinear0\.3330\.333282813/2513/257/137/1322ecoli\_like0\.8090\.809252557/5857/5853/5753/5755insurance0\.2390\.23910210239/5239/5216/3916/394141petshop\_high0\.3900\.390727227/4227/4223/2723/274242petshop\_low0\.3550\.355696924/4324/4319/2419/243737petshop\_temp10\.2310\.231808022/4322/4312/2212/223939petshop\_temp20\.2310\.231808022/4322/4312/2212/223939sachs\_full0\.3430\.34323237/177/176/76/799sachs\_nonlinear0\.2960\.29619196/176/174/64/633sachs\_obs0\.3450\.34519197/177/175/75/70Average0\.4560\.456––––Table 12:Per\-dataset results for PC withα=0\.01\\alpha=0\.01\.PC completes13/1513/15benchmark datasets and obtains average directed\-edgeF1=0\.464F\_\{1\}=0\.464over completed runs\. The two temporal Petshop variants time out after four days and are excluded from the average\.DatasetF1↑F\_\{1\}\\uparrowSHD↓\\downarrowFoundOrientedExtraalarm\_like0\.6670\.667333337/4637/4633/3733/3722asia0\.8420\.842338/88/88/88/80asia\_nonlinear0\.5710\.571665/85/84/54/50causal\_chambers0\.2140\.21444446/396/396/66/688child0\.7170\.717151520/2520/2519/2019/200child\_nonlinear0\.4350\.435262613/2513/2510/1310/1311ecoli\_like0\.5950\.595494943/5843/5836/4336/4355insurance0\.3060\.306595919/5219/5213/1913/1999petshop\_high0\.3820\.382555520/4220/4217/2017/202121petshop\_low0\.3110\.311626219/4319/4314/1914/192323petshop\_temp1N/A; timeout after four dayspetshop\_temp2N/A; timeout after four dayssachs\_full0\.4520\.45217177/177/177/77/733sachs\_nonlinear0\.2580\.25823237/177/174/74/744sachs\_obs0\.2760\.27621217/177/174/74/711Average0\.4640\.464––––Table 13:Per\-dataset results forFoundCauseon the real\-world benchmark suite\.We report directed\-edgeF1F\_\{1\}, skeletonF1F\_\{1\}, structural Hamming distance \(SHD\), normalized SHD, number of true edges found, number of found edges oriented correctly, number of extra predicted edges, and directed\-edge AUROC\. The final row reports macro\-averages across datasets where applicable\.DatasetddnnF1↑F\_\{1\}\\uparrowSkel\-F1↑F\_\{1\}\\uparrowSHD↓\\downarrowSHD/d↓/d\\downarrowFoundOrientedExtraAUROC↑\\uparrowalarm\_like3737200020000\.8330\.8330\.8330\.83318180\.490\.4945/4645/4645/4545/4517170\.9820\.982asia88100010001\.0001\.0001\.0001\.00000\.000\.008/88/88/88/800\.9790\.979asia\_nonlinear88100010000\.9330\.9330\.9330\.933110\.120\.127/87/87/77/700\.9380\.938causal\_chambers202010000100000\.5880\.5880\.5880\.58828281\.401\.4020/3920/3920/2020/20990\.7400\.740child2020100010000\.6030\.6030\.7620\.76220201\.001\.0024/2524/2519/2419/2414140\.8430\.843child\_nonlinear2020100010000\.5710\.5710\.6670\.66716160\.800\.8014/2514/2512/1412/14330\.7280\.728ecoli\_like4646100010000\.9050\.9050\.9050\.90512120\.260\.2657/5857/5857/5757/5711110\.9880\.988insurance2727150015000\.4030\.4030\.6050\.60559592\.192\.1936/5236/5224/3624/3631310\.6910\.691petshop\_high39395895890\.2340\.2340\.4940\.49449491\.261\.2619/4219/429/199/1916160\.6560\.656petshop\_low41415895890\.1840\.1840\.3680\.36855551\.341\.3414/4314/437/147/1419190\.6070\.607petshop\_temp14141165216520\.2930\.2930\.5070\.50745451\.101\.1019/4319/4311/1911/1913130\.6660\.666petshop\_temp24141165216520\.3120\.3120\.5190\.51945451\.101\.1020/4320/4312/2012/2014140\.6660\.666sachs\_full1111746674660\.4800\.4800\.4800\.48013131\.181\.186/176/176/66/6220\.6500\.650sachs\_nonlinear1111100010000\.3570\.3570\.6430\.64314141\.271\.279/179/175/95/9220\.6310\.631sachs\_obs11118538530\.3200\.3200\.6400\.64013131\.181\.188/178/174/84/800\.5750\.575Average––0\.5350\.5350\.6630\.66325\.925\.90\.980\.98–––0\.7560\.756

## Appendix ITest\-Time Ablation Study

To isolate the contribution of individual architectural components, we perform test\-time ablations on the fully trained Causal Discovery Transformer \(CDT\) without any retraining\. Our protocol replaces specific intermediate activations with neutral values \(zero tensors or learned constants\) at inference, while keeping every other weight, dataflow, and precision setting unchanged\. Because the forward pass is preserved end\-to\-end, each ablation produces a valid DAG prediction𝑫^∈\[0,1\]D×D\\hat\{\\bm\{D\}\}\\in\[0,1\]^\{D\\times D\}; only the*information content*of a chosen component is suppressed\.

##### Scope and Caveats\.

Test\-time ablation measures how much the trained model*currently relies on*a component, not whether a freshly trained model*could have learned to work without it*\. Removing a component the model was trained to use introduces a distribution shift in its downstream inputs, so the measured drop is an upper bound on true necessity\. We mark each ablation as eitherclean\(the architecture was trained to support the ablated state, or the ablated value is exactly in\-distribution for downstream modules\) orshift\(the ablation produces activation statistics the model never saw during training\)\.

### I\.1Experimental Protocol

##### Datasets\.

For the purposes of Ablation study, we evaluate on five real\-world causal discovery benchmarks \(asia,child\_nonlinear,alarm\_like,sachs\_obsandcausal\_chambers\_lt\) plus the Tübingen cause–effect pairs described in Appendix[C](https://arxiv.org/html/2606.17516#A3)\. As mentioned before, none of these datasets appear in training; the model was trained exclusively on synthetic structural causal models generated on\-the\-fly\.

##### Implementation\.

Each ablation modifies a single component of the trained model \(e\.g\., a module or parameter tensor\) while leaving the remainder of the architecture unchanged\. For each run, the model is evaluated with the ablation applied, after which the original model state is restored before the next evaluation\.

To verify correct restoration, we compare the predicted edge probabilities𝑫^∈\[0,1\]D×D\\hat\{\\bm\{D\}\}\\in\[0,1\]^\{D\\times D\}on a fixed synthetic validation dataset before and after each ablation, ensuring they match the baseline predictions\. All evaluations \(one baseline and six ablations\) are performed using the same checkpoint and within a single execution environment, ensuring that differences in performance are attributable solely to the ablated component\.

##### Checkpoint\.

All ablations use a single trained model checkpoint corresponding to epoch 320 \(approximately139139M parameters\), selected based on the best skeletonF1F\_\{1\}score on a held\-out synthetic validation set\. The model is trained using Schedule\-Free AdamW, which maintains both training iterates and an exponential moving average of parameters\. We use the averaged \(evaluation\) parameters for all experiments\. Implementation\-specific naming artifacts from compilation are removed when loading the checkpoint, but do not affect the model architecture or weights\.

##### Inference Protocol\.

All ablations are evaluated using a fixed inference procedure\. Given an input dataset𝑿∈ℝN×D\\bm\{X\}\\in\\mathbb\{R\}^\{N\\times D\}, we performK=10K=10stochastic inference passes, each consisting of bootstrap resampling of rows and a random permutation of variable indices\. Predictions are mapped back to the original variable order via the inverse permutation, and the resulting logits are averaged across theKKruns\.

The averaged logits are calibrated using temperature scaling with temperatureT=0\.65T=0\.65, estimated on a held\-out synthetic validation set\. A threshold for edge prediction is then determined by fitting a two\-component Gaussian mixture model \(GMM\) to the logit distribution\. Finally, self\-consistency pruning is applied by retaining only edges that appear in at least50%50\\%of theKKruns\.

For computational efficiency, datasets with more than50005000samples are subsampled along the sample dimension using a fixed random seed\. All inference computations are performed without gradient tracking using mixed precision \(bfloat16\), with selected numerically sensitive operations executed in float32 for stability\.

### I\.2Ablations

We describe each ablation by specifying: \(i\) the component being modified, \(ii\) the replacement value or tensor, \(iii\) the downstream modules affected, and \(iv\) the ablation category\.

Let𝒁∈ℝD×dh\\bm\{Z\}\\in\\mathbb\{R\}^\{D\\times d\_\{h\}\}denote the per\-variable representations obtained after encoder processing and pooling,𝑺∈ℝD×D×ds\\bm\{S\}\\in\\mathbb\{R\}^\{D\\times D\\times d\_\{s\}\}the tensor of pairwise statistics withds=45d\_\{s\}=45, and𝑪^∈\[0,1\]D×D\\hat\{\\bm\{C\}\}\\in\[0,1\]^\{D\\times D\}the predicted confounding matrix obtained via noisy\-OR aggregation\.

The pairwise statistics𝑺\\bm\{S\}include symmetric features, antisymmetric features, V\-structure indicators, and reliability metadata, which together capture dependence, directional asymmetry, and statistical confidence\.

#### I\.2\.1Pairwise Statistics Disabled \(shift\)

Intervention\.We ablate the pairwise statistics by replacing the tensor𝑺∈ℝD×D×ds\\bm\{S\}\\in\\mathbb\{R\}^\{D\\times D\\times d\_\{s\}\}with zeros:

𝑺←𝟎,ds=45\.\\bm\{S\}\\leftarrow\\bm\{0\},\\qquad d\_\{s\}=45\.All subsequent computations that depend on𝑺\\bm\{S\}receive this zero input\.

Propagation\.Setting𝑺=𝟎\\bm\{S\}=\\bm\{0\}affects all modules that use pairwise statistics:

- •*Stat\-conditioned attention bias\.*The per\-layer projection from𝒔i​j\\bm\{s\}\_\{ij\}to attention bias reduces to a constant \(bias\-only\) term\. Since this constant is added uniformly across attention logits, it cancels under the softmax, resulting in effectively zero statistical bias\.
- •*Existence head\.*The projection of symmetric statistics receives𝟎\\bm\{0\}, producing a constant feature that does not depend on\(i,j\)\(i,j\)\.
- •*Direction head\.*Antisymmetric statistics are zero, so their projection is also constant\. Directional predictions therefore rely only on embedding differences𝒛i−𝒛j\\bm\{z\}\_\{i\}\-\\bm\{z\}\_\{j\}\.
- •*Feature gate\.*The gating network receives no statistical input, and depends solely on the encoder representations\.
- •*Confounder module\.*The pairwise feature projection used to bias confounding logits becomes constant across pairs\.

What remains active\.The encoder representations𝒁∈ℝD×dh\\bm\{Z\}\\in\\mathbb\{R\}^\{D\\times d\_\{h\}\}, the parent–child role score, the confounder token attention mechanism, the triangular refinement module, and the fusion MLP all remain unchanged\.

Cleanness status:shift\.The model is trained with informative pairwise statistics; replacing𝑺\\bm\{S\}with zeros induces a distribution shift in intermediate representations, as downstream modules receive inputs outside their training regime\.

#### I\.2\.2Triangular Edge Refinement Disabled \(clean\)

Intervention\.We disable the triangular refinement module by bypassing all refinement updates\. The model directly uses the base edge logitsℓbase∈ℝD×D\\bm\{\\ell\}^\{\\mathrm\{base\}\}\\in\\mathbb\{R\}^\{D\\times D\}produced by the edge predictor\.

Propagation\.Without refinement, the model skips theR=3R=3rounds of triangle\-based updates that incorporate higher\-order interactions across variable triples\. In particular, no information is aggregated over intermediate variables, and the learned blending between base and refined logits is not applied\. The final logits are therefore given by

ℓ=ℓbase\+λr​𝑹,\\bm\{\\ell\}=\\bm\{\\ell\}^\{\\mathrm\{base\}\}\+\\lambda\_\{r\}\\bm\{R\},where𝑹\\bm\{R\}denotes the parent–child role score andλr=0\.1\\lambda\_\{r\}=0\.1is its fixed weight\. The logits are subsequently clamped to the range±15\\pm 15before conversion to probabilities\.

What remains active\.The encoder representations𝒁\\bm\{Z\}, the factored edge predictor \(existence and direction heads\), the parent–child role score, the confounder module, and the fusion MLP remain unchanged\.

Cleanness status:clean\.The base predictor is trained jointly with the triangular module, and disabling refinement corresponds to relying solely on the base logits\. This setting is consistent with the training regime, as the model learns to combine base and refined predictions through blending weights, and includes configurations where the contribution of the refinement module is reduced\.

#### I\.2\.3PMA Pooling Replaced by Max Pooling \(clean\)

Intervention\.We replace attention\-based pooling \(PMA\) with max pooling over the sample dimension\. Let𝒛iPMA\\bm\{z\}\_\{i\}^\{\\mathrm\{PMA\}\}and𝒛imax\\bm\{z\}\_\{i\}^\{\\max\}denote the PMA and max\-pooled representations for variableii, respectively\. The original model combines these via a learned gateg∈ℝg\\in\\mathbb\{R\}:

𝒛i=σ​\(g\)​𝒛iPMA\+\(1−σ​\(g\)\)​𝒛imax\.\\bm\{z\}\_\{i\}=\\sigma\(g\)\\,\\bm\{z\}\_\{i\}^\{\\mathrm\{PMA\}\}\+\\bigl\(1\-\\sigma\(g\)\\bigr\)\\,\\bm\{z\}\_\{i\}^\{\\max\}\.In this ablation, we setσ​\(g\)≈0\\sigma\(g\)\\approx 0, yielding

𝒛i≈𝒛imax\.\\bm\{z\}\_\{i\}\\approx\\bm\{z\}\_\{i\}^\{\\max\}\.
Propagation\.The PMA sub\-network still computes𝒛iPMA\\bm\{z\}\_\{i\}^\{\\mathrm\{PMA\}\}, but its contribution is negligible due to the gating\. As a result, the encoder output consists effectively of max\-pooled representations over samples, removing the model’s ability to learn adaptive, attention\-based aggregation\.

What remains active\.All other components, including the encoder, edge predictor, confounder module, and triangular refinement, remain unchanged\.

Cleanness status:clean\.The gating parameter is learned during training and operates continuously in\[0,1\]\[0,1\]\. In the trained model,σ​\(g\)≈0\.03\\sigma\(g\)\\approx 0\.03, indicating that max pooling already dominates the aggregation\. Settingσ​\(g\)≈0\\sigma\(g\)\\approx 0therefore corresponds to a small extrapolation from the training regime\.

#### I\.2\.4Confounder Feedback Disabled \(clean\)

Intervention\.We disable the use of confounder predictions in edge prediction by setting the confounding matrix to zero:

𝑪^←𝟎,𝑪^∈\[0,1\]D×D\.\\hat\{\\bm\{C\}\}\\leftarrow\\bm\{0\},\\qquad\\hat\{\\bm\{C\}\}\\in\[0,1\]^\{D\\times D\}\.The internal confounder representations \(e\.g\., variable–token loadings\) are still computed, but are not used by downstream modules\.

Propagation\.The existence head of the edge predictor receivesC^i​j\\hat\{C\}\_\{ij\}as an input feature\. With𝑪^=𝟎\\hat\{\\bm\{C\}\}=\\bm\{0\}, this feature channel is identically zero for all pairs\(i,j\)\(i,j\), removing the direct influence of inferred confounding on edge existence\. All other input features and computations in the edge predictor remain unchanged\.

What remains active\.The encoder representations𝒁\\bm\{Z\}, the factored edge predictor \(excluding the confounder input\), the parent–child role score, and the triangular refinement module remain unchanged\. In particular, higher\-order structural patterns such as forks and colliders can still be captured indirectly through pairwise and triangular interactions\.

Cleanness status:clean\.The model is trained to handle the absence of confounder input, as the training pipeline includes configurations where𝑪^\\hat\{\\bm\{C\}\}is zero\. This ablation therefore corresponds to a setting within the training regime, isolating the contribution of confounder feedback without introducing distribution shift\.

#### I\.2\.5Encoder\-conditioned Feature Gate Replaced by Constant0\.50\.5\(shift\)

Intervention\.We replace the encoder\-conditioned feature gate with a constant vector\. In the original model, the gate𝒈i​j∈\[0,1\]da\\bm\{g\}\_\{ij\}\\in\[0,1\]^\{d\_\{a\}\}\(withda=20d\_\{a\}=20\) is computed from the variable embeddings𝒛i,𝒛j∈ℝdh\\bm\{z\}\_\{i\},\\bm\{z\}\_\{j\}\\in\\mathbb\{R\}^\{d\_\{h\}\}and pairwise statistics𝒔i​j∈ℝds\\bm\{s\}\_\{ij\}\\in\\mathbb\{R\}^\{d\_\{s\}\}as

𝒈i​j=σ​\(𝑾2​GELU⁡\(𝑾1​\[𝒛i−𝒛j;𝒔i​j\]\)\)\.\\bm\{g\}\_\{ij\}=\\sigma\\\!\\left\(\\bm\{W\}\_\{2\}\\,\\operatorname\{GELU\}\\\!\\left\(\\bm\{W\}\_\{1\}\[\\bm\{z\}\_\{i\}\-\\bm\{z\}\_\{j\};\\bm\{s\}\_\{ij\}\]\\right\)\\right\)\.In this ablation, we set

𝒈i​j=12​1da,\\bm\{g\}\_\{ij\}=\\tfrac\{1\}\{2\}\\,\\bm\{1\}\_\{d\_\{a\}\},removing dependence on both embeddings and statistics\.

Propagation\.The gate modulates antisymmetric features used by the direction head and confounder module\. Replacing𝒈i​j\\bm\{g\}\_\{ij\}with a constant scales all antisymmetric features uniformly by a factor of0\.50\.5, eliminating pair\-specific feature selection\. Downstream projections and predictors continue to operate on these uniformly scaled features\.

What remains active\.All other components, including the encoder representations, symmetric features, factored edge predictor, confounder module, and triangular refinement, remain unchanged\.

Cleanness status:shift\.In the trained model, the gate produces a non\-uniform distribution of values across pairs, reflecting data\-dependent feature selection\. Replacing it with a constant removes this adaptive behavior and introduces a distribution shift in the inputs to downstream modules\.

#### I\.2\.6Direction\-head Activations Zeroed \(shift\)

Intervention\.We ablate the direction head by replacing its hidden activation with zeros prior to fusion\. Let𝒉i​jexist,𝒉i​jdir∈ℝdh\\bm\{h\}^\{\\mathrm\{exist\}\}\_\{ij\},\\bm\{h\}^\{\\mathrm\{dir\}\}\_\{ij\}\\in\\mathbb\{R\}^\{d\_\{h\}\}denote the existence and direction representations for a pair\(i,j\)\(i,j\)\. In this ablation,

𝒉i​jdir←𝟎,ℓi​j=fusion​\(\[𝒉i​jexist;0\]\),\\bm\{h\}^\{\\mathrm\{dir\}\}\_\{ij\}\\leftarrow\\bm\{0\},\\qquad\\bm\{\\ell\}\_\{ij\}=\\mathrm\{fusion\}\\\!\\bigl\(\[\\bm\{h\}^\{\\mathrm\{exist\}\}\_\{ij\};\\,\\bm\{0\}\]\\bigr\),whereℓi​j\\bm\{\\ell\}\_\{ij\}is the resulting base logit\.

Propagation\.The fusion module receives only existence\-based features, as the contribution of the direction head is removed\. Consequently, the learned interactions between symmetric \(existence\) and antisymmetric \(direction\) signals are lost, and predictions rely solely on existence\-related information within the fusion pathway\.

What remains active\.Other sources of directional information remain available, including the parent–child role score and the triangular refinement module, which can still capture asymmetric patterns such as colliders and chains\. All upstream components, including the encoder and existence head, are unchanged\.

Cleanness status:shift\.During training, the fusion module receives both existence and direction representations jointly\. Replacing the direction input with zeros produces inputs outside this training distribution, introducing a shift\. However, alternative pathways for direction prediction remain active, so orientation is degraded but not eliminated\.

### I\.3Results from Ablations

Cleanablations target components for which the architecture was trained to support a zero or bypassed path; the resulting numbers reflect the contribution of the component in isolation\.Shiftablations suppress a signal the model was trained to rely on and therefore provide an*upper bound*on the component’s necessity: a large drop could indicate either that the component is genuinely important or that the downstream layers are operating out of distribution\. We report results for both classes but interpret them accordingly\.

##### Interpretation\.

The ablation results in Table[14](https://arxiv.org/html/2606.17516#A9.T14)indicate that the largest gains arise from components that encode explicit causal structure\. Disabling the pairwise statistics leads to the most significant degradation, reducing averageF1F\_\{1\}from0\.6040\.604to0\.380\.38, AUROC from0\.8760\.876to0\.750\.75, and Tübingen accuracy from63\.7%63\.7\\%to53%53\\%\. This confirms that the symmetric, antisymmetric, and collider\-oriented statistics provide a primary source of causal inductive bias rather than acting as auxiliary features\.

Ablating the direction head also causes a substantial drop, lowering averageF1F\_\{1\}to0\.470\.47and Tübingen accuracy to56%56\\%, while AUROC remains relatively high at0\.860\.86\. This suggests that the model retains the ability to rank likely adjacent pairs but loses much of its capacity to correctly orient edges when the dedicated direction pathway is removed\.

Other components contribute more targeted improvements\. Removing triangular refinement reduces averageF1F\_\{1\}from0\.6040\.604to0\.550\.55, with pronounced effects onchild\_nl,alarm\_like, andcausal\_chambers\_lt, supporting its role in capturing higher\-order structures such as chains and colliders\. Replacing the encoder\-conditioned feature gate with a constant also degrades performance, indicating that adaptive feature selection is important for modulating the reliability of pairwise statistics\. In contrast, replacing PMA pooling with max pooling produces negligible changes, suggesting that the learned model operates close to the max\-pooling regime or that the benchmark tasks place limited demand on adaptive sample aggregation\.

Disabling confounder feedback yields a modest but consistent drop, particularly onsachsandcausal\_chambers\_lt, indicating that explicit modeling of latent confounding is most beneficial in settings with hidden common causes or complex dependencies\. Overall, the ablations support the central architectural claim: performance is driven primarily by explicit pairwise causal statistics, a dedicated direction pathway, and graph\-level refinement, while pooling and confounder feedback provide secondary but stabilizing contributions\.

Configasiachild\_nlalarmsachscc\_ltAvgF1F\_\{1\}Avg AUROCTüb\. %Baseline \(full model\)0\.8750\.875/0\.9840\.9840\.4550\.455/0\.8750\.8750\.8220\.822/0\.9990\.9990\.2220\.222/0\.6800\.6800\.6480\.648/0\.8410\.8410\.6040\.6040\.8760\.87663\.763\.7Pairwise stats disabled \(set→\\to𝟎\)\\mathbf\{0\}\)0\.55​\(−0\.33\)0\.55\\ \(\\\!\-\\\!0\.33\)/0\.840\.840\.28​\(−0\.18\)0\.28\\ \(\\\!\-\\\!0\.18\)/0\.720\.720\.60​\(−0\.22\)0\.60\\ \(\\\!\-\\\!0\.22\)/0\.930\.930\.10​\(−0\.12\)0\.10\\ \(\\\!\-\\\!0\.12\)/0\.550\.550\.35​\(−0\.30\)0\.35\\ \(\\\!\-\\\!0\.30\)/0\.700\.700\.38​\(−0\.22\)0\.38\\ \(\\\!\-\\\!0\.22\)0\.75​\(−0\.13\)0\.75\\ \(\\\!\-\\\!0\.13\)53​\(−11\)53\\ \(\\\!\-\\\!11\)Triangular edge refinement disabled0\.85​\(−0\.02\)0\.85\\ \(\\\!\-\\\!0\.02\)/0\.980\.980\.38​\(−0\.07\)0\.38\\ \(\\\!\-\\\!0\.07\)/0\.850\.850\.76​\(−0\.06\)0\.76\\ \(\\\!\-\\\!0\.06\)/0\.990\.990\.20​\(−0\.02\)0\.20\\ \(\\\!\-\\\!0\.02\)/0\.670\.670\.58​\(−0\.07\)0\.58\\ \(\\\!\-\\\!0\.07\)/0\.820\.820\.55​\(−0\.05\)0\.55\\ \(\\\!\-\\\!0\.05\)0\.86​\(−0\.02\)0\.86\\ \(\\\!\-\\\!0\.02\)63\.0​\(−0\.7\)63\.0\\ \(\\\!\-\\\!0\.7\)PMA Pooling replaced by Max\-Pool0\.87​\(−0\.00\)0\.87\\ \(\\\!\-\\\!0\.00\)/0\.980\.980\.45​\(−0\.00\)0\.45\\ \(\\\!\-\\\!0\.00\)/0\.870\.870\.82​\(−0\.00\)0\.82\\ \(\\\!\-\\\!0\.00\)/0\.990\.990\.22​\(−0\.00\)0\.22\\ \(\\\!\-\\\!0\.00\)/0\.680\.680\.64​\(−0\.01\)0\.64\\ \(\\\!\-\\\!0\.01\)/0\.840\.840\.60​\(−0\.00\)0\.60\\ \(\\\!\-\\\!0\.00\)0\.87​\(−0\.00\)0\.87\\ \(\\\!\-\\\!0\.00\)63\.7​\(−0\.0\)63\.7\\ \(\\\!\-\\\!0\.0\)Confounder feedback disabled0\.87​\(−0\.00\)0\.87\\ \(\\\!\-\\\!0\.00\)/0\.980\.980\.44​\(−0\.02\)0\.44\\ \(\\\!\-\\\!0\.02\)/0\.870\.870\.82​\(−0\.00\)0\.82\\ \(\\\!\-\\\!0\.00\)/0\.990\.990\.19​\(−0\.03\)0\.19\\ \(\\\!\-\\\!0\.03\)/0\.660\.660\.62​\(−0\.03\)0\.62\\ \(\\\!\-\\\!0\.03\)/0\.830\.830\.59​\(−0\.02\)0\.59\\ \(\\\!\-\\\!0\.02\)0\.87​\(−0\.01\)0\.87\\ \(\\\!\-\\\!0\.01\)63\.7​\(−0\.0\)63\.7\\ \(\\\!\-\\\!0\.0\)Encoder\-conditioned Geature Gatereplaced by 0\.50\.83​\(−0\.05\)0\.83\\ \(\\\!\-\\\!0\.05\)/0\.970\.970\.42​\(−0\.03\)0\.42\\ \(\\\!\-\\\!0\.03\)/0\.860\.860\.80​\(−0\.02\)0\.80\\ \(\\\!\-\\\!0\.02\)/0\.990\.990\.20​\(−0\.02\)0\.20\\ \(\\\!\-\\\!0\.02\)/0\.660\.660\.60​\(−0\.05\)0\.60\\ \(\\\!\-\\\!0\.05\)/0\.820\.820\.57​\(−0\.03\)0\.57\\ \(\\\!\-\\\!0\.03\)0\.86​\(−0\.01\)0\.86\\ \(\\\!\-\\\!0\.01\)60​\(−4\)60\\ \(\\\!\-\\\!4\)Direction\-head Activations zeroed0\.70​\(−0\.18\)0\.70\\ \(\\\!\-\\\!0\.18\)/0\.970\.970\.32​\(−0\.13\)0\.32\\ \(\\\!\-\\\!0\.13\)/0\.850\.850\.68​\(−0\.14\)0\.68\\ \(\\\!\-\\\!0\.14\)/0\.990\.990\.15​\(−0\.07\)0\.15\\ \(\\\!\-\\\!0\.07\)/0\.660\.660\.48​\(−0\.17\)0\.48\\ \(\\\!\-\\\!0\.17\)/0\.820\.820\.47​\(−0\.13\)0\.47\\ \(\\\!\-\\\!0\.13\)0\.86​\(−0\.02\)0\.86\\ \(\\\!\-\\\!0\.02\)56​\(−8\)56\\ \(\\\!\-\\\!8\)

Table 14:Test\-time ablation results\. Numbers for the baseline row are taken from the Epoch\-320 training log; ablation rows are architectural predictions\.Δ\\Deltavalues \(in parentheses\) are relative to the baseline row; negativeΔ\\Deltaindicates the model relies on the ablated component\. Forcleanablations the architecture was trained to support the ablated state; forshiftablations theirΔ\\Deltavalues should be read as upper bounds on true necessity\.

## Appendix JRobustness to Missing Data

A practical causal discovery method must handle incomplete observations, which are ubiquitous in real\-world data \(e\.g\., missing medical records, failed sensors, detection limits\)\.FoundCausenatively incorporates missingness via a two\-channel input representation\[xn​i​mn​i,mn​i\]\[x\_\{ni\}\\,m\_\{ni\},\\,m\_\{ni\}\], wherexn​ix\_\{ni\}is the observed value of variableiiin samplennandmn​i∈\{0,1\}m\_\{ni\}\\in\\\{0,1\\\}is the observation mask\. The mask is propagated throughout the model, including attention layers, pooling, and normalization steps \(Sec\.[3](https://arxiv.org/html/2606.17516#S3)\)\. To evaluate robustness, we re\-run the full benchmark suite under*Missing At Random*\(MAR\) corruption at two levels and compare against the same baselines as in Table[1](https://arxiv.org/html/2606.17516#S5.T1)\.

##### MAR Generation\.

For each dataset with variables𝒱\\mathcal\{V\}, we randomly partition variables into a*driver*set𝒟\\mathcal\{D\}and a*target*set𝒯=𝒱∖𝒟\\mathcal\{T\}=\\mathcal\{V\}\\setminus\\mathcal\{D\}, with\|𝒟\|=\|𝒯\|\|\\mathcal\{D\}\|=\|\\mathcal\{T\}\|\. Driver variables are always observed, while target variables are subject to missingness\. For each target variablej∈𝒯j\\in\\mathcal\{T\}and samplenn, the observation mask is drawn as

Mn,j∼Bernoulli​\(1−σ​\(αj\+𝜷j⊤​𝑿n,𝒟std\)\),\\textbf\{M\}\_\{n,j\}\\sim\\mathrm\{Bernoulli\}\\\!\\left\(1\-\\sigma\\\!\\left\(\\alpha\_\{j\}\+\\bm\{\\beta\}\_\{j\}^\{\\top\}\\bm\{X\}^\{\\mathrm\{std\}\}\_\{n,\\mathcal\{D\}\}\\right\)\\right\),where𝑿n,𝒟std\\bm\{X\}^\{\\mathrm\{std\}\}\_\{n,\\mathcal\{D\}\}denotes standardized driver values \(zero mean, unit variance\),𝜷j∼𝒩​\(𝟎,𝑰\)\\bm\{\\beta\}\_\{j\}\\sim\\mathcal\{N\}\(\\bm\{0\},\\bm\{I\}\)is sampled independently for each target variable, andαj\\alpha\_\{j\}is chosen to match a desired marginal missing rate\. This construction satisfies the MAR assumption, as missingness depends only on observed variables\. The driver/target split is fixed per \(dataset, missing level\) pair for reproducibility\.

We evaluate two missingness levels:10%10\\%, which lies within the22–12%12\\%range seen during training, and30%30\\%, which represents a substantial out\-of\-distribution regime\.

##### Inference Protocol\.

FoundCausedirectly consumes masked inputs without imputation\. For methods requiring complete data \(PC, FCI, GES, GRaSP, DirectLiNGAM, ICA\-LiNGAM, SCORE, and correlation\-based baselines\), we apply per\-variable mean imputation\. For continuous optimization methods that support incomplete observations \(e\.g\., NOTEARS\-Linear, DAGMA\), we use pairwise\-deletion covariance estimation\. For amortized baselines \(AVICI, CSIvA, SEA, Cond\_FIP\), we provide masked inputs when supported and otherwise use mean imputation\. All other inference settings \(temperature scaling,K=10K=10bootstrap runs, GMM thresholding, and self\-consistency pruning\) are identical to the clean\-data evaluation\.

##### Results at10%10\\%MAR\.

Table[15](https://arxiv.org/html/2606.17516#A10.T15)reports performance under10%10\\%MAR missingness\.FoundCauseexhibits only a minor degradation \(F1:−0\.009F\_\{1\}:\-0\.009, AUROC:−0\.007\-0\.007\), consistent with this level being within the training distribution\. Amortized baselines that natively handle missingness \(e\.g\., AVICI\) show moderate relative drops of approximately22–3%3\\%inF1F\_\{1\}\. In contrast, classical methods relying on mean imputation degrade more substantially \(33–8%8\\%absolute inF1F\_\{1\}\), with PC and FCI additionally losing coverage on several datasets\.

##### Results at30%30\\%MAR\.

Table[16](https://arxiv.org/html/2606.17516#A10.T16)reports results at30%30\\%MAR, a regime well beyond the training distribution and challenging for imputation\-based approaches\. PC and FCI lose additional datasets due to rank\-deficient partial correlation estimates under high missingness, particularly on densed=41d=41PetShop graphs\. Classical methods using mean imputation degrade sharply \(0\.090\.09–0\.130\.13absolute inF1F\_\{1\}\), reflecting bias introduced under MAR\. In contrast,FoundCauseretains approximately92%92\\%of its clean\-dataF1F\_\{1\}\(0\.535→0\.4910\.535\\rightarrow 0\.491\), demonstrating graceful degradation\. This robustness suggests that the mask\-aware architecture effectively leverages missingness patterns, allowing information carried by observed values to be partially recovered through the mask channel\.

MethodCov\.AUROC↑\\uparrow\(Δ\\Delta\)F1F\_\{1\}↑\\uparrow\(Δ\\Delta\)Skel\-F1F\_\{1\}↑\\uparrow\(Δ\\Delta\)SHD/d/d↓\\downarrow\(Δ\\Delta\)*Classical baselines \(mean imputation\)*PC \(α=0\.01\\alpha\{=\}0\.01\)12/1512/150\.688​\(−0\.022\)0\.688\\,\(\-0\.022\)0\.435​\(−0\.029\)0\.435\\,\(\-0\.029\)0\.618​\(−0\.024\)0\.618\\,\(\-0\.024\)1\.471​\(\+0\.088\)1\.471\\,\(\+0\.088\)PC \(α=0\.05\\alpha\{=\}0\.05\)12/1512/150\.669​\(−0\.023\)0\.669\\,\(\-0\.023\)0\.401​\(−0\.029\)0\.401\\,\(\-0\.029\)0\.620​\(−0\.026\)0\.620\\,\(\-0\.026\)1\.586​\(\+0\.091\)1\.586\\,\(\+0\.091\)FCI \(α=0\.05\\alpha\{=\}0\.05\)12/1512/150\.614​\(−0\.023\)0\.614\\,\(\-0\.023\)0\.352​\(−0\.028\)0\.352\\,\(\-0\.028\)0\.459​\(−0\.024\)0\.459\\,\(\-0\.024\)1\.294​\(\+0\.085\)1\.294\\,\(\+0\.085\)GES15/1515/150\.702​\(−0\.023\)0\.702\\,\(\-0\.023\)0\.427​\(−0\.029\)0\.427\\,\(\-0\.029\)0\.607​\(−0\.024\)0\.607\\,\(\-0\.024\)1\.682​\(\+0\.086\)1\.682\\,\(\+0\.086\)GRaSP15/1515/150\.720​\(−0\.027\)0\.720\\,\(\-0\.027\)0\.458​\(−0\.030\)0\.458\\,\(\-0\.030\)0\.651​\(−0\.025\)0\.651\\,\(\-0\.025\)1\.541​\(\+0\.087\)1\.541\\,\(\+0\.087\)DirectLiNGAM15/1515/150\.641​\(−0\.034\)0\.641\\,\(\-0\.034\)0\.280​\(−0\.034\)0\.280\\,\(\-0\.034\)0\.478​\(−0\.032\)0\.478\\,\(\-0\.032\)2\.814​\(\+0\.133\)2\.814\\,\(\+0\.133\)ICA\-LiNGAM15/1515/150\.617​\(−0\.036\)0\.617\\,\(\-0\.036\)0\.258​\(−0\.034\)0\.258\\,\(\-0\.034\)0\.466​\(−0\.034\)0\.466\\,\(\-0\.034\)3\.159​\(\+0\.151\)3\.159\\,\(\+0\.151\)NOTEARS\-Linear15/1515/150\.570​\(−0\.018\)0\.570\\,\(\-0\.018\)0\.233​\(−0\.018\)0\.233\\,\(\-0\.018\)0\.493​\(−0\.016\)0\.493\\,\(\-0\.016\)1\.523​\(\+0\.053\)1\.523\\,\(\+0\.053\)DAGMA \(linear\)15/1515/150\.580​\(−0\.017\)0\.580\\,\(\-0\.017\)0\.250​\(−0\.017\)0\.250\\,\(\-0\.017\)0\.523​\(−0\.017\)0\.523\\,\(\-0\.017\)1\.597​\(\+0\.051\)1\.597\\,\(\+0\.051\)SCORE15/1515/150\.521​\(−0\.033\)0\.521\\,\(\-0\.033\)0\.112​\(−0\.031\)0\.112\\,\(\-0\.031\)0\.331​\(−0\.035\)0\.331\\,\(\-0\.035\)4\.790​\(\+0\.252\)4\.790\\,\(\+0\.252\)Naive correlation15/1515/150\.559​\(−0\.025\)0\.559\\,\(\-0\.025\)0\.164​\(−0\.026\)0\.164\\,\(\-0\.026\)0\.401​\(−0\.026\)0\.401\\,\(\-0\.026\)4\.832​\(\+0\.179\)4\.832\\,\(\+0\.179\)*Amortised / foundation\-model methods*AVICI15/1515/150\.719​\(−0\.013\)0\.719\\,\(\-0\.013\)0\.479​\(−0\.019\)0\.479\\,\(\-0\.019\)0\.636​\(−0\.015\)0\.636\\,\(\-0\.015\)1\.194​\(\+0\.052\)1\.194\\,\(\+0\.052\)CSIvA15/1515/150\.717​\(−0\.022\)0\.717\\,\(\-0\.022\)0\.481​\(−0\.026\)0\.481\\,\(\-0\.026\)0\.638​\(−0\.020\)0\.638\\,\(\-0\.020\)1\.175​\(\+0\.067\)1\.175\\,\(\+0\.067\)SEA15/1515/150\.725​\(−0\.019\)0\.725\\,\(\-0\.019\)0\.490​\(−0\.022\)0\.490\\,\(\-0\.022\)0\.650​\(−0\.017\)0\.650\\,\(\-0\.017\)1\.113​\(\+0\.061\)1\.113\\,\(\+0\.061\)Cond\_FIP15/1515/150\.706​\(−0\.022\)0\.706\\,\(\-0\.022\)0\.464​\(−0\.025\)0\.464\\,\(\-0\.025\)0\.627​\(−0\.019\)0\.627\\,\(\-0\.019\)1\.246​\(\+0\.070\)1\.246\\,\(\+0\.070\)FoundCause\(ours\)15/1515/150\.749​\(−0\.007\)\\bm\{0\.749\}\\,\(\\bm\{\-0\.007\}\)0\.526​\(−0\.009\)\\bm\{0\.526\}\\,\(\\bm\{\-0\.009\}\)0\.656​\(−0\.007\)\\bm\{0\.656\}\\,\(\\bm\{\-0\.007\}\)1\.013​\(\+0\.033\)\\bm\{1\.013\}\\,\(\\bm\{\+0\.033\}\)Table 15:Results under𝟏𝟎%\\bm\{10\\%\}MAR missingness\. Parenthesized values denote changes relative to the clean\-data baseline \(Table[1](https://arxiv.org/html/2606.17516#S5.T1)\)\. Negative deltas for AUROC,F1F\_\{1\}, and Skel\-F1F\_\{1\}indicate performance degradation, while positive deltas forSHD/d\\mathrm\{SHD\}/dindicate degradation\.FoundCauseexhibits the smallest drop across all metrics\.MethodCov\.AUROC↑\\uparrow\(Δ\\Delta\)F1F\_\{1\}↑\\uparrow\(Δ\\Delta\)Skel\-F1F\_\{1\}↑\\uparrow\(Δ\\Delta\)SHD/d/d↓\\downarrow\(Δ\\Delta\)*Classical baselines \(mean imputation\)*PC \(α=0\.01\\alpha\{=\}0\.01\)10/1510/150\.624​\(−0\.086\)0\.624\\,\(\-0\.086\)0\.362​\(−0\.102\)0\.362\\,\(\-0\.102\)0\.557​\(−0\.085\)0\.557\\,\(\-0\.085\)1\.704​\(\+0\.321\)1\.704\\,\(\+0\.321\)PC \(α=0\.05\\alpha\{=\}0\.05\)10/1510/150\.607​\(−0\.085\)0\.607\\,\(\-0\.085\)0\.330​\(−0\.100\)0\.330\\,\(\-0\.100\)0\.561​\(−0\.085\)0\.561\\,\(\-0\.085\)1\.834​\(\+0\.339\)1\.834\\,\(\+0\.339\)FCI \(α=0\.05\\alpha\{=\}0\.05\)10/1510/150\.550​\(−0\.087\)0\.550\\,\(\-0\.087\)0\.286​\(−0\.094\)0\.286\\,\(\-0\.094\)0\.402​\(−0\.081\)0\.402\\,\(\-0\.081\)1\.498​\(\+0\.289\)1\.498\\,\(\+0\.289\)GES15/1515/150\.629​\(−0\.096\)0\.629\\,\(\-0\.096\)0\.344​\(−0\.112\)0\.344\\,\(\-0\.112\)0\.543​\(−0\.088\)0\.543\\,\(\-0\.088\)1\.972​\(\+0\.376\)1\.972\\,\(\+0\.376\)GRaSP15/1515/150\.641​\(−0\.106\)0\.641\\,\(\-0\.106\)0\.370​\(−0\.118\)0\.370\\,\(\-0\.118\)0\.581​\(−0\.095\)0\.581\\,\(\-0\.095\)1\.798​\(\+0\.344\)1\.798\\,\(\+0\.344\)DirectLiNGAM15/1515/150\.558​\(−0\.117\)0\.558\\,\(\-0\.117\)0\.193​\(−0\.121\)0\.193\\,\(\-0\.121\)0\.401​\(−0\.109\)0\.401\\,\(\-0\.109\)3\.215​\(\+0\.534\)3\.215\\,\(\+0\.534\)ICA\-LiNGAM15/1515/150\.536​\(−0\.117\)0\.536\\,\(\-0\.117\)0\.174​\(−0\.118\)0\.174\\,\(\-0\.118\)0\.391​\(−0\.109\)0\.391\\,\(\-0\.109\)3\.528​\(\+0\.520\)3\.528\\,\(\+0\.520\)NOTEARS\-Linear15/1515/150\.491​\(−0\.097\)0\.491\\,\(\-0\.097\)0\.160​\(−0\.091\)0\.160\\,\(\-0\.091\)0\.417​\(−0\.092\)0\.417\\,\(\-0\.092\)1\.716​\(\+0\.246\)1\.716\\,\(\+0\.246\)DAGMA \(linear\)15/1515/150\.503​\(−0\.094\)0\.503\\,\(\-0\.094\)0\.176​\(−0\.091\)0\.176\\,\(\-0\.091\)0\.447​\(−0\.093\)0\.447\\,\(\-0\.093\)1\.798​\(\+0\.252\)1\.798\\,\(\+0\.252\)SCORE15/1515/150\.424​\(−0\.130\)0\.424\\,\(\-0\.130\)0\.051​\(−0\.092\)0\.051\\,\(\-0\.092\)0\.237​\(−0\.129\)0\.237\\,\(\-0\.129\)5\.341​\(\+0\.803\)5\.341\\,\(\+0\.803\)Naive correlation15/1515/150\.481​\(−0\.103\)0\.481\\,\(\-0\.103\)0\.100​\(−0\.090\)0\.100\\,\(\-0\.090\)0\.333​\(−0\.094\)0\.333\\,\(\-0\.094\)5\.427​\(\+0\.774\)5\.427\\,\(\+0\.774\)*Amortised / foundation\-model methods*AVICI15/1515/150\.685​\(−0\.047\)0\.685\\,\(\-0\.047\)0\.432​\(−0\.066\)0\.432\\,\(\-0\.066\)0\.596​\(−0\.055\)0\.596\\,\(\-0\.055\)1\.318​\(\+0\.176\)1\.318\\,\(\+0\.176\)CSIvA15/1515/150\.672​\(−0\.067\)0\.672\\,\(\-0\.067\)0\.418​\(−0\.089\)0\.418\\,\(\-0\.089\)0\.586​\(−0\.072\)0\.586\\,\(\-0\.072\)1\.334​\(\+0\.226\)1\.334\\,\(\+0\.226\)SEA15/1515/150\.681​\(−0\.063\)0\.681\\,\(\-0\.063\)0\.428​\(−0\.084\)0\.428\\,\(\-0\.084\)0\.603​\(−0\.064\)0\.603\\,\(\-0\.064\)1\.286​\(\+0\.234\)1\.286\\,\(\+0\.234\)Cond\_FIP15/1515/150\.649​\(−0\.079\)0\.649\\,\(\-0\.079\)0\.388​\(−0\.101\)0\.388\\,\(\-0\.101\)0\.570​\(−0\.076\)0\.570\\,\(\-0\.076\)1\.451​\(\+0\.275\)1\.451\\,\(\+0\.275\)FoundCause\(ours\)15/1515/150\.717​\(−0\.039\)\\bm\{0\.717\}\\,\(\\bm\{\-0\.039\}\)0\.491​\(−0\.044\)\\bm\{0\.491\}\\,\(\\bm\{\-0\.044\}\)0\.631​\(−0\.032\)\\bm\{0\.631\}\\,\(\\bm\{\-0\.032\}\)1\.132​\(\+0\.152\)\\bm\{1\.132\}\\,\(\\bm\{\+0\.152\}\)Table 16:Results under𝟑𝟎%\\bm\{30\\%\}MAR missingness\. This setting represents heavy corruption, exceeding by3×3\\timesthe maximum missingness seen duringFoundCausetraining\. Parenthesized values denote changes relative to the clean\-data baseline \(Table[1](https://arxiv.org/html/2606.17516#S5.T1)\)\. Methods requiring complete data \(e\.g\., PC, FCI\) lose additional coverage, and classical imputation\-based methods degrade substantially\. In contrast,FoundCauseexhibits a modest relative drop of8\.2%8\.2\\%inF1F\_\{1\}\(0\.535→0\.4910\.535\\rightarrow 0\.491\)\.
##### Discussion\.

Three observations emerge from Tables[15](https://arxiv.org/html/2606.17516#A10.T15)and[16](https://arxiv.org/html/2606.17516#A10.T16)\. First,FoundCause’s native handling of missingness provides a consistent advantage across all metrics and both missingness levels: itsF1F\_\{1\}drop is−0\.009\-0\.009at10%10\\%MAR and−0\.044\-0\.044at30%30\\%MAR, roughly half that of the next best method \(AVICI\), which also supports mask\-aware inputs\.

Second, classical methods relying on mean imputation exhibit increasing bias as missingness grows\. For example, the strongest classical baseline \(GRaSP\) degrades fromF1=0\.488→0\.458→0\.370F\_\{1\}=0\.488\\rightarrow 0\.458\\rightarrow 0\.370across clean,10%10\\%, and30%30\\%MAR settings, corresponding to a24%24\\%relative decline\.

Third,FoundCausemaintains top rank on all metrics while achieving full coverage \(15/1515/15datasets\) at both missingness levels, making it the only method that remains consistently usable across a wide range of data\-quality conditions\. Notably, the30%30\\%MAR setting lies well beyond the12%12\\%maximum missingness seen during training, indicating that the model’s mask\-aware design generalizes beyond the training distribution—a property that imputation\-based methods do not possess by construction\.

## Appendix KPerformance of SEA and CSIvA on \>50 Dimensional Datasets

For comparison, we report results of SEA and CSIvA on the same four datasets used to evaluate the dimension generalization ofFoundCause\.

DatasetDDAUROC↑\\uparrowF1↑F\_\{1\}\\uparrowSkel\-F1↑F\_\{1\}\\uparrowSHD/d↓/d\\downarrowHailfinder\[[1](https://arxiv.org/html/2606.17516#bib.bib30)\]56560\.7210\.7210\.4940\.4940\.5890\.5891\.281\.28HEPAR II\[[39](https://arxiv.org/html/2606.17516#bib.bib25)\]70700\.7280\.7280\.4380\.4380\.5960\.5961\.511\.51WIN95PT\[[46](https://arxiv.org/html/2606.17516#bib.bib29)\]76760\.7710\.7710\.4820\.4820\.5630\.5631\.631\.63DREAM4\-100\[[30](https://arxiv.org/html/2606.17516#bib.bib24),[29](https://arxiv.org/html/2606.17516#bib.bib32)\]1001000\.7020\.7020\.4060\.4060\.5680\.5681\.741\.74Table 17:Dimension\-generalization of SEA on additional real\-world datasets\.DatasetDDAUROC↑\\uparrowF1↑F\_\{1\}\\uparrowSkel\-F1↑F\_\{1\}\\uparrowSHD/d↓/d\\downarrowHailfinder\[[1](https://arxiv.org/html/2606.17516#bib.bib30)\]56560\.7580\.7580\.5330\.5330\.6150\.6151\.191\.19HEPAR II\[[39](https://arxiv.org/html/2606.17516#bib.bib25)\]70700\.7260\.7260\.4390\.4390\.5940\.5941\.551\.55WIN95PT\[[46](https://arxiv.org/html/2606.17516#bib.bib29)\]76760\.7960\.7960\.4910\.4910\.5710\.5711\.651\.65DREAM4\-100\[[30](https://arxiv.org/html/2606.17516#bib.bib24),[29](https://arxiv.org/html/2606.17516#bib.bib32)\]1001000\.5890\.5890\.2760\.2760\.4620\.4622\.412\.41Table 18:Dimension\-generalization of CSIvA on additional real\-world datasets\.Tables[17](https://arxiv.org/html/2606.17516#A11.T17)and[18](https://arxiv.org/html/2606.17516#A11.T18)present results on four Bayesian networks withD∈\{56,70,76,100\}D\\in\\\{56,70,76,100\\\}, all exceedingFoundCause’s training range \(D≤50D\\leq 50\)\. SEA is trained up toD=100D=100, while CSIvA is trained up toD=80D=80\. As expected, each method performs best on datasets closest to its training distribution: CSIvA leads on discrete benchmarks within its range \(e\.g\.,hailfinder,win95pts\), while SEA performs best ondream4\-100at its training ceiling\.

From Table[3](https://arxiv.org/html/2606.17516#S5.T3),FoundCause, despite extrapolating1\.51\.5–22x beyond its training dimension, remains competitive, within0\.020\.02–0\.050\.05inF1F\_\{1\}of the best method on three of the four datasets, and achieves the highest AUROC onwin95pts\. The largest gap occurs ondream4\-100, reflecting the difficulty of22x extrapolation; closing this gap would require extending the training range, which the architecture supports without modification\.