Identifying Causal Direction via Variational Bayesian Compression
Authors: Quang-Duy Tran
Bao Duong
Phuoc Nguyen
Thin Nguyen
Venue: ICML 2025 spotlightposter (2025)
Keywords: causal discovery
neural networks
variational Bayesian code
Primary area: general_machine_learning->causality
TL;DR: A framework leveraging variational Bayesian compression of neural networks for the cause-effect identification task

Abstract: Telling apart the cause and effect between two random variables with purely observational data is a challenging problem that finds applications in various scientific disciplines. A key principle utilized in this task is the algorithmic Markov condition, which postulates that the joint distribution, when factorized according to the causal direction, yields a more succinct codelength compared to the anti-causal direction. Previous approaches approximate these codelengths by relying on simple functions or Gaussian processes (GPs) with easily evaluable complexity, compromising between model fitness and computational complexity. To address these limitations, we propose leveraging the variational Bayesian learning of neural networks as an interpretation of the codelengths. This allows the improvement of model fitness, while maintaining the succinctness of the codelengths, and the avoidance of the significant computational complexity of the GP-based approaches. Extensive experiments on both synthetic and real-world benchmarks in cause-effect identification demonstrate the effectiveness of our proposed method, showing promising performance enhancements on several datasets in comparison to most related methods.

Identifying Causal Direction via Variational Bayesian Compression
Quang-Duy Tran
Bao Duong'
Phuoe Nguyen' Thin Nguyen'
Abstract
medical research (Guyon et al., 2019), they often imprac- tical due to resource constraints and ethical considerations. Studying passive observations offers viable alternative for causal inference, despite requiring assumptions on the data generating process to detect the asymmetry between the causal and anti-causal directions. One intuitive interpreta- tion of the causal asymmetry from an information -theorettc perspective is the independence postulate of algorithmic Markov kernels (Janzing & Schölkopf, 2010), which states that the true causal direction must yield the lowest algo- rithmic Kolmogorov complexity factorization of the joint distribution.
Telling apart the cause and effect between two random variables with purely observational data is a challenging problem that finds applications in various scientific disciplines. A key principle utilized in this task the algorithmic Markov con- dition, which postulates that the joint distribution, when factorized according the causal direction, yields more succinet codelength compared to the anti-causa direction. Previous approaches ap- proximate these codelengths by relying on simple functions or Gaussian processes (GPs) with eas- ily evaluable complexity, compromising between model fitness and computational complexity. To address these limitations, we propose leveraging the variational Bayesian learning of neural net- works as an interpretation of the codelengths. This allows the improvement of model fitness, while maintaining the succinetness of the codelengths, and the avoidance of the significant computational complexity of the GP-based approaches. Exten- sive experiments both synthetic and real-world benchmarks in causc-effect identification demon- strate the effectiveness of our proposed method, showing promising performance enhancements on several datasets in comparison to most related methods.
However, because of the incomputability of the Kolmogorov complexity (Li & Vitányi. 2019). approximation meth- ods via the principle of Minimum Description Length (MDL, Grünwald, 2007: Marx & Vreeken, 2017; 2019a;b) are proposed to instantiate this complexity empirically. The two-part MDL aims to find the model that minimize two criteria: (1) the complexity of the data given that model and (2) the model complexity of the model used for modeling the data. The former complexity measures the model fitness, which is commonly evaluated through the log-likelihood of the data given the model. The options for estimating the latter model complexity are more diverse, which usually involve the number of available models and the parameters in each model. As these methods attempt to minimize the total codelength, they are also called compression- based methods.
While these approaches have shown promising results, the estimation of the conditional distributions in these methods relies on traditional regression methods that offer easily evaluable model complexity. If the ground truth conditional models are more complex or too distinct from the predefined ones, this implementation for the model classes can result in lower fitness and higher complexity of data given the model, leading to suboptimal approximations for the algorithmic complexity. Dhir al. (2024a) overcome this restriction by leveraging Gaussian processes (GPs), though this involves significant trade-off between flexibility and computational complexity.
1. Introduction
Cause-effect identification or bivariate causal discovery- the task of telling apart the cause and the effect between two random variables-is a critical task across various sci- entific disciplines. including biology, economics, and so- ciology (Pearl, 2009). While randomized controlled tri- als (RCTs) are considered the most accurate method for identifying these types of causal relationships, especially in
'Dekkin Applied Artificial Intelligence Initiative, Deakin Uni- versity, Geelong, Australia. Correspondence to: Quang- Duy Tran <q.tran @ edu..
To address the limitations of these previous compression- based methods, we propose the leverage of neural networks for learning the conditional models, which are regarded
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Identifying Causal Direction via Variational Bayesian Compression
as universal approximators (Hornik al., 1989). More- over, the algorithmic complexity of the networks can be approximated empirically via the concepts of bits-back cod- ing (Hinton & van Camp, 1993; Wallace, 1990) and varia- tional Bayesian coding (Honkela & Valpola, 2004; Louizos et al., 2017; Blier & Ollivier, 2018). In this study, we intro- duce COMIC- a Bayesian C0Mpression-basee approach to Identifying the Causal direction that improves the model fitness of compression-basee based methods without compromis- ing the computability of model complexity and avoids the higher computational complexity of GP-based modeling. Correspondingly, through extensive empirical evaluation on benchmarks for the cause-effect identification task, our ap- proach demonstrates promising results with performance im- provements compared to a majority of complexity-bee and maximum likelihood-based approaches on multiple bench- marks.
Bühlmann et al., 2014; Peters al., 2014), where the ef- feet Y is generated from function of the cause X and independent noise By (X TT Ey) by adding them as Y = f (X) + Ey. In this model, the cause is assumed to only contribute to the mean, which is can be estimated by mean regression methods (Shimizu et al., 2006; Hoyer et al., 2008; Bühlmann et al., 2014; Peters et al., 2014). From the estimated models, CAM (Bühhmann et al., 2014) de- termines the causal direction by selecting the one with the greater maximum likelihood, whereas RESIT (Peters et al., 2014) quantifies the independence between the cause and the estimated noise with the Hilbert-Schmidt Independence Criterion (HSIC, Gretton et al., 2005).
Due to the striet assumption on model classes of ANMs, generalization approaches for the ANMs have been in- troduced to allow for more complex functions, including post nonlinear models (PNL.s, Zhang & Hyvärinen, 2009) and heteroscedassic/llcattino Immer et al., 2023). Specifically, LSNMs assume that the cause not only contributes to the mean but also the scale through another function g (X), resulting the effect Y == f (X) + g (X) x Ey. These models are more flexible than the ANMs due to their ability to also cover multiplica- tive noise models with the scale functions. LOCI (Immer et al., 2023) chooses the Gaussian likelihood to estimate the mean and scale functions, and recovers the noise from the fitted models. Similar to RESIT, LOCI also consider HSIC as a criterion in addition to the likelihood for predicting the causal direction. Since the Gaussian likelihood may not be robust to epistemic uncertainty, ROCHE (Tran et al., 2024a) suggested more robust estimation for LSNMs by replacing the Gaussian likelihood with a likelihood based on Student's t- t-distribution.
noise models (LSNMs,
Contributions The key contributions of this work can be outlined as follows:
1. We propose the utilization of Bayesian neural networks for modeling the conditional distributions to address the challenge of balancing flexibility and scalability, which hindering previous complexity/oompresion based methods. By minimizing the variational Bayesian codelength, can approximate the algo- rithmic complexity of neural networks.
2. From our proposed encoding scheme of the data given each causal direction, the causal identifiabil- ity can be proven. In particular, our models are non- eeparable-compatibbe which implies that given a suffi cient amount of data, the causal direction can be iden- tified via the complexity of our models.
Principle of Independent Causal Mechanisms Beside the algorithmic complexity. there are other approaches for interpreting the principle of independent causal mechanisms (ICMs, Peters et al., 2017, Sec. 2. which assumes the independence between the marginal distribution of the cause and the conditional distribution of the effect given the cause. With the assumption of low noise levels and invertible causal mechanisms, IGCI (Daniušis et al., 2010) formulates this independence using orthogonality in information space to distinguish cause and effect, which is implemented by the relative entropy distances. CDCI (Duong & Nguyen, 2022) postulates that due to this ICM principle, the shape of the conditional distributions will be invariant, and compute the variations in shape to find the causal directions.
3. The capability of our approach is assessed on both syn- thetic and real-world bivariate causal discovery bench- marks. In comparison to most related approaches, our method achieves performance improvements in several benchmarks, demonstrating the effectiveness of our approach in identifying the causal direction.
2. Related Works
Despite being a well-defined task, the amount of information for determining the causal direction in the bivariate setting is limited, hindering further improvements in both theoretical and empirical results. In the following section, we provide an overview of recent related publications about this task, which includes two popular approaches.
The algorithmic complexity
methods (Marx & Vreeken, 2019a;b; Tagasovska et al., 2020) do not only consider the model fitness objective as in FCM-based meth- ods but also examine the complexity of the model. SLOPE and SL.OEER (Marx & Vreeken, 2017; 2019b) use a set of
Functional Causal Models The earliest functional causal model (FCM, Pearl, 2009) being proposed is the additive noise model (ANM, Shimizu et al., 2006; Hoyer et al., 2008;
2

---

Identifying Causal Direction via Variational Bayesian Compression
basis functions regress the data globally and locally to account for both deterministic and non-deterministic fune- tions, and compute the codelengths for encoding the data
where K K() is the Kolmogorov complexity, PAG., denotes the parents of X, in the causal hypothesis G, and = denotes equality up to a constant, which is independent of p (-).
with two-part codes. SLOPPY (Marx & Vreeken, 2019a) is an improvement of RECI (Blöbaum et al., 2018) that utilizes regularized regressions of Identifiable Regression- based Scoring Functions and find the one with the lowest regularized score to find the minimal model in each direc- tion. QCCD (Tagasovska et al., 2020) uses non-parametric conditional quantile regression methods and encodes the data via each quantile model. Our (inK--COMC-ala be- longs to this category where the conditional distributions of the data are modeled by neural networks and encoded by the variational Bayesian coding scheme. This approach allows for more flexibility compared to previous methods while allowing for the balance between model fitness and model complexity. Another study by Dhir et al. (2024a) interprets the principle of ICMs from the view of the Bayesian model selection and proposes using the marginal likelihoods esti- mated by latent variable Gaussian processes (GPLVM, Tit- sias & Lawrence, 2010) for identifying the causal direction.
In the bivariate setting with two random variables X and Y, if X causes Y (denoted as X + Y), Eq. (1) will become
K (o = K (p (x)) + K (p (Y X)).
(2)
Following the algorithmically independent conditionals, the algorithmic independence of Markov kernels has also been postulated by Janzing & Schölkopf (2010); Mooij et al. (2010) as follows
Postulate 2 (Algorithmic Independence of Markov Kemels, Janzing & Schölkopt, 2010). If + YY, the marginal dis- tribution of the cause p (X) and the conditional distribution of the effect given the cause p (Y) | X) are algorithmically independent of each other. In other words, their algorith- mic mutual information (Ix) will be equal zero up to an additive constant.
IA (p PY xx)=
3. Preliminaries
and this independence does not hold in the other direction.
In this work, the assumption of causal sufficiency is adopted, which is similar to previous publications (Immer et al., 2023; Marx& & Vreeken, 2017; 2019a;b; Mooij et al., 2016; Peters et al., 2014; Tagasovska et al., 2020; Tran et al., 2024a). This means that we assume that there is no hidden confounder between the two random variables. In other words, given two variables X and Y, if they not independent, either X or Y will be the cause of the other.
From these postulates, Mooij et al. (2010, Thm. 1) induce a rule for identifying the causal direction.
Theorem 3.1 (Asymmetry in Complexities of Markov Ker- nels, Mooij et al., 2010). If X is the cause of Y and Pos. 2 holds, the description of the joint distribution K (p Y)) via the description of the marginal distribution of the cause K (p (X)) and the description of the conditional distribu- tion of the effect given the cause K (p | [X) is the most succinct one, or formally,
As mentioned in Sec. 2, interpreting the principle of indepen- dent causal mechanisms is one of two major approaches for determining the direction of a causal relation. The algorith- mie interpretation of this independence via the Kolmogorov complexity originates from the postulate about the algorith- mic independence of conditionals (Janzing & Schölkopf, 2010, Eq. 26) as follows:
K (x)())) + K X)) KWx
As a consequence of this rule, causal indicator score can be obtained for the assumed causal directions of X + Y as follows
Postulate (Algorithmic Independence of Conditionals, Janzing& & Schölkopf, 2010). Let G be causal hypothe- sis, represented by a directed acyclic graph (DAG), over a sei of d variables X1,...,Xa with joint density (XI..... Xa). which is lower semi-computable, that is, K (p (XI,..., Xa)) < co. The causal hypothesis is only acceptable if the shortest description (i.e., the Kolmogorov complexity) of the joint density K (p (X1,...,xa) is equal to concatenation of the shortest description of the Markov kernels up to independent additive constant. This postu- late can be described formally as
Ax-Y = K (x))) + K (Y Xx),
(5)
and vice versa for the remaining direction of V - From these indicator scores, we can infer that X + Yfff A > and Y + X if Ay - Ax- <
The Kolmogorov complexity is not computable in prac- tice (Li & Vitányi, 2019). Hence, the causal indicators in the previous section are substituted by approximating approaches such as Minimum Message Length (MML, Wal- lace& & Freeman, 1987) or Minimum Description Length (MDL, Rissanen, 1978) in previous ceammation-theorete methods (Mooij et al., 2010; Marx & Vreeken, 2017:
d
K (X (X....., ..). +
(1) Σ K (x PAaa)),,
1=1
3

---

Identifying Causal Direction via Variational Bayesian Compression
2019b;a). MML and MDL share a two-part coding principle where the complexity or codelength' T.2PP of the data D is computed via model M E . by combining the complex- ity (fitness) of the data given that model L1 (D) M) and the complexity of the model L2 (M) as
models for computing the codelengths of conditional distri- butions. Moreover, due to their flexibility and capability of approximating wide range of functions, neural networks allow us to only focus on this class of models where com- plexity scores are determined solely by their parameters.
(D) = L1 (D M) + L2 (M).
(6)
The prequential code (Dawid, 1984) and the variational Bayesian code (Honkela & Valpola. 2004) are two effec- tive MDL approaches for encoding the conditional distribu- tion modeled by neural networks (Blier & Ollivier, 2018; Grünwald, 2007; Voita & Titov, 2020). The former encodes the model implicitly through sequential data transmission, as seen in applications such as time series. The latter involves predefining the priors over the parameters and utilizes vari- ational inference to learn the posteriors from the samples. In contrast to the online prequential code, this approach is aligned with the wo-ppart code. Despite these differences in coding strategies, both methods yield consistent results, as noted by Voita & Titov (2020). We opted for the variational Bayesian code to assess the codelengths because it explicitly captures the model complexity and does not necessitate a specific order of data transmission.
The codelength with a model M* that minimizes this equa- tion is appointed as an instantiation for the algorithmic com- plexity. If there are multiple solutions for M the one with the smallest model complexity L, is selected to model the data. From this two-part code, we can attain an approxima- tion for the causal indicator score in Eq. (5) as follows
A 2 A = LM (X) -
2p M 1x),
where M and My x are models that minimize 1 (X) (Y X), respectively.
and 2 x/x
The definitions of Kolmogorov complexity and algorith- mic mutual information, well as the discussion of the MDL-based instantiation of Kolmogorov complexity in the context of causal discovery, as referenced in this section, are provided in App. A.
4.2. Variational Bayesian Code for Evaluating Complexity of Neural Networks
The problem of encoding a conditional distribution P (Y X) is often defined via transmitting perspective.
4. COMIC: Bayesian Compression for Identifying Causal Direction
Alice has dataset DV = { (-²), yy(") } that needs to be transported to Bob, who has already got the input sam-
4.1. Classes of Models for the Conditionals
plcs of this dataset {x(1) N },m The most efficient method is to encode the conditional distribution p pYY X) so that
For the conditional distribution, most previous MDL-related studies choose the set M of candidate functions by predefin- ing a list of basis functions (Marx & Vreeken, 2017; 2019b) or using more advanced regression methods such as cubic spline regression (Marx & Vreeken, 2019a). Although these classes of functions have the model codelengths that are easily computable (e.g., through the number of bits of linear parameters in polynomial regressions), the fitness can be compromised when the ground truth models are more com- plex, leading to suboptimal codelengths. Dhir et al. (2024a) utilize GPLVM (Titsias& & Lawrence, 2010) improve the fitness; however, this substantially increases the computa- tional complexity due to the poor scalability of GPs.
Bob can predict the remaining output part {)()) W fi=1 of DW The variational Bayesian code a two-part code in MDL where both Alice and Bob first designate a class of model M {p x; 0) O E @). where I represents the parameters of the conditional probability density function (y xr, 0), and a prior distribution for the parameters with the probability density function p (0). The corresponding two-part codelength? for this setting can be computed as
p()) (ya:N)
| (o)a))))) = log (y | x (1:N) 0
log (0),
(8)
Neural networks are one class of models that can overcome these limitations thanks to their universality in approxima- tions (Hornik et al., 1989) and better scalability. Addi- tionally, the complexity of neural networks has also been well studied and implemented with different encoding meth- ods (Blier& & Ollivier, 2018; Louizos et al., 2017; Voita & Titov, 2020). Hence, neural network can be viable class of
where the former term corresponds to L1, the latter cor- responds La in Eq. (6), and p (y(1:N) | (,(:.)) 0) == II P (1(6) 2000 0). The variational Bayesian code is
"In most information thhorrettttt literature, the amount of infor- mation measured in bits. As result, the logarithm with the base of 2 is commonly chosen. In this work, for convenience, we will use the natural logarithm and measure the amount of information in nats (natural units of information, Grünwald, 2007). The amount in nats can easily be converted to the corresponding value in bits by dividing it by log 2.
Since algorithmic complexity is defined based on the length of the code, or codelength, we use the terms of "complexity" and "codelength" " interchangeably in this work.
4

---

Identifying Causal Direction via Variational Bayesian Compression
based on the bits-back coding scheme (Wallace, 1990; Hin- ton & van Camp, 1993). In this scheme, Alice employs codelength with redundant code and computes the code- length with respect to that auxiliary information (i.e., the redundant code). After that, Bob will perform the same learning process and observe the choice made by this pro- cess to retrieve the auxiliary information and the transmitted data (Honkela & Valpola. 2004).
The selection of the priors p (0) and the variational dis- tributions 2o (0) over the parameters is necessary step in Bayesian learning. In this work, we employ Gaussian distributions as the family for both the priors p (0) and the mean-field variational posteriors Io (0). The details on these priors and variational posteriors are provided in App. B.
4.3. Identifying Causal Direction via the Codelengths
By applying this scheme, instead of finding point estimate of the parameters 0" that minimize Eq. (8), we introduce re- dundant code by choosing the parameters from a variational distribution Go (0) and compute the expected codelength on this distribution. In this scenario, the amount of redun- dant code to encode 9o (0) is its entropy H (1s (0)) = (0) [- log 9 (0) (Honkela & Valpola. 2004). We can obtain the amount of original information being transmit- ted by deducting the excessive entropy H (go (0)) from the expectation codelength as follows
From the formulation of variational Bayesian codelength of the conditional distribution described in the previous section, we can approximate the conditional codelength for the assumed causal direction X - Y using Eq. (10) with the chosen likelihood distribution being the Gaussian distribution as
P (y(LN)) ( .0)
N III N (y 1 (200) o² (2() 1=1
T, (0)(0) (yOLM) (x(1.)))
(13)
= -Bgo(w) [ (W(UN) (oUUN))))) [0] -H (gs (0)) (9)
where 1s (:: () and (-: are modeled via a neural network fy R θ - IR2 with parameters 0 E O. In particular, we implement neural network with one hidden layer and two output nodes fy, (:)) and fy.2 (: 0). and compute the parameters of the likelihood as follows
= =-Bqs(m) [lggp ( (yaan)) 2..N)), .0]
+ +KL (0) (()).
(10)
where the first term corresponds to the fitness of the data given the model and the second term corresponds to the complexity of the model of the two-part MDL prin- ciple (Honkela & Valpola, 2004; Louizos et al., 2017). -I. -0) (,(1: (.((:))) is known as the evidence lower bound (ELBO) in the variational inference problem (Blei et al., 2017). From this perspective, by minimizing Eq. (10), we can expect to achieve the negative logarithm of the evidence in the Bayesian inference problem, which is also known as the Bayesian codelength in MDL cod- ing (Grünwald, 2007) as follows
(:; fy, (::)) and (::)) = (fr.2 (=))), (14)
where Ç (-) is positive link function, such as the exponen- tial or softplus function, for ensuring the standard deviation values being positive.
Both c(::N) and (11-N) are standardized with respect to their corresponding sample means and standard deviations to avoid the scale-based bias on the identifiability (Reisach et al., 2021). Regarding the codelength of the assumed cause X, we adopt standard Gaussian distribution N' (= 0, 1) to encode the data using the marginal codelength, which is a common choice in previous methods (Mooij et al., 2016; Immer et al., 2023). We denote this marginal codelength Lx (I(1:))).
L Bayes
(2(1:)) 1 ()()))))
= -log/ / (p((i)) (o(li)) 0) p (0) do.
(11)
The approximated causal indicator score for this direction is the sum of the two codelengths
The I. Bayes (y(I:N)
between the optimal codelength (x(1) and its upper bound a.(1:N)) can also be formulated from the variational perspective as follows
gap
yar
(y(I:N)
Å A (DN) = -L (2(2:N)
(0)
+ +or (0) (y(liN) aellin)) (15)
(o) (y(IN)) =((IN))) Bayes (D(LNN) =((NN))) p(o) = KT. 2o (0) ll (0) a(1:,)) v(an))
where gø- (0) is the model that minimizes the variational Bayesian codelength. The corresponding score Ây for the reversed direction Y + X is estimated by a similar procedure. Once the scores in both directions are obtained, the difference between them provides final score:
(12)
where p (0) (1:N) (1:N) is the posterior distribution of the parameters given the observed samples. From this for- mulation, it is obvious that we can retrieve the Bayesian
codelength iff. qo (0) converges to p (0) [<(IN)) (((:,,) (L..))).
A var (DN) = yar (D A AvL (DM), (16)
5

---

of Our Identifying Causal Direction via Variational Bayesian Compression
which indicates the inferred causal direction with its abso- lute value reflecting the confidence of the inference.
Proposition 4.2
Causal Models). Let the two Bayesian causal models Mx-Y and My respectively be factorized into the marginal densities-pl (x Mx-x) and p (y My-x)) and the conditional densities (y x By, Mx-xy and p (x y, 0x, XxxMy-x) The marginal densities are standard Gaussian (. 0,1). The conditional densities are Gaus- sian likelihoods (N (=: (;;0),000) o² (-:)), where i ((:0) and a (::0) are parametrized by neural networks as de- scribed in Eq. (14), Sec. 4.3, and the parameters O follow Gaussian priors specified in App. B. Under these conditions, the Bayesian causal models are not separable-compatte
If the optimized variational distribution 2o- (0) converges to the posterior distribution p (Ø) x=(1:M) yyl:))) and V (0, is the ground truth distribution of p (X), we can expect À X to converge to the Bayesian causal indicator score
Bayes X-Y (DM)= = LN (2((1.))
A
+ L p(0)
(ya:N)) | a
(17)
= - -ogog (DN)) M
(18)
Corollary 4.3 (Causal Identifiability of Our Causal Mod- els). Let p (- Mx-x) ) and (- My respectively be the marginal likelihoods of the data given the causal mod- els Mx-Y and - Because of the non-separable- compatibility in Prop. 4.2, there exists a dataset DN whose causal direction is identifiable via the difference between Bayesian indicator scores:
where p (DN .x--r) is the marginal likelihood of the dataset DV factorized in accordance with the causal model M Conversely, Blayss (D*) can be achieved if Y N (0, 1) and the minimized variational codelength in this case also converges to the Bayesian codelength.
4.4. Causal Identifiability
Bayes (DN) == A
Bayes -X (DM)
- -AX-Y (DN). (19)
The identifiability of our approach is closely related to the identifiability of Bayesian causal models via marginal like- lihoods. First, we introduce the definition of separable- compatibility by Dhir et al. (2024a), which is a necessary condition of two Bayesian causal models being unidentifi- able via their marginal likelihoods, regardless of the dataset DAN
Cor. 4.3 implies that in large sample limit, the Bayesian causal indicator scores defined in Eq. (17) can distinguish the cause and effect from observational data. Details on the causal identifiability via marginal likelihoods and the proof of Prop. 4.2 are further discussed in App. C.
Definition 4.1 (Separable-Compatibin) of Bayesian Causal Models, informally restated from Dhir et al., 2024a). Two causal models Mx-+ and My-x are separable- compatible if the anti-causal factorizations of Mx-x and My-x respectively belong to the same classes of distri- butions as Mx and Mx-y, and the priors of these anti-causal factorizations can also be factorized with respect to the priors of My and Mx-x-
5. Experiments
Throughout this section, the empirical performance of our COMIC approach is evaluated in comparison with state- of-the-art complexity- based and regression- methods for bivariate causal discovery. Implementation details and descriptions of the benchmarks related to this section are provided in App. D and E, respectively.
In this definition, the anti-causal factorization of the causal model Mx-x refers to the factorization of the joint dis- tribution p (X, Y) into p (Y) and p (XX Y) with respect to Mx-Y. Similarly, the anti-causal factorization of the causal model My involves factorizing p (X, Y) into p (X) and p (Y | X) according to My If the anti- causal factorization the causal model Mx-+ results in the same distributions as the causal factorization o My or vice versa, the two causal models are said to be separable- compatible.
5.1. Experimental Settings
Benchmarks Following previous works (Immer et al., 2023; Marx & Vreeken, 2019a; Tagasovska et al., 2020; Tran et al., 2024a), the experiments in this section utilizes both synthetic and real-world data for evaluation. The synthetic data consists of 12 common benchmarks. The first collection of synthetic datasets, including AN, AN-s, LS, LS-s, and MN-U, are proposed by Tagasovska et al. (2020). The next group of simulated benchmarks, compris- ing SIM, SIM-e, SIM-G, and SIM-In, is introduced by Mooij et al. (2016). The remaining synthetic benchmarks consist of the CE-Multi, CE-Net, and CE-Cha datasets described by Goudet et al. (2018). For real-world data, we choose the Tübingen cause-effect pairs (Mooij et al., 2016).
As a result of our Bayesian coding scheme, the identifi- ability of our method is verifiable through an orthogonal perspective of marginal likelihoods (Dhir et al., 2024a). Let us assume that the approximated scores in Eq. (15) would converge to the Bayesian scores in Eq. (17), our results on the non- employed in our method can be presented as follows:
of Bayesian causal models
Baselines
We compare our work against various base- line ICM-based and -ased methods. Methods based
6

---

Identifying Causal Direction via Variational Bayesian Compression
Accuracy
Accuracy AN
Accuracy
Accuracy
AN
AN
AN
AN-
Tülinngnn
AN
Takingen
AN4
1.00 Tültingen
AN-
1.00
Tübinggn
15
D'25
CL Cba
LS
D575
CB Cha
1S
D:25
CD Cha
LS
CE Cha
0.50
D. EV
U.5O
0.50
0,25
15-5
CE-Ne
CB-NNt
MI-U
CE-Mulii
MN-U
CE-Mulii
MINU
CE.Mulli
MN-U
CE-Mulii
SIM-In
c SIM-G
SIM
SIM-G
SIM SIM
SIM SIM-G
Bi-AUROC
Bi-AUROC
Bi-AUROC
Bi-AUROC
AN
AN
AN
AN
AN
1.00 Tübingen
AN
.00
Tübingen
AN
1.00 Tübingen
AN
1.00 Tübingen
CD
LS
Cha
CE Cha
8.00
B.d0
MN-U
CE-Muli
MN-U
CE-Mulii
CF-Muli
MN-E
CE-Mulli
SIM
SIM In
SIM
SIM
SIM
SIM
SIM
SIM In
(a) COMIC (Ours)
(b) SL.PPY
(c) SLOPE
(d) QCCD
(solid: AIC, dashed: BIC)
(solid: SLOPER, dashed: SLOPE)
Accuracy
Accuracy
Accuracy
Accuracy
1.00 Tübingen
Tubingen
Tübingen
1.00 Tühingen
D.50
D,25
b.2b
1.5
BLO
CE-Net
LS
6.00
CE-Ne: LS
6.00
CE-Net LS-s
d.on
CE-Net
CE Mulii
CE Muli
MN
CL Multi
CE Multi
SIM
SIM
SIM
SIM
SIM
SIM-c
SIM-G
SiM-c
SIM-G
SIMe
SIM-G
SIM-c
SIM-G
Bi-AUROC
Bi-AUROC
Bi-AUROC
Bi-AUROC AN
AN
AN
AN
Tühingm
Tuthingen
Tishngen
Tabingen
CF-Chi
1.S
CF-Chi
IS
CF-Clo
1.5
CF-Cla
D,50
D.5O
D.50
DLSO
D,25
0,25
1.5
CE Nei
LS
6.00
CC-Ne:
LS
6-00
CE-Net LS
Net
MN U
CB Multi
MN
CE Muldi
MN U
CE Mula
MN U
CE Multi
STM-In
SIM
STM-In
SIM
STM-h
SIM
SIM-:
SIM-G
SIM-
SIM-G
SIMe
SIM-G
SIMcc
SIM-G
(solid: (e) IGCI
(I) GPLVM
(g) LOCI
(h) CAM
Gaussian)
Figure 1. Performance on synthetic and real- world benchmarks. Higher accuracy and bidirectional AUROC (Bi-AUROC) values are preferable. Our COMIC method is compared against SLOPPY (with the AIC and BIC variants, Marx & Vreeken, 2019a), SLOPER (Marx Vreeken, 2019b) and SLOPE (Marx & Vreeken, 2017). QCCD (Tagasovska et al., 2020), IGCI (with uniform and Gaussian reference measures, Daniušis al., 2010), GPLVM (Dhir al., 2024a), LOCI (Immer al., 2023), and CAM (Bühlmann al., 2014). The marginal likelihood-hee haathoood-haaddddoo objective in COMIC achieves promising results with enhanced performance compared to methods with the maximum likelihood complexity based methods, COMIC also achieves better performance on multiple synthetic datasets and obtains comparable results on real-world data.
objectives LOCI and CAM, especially on real-world Tübingen benchmark. In comparison to most

---

Identifying Causal Direction via Variational Bayesian Compression
1.00
5.2. Results & Discussion
0.75
The experimental results on all aforementioned benchmarks are visualized in Fig. 1 and 2. The overall results in Fig. 2 demonstrate that our COMIC approach achieves second- best performance in comparison to complexity-basee -baeed and maximum rikelihood-basee methods. Our method is also more consistent across all benchmarks, as indicated by the smaller standard deviation of the accuracy and Bi-AUROC scores compared to most baseline methods. Furthermore, COMIC outperforms the maximum likelihood- meth- ods, i.e.. LOCI and CAM, suggesting that its marginal likelihood- based objective delivers promising results with improved performance relative to methods with the maxi- mum likelihood- -ased objectives.
Accuracy
0.50
0.25
0.00
COMIC (Ours)
SLOPPY-AIC
SLOrPy- BIC
SLOPPER
StoPF
IGCI-Uanton Lnttooi
QCCD
IGCT Gaussian
GPLLMM
LOC1
CAM
1.00
ST-AURCC +
0.50
GPLVM is the baseline that attains the best overall results. However, it important to note that GPLVM also requires the most substantial computational resources, with execu- tion GPUs (as noted in App. D), compared to COMIC and other baselines, which can be executed on CPUs. This high resource demand stems from the computational com- plexity of Gaussian processes with latent variables and the extensive use of random restarts for hyperparameter selec- tion.
0.25
0.00
COMIC COMIC (Ohrs)
SLOPY--ACC
Stopfy- BIC
St. SLOPPER
SLOVE
IGCL Unifemm
QCCD
IGCI Garsion
GPLVM
LOCT
CAM
Figure Overall performance of COMIC and the baseline meth- ods on all benchmarks. The accuracy and bidircctional AUROC (Bi-AUROC) scores are averaged over all datasets. The results indicate that our COMIC approach achieves better and more con- sistent overall results comparison to majority baselines.
Performance on Synthetic Benchmarks On synthetic AN, AN-s, LS, LS-s, and MN-U benchmarks. COMIC perfectly determines the causal directions in every case. QCCD, GPLVM, and LOCI also achieve equivalent ac- curacy and Bi-AUROC scores on these datasets. These findings align with theoretical expectation that models iden- tifiable via maximum likelihood are also identifiable via marginal likelihood (Dhir et al., 2024a). IGCI with Gaus- sian reference measure is another method that excels in this group of datasets. The remaining baseline approaches do not perform as well on the more challenging LS, LS-s, and MN-U sets. In particular, the predictions of SL.PPR, SLoPe, IGCI with uniform reference measure, and CAM on these datasets are noticeably suboptimal compared to other baselines.
on the principle of ICMs include SL.PPY (Marx & Vreeken, 2019a), SLOPE (Marx & Vreeken, 2017; 2019b), QCCD (Tagasovska et al., 2020), and IGCI (Daniušis et al., 2010), and GPLVM (Dhir et al., 2024a). For FCM-based methods, we choose CAM (Bühlmann et al., 2014) as a representative for ANM-based methods, and LOCI (Immer al., 2023) with maximum likelihood scoring to represent LSSNM-base methods. Some of the ICM-based methods have different variants, which also include in the evalua- tions.
On the more challenging SIM and SIM-c datasets, COMIC approach demonstrates performance comparable to complexity-basee methods, such as SLOPPY with AIC, SLOPER, and SLOPE. A majority of baseline methods perform decently on SIM-G and SIM-In, which are more manageable compared to the previous two sets. COMIC achieves the accondhiighest scores on SIM-In, with slightly lower accuracy and Bi-AUROC than GPLVM. One unex- pected outlier among the baseline models is IGCI, which exhibits the lowest accuracy and Bi-AUROC on SIM-In, de- spite its intended focus on telling apart causal directions in low- scenarios. Moreover, SLOPPY with AIC, SLOPER, Slope, and QCCD underperform relative to COMIC on the
Evaluation Metrics The identification result of each pair in a dataset considered as a sample in a binary classifi- cation problem of that dataset. Hence, the accuracy score and the area under receiver operating characteristic curve (AUROC) are commonly utilized as evaluation metries for the task of bivariate causal discovery. As the ground truth directions in the benchmarks are imbalanced, we compute the bidirectional AUROC (Bi-AUROC, Guyon et al., 2019, Sec. 2.4.3), which is the average of the forward AUROC and the backward AUROC corresponding to X + Y and y X, respectively.
8

---

Identifying Causal Direction via Variational Bayesian Compression
SIM-G dataset with nearly Gaussian distributions over the causes.
Limitations & Future Work One persistent limitation of our current work is the non-convexity of the learning objec- tive, which may require further investigation into the con- vergence and consistency. Additionally, the use of standard Gaussian codelength to encode the marginal distribution of the cause can introduce a bias toward "more Gaussian" causes, potentially posing a hindrance in intricate settings, In future work, we plan to adapt our approach to multivari- ate settings, explore alternative priors and likelihoods, and consider other marginal likelihood estimation techniques for neural networks, such as Laplace approximation, to en- hance its applicability. Regarding the multivariate extension, we also provide a brief discussion in App. H on potential adaptations of the bivariate methods, including ours, to mul- tivariate data.
The diverse causal modeling of CE-Multi results in varying performance across benchmarks. Information -hheorrttint ap- proaches, including COMIC, SLOPPY with BIC, SLOPER, SLOPE, IGCI with Gaussian reference measure, and GPLVM, tend to perform well on this benchmark. QCCD and uniform-refereeced IGCI, while also belonging to this category, perform less effectively. Since CAM is designed for additive noise models, it exhibits inefficacious perfor- mance on this dataset. A majority of methods featuring regressions in the learning process can adequately predict causal relations of the CE-Net benchmark. Hence, IGCI is the only baseline tested that underperforms on this bench- mark. Due to the difficulty of CE-Cha, the results on this set resemble those on SIM, where most methods struggle, except for GPLVM.
Impact Statement
This paper introduces a novel method whose goal is to ad- vance the field of causal discovery, with potential applica- tions across various scientific disciplines. While our work potentially carries broad societal implications, we do not believe that there are any specific consequences that require partieular emphasis here.
Performance on Real-World Benchmark On the real- world Tübingen benchmark, complexity-based methods, including our COMIC approach, GPLVM, SLOPPY with AIC, SLOPE, and QCCD, achieve comparable accuracy and Bi-AUROC scores. This highlights the effectiveness of compression-based methods, whose Bi-AUROC scores are notably higher than the remaining maximum likelihood- based methods that focus solely model fitness. As noted by Marx & Vreeken (2019a), SLOPPY with AIC offers greater flexibility for more complex datasets, such as this one, and performs better than the BIC variant. A simi- lar pattern appears in IGCI, where the uniform reference measure variant obtains higher accuracy and Bi-AUROC. Although SLOPER outperforms its respective variant SLOPE on previous benchmarks, it is surpassed by SLOPE on this benchmark.
References
Blei, D. M. and Jordan, M. I. Variational inference for Dirichlet process mixtures. Bayesian Analysis, 1(1):121- 143, 2006.
Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. Varia- tional inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859-89)) 877, 2017.
Blier, L. and Ollivier, Y. The description length of deep learning models. In Advances in Neural Information Processing Systems (NeurIPS), volume 31, 2018
6. Conclusion
In this work, we have proposed COMIC-a neural network compression-based approach for determining the cause and effect via the variational Bayesian code-where more uni- versal and scalable class of neural networks is utilized for modeling the conditionals to improve fitness and induce better codelengths. With the variational Bayesian coding scheme, the algorithmic complexity of these networks can be assessed empirically to approximate the theoretical Kol- mogorov complexity, and its identifiability can also be scru- tinized from a marginal likelihood- based perspective. The effectiveness of COMIC has been validated through com- prehensive experiments, delivering promising results and demonstrating enhanced performance compared to most related methods based on the compression and FCM regres- sion objectives across multiple benchmarks.
Blöbaum, P., Janzing, D., Washio, T., Shimizu, S., and Schölkopf, B. Cause-effect inference by comparing re- gression errors. In Proceedings of the International Con- ference on Artificial Intelligence and Statistics (AISTATS), volume 84 of Proceedings of Machine Learning Research, pp. 900-909,10 2018
Bornschein, J., Chiappa. S., Malek, A., and Ke, R. N. Pre- quential MDL for causal structure learning with neural networks. In Proceedings of the Conference on Neural Informational Processing Systems (NeurIPS) WHY-21 Workshop, 2021.
Bühlmann, P., Peters, J., and Ernest, J. CAM: Causal addi-
tive models, digh-dimensiooal order search and penalized
9

---

Identifying Causal Direction via Variational Bayesian Compression
regression. The Annals of Statistics, 42(6))2526-2556, 2014.
Duong, B. and Nguyen, T. Heteroscedastic causal structure learning. In Proceedings of the European Conference on Artificial Intelligence (ECAT), volume 372 of Frontiers in Artificial Intelligence and Applications, pp. 598-605, 2023.
Chickering. D. M. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3 (Nov):507-2) 5544 2002.
Goudet, O., Kalainathan, D., Caillou, P., Guyon, I., Lopez- Paz, D., and Sebag. M. Learning functional causal mod- els with generative neural networks. Explainable and Interpretable Models in Computer Vision and Machine Learning, pp. 39-80, 2018.
Coker, B., Bruinsma, W. P., Burt, D. R., Pan, W., and Doshi- Velez, F. Wide mean-field Bayesian neural networks ignore the data. In Proceedings of the International Con- ference on Artificial Intelligence and Statistics (AISTATS), volume 151 of Proceedings of Machine Learning Re- search, pp. 5276-5333, 2022.
Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. Measuring statistical dependence with Hilbert-Schmidt norms. In Proceedings of the International Conference on Algorithmic Learning Theory (ALT), pp. 63-77, 2005.
Cundy, C., Grover, A., and Ermon, S. BCD Nets: Scal- able variational approaches for Bayesian causal discovery. In Advances in Neural Information Processing Systems (NeurIPS), volume 34, pp. 7095-7110, 2021.
Grünwald, P. D. The minimum description length principle. Adaptive Computation and Machine Learning. The MIT Press, 2007.
Daniušis, P., Janzing, D., Mooij. J., Zscheischler, J., Steudel, B., Zhang, K., and Schölkopf, B. Inferring deterministic causal relations. In Proceedings of the Conference on Uncertainty in Artificial Intelligence (UAI), 2010.
Guyon, I., Statnikov, A., and Batu, B. B. Cause effect pairs in machine learning. Springer, 2019.
Dawid, A. P. Present position and potential developments: Some personal views: Statistical theory: The prequential approach. Journal of the Royal Statistical Sociery. Series (General), 147(2):278 292. 1984.
Hinton, G. E. and van Camp, D. Keeping the neural net- works simple by minimizing the description length of the weights. In Proceedings of the Annual Conference on Computational Learning Theory (COLT), pp. 5-13, 1993.
Honkela, A. and Valpola, H. Variational learning and bits- back coding: An infformation-theorett view to Bayesian learning. IEEE Transactions on Neural Networks, 15(4): 800-810, 2004.
Deleu, T., Góis, A., Emezue, C., Rankawat, M., Lacoste- Julien, S., Bauer, S., and Bengio, Y. Bayesian structure learning with generative flow networks. In Proceedings of the Conference on Uncertainty in Artificial Intelligence (UAI), volume 180 of Proceedings of Machine Learning Research, pp. 518-528, 2022.
Horiik, K., Stinchcombe, M., and White, H. Multilayer feedforward networks are universal approximators. Neu- ral Networks, 2(5):359- 366, 1989.
Dhir, A., Power, S., and der Wilk, M. Bivariate causal discovery using Bayesian model selection. In Proceed- ings of the International Conference on Machine Learn- ing (ICML), volume 235 of Proceedings of Machine Learning Research, pp. 10710-10735, 2024a.
Hoyer, P., Janzing, D., Mooij, J. M., Peters, J., and Schölkopf, B. Nonlinear causal discovery with addi- tive noise models. In Advances Neural Information Processing Systems (NIPS), volume 21, 2008
Dhir, A.. Sedgwick, R., Kori, A., Glocker, B., and van der Wilk, M. Continuous Bayesian model selec- tion for multivariate causal discovery. arXiv preprint mrXiv:2411.0015: 2024b.
Immer, A., Schultheiss, C., Vogt, J. E., Schölkopf, B., Bühlmann, P., and Marx, A. On the identifiability and estimation causal location-scale noise models. In Pro- ceedings of the International Conference on Machine Learning (ICML), volume 202 of Proceedings of Ma- chine Learning Research, pp. 14316-14332, 2023.
Dhir, A., Ashman, M., Requeima, J., and van der Wilk, M. A meta-learning approach to Bayesian causal discovery. In Proceedings of the International Conference on Learning Representations (ICLR), 2025.
Janzing, D. and Schölkopf, B. Causal inference using the algorithmic Markov condition. IEEE Trunsactions on Information Theory, 56(10):5188-5144 2010.
Duong, B. and Nguyen, T. Bivariate causal discovery via conditional divergence. In Proceedings of the Conference on Causal Learning and Reasoning (CLeaR), volume 177 of Proceedings of Machine Learning Research, pp. 336-252, 2022.
Janzing, D., Mooij, J., Zhang, K., Lemeire, J., Zscheis- chler, J., Daniušis, P., Steudel, B., and Schölkopf, B. Information infornation-geomeetic approach to inferring causal direc- tions. Artificial Intelligence, 182-183:1-33 2012.
10

---

Identifying Causal Direction via Variational Bayesian Compression
Kalainathan, D., Goudet, O., and Dutta, R. Causal discovery toolbox: Uncovering causal relationships in Python. Jour- nal of Machine Learning Research, 21(37)11-5, 2020.
Marx, A. and Vreeken, J. Telling cause from effect by local and global regression. Knowledge and Information Systems, 60(3):1277-1305. 2019b.
Kaltenpoth, D. Don't Confound Yourself: Causality from Biased Data, PhD thesis, Saarland University, 2024.
Marx, A. and Vreeken, Formally justifying MDL-based inference of cause and effect. In Proceedings of the AAAI Workshop on Information- anformation-Theorrticc Causal Inference and Discovery (ITCI'22), 2022.
Kaltenpoth, D. and Vreeken, J. Causal discovery with hid- den confounders using the algorithmic Markov condition. In Proceedings of the Uncertainty in Artificial Intelli- gence (UAI), volume 216 of Proceedings of Machine Learning Research, pp. 1016-1026, 2023.
Mian, O. A., Marx, A., and Vreeken, J. Discovering fully oriented causal networks. In Proceedings of the AAAI Conference on Artificial Intelligence (AAAT), volume 35, 2021.
pp.
Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. In Proceedings of the International Confer- ence on Learning Representations (ICLR), 2015.
Mooij, J. M., Stegle, O., Janzing, D., Zhang, K., and Schölkopf. Probabilistic latent variable models for distin- guishing between cause and effect. In Advances in Neural Information Processing Systems (NeurIPS), volume 23, 2010.
Kolmogorov, A. N. On tables of random numbers. Sankhyā: The Indian Journal of Statistics, Series A (1961-2002), (25):369-376, 1963.
Mooij, J. M., Peters, J., Janzing, D., Zscheischler, J., and Schölkopf, B. Distinguishing cause from effect using observational data: Methods and benchmarks. Journal of Machine Learning Research, 17(32)11-12,, 2016
Li, M. and Vitányi, P. An Introduction to Kolmogorov Com- plexity and Its Applications. Texts in Computer Science. Springer, 4th edition, 2019.
Lin, Y., Huang, Y., Liu, W., Deng, H., Ng. I., Zhang, K., Gong, M., Ma, Y.AA., and Huang, B. A skewness-ben -ased criterion for addressing heteroscedastic noise in causal discovery. Proceedings of the International Conference on Learning Reprexentations (ICLR), 2025.
Neal, R. M. Priors for infinite networks. In Bayesian Learning for Neural Networks, volume 118 of Lecture Notes in Statistics, chapter 2, pp. 29-53. Springer, 1996.
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A.,, Yang. E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. PyTorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems (NeurIPS), volume 32, 2019.
Lorch, L.,, Rothfuss, J.. Schölkopf, B., and Krause, A. DIBS: Differentiable Bayesian structure learning. In Advances in Neural Information Processing Systems (NeurlPS), volume 34, pp. 14111-24123, 2021.
Lorch, L., Sussex, S., Rothfuss, J., Krause, A., and Schölkopf, B. Amortized inference for causal structure learning. In Advances in Neural Information Processing Systems (NeurIPS), volume 35, pp. 13104-13118, 2022.
Pearl, J. Causality. Cambridge University Press, 2009.
Peters, J., Mooij, M., Janzing, D., and Schölkopf, B. Causal discovery with continuous additive noise mod- els. Journal of Machine Learning Research, 15(8))22009- 2053, 2014.
Louizos, C., Ullrich, K., and Welling, M. Bayesian compres- sion for deep learning. In Advances in Neural Information Processing Systems (NeurIPS), volume 30, 2017.
MacKay, D. J. C. Comparison of approximate methods for handling hyperparameters. Neural Computation, 11(5): 1035-1068, 1999.
Peters, J., Janzing, D., and Schölkopf, B. Elements of Causal Inferecee: Foundations and Learning Algorithmss. Adap- tive Computation and Machine Learning. The MIT press, 2017.
Marx, A. and Vreeken, J. Telling cause from effect using MDL-based local and global regression. In Proceedings of the IEFE International Conference on Data Mining (ICDM), pp. 307-316, 2017.
Reisach, A., Seiler, C., and Weichwald, S. Beware of the simulated DAG! Causal discovery benchmarks may be easy to game. In Advances in Neural Information Pro- cessing Systems (NeurIPS), volume 34, pp. 27772-2774, 2021.
Marx, A. and Vreeken, J. Identifiability of cause and ef- feet using regularized regression. In Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD), pp. 852- 2019a.
Rissanen, J. Modeling by shortest data description. Auto- matica, 14(5):465-477 471, 1978.

---

Identifying Causal Direction via Variational Bayesian Compression
Shimizu, S., Hoyer, P. O., Hyvärinen, A., Kerminen, A., and Jordan, M. A linear Gan-Gansiian acyclic model for causal discovery. Journal of Machine Learning Research, 7(72), 2006.
Williams, C. Computing with infinite networks. In Ad- vances in Neural Information Processing Systems (NIPS), volume 9, 1996.
Zhang, and Hyvärinen, On the identifiability of the post-nonlinear causal model. In Proceedings of the Con- ference Uncertainty in Artificial Intelligence (UAI), pp- 647-655, 2009.
Solomonoff, R. A formal theory of inductive inference. Part I. Information and Control, 7(1):1-22, 1964a.
Solomonoff, J. A formal theory inductive inference. Part II. Information and Control, 7(2):224-254, 1964b.
Sønderby, C. K., Raiko, T., Maaløe, I.., Sønderby, S. K., and Winther, O. Ladder variational autoencoders. In Ad- vances in Neural Information Processing Systems (NIPS), volume 29, 2016.
Tagasovska, N., Chaoez-Demoulin, V., and Vatter, T. Dis- tinguishing cause from effect using quantiles: Bivariate quantile causal discovery. In Proceedings of the Interna- tional Conference on Machine Learning (ICML), volume 119 of Proceedings of Machine Learning Research, pp. 9311-9323, 2020.
Titsias, M. and Lawrence, N. D. Bayesian Gaussian process latent variable model. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AIS- TATS), volume 9 of Proceedings of Machine Learning Research, pp. 844 851, 2010.
Tran, Q.-D,, Nguyen, P., Duong, B., and Nguyen, T. Dif- ferentiable Bayesian structure learning with acyclicity assurance. In Proceedings of the IEEE International Conference Data Mining (ICDM), pp. 598-607, 2023.
Tran, Q.-D., Duong, B., Nguyen, P., and Nguyen, T. Robust estimation of causal heteroscedastie noise models. Proceedings of the SIAM International Conference on Data Mining (SDM), pp. 788-796. 2024a.
Tran, Q-D., Nguyen, P., Duong, B., and Nguyen, T. Con- straining acyclicity of diffcrentiable Bayesian structure learning with topological ordering. Knowledge and Infor- mation Systems, 66(9):5605-560 5630, 2024b.
Voita, E. and Titov, I. Information-thhereti probing with minimum description length. In Proceedings of the Con- ference on Empirical Methods Natural Language Pro- cessing (EMNLP), pp. 183-196, 2020.
Wallace, C. S. Classification by minimum--messge--eee inference. In Advances in Computing and Information Proceedings of the International Conference on Com- puting and Information (ICCI), volume 468 of Lecture Notes in Computer Science, pp. 72-81, 1990.
Wallace, C. S. and Freeman, P. R. Estimation and infer- ence by compact coding. Journal of the Royal Statistical Society. Series (Methodological), 49(3):240-265 1987.
12

---

Identifying Causal Direction via Variational Bayesian Compression
A. Kolmogorov Complexity, Algorithmic Mutual Information, & Minimumm Description Length
Kolmogorov Complexity In this section, we will introduce some related definitions and notations regarding the algorithmic complexity, which is also known as the Kolmogorov complexity. This complexity represents the length of the ultimate lossless compression of the data (Marx & Vreeken, 2019b).
Definition A.1 (Kolmogorov Complexity, Kolmogorov, 1963; Li & Vitányi, 2019). With a universal Turing machine U4 for every program p {0, 1)" that generates æ and halts from an input y, the Kolmogorov complexity the length of the shortest program. Formally, the Kolmogorov complexity is defined as
K (x y) min (length (p) 1 (w.p)) =]-
(20)
When the input y is an empty string e, we have K (I) K (æ e).
In various causal discovery literature, it is necessary to evaluate the complexity of a distribution with a probability density function. As a consequence, we consider the Kolmogorov complexity of a function (x).
Definition A.2 (Kolmogorov Complexity of a Function, Li & Vitányi, 2019). To compute the complexity of a function f (x), we need to find the shortest program that generates o (a:) from x with precision g as follows
K (x)) = min Iength (p) I (fa, (a,)))) - (x)1 < #
(21)
Algorithmic Mutual Information From the Kolmogorov complexity, the algorithmic mutual information between two binary strings is also defined quantify the amount of overlapping information.
Definition A.3 (Algorithmic Mutual Information, Li & Vitányi, 2019). The algorithmic mutual information between two binary strings æ and y
LA (x: v) = K (v) K (y x
K (x) + K (v) K (x,
(22)
where x² is the shortest description of the string x and denotes equality up to an additive constant, which only depends on U and is independent of x and y-
Minimum Description Length in Causal Discovery Literature In the context of causal discovery task, there are two major hindrances (Kaltenpoth & Vreeken, 2023) when estimating the Kolmogorov complexity K (p (X)) of a distribution (X) from the data samples x: (I) there is no knowledge about the true underlying distribution (X), and (2) the Kolmogorov complexity is incomputable. The former problem can be resolved by estimating the model through the joint complexity K (xx, p (X)), which on expectation over p (X) will yield K (p (X)) + (p (X)) up to an additive constant (Marx & Vreeken, 2022). The latter problem requires an approximating codelength L (z. p (X)) that mirrors K (x, p (X)), commonly selected according the minimum description length (MDL, Grünwald, 2007) principles.
The MDL codelengths are computed by limiting the finding of the shortest program/model that generate the data and halt from a predefined set instead of the universal Solomonoff prior (Grünwald, 2007: Li & Vitányi, 2019; Solomonoff, 1964a;b), By employing sufficiently broad class of models, designed to maintain independence from their conditioning variables, and using a large enough number of samples, the gap between the approximated codelength and the Kolmogorov complexity can be expected to decrease (Kaltenpoth, 2024; Marx & Vreeken, 2022). While approximation gaps may still exist in our method, they should not affect the identifiability result, which is assessed through marginal likelihoods and remains independent of the Kolmogorov complexity-bee postulates in Sec. 3.
B. Prior and Variational Distributions
Let us consider fully-connected layer with the weight matrix W = [Wey]] and the bias vector b = [by]. The parameters W and b in each layer form the set of parameters in Sec. 4.
Priors over Parameters As the parameters are continuous, we consider Gaussian priors with zero means and variances as hyperparameters. For simplicity, we utilize the standard Gaussian distribution for the biases. The priors of the weights and
13

---

Identifying Causal Direction via Variational Bayesian Compression
biases be formulated
P (W) = II (w)), p(wsy) = N (wyy) |0, (0. 1
(23)
p (b) = I p (b,).
p(b,) N (by 10,0)
(24)
Hyperparameters Selection Let z = [zij, Øb] denote the scale hyperparameters of the priors as in Eq. (23) and (24). The marginal likelihood can be acquired by marginalizing over z as
log (y(:)) (x(l:))) log / p (yJiN) (,(:N), 2)
,2) p (2) dz,
(25)
where p (z) is the prior over the hyperparameters z. Following previous work by Dhir et al. (2024a), the integral in the marginal likelihood above is approximated with Laplace approximation (LA, MacKay, 1999). In this approximation method, we assume a Gaussian distribution around the maximum a posteriori (MAP) solution and consider the marginal likelihood as the normalizing constant of this Gaussian distribution. The formulation of the normalizing constant computed around the MAP solution in LA is as follows
log (D(LN)) lollar)) so log (U(LN)) ( 2)) (2)
A2
(26)
2T
where 2 = arg max loogp (y(LIN) 2 ==(aa)))] and
(27)
x
Ax == - V2 logo (y(:N)) 2 allla)
(28)
We also follow Dhir et al. (2024a) by choosing uniform hyperprior over the hyperparameters z and assuming a sufficient amount of samples so that the distribution concentrates around a single point. Given these conditions, we can simply approximate log p (y(1: (.(::))) ~ log (y(2: (1:1.) 2). From an implementation perspective, this means that we can regard the priors z as additional parameters to be optimized in conjunction with the following variational parameters.
Variational Posteriors over Parameters To evaluate the posteriors of W and b, we employ mean-field variational inference (MF-VI, Blei et al., 2017) where the posterior of each weight and bias is assumed to be independent. The posterior distribution of each weight and bias is modeled Gaussian distribution, given by
ls (W)=
II 4s (wx)), gø (wij) N (wiy jw.ej, w.ss))
(29)
ao (b) II as (6,),
as (b,) N (by |b.,, ob.))
(30)
Forward Pass With Gaussian Variational Posteriors The forward pass of a fully-connected Bayesian layer with the Gaussian variational posteriors is presented in Alg. 1. H. E NxA where N is the number of samples and A is the input dimension. Let Mw = [pw, and Vw = (ow.;) represent the means and variances of the weight matrix W,
Pb = [mn_s] and oî [b]]] be the means and variances of the bias vector b. @ be the Hadamard (element-wise product, and a2 and o respectively denote the Hadamard dieementwwise square and square root. Corresponding to this algorithm, the output matrix will include N samples of B nodes. Given that the cost of sampling one instance from the standard Gaussian distribution N (0, 1) is O (1), the computational complexity of forwarding through one layer using local reparametrization is O (NA
C. Causal Identification via Bayesian Model Selection
C.1. Bayesian Model Selection Criterion
Let X and Y not be independent of each other, and assume there are no hidden confounders between these two random variables (i.e., the assumption of causal sufficiency). Dhir al. (2024a) have introduced Bayesian framework for inferring
14

---

Identifying Causal Direction via Variational Bayesian Compression
Algorithm 1 Forward pass of fully-connected Bayesian layer with Gaussian variational posteriors
Require: input: Hi RRNAA; parameters: Mw RRABB Vw RAXB Mb RR o R
1: Mout < My + Mb
2: Vout < V + tog
3: E ~ N 1) {dims: N x B}
4: Hour f Mom + Vol 2 E
5: return Hout
the causal direction between two possible choices Mx-x and My from an observational dataset DV with N samples via the posterior
p (Mi DV) p (DN M.) (M)
(31)
(DN)
where Mj E M = [Mx and p (D*) = EMMEMP (DN | M.)) p (M.). The likelihood p (D M.) has a prior p (T] M,) and can be retrieved by marginalizing p (DN ग. M&) over т. Hence, this likelihood is also referred as the "marginal likelihood".. As we do not possess any knowledge about each causal direction, the prior probabilities are set uniformly to p (M,) = [M] = 0.5.
With this prior, the log-ratio between two posteriors of the causal directions are then computed balance the evidence p (D) and achieve the following equations
p (M DN
p (DN Mx-x) (M--x))>(XX>))) (Mx+r)
log (M DM)
log
(32)
p (DN D
My (My
(DN
Mx-v).
log
(as (Mx-+y) = (My- = 0.5)
(33)
p (DN My+x)
Eq. (33) shows that choosing a model M, so that log p (M, DN) > logp (M; DM) is equivalent to log (DN N M.) > log (DN) | MM),, where M, is the inverse direction of Mi, therefore, we can select the causal direction depending on the log-ratio of the marginal likelihoods as follows
Mx-.Y if logp (DN
x log (DN My > 0,
M =
MY indecisive
if log p D Mx-x) - -log (DN
My < 0,
(34)
if log (DN Mx-r)= - llgg (DN My = 0.
C.2. Causal Identifiability of Marginal Likelihoods
From Eq. (34), Dhir et al. (2024a) study the causal identifiability based on the separability of the two marginal likelihoods of the data given the models. First, the distribution-equivalenn causal models and Bayesian distribution-equivalent causal models are defined, which are the cases where the causal direction can not be determined.
Definition C.1 (Distribution--Euuivalleee of Causal Models, restated from Def. 2.2 of Dhir et al., 2024a). Let two causal models Mx-Y and My respectively have (mx,cyix) E Cx x Cylx and (my, Cy denote the sets of marginal distributions of X and Y, and Cyix and Cxly denote the sets of conditional distributions of given and X given Y. If Mx-y and My are distribution-equivalee
E Cy x Cxjy. where Cx and
there exists a unique bijective map Y Cx x Cy1x + Cy x Cxjy such that for every (mxxey|x) E Cx x Cyjx, there holds an equality of joint likelihoods
mx (z) myjx (y x) my (y) cxjy (x v) Vx. V.
(35)
where (my,exiv) (mxxcypx).
If the two causal models are distribbttion-equivllee the maximum likelihoods cannot distinguish the causal direction. Otherwise, there exists some dataset D (N sufficiently large) such that the maximum likelihood can determine the cause direction. For Bayesian causal models, which include auxiliary prior distributions over the marginal distributions (Cx and Cy) and the conditional distributions (Cyjx and Cxly)), we have a definition of Bayesian distribution -equvvalenn causal models as follows:
15

---

Identifying Causal Direction via Variational Bayesian Compression
Definition C.2 (Bayesian Distributtion--quivall
of Causal Models, restated from Def. 4.4 of Dhir et al., 2024a). Given two Bayesian causal models Mx-+ and My-x, they are Bayesian distrlbution-equiialeet if for all N E N and for all DV, the marginal likelihoods of the data given the models are equal, i... p (D Mx-r) p (DN My-x).
If two possible causal models Mx-y and -y are Bayesian distribution- -euibuuiiiin-uivvenntteee according to this definition, regardless of the observational dataset D> the log-ratio of marginal likelihoods cannot be utilized as criterion for distinguishing the causal direction (third case in Eq. (34)). Bayesian iettribution-qqivvllnnt Ietributiionequivaleent since it not only implies the equality of maximum likelihoods but also imposes the equality of marginal likelihoods, which is the expectations of the likelihoods over their prior distributions. Hence, if two causal models are not distribution-equiaalenn (i.e.. identifiable via maximum likelihoods), they are not Bayesian distribution-equivaleee (i.e., identifiable via marginal likelihoods). In cases of distribution equivvllnt causal models, these models need a necessary condition of epparabl-commattbllii for them to be Bayesian distribution -atattioon-ivvvleeeeeee which is also proposed by Dhir et al. (2024a) determine the set of distributions C that can be utilized for modeling the marginal likelihoods. The definition of this condition is as follows:
is a more specific case of
Definition C.3 (Separable- (Separabee-oompatbblle Compatibility of Bayesian Causal Models, restated from Def. 4.6 of Dhir et al.. 2024a). Let two Bayesian causal models Mx-+y and My-x include (rix) E Cx Crix and (my,,xir) E Cy x Cxpy with their corresponding priors being Tx-xY xxIxx = ax (mx) πyix (evix) and Ty- (my,Cxir) = ay (my) x (exly), respectively. Let us denote as in Def. C.1, if the push-forward x is separable with respect to Cy x Cxly, i.e., -y- (7 (mx, mylx)) )) Ky (my) Txly (ex(r), and ny my is separable with respect to Cx x Cy then the two causal models are separable-compatin
This definition is similar to the principle of independent causal mechanisms where the joint density can only be factorized (separable) with respect to the causal direction and the anti-causal factorization do not satisfy separability. This means that even the prior of the anti-causal models can be represented the product of two priors of distributions, these two priors are not independent of each other. As the separable-compatiiiii is the necessary condition of the Bayesian distribution- equivalence, if the two causal models are not separable-compatee they are not Bayesian cestriuttinneqquialennr implying that they are identifiable via the marginal likelihoods.
C.3. Analysis of Our Models (Proof of Prop. 4.2, Sec. 4.4)
To validate the identifiability of the causal models of COMIC, we need to analyze their separable-compptiillty Let us recall the non-separable--compiiiii of our proposed causal models in Prop. 4.2 as follows:
of Our Causal Models). Let the two Bayesian causal models Mx-+ and My-x respectively be factorized into the marginal densities-pi (x Mx-Y) and (y My- the conditional densities-p(y (y| z. , Mx-x) and (x y, My My)x).. The marginal densities are standard Gaussian N 0, 1). The conditional densities are Gaussian likelihoods N (.1 (- (=; (²;;;);)))))))))) ))) where (:; 0) and a (-:0) are parametrized by neural networks as described in Eq. (14), Sec. 4.3, and the parameters O follow Gaussian priors specified in App. B.
Proposition 4.2
Celeparable-Coppatiile
Under these conditions, the Bayesian causal models are not separuble-compuibble
The proof of this pon-sepprrablc-cmmppbbiiiii is provided below. For ease of proof, neural networks with one hidden layer can be considered as Gaussian processes under Central Limit Theorem when the number of hidden units is large enough (Neal, 1996; Williams, 1996).
Additive Noise Models (ANMs) First, we study the case of an ANM where the mean of y is modeled by a neural network fy R + IR with one hidden layer. With an input x, the hidden layer is computed as follows:
by = by (x) = a (wITy x + bby),
(36)
where Why = [wny..n 1 E RRIDD, Why is sampled from a Gaussian prior A (0, shy. bby = [bn E R RR,,,
bhy-in in ~ N (0, o and a (.) is an activation function such as the hyperbolic tangent (tanh). From the features by in the hidden layer, the output function fy is then linearly represented as
fr (x) = w J. h (z) + byy.
(37)
16

---

Identifying Causal Direction via Variational Bayesian Compression
where Wy (t E Dbx1, Ws is also sampled from a Gaussian distribution N (0, 10.27r.0 2)rren) and bsr N (0.03,,) ). Let us denote that fr (W/r, byy Ony (Way bby ). 2/y iiviin::77 o/v] and zn. the Gaussian process of IY is GP (my (ar), K (2r, a')) where
2
by in hy
myy (x) Eo,, Ony [fr (x)] 0,
(38)
and
Ky (2,z') = Eo Oby
[fr (x) fy (z')]
Eosy Ony
[(w]. J. by + byy) (w) [w byy b (by= - (x¹))
Easy Bby
[[W], hy) (w fr h(y) -w wirhrbsr + byv w hy + b
Er Pay
( T fr hr) (w]y h1y)] + Ensy [
]]
EE Oay
Σ W W hy Σ Sr So hyy
+ Bosr [6,r]
Dn Dr.
EE ty Ony
Wr fr on My Ny
+ Ebrry
[b,y
in
Jn
Dr
Ew Ony
Σ
w Jy An by Nlysen
+ Ensy [67r]
in
Dn
EBny
Σ Eorr by, Hysiu
+ Bbbr 10
BBuy
[hy diag (zit)) hly] + obr
2 by (x) diag (2?,) hy (x') + +obyy
Bouy
(39)
With fr (x) + ey and ey ~ N (0, o2 o ), the generative processes are
fy I, Zy any (0, (x,x)),
(40)
y v fy N of
(41)
y x 2, 2 ~ N (0, (0,Kiy (2,z) + o²
(42)
For the causal model Mx-Y,, the marginal distribution p (x) and the conditional distribution p (y ( with respect to the Gaussian process above are as follows
(x Mx-y)==N = | ( x, 2y 2hy Mx-Y)= = (y 0, Ky (I, +
(43)
(44)
where 2 and 2iy are hyperparameters. The selection method for these hyperparameters is discussed in App. B. The corresponding distributions for the causal model My are
p My+ (1 00,1),
(45)
p (x y, 2 2
N (s 0, (v.y) + +0
(46)
To analyze the separable-compatinnn we need to compare the causal factorization of Mx-,Y with the anti-causal factorization of My-x, which is formulated as
/ (x v, 2 tix My pvv My dy
p (x 2)x 2ix My
(47)
17

---

Identifying Causal Direction via Variational Bayesian Compression
p (y x, 23x Tix My-x) (
y, fx 2nx My
p (y
(48)
p x "/ nxx My++x)
As f (y.y) y) is non-linearly dependent on y, the marginal distribution of 2²xx Thx My is non-Gaussian, which makes it not belong the same class of distribution p (z: Mx-+)). Moreover, because the Gaussian distribution is not a conjugate prior if the likelihood is a aariance-parametrizes Gaussian distribution, the posterior distribution 2, 2]x hx -y will not be Gaussian. Hence, these distributions are not in the same classes of distributions as those of Mx-Y. In addition, let us examine the variance x with respect to the anti-causal marginal
[x-]] / N 11110 (f (22N(e (xr 0, Fxx (v.u) + o²xx dx dy
E "x(leex My 1x)
(49)
= / No (K (v.v) + oxxx) dy-
(50)
Although this integral does not have a closed form, is to be expected that the variance of xt x)x Thx -y is dependent on 2), and L,x. which should also be the case for the variance of y x. 22,x 2 My Hence, the distributions in the anti-causal factorization share identical priors, which make them not independent of each other and violate the principle of independent causal mechanisms (ICMs, Peters et al., 2017, Sec. 2.1). From these aspects, we can accept that the models Mx-y and - are not separable-compatible.
Location-Scale Noise Models (LSNMs) In LSNM, not only the mean but also the scale of the noise is parametrized. As mentioned in Eq. (14), Sec. 4.3, we consider the case where y is generated from x via a neural network fy : IR - R2 with two output nodes fY,1 and fr.2. The generating process of can be presented as follows
y fy N (fr.1, (frz)).
(51)
where C (-) is a positive link function, such as exponential or softplus function. In our work, as fy, and fy. are two output nodes of the same network, the Gaussian processes of fy,1 and fy. from x can be described as
fY, x, F. , zhy ~ ~N (0, Krr, (x,x))
2/y,
(52)
where K Kr, (x, x²) = Eeny by (z)
T diag (2)x3)) hy (2') + o and i E {1,2}.
The conditional distributions p y x, 27, fr. 2/y2' 2/y Mx Mx- and x y. fx. 2 fx.=' ; will be more complicated in this setting, however, we can expect similar ceparable-compatibiie as in the location-only setting above.
Since the property of separable-compattiiity does not hold in our causal models, following the definition of separable- compatibility in Def. C.3 and Bayesian serriiutiooneuuvvlleneee in Def. C.2, our causal models are not distribution-equivalleee which indicates that there exists dataset where the causal direction is identifiable via marginal likelihoods. As result, we propose Cor. 4.3, which is recalled as follows
Corollary 4.3 (Causal Identifiability of Our Causal Models). Let ( Mx-xy) and (- My-x) respectively be the marginal likelihoods of the data given the causal models Mx-Y and My. Because of the in Prop. 4.2, there exists dataset D whose causal direction is identifiable via the difference between Bayesian indicator scores:
Bayes (DN)
Bayes (DN) = A -
- 4 X- (DN).
Bayes
(19)
D. Implementations & Hyperparameters
All the experiments are conducted on a workstation with an Intel Core i7 processor, 64 GB of memory, and 3 TB of storage, except for those involving GPLVM, which are executed on NVIDIA Tesla V100 GPUs.
D.1. COMIC
Our proposed Caaohod-COMIC--C implemented in Python with the PyTorch library (Paszke et al., 2019). The implemen- lation of the fully-connected Bayesian layers with Gaussian priors is adapted from the publicly available code of Louizos
Bayesian layers: https: //gi thub. com http:///weeeeaaaaooooooaaaaaaaaa
18

---

Identifying Causal Direction via Variational Bayesian Compression
al. (2017). As mentioned in App. C, Bayesian neural networks are employed to model the mean and the standard deviation of each assumed direction, i.e., X Y and Y + X. Each neural network includes one hidden layer with Dn 50 nodes with the hyperbolic tangent (tanh) as the activation function and fully-connected output layer. The natural exponential function is chosen as the positive link function for ensuring the positivity of the predicted standard deviation values. The neural networks of COMIC are optimized to minimize the variational Bayesian objective with the Adam optimizer (Kingma & Ba, 2015) and the cosine learning rate scheduler with a maximum learning rate of 10-?, a minimum learning rate of 10-6, and T = 2, 500 training epochs. To avoid the bad local optima of the variational objective, Louizos et al. (2017) suggest the adoption of the "warm-up" scheme from Sønderby et al. (2016), where the model complexity term KL (ge (w) II pw (w)) is weighted by a hyperparameter |3 annealed linearly from 0 to 1 for the first Twu training epochs. In our implementation, we choose Twu- = 250. In addition, to avoid overfitting the hyperparameters of the priors, we pretrain the models and the hyperparameters with maximum a posteriori (MAP) estimation before the variational optimization for Twu 2, 500 epochs, where the parameters and hyperparameters are optimized to minimize the negative joint log- logllikelihood of the data and the parameters. All samples of each pair are utilized in the training process, i.e., we set the batch size to be equal to the number of available samples of each dataset.
D.2. Baselines
For information theory-based methods including SLOPPY (Marx & Vreeken, 2019a), SLOPE (Marx & Vreeken, 2017; 2019b), and QCCD6 (Tagasovska et al., 2020), we utilize their original repositories R. The rpy2 interface is employed to incorporate these methods into our source code. As a Python version of IGCI (Daniušis et al., 2010; Janzing al., 2012) is available in the Causal Discovery Toolkit (CDT?, Kalainathan et al., 2020), we employ this version for our experiments. With GPLVM we utilize the Python implementation from the repository of Dhir et al. (2024a). Following to the recommendation of Dhir et al. (2024a), we choose the closed form GPLVM for all datasets except for the Tübingen dataset where the stochastic GPLVM is employed. For regression-based methods, the Python source code of LOCI' is obtainable from the repository of Immer et al. (2023). The R version of CAM (Bühlmann et al., 2014) can be accessed via the Comprehensive R Archive Network (CRAN).
E. Descriptions of the Benchmarks
In the upcoming descriptions, we denote X as the cause, Y as the effect, and Ey as the noise term (Ey 11 X).
AN, AN-s, LS, LS-s, & MN-U (Tagasovska et al., 2020) As their names are the abbreviations of their generating models, the AN and AN-s are simulated from additive noise models Y = fX) + Ey. Correspondingly, the LS and LS-s are generated from the location-scale (heteroscedastic) noise models Y = (X) + (X) x Ey, and the MN-U pairs are sampled with the multiplicative noise models Y == f (X) x Ey as the data generating processes. Injective sigmoid-type functions from Bühlmann et al. (2014) are chosen in contrast to Gaussian processes for the model functions f and g in the MN-U set and the ones with the "-s" suffix to produce more difficult settings. Except for the MN-U benchmark with uniform noises, other datasets are sampled with the Gaussian distribution for Ey. Each dataset contains 100 pairs with 1, 0000 samples for each pair.
SIM, SIM-c, SIM-G, & SIM-In (Mooij et al., 2016) In each dataset, similar to the five benchmarks mentioned above, there are also 100 pairs where each pair has the number of samples of 1, 000. The cause-and-effect samples are generated from a more complicated process as X == fx (Ex) + Mx and Y == fr Ex, + My, where Ex and Ey are sampled from randomly generated distributions P (Ex) and P (Ey). fx and fr are drawn from Gaussian processes, and Mx and My are Gaussian measurement noises. The SIM dataset has this default configuration. In SIM-c, there is a one-dimensional confounder Z included as an additional input of both fx and fy. The measurement noise levels are reduced in the SIM-In pairs to create approximately deterministic relations. The cause X in each SIM-G pair follows distribution that is close
*SLOPe: httpss//eeaa.eggccccpps...pppppp "QCCD: CCCCD//ttpp///tttoooc/aaaaaa//// CCC.://ww/.ttttaooo///.e/oootoss
WGt.W//a.rre/.gaauc.cccaaaaaaaaa aiscovery baypyaaaaannnnn laseesianmmoeeeee selection
"CAM: https////ggttttb.ccc//rraa//////
19

---

Identifying Causal Direction via Variational Bayesian Compression
the Gaussian distribution and the corresponding effect Y iimitate a nonlinear additive noise model with Gaussian noise.
CE-Multi, CE-Net, & CE-Cha (Goudet et al., 2018) CE-Multi encompasses 300 pairs generated from four families of noise models: pre-additive noise model Y := (X + Ey), post-additive noise model Y == f (X) + Ey. pre-multiplicative noise model Y = f (X F and post-multiplicative noise model Y = f (X) Ey, where f can be both linear and nonlinear functions. CE-Net also involves 300 pairs, whose causes are sampled from random distributions and analogous effects synthesized from a neural network with random weights. The CE-Cha set originates from the Chal.ean Cause- Pairs Challenge (Guyon et al., 2019) where the selected 300 pairs only contain either + or Y - X causal relations with X and Y being continuous variables. Each pair of these datasets includes 1, 500 samples.
Tübingen Cause-Effect Pairs (Mooij al., 2016) This real-world benchmark involves 108 cause-effect pairs col- lected from various domains, such as metcorology, biology, economy, engineering, medicine, etc. As most related approaches (Bühlmann al., 2014; Immer et al., 2023; Marx & Vreeken, 2017; 2019a;;; Peters al., 2014) and COMIC are designed for univariate numeric cause-effect pairs, we adopt existing experimental setups and consider 99 pairs that have one-dimensional continuous causes and effects. Every pair in this dataset is weighted to avoid biased results toward related and similar examples of causal relations.
Accessing the Benchmarks
All the benchmarks are available in the repository of LOCI (Immer et al., 2023)!!
F. Detailed Experimental Results
The detailed experimental results in Fig. 1 (Sec. 5) are listed in Tab. 1.
U'See 9
20

---

Identifying Causal Direction via Variational Bayesian Compression
Table Detailed accuracy and bidirectional AUROC scores COMIC approach and baseline methods all henchmarks. Higher accuracy and bidirectional AUROC (Bi-AUROC) are preferable. The highest scores are listed in bold, the second-best overall results are listed in italic. Our COMIC method is compared against SLOPPY (including AIC and BIC variants, Marx Vreeken, 2019a). SLOPER (Marx & Vreeken, 2019h) and SLOPE (Marx & Vreeken, 2017), QCCD (Tagasovska et al., 2020), IGCI (with uniform and Gaussian reference measures, Daniušis et al., 2010), GPLVM (Dhir et al., 2024a), LOCI (Immer et al., 2023), and CAM (Bühlmann al., 2014). Our COMIC method ranks second overall while maintaining decently consistent performance, as indicated by its lower standard deviations compared to most baseline methods.
Method
AN
AN-s
LS
LS-s
MN-U
SIM
SIM-c
SIM-G
SIM-In
CE-Multi
CE-Net
CE-Cha
Tübingen
Overall
Accuracyt
COMIC (Ours)
1.00
1.00
1.00
1.00
1.00
0.57
0.54
0.78
0.89
0.79
0.78
0.43
0.67
0.804 + 0.200
SLOPPY-AIC
0.75
0.31
0.71
0.12
0.04
0.54
0.64
0.52
0.57
0.90
0.58
0.71
0.511 0.249
SLOPPY-BIC
1.00
0.56
0.77
0.45
SLOPER
0.91
0.56
0.80
0.06
0.10
0.57
0.66
0.60
0.81
0.86
0.73
0.59
0.67
0.609 0.260
SLOPE
0.18
0.28
0.20
0.12
0.07
0.46
0.54
0.46
0.47
0.89
0.62
0.56
0.72
0.429 1 0.246
QCCD
1.00
0.85
1.00
0.99
0.99
0.64
0.72
0.66
0.85
0.52
0.82
0.54
0.77
0.796 0.171
IGCI-Uniform
0.20
0.35
0.46
0.34
0.11
0.37
0.45
0.53
0.51
0.92
0.55
0.464
IGCI-Grusiaa
0.89
0.97
0.95
0.94
0.36
0.86
0.68
0.55
0.712 0.210
GPLVM
1.00
1.00
1.00
1.00
1.00
0.83
0.79
0.92
0.90
0.96
0.96
0.75
0.70
0.909 = 0.105
LOCT
1.00
1.00
1.00
1.00
1.00
0.49
0.50
0.78
0.79
0.72
0.77
0.43
0.56
0.773 0.219
CAM
1.00
1.00
0.98
0.52
0.56
0.62
0.82
0.87
0.35
0.78
0.53
0.52
0.725 0.217
Bidirectional AUROC"
COMIC (Ours)
0.59
0.95
0.91
0.47
0.843 0.
SLOPPY-AIC
0.84
0.27
0.80
0.03
0.00
0.60
0.68
0.58
0.67
0.97
0.78
0.62
0.79
0.587 0.301
SLOPPY-BIC
1.00
1.00
1.00
0.58
1.00
0.74
0.74
0,90
0.90
0.59
0.88
0.51
0.59
0.805 = 0.182
SLOPER
0.97
0.60
0.87
0.01
0.03
0.67
0.71
0.68
0.93
0.96
0.84
0.62
0.79
0.668 = 0.313
SLOPE
0.09
0.23
0.12
0.03
0.01
0.48
0.57
0.45
0.46
0.97
0.68
0.59
0.81
0.422 0.307
QCCD
0.91
IGCI-Unifnne Uniform
0.17
0.37
0.51
0.42
0.02
0.38
0.43
0.57
0.51
0.98
0.59
0.58
0.69
0.478 0.235
IGCI-Gmussian
0.98
1.00
0.99
0.99
0.94
0.32
0.38
0.95
0.66
0.78
0.57
0.56
0.64
0,750 0.244
GPLVM
1.00
1.00
1.00
1.00
1.00
0.92
0.91
0.99
0.97
0.98
0.99
0.82
0.80
0.951 + 0.069
LOCI
1.00
1.00
1.00
1.00
1.00
0.53
0.56
0.85
0.91
0.88
0.82
0.45
0.65
0.819 0.202
0.38
0.719 0.240
G. Ablation Studies
Layer Width of Bayesian Neural Networks Coker et al. (2022) have indicated that mean-field variational inference (MF-VI) of wide Bayesian neural networks can cause the network to ignore the data, where the variational posteriors will collapse to the priors instead of the true posteriors. In this part, examine the behavior of our models with different hidden layer widths. The result is depicted in Tab. 2. From these results, there are impacts on performance when the width of the hidden layer increases. When the width reaches 100, the reductions in accuracy scores become noticeable. Nonetheless, the Bi-AUROC scores still remain satisfactory. At 200 hidden nodes, we begin to experience severer declines in performance. However, it is important to note that such excessively large width is not recommended in the bivariate setting, where there is only one input dimension and up to two output dimensions. Thus, our approach will not be substantially affected by the ignorance of the data (Coker al., 2022), if the width of the hidden layer stays within appropriate range (S 100).
21

---

Identifying Causal Direction via Variational Bayesian Compression
Table The effect of hidden layer width the performance of COMIC. Higher accuracy and bidirectional AUROC (Bi-AUROC) preferable. The impact of the width on the Bi-AUROC is tolerable up until 100 nodes. When the width reaches 200 discouraged choice bivariate setting the decline in performance is severe.
Width
AN
AN-s
LS
LS-s
MN-U
SIM
SIM-c
SIM-G
SIM-In
CE-Multi
CE-Net
CE-Cha
Tübingen
Accuracyt
10
1.00
1.00
1.00
1.00
1.00
0.53
0.53
0.83
0.81
0.76
0.77
0.43
0.45
20
1.00
1.00
1.00
1.00
1.00
0.54
0.53
0.81
0.89
0.77
0.77
0.43
0.60
50
1.00
1.00
1.00
1.00
1.00
0.57
0.54
0.78
0.89
0.79
0.78
0.43
0.67
100
0.98
0.94
0.95
0.98
0.96
0.54
0.55
0.76
0.88
0.79
0.72
0.48
0.63
200
0.97
0.91
0.89
0.98
0.87
0.54
0.53
0.74
0.84
0.76
0.69
0.47
0.57
Bidirectional AUROC1
10
1.00
1.00
1.00
1.00
1.00
0.58
0.58
0.90
0.92
0.90
0.83
0.46
0.53
20
1.00
1.00
1.00
1.00
1.00
0.58
0.59
0.89
0.90
0.83
0.46
0.68
50
1.00
1.00
1.00
1.00
1.00
0.58
0.59
0.85
0.95
0.91
0.84
0.47
0.75
100
1.00
0.99
0.99
1.00
0.99
0.62
0.60
0.84
0.95
0.93
0.80
0.49
0.70
200
0.99
0.97
0.97
0.99
0.96
0.62
0.57
0.81
0.93
0.92
0.76
0.48
0.69
Sparsity-Inducing Priors In addition to the MAP approach for selecting the hyperparameters in App. B, the sparsity- inducing prior proposed by Louizos et al. (2017) can be used as an alternative hyperprior. In this prior, the scale hyperparam- eters follow the -nnimmfoorative Jeffreys hyperprior (2) x |=] rather than a uniform hyperprior. With this choice of prior, the posteriors the hyperparameters are inferred using MF-VI with Gaussian distributions adopted as the variational posteriors (Louizos et al., 2017). The Gaussian posteriors of the hyperparameters enable pruning mechanism in which the pruned from the network. We compare the results between these two hyperpriors
weights satisfying log o log 2 in Tab. 3. The results indicate that the uniform hyperprior is generally more effective in comparison to the sparsity-induchie Jeffreys hyperprior most benchmarks, except for SIM-c. Moreover, pruning yields no significant differences in scores compared to the unpruned case of the Jeffreys hyperprior.
Table 3, Performance of COMIC with uniform and Jeffreys hyperpriors (with and without the post-hoc pruning, Louizos et al.. 2017). Higher accuracy and bidirectional AUROC (Bi-AUROC) are preferable. The uniform hyperpriors provide better performance most datasets, except for SIM-c.
Hyperprior
AN
AN.s
LS
L.S.s
MN-U
SIM
SIM-c
SIM-G
STM-In
CE-Multi
CE-Net
CF-Cha
Tübhingnn
Accuracyt
Uniform
1.00
1.00
1.00
1.00
1.00
0.57
0.54
0.78
0.89
0.79
0.78
0.43
0.67
Jeffreys w/o Pruning
1.00
1.00
1.00
1.00
1.00
0.54
0.62
0.66
0.85
0.84
0.70
0.47
0.58
Jeffreys w/ Pruning
1.00
1.00
1.00
1.00
1.00
0.55
0.62
0.67
0.85
0.84
0.70
0.47
0.58
Bidirectional AUROC1
Uniform Jeffreys w/o Pruning
1.00
1.00
1.00
1.00
1.00
0.58
0.59
0.85
0.95
0.91
0.84
0.47
0.75
1.00
1.00
1.00
1.00
1.00
0.58
0.63
0.80
0.93
0.94
0.80
0.47
0.63
Jeffreys w/ Pruning
1.00
1.00
1.00
1.00
1.00
0.57
0.62
0.50
0.92
0.94
0.80
0.47
0.63
Choices for Encoding the Causes While the standard Gaussian distribution is a common choice previous methods for encoding the marginal distributions of the causes, discussed Sec. 4.3, we also investigate more advanced alternative for modeling the marginal distribution of the causes, which is the variational Bayesian Gaussian mixture model (VB- B-MMMM Blei & Jordan, 2006). The VB-GMM is flexible model that can capture multimodal distributions, which may be beneficial in scenarios where the causes exhibit complex distributions. Despite that, the identifiability of the causal direction is more challenging to evaluate theoretically with the VB-GMM. As shown in Tab. 4, VB-GMM substantially improves prediction accuracy and Bi-AUROC over the standard Gaussian on certain synthetic benchmarks. such as SIM, SIM-c, and CE-Net. This improvement, however, not uniform across datasets: on SIM-G and CE-Multi,, performance drops notably, while on
22

---

Identifying Causal Direction via Variational Bayesian Compression
the real-world Tübingen dataset, VB-GMM increases accuracy but slightly reduces Bi-AUROC.
Table The effect of different cause encoding choices on the performance of COMIC, including the standard Gaussian and variational Bayesian Gaussian mixture models (VB-GMM) encoding methods. Higher accuracy and hidirectional AUROC (Bi-AUROC) are preferable. Using the more advanced VB-GMM method for encoding the causes yields mixed results, showing significant improvements in performance on some datasets, but reducing some others, including the real-world Tübingen benchmark.
Cause Encoding
AN
AN-s
LS-s
MN-U
SIM
SIM-e
SIM-G
SIM-In
CE-Multi
CE-Net
CE-Cha
Tübingen
Accuracyt
Standard Gaussian
1.00
1.00
1.00
1.00
1.00
0.57
0.54
0.78
0.89
0.79
0.78
0.43
0.67
VB-GMM
1.00
1.00
1.00
1.00
1.00
0.77
0.71
0.57
0.90
0.74
0.86
0.48
0.73
Bidirectional AUROC1
Standard Gaussian
1.00
1.00
1.00
1.00
1.00
0.58
0.59
0.85
0.95
0.91
0.84
0.47
0.75
VB-GMM
1.00
1.00
1.00
1.00
1.00
0.80
0.79
0.67
0.95
0.81
0.96
0.47
0.73
H. From Bivariate to Multivariate Settings
There are some probable approaches to adapting this work to multivariate settings, which we will explore in this section.
Order-Based Causal Discovery The order-based approach is applied by FCM-based methods (Zhang & Hyvärinen, 2009; Peters et al., 2014; Duong & Nguyen, 2023; Lin et al., 2025), which includes choosing a criterion for discovering the topological order of the variables to construct a corresponding full directed acyclic graph (DAG) and then pruning the excessive edges from this full DAG. For example, the multivariate version of RESIT by Peters et al. (2014) utilizes the HSIC (Gretton al., 2005) score between each variable with the remaining one as a sorting criterion, whose value will be highest in sink node (the last node in topological order). In HOST by Duong & Nguyen (2023), with the assumption of standard Gaussian noise in heteroscedastic noise models (HNMs), the normality of the residuals is evaluated with the Wilk test after regression, which is expected to be the highest in sink nodes. Lin et al. (2025) present a more generalized criterion based on minimal skewness in sink nodes, which enables the identification of HNMs with symmetric noise distributions.
As the joint distribution of sink node X, and its set of non-descendants ND6 (X₃) in the causal graph G can be factorized as y (X,, NDa (X,)) = p (X; NDo (X,)) p (NDa (Xz)) under the principle of independent causal mechanisms, the Kolmogorov complexity can be computed from the Postulate of Algorithmic Independence of Conditionals (Pos. 1, Sec. 3) as K (p (X,, NDG.,)) $ K (p (X, NDG.))) + K (p (NDg.)). We can expect the factorization of the codelengths corresponding to the sink node as the most succinct. As a consequence, we can choose a node yielding the shortest sum of codelengths as the sink node.
Score-Based Causal Discovery MDL-based scores have also appeared in score-based causal discovery literature (Born- schein et al., 2021; Mian et al., 2021) evaluation metrics for candidate causal graphs. GLOBE from Mian et al. (2021) is an extension of SLOPE (Marx & Vreeken, 2017), which employs the algorithm of greedy equivalence search (GES, Chickering, 2002) for searching in the graph space. Bornschein et al. (2021) also propose a prequential MDL--basd score for ranking graphs via their conditional probability distributions of a causal graph. In these approaches, the graphs maximizing the criteria are chosen as the causal graphs. Following this learning scheme, an extension of GPLVM (Dhir et al., 2024a) have also been introduced by Dhir al. (2024b), which adapt the posterior-based model selection criterion in App. C to multivariate formulation.
In addition these point-estimation approaches, Bayesian causal discovery approach can also be utilized infer the full posterior distribution over the causal graphs given the data (G D). This posterior is proportional to the product of the marginal likelihood p (DI G) = I'p (D) G, () p (0) d0 and a prior over the structures (G). The posterior (G DD) can either be estimated in an unsupervised manner, such as those proposed by Cundy et al. (2021); Lorch et al. (2021); Deleu et al. (2022); Tran et al. (2023; 2024b), or learned via supervised training, as demonstrated by Lorch et al. (2022); Dhir et al. (2025). Since the marginal likelihood p (D) G)) is equivalent to the Bayesian codelength, which can be effectively
23

---

Identifying Causal Direction via Variational Bayesian Compression
evaluated via the variational Bayesian framework. Hence, approach is naturally well-suited for this Bayesian approach to multivariate causal discovery.
24


Reviews:
[{"review_id": "RFe8bigT6s", "reviewer": "Reviewer_NJ5G", "summary": "The paper proposes a new method, called ‘COMIC’, for bivariate causal direction identification under the causal sufficiency assumption. It is based on the familiar ICM + MDL principles, utilising a variational Bayesian learning of the complexity of neural network approximations to the marginal/conditionals implied the two different decompositions X-> Y vs. Y-> X. It is intended to offer the same flexibility in model fitness obtained via GPs, but at a significantly lower computational cost / better scalability. The end result is evaluated on a range of benchmark data sets against a range of alternatives, and found to provide near universal improvement over existing methods on nearly all data sets.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank Reviewer NJ5G for recognizing the positive aspects of our work. We expect the following rebuttal will address your concerns.\n\n**1. Underlying model assumptions**\n\nWe summarize three key definitions related to the identifiability of our Bayesian causal models (BCMs), introduced in [1], as follows:\n1. *Distribution-equivalence:* causal model selection via **maximum** likelihood will not be able to identify the causal direction;\n2. *Bayesian distribution-equivalence:* Bayesian causal model selection via **marginal** likelihood will not be able to identify the causal direction, which is a stricter condition than distribution-equivalence;\n3. *Separable-compatibility:* the necessary condition for two distribution-equivalent BCMs being Bayesian distribution-equivalent.\n\nIf two BCMs are distribution-equivalent, as long as they are not separable-compatible, they will not be Bayesian distribution-equivalent [1]. Similar to previous methods based on functional causal models (FCMs), the BCMs as in GPLVM [1] and our work are assumed to be the *underlying generating processes* of the data. With non-separable-compatible BCMs, given enough data, we can estimate the underlying BCMs and identify the causal direction with the marginal likelihoods. The benefit of this approach is that BCMs can be designed to be flexible enough (while keeping non-separable-compatibility) to cover a broad spectrum of generating models, such as GPLVM [1] or Bayesian neural networks in our work. Therefore, your example of normalized linear Gaussian models do not lie into the scope of our BCMs which assume that the effect is generated via nonlinear Bayesian neural networks.\n\n**2. Performance over scalability results**\n\nSince we use the official implementations of the baselines, we do not believe a direct comparison of running time would be fair. Instead, we have opted to evaluate the computational complexity of our forward pass in App. B. However, for a rough comparison, on the AN dataset, our work only took around 3 minutes on a CPU configuration, whereas GPLVM required on average about 2 days 11 hours on 10 GPUs.\n\n**3. Causal sufficiency assumption**\n\nThank you for your suggestion. Since our performance on the SIM-c (with confounders) is promising, relaxing the causal sufficiency assumption will be a potential direction of study for our future work.\n\n[1] Dhir, A., Power, S., & Van Der Wilk, M. (2024). Bivariate Causal Discovery using Bayesian Model Selection. In ICML.", "claims_and_evidence": "Claims are properly supported by theoretical derivations and experimental results.", "methods_and_evaluation_criteria": "Well-known, challenging problem. New approach that seems both efficient and effective, with lots of potential for further extensions. \nMain criticism is that the paper keeps following the pervasive simplifying assumption of ‘no confounding’. Yes, everyone does it, and it makes life so much easier for writing papers, but it is unrealistic and unhelpful in real-world contexts, making the resulting methods, however fast&fancy, often largely useless in practice … maybe something to tackle next?", "theoretical_claims": "Yes, at least to some degree (see Q at 12)", "experimental_designs_or_analyses": "Extensive evaluation (incl. Appendix F). \nOnly thing surprisingly missing from the experimental results is an indication of performance / scalability of COMIC vs. e.g. GPLVM, as that aspect was one of the main motivations behind the current approach.", "relation_to_broader_scientific_literature": "Good: all relevant bases seem to be covered.", "essential_references_not_discussed": "None that come to mind.", "other_strengths_and_weaknesses": "Well written paper with interesting and meaningful contribution to a well-established problem. Fairly technical, but readable. Would have liked to get some more insight/intuition for the resulting computational complexity and stability of the result.\nBut overall solid paper, so for now: clear accept.", "other_comments_or_suggestions": "No.", "questions_for_authors": "Perhaps I am misunderstanding part of Prop.4.2, but it reads as a claim about your (neural network) causal model approximations, stating that they are not separable-compatible, and hence capable of distinguishing between X -> Y vs. Y -> X in the large sample limit, without reference to the actual underlying causal model? But if the underlying model responsible for generating the data is e.g. linear Gaussian (so separable-compatible), then this should be impossible, right? Or is this implicitly excluded in the assumption that the underlying model itself must be a neural network with the  stated parameters? \nTo be clear: I would understand the neural network model as a universal approximation to the real underlying model, and so a Prop4.2/Cor4.3 along the form ‘If the underlying model for X->Y is not distributionally equivalent to a model with Y->X, then our neural network approximation would yield different Bayesian codelengths for the two alternatives, i.e. they would not be separable-compatible.’ \n\nCould you clarify what exactly is meant here?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "JKk52Va1OT", "reviewer": "Reviewer_3KpH", "summary": "This work does not focus on traditional or recently popular functional class-based methods. Instead, it aims to study more general models, where the asymmetry in determining causal direction is assumed based on the Kolmogorov complexity. However, due to the incomputability of Kolmogorov complexity, the Minimum Description Length (MDL) principle serves as a practical proxy. Building on this, the authors explore Bayesian neural networks for causal discovery, leveraging the natural connection between MDL and likelihood. Identifiability analysis for causal direction is provided, and experimental results demonstrate strong performance on both synthetic and real-world benchmark datasets.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer 3KpH for your valuable insights. We provide our responses below to clarify and resolve your concerns.\n\n**1. Kolmogorov complexity and MDL**\n\nThere are two problems [1] when directly applying the Kolmogorov criterion in Eq. (4): \n1. We do not have access to the true generating models $P_{X}$ and $P_{Y \\mid X}$, and \n2. The Kolmogorov complexity is not computable in practice.\n\nThe former problem can be resolved by estimating the model through the joint complexity $K(x, P)$, which on expectation over $P(x)$ will yield $K(P) + H(P)$ up to an additive constant [2]. The latter problem requires an approximating codelength $L(x, P)$ that mirrors $K(x, P)$, commonly selected according to MDL or MML principles. Despite the performance achieved, the gap between the practical $L(x, P)$ and the theoretical $K(x, P)$ is still an open problem in MDL-based causal discovery methods [1, 3], which we have not been able to rectify at the moment.\n\nHowever, although our method begins from the Kolmogorov formulation, it is important to note that the identifiability of our models are studied from an orthogonal perspective of marginal likelihoods, which should not be theoretically impacted by the MDL approximation.\n\n**2. MDL of Bayesian Neural Networks**\n\nWe would like to clarify that the prior over the parameters in our method are Gaussian distributions with zero means and learned variances (as described in App. B and C.3). We have also analyzed the normal-Jeffreys sparsity-inducing prior in Tab. 3, App. F.2, which indicated no significant improvements in performance. Additionally, Bayesian learning has already been shown to follow Occam's razor in model selection [4]. Hence, the sparsity-inducing prior is not necessary in our method.\n\nRegarding the number of layers, as our model with one hidden layer has already achieved adequate performance, we did not consider including additional layer(s). We have performed additional experiments with two hidden layers. However, there are no substantial differences in performance when the number of layer is increased. Moreover, with one more layer, the AUROC scores on some datasets such as SIM, SIM-c, and Tübingen are slightly decreased.\n\nAccuracy:\n\n| Hidden nodes | AN | AN-s | LS | LS-s | MN-U | SIM | SIM-c | SIM-G | SIM-ln | CE-Multi | CE-Net | CE-Cha | Tübingen |\n|:-|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|\n| 20 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.90 | 0.93 | 0.99 | 1.00 | 0.89 | 0.97 | 0.90 | 0.87 |\n| 10, 5 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.87 | 0.92 | 0.99 | 1.00 | 0.91 | 0.96 | 0.91 | 0.88 |\n| 10, 10 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.89 | 0.93 | 0.98 | 1.00 | 0.90 | 0.96 | 0.91 | 0.87 |\n| 20, 5 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.88 | 0.92 | 0.99 | 1.00 | 0.90 | 0.97 | 0.89 | 0.87 |\n| 20, 10 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.87 | 0.91 | 0.99 | 1.00 | 0.91 | 0.97 | 0.90 | 0.89 |\n\nAUROC:\n\n| Hidden nodes | AN | AN-s | LS | LS-s | MN-U | SIM | SIM-c | SIM-G | SIM-ln | CE-Multi | CE-Net | CE-Cha | Tübingen |\n|:-|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|\n| 20 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.97 | 0.99 | 1.00 | 1.00 | 0.98 | 1.00 | 0.97 | 0.97 |\n| 10, 5 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.96 | 0.99 | 1.00 | 1.00 | 0.98 | 1.00 | 0.97 | 0.96 |\n| 10, 10 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.97 | 0.99 | 1.00 | 1.00 | 0.98 | 1.00 | 0.98 | 0.96 |\n| 20, 5 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.96 | 0.98 | 1.00 | 1.00 | 0.98 | 1.00 | 0.97 | 0.96 |\n| 20, 10 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.96 | 0.98 | 1.00 | 1.00 | 0.98 | 1.00 | 0.97 | 0.96 |\n\n**3. References to be discussed**\n\nThank you for your recommendations. We will include additional related works in multivariate causal discovery in our multivariate discussion in App. G. As the aim of our work is the evaluation the codelengths/complexity of neural networks rather than estimating Bayesian neural networks (BNNs), we did not include a literature review on the BNNs, which we will examine in our future work. However, we do discuss related choices for computing the complexity of neural networks in Sec. 4.1.\n\n[1] Kaltenpoth, D., & Vreeken, J. (2023). Causal discovery with hidden confounders using the algorithmic Markov condition. In UAI.\n\n[2] Marx, A., & Vreeken, J. (2022). Formally Justifying MDL-based Inference of Cause and Effect. In AAAI ITCI'22 Workshop.\n\n[3] Marx, A., & Vreeken, J. (2019). Telling cause from effect by local and global regression. Knowl. Inf. Syst., 60.\n\n[4] Rasmussen, C., & Ghahramani, Z. (2000). Occam's razor. In NIPS.", "claims_and_evidence": "The proposed method aligns with the claims made in the paper.", "methods_and_evaluation_criteria": "Empirical evaluations on both synthetic and real-world datasets validate the effectiveness of the proposed method.", "theoretical_claims": "I reviewed the theoretical aspects at a high level but did not rigorously verify the correctness", "experimental_designs_or_analyses": "The provided analyses are generally well-conducted.", "relation_to_broader_scientific_literature": "The paper addresses an important problem, which has received limited attention in prior research.", "essential_references_not_discussed": "The authors may discuss related works that address general non-linear relationships by leveraging the non-stationarity of observations, such as Monti et al., \"Causal Discovery with General Non-Linear Relationships Using Non-Linear ICA\", and Huang et al., \"Causal Discovery from Heterogeneous/Nonstationary Data\" (JMLR, 2020). Additionally, since Bayesian neural networks are mentioned, it would be helpful to include a discussion in the related work section on their relevance, particularly in the context of causal discovery and uncertainty quantification, e.g., Bayesian causal discovery mentioneds.", "other_strengths_and_weaknesses": "It is expected that some methods consider non-parametric functions beyond functional class-based approaches for causal discovery, as the latter are often difficult to validate in real applications. While I agree with Postulates 1 and 2 regarding the causal asymmetry from the Kolmogorov complexity (though I am unsure how widely accepted they are in the community), I have the following concerns:\n\n1) The gap between Kolmogorov complexity and MDL:  It is well known that Kolmogorov complexity is incomputable, meaning there is no general algorithm to compute the exact value of Kolmogorov complexity for arbitrary data. In this case, how can an approximation method be used to provide a solution? In other words, if an approximation exists, does that mean the problem is constrained in a way that makes it well-defined and solvable? Furthermore, does such a constrained Kolmogorov complexity remain consistent with the ability to determine causal asymmetry, similar to standard Kolmogorov complexity? In other words, do the Postulates still hold under this constrained Kolmogorov complexity?\n\n2) The gap between MDL and Bayesian Neural Networks: Given that this work adopts an MDL approach, it is reasonable to use a Bayesian framework for model selection. However, there has been significant recent progress in Bayesian neural networks, particularly regarding compact model selection, such as using sparsity-inductive priors. In other words, when designing a model, one may focus on how to learn a compact model that aligns with the prior knowledge inherent in the MDL framework. However, this work appears to introduce a relatively simple sample prior, a Gaussian distribution with zero mean and unit variance. How can this model effectively learn a sparse representation? Furthermore, how should the number of layers in the neural network be determined to maximize the ability to learn a sparse model from the data?\n\nOverall, the gap between Kolmogorov complexity and MDL raises concerns about whether causal asymmetry holds under the MDL approximation. Furthermore, the implementation of methods—such as using a relatively simple sample prior, specifically a Gaussian distribution with zero mean and unit variance—further deepens these concerns. \n\nI find this work compelling, particularly its exploration of a practical direction for causal discovery. However, the concerns raised above have led me to question the significance of the contribution. I would be happy to revise my rating if these issues are addressed effectively in a revised version of the work.", "other_comments_or_suggestions": "none", "questions_for_authors": "none", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "aK4PMKzKrj", "reviewer": "Reviewer_YARz", "summary": "The paper introduces COMIC (Causal direction identification via Bayesian COMpression), a novel method for determining causal relationships between pairs of variables using variational Bayesian compression with neural networks. This approach improves upon existing compression-based methods by balancing model fitness and complexity through variational inference, demonstrating superior performance across multiple synthetic and real-world benchmarks. COMIC effectively models conditional distributions while considering both model accuracy and complexity, achieving state-of-the-art results.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer YARz for your novelty and favorable qualities of our work. Our responses below aim to rectify the issues that you mentioned.\n\n**1. Formulating the capability of the neural networks in Prop. 4.2**\n\nThank you for your recommendation, we will specifically formalize the function of the \"single-hidden-layer neural networks\" as described in App. C.3 to make our proposition more mathematically rigorous.\n\n**2. Theorem 3.1**\n\nSince Theorem 3.1 is directly related to Eq. (5) and (7) and is the criterion for determining the causal direction, it was presented to clarify the origin of our causal scores and make our paper self-contained and easier to follow.\n\n**3. Hyperparameter configurations in location-scale models**\n\nSimilar to the additive noise models, our location-scale models are also Bayesian neural networks with one hidden layer, whose details on hyperparameters are presented in App. B and C.3. The only difference between these types of models is that the location-scale ones have an additional output node to regress the log scale parameters for the Gaussian likelihood for enhancing the ability to model aleatoric uncertainty of our models.\n\n**4. Putting more experiment results in the main paper**\n\nDue to the current limit of space, we can only accommodate the most crucial results in the main paper. We will consider present directly or including more specific discussions and references to these key results in the next version of our manuscript. \n\n**5. Choice of prior for Bayesian neural networks**\n\nOur selection of Gaussian priors is influenced by previous works on evaluating neural network's complexity, especially [1-3]. As a result, we utilize Gaussian priors, which is a common choice, for variational Bayesian learning of neural networks. These studies have shown that variational learning with the Gaussian priors and hyperpriors on the variances yields adequate and consistent results with respect to the implicit prequential coding, which has been demonstrated to achieve good compression bounds in [2] and [3].\n\n[1] Louizos, C., Ullrich, K., & Welling, M. (2017). Bayesian compression for deep learning. In NeurIPS.\n\n[2] Blier, L., & Ollivier, Y. (2018). The description length of deep learning models. In NeurIPS.\n\n[3] Voita, E., & Titov, I. (2020). Information-Theoretic Probing with Minimum Description Length. In EMNLP.", "claims_and_evidence": "I list several core claims of the paper\n1. COMIC improves model fitness while promoting succinct codelengths: The method achieves superior performance on both synthetic and real-world benchmarks, outperforming related complexity-based and structural causal model regression-based approaches.\n2. The variational Bayesian coding scheme effectively approximates the algorithmic complexity of neural networks Theoretical derivation shows that variational coding length can be decomposed into model fitness (ELBO) and complexity terms (KL divergence), with experimental validation demonstrating its effectiveness in practice.\n3.  COMIC is identifiable for causal direction: Based on Bayesian model selection theory, the authors prove that their models are non-separable-compatible, ensuring causal direction can be distinguished by marginal likelihood differences.", "methods_and_evaluation_criteria": "The paper utilizes neural networks to learn conditional distributions, overcoming limitations of traditional regression methods. It employs variational inference to assess neural network complexity, explicitly capturing model complexity while avoiding high computational costs.\nThe causality is determined by comparing total coding lengths (marginal and conditional) for different directions.\nEvaluation criteria: mainly on proportion of correctly identified causal directions (auruc). This criteria makes sense.", "theoretical_claims": "I find no obvious problems.", "experimental_designs_or_analyses": "The method is compared with a set of other ones that at least relates to the complexity based principles. Also, they conduct location-scale modeling vs. location-only modeling: The former shows superior performance on complex datasets (e.g., 97% Bi-AUROC on Tübingen vs. significant drops for the latter), and the impact of model complexity: Ignoring KL divergence (optimizing only likelihood) leads to performance degradation (e.g., accuracy drops from 90% to 83% on SIM dataset). Hidden layer width is also studied. This is a factor that relates to \"how well\" the complexity can be approximated by the neural statistical approaches. I thus thing the evaluations are comprehensive.", "relation_to_broader_scientific_literature": "This method contributes a practical approach for causal discovery. It is based on the statistical approximation of KM complexity, and the ICM principle. Although the princicple is already proposed for years, this new method provides a bridge connecting the uncomputable KM to a statistical approximations, which makes causal discovery from observational data tractable. Its perspective of approximation KM using neural nets with appropriate optimization method is, to me, novel.", "essential_references_not_discussed": "The citations are appropriate.", "other_strengths_and_weaknesses": "The method is practical. On the other hand, the approximation limits are not taht clear, although it does not affect the main conclusions of the paper.", "other_comments_or_suggestions": "It would be better to come out with a compehensive approach on how to select to priror of the Bayes method.", "questions_for_authors": "1. Prop 4.2: I do not think \"modeled by single-hidden-layer neural\" should appear in  a math claim. This should be formulated by the capability of the neural nets or function to be more rigorous.\n2. Theorem 3.1: This better appears in appendix since it is already proven.\n3. Appendix C.3: concerning location-scale model, is there some specific Bayes hyperparameter configuations?\n4. You may also consider putting several core experiments regarding the key aspects of the paper, say, to what extend can the Bayes method approximates the uncomputable complexity, in main paper, rather than on appendix.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "XoU96jUBci", "reviewer": "Reviewer_3osf", "summary": "This work proposes using Bayesian neural networks with variational inference (I will refer to this as variational BNNs) for bivariate causal discovery via the MDL principle (COMIC). The causal direction is determined by the direction that best trades off model complexity with model fit under the factorization $P = P_{\\text{cause}} P_{\\text{effect}\\mid\\text{cause}}$. COMIC parametrizes the conditional model as a variational BNN and takes advantage of the fact that minimizing the ELBO can be directly interpreted as description length, and thus variational inference in this case performs MDL. The authors adapt recent identifiability theory from a related paper that also uses Bayesian model selection in this setting (Dhir et al., 2024) to show identifiability. Experimental results show that COMIC outperforms most competitors on a large set of benchmark datasets.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate Reviewer 3osf's comments to further elucidate and enhance our work. The following responses address the major points that you raised.\n\n**1. Relation to Appendix B.2 of [1]**\n\nOur work originates from the causal discovery via MDL using neural networks, with variational Bayesian codelength being the encoding method. This approach latter converges to the marginal likelihood-based method outlined in GPLVM [1]. Although the two-part MDL formulation is briefly discussed in [1, App. B.2], to the best of our knowledge, our work is the first one to *explicitly* integrate this variational MDL codelengths into the learning process of neural networks for the downstream task of causal discovery. This approach fundamentally differs from GPLVM, which originally adopts a non-parametric approach to model inference and only considers the formulation with respect to the model complexity in a post-hoc manner. While both methods employ the variational learning, GPLVM uses the KL divergence terms to infer the distributions over latent variables and inducing points of sparse GPs [1, App. F.1], rather than optimizing the conventional model parameters. In contrast, our approach *directly* incorporates the model's complexity into the learning process by considering the complexity via the distributions over the parameters.\n\n**2. $N(0,1)$ for the marginals**\n\nAs encoding the conditionals is our main focus in this work, we choose $N(0,1)$ as an uninformative choice to encode the marginals. In fact, excepts for the AN, AN-s, LS, LS-s, MN-U, and SIM-G pairs, the remaining benchmarks consist of pairs with a non-Gaussian cause, including the real-world Tübingen dataset. We have run additional analysis on different choices for the marginals including $N(0, 1)$, $U(X_{min}, X_{max})$, and Variational Bayesian Gaussian mixture model (GMM) [2]. The uniform encoding negatively affects the performance on most datasets. GMMs provide advantages on SIM, CE-Multi, and CE-Net synthetic benchmarks, whereas the performance on SIM-G, SIM-ln, and CE-Cha benchmarks are noticeably reduced. Although there is a slight increase on the Tübingen prediction accuracy, the AUROC score is also substantially decreased. As a result, we believe the standard Gaussian remains a reasonable choice for our method. Additionally, using a more flexible model for encoding the marginals will also challenge the theoretical analysis of our models' separable-compatibility.\n\nAccuracy:\n\n| Marginal | AN | AN-s | LS | LS-s | MN-U | SIM | SIM-c | SIM-G | SIM-ln | CE-Multi | CE-Net | CE-Cha | Tübingen |\n|:-|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|\n| Gaussian | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.90 | 0.93 | 0.99 | 1.00 | 0.89 | 0.97 | 0.90 | 0.87 |\n| Uniform | 1.00 | 0.93 | 1.00 | 1.00 | 1.00 | 0.88 | 0.84 | 0.83 | 0.99 | 0.94 | 0.94 | 0.80 | 0.86 |\n| GMM | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.97 | 0.94 | 0.92 | 0.94 | 0.98 | 1.00 | 0.80 | 0.89 |\n\nAUROC:\n\n| Marginal | AN | AN-s | LS | LS-s | MN-U | SIM | SIM-c | SIM-G | SIM-ln | CE-Multi | CE-Net | CE-Cha | Tübingen |\n|:-|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|-:|\n| Gaussian | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.97 | 0.99 | 1.00 | 1.00 | 0.98 | 1.00 | 0.97 | 0.97 |\n| Uniform | 1.00 | 0.98 | 1.00 | 1.00 | 1.00 | 0.96 | 0.91 | 0.93 | 1.00 | 0.99 | 0.98 | 0.89 | 0.95 |\n| GMM | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.98 | 0.98 | 0.91 | 1.00 | 0.99 | 0.81 | 0.87 |\n\n[1] Dhir, A., Power, S., & Van Der Wilk, M. (2024). Bivariate Causal Discovery using Bayesian Model Selection. In ICML.\n\n[2] Blei, D. M., & Jordan, M. I. (2006). Variational Inference for Dirichlet Process Mixtures. Bayesian Anal., 1(1).", "claims_and_evidence": "To my judgement, the claims are well-supported.", "methods_and_evaluation_criteria": "The proposed method builds on a large body of literature using the MDL principle for causal discovery. Given the connection between variational BNNs and MDL, and given recent work that Bayesian methods provably give sharper identifiability (Dhir et al., 2024), the method proposed is timely and makes a lot of sense. The benchmarks used are standard in this area and go beyond what was done in the previous work by Dhir et al.", "theoretical_claims": "The claim is that the variational BNN model with Gaussian likelihoods is not separable--compatible, which yields causal direction identifiability via the marginal likelihoods in the sense of Dhir et al., 2024. \n\nI found this to be a weak part of the paper. To my understanding, the model classes in Proposition 4.2 represent the class of Gaussian LSNMs, but 1) the proof is only given rigourously for ANMs, which are essentially identifiable in the non-Bayesian sense, which according to Dhir et al. implies identifiability in the Bayesian sense, and the same for LSNMs (Immer et al., 2023). Also, I did not check the proof closely (I am not that familiar with the idea of separable--compatible), but I couldn't actually tell from the proof why it was necessary to analyze the BNNs as GPs (which only holds in the infinite--width limit), it seems the Gaussianity of the observation noise would be enough (if the identifiability of ANMs with equal noise variance in general were not enough already!)", "experimental_designs_or_analyses": "I think the benchmarks/experimental design is pretty standard for bivariate causal discovery.", "relation_to_broader_scientific_literature": "I think the relevance of bivariate causal discovery to the broader scientific literature is quite clear at this point, but it may be worthwhile to point out the relation to multivariate causal discovery (e.g., how to use this method directly or that one can first use bivariate causal discovery to orient edges in CPDAGs).", "essential_references_not_discussed": "Appendix B.2 of Dhir et al., 2024 also connect the MDL principle to evidence maximization in Bayesian models.", "other_strengths_and_weaknesses": "### Strengths\n\nI mentioned above that I think this is generally a timely paper that spells out the known connection between variational BNNs and MDL to connect MDL-based causal discovery (older idea in this area) to Bayesian model selection-based causal discovery (a new idea with new theoretical guarantees). It seems that this connection may have been noted by previous work (Dhir et al., 2024), but it was brief and buried in the appendix---my judgement is that it is worthwhile that the current submission expands on this. The other strength is that the experimental results really appear to be SOTA on these common benchmarks, improving already impressive performance by GPLVM. \n\n### Weaknesses\n\nAs a methodology, the novelty is limited as it can be seen as a version of the GPLVM-based method (Dhir et al., 2024) replacing the marginal likelihood with the ELBOs (a lower bound), and replacing the GP with a BNN. The identifiability theory also closely follows the previous paper.", "other_comments_or_suggestions": "- I would have preferred more background in Section 4.4, e.g., it was not clear what \"priors over anti-causal factorizations\" were in Definition 4.1 without reading the reference Dhir et al., 2024.\n- Under Corollary 4.3, I would also include the qualifier that we need the marginal model to be correctly specified as $P_{\\text{cause}} = N(0,1)$.", "questions_for_authors": "For me, a major conceptual contribution of the paper is that it spells out this tight connection between the MDL principle and Bayesian inference. However, I am not an expert in these fields, so one main question that I have which will influence my score is:\n\n1. Can you clarify to what extent your description of the connection between Bayesian causal discovery and the MDL principle extends that of Appendix B.2. of Dhir et al., 2024? \n\nOther questions (ordered in terms of importance):\n\n2. I understand using the same marginal model for the cause variable to be fair to both directions (Eqn. 14), but would this not just prefer the \"more Gaussian\" variable for the cause? Is this a reasonable bias for real-world data? \n3. Can you clarify the contribution of your identifiability proof compared to the known existing results of (non-Bayesian) identifiability of ANM/LSNM? \n\nI am not particularly picky about (3), but if the bias in (2) is prominent I think it would be better to probe the sensitivity of the method to this (e.g., what does it imply for non-Gaussian causes)?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper views bi-variate causal discovery as a model selection problem under Kolmogorov complexity, and proposes a practical approximation to it via a variational Bayesian approach based on the minimum description length principle. Specifically, the idea is to compare the compression code length of L(p(x)) + L(p(y|x)) vs L(p(y)) + L(p(x|y)), where the code for the marginal distribution is computed using a standard Gaussian distribution for entropic encoding, and the code length for the conditional distribution is computed via the variational lower-bound of a Bayesian neural network fitted to that conditional distribution. Experiments on various bi-variate causal discovery benchmarks show the proposed method's effectiveness when compared with a number of baseline approaches.

Reviewers found the proposed approach novel and the presentation is sound. They have a few questions regarding the gap between Kolmogorov complexity and MDL, as well as questions regarding the Bayesian neural network methodological designs. Author rebuttal has addressed most of these concerns.

I briefly read the paper and I personally think that the proposed is very interesting and should be presented to the causal ML community at ICML. The only major concern I have is regarding the entropic coding of the marginal distribution under standard Gaussian, which sounds like a crude approximation when the cause is non-Gaussian. However bi-variate causal discovery becomes easier when having non-Gaussian variables, and I guess this partly connects to the empirical observation that Gaussian entropic encoding seems already enough?

Also in Bayesian neural network literature, Laplace methods are empirically shown to perform better regarding model selection and hyper-parameter optimisation. It would be great to comment on this point since the authors choose variational lower-bound as the code length for the conditional distribution.
Decision: Accept (spotlight poster)