Demystifying amortized causal discovery with transformers

## Abstract

Abstract Supervised learning approaches for causal discovery from observational data often achieve competitive performance despite seemingly avoiding explicit assumptions that traditional methods make for identifiability. In this work, we investigate CSIvA [ 1 ] , a transformer-based model promising to train on synthetic data and transfer to real data. First, we bridge the gap with existing identifiability theory and show that constraints on the training data distribution implicitly define a prior on the test observations. Consistent with classical approaches, good performance is achieved when we have a good prior on the test data, and the underlying model is identifiable. At the same time, we find new trade-offs. Training on datasets generated from different classes of causal models, unambiguously identifiable in isolation, improves the test generalization. Performance is still guaranteed, as the ambiguous cases resulting from the mixture of identifiable causal models are unlikely to occur (which we formally prove).
Overall, our study finds that amortized causal discovery still needs to obey identifiability theory, but it also differs from classical methods in how the assumptions are formulated, trading more reliance on assumptions on the noise type for fewer hypotheses on the mechanisms.

## 1 Introduction

Causal discovery aims to uncover the underlying causal relationships between variables of a system from pure observations, which is crucial for answering interventional and counterfactual queries when experimentation is impractical or unfeasible [, , ]. Unfortunately, causal discovery is inherently ill-posed []: unique identification of causal directions requires restrictive assumptions on the class of structural causal models (SCMs) that generated the data [, , ]. These theoretical limitations often render existing methods inapplicable, as the underlying assumptions are usually untestable or difficult to verify in practice [].

Recently, supervised learning algorithms trained on synthetic data have been proposed to overcome the need for specific hypotheses, which restrains the application of classical causal discovery methods to real-world problems [, , , , ]. Seminal work from Lopez-Paz et al. [] argues that this learning-based approach to causal discovery would allow dealing with complex
data-generating processes and would greatly reduce the
need for explicitly crafting identifiability conditions a-priori: despite this ambitious goal, the output of these methods is generally considered unreliable, as no theoretical guarantee is provided. A pair of non-identifiable structural causal models can be associated with different causal graphs $\mathcal{G}\neq\tilde{\mathcal{G}}$, while entailing the same joint distribution $p$ on the system’s variables. It is thus unclear how a learning algorithm presented with observational data generated from $p$ would be able to overcome these theoretical limits and correctly identify a unique causal structure. However, the available empirical evidence seems not to care about impossibility results, as these methods yield surprising generalization results on several synthetic benchmarks. Our work aims to bridge this gap by studying the performance of a transformer architecture for causal discovery through the lens of the theory of identifiability from observational data.
Specifically, we analyze the CSIvA (Causal Structure Induction via Attention) model for causal discovery [], focusing on bivariate graphs, as they offer a controlled yet non-trivial setting for the investigation.
As our starting point, we provide closed-form examples that identify the limitations of CSIvA in recovering causal structures of linear non-Gaussian and nonlinear additive noise models, which are notably identifiable, and demonstrate the expected failures through empirical evidence.
These findings suggest that the class of structural causal models that can be identified by CSIvA is inherently dependent on the specific class of SCMs observed during training. Thus, the need for restrictive hypotheses on the data-generating process is intrinsic to causal discovery, both in the traditional and modern learning-based approaches: assumptions on the test distribution either are posited when selecting the algorithm (traditional methods) or in the choice of the training data (learning-based methods). To address this limitation, we theoretically and empirically analyze *when* training CSIvA on datasets generated by multiple identifiable SCMs with different structural assumptions improves its generalization at test time.
In summary:

- •
We show that the class of structural causal models that CSIvA can identify is defined by the class of SCMs observed through samples during the training. We reinforce the notion that identifiability in causal discovery inherently requires assumptions, which must be encoded in the training data in the case of learning-based approaches.
- •
To overcome this limitation, we study the benefits of CSIvA training on mixtures of causal models. We analyze when algorithms learned on multiple models are expected to identify broad classes of SCMs (unlike many classical methods). Empirically, we show that training on samples generated by multiple identifiable causal models with different assumptions on mechanisms and noise distribution results in significantly improved generalization abilities.

#### Closely related works and their relation with CSIvA.

In this paper, we study amortized inference of causal graphs, i.e. optimization of an inference model to directly predict a causal structure from newly provided data. This is the first work that attempts to understand the connection between identifiability theory and amortized inference, while several algorithms have been proposed. In the context of purely observational data, Lopez-Paz et al. [] defines a distribution regression problem [] mapping the kernel mean embedding of the data distribution to a causal graph, while Li et al. [] relies on equivariant neural network architectures. More recently, Lippe et al. [] and Lorch et al. [] proposed learning on interventional data, in addition to observations (in the same spirit as CSIvA). Despite different algorithmic implementations, the target object of estimation of most of these methods is the distribution over the space of all possible graphs, conditional on the input dataset (similarly, the ENCO algorithm in Lippe et al. [] models the conditional distribution of individual edges). This justifies our choice of restricting our study to the CSIvA architecture (despite this being a clear limitation), as in the infinite observational sample limit, these methods approximate the same distribution. Methods necessarily requiring interventional data [, , ], and learning-based algorithms unsuitable for amortized inference [, , , , ] are out of the scope of this work.

## 2 Background and motivation

We start introducing structural causal models (SCMs), an intuitive framework that formalizes causal relations. Let $X$ be a set of random variables in $\mathbb{R}$ defined according to the set of structural equations:

$$ $X_{i}\coloneqq f_{i}(X_{\operatorname{PA}^{\mathcal{G}}_{i}},N_{i}),\hskip 5.6 9054pt\forall i=1,\ldots,k.$ (1) $$

$N_{i}\in\mathbb{R}$ are noise random variables. The function $f_{i}$ is the causal mechanism mapping the set of direct causes $X_{\operatorname{PA}^{\mathcal{G}}_{i}}$ of $X_{i}$ and the noise term $N_{i}$, to $X_{i}$’s value. The causal graph $\mathcal{G}$ is a directed acyclic graph (DAG) with nodes $X=\{X_{1},\ldots,X_{k}\}$, and edges $\{X_{j}\rightarrow X_{i}:X_{j}\in X_{\operatorname{PA}^{\mathcal{G}}_{i}}\}$, with $\operatorname{PA}^{\mathcal{G}}_{i}$ indices of the parent nodes of $X_{i}$ in $\mathcal{G}$. The causal model induces a density $p_{X}$ over the vector $X$.

### 2.1 Causal discovery from observational data

Causal discovery from observational data is the inference of the causal graph $\mathcal{G}$ from a dataset of i.i.d. observations of the random vector $X$. In general, without restrictive assumptions on the mechanisms and the noise distributions, the direction of edges in the graph $\mathcal{G}$ is not identifiable, i.e. it can not be found from the population density $p_{X}$. In particular, it is possible to identify only a Markov equivalence class, which is the set of graphs encoding the same conditional independencies as the density $p_{X}$.
To clarify with an example, consider the causal graph $X_{1}\rightarrow X_{2}$ associated with a structural causal model inducing a density $p_{X_{1},X_{2}}$. If the model is not identifiable, there exists an SCM with causal graph $X_{2}\to X_{1}$ that entails the same joint density $p_{X_{1},X_{2}}$. The set $\{X_{1}\rightarrow X_{2},X_{2}\rightarrow X_{1}\}$ is the Markov equivalence class of the graph $X_{1}\rightarrow X_{2}$, i.e. the set of all graphs with $X_{1},X_{2}$ mutually dependent. Clearly, in this setting, even the exact knowledge of $p_{X_{1},X_{2}}$ cannot inform us about the correct causal direction.

###### Definition 1 (Identifiable causal model) .

Consider a structural causal model with underlying graph $\mathcal{G}$ and $p_{X}$ joint density of the causal variables. We say that the model is identifiable from observational data if the density $p_{X}$ can not be entailed by a structural causal model with graph $\tilde{\mathcal{G}}\neq\mathcal{G}$.

We define the post-additive noise model (post-ANM) as the causal model with the set of equations:

$$ $X_{i}\coloneqq f_{2,i}(f_{1,i}(X_{\operatorname{PA}^{\mathcal{G}}_{i}})+N_{i}) ,\hskip 2.84526pt\forall i=1,\ldots,d,$ (2) $$

with $f_{2,i}$ invertible map and mutually independent noise terms. When $f_{2,i}$ is a nonlinear function, the post-ANM amounts to the identifiable post-nonlinear model (PNL) []. When $f_{2,i}$ is the identity function and $f_{1,i}$ nonlinear, it simplifies to the nonlinear additive noise model (ANM)[, ], which is known to be identifiable, and is described by the set of structural equations:

$$ $X_{i}\coloneqq f_{1,i}(X_{\operatorname{PA}^{\mathcal{G}}_{i}})+N_{i}.$ (3) $$

If, additionally, we restrict the mechanisms $f_{1,i}$ to be linear and the noise terms $N_{i}$ to a non-Gaussian distribution, we recover the identifiable linear non-Gaussian additive model or LiNGAM []:

$$ $X_{i}=\sum_{j\in\operatorname{PA}_{i}^{\mathcal{G}}}\alpha_{j}X_{j}+N_{i}, \hskip 10.00002pt\alpha_{j}\in\mathbb{R}.$ (4) $$

###### Definition 1 (Identifiable causal model) .

### 2.2 Motivation and problem definition

Causal discovery from observational data relies on specific assumptions, which can be challenging to verify in practice []. To address this, recent methods leverage supervised learning for the amortized inference of causal graphs [, , , , , , ], optimizing an inference model to directly predict a causal structure from a provided dataset. While these approaches aim to reduce reliance on explicit identifiability assumptions, they often lack a clear connection to the existing causal discovery theory, making their outputs generally unreliable. We illustrate this limitation through an example.

###### Example 1 .

We consider the CSIvA transformer architecture proposed by Ke et al. [], which can learn a map from observational data to a causal graph. The authors of the paper show that, in the infinite sample regime, the CSIvA architecture exactly approximates the conditional distribution $p(\cdot|\mathcal{D})$ over the space of possible graphs, given a dataset $\mathcal{D}$. Identifiability theory in causal discovery tells us that if the class of structural causal models that generated the observations is sufficiently constrained, then there is only one graph that can fit the data within that class. For example, consider the case of a dataset that is known to be generated by a nonlinear additive noise model, and let $p(\cdot|\mathcal{D},\textnormal{ANM})$ be the conditional distribution that incorporates this prior knowledge on the SCM: then $p(\cdot|\mathcal{D},\textnormal{ANM})$ concentrates all the mass on a single point $\mathcal{G}^{*}$, the true graph underlying the $\mathcal{D}$ observations.
Instead, in the absence of restrictions on the structural causal model, all the graphs in a Markov equivalence class are equally likely to be the correct solution given the data. Hence, $p(\cdot|\mathcal{D})$, the distribution learned by CSIvA, assigns equal probability to each graph in the Markov equivalence class of $\mathcal{G}^{*}$.

Our arguments of Example 1 are valid for all learning methods that approximate the conditional distribution over the space of graphs given the input data [, , , , ], and suggest that these algorithms are at most informative about the equivalence class of the causal graph underlying the observations. However, the available empirical evidence does not seem to highlight these limitations, as in practice these methods can infer the true causal DAG on several synthetic benchmarks. Thus, further investigation is necessary if we want to rely on their output in any meaningful sense.
In this work, we analyze these "black-box" approaches through the lens of established theory of causal discovery from observational data (causal inference often lacks experimental data, which we do not consider).
We study in detail the CSIvA architecture [] (see Appendix A), a variation of the transformer neural network [] for the supervised learning of algorithms for amortized causal discovery. This model is optimized via maximum likelihood estimation, i.e. finding $\Theta$ that minimizes $-\mathbf{E}_{\mathcal{G},\mathcal{D}}[\ln\hat{p}(\mathcal{G}|\mathcal{D};
\Theta)]$, where $\hat{p}(\mathcal{G}|\mathcal{D};\Theta)$ is the conditional distribution of a graph $\mathcal{G}$ given a dataset $\mathcal{D}$ parametrized by $\Theta$.
We limit the analysis to CSIvA as it is a simple yet competitive end-to-end approach to learning causal models.
While this is clearly a limitation of the paper, our theoretical and empirical conclusions exemplify both the role of theoretical identifiability in modern approaches and the new opportunities they provide. Additionally, it fits well within a line of works arguing that specifically transformers can learn causal concepts [, , ] and identify different assumptions in context [].

###### Example 1 .

## 3 Experimental results through the lens of theory

In this section, we present a comprehensive analysis of causal discovery with transformers and its relation to the theoretical boundaries of causal discovery from observational data. We show that suitable assumptions must be encoded in the training distribution to ensure the identifiability of the test data, and we additionally study the effectiveness of training on mixtures of causal models to overcome these limitations, improving generalization abilities.

### 3.1 Experimental design

We concentrate our research on causal models of two variables, causally related according to one of the two graphs $X\rightarrow Y$, $Y\rightarrow X$.
Bivariate models are the simplest non-trivial setting with a well-known theory of causality inference [, , ], but also amenable to manipulation. This allows for comprehensive training and analysis of diverse SCMs and facilitates a clear interpretation of the results.

#### Datasets.

Unless otherwise specified, in our experiments we train CSIvA on a sample of $15000$ synthetically generated datasets, consisting of $1500$ i.i.d. observations. Each dataset is generated according to a single class of SCMs, defined by the mechanism type and the noise terms distribution. The coefficients of the linear mechanisms are sampled in the range $[-3,-0.5]\cup[0.5,3]$, removing small coefficients to avoid close-to-unfaithful effects []. Nonlinear mechanisms are parametrized according to a neural network with random weights, a strategy commonly adopted in the literature of causal discovery [, ]. The post-nonlinearity of the PNL model consists of a simple map $z\mapsto z^{3}$. Noise terms are sampled from common distributions and a randomly generated density that we call mlp, previously adopted in Montagna et al. [], defined by a standard Gaussian transformed by a multilayer perceptron (MLP) (Appendix B.2). We name these datasets mechanism-noise to refer to their underlying causal model. For example, data sampled from a nonlinear ANM with Gaussian noise are named nonlinear-gaussian. More details on the synthetic data generation schema are found in Appendix B.2. All data are standardized by their empirical variance to remove opportunities to learn shortcuts [, , ].

#### Metric and random baseline.

As our metric we use the structural Hamming distance (SHD), which is the number of edge removals, insertions or flips to transform one graph to another. In the context of bivariate causal graphs with a single edge, this is simply an error count, so correct inference corresponds to $\textnormal{SHD}=0$, and an incorrect prediction gives $\textnormal{SHD}=1$. Additionally, we define a reference random baseline, which assigns a causal direction according to a fair coin, achieving $\textnormal{SHD}=0.5$ in expectation. Each architecture we analyze in the experiments is trained $3$ times, with different parameter initialization and training samples: the SHD presented in the plots is the average of each of the $3$ models on $1500$ distinct test datasets of $1500$ points each, and the error bars are $95\%$ confidence intervals.

We detail the training hyperparameters in Appendix B.1. Next, we analyze our experimental results, starting by investigating how well CSIvA generalizes on distributions unseen during training.

### 3.2 Warm up: is CSIvA capable of in and out-of-distribution generalization?

#### In-distribution generalization.

First, we investigate the generalization of CSIvA on datasets sampled from the structural casual model that generates the train distribution, with mechanisms and noise distributions fixed between training and testing. We call this in-distribution generalization. As a benchmark, we present the performance of several state-of-the-art approaches from the literature on causal discovery: we consider the DirectLiNGAM, and NoGAM algorithms [, ], respectively designed for the inference on LiNGAM and nonlinear ANM generated data(^1^11The [causal-learn](https://causal-learn.readthedocs.io/en/latest/) implementation of the PNL algorithm could not perform better than random on our synthetic post-nonlinear data, and we observed that this was due to the sensitivity of the algorithm to the variance scale. So we report the plot of Figure 1(c) without benchmark comparison. We remark that the point of this experiment is not to make any claims on CSIvA being state-of-the-art but to validate that the performance we obtain in our re-implementation is non-trivial. This is clear for PNL, even without comparison.). The results of Figure 1 show that CSIvA can properly generalize to unseen samples from the training distribution: the majority of the trained models present SHD close to zero and comparable to the relative benchmark algorithm.

Figure: (a) Linear
Refer to caption: x1.png

#### Out-of-distribution generalization.

In practice, we generally do not know the SCM defining the test distribution, so we are interested in CSIvA’s ability to generalize to data sampled from a class of causal models that is unobserved during training. We call this out-of-distribution generalization (OOD). We study OOD generalization to different noise terms, analyzing the network performance on datasets generated from causal models where the mechanisms are fixed with respect to the training, while the noise distribution varies (e.g., given linear-mlp training samples, testing occurs on linear-uniform data). Orthogonally to these experiments, we empirically validate CSIvA’s OOD generalization over different mechanism types (linear, nonlinear, post-nonlinear), while leaving the noise distribution (mlp) fixed across test and training. In Figure 2(a), we observe that CSIvA cannot generalize across the different mechanisms, as the SHD of a network tested on unseen causal mechanisms approximates that of the random baseline. Further, Figure 2(b) shows that out-of-distribution generalization across noise terms does not work reliably, and it is hard to predict when it might occur.

#### Implications.

CSIvA generalizes well to test data generated by the same class of SCMs used for training, in line with the findings in Ke et al. [], which validates our implementation and training procedure. However, it struggles when the test data are out-of-distribution, not generated by causal models with the same mechanisms and noise terms it was trained on. While training on a wider class of SCMs might overcome this limitation, it requires caution. The identifiability of causal graphs indeed results from the interplay between the data-generating mechanisms and noise distribution. However, as we argue in our Example 1, the class of causal models that a supervised learning algorithm can identify is generally not clear.
In what follows, we investigate this point and its implications for CSIvA, showing that the identifiability of the test samples can be ensured by imposing suitable assumptions on the class of SCMs generating the training distribution.

Figure: (a) Mechanisms
Refer to caption: x5.png

### 3.3 How does CSIvA relate to identifiability theory for causal graphs?

The CSIvA algorithm does not make structural assumptions about the causal model underlying the input data. This implies that the output of this method is unclear: as CSIvA targets the conditional distribution $p(\cdot|\mathcal{D})$ over the space of graphs, in the absence of restrictions on the functional mechanisms and the distribution of the noise terms, the causal graph $X\to Y$ is indistinguishable from $Y\to X$, as they are both equally likely to underlie the joint density $p_{X,Y}$ generating the data.
As we discuss in Example 1, the graphical output of the trained architecture could at most identify the equivalence class of the true causal graph. Yet, our experiments of Section 3.2 show that CSIvA is capable of good in-distribution generalization, often inferring the correct DAG at test time.
We explain this seeming contradiction with the following hypothesis, which motivates the analysis in the remainder of this section.

###### Hypothesis 1 .

The class of structural causal models that can be identified by CSIvA is defined by the class of structural causal models underlying the generation of the training data.

To support and clarify our statement, we present the following example, adapted from Hoyer et al. [].

###### Example 2 .

Consider the causal model $Y=f(X)+N,$ where $f(X)=-X$ and $p_{X},p_{N}$ are Gumbel densities $p_{X}(x)=\operatorname{exp}(-x-\operatorname{exp}(-x))$ and $p_{N}(n)=\operatorname{exp}(-n-\operatorname{exp}(-n))$. This model satisfies the assumptions of the LiNGAM, so it is identifiable, in the sense that a backward linear model with the same distribution does not exist. However, in this special case, we can build a backward nonlinear additive noise model $X=g(Y)+\tilde{N}$ with independent noise terms: taking $p_{Y}(y)=\operatorname{exp}(-y-2\log(1+\operatorname{exp}(-y)))$ to be the density of a logistic distribution, $p_{\tilde{N}}(\tilde{n})=\operatorname{exp}(-2\tilde{n}-\operatorname{exp}(-
\tilde{n}))$ and $g(y)=\log(1+\operatorname{exp}(-y))$; we see that $p_{X,Y}$ can factorize according to two opposite causal directions, as $p_{X,Y}(x,y)=p_{N}(y-f(x))p_{X}(x)=p_{\tilde{N}}(x-g(y))p_{Y}(y)$.
Given a dataset $\mathcal{D}$ of observations from the forward linear model, causal discovery methods like DirectLiNGAM [] can provably identify the correct causal direction $X\rightarrow Y$, assuming that sufficient samples are provided.
Instead, the behavior of CSIvA seems hard to predict: given that the network approximates the conditional distribution $p(\cdot|\mathcal{D})$ over the possible graphs, for $\mathcal{D}$ with arbitrary many samples we have $p(X\rightarrow Y|\mathcal{D})=p(Y\rightarrow X|\mathcal{D})=0.5$. On the other hand, given the prior knowledge that the data-generating SCM is a linear non-gaussian additive noise model, we have $p(X\rightarrow Y|\mathcal{D},\textnormal{LiNGAM})=1$, because the LiNGAM is identifiable. In this sense, the class of structural causal models that CSIvA correctly infers appears to be determined by the structural causal models underlying the generation of the training data. Under our Hypothesis 1, training CSIvA exclusively on LiNGAM-generated data is equivalent to learning the distribution $p(\cdot|\mathcal{D},\textnormal{LiNGAM})$, such that the network should be able to identify the forward linear model, whereas it could only infer the equivalence class of the causal graph if its training datasets include observations from a nonlinear additive noise model.

The empirical results of Figure 3(a) show that CSIvA behaves according to our hypothesis: when training exclusively occurs on datasets $\{\mathcal{D}_{i,\to}\}_{i}$ generated by the forward linear-gumbel model of Example 2, the network can identify the causal direction of test data generated according to the same SCM. Similarly, the transformer trained on datasets
$\{\mathcal{D}_{i,\leftarrow}\}_{i}$ from the backward nonlinear model of the example can generalize to test data coming from the same distribution. According to our claim, instead, the network that is trained on the union of the training samples $\{\mathcal{D}_{i,\to}\}_{i}\cup\{\mathcal{D}_{i,\leftarrow}\}_{i}$ from the forward and backward models (50:50 ratio in Figure 3(a)) displays the same test SHD (around $0.5$) as a random classifier assigning the causal direction with equal probability.

Further, we investigate CSIvA’s relation with known identifiability theory by training and testing the architecture on data from a linear Gaussian model, which is well-known to be unidentifiable. Not surprisingly, the results of Figure 3(b) show that none of the algorithms that we learn can infer the causal order of linear Gaussian models with test SHD any better than a random baseline.

Figure: (a) Invertible data
Refer to caption: x7.png

###### Hypothesis 1 .

###### Example 2 .

#### Implications.

Our experiments show that CSIvA learns algorithms that closely follow identifiability theory for causal discovery. In particular, while the method itself does not require explicit assumptions on the data-generating process, the chosen training data ultimately determines the class of causal models identifiable during inference.
Notably, previous work has argued that supervised learning approaches in causal discovery would help with "dealing with complex data-generating processes and greatly reduce the need of explicitly crafting identifiability conditions a-priori", Lopez-Paz et al. []. In the case of CSIvA, this expectation does not appear to be fulfilled, as the assumptions still need to be encoded explicitly in the training data. However, this observation opens two new important questions: (1) Can we train a single network to encompass multiple (or even all) identifiable causal structures? (2) How much ambiguity might exist between these identifiable models?

### 3.4 A low-dimensions argument in favor of learning from multiple causal models

Example 2 of the previous section shows that elements of distinct classes of identifiable structural causal models, such as LiNGAM and nonlinear ANM, may become non-identifiable when we consider their union. In this section, we show that in the class of post-additive noise models given by equation (2) (obtained as the union of the LiNGAM, the nonlinear ANM, and the post-nonlinear model), the set of distributions that is non-identifiable is negligible. Our proposition extends the results of Hoyer et al. [], which are limited to the case of linear and nonlinear additive noise models, and Zhang and Hyvärinen [], which provides the conditions of identifiability of the post-ANM without bounding the set of non-identifiable distributions.

Let $X,Y$ be a pair of random variables generated according to the causal direction $X\rightarrow Y$ and the post-additive noise model structural equation:

$$ $Y=f_{2}(f_{1}(X)+N_{Y}),$ (5) $$

where $N_{Y}$ and $X$ are independent random variables, and $f_{2}$ is invertible.
If the SCM is non-identifiable, the data-generating process can be described by a backward model with the structural equation:

$$ $X=g_{2}(g_{1}(Y)+N_{X}),$ (6) $$

$N_{X}$ independent from $Y$, and $g_{2}$ invertible.
We introduce the random variables $\tilde{X},\tilde{Y}$, such that the forward and backward equations can be rewritten as

$$ $\begin{split}Y=f_{2}(\tilde{Y}),\hskip 10.00002pt\tilde{Y}\coloneqq f_{1}(X)+N _{Y},\\ X=g_{2}(\tilde{X}),\hskip 10.00002pt\tilde{X}\coloneqq g_{1}(Y)+N_{X}.\\ \end{split}$ $$

We note that this implies that the following invertible additive noise models on $\tilde{X},\tilde{Y}$ hold:

$$ $\displaystyle\tilde{Y}=h_{Y}(\tilde{X})+N_{Y},\hskip 10.00002pth_{Y}\coloneqq f _{1}\circ g_{2},$ (7) $\displaystyle\tilde{X}=h_{X}(\tilde{Y})+N_{X},\hskip 10.00002pth_{X}\coloneqq g _{1}\circ f_{2}.$ (8) $$

###### Proposition 1 (Adapted from Hoyer et al. [ 7 ] ) .

Let $p_{N_{Y}},h_{X},h_{Y}$ be fixed, and define $\nu_{Y}\coloneqq\log p_{N_{Y}}$, $\xi\coloneqq\log p_{\tilde{X}}$. Suppose that $p_{N_{Y}}$ and $p_{\tilde{X}}$ are strictly positive densities, and that $\nu_{Y},\xi,f_{1},f_{2},g_{1},$ and $g_{2}$ are three times differentiable. Further, assume that for a fixed pair $h_{Y},\nu_{Y}$ exists $\tilde{y}\in\mathbb{R}$ s.t. $\nu^{\prime\prime}_{Y}(\tilde{y}-h_{Y}(\tilde{x}))h^{\prime}_{Y}(\tilde{x})\neq
0$ is satisfied for all but a countable set of points $\tilde{x}\in\mathbb{R}$. Then, the set of all densities $p_{\tilde{X}}$ of $\tilde{X}$ such that both equations (5) and (6) are satisfied is contained in a 2-dimensional space.

###### Proposition 1 (Adapted from Hoyer et al. [ 7 ] ) .

#### Implications.

Our result is closely related to Theorem 1 of Hoyer et al. [], which we simply generalize to the post-ANM. Intuitively, it says that the space of all continuous distributions such that the bivariate post-ANM is non-identifiable is contained in a 2-dimensional space. As the space of continuous distributions of random variables is infinite-dimensional, we conclude that the post-ANM is generally identifiable, which suggests that the setting of Example 2 is rather artificial. Our results provide a theoretical ground for training causal discovery algorithms on datasets generated from multiple identifiable SCMs.
This is particularly appealing in the case of CSIvA, given the poor OOD generalization ability observed in our experiments of Section 3.2.

### 3.5 Can we train CSIvA on multiple causal models for better generalization?

In this section, we investigate the benefits of training over multiple causal models, i.e. on samples generated by a combination of classes of identifiable SCMs characterized by different mechanisms and noise terms distribution. Our motivation is as follows: given that our empirical evidence shows that CSIvA is capable of in-distribution generalization, whereas dramatically degrades the performance when testing occurs out-of-distribution, it is thus desirable to increase the class of causal models represented in the training datasets. We separately study the effects of training over multiple mechanisms and multiple noise distributions and compare the testing performance against architectures trained on samples of a single SCM.

Figure: (a) Linear test data
Refer to caption: x9.png

#### Mixture of causal mechanisms.

We consider four networks optimized by training of CSIvA on datasets generated from pairs (or triples) of distinct SCMs, with fixed mlp noise and which differ in terms of their mechanisms type: linear and nonlinear; nonlinear and post-nonlinear; linear and post-nonlinear; linear, nonlinear and post-nonlinear. The number of training datasets for each architecture is fixed ($15000$) and equally split between the causal models with different mechanism types.
The results of Figure 4 show that the networks trained on mixtures of mechanisms all present significantly better test SHD compared to CSIvA models trained on a single mechanism type. We find that learning on multiple SCMs improves the SHD from $\sim 0.5$ to $\sim 0.2$ both on linear and nonlinear test data (Figures 4(a) and 4(b)), and even better accuracy is achieved on post-nonlinear samples, as shown in Figure 4(c).

#### Mixture of noise distributions.

Next, we analyze the test performance of three CSIvA networks optimized on samples from structural causal models that have different distributions for their noise terms, while keeping the mechanism types fixed. Figure 5 shows that training over different noises (beta, gamma, gumbel, exponential, mlp, uniform) always results in a network that is agnostic with respect to the noise distributions of the SCM generating the test samples, always achieving SHD $<0.1$, with the exception of datasets with mlp error terms ($0.2$ average SHD on nonlinear and pnl data).

Figure: (a) Linear test data
Refer to caption: x13.png

#### Implications.

We have shown that learning on mixtures of SCMs with different noise term distributions and mechanism types leads to models generalizing to a much broader class of structural causal models during testing. Hence, combining datasets generated from multiple models looks like a promising framework to overcome the limited out-of-distribution generalization abilities of CSIvA observed in Section 3.2.
However, it is easier to incorporate prior assumptions on the class of causal mechanisms (linear, non-linear, post-non-linear) compared to the noise distributions (which are potentially infinite). This introduces a trade-off between amortized inference and classical methods for causal discovery: for example, RESIT, NoGAM, and CAM [, , ] algorithms require no assumptions on the noise type, but only work for a limited class of mechanisms (nonlinear).

## 4 Conclusion

In this work, we investigate the interplay between identifiability theory and supervised learning for amortized inference of causal graphs, using CSIvA as the ground of our study. Consistent with classical algorithms, we demonstrate that good performance can be achieved if (i) we have a good prior on the structural causal model generating the test data (ii) the setting is identifiable. In particular, prior knowledge of the test distribution is encoded in the training data in the form of constraints on the structural causal model underlying their generation. With these results, we highlight the need for identifiability theory in modern learning-based approaches to causality, while past works have mostly disregarded this connection. Further, our findings provide the theoretical ground for training on observations sampled from multiple classes of identifiable SCMs, a strategy that improves test generalization to a broad class of causal models. Finally, we highlight an interesting new trade-off regarding identifiability: traditional methods like LiNGAM, RESIT, and PNL require strong restrictions on the structural mechanisms underlying the data generation (linear, nonlinear or post-nonlinear) while generally being agnostic relative to the noise terms distribution. Training on mixtures of causal models instead offers an alternative that is less reliant on assumptions on the mechanisms, while incorporating knowledge about all possible noise distributions in the training data is practically impossible to achieve. We leave it to future work to reproduce our analysis on a wider class of architectures, as well as extending our study to interventional data with more than two nodes.

## Appendix A Learning to induce: causal discovery with transformers

### A.1 A supervised learning approach to causal discovery

First, we describe the training procedure for the CSIvA architecture, which aims to learn the distribution of causal graphs conditioned on observational and/or interventional datasets. We omit interventional datasets from the discussion as they are not of interest to our work. Training data are generated from the joint distribution $p_{\mathcal{G},\mathcal{D}}$ between a graph $\mathcal{G}$ and a dataset $\mathcal{D}$. First, we sample a set of directed acyclic graphs $\{\mathcal{G}^{i}\}_{i=1}^{n}$ with nodes $X_{1},\ldots,X_{d}$, from a distribution $p_{\mathcal{G}}$. Then, for each graph we sample a dataset of $m$ observations of the graph nodes $\mathcal{D}^{i}=\{x_{1}^{j},\ldots,x_{d}^{j}\}_{j=1}^{m}$, $i=1,\ldots,n$. Hence, we build a training dataset $\{\mathcal{G}^{i},\mathcal{D}^{i}\}_{i=1}^{n}$.

The CSIvA model defines a distribution $\hat{p}_{\mathcal{G}|\mathcal{D}}(\cdot;\Theta)$ of graphs conditioned on the observational data and parametrized by $\Theta$. Given an invertible map $\mathcal{G}\mapsto A$ from a graph to its binary adjacency matrix representation of $d\times d$ entries (where $A_{ij}=1$ iff $X_{i}\rightarrow X_{j}$ in $\mathcal{G}$), we consider an equivalent estimated distribution $\hat{p}_{A|\mathcal{D}}(\cdot;\Theta)$, which has the following autoregressive form:

$$ $\hat{p}_{A,\mathcal{D}}(A|\mathcal{D};\Theta)=\prod_{l=1}^{d^{2}}\sigma(A_{l}; \rho=f_{\Theta}(A_{1},\ldots,A_{l-1},\mathcal{D})),$ $$

where $\sigma(\cdot;\rho)$ is a Bernoulli distribution parametrized by $\rho$. $\rho$ itself is a function of $f_{\Theta}$ defined by the encoder-decoder transformer architecture, taking as input previous elements of the matrix $A$ (here represented as a vector of $d^{2}$ entries) and the dataset $\mathcal{D}$. $\Theta$ is optimized via maximum likelihood estimation, i.e. $\Theta^{*}=\operatorname*{argmin}_{\Theta}-\mathbf{E}_{\mathcal{G},\mathcal{D}
}[\ln\hat{p}(\mathcal{G}|\mathcal{D};\Theta)]$, which corresponds to the usual cross-entropy loss for the Bernoulli distribution. Training is achieved using stochastic gradient descent, in which each gradient update is performed using a pair $(\mathcal{D}^{i},A^{i})$, $i=1\ldots,d$. In the infinite sample limit, we have $\hat{p}_{\mathcal{G}|\mathcal{D}}(\cdot;\Theta^{*})=p_{\mathcal{G}|\mathcal{D}
}(\cdot)$, while in the finite-capacity case, it is only an approximation of the target distribution.

### A.2 CSIvA architecture

In this section, we summarize the architecture of CSIvA, a transformer neural network that can learn a map from data to causally interpreted graphs, under supervised training.

#### Transformer neural network.

Transformers [] are a popular neural network architecture for modeling structured, sequential data data. They consist of an encoder, a stack of layers that learns a representation of each element in the input sequence based on its relation with all the other sequence’s elements, through the mechanism of self-attention, and a decoder, which maps the learned representation to the target of interest. Note that data for causal discovery are not sequential in their nature, which motivates the adaptations introduced by Ke et al. [] in their CSIvA architecture.

#### CSIvA embeddings.

Each element $x_{i}^{j}$ of an input dataset is embedded into a vector of dimensionality $E$. Half of this vector is allocated to embed the value $x_{i}^{j}$ itself, while the other half is allocated to embed the unique identity for the node $X_{i}$. We use a node-specific embedding because the values of each node may have very different interpretations and meanings. The node identity embedding is obtained using a standard 1D transformer positional embedding over node indices. The value embedding is obtained by passing $x_{i}^{j}$, through a multi-layer perceptron (MLP).

#### CSIvA alternating attention.

Similarly to the transformer’s encoder, CSIvA stacks a number of identical layers, performing self-attention followed by a nonlinear mapping, most commonly an MLP layer. The main difference relative to the standard encoder is in the implementation of the self-attention layer: as transformers are in their nature suitable for the representation of sequences, given an input sample of $D$ elements, self-attention is usually run across all elements of the sequence. However, data for causal discovery are tabular, rather than sequential: one option would be to unravel the $n\times d$ matrix of the data, where $n$ is the number of observations and $d$ the number of variables, into a vector of $n\cdot d$ elements, and let this be the input sequence of the encoder. CSIvA adopts a different strategy: the self-attention in each encoder layer consists of alternate passes over the attribute and the sample dimensions, known as alternating attention []. As a clarifying example, consider a dataset $\{(x_{1}^{i},x_{2}^{i})\}_{i=1}^{n}$ of $n$ i.i.d. samples from the joint distribution of the pair of random variables $X_{1},X_{2}$. For each layer of the encoder, in the first step (known as attention between attributes), attention operates across all nodes of a single sample $(x_{1}^{i},x_{2}^{i})$ to encode the relationships between the two nodes. In the second step (attention between samples), attention operates across all samples $(x_{k}^{1},\ldots,x_{k}^{n}),k\in\{1,2\}$ of a given node, to encode information about the distribution of single node values.

#### CSIvA encoder summary.

The encoder produces a summary vector $s_{i}$ with $H$ elements for each node $X_{i}$, which captures essential information about the node’s behavior and its interactions with other nodes. The summary representation is formed independently for each node and involves combining information across the $n$ samples. This is achieved with a method often used with transformers that involves a weighted average based on how informative each sample is. The weighting is obtained using the embeddings of a summary "sample" $n+1$ to form queries, and embeddings of node’s samples $\{x_{i}^{j}\}_{j=1}^{n}$ to provide keys and values, and then using standard key-value attention.

#### CSIvA decoder.

The decoder uses the summary information from the encoder to generate a prediction of the adjacency matrix $A$ of the underlying $\mathcal{G}$. It operates sequentially, at each step producing a binary output indicating the prediction $\hat{A}_{i,j}$ of $A_{i,j}$, proceeding row by row. The decoder is an autoregressive transformer, meaning that each prediction $\hat{A}_{i,j}$ is obtained based on all elements of $A$ previously predicted, as well as the summary produced by the encoder. The method does not enforce acyclicity, although Ke et al. [] shows that in cyclic outputs genereally don’t occur, in practice.

## Appendix B Training details

### B.1 Hyperparameters

In Table 1 we detail the hyperparameters of the training of the network of the experiments. We define an iteration as a gradient update over a batch of $5$ datasets. Models are trained until convergence, using a patience of $5$ (training until five consecutive epochs without improvement) on the validation loss - this always occurs before the $25$-th epoch (corresponding to $\approx 150000$ iterations). The batch size is limited to $5$ due to memory constraints.

**Table 1: Hyperparameters for the training of the CSIvA models of the experiments in Section 3.**
| Hypeparameter | Value |
| --- | --- |
| Hidden state dimension | 64 |
| Encoder transformer layers | 8 |
| Decoder transformer layers | 8 |
| Num. attention heads | 8 |
| Optimizer | Adam |
| Learning rate | $10^{-4}$ |
| Samples per dataset ($n$) | $1500$ |
| Num. training datasets | $15000$ |
| Num. iterations | $<150000$ |
| Batch size | $5$ |

### B.2 Synthetic data

In this section, we provide additional details on the synthetic data generation, which was performed with the causally(^2^22[https://causally.readthedocs.io/en/latest/](https://causally.readthedocs.io/en/latest/)) Python library []. Our data-generating framework follows that of Montagna et al. [], an extensive benchmark of causal discovery methods on different classes of SCMs.

#### Causal mechanisms.

The nonlinear mechanisms of the PNL model and the nonlinear ANM model are generated by a neural network with one hidden layer with $10$ hidden units, with a parametric ReLU activation function. The network weights are randomly sampled according to a standard Gaussian distribution. The linear mechanisms are generated by sampling the regression coefficients in the range $[-3,-0.5]\cup[0.5,3]$.

#### Distribution of the noise terms.

We generated datasets from structural causal models with the following distribution of the noise terms: Beta, Gamma, Gaussian (for nonlinear data), Gumbel, Exponential, and Uniform. Additionally, we define the mlp distribution by nonlinear transformations of gaussian samples from a guassian distribution centered at zero and with standard deviation $\sigma$ uniformly sampled between $0.5$ and $1$. The nonlinear transformation is parametrized by a neural network with one hidden layer with $100$ units, and sigmoid activation function. The weights of the network are uniformly sampled in the range $[-1.5,1.5]$. We additionally standardized the output of each mlp sample by the empirical variance computed over all samples.

Data are standardized with their empirical variance, which removes the presence of shortcuts which could be learned by the network, notably varsortability [] and score-sortability [].

### B.3 Computer resources

Our experiments were run on a local computing cluster, using any and all available GPUs (all NVIDIA). For replication purposes, GTX 1080 Ti’s are entirely suitable, as the batch size was set to match their memory capacity, when working with bivariate graphs. All jobs ran with $10\mathrm{GB}$ of RAM and $4$ CPU cores. The results presented in this paper were produced after $145$ days of GPU time, of which $68$ were on GTX 1080 Ti’s, $13$ on RTX 2080 Ti’s, $11$ on A10s, $19$ on A40s, and $35$ on RTX 3090s. Together with previous experiments, while developing our code and experimental design, we used 376 days of GPU time (for reference, at a total cost of $492.14$ Euros), similarly split across whichever GPUs were available at the time: $219$ on GTX 1080 Ti’s, $38$ on RTX 2080 Ti’s, $18$ on A10s, $63$ on RTX 3090s, $31$ on A40s, and $6$ on A100s.

## Appendix C Further experiments

We present our experimental results on one further question, to help clarify the results in the main text of the paper. Our aim is to understand when to make tradeoffs between computational resources, and having models that have been trained on a wider variety of SCMs. We compare training on multiple SCMs to single-SCM training, when all models see the same amount of training data from each SCM type as a non-mixed model (i.e. a mixed network trains on $15,000$ linear datasets and $15,000$ PNL datasets, instead of $15,000$ divided between the two SCM types).

In the main text of this paper, we compare neural networks trained on a mix of structural causal models (e.g. noise distributions, or mechanism types), to models trained on a single mechanism-noise combination, where all models have the same amount of training data, $15,000$ datasets. In mixed training, we split these evenly, so a "lin, nl" model is trained on $7,500$ datasets from linear SCMs, and $7,500$ from nonlinear SCMs. Our results in this framework are promising, and show that for many combinations of SCM types, we can train one model instead of two, and achieve good progress, while making a $50\%$ savings on training costs. However, if our training budget is high/unlimited, we should also ask whether there is a downside to mixed training - can we achieve the same performance as a model trained on a single SCM type? Fig. 6 shows good results in this direction - the models trained with the same number of datasets per SCM type as an unmixed model had similar (or even better, for PNL data) performance as the un-mixed model trained on the same SCM type as the test data. These mixed models are also significantly more useful than having 2 or 3 separate models per SCM type, as they have good across-the-board performance. However, if we used the same computational resources to train 3 separate networks (one for each mechanism type) and wanted to use them for causal discovery on a dataset with unknown assumptions, we would be left with the rather difficult task of deciding which model to trust.

Figure: (a) Linear test data
Refer to caption: x17.png

## Appendix D Theoretical results and proofs

Before stating the proof of Proposition 1, we show under which condition the pair of random variables $X,Y$ satisfies the forward and backward models of equations (5), (6): this is relevant for our discussion, as the proof of Proposition 1 consists of showing that this condition is almost never satisfied.

#### Notation.

We adopt the following notation: $\nu_{X}\coloneqq\log p_{N_{X}}$, $\nu_{Y}\coloneqq\log p_{N_{Y}}$, $\xi\coloneqq\log p_{\tilde{X}}$, $\eta\coloneqq\log p_{\tilde{Y}}$, and $\pi\coloneqq\log p_{\tilde{X},\tilde{Y}}$.

###### Theorem 1 (Theorem 1 of Zhang and Hyvärinen [ 8 ] ) .

Assume that $X,Y$ satisfies both causal relations of equations (5) and (6). Further, suppose that $p_{N_{Y}}$ and $p_{\tilde{X}}$ are positive densities on the support of $N_{Y}$ and $\tilde{X}$ respectively, and that $\nu_{Y},\xi,f_{1},f_{2},g_{1},$ and $g_{2}$ are third order differentiable. Then, for each pair $(\tilde{x},\tilde{y})$ satisfying $\nu^{\prime\prime}_{Y}(\tilde{y}-h_{Y}(\tilde{x}))h_{Y}(\tilde{x})\neq 0$, the following differential equation holds:

$$ $\xi^{\prime\prime\prime}=\xi^{\prime\prime}\left(\frac{h^{\prime\prime}_{Y}}{h ^{\prime}_{Y}}-\frac{\nu^{\prime\prime\prime}_{Y}h^{\prime}_{Y}}{\nu^{\prime \prime}_{Y}}\right)+\frac{\nu^{\prime\prime\prime}_{Y}\nu^{\prime}_{Y}h^{ \prime\prime}_{Y}h^{\prime}_{Y}}{\nu^{\prime\prime}_{Y}}-\frac{\nu^{\prime}_{Y }(h^{\prime\prime}_{Y})^{2}}{h^{\prime}_{Y}}-2\nu^{\prime\prime}_{Y}h^{\prime \prime}_{Y}h^{\prime}_{Y}+\nu^{\prime}_{Y}h^{\prime\prime\prime}_{Y},$ $$

and $h_{X}$ is constrained in the following way:

$$ $\frac{1}{h^{\prime}_{X}}=\frac{\xi^{\prime\prime}+\nu^{\prime\prime}_{Y}(h^{ \prime}_{Y})^{2}-\nu^{\prime}_{Y}h^{\prime\prime}_{Y}}{\nu^{\prime\prime}_{Y}h ^{\prime}_{Y}},$ (9) $$

where the arguments of the functions have been left out for clarity.

###### Proof of Theorem 1 .

We demonstrate separately the two statements of the theorem.

###### Theorem 1 (Theorem 1 of Zhang and Hyvärinen [ 8 ] ) .

###### Proof of Theorem 1 .

#### Part 1.

Given that equations (5) and (6) hold, this implies that the forward and backward models on $\tilde{X},\tilde{Y}$ of equations (7) and (8) are also valid, namely that:

$$ $\displaystyle\tilde{Y}=h_{Y}(\tilde{X})+N_{Y},$ $\displaystyle\tilde{X}=h_{X}(\tilde{Y})+N_{X}.$ $$

These are the structural equations of two causal models, associated with the forward $\tilde{X}\rightarrow\tilde{Y}$ and backward $\tilde{Y}\rightarrow\tilde{X}$ graphs, respectively. Applying the Markov factorization of the distribution according to the forward direction, we get:

$$ $p_{{\tilde{X},\tilde{Y}}}(\tilde{x},\tilde{y})=p_{{\tilde{Y}|\tilde{X}}}( \tilde{y}|\tilde{x})p_{{\tilde{X}}}(\tilde{x})=p_{{N_{Y}}}(\tilde{y}-h_{Y}( \tilde{x}))p_{{\tilde{X}}}(\tilde{x}),$ $$

which implies

$$ $\pi(\tilde{x},\tilde{y})=\nu_{Y}(\tilde{y}-h_{Y}(\tilde{x}))+\xi(\tilde{x}),$ (10) $$

for any $\tilde{x},\tilde{y}$.
Similarly, the Markov factorization on the backward model implies:

$$ $\pi(\tilde{x},\tilde{y})=\nu_{X}(\tilde{x}-h_{X}(\tilde{y}))+\eta(\tilde{y}).$ (11) $$

From (11), we have that:

$$ $\displaystyle\frac{\partial^{2}}{\partial\tilde{x}^{2}}\pi(\tilde{x},\tilde{y} )=\nu^{\prime\prime}_{X}(\tilde{x}-h_{X}(\tilde{y}))$ $\displaystyle\frac{\partial^{2}}{\partial\tilde{x}\partial\tilde{y}}\pi(\tilde {x},\tilde{y})=-\nu^{\prime\prime}_{X}(\tilde{x}-h_{X}(\tilde{y}))h^{\prime}_{ X}(\tilde{y}),$ $$

which implies

$$ $\frac{\partial}{\partial\tilde{x}}\left(\frac{\frac{\partial^{2}}{\partial \tilde{x}^{2}}\pi(\tilde{x},\tilde{y})}{\frac{\partial^{2}}{\partial\tilde{x} \partial\tilde{y}}\pi(\tilde{x},\tilde{y})}\right)=0.$ (12) $$

Computing the same set of partial derivatives from (10), we find:

$$ $\displaystyle\frac{\partial^{2}}{\partial\tilde{x}^{2}}\pi(\tilde{x},\tilde{y} )=\nu^{\prime\prime}_{Y}(\tilde{y}-h_{Y}(\tilde{x}))(h^{\prime}_{Y}(\tilde{x}) )^{2}-\nu^{\prime}_{Y}(\tilde{y}-h_{Y}(\tilde{x}))h^{\prime\prime}_{Y}(\tilde{ x})+\xi^{\prime\prime}(\tilde{x})$ $\displaystyle\frac{\partial^{2}}{\partial\tilde{x}\partial\tilde{y}}\pi(\tilde {x},\tilde{y})=-\nu^{\prime\prime}_{Y}(\tilde{y}-h_{Y}(\tilde{x}))h^{\prime}_{ Y}(\tilde{x}).$ $$

from which follows:

$$ $\begin{split}\frac{\partial}{\partial\tilde{x}}\left(\frac{\frac{\partial^{2}} {\partial\tilde{x}^{2}}\pi(\tilde{x},\tilde{y})}{\frac{\partial^{2}}{\partial \tilde{x}\partial\tilde{y}}\pi(\tilde{x},\tilde{y})}\right)&=-2h^{\prime\prime }_{Y}+\frac{\nu^{\prime}_{Y}h^{\prime\prime\prime}_{Y}}{\nu^{\prime\prime}_{Y} h^{\prime}_{Y}}-\frac{\xi^{\prime\prime\prime}}{\nu^{\prime\prime}_{Y}h^{ \prime}_{Y}}+\frac{\nu^{\prime\prime\prime}_{Y}\nu^{\prime}_{Y}h^{\prime\prime }_{Y}}{(\nu^{\prime\prime}_{Y})^{2}}-\frac{\nu^{\prime}_{Y}(h^{\prime\prime}_{ Y})^{2}}{\nu^{\prime\prime}_{Y}(h^{\prime}_{Y})^{2}}+\frac{\xi^{\prime\prime} \nu^{\prime\prime\prime}_{Y}h^{\prime\prime}_{Y}}{(\nu^{\prime\prime}_{Y})^{2} \nu^{\prime\prime}_{Y}(h^{\prime}_{Y})^{2}}\\ &=0.\end{split}$ $$

where we drop the input arguments for conciseness. The equality with 0 0 is given by the equality with (12). Manipulating the above expression, the first claim follows.

#### Part 2.

Next, we prove the constraint derived on $h_{X}$. To do this, we exploit the fact that $\tilde{Y}$ is independent of $N_{X}$, which implies the following condition []:

$$ $\frac{\partial^{2}}{\partial\tilde{y}\partial n_{x}}\log p(\tilde{y},n_{x})=0,$ (13) $$

for any $(\tilde{y},n_{x})$. According to equations (7), (8), we have that:

$$ $\displaystyle\tilde{Y}=h_{Y}(\tilde{X})+N_{Y},$ $\displaystyle N_{X}=\tilde{X}-h_{X}(\tilde{Y}),$ $$

such that we can define an invertible map $\Phi:(\tilde{y},n_{x})\mapsto(\tilde{x},n_{Y})$. It is easy to show that the Jacobian of the transformation has determinant $|J_{\Phi}|=1$, such that

$$ $p(\tilde{y},n_{Y})=p(\tilde{x},n_{Y}),$ $$

where $(\tilde{x},n_{Y})=\Phi^{-1}(\tilde{y},n_{X})$. Thus, being $\tilde{X},N_{Y}$ independent random variables, we have that:

$$ $\log p(\tilde{y},n_{X})=\log p(\tilde{x})+\log p(n_{Y})=\xi(\tilde{x})+\nu_{Y} (n_{Y}).$ $$

Given that $\tilde{X}=h_{X}(\tilde{Y})+N_{X}$, we have that

$$ $\frac{\partial^{2}}{\partial\tilde{y}\partial\tilde{n}_{X}}\log p(\tilde{x})= \xi^{\prime\prime}h^{\prime}_{X},$ $$

while $N_{Y}=\tilde{Y}-h_{Y}(\tilde{X})$ implies

$$ $\frac{\partial^{2}}{\partial\tilde{y}\partial\tilde{n}_{X}}\log p(n_{Y})=-\nu^ {\prime\prime}_{Y}h^{\prime}_{Y}+\nu^{\prime\prime}_{Y}h^{\prime}_{X}(h^{ \prime}_{Y})^{2}-\nu^{\prime}_{Y}h^{\prime}_{X}h^{\prime\prime}_{Y},$ $$

such that

$$ $\log p(\tilde{x},n_{Y})=\xi^{\prime\prime}h^{\prime}_{X}+-\nu^{\prime\prime}_{ Y}h^{\prime}_{Y}+\nu^{\prime\prime}_{Y}h^{\prime}_{X}(h^{\prime}_{Y})^{2}-\nu^ {\prime}_{Y}h^{\prime}_{X}h^{\prime\prime}_{Y},$ $$

which must be equal to zero, being equal to the LHS of (13). Thus, we conclude that

$$ $\frac{1}{h^{\prime}_{X}}=\frac{\xi^{\prime\prime}+\nu^{\prime\prime}_{Y}(h^{ \prime}_{Y})^{2}-\nu^{\prime}_{Y}h^{\prime\prime}_{Y}}{\nu^{\prime\prime}_{Y}h ^{\prime}_{Y}},$ $$

proving the claim.
∎

### D.1 Proof of Proposition 1

###### Proof.

Under the hypothesis that equations (5), (6) hold, i.e. when the data generating process satisfy both a forward and a backward model, by Theorem 1 we have that:

$$ $\xi^{\prime\prime\prime}(\tilde{x})=\xi^{\prime\prime}(\tilde{x})G(\tilde{x}, \tilde{y})+H(\tilde{x},\tilde{y}),$ (14) $$

where

$$ $\displaystyle G(\tilde{x},\tilde{y})=\left(\frac{h^{\prime\prime}_{Y}}{h^{ \prime}_{Y}}-\frac{\nu^{\prime\prime\prime}_{Y}h^{\prime}_{Y}}{\nu^{\prime \prime}_{Y}}\right),$ $\displaystyle H(\tilde{x},\tilde{y})=\frac{\nu^{\prime\prime\prime}_{Y}\nu^{ \prime}_{Y}h^{\prime\prime}_{Y}h^{\prime}_{Y}}{\nu^{\prime\prime}_{Y}}-\frac{ \nu^{\prime}_{Y}(h^{\prime\prime}_{Y})^{2}}{h^{\prime}_{Y}}-2\nu^{\prime\prime }_{Y}h^{\prime\prime}_{Y}h^{\prime}_{Y}+\nu^{\prime}_{Y}h^{\prime\prime\prime} _{Y}.$ $$

Define $z\coloneqq\xi^{\prime\prime\prime}$, such that the above equation can be written as $z^{\prime}(\tilde{x})=z(\tilde{x})G(\tilde{x},\tilde{y})+H(\tilde{x},\tilde{y})$. given that such function $z$ exists, it is given by:

$$ $z(\tilde{x})=z(\tilde{x}_{0})e^{\int_{\tilde{x}_{0}}^{\tilde{x}}G(t,y)dt}+\int _{\tilde{x}_{0}}^{\tilde{x}}e^{\int_{\hat{t}}^{\tilde{x}}G(t,y)dt}H(\hat{t},y) d\hat{t}.$ (15) $$

Let $\tilde{y}$ such that $\nu^{\prime\prime}_{Y}(\tilde{y}-h_{Y}(\tilde{x}))h^{\prime}_{Y}(\tilde{x})\neq
0$ holds for all but countable values of $\tilde{x}$. Then, $z$ is determined by $z(\tilde{x}_{0})$, as we can extend equation (15) to all the remaining points. The set of all functions $\xi$ satisfying the differential equation (14) is a 3-dimensional affine space, as fixing $\xi(\tilde{x}_{0}),\xi^{\prime\prime}(\tilde{x}_{0}),\xi^{\prime\prime}(\tilde
{x}_{0})$ for some point $\tilde{x}_{0}$ completely determines the solution $\xi$. Moreover, given $\nu_{Y},h_{X},h_{Y}$ fixed, $\xi^{\prime\prime}$ is specified by (9) of theorem 1, which implies:

$$ $\xi^{\prime\prime}=\frac{\nu^{\prime\prime}_{Y}h^{\prime}_{Y}}{h^{\prime}_{X}} +\nu^{\prime}_{Y}h^{\prime\prime}_{Y}-\nu^{\prime\prime}_{Y}(h^{\prime}_{Y})^{ 2},$ $$

which confines $\xi$ solutions of (14) to a 2-dimensional affine space.
∎

###### Proof.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - Have you explored how your results might extend to larger graphs? What are the main challenges you foresee?
- Can you provide more concrete guidance on choosing optimal ratios when mixing different types of causal models during training?
soundness: 4
strengths: - Establishes a clear connection between how assumptions in classical causal discovery are translated to training data assumptions in amortized approaches.
- Provides strong empirical validation and theoretical justification of generalization in amortized learning improves with mixed training, something that previously has not been studied to this extent, providing an important result and discussion point that will benefit the causal discovery community. 
- Establishes a testing methodology for amortized causal discovery that if adopted stands to benefit the field in improving the methods.
summary: This paper analyzes amortized causal discovery, where synthetic data allows supervised learning of causal discovery. The investigation focuses on whether learning-based methods can overcome traditional identifiability constraints in causal discovery. Extensive out-of-distribution experiments in a bivariate setting with a CSIvA transformer demonstrate that the method cannot generalize well to unseen mechanisms or noise families. This finding establishes an important link to traditional causal discovery literature, revealing how identifiability assumptions are transferred from algorithm design to training data selection. The work further demonstrates empirically that model generalization can be significantly improved by training on diverse families of mechanisms and noise distributions. A theoretical justification for this approach is provided by proving that the set of non-identifiable structures arising from such unions is negligible, particularly extending known results to post-additive noise models.
weaknesses: - Limited Scope:
  - Analysis restricted to bivariate cases
  - Not clear how well the results extend to more complex graphs, nor is there a clear path for extending results to larger graphs
  - Doesn't address practical challenges like confounding or selection bias
- Experimental Limitations:
  - Only evaluated on purely synthetic data
  - Limited exploration of how results scale with dataset size

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: See above weakness.
soundness: 2
strengths: The problem that the authors are trying to attack is actually important. Identifiability is central to causal discovery because it determines whether the true causal structure can be uniquely recovered from the data. Without identifiability, any conclusions about causality are inherently uncertain, making the validation and application of causal models problematic. Meanwhile, models such as CSIvA does not provide any identifiability results, but achieves good performance in practice. Thus to ensure the reliability of such models, identifiability results would be required.
summary: The paper tries to analysis the CSIvA model on bivariate causal models. The authors proposed an hypothesis to attempt to bridge the gap between model free (i.e. end-to-end learning without given knowledge of causal generating process) causal discovery algorithms and the identifiability theory, but only use some examples, without a proper proof.
weaknesses: 1. The hypothesis in the paper is not proved.
2. Actually, for a causal discovery algorithm, if it is differentiable at each steps, then it would be trivial to use a RNN to mimic  the action of the algorithm exactly. Then a large enough transformer can be used to fit the RNN easily. Thus it would be trivial that a transformer would perform good in identifiable causal discovery cases if such algorithm exists, For example, for linear gaussian with equal variance, each steps of NOTEARS is differentiable if the objective is differentiable. Thus we can use a differentiable neural network to mimic NOTEARS. In other cases, such explanation may also holds.

3. The experimental on 2-variables dataset are too simple.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: See above.
soundness: 2
strengths: Some insights from the experiments are provided.


1) For in-distribution generalization, CSIvA effectively generalizes to unseen samples from the training distribution.

2) For out-of-distribution generalization, it is unsurprising that CSIvA performs poorly.

3) CSIvA's performance depends on whether the data is sampled from an identifiable SCM.

4) Training CSIvA on mixtures of SCMs leads to improved generalization.
summary: This work explores supervised learning for causal discovery, a promising direction to potentially replace traditional score-based or constraint-based methods. However, many gaps remain unclear when considering the supervised learning paradigm for this task. To clarify these gaps, we conduct extensive experiments using a specialized method, CSIvA, on a variety of synthetic datasets. The results of these experiments are presented, along with a detailed discussion of their implications.
weaknesses: I am generally open to, encouraging, and expecting new paradigms that go beyond traditional score-based or constraint-based methods for causal discovery. However, this work raises several key concerns, prompting me to review this draft with caution.

1) As an experimental paper, it is expected to provide a wide range of existing methods applied to various datasets, not just simulated data. However, the experimental settings in this work are relatively simple and do not appear to be robust enough for a comprehensive experimental study.

2) The insights from the experiments do not yield particularly exciting results and seem somewhat predictable. I would expect a deeper discussion of these experiments or even some theoretical progress. For example, it would be helpful to explore why training data sampled from identifiable SCMs performs well or how the learned model (e.g., a transformer) leverages this information. Additionally, while the learned model may not explicitly satisfy identifiable conditions (such as those in linear non-Gaussian additive noise models), it would be interesting to investigate whether the model implicitly learns such information.

3) For data sampled from a non-identifiable SCM, there may be multiple graphs within a single Markov equivalence class. In this case, do supervised learning paradigms truly offer no advantage? For instance, could we apply a Bayesian approach to the model to yield the most probable result? Additionally, given the existence of multiple possible graphs, could we frame this as a one-to-many problem, perhaps using techniques like Mixture Density Networks?

Overall, although the insights from this work are not surprising, it is still interesting and valuable, and I believe it contributes to the community. However, the current version is not yet ready for publication. I encourage the authors to incorporate more robust experimental results, provide deeper insights, and offer rigorous justifications.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: I would appreciate clarification on how the method handles independent variable pairs. Does it require a training dataset of independent variables, or is it primarily designed to identify a DAG from the Markov equivalent class, meaning independent cases don't need to be considered? Addressing this point explicitly could be helpful for readers with varying levels of expertise in the field.
soundness: 3
strengths: The paper is well-organized and clearly presented. The authors articulate their points effectively, with implications interwoven throughout the sections and the corresponding background provided in the appendix. The experiments on synthetic datasets are comprehensive, demonstrating that the performance of CSIvA is non-trivial. The empirical results support the proposed hypothesis. The experiment demonstrates that CSIvA’s generalization ability can be significantly enhanced by incorporating multiple distributions into the training dataset. 

This paper potentially addresses the limitations of supervised learning in causal discovery, particularly the challenge of applying pre-trained models to datasets that may be out-of-distribution. Specifically,  the highlight of this paper is that it offers a possible way to make supervised learning-based causal discovery practical: if the training dataset is sufficiently comprehensive, the model’s predictions are somewhat reliable.
summary: The paper presents a comprehensive analysis of CSIvA in the context of bivariate causal models. The authors first demonstrate that while CSIvA generalizes well in distribution, it fails to generalize out-of-distribution, regardless of changes in mechanisms or noise distributions. By hypothesizing that CSIvA can identify the SCM underlying the data generation process, the authors challenge the previous claim that “supervised learning methods can address complex data generation and reduce the need for explicit assumptions.” Additionally, the paper takes an empirical step forward by showing that CSIvA has the potential to capture multiple identifiable causal structures within a single network, supported by theoretical claims that post-ANM models are generally identifiable.
weaknesses: In my view, the largest weakness is the lack of evaluation on real-world datasets. Since the experiments already demonstrate the potential of achieving a relatively robust predictor by incorporating multiple distributions in the training dataset, an investigation and comparison with traditional methods on any real-world datasets would be valuable and appreciated. For example, https://www.cmu.edu/dietrich/causality/projects/causal_learn_benchmarks/. This is just an example, the author does not need to be constrained by this and can use any dataset that is valid for causal discovery. If there are any reasons why experiments on real-world datasets may not be feasible or practical, I would be very interested in hearing them and would welcome the explanation with an open mind.



Another limitation, which is also discussed by the authors, is the scope of analysis is rather limited in a bivariate model without considering intervention.

Meta-review: While the manuscript targets an important question that why amortized causal discovery with transformers works, major concerns raised by the reviewers—particularly regarding the limited scope of experiments (mostly bivariate) and the lack of theoretical claims—still remain. 

Two reviewers were marginally positive but did not fully endorse the paper, while two others maintained their skepticism. Given the discussion, the additional experiments, and the persistent reservations, the balance of arguments weighs slightly against acceptance at this time. The paper provides a valuable perspective and represents a step forward in understanding amortized causal discovery, but it may require further theoretical development or more robust multivariate evidence to be fully convincing. Therefore, the final recommendation is to reject the paper in its current form.

Decision: Reject