Identifying Causal Direction via Dense Functional Classes

TL;DR

This paper introduces LCUBE, a fast and interpretable method for causal inference between two variables using the MDL principle and cubic splines, outperforming existing methods on real-world datasets.

Contribution

The paper proposes a novel MDL-based causal score using cubic splines, with proven identifiability and practical advantages over prior approaches.

Findings

LCUBE achieves superior precision on Tuebingen cause-effect dataset.

It outperforms state-of-the-art methods on multiple benchmarks.

The method is simple, interpretable, and fast with only one hyperparameter.

Abstract

We address the problem of determining the causal direction between two univariate, continuous-valued variables, X and Y, under the assumption of no hidden confounders. In general, it is not possible to make definitive statements about causality without some assumptions on the underlying model. To distinguish between cause and effect, we propose a bivariate causal score based on the Minimum Description Length (MDL) principle, using functions that possess the density property on a compact real interval. We prove the identifiability of these causal scores under specific conditions. These conditions can be easily tested. Gaussianity of the noise in the causal model equations is not assumed, only that the noise is low. The well-studied class of cubic splines possesses the density property on a compact real interval. We propose LCUBE as an instantiation of the MDL-based causal score utilizing cubic regression splines. LCUBE is an identifiable method that is also interpretable, simple, and very fast. It has only one hyperparameter. Empirical evaluations compared to state-of-the-art methods demonstrate that LCUBE achieves superior precision in terms of AUDRC on the real-world Tuebingen cause-effect pairs dataset. It also shows superior average precision across common 10 benchmark datasets and achieves above average precision on 13 datasets.

Tables

Table 1. TABLE I : 1. column: AUDRC on real data of Tü pairs. 2. column: Average AUDRC/ACC over 10 datasets, 2. column: Average AUDRC/ACC over 13 datasets

Method	AUDRC $\uparrow$	Aver.AUDRC/ACC $\uparrow$	Aver. AUDRC/ACC $\uparrow$
LCUBE	87	91.5    87.5	85.5    81.4
Sloppy	75	83.7    81.1	78.9    73.8
CAM	70	79.4    78.2	73.8    72.5
QCCD	84	89.7    83.2	85.7    78.2
RESIT	71	65.7    65.0	68.4    64.4
IGCI	74	42.6    40.0	49.5    44.7
LOCIM	45	80.3    78.1	78.2    74.4
LOCIH	47	82.9    78.6	80.1    77.3
HECI	78	71.4    62.9	73.4    65.3
NN-LOCIM	66	87.3    81.2	84.5    77.1
NN-LOCIH	56	91.3    85.6	89.1    84.1
ROCHE	74	85.6    83.1	84.2    82.3

Table 2. TABLE II : AUDRC in % \% for concrete datasets and methods. CAM to NN-LOCI H are from [ 23 ] .

	Tü	AN	ANs	LS	LSs	MNU	SIM	SIMc	SIMln	SIMG	Multi	Net	Cha
LCUBE	87	100	100	100	99	100	75	76	93	85	50	88	58
Sloppy	75	100	100	100	70	99	65	63	80	85	70	61	58
CAM	70	100	100	100	36	76	68	69	87	88	85	40	41
QCCD	84	100	91	100	100	100	71	83	92	76	63	94	61
RESIT	71	100	100	69	4	0	75	84	83	71	68	81	83
IGCI	74	18	35	60	49	1	34	41	51	63	99	61	58
LOCIM	45	98	95	88	86	90	68	56	90	87	75	84	55
LOCIH	47	100	96	92	54	92	84	85	88	91	51	93	68
HECI	78	100	63	99	72	20	59	64	86	73	99	84	57
NN-LOCIM	66	100	100	100	100	100	60	63	95	89	93	86	47
NN-LOCIH	56	100	100	99	97	100	89	93	86	93	77	97	71
ROCHE	74	100	100	100	100	65	76	85	78	78	80	81	78

Table 3. TABLE III : Accuracy in % \% for concrete datasets and methods. CAM to NN-LOCI H are from [ 23 ] .

	Tü	AN	ANs	LS	LSs	MNU	SIM	SIMc	SIMln	SIMG	Multi	Net	Cha
LCUBE	72	100	100	100	98	100	67	68	90	80	45	81	57
Sloppy	76	100	100	100	56	96	64	62	77	81	50	46	52
CAM	58	100	100	100	53	86	57	60	87	81	35	78	47
QCCD	77	100	82	100	96	99	62	72	80	64	51	80	54
RESIT	57	100	100	61	6	2	77	82	87	78	37	78	72
IGCI	68	20	35	46	34	11	37	45	51	53	92	35	55
LOCIM	52	99	98	94	94	93	52	48	77	74	66	75	46
LOCIH	56	99	98	85	53	90	75	76	73	81	66	84	70
HECI	71	98	55	92	55	33	49	55	65	56	91	72	57
NN-LOCIM	57	100	100	100	100	100	48	50	79	78	72	76	43
NN-LOCIH	60	100	100	95	89	100	79	83	72	78	78	87	72
ROCHE	70	100	100	100	100	62	71	80	74	74	78	81	80

Equations

$$Y_{\alpha }:=f(X)+\alpha N,$$

$${\tilde{Y}}_{\alpha }:=\frac{Y_{\alpha }-\alpha r_{-}}{1+\alpha r_{+}-\alpha r_{-}}.$$

$$Y_n:=f(X)+\frac{1}{n}N$$

$$Y_n\approx h(X)+\frac{1}{n}N.$$

$$y_i\approx f(x_i)\approx h(x_i,\mathbf{\theta }),i=1,\dots ,n.$$

$$K(P(X))+K(P(Y|X))\stackrel{+}{\le }K(P(Y))+K(P(X|Y))$$

$$M^{\ast }=argmi{n}_{M\in \mathcal{M}}L(M)+L(D|M)$$

$$\underset{{\mathbf{\theta }}_X}{\min }(L(X)+L(Y|X,{\mathbf{\theta }}_X))\stackrel{+}{\le }\underset{{\mathbf{\theta }}_Y}{\min }(L(Y)+L(X|Y,{\mathbf{\theta }}_Y))$$

$$\underset{{\mathbf{\theta }}_X}{\min }L(Y|X,{\mathbf{\theta }}_X)\stackrel{+}{\le }\underset{{\mathbf{\theta }}_Y}{\min }L(X|Y,{\mathbf{\theta }}_Y).$$

$$\underset{n\to \infty }{\lim }\frac{L({\tilde{Y}}_n|X,{\hat{\mathbf{\theta }}}_X)}{L(X|{\tilde{Y}}_n,{\hat{\mathbf{\theta }}}_Y)}\le 1,$$

$$S(Y\mid X,\phi ):=\gamma (E[(Y-\phi (X){)}^2])+\lambda (\Vert {\beta }_{\phi }{\Vert }_0)$$

$$\underset{n\to \infty }{\lim }\frac{S(E[({\tilde{Y}}_n-\phi (X){)}^2],\Vert {\beta }_{\phi }{\Vert }_0)}{S(E[(X-\psi ({\tilde{Y}}_n){)}^2],\Vert {\beta }_{\psi }{\Vert }_0)}\le 1,$$

$$E[({\tilde{Y}}_n-\phi (X){)}^2]\approx \frac{RSS({\hat{\mathbf{\theta }}}_X)}{n}$$

$$\text{mbox}andE[(X-\psi ({\tilde{Y}}_n){)}^2]\approx \frac{RSS({\hat{\mathbf{\theta }}}_Y)}{n}$$

$$f(x)\approx s_f(x)=b_0+b_{1}x+b_2{x}^2+b_3{x}^3+\sum _{j=1}^m{\beta }_j(x-k_j{)}_{+}^3.$$

$$\mathbf{\theta }=(m,\mathbf{k},\mathbf{b},\mathbf{\beta })$$

$$\mathbf{X}:=(1,\mathbf{x},{\mathbf{x}}^2,{\mathbf{x}}^3,(\mathbf{x}-k_{1}1{)}_{+}^3,\dots ,(\mathbf{x}-k_{m}1{)}_{+}^3)$$

$$({\hat{\mathbf{b}}}^{\top },{\hat{\mathbf{\beta }}}^{\top })=({\mathbf{X}}^{\top }\mathbf{X}{)}^{+}{\mathbf{X}}^{\top }\mathbf{y},$$

$$Y\approx s_f(X)$$

$$X\approx s_g(Y)$$

$$L(Y|X,\mathbf{\theta }):=L(\mathbf{\theta })+L(Y|\mathbf{\theta }).$$

$$L(\mathbf{\theta })=L(m)+L(\mathbf{k}|m)+L(\mathbf{b},\mathbf{\beta },|m,\mathbf{k})$$

$$L(\mathbf{\theta })=\log m+\sum _{j=1}^m\log u_j+\frac{m+4}{2}\log n$$

$$u_j:=\text{mbox}thenumberof{x}_i^{\prime }s\text{mbox}sothat{k}_{j-1}\le x_i\le k_j$$

$$L(Y|\mathbf{\theta })=\frac{n}{2}\log (\frac{RSS(\mathbf{\theta })}{n})+C$$

$$L(Y|X,\mathbf{\theta })=L(\mathbf{\theta })+L(Y|\mathbf{\theta })=\log m+\sum _{j=1}^m\log u_j+\frac{m+4}{2}\log n+\frac{n}{2}\log (\frac{RSS(\mathbf{\theta })}{n})$$

$$\underset{n\to \infty }{\lim }\frac{{\hat{\Delta }}_n(X\to Y)}{{\hat{\Delta }}_n(Y\to X)}\le 1,$$

$$\lambda (\Vert {\beta }_{\phi }{\Vert }_0):=\log m+\sum _{j=1}^m\log u_j^X+\frac{m+4}{2}\log n=\log (m\times \prod _{j=1}^m{u}_j^X\times n^{\frac{m+4}{2}}).$$

$$\underset{\alpha \to 0}{\lim }\frac{\mathbb{E}[({\tilde{Y}}_{\alpha }-\phi (X){)}^2]}{\mathbb{E}[(X-\psi ({\tilde{Y}}_{\alpha }){)}^2]}\le 1$$

$$\underset{n\to \infty }{\lim }\frac{\mathbb{E}[({\tilde{Y}}_n-\phi (X){)}^2]}{\mathbb{E}[(X-\psi ({\tilde{Y}}_n){)}^2]}\le 1.$$

Full text

I Introduction

We deal with a problem to decide about two univariate continuous-valued variables $X$ and $Y$ without hidden confounders, which variable is causing another one. In general, it is not possible to make definite statements about causality without some assumptions on the underlying model [33]. Recently, Blöbaum et al. in [31] formulated a set of assumptions under which the true causal direction is identifiable via regression. They showed that it is possible to identify cause from effect by selecting the direction with smaller residual error. The idea is to fit both $Y=f(X)+N_{X}$ and $X=g(X)+N_{Y}$ by minimizing the respective residual errors $N_{X}$ and $N_{Y}$, and then to compare the sum of squared errors. They showed that, among models with complexity matching that of the true model, the model in the correct causal direction yields a lower error. However, since the complexity of the true model is not know, RECI may compare residuals of arbitrarily underfit or overfit models. Paper [30] builds on the results of [31], and formulates an additional condition under which cause and effect are identifiable via so-called regularized regression. The so-called Identifiable Regression-based Scoring Function (= IRSF) is proposed and a general framework for the identifiability of these functions is presented, including an instantiation by method SLOPPY, using splines and AIC and BIC as scoring function. Method Slope [28] also uses spline regression, however without identifiability conditions. Methods based on regression error are able to decide between Markov equivalent DAGs under the assumption of having a non-linear function and low noise [31]. The method we present in this paper will also adopt these assumptions.

In this paper, we propose a causal score using the Minimum Description Length (MDL) of regression models via functions from dense functional classes. A functional class is dense in the set of all continuous functions $C([a,b])$ on a compact interval $[a,b]$ if it can approximate any continuous function on $[a,b]$ within arbitrary precision. Under certain assumptions on the causal models, we prove identifiability of these causal scores using MDL of functions from dense classes. We propose method LCUBE as an instantiation of this MDL-based score for the cubic regression splines which possess the density property. LCUBE uses MDL encoding that allows more detailed characterization of cubic splines models than the encoding via AIC or BIC criteria.

The key contributions of this work can be outlined as follows:

• •

We introduce a bivariate causal score based on the MDL principle, using functions that possess the density property on a compact real interval. We prove the identifiability of these causal scores under specific conditions, which can be easily tested in practice.

• •

We propose LCUBE as an instantiation of the MDL-based causal score using cubic regression splines. LCUBE is an identifiable method that is also interpretable, simple and very fast.

• •

Empirical evaluations compared to state-of-the-art methods demonstrate that LCUBE achieves superior precision in terms of AUDRC on the real-world Tübingen cause-effect pairs dataset. It also outperforms the methods in terms of average precision across 10 common benchmark datasets and achieves above-average precision on 13 datasets.

II Preliminaries

Recent literature on causal inference imposes certain model restrictions in order to to identify the cause and effect reliably. For example, [41], [34] provide theoretical foundations to establish identifiability of non-linear additive noise models. In our work, we assume that causal relations between $X$ and $Y$ can be expressed by a non-linear bivariate functional model, an instance of the Additive Noise Models (ANM), [41], and that there are no other causes of $Y.$ Our models further rely on the assumptions on regression from [31]. Blöbaum et al. in [31] formulated the following three assumptions under which the true causal direction is identifiable via regression.

Assumption 1**.**

[31]** (Causal model). The effect can be written as

[TABLE]

with noise term $N$ and parameter $\alpha$ restricting the noise level.

Assumption 2**.**

[31]** (Unbiased noise). The noise term $N$ is unbiased and has unit variance.

Assumption 3**.**

[31]** (Compact supports). The distribution of $X$ has compact support and w.l.o.g. $X$ attains values between 0 and 1, which can be achieved by normalizing $X$. Further, the distribution of $N$ has compact support and there exist values $r_{+}>0>r_{-}$ such that for each value $x\in X$, $[r_{-},r_{+}]$ is the smallest interval containing the support for the conditional density of $N$ given $x$. Hence, we know that $[\alpha r_{-},1+\alpha r_{+}]$ is the smallest interval containing the support of the density of $Y_{\alpha}$ and rescale it to

[TABLE]

$\tilde{Y}_{\alpha}$* has the same scale as $X$ and has values between 0 and 1.*

Now consider an integer $n\geq 2$ integer and $\alpha:=\frac{1}{n}$ Eq. (1). With a slight abuse of notation, we write $Y_{n}$ instead of $Y_{\frac{1}{n}}$, and define the following ANM model

[TABLE]

with noise $N$. One can see that (2) is a special case of (1). We stress here that the noise variable $N$ in Eq. (1). does not have to be necessarily Gaussian.

II-A Identifiability of likelihood score-based methods in causal discovery

Schultheiss and Bühlmann in [21] discuss causal discovery in structural causal models using Gaussian likelihood scoring and analyze the effect of model misspecification in a general directed acyclic graphs. They came to the following conclusions in the case where the data-generating distribution comes from a linear structural equation model and linear regression functions are used for estimation. When the true error distribution is Gaussian, one can only identify the Markov equivalence class of the underlying data-generating DAG. The same holds true when the error distribution is non-Gaussian, but one wrongly relies on a Gaussian error distribution for estimation. When likelihood scores in a linear structural equation model and linear regression are used for comparison of different DAGS and Gaussianity of errors is not assumed, the multivariate density of the DAGs factorize only for independent error functions. The true DAG maximizes this likelihood score due to the properties of the Kullback-Leibler divergence. In practice, this comes with the additional difficulty of estimating the densities and one needs to add some penalization to prefer simpler graphs and avoid selecting complete graphs [22].

In this paper, we face the pitfalls associated with likelihood-based causal scores by employing minimum description length (MDL)-based scores. The MDL-based scores inherently penalize more complex graph structures and do not assume Gaussianity of the error terms. As we demonstrate in our experiments using concrete MDL scores, this approach yields high causal accuracy in practice.

III Functional model with density property and MDL

III-A Density property of a functional class

In general, a subset $A$ of a topological space $X$ is said to be dense in $X$ if every point of $X$ either belongs to $A$ or is arbitrarily ”close” to a point in $A$. In functional approximation theory, the Weierstrass Approximation Theorem states that any given complex-valued continuous function defined on a closed interval $[a,b]$ can be uniformly approximated as closely as desired by a polynomial function. In other words, polynomial functions are dense in the space of continuous complex-valued functions on $[a,b]$, when equipped with the supremum norm. For the purposes of our paper we restrict ourselves to real-valued continuous functions on a closed real interval $[a,b]$, which is denoted as $C([a,b])$.

Assume now a functional causal model from Eq. (2) with error $\frac{1}{2n}N$ where we replace $f$ by its dense approximation $h$ within error $\frac{1}{2n}N$ (for the supremum norm on $C[a,b]$). So for a sufficiently large $n$, we get from Eq. (2) model

[TABLE]

Since $Y_{n}$ can be approximated by $h$ on the compact interval $[a,b]$ within error $\frac{1}{n}N,$ the guarantee of this error holds also for every subset of $y_{i}=f(x_{i}),i=1,\dots,n$, drawn from continuous $X$ on a compact interval. Therefore given an approximation of $f$ by $h$ from the dense functional class $H$, we can express the regression model on finite samples as

[TABLE]

III-B Kolmogorov complexity and MDL

To be able to distinguish between Markov equivalence classes, Janzing and Schölkopf in [2] postulated the algorithmic equivalent of the principle of independent mechanisms. Under the algorithmic model of causality, they formulated the postulate of algorithmic independence of conditionals using Kolmogorov complexity. This led to the following inference rule, as described in [3]: If $X\to Y$ is the true causal direction, then for Kolmogorov complexity $K$ of a cause $X$ and effect $Y$ given cause holds

[TABLE]

where $\overset{+}{\leq}$ denotes inequality up to additive constants. However, the Kolmogorov complexity is not computable [4] and the true distribution is not known, only a data sample is usually available. A principled way to solve at least the first part of the problem is to approximate Kolmogorov complexity via the Minimum Description Length principle (MDL). The Minimum Description Length principle is a practical variant of Kolmogorov complexity that approximates $K$ from above, see e.g., [10], [11]. The MDL principle has been applied in bivariate causal discovery for classes of probability functions, as in e.g. [5], [7]. Given data $D$, which may represent a sample from a distribution $P(X)$, the idea is to find that model $M^{*}\in{\cal M}$, such that

[TABLE]

where $L(M)$ is the length in bits required to describe the model $M$ within the model class ${\cal M}$, and $L(D|M)$ is the length in bits of the description of data $D$ given $M$. For more details to MDL and its coding, we refer the reader to [11].

In this paper we consider the model class ${\cal M}$ to be a class of functional models. We assume that each $f$ in the functional class ${\cal M}$ can be parameterized by a vector $\boldsymbol{\theta}$ and each function $f$ can be fitted to the data sample $D$. Accordingly, in (6), $L(M)$ measures the complexity of describing the parameter vector $\boldsymbol{\theta}$, and $L(D|M)$ quantifies the goodness of fit of the function $f$ given by the parameter vector $\boldsymbol{\theta}$ to on the data $D$. In this way, the causal direction rule from Eq. (5) becomes

[TABLE]

where $\boldsymbol{\theta}_{X},\boldsymbol{\theta}_{Y}$ are parameterizations of each causal model direction. If we assume for cause that it is normalized and has a uniform prior, then the codes $L(X),L(Y)$ do not depend on the parameterizations $\boldsymbol{\theta}_{X},\boldsymbol{\theta}_{Y}$, and we can instead use the rule

[TABLE]

We use the mean residual sum of squares (MRSS) as the goodness of fit. For the expected least-squared error $E(.)$, we will use the commonly adopted abbreviation ELSE.

IV

Identifiability of the MDL Causal Score via Dense Functional Classes

Let $\beta_{f}$ denote the set of parameters of a function $f$ and let $\|\beta_{f}\|_{0}$ denote its number of non-zero parameters. We will further assume for cause that it is normalized and has a uniform prior.

Assumption 4**.**

[30]** (Simplicity). Let $Y_{\alpha}$ be generated as in Assumption 1. Let $\varphi$ be the function minimizing the ELSE for predicting the effect $Y$ from the cause $X$ and $\psi$ be the function minimizing the ELSE in the anti-causal direction. We assume that $\psi$ has at least as many parameters as $\varphi$, i.e. $\|\beta_{\varphi}\|_{0}\leq\|\beta_{\psi}\|_{0}$.

And now we can present our result:

Theorem 1**.**

Let $H$ be a dense functional model class where each function from $H$ is parameterized by vector $\boldsymbol{\theta}$. Assume the Assumptions 1-4 for causal models given by Eq. (3) between cause and effect are satisfied. Let the encoding of the parameter vector $\boldsymbol{\theta}$ in the MDL description, for a sample set of size $n$, be given by the logarithm of a positive integer function that is strictly increasing with $n$. Assume that the goodness-of-fit measure is either ELSE or MRSS. Then, for the score function $L(\tilde{Y}_{n}|X,\hat{\boldsymbol{\theta}}_{X})$ from Eq. (8) to infer the causal direction ${X\to Y}$, holds

[TABLE]

with equality if and only if the function $f$ in Eq. (3) is linear.

Proof: To prove the theorem, we need the following definition from [30] and our Lemma 1.

Definition 1**.**

(Identifiable Regression-based Scoring Functions) [30]. Given two random variables $X$ and $Y$ and a regression function $\varphi$ that maps $X$ to $Y$, and a scoring function $S:\mathbb{R}_{\geq 0}\times\mathbb{N}\to\mathbb{R}$ which takes as input the expected least-squared error $E[(Y-\varphi(X))^{2}]$ and the number of parameters of $\varphi$, $\|\beta_{\varphi}\|_{0}$. A scoring function

[TABLE]

is called an Identifiable Regression-based Scoring Function (IRSF), if both $\gamma:\mathbb{R}_{\geq 0}\to\mathbb{R}$ and $\lambda:\mathbb{N}\to\mathbb{R}$ are strictly increasing.

Lemma 1**.**

Let the causal model be given by Eq. (3) where $Y_{n}$ denotes the effect $Y$ estimated via regression using $n$ samples. Assume that Assumptions 2–4 hold for this model. Let $\varphi$ denote the function that minimizes the ELSE when predicting the effect $Y$ from the cause $X$, and let $\psi$ be the function minimizing the ELSE for predicting $X$ from $Y$, i.e. $\varphi(x)=E[Y\mid x]$ and $\psi(y)=E[X\mid y]$. Let $S$ be an IRSF as defined in Definition 1. Then the following limit always holds

[TABLE]

with equality if and only if $\varphi$ is linear.

Proof: We have moved the proof of Lemma 1 to Appendix -A due to its technical nature. We note that a similar lemma was proven as Theorem 1 in [30], although for the model given by Eq. (1) and $\alpha\to 0$. Regarding our Theorem 1, we will show that $L(\tilde{Y}_{n}|X,\hat{\boldsymbol{\theta}}_{X})$ satisfies the conditions of Lemma 1. In the nominator an denominator of Eq. (9) the parameterizations $\hat{\boldsymbol{\theta}}_{X},\hat{\boldsymbol{\theta}}_{Y}$ are fixed, and therefore do not affect the limit in Eq.  (9). The fulfillment of the IRSF conditions for the function $\lambda:\mathbb{N}\to\mathbb{R}$ is straightforward: Since $b(n)$ is strictly increasing in $n$, it follows that also $\log(b(n))$ is strictly increasing, thus it holds for $\lambda:\mathbb{N}\to\mathbb{R}$, too. Secondly we show the fulfillment of the conditions for the function $\gamma:\mathbb{R}_{\geq 0}\to\mathbb{R}$. The statement is obvious for the goodness of fit ELSE, since $\log$ is a strictly increasing function. For MRSS: In a discrete finite scenario, MRSS reflects empirical variations but aligns with ELSE for large $n$. As the number of data samples $n$ converges to infinity, MRSS converges to ELSE, assuming unbiased estimation, see [18]. Functional spaces with the density property on a compact interval provide unbiased estimation for a sufficiently large $n$. So for sufficiently large $n$ we can approximate

[TABLE]

[TABLE]

with $\hat{\boldsymbol{\theta}}_{X}$, $\hat{\boldsymbol{\theta}}_{Y}$ parameterizations of models corresponding to $X\to Y$ and $Y\to X$, respectively. Thus for function $\gamma:=\frac{n}{2}\log\frac{RSS(\hat{\boldsymbol{\theta}}_{X})}{n}$ holds $\gamma:\mathbb{R}_{\geq 0}\to\mathbb{R}$ and this function is strictly increasing for $n$. Therefore the MDL score is IRSF according to Definition 1 and the statement of Theorem 1 follows. ∎

In the following, we use our instantiation of the cubic spline functional class for the set $H$. This class has the density property on $C([0,1])$ and we construct the MDL-based causal score about which we show that it is IRSF. In this way, the identifiability of the score under the given mild conditions is proven.

V Cubic spline regression for bivariate causal discovery

V-A The class of cubic splines and its density property

Cubic spline functions have been well studied over decades as they provide reliable approximation of continuous functions. The class of cubic splines with variable number of knots is dense, i.e., any continuous function having fourth derivative on a compact interval can be approximated by a cubic spline within an arbitrary precision, having an arbitrary finite number of knots, [40]. However, the parametric representation of any continuous function by cubic spline does not have to be unique. The density property holds even more general for degree $r\geq 3$ as we sketch in the following. In general, on a compact interval, for the subclass of continuous functions, splines of best approximation do have uniqueness properties in the case when they are tied together from smooth functions to form a Chebyshev system on $[a,b]$ [43]. [42]. Theorem 3 in [42] shows that if one approximates a function of smoothness $r-1$ by splines $s$ of degree $r$ with variable knots and the number at most $m$, then at the expense of the choice of knots, the order of approximation, i.e. $\|f-s\|_{L_{p}[a,b]}$ is equal to $m^{-r-1+i}$ where $i\leq r-1$ is integer and $m$ is the number of knots and $f:[a,b]\to\mathbb{R}$ is a given continuous function. Based on the above results of Subbotin in [42], with the increasing number of knots $m$, a continuous function can be approximated by a cubic spline function within arbitrary precision. For interpolation on a finite set, spline interpolant exists and is unique for a given set of data points, provided the knots and interpolation points satisfy certain geometric conditions, see the Schoenberg-Whitney conditions [40].

Consider a curve to be a mapping associated by parameters $\beta_{j}$ with a fixed degree and fixed knots. A spline of degree $r$ with knots at $k_{j},j=1,\dots,m$ is a piecewise polynomial with continuous derivatives up to order $r-1$ at each knot. If the degree and the knots are fixed, then we have a vector space of piecewise polynomials. The B-spline functions form a basis for this vector space; Thus the coefficients $\beta_{j}$ of any given curve in this vector space are uniquely determined [6].

V-B Representation of $(X,Y)$ by a cubic regression spline

We will model the variable $Y$ as a cubic spline function of the variable $X$, both of them given by datasets and affected by independent noise. The spline is given by the dataset and by the parametrization as in Eq. (4) and we will approximate $f$ by a cubic regression spline with $m$ knots

[TABLE]

$k_{j}$ is the position of the j-th knot, $\{b_{0},b_{1},b_{2},,b_{3},\beta_{1},\dots,\beta_{m}\}$ is the set of coefficients and $(c)_{+}=\max(0,c)$. We further assume that $\min(x_{i})<k_{1}<\dots<k_{m}<\max(x_{i})$ and $\{k_{1},\dots,k_{m}\}\subset\{x_{1},\dots,x_{n}\}$. Denote $\bf{k}$ $=(k_{1},\dots,k_{m})^{\top}$, $\boldsymbol{b}=(b_{0},b_{1},b_{2},b_{3})^{\top}$ and $\boldsymbol{\beta}=(\beta_{1},\dots,\beta_{m})^{\top}$. If function $f$ can be approximated by the spline function $s_{f}$ from (13), then the estimate of $f$ can be obtained by estimating $m,\bf{k}$, $\bf{b}$ and $\bf{\boldsymbol{\beta}}$. Thus the problem of estimation of $f$ becomes a model selection problem. We restrict ourselves to cubic splines with equidistant knot sequences, since this class also has the density property in a compact interval, see e.g.[17].

V-C Cubic spline parametric approximation of a continuous function

In the following we will construct the MDL description of a cubic regression spline function. The problem of estimating the function $s_{f}$ in (13) can be seen as the model selection problem where the model is characterized by a parameter vector

[TABLE]

where $m$ the number of knots, $\boldsymbol{b},\boldsymbol{\beta}$ the vectors of coefficients and $\boldsymbol{k}$ the vector of $m$ knot points. Assume now the estimates $\hat{m},\hat{\boldsymbol{k}}$ are known. Denote vectors $\boldsymbol{x}=(x_{1},\dots,x_{n})^{\top}$, $\boldsymbol{y}=(y_{1},\dots,y_{n})^{\top}$ and design matrix

[TABLE]

where ${\bf 1}$ is $n\times 1$ vector of ones. Then if $m,\boldsymbol{k}$, are already specified, by using normal equations for multivariate regression we can get the estimates of $\boldsymbol{b},\boldsymbol{\beta}$ as

[TABLE]

where $A^{+}$ is the inverse or pseudo-inverse of matrix $A$. (Since knots are not fixed, we do not have to find matrix $\beta$ uniquely, therefore pseudoinverse.) Eq. (16) is the least squares estimate of $\bf{b},\boldsymbol{\beta}$ given estimates of $m$, $\boldsymbol{k}$. If for the error term holds $N\sim N(0,1)$ in Eq. (1), the estimate of $\bf{b},\boldsymbol{\beta}$ in Eq. (16) is also the maximum likelihood estimate.

V-D Causal rule by cubic spline approximation

We will compute our scoring function by MDL and denote it for both causal directions as $\hat{\Delta}({X\to Y})$ and $\hat{\Delta}({Y\to X})$. To be able to compute these scores, we need to compute $L(Y|X,\boldsymbol{\theta})$ and $L(X|Y,\boldsymbol{\theta})$. $L(Y|X,\boldsymbol{\theta}),L(X|Y,\boldsymbol{\theta})$ will be computed as the MDL two-part code for the cubic spline regression functions $\hat{s}_{f}$ and $\hat{s}_{g}$, where $f$ and $g$ are functions for each causal directions $X\to Y$ and $Y\to X$, respectively, as in Eq. (13). The code lengths of $L(Y|X,\boldsymbol{\theta}),L(X|Y,\boldsymbol{\theta})$ are explained in the following section.

V-E Code length $L(Y|X,\boldsymbol{\theta})$ and $L(X|Y,\boldsymbol{\theta})$

Now we will compute codes $L(Y|X,\boldsymbol{\theta})$ and $L(X|Y,\boldsymbol{\theta})$. With a slight abuse of notation and for simplicity we do not write $\tilde{Y}_{n}$ but $Y=\{y_{1},\dots,y_{n}\}$ for a cubic spline regression on set $X=\{x_{1},\dots,x_{n}\}$, i.e.

[TABLE]

where $s_{f}$ is defined by (13) with $m$ knots and analogously

[TABLE]

where $s_{g}$ a cubic spline function for the opposite direction. Without the lack of generality, we explain only $L(Y|X)$. The estimation of $Y$ by $s_{f}(X)$ is a model selection problem where each model is defined by parameter $\boldsymbol{\theta}$. Having parameter vector $\boldsymbol{\theta}=(m,\boldsymbol{k},\bf{b},\boldsymbol{\beta})$ from Eq. (14) corresponding to a fitted spline model, we compute the two-part MDL code as

[TABLE]

$L(\boldsymbol{\theta})$ can be still decomposed into

[TABLE]

and $L(Y|\boldsymbol{\theta})$ is a goodness of fit. We derive

[TABLE]

where $u_{j}$ are defined as follows: Let $k_{0}=\min(x_{i})$ and $k_{m+1}=\max(x_{i})$ and similarly as in [27] let

[TABLE]

for $j=1,\dots,m$. Each $u_{j}$ is the $j-th$ successive (i.e. next) index difference. If we know the whole sequence $(u_{1},\dots,u_{m})$, we also know the whole knot sequence $\boldsymbol{k}$. For fluency in reading we moved the full derivation of $L(\boldsymbol{\theta})$ to Appendix -B.

V-F Code length $L(Y|\boldsymbol{\theta})$

Now we code the goodness of fit, using the code from [39], as the code of $Y$ given $\boldsymbol{\theta}$ by

[TABLE]

where the residual sum of squares $RSS(\boldsymbol{\theta}):=\sum_{i=1}^{n}(y_{i}-s_{f}(x_{i}|\boldsymbol{\theta}))^{2}$ and C is a constant which is for all models the same and thus we can be omitted from MDL for the model selection problem.

V-G The full code $L(Y|X,\boldsymbol{\theta})$

The full code for conditionals is

[TABLE]

where $RSS(\boldsymbol{\theta})=\sum_{i=1}^{n}(y_{i}-s_{f}^{m,\boldsymbol{k}}(x_{i}))^{2}$ and $s_{f}^{m,\boldsymbol{k}}(x_{i})$ is a cubic regression spline with knot sequence $\boldsymbol{k}$ with $m$ knots. Algorithm 1 summarizes these computations.

To compute $\hat{\Delta}_{n}({X\to Y})$ in Algorithm 1, we take the maximum number of equidistant knots $m_{max}$ of the cubic regression splines as a hyperparameter. Based on the observation from [42] mentioned in Section V-A, the parameter $m$ controls the precision of the uniform approximation of the function $f$ by the spline $s_{f}$. Algorithm 2 determines the causal direction between the scores $\hat{\Delta}_{n}(X\to Y)$ and $\hat{\Delta}_{n}(Y\to X)$.

V-H Identifiability of the MDL cubic spline regression score

Corollary 1**.**

Assume that the causal direction $X\to Y$ for the dataset $(X,Y)$ with $n$ samples is modeled by Eq.  (3) and the score $\hat{\Delta}_{n}({X\to Y})$ was computed by Algorithm 1 for even $m$. Then

[TABLE]

with equality if and only if the function $h$ in Eq. (3) is linear.

Proof: We prove Corollary 1 by showing that $L(Y|X,\boldsymbol{\theta})$ code from Eq. (24) is IRSF according to Definition 1. Assumption 1 is satisfied, since both Eq. (2) and Eq. (3) are special cases of Eq. (1). Assumption 2 is satisfied for both models from Eq. (2) and Eq. (3) due to the density of the class of cubic splines in $C([0,1])$ with an arbitrary number of knots $m$. Assumption 3: Using spline regression, compact supports are assumed for numerical stability. Assumption 4: It holds that if we know that $\varphi$ from Definition 1 consists of a linear combination of basis functions that are linearly independent of each other, we cannot find an inverse function that has fewer degrees of freedom, see [30]. Concretely, Kilbertus et al. in [29] showed that for low degree polynomials, it is not possible to formulate an inverse with less parameters as the original function. We can therefore assume that $\|\beta_{\varphi}\|_{0}\leq\|\beta_{\psi}\|_{0}$, i.e. Assumption 4 is satisfied for the MDL cubic spline regression score. In Definition 1, if $\varphi$ is the function minimizing the expected least-squared error (ELSE) for predicting the effect $Y$ from the cause $X$, then in case of MDL cubic regression score for a fixed $m$, it holds $\varphi:=\hat{s}_{f}$ and if $\psi$ is the function minimizing the ELSE in the anti-causal direction, which is in case of MDL cubic spline regression score for a fixed $m$ then it holds $\psi:=\hat{s}_{g}$. We define

[TABLE]

The value $m$ depends on $n$. The product $m\times\prod_{j=1}^{m}u_{j}^{X}\times n^{\frac{m+4}{2}}$ has for even $m$ values in $\mathbb{N}$ and thus for $\lambda:=\log(...)$ holds that $\lambda:\mathbb{N}\to\mathbb{R}$; Function $\lambda$ is strictly increasing function for $n$. Analogically holds the same also for $\lambda(\|\beta_{\psi}\|_{0})$ for the direction $Y\to X$. The relationship between ELSE and MRSS holds also for the class of cubic splines, as this class is an instance of a dense functional class in Theorem 1. As the number of knots $m$ increases and the data sample size $n$ grows, due to the density property of cubic splines, the MRSS converges to ELSE. ∎

Our method for computation of causal score which is performed by Algorithm 1 and Algorithm 2 is called LCUBE, as an abbreviation for mdL CUBic splinE score.

Complexity of the LCUBE algorithm: The complexity of computing the regression spline $s_{f}^{m}$ in Algorithm 1 is $O(nm+m^{3})$. The complexity of LCUBE by Algorithm 1 for $m\leq m_{max}$ together with Algorithm 2 is $O(2m_{max}(nm_{max}+m_{max}^{3}))$. Further details on the computation are provided in Appendix -C. In our experiments, $m_{max}$ was at most 10.

V-I Differences of LCUBE to Slope and Sloppy

There are two methods, Slope [28] and Sloppy [30], which use spline regression, and therefore are related to LCUBE. Slope employs a two-part MDL to approximate the algorithmic Markov condition for continuous data. It also assumes that the error is Gaussian distributed. The score used in Slope can be written as $\gamma(E[(\tilde{Y}_{\alpha}-\varphi(X))^{2}])$,where $\gamma$ is based on the negative log likelihood, plus a function $\rho$ over the parameters. Since $\rho$ does not purely consider the number of parameters, but assigns different weights according to the value, of the parameter, the corresponding scoring function is not an IRSF. No identifiability results are known for Slope. Sloppy is is the first instantiation of the IRSF framework. Method Sloppy (Version SloppyS) fits a cubic spline, where the degrees of freedom is controlled and that selection of splines is found, for which S is minimal. As a scoring function $S$, AIC and BIC are used.

The scoring function of LCUBE is the $L(Y|X,\boldsymbol{\theta})$ score from Theorem 1 and it satisfies Assumptions 1-4 for cubic splines. Moreover, our MDL encoding allows more detailed characterization of cubic spline models than the encoding via AIC or BIC applied by Sloppy.

VI Approximation by dense functional classes and identifiability of the LCUBE score

VI-A The classes of cubic splines and harmonic functions and their density property

For some functional classes that possess the density property, the upper bounds on the approximation errors are known. In approximation theory, the rate of approximation on a compact interval describes how quickly an approximating function converges to the target function as the approximation parameters, such as polynomial degree or number of basis functions, improve. The density property of cubic splines was discussed in Section V-A. A well known upper bound from [13] states that the class of cubic splines approximates smooth functions with an error rate of $O(h^{4})$ where $h$ is the spacing between knots. In Fourier approximation, the approximation error depends on the function’s smoothness and decreases at a rate of order $O(1/n^{2})$ for differentiable functions, where $n$ is the number of terms included in the Fourier series expansion, see [16]. In area of neural network approximation, the rates of approximation refer to how efficiently neural networks can approximate functions. These are measured in terms of the number of neurons or layers required to achieve a certain level of accuracy. For an overview, see e.g. [15]. Regarding cubic splines with the equidistant knot sequence $\boldsymbol{k}$ and values $\boldsymbol{b},\boldsymbol{\beta}$, we assume that they satisfy the Schoenberg-Whitney conditions [40]. These ensure the existence of a unique cubic spline that interpolates the given data points. As we will see from the experiments with LCUBE in Section VIII, these properties of cubic splines contribute to its high causal score, particularly for datasets with low noise.

We showed above that the set of cubic splines satisfies Assumption 4. Intuitively, if $X$ causes $Y$, the functional causal model $Y=f(X)$ may exhibit a smoother relationship compared to the anti-causal model $X=g(Y)$. This could mean that the functional dependence of $Y$ on $X$ is more regular and smoother, potentially requiring fewer knots in spline approximation than for the functional dependence of $X$ on $Y.$ The mapping from $Y$ to $X$ may be more complex, e.g., nonlinear or multi-valued, thus requiring more spline knots to adequately capture its variations.

VI-B The class of shallow feed-forward neural networks and its density property

The class of feed-forward neural networks with one hidden layer and $d$ input units defined on $C([a,b]^{d})$ (i.e., shallow feed-forward networks), has the universal approximation property [19]. This property states that for any sufficiently smooth function $f$ on a compact set with finitely many discontinuities, there exists a feedforward network with one hidden layer (abbreviated as FNN), denoted $F$ that can approximate $f$ arbitrarily well if the number of hidden units $m$ in $F$ is sufficiently large and the activation function $\sigma$ is not polynomial. In other words, the FNN class is dense in $C([a,b]^{d})$. The literature provides upper bounds on the number of hidden units required to achieve a given approximation error on a compact interval for various activation functions, e.g. [24], [14]. These upper bounds can serve as an upper bound on $m$ while keeping the the approximation error in a compact interval under some value. However, in a finite sample scenarios, these bounds are too high in practice. However, directly constructing a FNN and estimating its approximation error on a finite dataset is challenging. Such a FNN would have the generalization ability without being trained. It is well known that the output of an FNN depends heavily on the chosen training algorithm and its parameter settings. In short, to find encodings that characterize the structure and learning of FNNs for both causal and anti-causal models, such that they satisfy Assumption 4, is a difficult task.

VII Related Work

We focus on state-of-the-art bivariate causal discovery methods for continuous data that are closely related to our work and highlight their identifiability results. The first identifiable methods are those based on the Additive Noise Model (ANM) [41], where $Y=f(X)+N_{X}$ and the noise $N_{X}$ and cause $X$ are independent. This model is identifiable when there does not exist an ANM in the anti-causal direction; This holds for instance for linear functions $f$ and non-Gaussian noise $N_{X}$ [41] or for nonlinear functions and additive noise [20] or for post-nonlinear models [36]. Although identifiability conditions are provided in these publications, they are often difficult to verify in practice. Furthermore, the results regarding causal direction heavily depend on the chosen independence test and the fitting algorithm [32].

The most related methods to our work are RECI from Blöbaum et al. [31], Slope [28] and Sloppy [30], the last two from Marx and Vreeken, as these approaches base their inference rules on the regression error. Methods based on regression error are able to decide between Markov equivalent DAGs under the assumption of having a non-linear function and low noise, see [31]. As mentioned earlier, Slope does not provide identifiability results. Method RECI does not employ model selection. A third method based on regression error is CAM from [1], which was designed to infer general causal graphs. For the bivariate case, CAM determines for the causal direction using regularized log-likelihood scores.

The IGCI method from Janzing et al. [26] was also proposed for low noise settings. The causal direction is inferred using the Shannon entropy of the marginals. However, entropy estimation for continuous data is challenging. Method QCCD [44] approximates the algorithmic Markov condition using non-parametric conditional quantile estimation. Although QCCD performs well in practice, its lacks strong identifiability guarantees. To the group of methods with low noise we count also methods which assume a homoscedastic noise. For homoscedastic noise models, regression methods are used as model estimators in many methods, see Shimizu et al.[41], Peters et al. [35], Bühlmann et al. [1]. Method CAM from [1] also belongs among them, using the maximum likelihood for inferring the cause and effect. Similarly, method RESIT from Peters et al. [35] executes the mean estimation, but an additional step of independence testing is done subsequently.

Methods which do not assume a low noise use heteroscedastic noise models, e.g. method HECI from Xu et al. [46]. HECI partitions the domain of the cause into multiple segments where the noise is assumed to be homoscedastic. An estimator for maximizing the Gaussian log-likelihood is presented in method LOCI [23]. This method utilizes the log-likelihood from the mean and the variance parameters of the Gaussian distribution as a criterion (version LOCIM) and the residual is recovered for the next step of testing the independence with Hilbert-Schmidt Independence Criterion (HSIC), in version LOCIH. The conditions on identifiability of these LOCI versions are given, however, similarly as for ANMs, they are difficult to test in practice. There are also two neural-network estimators of LOCI [23], these without identifiability guarantees. Method ROCHE [45], also without identifiability guarantees, follows the same procedure as LOCIH. ROCHE uses a neural network to estimate parameters of an objective function which fits a log-likelihood of Student t-distribution on a data. In general, the neural network estimators of causal direction have several advantages over the likelihood-based approach, since they can capture non-linear relationships between variables, and can be more robust to noise and outliers, leading to more reliable causal direction estimates. Being aware of this advantage, we discuss the comparison of LCUBE to the methods in this group separately.

VIII Experiments

In our experiments, we used these three groups of methods for comparison. The regression based methods are represented by Sloppy version for splines, since it is most similar method to LCUBE. We do not use Slope and RECI, since they had similar or worse results than Sloppy on all benchmark datasets we use (see their performance in [30]). Instead of RECI, we compare to RESIT [35]. We used both synthetic and real-world causal discovery benchmark datasets to evaluate LCUBE, and compared it to state-of-the-art methods. As an evaluation metrics we used accuracy (ACC) and AUDRC introduced in [23]. ACC measures the fraction of correctly inferred cause-effect relationships and the AUDRC measures the area under the decision rate curve. [23] justifies the selection of AUDRC over AUROC, since it weights correctly identified $X\to Y$ pairs in the same way as correctly identified $Y\to X$ pairs and thus avoids an arbitrary selection of true positives and true negatives. (More in Appendix -D1.) For Tübingen pairs, for which a special weighting of each pair is given, we used the so called forced decision from [32]. Our code, its description and parameter setting are publicly available under https://github.com/suzi216/LCUBE.

VIII-A Benchmark datasets

We used the same benchmark datasets as in recent papers [45] and [23]. The first subroup of synthetic datasets are from [44] including AN, AN-s, LS, LS-s, and MN-U. They consist of nonlinear functions with additive noise (AN), sigmoidal functions with additive noise (AN-s), nonlinear and sigmoidal location scale functions (LS and LS-s), where $Y=f(X)+g(X).N_{Y}$ and sigmoid functions with multiplicative uniform noise (MN-U) with the effect Y is the form $Y=f(X).N_{Y}$. The datasets denoted with “-s” and the MN-U dataset use invertible sigmoid-type functions for the generative functions which can make the identification of the causal direction more complicated. All these simulated datasets consist of 100 cause-effect pairs with 1 000 samples per pair. Further we consider datasets SIM, SIM-c, SIM-ln, and SIM-G from [32], each having 300 pairs, with 500 samples per pair. These datasets are simulated with the functions $X=f_{X}(\varepsilon_{X})$ and $Y=f_{Y}(X,\varepsilon_{Y})$ in the cases without a confounder (in SIM, SIM-ln, and SIM-G) and the functions $Z=f_{Z}(\varepsilon_{Z})$, $X=f_{X}(\varepsilon_{X},\varepsilon_{Z})$, and $Y=f_{Y}(X,\varepsilon_{Y},\varepsilon_{Z})$ in the cases with a confounder Z (in SIM-c). In the SIM-ln dataset, low levels of noise are applied in the models, and the SIM-G dataset has approximations of Gaussian distributions for the cause $X$ and approximately Gaussian non-linear additive noise generative models. The third group of synthetic datasets are Multi, Net from [9], and Cha from [12]. The data pairs in Multi dataset are generated with pre-additive noise as $Y=f(X+\varepsilon_{Y})$, post-additive noise (similar to conventional additive noise models), pre-multiplicative noise $Y=f(X\varepsilon_{Y})$, and post-multiplicative noise $Y=f(X)\varepsilon_{Y}$. The pairs in the Net dataset are generated by neural networks with random distribution for $X$, such as exponential, gamma, log-normal, or Laplace distribution. The data pairs in the Cha benchmark are chosen from the ChaLearn Cause-Effect Pairs Challenge [12]. All three datasets have 300 continuous variable pairs with 500 samples per pair. And finally, we used the Tübingen dataset, consisting of 99 real-world causal pairs [32]. The sample size of Tübingen pairs varies, ranging from a few hundred to five thousand.

VIII-B Performance and comparison with the state-of-art methods

We evaluate LCUBE with a focus on benchmark methods:

Methods with similar assumptions, i.e. regression with low noise - Sloppy, CAM, QCCD, RESIT and IGCI;

• 2.

Methods using regression with heteroscedastic noise - LOCIM, LOCIH and HECI;

• 3.

Methods using regression with heteroscedastic noise and a neural network estimator - NN-LOCIM, NN-LOCIH and ROCHE.

The ACC and AUDRC values for CAM, HECI, QCCD, RESIT, IGCI, and LOCI are adopted from [23] (Table 1 and 2). We computed ACC and AUDRC values for ROCHE using results from their Github repository [45]. For Sloppy, we used the provided R code and ran it on the Cha, Net, and Multi datasets. The ACC and AUDRC values for all methods are reported in our Table II and Table III in Appendix -E, where also the performance of all methods on all 13 datasets is reported. The summary of the AUDRC values for LCUBE and a representative method from each category of benchmark methods is shown in Figure 1. The precisions in AUDRC and ACC for each method in average over 10 and 13 datasets are in Table I.

Regarding Figure 1, no method performs perfectly across all datasets. However, LCUBE achieves the highest AUDRC among the four methods on the AN, ANs, LS, LSs, MN-U, and Tübingen pairs datasets. AN, ANs, LS, LSs, MN-U have a small noise (not necessarily Gaussian). One can see that also Sloppy, LOCIH and ROCHE achieve a high precision on the Gaussian AN, ANs, LS. However, AUDRC values of Sloppy and LOCIH drop on LSs due to the higher complexity of this dataset. Sloppy’s AUDRC is lower than LCUBE’s on 12 datasets and it is higher only on set Multi. LOCIH ad ROCHE show overall high AUDRC on synthetic datasets, but significantly lower AUDRC compared to LCUBE and Sloppy on real Tübingen data pairs.

LCUBE achieves perfect AUDRC and accuracy on the AN, ANs, LS, LSs and MNU datasets. In general, the precision of the methods depends on the precision of two main steps: the precision of an estimator for the model and the precision of the causal score method. LCUBE allows very precise estimators due to the density of the set of cubic splines. Interestingly, the NN method ROCHE demonstrates high precision across all synthetic sets, except for MN-U, where it has ACC of 62 and AUDRC of 65. In contrast, LCUBE has perfect precision (100) on MN-U. While all NN methods perform well on synthetic datasets, they lack identifiability guarantees. They also perform poorly on real data of Tübingen pairs, with AUDRC of 56 for NN-LOCIH and 66 (NN-LOCIH), and 74 (ROCHE), while AUDRC of LCUBE achieves score 87. To the best of our knowledge, this is the highest score achieved by any state-of-the-art method. As stated before, LCUBE is a scoring method with identifiability guarantees.

LCUBE has only one hyper-parameter: the upper bound for number of knots $m_{max}$, which is also used by Sloppy. In contrast, Roche, NN-LOCIH and NN-LOCIM require multiple hyperparameters for their NN architectures and Adam optimizer. Methods that use HSIC score require additional hyperparameters, including the kernel type, bandwidth and regularization parameter. To accelerate the execution of the computational models, we parallelized both the code and the data. We used the NVIDIA RTX 4090 GPU and 64 GB of RAM. Regarding the real running time on our hardware, LCUBE is significantly faster than both Sloppy and ROCHE. For example, on the Tübingen pairs dataset, LCUBE takes approximately 3 minutes, while Sloppy takes 20 minutes and ROCHE around 6 hours. On the Cha, Net or Multi datasets, LCUBE takes 5 minutes and 2 minutes on the other datasets. In comparison, Sloppy takes 15 min on Cha, Net, Multi, and 10 min on AN, MN-U, SIM. Roche requires around 8 hours on Cha, Net, and Multi and 5 hours on AN, MN-U, and SIM.

IX Conclusions

In this work, we proposed a bivariate causal score based on the the Minimum Description Length (MDL) principle, using functions with the density property and proved the identifiability of these scores under specific conditions. These conditions can be easily tested. Gaussianity of the noise in the causal model equations is not assumed, only that the noise is low. We also introduced LCUBE as an instantiation of the MDL-based causal score utilizing well-studied cubic regression splines. LCUBE is identifiable, interpretable, simple and very fast. LCUBE has only one hyperparameter. It achieves superior precision in terms of AUDRC on the Tübingen cause-effect pairs dataset compared to the state-of-the-art methods. It also shows superior average precision across common 10 benchmark datasets and above average precision on 13 datasets. In future work, we aim to explore causal scores via other dense functional classes.

**Acknowledgements

**We thank Dr. Negar Safinia naini and Dr. Amir Rahnama for valuable discussions in the early phase of the paper and Prof. Jilles Vreeken for his comments in the final phase.

-A Proof of Lemma 1:

Under Assumptions 1-3 from Blöbaum et al. in [31], it holds

[TABLE]

with equality only if the $\phi$ is linear. (Note: The fraction on the left side of Eq. (26) is in Theorem 1 in [31] is formulated as $\lim_{\alpha\to 0}\frac{\mathbb{E}[(X-\psi(\tilde{Y}_{\alpha})^{2}]}{\mathbb{E}[(\tilde{Y}_{\alpha}-\phi(X)^{2}]}\geq 1$ and thus this inequality is equivalent to Eq.  (26).) Since Eq. (26) holds for any $\alpha\to 0$, it holds also for model given by Eq. (3) and thus we can write

[TABLE]

Similarly as in proof of Theorem 1 in [30], as $S$ is an IRSF, we can write it as $S(a,b):=\gamma(a)+\lambda(b)$ where $\gamma$ is a strictly increasing function. The statement does not change by applying $\gamma$ to the numerator and denominator in Eq. (27).

Based on Assumption 4, i.e. $\|\beta_{\varphi}\|_{0}\leq\|\beta_{\psi}\|_{0}$ holds and thus

[TABLE]

with equality if and only if $\|\beta_{\varphi}\|_{0}=\|\beta_{\psi}\|_{0}$. As $\lambda$ is strictly increasing, applying it to $\|\beta_{\varphi}\|_{0}\text{ and }\|\beta_{\psi}\|_{0}$ will not change this statement.

-B Code length $L(\boldsymbol{\theta})$ for a cubic regression spline

For the moment we select a fixed number of (internal) knots $m=|\boldsymbol{k}|$ where $\boldsymbol{k}=\{k_{1},\dots,k_{m}\}\subset X$.

-B1 Code $L(m)$

$L(m)$ can be approximated by $\log m$ when $m$ is reasonably large, see paper [38], Section 2.2.4. Thus based on this, $L(m)\approx\log_{2}m.$

-B2 Code $L(\boldsymbol{k}|m)$

First we will compute the code $L(\boldsymbol{k}|m)$. Since $\boldsymbol{k}=\{k_{1},\dots,k_{m}\}\subset X$, the sequence of indices $\boldsymbol{k}$ can be specified by the indices of those $x_{i}$’s where a knot is placed. This set of sorted indices can be compactly specified by their successive differences. Define $k_{0}=\min(x_{i})$ and $k_{m+1}=\max(x_{i})$ and similarly as in [27] let

[TABLE]

for $j=1,\dots,m$. Each $u_{j}$ is the $j-th$ successive (i.e. next) index difference. If we know the whole sequence $(u_{1},\dots,u_{m})$, we also know the whole knot sequence $\boldsymbol{k}$. If the knots $\boldsymbol{k}$ satisfy the Schoenberg-Whitney conditions, is each $u_{j}$ a positive integer and the sequence $\boldsymbol{k}$ can be coded by

[TABLE]

where $L(u_{j})=\log u_{j}$. Since $u_{j}$ are integer values, holds

[TABLE]

-B3 Code for $L(\boldsymbol{b},\boldsymbol{\beta}|m,\boldsymbol{k})$

Now we compute $L(\boldsymbol{b},\boldsymbol{\beta}|m,\boldsymbol{k})$. Given $m$, $\boldsymbol{k}$, the values $\boldsymbol{b},\boldsymbol{\beta}$ can be computed from Eq. (16) and the result is the conditional maximum likelihood estimates of $\boldsymbol{b}$ and $\boldsymbol{\beta}$ (or the least square estimates if the assumption of normal errors is not given). Based on [38], pp. 55-56, if a conditional max. lik. estimate is estimated from $n$ points, then it can be efficiently encoded with $\frac{1}{2}\log n$. Since $r=3$, one can see that each $b_{j}^{\prime}s$ and $\beta_{j}^{\prime}s$ is estimated from all $n$ x-values. Hence

[TABLE]

Since a cubic spline is defined by $4$ coefficients then $L(\boldsymbol{b},\boldsymbol{\beta}|m,\boldsymbol{k})=\frac{m+4}{2}\log n.$ To summarize,

[TABLE]

-C Complexity of LCUBE

The placement of $m$ knots has a comp. cost of $O(1)$. Constructing the spline basis in Eq. (13) requires $O(m)$ operations. Solving the least squares for $n$ samples and $m$ basis functions, which involves solving an $m\times m$ linear system with $n\times m$ design matrix, results to total costs $O(nm+m^{3})$.

-D Evaluation measures

We use evaluation metrics accuracy and AUDRC from [23]. For real data of Tübingen pairs, we use the so called forced decision as defined in [32].

-D1 AUDRC

The area under the decision rate curve (AUDRC). [23] justifies the selection of AUDRC over AUROC, since it weights correctly identified $X\to Y$ pairs in the same way as correctly identified $Y\to X$ pairs and thus avoids an arbitrary selection of true positives and true negatives, which can lead to non-interpretable results on unbalanced datasets, see [28], Table 1. The accuracy (ACC) measures the fraction of correctly inferred cause-effect relationships and the AUDRC measures how well the decision certainty indicates accuracy. The certainty is, for example, indicated by the likelihood or $p$-value difference in both directions. Thus, a high AUDRC indicates that an estimator tends to be correct when it is certain and only incorrect when it is uncertain. Given the ground truth direction $t(.)$ and a causal estimator for the direction $f(.)$ on $Q$ pairs with ordering $\pi(.)$ according to the estimator’s certainty,

[TABLE]

the indicator function $\mathbf{1}_{A}(x)$ is defined as 1 if $x\in A$ otherwise 0. I.e., we average the accuracy when iteratively adding the pair the estimator is the most certain about. Higher AUDRC means that a method predicts the correct causal directions if the confidence of the scores is high.

-D2 Forced decision: Evaluation of accuracy ACC

Given $(X,Y)$, the method has to decide either $X\to Y$ or $Y\to X$; We evaluate the accuracy of these decisions as

[TABLE]

where $d_{q}$ is the true causal direction for the q’th Tübingen pair (either $\leftarrow$ or $\rightarrow$), and $\hat{d}_{q}$ is the estimated direction (i.e. $\leftarrow$ or $\rightarrow$ or ? = undecided) and $w_{q}$ is the weight of the pair, see [32].

-E Performance of all methods

AUDCRC and accuracy in $\%$ for each method and dataset can be found in Table II and Table III.