Maximal Correlation-Based Post-Nonlinear Learning for Bivariate Causal Discovery
Authors: Tianjian Zhang
Feng Yin
Zhi-Quan Luo
Venue: Submitted to ICLR 2023 (2023)

Abstract: Bivariate causal discovery aims to determine the causal relationship between two random variables from passive observational data (as intervention is not affordable in many scientific fields), which is considered fundamental and challenging. Designing algorithms based on the post-nonlinear (PNL) model has aroused much attention for its generality. However, the state-of-the-art (SOTA) PNL-based algorithms involve highly non-convex objectives for neural network training, which are time-consuming and unable to produce meaningful solutions with finite samples. In this paper, we propose a novel method that incorporates maximal correlation into the PNL model learning (short as MC-PNL) such that the underlying nonlinearities can be accurately recovered. Owing to the benign structure of our objective function when modeling the nonlinearities with linear combinations of random Fourier features, the target optimization problem can be solved rather efficiently and rapidly via the block coordinate descent. We also compare the MC-PNL with SOTA methods on the downstream synthetic and real causal discovery tasks to show its superiority in time and accuracy. Our code is available at https://anonymous.4open.science/r/MC-PNL-E446/.

Under review as a conference paper at ICLR 2023
MAXIMAL CORRELATION-BASED POST-NONLINEAR LEARNING FOR BIVARIATE CAUSAL DISCOVERY
Anonymous authors Paper under double-blind review
ABSTRACT
Bivariate causal discovery aims to determine the causal relationship between two random variables from passive observational data (as intervention is not affordable in many scientific fields), which is considered fundamental and challenging. De- signing algorithms based on the post-nonlinear (PNL) model has aroused much attention for its generality. However, the state-of-the-art (SOTA) PNL-bbssed al- gorithms involve highly non-convex objectives due to the use of neural networks and non-convex losses, thus optimizing such objectives is often time-consuming and unable to produce meaningful solutions with finite samples. In this paper, we propose a novel method that incorporates maximal correlation into the PNL model learning (short as MC-PNL) such that the underlying nonlinearities can be accurately recovered. Owing to the benign structure of our objective function, when modeling the nonlinearities with linear combinations of random Fourier features, the target optimization problem can be solved rather efficiently and rapidly via the block coordinate descent. We also compare the MC-PNL with SOTA methods on the downstream synthetic and real causal discovery tasks to show its superiority in time and accuracy. Our code is available at httts::/anooyyoos.4opppe.. Atp.eeepneeccnce//MCCCCCC PNL.FE466/
INTRODUCTION AND RELATED WORKS
Causal discovery is an old and new topic to the machine learning community, which aims to find causal relationships among variables. Many recent attempts at application have emerged in various scientific domains, such as climate science JEbert-Uphoti & Deng 2012: Runge et al. 2019). bioinformatics (Choi al. 2020, Foraita al. 2020) Shen et al.) 2020). etc. The gold standard for causal discovery is to conduct randomized experiments (via interventions), however, interventions are often expensive, unethical, and impractical. It highly demanded to discover causal relationships purely from passive observational data. In the past three decades, many pioneer algorithms for directed acyelic graph (DAG) searching have been developed for multi-variate causal discovery to reduce the computational complexity and improve the accuracy. For example, there constraint/ndeppeddncee based algorithms such as IC, PC, FCI (Pearl 2009 Spirtes et al. 2000), RFCI (Colombo et al.] 2012) (too many to be listed), as well score-based methods such as GES (Chickering) [2002). NOTEARS (Zheng et al., 2018), etc. However, the algorithms mentioned above can merely return Markov equivalence class (MEC) that encodes the same set of conditional independencies, with many undetermined edge directions; moreover, the discovered DAG may not necessarily be causal. In this paper, we will focus on a fundamental problem, namely bivariate causal discovery, which aims to determine the causal direction between two random variables X and Y. Bivariate causal discovery one promising routine for further identification of the underlying causal DAG (Peters et al., (017).
Bivariate causal discovery is a challenging task, which cannot be directly solved using the existing methodologies for the multivariate case, as the two candidate DAGs, X -> Y and X <- Y, are in the same MEC. More assumptions should be imposed to make bivariate causal discovery feasible, as summarized by Peters et al. (0117). One assumption is on the a priori model class restriction, e.g., linear non-Gaussiaa acyclic model (LiNGAM) (Shimizu et al. (2006), nonlinear additive noise model (ANM) (Mooij et al.) 2016). post-nonlinear (PNL) model (Zhang & Hyvärinen] 2009), etc. the "independence of cause and mechanism" leading the algorithms of trace condition (Janzing et al., 20102,1)) IGCI (Janzing et al., 2012), distance correlations (Liu & Chan 2016), meta-transfer (Bengio al. 2020), CDCI (Duong & Nguyen 2022), etc. There
The other assumption

---

Under review as a conference paper at ICLR 2023
also seminal works focusing on causal discovery in linear/nonlinear dynamic systems, which out of the scope of this paper, and the corresponding representatives are Granger causality test (Granger 1969) and convergent cross mapping (Sugihara et al. 2012 Ye et al. 2015).
In this work, we focus on the PNL model, which is more general than LINGAM and ANM. The existing works merely show the identifiability results with infinite data samples (i.e. known joint distribution), while practical issues with finite sample size are seldom discussed. We reveal the diffi- culties with the current PNL-based algorithms in the finite sample regime, such as insufficient model fitting, slow training progress, and unsatisfactory independent test performance, and correspondingly propose novel and practical solutions.
The main contributions of this work are as follows.
We point out various practical training issues with the existing PNL model learning algo- rithms, in particular PNL-MLP and AbPNL., and propose a new algorithm called MC-PNL (specifically the maximal correlation- baaeed algorithm with independence regularization), which can achieve a better recovery of the underlying nonlinear transformations.
2. We suggest using the randomized dependence coefficient (RDC) instead of the Hilbert- Schmidt independence criterion (HSIC) for the independent test and give a universal view of some widely used dependence measures.
3. We use MC-PNL for model learning in bivariate causal discovery and show that our method outperforms other SOTA independence test-based methods on various benchmark datasets.
2
PRELIMINARIES
In this section, we will introduce the HSIC as a dependence measure, the current HSIC-based causal discovery methods for PNL model, and other relevant learning methods based on the Hirschfeld- Gébelein-Rényi (HGR) correlation. Our proposed MC-PNL method exploits all these ingredients.
2.1 HSIC SCORE AND HSIC-BASHE REGRESSION
Regression by dependence minimization MMoijj et al. 2009) has attracted lots of attention recently. Greenfeld & Shalit (2020) has shown its power for robust learning, in particular the unsupervised covariate shift task. Let us consider the following regression model,
(1) where the additive noise is independent (symbolized by ILI with the input variable X, and the selected regression model fo is to be learned via minimizing the dependency between the input variable X and the residual Y - (X)) A widely used dependence measure is the Hilbert-Schmidt independence criterion (HSIC) (Gretton et al. 2005 (007).
Y = f(X) + e, Il X,
Definition 1 (HSIC). Let X, Z Pxz be jointly distributed random variables. and F, G be the reproduced kernel Hilbert spaces with kernel functions is and I, the HSIC can be expressed as,
HSIC(X, (;U F,),) =ExzEx (2, x I (z, x²) - ExExx k (x, x') EzlEzi (2,
2 Exc (z, x²) Ezi (2,
where x' and z'' denote independent copies of x and 2, respectively.
Remark 2.1. We can conclude that: (a) X 11 7 = HSTC O: (b) with a proper universal kernel (e.g. Gaussian kernel), X IL + HSIC(X, Z) Gretton et al. 2005
When the joint distribution Pxz is unknown, given dataset with samples (x [01, .2...., xx E RR" = [21, 12,22, 22,..., 2.1TT E R²), biased HSIC estimate can be constructed as,
HSTC (x, 2F F.9
tr(KHHH),
n²
where Ki,; = k (z), Lij = 1 (zi, zj). and Hl = -117 E x is centering matrix. The Gaussian kernel ls (xi, exp (-(x:; 20-) is commonly used, and the same for 1. One
2

---

Under review as a conference paper at ICLR 2023
can intuitively interpret this empirical HSIC the inner- inner-rrouct of two centralized kernel matrices (HKH, HLH), where the kernel matrices summarize the sample similarities.
Mooij al. (2009) first proposed to the above defined empirical HSIC for model learning. Con- cretely, the regression model a linear combination of the basis functions, fa( > (2))(x). and the parameters are learned from:
A 11112 i E arg min H fo(x)) + 2 BERP
(4)
where fe is applied elementwisely to the data points, and À >> 0 is a penalty parameter (we will keep using penalty parameter under different contexts). One key advantage of this formulation that it requires no assumption on the noise distribution. (Greenfeld & Shalit (2020) implemented fo using neural networks, and showed the learnability of the HSIC loss theoretically.
2.2 CAUSAL DISCOVERY WITH DOST-NONLI MODEL
The bi-variate post-nonlinear model is expressed as, Y (X) + <), where fi denotes the nonlinear effect of the cause, é is the independent noise, and f2 denotes the invertible post-nonlinear distortion from the sensor or measurement side. The goal is to find the causal direction X - Y from a set of passive observations on X and Y. Note that from the data generating process, is independent with X but not Y. Taking this asymmetry
prior, one test the causal direction by first learning the underlying transformations, fz" and f1, and then checking the independence between the residual ((-) = 12 (V) - (i(X) and the input X.
The PNL-MLP algorithm proposed by [Zhang & Hyvärinen (2009) tests between hypotheses (X + Y and X + Y) as follows. Under the hypothesis X y, one can parameterize fi and fz by two multi-layer perceptrons (MLPs) f(-+) and 9(-)). and learn them via minimizing the mutual information (MI):
( E arg min MI = - f(-))x);xx,
(5)
ft 109
where 9(-)) f(-) are applied elementwisely. The estimated residual is f(-)) = (y) (x)) Similarly, under the hypothesis X +- Y, one obtain estimate of f(()) ((-)(x) - (x-)
via minimizing MI(r(-)iy). The causal direction is determined by comparing HSIC (f(-)): x) and HSIC (f) If HSIC (fz-). < HSIC (f(+).)). the hypothesis X + Y is endorsed; otherwise, the hypothesis X Y is endorsed.
However, the MI between random variables is often difficult to calculate (see supplement Aj. and tuning the MLPs requires many tricks mentioned in Zhang Hyvärinen (2009), altogether bringing huge difficulties to handle large-scale datasets with many variable pairs. Uemura & Shimizu (2020)
proposed AbPNL method that uses HSIC instead MI, and imposes the invertibility restriction fz via auto unto-enccoer to eliminate nonsense solutions,
min HSIC (x g(y) f(x)) A - g'(o())) F.g.9'
where 9.g' g are encoder and decoder MIPss. The subscript (+) is omitted for conciseness here.
We summarize the architectures of the above-mentioned two methods in Figure I inherent issues exist concerning the cost function and the neural network training procedure when dealing with finite sample datasets, see in Section 3.1
Nevertheless,
2.3 PNL LEARNING THROUGH MAXIMAL CORRELATION
Another routine to learn the nonlinear transformations ff and g is through the HGR maximal correlation
(Hirschfeld 1935 Gebelein 194] Rényi 1959).
Definition 2 (HGR maximal correlation). Let X, Y be jointly distributed random variables. Then,
3

---

Under review as a conference paper at ICLR 2023
HGR(X;Y) =
sup F.x-R.gy= g -[(
([
is the HGR maximal correlation between X and Y, and f, 9 the associated maximal correlation functions.
Remark 2.2. The HGR maximal correlation p² is attractive as a measure of dependency due to some 11; (2) X and Y are independent if and only if p* 0; useful properties: (1) Bounded p² :0 p² 1;
(3) there exists f and g such that f(X) = g(Y) with probability if and only if pf =1
The optimal unit-variance feature transformations, f and g², can be found by iteratively updating and g in 77. However, for causal discovery applications, one fatal issue is that the learned f and 9 are constrained to have unit-variance, thus being unable to reflect the true magnitudes of the underlying functions f and g- As consequence, the resulting residual can be incorrect for the independence tests in the next stage. We found two possible remedies in the literature, namely the alternating conditional expectation (ACE) algorithm (Breiman & Friedman 1985) and a soft version of (7) (Soft-HGR) ){ (Wang et al. 2019).
The ACE algorithm solves the regression problem (8) by computing the conditional mean alternatively,
min (f g(V)) s.t. [[f(X) = E[g(Y)) = 0, [Y)] 1, fi
which only retains the unit-variance constraint on g. The equivalence 7 was established, and the regression optimal transformation
equals (p" f"; (",g") :,),, see Theorem 5. in Breiman & Friedman 1985). The other formulation, Soft-HGR, relaxes the unit-variance constraints as follows,
max EE [ (F())((Y) var(f(X)) var(g(Y)), s.t. f.g
0.
allows certain linear transformations (af" a "=") where a E R\(0) can produce infinitely many equivalent local minima. This scale ambiguity results in enormous useless solutions for causal discovery, and the desired one should make the estimated residual independent with the input. We will show how our proposed method is able eliminate those undesired solutions in Section 4
We notice that the recent proposed Variance-Invariancee Coviaannn-o Regu- larization (VICReg) (Bardes al. 2022) shares similar intuitions with the HGR maximal correlation. When the dimension of representation vectors (i.e., f and g) reduces to one, the covariance term disappears, and the VICReg objective becomes,
Connections to VICReg,
min ( f.g
[max(0, - var(f(x)) + (0, var (ar(g(()) (g())))
(10)
invariance term
variance term
where X, 7 > O are the paramm-srraeetrss that need to be tuned. The invariance term encourages the alignment of the learned features; and the variance term encourages 7-bbunded variation to avoid trivial solutions like f(X) = g(Y) constant. To see the connections, we rewrite Soft-HGR (9 as,
min [(()) o( Y)) + var var(s(Y)) var(s(X)) vYr(g(Y)),
(1I)
invariance [[] EYg( = 0,
s.t.
in which the variance of ff and g are also encouraged but not allowed to grow simultancously.
PRACTICAL ISSUES WITH EXISTING ALGORITHMS
In this section, we summarize several practical issues of the existing algorithms for PNL learning. including among others PNL-MLP (Zhang & Hyvärinen (010)) and AbPNL (Uemura & Shimizul 2020). These issues motivate our novel MC-PNL method to be introduced Section 4 see the comparisons in terms of their architectures in Figure
4

---

Under review as a conference paper at ICLR 2023
3.1 ISSUES ON MODEL LEARNING
Over-fitting issue. The general idea of the PNL model learning, according to Section 2.2 is to encourage the independence between the input and the estimated residual. Both PNI--MLP and AbPNL. use neural networks to parameterize f and g. But it is skeptical that meaningful representations can really be learned with finite samples. Let us review the dependence minimization problem below,
1 tr Len where r = f(x) -(y)). nT2
min HSIC(X. Fig
(12)
We argue that it is utmost difficult learn meaningful representations of and g via minimizing solely the HSIC score, due to the enormous degrees of freedom for f and g to learn arbitrary random noise. We conducted experiments using both wide ever-parameterized and narrow deep neural networks with sufficient representation power. In our simulation results (see supplement [). for both network architectures, the objective values can reach zero but unfortunately produce meaningless estimates. This is unsurprising though, as one can force to be samples from arbitrary independent random noise (Yun et al.) 2019; Zhang et al. 202I). To aid with that, we propose to cooperate dependence minimization with maximal correlation, which helps to obtain desired solutions, see Figure (e) for illustration and Section for details.
Optimization issue. The optimization of neural networks is long-standing problem, and yet there is not any study on the optimization landscape of the HSIC loss with neural networks. Typically, first-order methods such as stochastic gradient descent used in the existing causal discovery methods, and initialization is crucial to the causal discovery accuracy, see in supplement [ this paper, suggest parameterizing both and g linear combination of random Fourier features and using a linear kernel for HSIC, which admit a benign landscape with symmetry (see Chapter 7 in Wright & Ma (002)) for the non-convex optimization.
X
X
X
X
X
USTTC(: Res)
DISICI Res)
min MI(X;Ra Ros)
Y-
9
Res
Y-
Ros
Y
Res
ree(Y,y" y Y
corr (f (XX) (x)))
+A +X-mo(YY)")
min AUSSCC comm (XX))
(a) PNL-MLP
(b) AbPNL
(c) (MCrs)
Figure 1: Architectures of PNL learning frameworks.
3.2 ISSUES ON INDEPENDENCE TEST
As the independence test is critical to the accuracy of causal discovery, we have to cautiously choose the dependence measure. A Although HSIC is widely used, there several drawbacks of HSIC (e.g., the choice of kernel and corresponding hyper-parameter are user-defined, the values of HSIC depends on the scale of the random variables). In this section, we show experimentally that the HSIC score is not the best choice, and we favor randomized dependence coefficient (RDC) (opez--Paz al. 2013) particularly for finite samples.
We generated various synthetic datasets following the PNL models, see supplement C in which we know in advance that the injected noise IL X and _1 Y. Thus, we are able compare various dependence measures, by checking whether Dep(x, E) < Dep(y, E) on various datasets. In this section, the compared dependence measures are HSIC (Gretton al. 2005), its normalized variant (NOCCO) (Fukumizu et al.) (007). and RDC (Lopez-Paz et al. 2013). Besides, we also study the impact of different choices of linear, Gaussian radial basis function (RBF), and rational quadratic (RQ) kernels. We note here that RDC is computational tractable estimator inspired by the HGR maximal correlation. It shows that RDC outperforms other dependence measures especially when the sample size is small, see Table I Thus, advocate to use RDC to measure dependency with finite samples, Finally, we give a universal view the aforementioned dependence measures in supplement D
5

---

Under review as a conference paper at ICLR 2023
Table 1: The independence test accuracy (%) with known injected noise
samples noise level de
1000
n
1
10
0.1
10
IISIC-linear
88
100
85
100
HSIC-RBF
78
85
99
100
86
97
100
HSIC-RQ
79
100
86
99
100
91
96
100
NOCCO-RBF
87
83
75
83
85
72
90
92
90
NOCCO-RQ
87
82
73
81
86
65
85
94
87
RDC
94
100
100
98
100
100
100
100
100
PROPOSED METHOD
In this section, we propose maximal correlation- post-nonlinear model learning framework, called MC-PNL, to accurately estimate the nonlinear functions and compute the corresponding residuals. After then, independence tests will be conducted determine the causal direction.
4.1 MAXIMAL CORRELATION- PNL MODEL LEARNING
As we can see in the previous sections, minimizing HSIC T2 requires no assumption on the noise distribution and encourages the independence of the residual, but it can easily get stuck at meaningless local minima. Maximal HGR correlation based methods can learn meaningful transformations as its name suggested, not necessarily produce independent residual.
To combine their strengths, we propose the following MC-PNL objective, (X + var(f(X)) vY)) + f(X) Bb,
min f-9
(13)
s.t.
where Dep(.,-) 2 is dependence measure (e.g., HSIC with different kernel functions), and > hypeer-parameter that penalizes the dependence between the input variable X and the estimated residuals f(X) g(Y). This novel objective can learn meaningful feature transformations with the Soft-HGR term, and resolve the scale ambiguity via the dependence minimization term. The objective 1 is consistent with minimizing MI principle, under the assumptions of invertible PNL generating functions and Gaussian noise, see details in supplement
Parameterization with Random Features
For ease of optimization, we parameterize the transformation functions as the linear combination of the random features, namely fx; ai) a (xx) and g(y; (3) := 33 <b(v)) where the random features (xx) E 1 (b) E k2 are nonlinear projections as described in López-Paz et al. 2013)).
see supplement F For given dataset {(x;, 3/4 H the corresponding feature matrices are denoted + : (p(x1) ), &
(xz2).... (fxx)) E IR x and V (v()), (8y))....... (by)) E R2 We
further denote the residual vector as r == T T B.
a
Consequently, 13) can written as the following non-convex programming problem, - + 2. o कH a V B ADep(x,
min J(a, - T 0.33
(14)
s.t.
a' 31'v1 11
where 1 is an all-ones vector, and the dependence measure Dep(x, r) can be specially set to the HSIC score with linear kernel, namely,
HSIC Nin (x,r) =
n² tr(KaoH H) tr tr(KxxHrrr 1
n²
no ir( (((KxxH1(0)) - xT
HVT (a HH* HTT 3d V 20- () (15)
n²
Remark: We adopt the HSIC with linear kernel L mainly for favorable optimization structure, as the resulting HSIC score admits a quadratic form w.r.t both a and B. Note that the penalty HSIC term is always non-negative, but the Soft-HGR objective can be negative.

---

Under review as a conference paper at ICLR 2023
The above problem be solved via simple block coordinate descent (BCD) algorithm that updates a and /3 iteratively, see Algorithm [ Essentially, (14) is multi-convex (Xu & Yin 2013). and each
update (line 3 or 4 in Algorithm[] D the sub-problem a linearly constraimed quadratic programming. When the sub-problem is strictly convex, one can obtain the unique minimum in closed-form each update, which admits convergence guarantee a critical point (Grippo & Sciandrone 2000). More details on the subproblem optimization and the landscape study can be found in supplement
Algorithm 1 BCD for problem 14
1: Initialize (0) and (3 I| Use random initializatiin
2: for t=1:T do
a 0. 4: Update 3 c- arg ming ((xl"),)) s.t. 3² v1 0.
3: Update (+ arg min, J(a, ((t,))
5: if stopping creteria is met then
6: return ,
7: end if
8: end for
Remark: We can also impose the invertability of g by limiting the derivative
dy (V) to be positive dy (b22).......
(or negative) in line 4, i.e., T > O, where T = dy (b1),, d
(yn)) E TR
Fine-tune: Algorithm| 1 may produce solutions with distortions, see Figure probably due to the use of the linear kernel. To cope with that, one can enlarge the penalty of dependence X, and use HSIC with universal kernels other dependence measures. Besides, propose a banded loss to reinforce a banded residual plot, see in supplement[ H
4.2 DISTINGUISH CAUSE FROM EFFECT VIA INDEPENDENCE TEST
Following the framework proposed by Zhang Hyvärinen 2009). we distinguish cause from effect according Algorithm 2 We first fit nonlinear models f(-)) 9(-) under hypothesis X + Y. and
f(-)) under hypothesis + Y. After the learning iterations, we conduct independence tests.
If Dep (f) Dep (f(-)) y), the hypothesis X + is supported: otherwise. the hypothesis X + Y is supported. We use the RDC for the independent test, as introduced in Section 3.2
Algorithm 2 The MC-PNL method for causal direction prediction.
Input: The standardized data x, E
I".
Output: The causal score Cx and direction.
1. Fit PNL models via Algorithm aand estimate residuals under hypotheses, X Y and X +- Y.
Under hypothesis X + Y: ((-+) (y-) - (x)
Under hypothesis X Y: f(+) (x-)(x) - (y)-
2. Calculate the causal score Cx- == Dep (f Dop (f)
3. Output the causal score C and
X +Y, Cx > 0,
direction
x Y, if if O,
Towards trustworthy decisions, Liu & Chan 2016 proposed to make no decision when |x is less than a threshold 8 >> 0. Besides, bootstrap [Efron 1992 Zoubir & Boashash 1998) can also be used for uncertainty quantification, see supplement]
5 EXPERIMENTS
In the following, we show the performance of MC-PNL in model learning and its application to bivariate causal discovery.
7

---

Under review as a conference paper at ICLR 2023
5.1 NONLINEAR FUNCTION FITTING
For better demonstration, we generated two synthetic datasets from the PNL model, Y fa (X) +e), and each contains 1000 samples. The data generation mechanisms are as follows,
Synn-: - 10 f2(2) z2, U(0.1, U).
= X
. Syn-2: fx()) sin(7X), = exp(2). ~ U1,, 1), 0.3...
We apply Algorithm to both datasets and show the learned nonlinear transformations well the corresponding residual plots Figure [ The underlying nonlinear functions are correctly learned under the true hypothesis but with certain distortions. We also show that, after fine-tuning with our proposed banded loss or HSIC-RBF loss, such distortion can be fixed up, see supplement H
ADC:0.40
RDC
Essmaadd Modiona
Normalttod
0.5
0.3-
bamaaed
-0.5
-1
Normakize..
Vermsizzex
Norma
Normlloed
1.5
ROCO.82
1.5
RDC.0.99
1.0
Eatmated Rendinn
Nama lood
0.5
0.0
Extmattdd
Nora
Nomrattzdd
Noral
(a) Syn-1
(b) Syn-2
Figure 2: The sub-figures (a) and (b) show the nonlinear function fitting of the two datasets. In each sub-figure, the top row shows the learned (z) (red line) and the residual plot under the correct hypothesis + Y, which has lower RDC value; the bottom row under the opposite X Y.
Convergence Results. We demonstrate the convergence profile of our algorithm with Syn-2, see Figure 3 Results for Syn-1 can be found in the supplement ! The top row shows the snapshots of the learned representations, where we do impose independence regularization (A 0). The algorithm, starting from different random initializations, convergences quickly to the local minimizers sharing the same objective value. The bottom row is with independence regularization À where the solutions have sign symmetry.
5.2 BIVARIATE CAUSAL DISCOVERY
We evaluated the causal discovery accuracy on both synthetic and real datasets.
The generated synthetic datasets all follow the PNL model. And we considered the following two settings: 1) PNL-A: fi are general nonlinear functions generated by polynomials with random coefficients; and fz are monotonic nonlinear functions generated by unconstrained monotonic neural networks (UMNN) (Wehenkel & Louppe 2019); 2) PNL-B: Both fi and fa are monotonic generated by UMNN. The variances of /i, // are rescaled to 1. The input variable X is sampled either from Gaussian mixture (mixG) or uniform (unif) distribution, and the injected noise is generated from normal distributions N nss),, where ns {0.2, 0.4, 0.6, 0.8, 1}. Each configuration contains 100 data pairs, and each data pair has 1000 samples.
Synthetic Datasets:
Gene Data: Discovering gene-gene causal relationships is important application. We used the data in DREAM4 competition the scRNA-seg data (GSE57872) (Han et al. 2017), supplement
D4-S2C) (Marbach et al.) 2009 2010) and
Baselines & Evaluation: Thanks to the implementation by Kalainathan et al. (2020). we can easily compare our proposed method with various existing algorithms. In this paper, compared our
8

---

Under review as a conference paper at ICLR 2023
epcch -0
0.6
spochin
0.4
spoch=10
0.4
epoch=:0
0.2
4020
6.0
2026
00
2020
-04
-1.-5.6-- 64555 1.2 15
15 0.0000033000 05 2.0
1.3.000..5. 05 1.
1.3-1 0,5,5 05 12
epoch o
2.0
1.5
3.5-
epoch-20
1020
1.0
1.0-
4020
0.0
2010
-65
05
-05
-2000
20
15.100.55500 05 15
-.5-1... 50. O5 133
-15-16.4556
Figure 3: The Algorithm converges Syn-2. We plot the snapshots of the feature transformations at training epochs [0, 5. 10, 20]. under 15 random initializations (indicated by colors). Upper: À 0, most iniiaalizations converge to local minimizers (symmetry: (a.3) + (aa, B3). Lower:
5, most initializations converge to two local minimizers (symmetry: (oz, 3) + - (),
proposed algorithm on both synthetic datasets and real datasets with several baseline algorithms, including ANM (Hoyer et al. 2008), CDS (Fonollosa 2019), IGCI (Janzing et al. 2012), RECI (Blöbaum et al. 2018), CDCI (Duong Nguyen 2022], OT-PNL (Tu al. 2022). AbPNL (Uemura & Shimizu 2020p. Our implementation of MC-PNL follows Algorithm (without fine-tuninge and we empirically set À = 5 (the choice of À is briefly discussed supplement D. We also conducted causal discovery the PNI. learned by the ACE algorithm. The ROC-AUC score is used for the evaluation.
Table 2: Comparison of bivaraite causal discovery ROC-AUC on synthetic and real datasets
Dataset
ANM
CDS
IGCI
RECI
CDCI
OT-PNL
AbPNL
ACE!
MC-PNL
PNL-A-mixC
0.256
0.207
0.932
0.537
0.410
0.431
0.645
0.580
0.708
PNL-B-mixG
0.150
0.160
0.908
0.462
0.304
0.309
0.672
0.536
0.771
PNL-A-unif
0.203
0.390
0.681
0.879
0.544
0.711
0.517
0.514
0.617
PNL-B-unif
0.094
0.311
0.866
0.929
0.535
0.536
0.599
0.418
0.608
D4-S1
0.604
0.582
0.380
0.550
0.651
0.474
0.408
0.592
0.646
D4-S2A
0.616
0.580
0.447
0.592
0.673
0.472
0.519
0.558
0.626
D4-S2B
0.521
0.529
0.450
0.491
0.614
0.517
0.501
0.495
0.519
D4-S2C
0.556
0.564
0.441
0.521
0.590
0.490
0.145
0.538
0.576
GSES7872
0.493
0.457
0.599
0.474
0.478
0.538
0.499
Average time? (s)
20.11
7.67
0.50
0.27
0.26
7220
~ 9300
21.68
30.31
Independence test-based methods.
Average running time evaluated on synthetic data containing 100 pairs, and each pair has 1000 samples.
We report the comparison of ROC-AUCs in Table The results averaged five different noise scales for the synthetic datasets, Our proposed MC-PNL consistently outperforms other independence test-based methods the synthetic PNL data. Especially compared with AbPNL, our MC-PNL is not sensitive to the initializations and is much more efficient (w.rt. training time); compared with ACE (without independence regularizer). MC-PNL has better causal discovery accuracy. And for real datasets, our methods is quite competitive.
CONCLUSIONS
In this paper, we focus on the PNL model learning and propose a maximal correlation-based method, which can recover the nonlinear transformations accurately and swiftly in an iterative manner. The key is to incorporate with maximal correlation to avoid learning arbitrary independent noise, and the proposed MC-PNL is more reliable than previous methods that are solely based on the independence loss. Besides the PNL model learning, we conduct experiments on the downstream causal discovery task where MC-PNL is superior the SOTA independence test-based methods.

---

Under review as a conference paper at ICLR 2023
REFERENCES
Adrien Bardes, Jean Ponce, and Yann LeCun. VICReg: Variance-invaaiann--ovvai regularization for self-supervised learning. In International Conference on Learning Representations, 2022. URL
leaming
https://openieviee net forum derrumiid-vae iid xm6 www.YDDDDDDDD 20DDDbb
Yoshua Bengio. Tristan Deleu, Nasim Rahaman, Nan Rosemary Ke, Sebastien Lachapelle, Olexa Bilaniuk, Anirudh Goyal, and Christopher Pal. A meta-transfer objective for learning to disentangle causal mechanisms. In International Conference on Learning Representations, 2020. URL
https::/openreviie net/
wtr.mmirrmmmiiiooo wtt.://www.etttoouudddddd... rYFMI 1
Patrick Blöbaum, Dominik Janzing, Takashi Washio, Shohei Shimizu, and Bernhard Schölkopf. Causc-effect inference by comparing regression errors. In International Conference Artificial Intelligence and Statistics, pp. 900-909. PMLR, 2018.
Leo Breiman and Jerome H Friedman. Estimating optimal transformations for multiple regression and correlation. Journal of the American stafistical Association, 80(391):580-555) 1985.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(Nov):507-554, 2002.
Junsouk Choi, Robert Chapkin, and Yang Ni. Bayesian causal structural learning with zero-inflated poisson bayesian networks. In H. Larochelle, M. Ranzato, Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Neural Information P Processing Systems, volume 33, pp. 5887-58877 Curran As- sociates, Inc., 2020, URL Bttp:tttttt:////ooeeeeiii wttp://www...nnooooooo.....ooo..
neurips /ttpa//ppppcc/ppppp/////////eeee
Diego Colombo, Marloes H Maathuis, Markus Kalisch, and Thomas S Richardson. Learning high- dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, pp. 2012.
Povilas Daniušis, Dominik Janzing, Joris Mooij, Jakob Zscheischler, Bastian Steudel, Kun Zhang, and Bernhard Schölkopf. Inferring deterministic causal relations. In Proceedings of the Twenty-Sixth Conference on Uncertainty Artificial Intelligence, pp. 143-150, 2010.
Bao Duong and Thin Nguyen. Bivariate causal discovery via conditional divergence. In First Conference on Causal Learning and Reasoning, 2022. URL
Forum?1d=8X6cWIIe
2v
Imme Ebert- Uphoff and Yi Deng. Causal discovery for climate research using graphical models. Journal of Climate, 25147):5648-5665 2012.
Bradley Efron. Bootstrap methods: another look at the jackknife. In Breakthroughs in statistics, pp. 569-593. Springer, 1992.
José AR Fonollosa.
Conditional distribution variability measures for causality detection. In Cause Effect Pairs in Machine Learning, pp. 339-3439-493 347. Springer, 2019.
Ronja Foraita, Juliane Friemel, Kathrin Günther, Thomas Behrens, Jörn Bullerdiek, Rolf Nimzyk, Wolfgang Ahrens, and Vanessa Didelez. Causal discovery of gene regulation with incomplete data. Journal of the Royal Statistical Society: Series A (Statistics in Society), 183(4):1747-1775 2020.
Kenji Fukumizu, Arthur Gretton, Xiaohai Sun, and Bernhard Schölkopf. Kernel mea- sures of conditional dependence. In J. Platt. D. Koller, Y. Singer, and S. Roweis (eds.), Advances in Neural Information Processing Systems, volume 20. Curran Asso- ciates, Inc., 2007. URL httpsttt///:::/reeeeeee http://eeeeeeeee...eeeooee.eee/e ccc/pappe//000//11ee "3572443073311222253333628288
Paper.
Hans Gebelein. Das statistische problem der korrelation als variations-und eigenwertproblem und sein zusammenhang mit der ausgleichsrechnung. ZAMM-Journal of Applied Mathematics and Mechanics/Zeiischrjff für Angewandte Mathematik und Mechanik, 21(6):364-379, 1941.
Clive WJ Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: journal of the Econometric Society, pp. 424 438, 1969.
10

---

Under review as a conference paper at ICLR 2023
Daniel Greenfeld and Uri Shalit. Robust learning with the Filbeer-Schhidd independence criterion. In International Conference on Machine Learning, pp. 3759-3768. PMLR, 2020.
Arthur Gretton, Ralf Herbrich, Alexander Smola, Olivier Bousquet, Bernhard Schölkopf, et al. Kernel methods for measuring independence. 2005.
Arthur Gretton, Kenji Fukumizu, Choon Teo, Le Song, Bernhard Schölkopf, and Alex Smola. A kernel statistical test of independence. In Platt, D. Koller, Y. Singer, and S. Roweis (eds.), Advances in Neural Information Processing Systems, volume 20. Curran Asso- ciates, Inc., 2007. URL /poptttts////ooceennnsssss. neurips...eu.p....e.. http://www..aaaaoaooooo////////d
pdf
Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf. and Alexander Smola. A kernel two-sample test. The Journal of Machine Learning Research, 13(1):723-773 2012.
L. Grippo and Marco Sciandrone. On the convergence of the block nonlinear Gauss-Seidel method under convex constraints. Oper. Res. Lett., 26(3):127-136, 2000. URL http:
htm ht
Heonjong Han, Won Cho, Sangyoung Lee, Ayoung Yun, Hyojin Kim, Dasom Bac, Sunmo Yang, Chan Ycong Kim, Muyoung Lee, Eunbeen Kim, Sungho Lec, Byunghee Kang, Dabin Jeong, Yaeji Kim, Hyeon-Nae Jeon, Haein Jung, Sunhwee Nam, Michael Chung, Jong-Hoon Kim, and Insuk Lee. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Research, 46(D1):D380-D386 October 2017. doi: 10.1093/ nar/gkx1013, URL httttp:://di..rggp/pttttp http:///.....///////////////
Hermann O Hirschfeld. connection between correlation and contingency. In Mathematical Pro- ceedings of the Cambridge Philosophical Society, volume 31, pp. 520-524, Cambridge University Press, 1935.
Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008
Dominik Janzing, Patrik O. Hoyer, and Bernhard Schölkopf. Telling cause from effect based on high-dimensional observations. In Johannes Fürnkranz and Thorsten Joachims (eds.), Proceedings of the 27th International Conference Machine Learning (ICML-10) pp. 479 Haifa, Israel, June 2010. Omnipress. URL htllhtp:/// www. wwwiicmm20000
Dominik Janzing, Joris Mooij, Kun Zhang, Jan Lemeire, Jakob Zscheischler, Povilas Daniušis, Bastian Steudel, and Bernhard Schölkopf. information-geommtrii approach to inferring causal directions. Artificial Intelligence, 182:1-31, 2012.
Diviyan Kalainathan, Olivier Goudet, and Ritik Dutta. Causal discovery toolbox: Uncovering causal relationships in python. Journal of Machine Learning Research, 21(37):1-5, 2020. URL http:// /mmpt//mmm..///////// org wttpg//appr///eeee hetppeersvv1/a11 /ttp///111111 187.html
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR (Poster), 2015. httptt/t///aai
org. /aaa////11... 1112. 6980
Joshua Ka-Wing Lee. Maximal Correlation Feature Selection and Suppression With Applications.
PhD thesis, Massachusetts Institute of Technology, 2021.
Furui Liu and Laiwan Chan. Causal inference discrete data via estimating distance correlations. Neural computation, 28(5):801-81,, 2016.
David Lopez-Paz, Philipp Hennig, and Bernhard Schölkopf. The randomized dependence coefficient. Advances in neural information processing systems, 26, 2013.
David Lopez-Paz, Robert Nishihara, Soumith Chintala, Bernhard Scholkopf, and Léon Bottou. Discovering causal signals in images. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 6979-6987, 2017.
11

---

Under review as a conference paper at ICLR 2023
David López-Paz. Philipp Hennig, and Bernhard Schölkopf. The randomized dependence coefficient. In Christopher J. C. Burges, Léon Bottou, Zoubin Ghahramani, and Kilian Q. Weinberger (eds.), NIPS, pp. 1-9, 2013. UtLhttp:/ sestp://dblpp wnw.uniiiittt... weeiierrreeeeeeeee conf/ http///oonn//m//n//s//111111 html #Lopez 2- # PaZHS13
Wan-Duo Kurt Ma, JP Lewis, and Bastiaan Kleijn. The HSIC bottleneck: Deep learning without back-propagation. In Proceedings of the AAAI Conference Artificial Intelligence, volume 34, pp. 5085-5092, 2020.
Daniel Marbach, Thomas Schaffter, Claudio Mattiussi, and Dario Floreano. Generating realistic in silico gene networks for performance assessment of reverse engineering methods. J. Comput. Biol., 16(2):229 29,,, February 2009.
Daniel Marbach, Robert J Prill, Thomas Schaffter, Claudio Mattiussi, Dario Floreano, and Gustavo Stolovitzky. Revealing strengths and weaknesses of methods for gene network inference. Proc. Natl. Acad. Sci. U. S. A., 107(14):6286-229) April 2010.
Joris Mooij. Dominik Janzing, Jonas Peters, and Bernhard Schölkopf. Regression by dependence minimization and its application to causal inference in additive noise models. In Proceedings of the 26th annual international conference machine learning, pp. 745-752, 2009.
Joris M Mooij. Jonas Peters, Dominik Janzing, Jakob Zscheischler, and Berhhrrd Schölkopf. Distin- guishing cause from effect using observational data: methods and benchmarks. The Journal of Machine Learning Research,
2016.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019.
Anoop P. Patel, Itay Tirosh, John J. Trombetta, Alex K. Shalek, Shawn M. Gillespie, Hiroaki Wakimoto, Daniel Cahill, Brian V. Nahed, William T. Curry, Robert I.. Martuza, David N. Louis, Orit Snrenblatt-Rosen Mario I.. Suvà, Aviv Regev, and Bradley F. Bernstein. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science, 14466190) 2014. doi: 10.1126/sciencc 1254257. URL https:/ httt:::::/www..ce. www.:wc.in.n.....mmwomwcm...eeee (aarrg/do///bb//// 1126/ 2766/ccence.. 1254257
Judea Pearl. Causality.
Cambridge University Press, September 2009.
doi: 10.1017/
cbo978051 1803161. htt.httss: wtt.//ooooog/...// 11t1://0.01177700077711111111111
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: foundations
and learning algorithms. The MIT Press, 2017.
Ali Rahimi and Benjamin Recht. Weighted sums of random kitchen sinks: Replacing minimiza- tion with randomization in learning. In D. Koller, D. Schuurmans, Y. Bengio. and L. Bot- tou (eds.), Advances Neural Information Processing Systems, volume 21. Curran Asso- ciates, Inc., 2008. URL htttps::/procceeedii hww.r/occeednnnn.............. nwwrips httpr/pppscccppppee/////////// [ccr2244923375555444444446666
-ddf
A. Rényi. On measures of dependence. Acta Mathematica Academiae Scientiaram Hungaricae, 10(3-4):441-551 September 1959. doi: 10.1007760020225507. URL
1007/bE02024507
Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Camps- Calps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D. Mahecha, Jordi Muñoz-Mari, Egbert H. van Nes, Jonas Peters, Rick Quax, Markus Reichstein, Marten Scheffer, Bernhard Schölkopf, Peter Spirtes, George Sugihara, Jie Sun, Kun Zhang, and Jakob 7scheischler. Inferring causation from time series in earth system sciences. Nature Communications, 10(1), June 2019. doi: 10. 1038//41467-011-10000
019 10105-2
Xinpeng Shen, Sisi Ma, Prashanthi Vemuri, and Gyorgy Simon. Challenges and opportunities with causal discovery algorithms: Application to Alzheimer's pathophysiology. Scientific Reports, 10 (1), February 2020. doi: 10. 10.1038//41598 020- 0200000005522000 59669-x. URL/https: //o.i. org/10.1038.... 841 941598-020-596555
12

---

Under review as a conference paper at ICLR 2023
Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A linear non- non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
George Sugihara, Robert May. Hao Ye, Chih-hao Hsich, Ethan Deyle, Michael Fogarty, and Stephan Munch. Detecting causality in complex ecosystems. science, (3886106):4996-1 2012.
Masashi Sugiyama and Makoto Yamada. On kemel parameter selection in Hilbert-Schmin inde- pendence criterion. IEICE TRANSACTIONS on Information and Systems, 95(10):2564-2555 2012.
Taiji Suzuki, Masashi Sugiyama, Takafumi Kanamori, and Jun Sese. Mutual information estima- tion reveals global associations between stimuli and biological processes. BMC Bioinform., 10 (S-1), 2009. URL http: http:/bblp..e/ uni-trier........ de/ /ttp://d/d//orrrl//sa/a//oo/ /tt/:/oorralls/mmm///ommooo html 4Suzuk
Sergios Theodoridis. Machine learning: a Bayesian and optimization perspective. Academic press, 2015.
Ruibo Tu, Kun Zhang, Hedvig Kjellstrom, and Cheng Zhang. Optimal transport for causal discovery. In International Conference on Learning Representationss 2022. URL /tllttts///ppeeeiiiee
net/ pret/forun: Ldmpumm3cc0
Kento Uemura and Shohei Shimizu. Estimation of post-nonlinear causal models using autoencoding structure, In /CASSP 2020 2020 IEEE International Conference on Acoustics, Speech and Signal doi: N0.1109/CASSP44076..6 000.000777..0009900 2020.9033468.
Processing (ICASSP), pp. 3312-3316, 2020.
Lichen Wang. Jiaxiang Wu, Shao-Lum Huang, Lizhong Zheng, Xiangxiang Xu, Lin Zhang, and Junzhou Huang. An efficient approach to informative feature extraction from multimodal data. In Proceedings of the AAAI Conference Artificial Intelligence, volume 33, pp. 5281-5288, 2019.
Antoine Wehenkel and Gilles Louppe. Unconstrained monotonic neural networks. In Advances in
Neural Information Processing Systems, pp. 1543-1553, 2019.
John Wright and Yi Ma. High--Dmensional Data Analysis with Low-Dimensional Models: Principles, Computation, and Applications. Cambridge University Press, 2022.
Yangyang Xu and Wotao Yin. A block coordinate descent method for regularized multiconvex optimization with applications nonnegative tensor factorization and completion. SIAM Journal on imaging sciences, 1738:1758-1888 2013.
Hao Ye, Ethan R Deyle, Luis J Gilarranz, and George Sugihara. Distinguishing time-delayed causal interactions using convergent cross mapping. Scientific reports, 5(1):1-9, 2015.
Chulhee Yun, Suvrit Sra, and Ali Jadbabaie. Small ReLU networks are powerful memorizers: a tight analysis of memorization capacity. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Bus Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 32. Cur- ran Associates, Inc., 2019. URL httpppps/pppccceeiings..... fi File/dbea3dde2220 7c1 70c4 77c777777777777777111133 wwuruinggs.euuupppp httpppsccccpppp/e1199.
.pff
Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals. Understanding deep learning (still) requires rethinking generalization. Communications of the ACM, 64(3):107-115, 2021.
K Zhang and A Hyvärinen. On the identifiability of the post-nonlinear causal model. In 25th Conference on Uncertainty in Artificial Intelligence (UAI 2009), pp. 647-655. AUAI Press, 2009.
Kun Zhang and Aapo Hyvärinen. Distinguishing causes from effects using nonlinear acyclic causal models. In Isabelle Guyon, Dominik Janzing, and Bernhard Schölkopf (eds.), Proceedings of Workshop Causality: Objectives and Assessment NIPS 2008, volume 6 of Proceedings of Machine Learning Research, pp. 157-164, Whistler, Canada, 12 Dec 2010. PMLR. URL https: Sitts::/prroese ngs ml www.sppessseeee pwwps//////eee..pppppp.. zhang 10..htm
13

---

Under review as a conference paper at ICLR 2023
Kun Zhang, Zhikun Wang, Jiji Zhang, and Bernhard Schölkopf. On estimation of functional causal models: general results and application to the post-nonlinear causal model. ACM Transactions on Intelligent Systems and Technology (TIST), 7(2):1-22, 2015.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. DAGs with NO TEARS:
Continuous Optimization for Structure Learning. In Advances in Neural Information Processing Systems, 2018.
Abdelhak M Zoubir and Boualem Boashash. The bootstrap and its application in signal processing. IEEE signal processing magazine, 15(1):56-76, 1998.
14

---

Under review as a conference paper at ICLR 2023
Supplementary Document
Maximal Correlation-Based Post-Nonlinear Learning for Bivariate Causal Discovery
A DISCUSSION ON MI MINIMIZATION
was shown that minimizing the MI in 5 is equivalent to maximizing E logg (r(-)) E log d ()) (w) Zhang & Hyvärinen 2009,, where is the assumed noise density, We find this objective interpretable, since the first term, I log p (r(-)), can be understood the data fitting be understood from information-geomeeni perspective DDaniušis et al. (010)) However, such equivalent form requires a known noise distribution calculate the log-likelihod. Some works (Ma al. 2020 Uemura Shimizu 2020) have been proposed to avoid this difficulty by using HSIC instead of MI.
term, and the second term, E log dy ((-) (v)
d
B EXPERIMENTS ON MINIMIZING HSIC
In this section, we show the PNL model learning result by minimizing L2 generated two synthetic datasets from PNL model, Y f2 (fi (X) ) + c), and each contains 1000 data samples. The data generation mechanisms are as follows (see Figure D
Syn-1: fi() X-1 10X, f2(Z) Z X (0)
. Syn-2: = f2))) = -,00
f
g"
scatter plot
Figure 4: The ground truth transformations of f" and g" of Syn-l (top) and Syn-2 (bottom).
We build different MIPs with the following configurations.
Narrow deep MLP: the input and output both one-dimensional; there are hidden layers, each with 5 neurons. The activation function is Meaky-ReLUU
Wide over-parameterized MIP: the input and output are both one-dimensional: there is only single hidden layer with 9000 neurons. The activation function is Leaky-ReLU.
We use the default initialization method in Py (Paskke et al. 2019) and make sure the exact same initial weights for narrow/wde MLPs are used (i.c., the initializations for different datasets the same).
Optimization Setup:
We set the batch size be 32. We use Adam (Kingma & Ba 2015 for the optimization (the learning rates are 10-3 and 10-6 6 for narrow deep and wide over parameterizeed MLPs, respectively, while all other parameters are set by default).
We report the learning results in Figure S The learned transformations (see row and row 4 in Figure 5) deviates far away from the underlying functions, and are quite similar across datasets. The
15

---

Under review as a conference paper at ICLR 2023
Wey Ove Po
Whd Outo-rraneeeted RtP
24
14
12
10:
286
364
214
62
opcc1
spon
34
660
535
316
28
-25-30-15
33.36-13
-16
-15
11
01 18
-3 11----111 03
-15-17-55
835
7.
665
Kn
Syn-1, narrow MLPs
Syn-1, wide MLPs
Syn-2, narrow MLPs
Syn-2, wide MLPs
Figure Visualization of the learned nonlinearities (trained solely with HSIC, under different datasets and MLP configurations). From top to bottom, the convergence results, residual plot, learned f, learned g. are plotted. Each column shows one specific configuration. None of them learns meaningful nonlinearities, and the learned transformations are quite similar across datasets.
possible reason is that, the solutions were started from the same initialization and trapped the local minima near the initializations.
To verify whether such HSIC-based PNL learning algorithm stable for causal discovery, we further evaluate the AbPNL on the following dataset. We build 100 data pairs with different random seeds, following the same mechanism, Syn- and each contains 1000 data samples. And we applied the AbPNL (Uemura & Shimizu (020) with different initializations on each of those data pairs. The results in Table show that the causal discovery stableness for AbPNL not satisfactory.
Table 3: Comparison of bivaraite causal discovery AUC on 100 realizations of Syn-1
Dataset
ANM
CDS
IGCI
RECI
CDCI
AbPNL
ACE
MC-PNL
Syn-1
0.495
1
0.528
0.281
1
C SYNTHETIC DATASETS FOR INDEPENDENCE TEST
In this section, we describe the synthetic data generation from PNL model for the independent test. The data were generated from the following model, Y fz (fi (X) + c) X ~ GMM, ~ N(0, o²). where f1, f2 are randomly initialized monotonic neural networks (Wehenkel & Louppe 2019) with layers and 100 integration steps, and each layer contains 100 units. The cause term X is sampled from a Gaussian mixture model as described in Lopez-Pazz et al. (0017). The datasets were configured with various noise levels and sample sizes. There three different injected noise levels, E {0.1, 1, 10). and three different sample sizes, N e {1000, 2000, 5000}. And under each configuration. we generated 100 data pairs for evaluating the independence test accuracy.
16

---

Under review as a conference paper at ICLR 2023
D A UNIVERSAL VIEW OF DEPENDENCE MEASURES
Actually the discussed dependence measures in Section 3.2 are all closely related to the mean squared contingency introduced by (Rényi 1959) and rediscovered due to its squared version called squared-lor mutual information (SMI) (Suzuki et al. 2009).
SMI /f p(o))())
p(z,z) [( p(r)p(u)
drdy // -1)²
))
p(a, v) p(z)p(()
p(a, y)dady -1 1.
(16)
When the density ratio DR(x, p(a)p(y) is constant (namely X and Y are independent), the SMI should be zero. To estimate the SMI, one can first approximate DR(x, y) by a surrogate function [z, y) parameterized by 0. The optimal parameter can be obtained via minimizing the following squared-error loss JDR
pla.v)
(0) l/ (DRo(z,)) (2, y) DR(a, v))² p(z))(y)dzdd
(17)
// (Ro( (())(())xxd - // DRo(E, v)((z, y)drdy + Const.
Then the empirical SMI can be calculated as, SMI 1/4 Σj= (Rg(2,
We show that, with different parameterizations of the density ratio, the resulting SMI will be equivalent to different dependence measures, see Table
Table 4: Connections between DR parameterization and dependence measure
Density ratio surrogate function (R)(1,))
Corresponding dependence measure
DR (z, y) = 1+
0,K (a, x.)) L (y, 34)
variant of LSMI (Sugiyama & Yamada 2012
DR v)
Ei=1 1= K
(z, 2x)
L (y. V1)
HSIC (Gretoo et al. 2005)
DRa (z, v) = 1 + D
1 (1( (Ro v) = 1 (x)g(())
m-mode HGR correlation (Wang et al.
2019
HGR correlation (Rényl 1959)
When fig g are the linear combinations of random features. [(x) = a² ((x).9()) = 3 (vy) the corre- sponding dependence measure will be RDC [Loez-PPaz et al. 2013).
Sugiyama & Yamada 2012) proposed to approximate the density ratio by DR y)
En 0₂K (x, x4) L (y, 3.). where θ has closed-form solution via minimizing 17 After then, they approximated the SMI using the empirical average of Equation (16). Ej= DRa (T,
- Σラ 0,K (I. 2;) L (v. 1. It is shown that, the first term is actually the empirical HSIC, when (6: = 1 We argue that there is a flaw above, when X and Y are independent, both the SMI and HSIC score should be zero. A simple modification is to model the density ratio by DR 3) 1 + >x O,K( (x. z;) L (v, y:). The constant 1 here is to exclude all the independence terms, and the rest ones should model the dependency only. This modification will not hurt the quadratic form of J DD (0), and maintains good interpretation. The SMI reduced to HSIC score, when (0.)2 A.
We extend this idea to approximate the density ratio by DRo(z,y) = = f(x)g(y), where f. g zero mean and unit variance functions parameterized by 9, the resulting SMI will be equal to the HGR maximal correlation. Similarly, the constant 1 will capture the independence part, and f())g(v) will capture the dependencies.
Proposition 1. The density ratio estimation problem (17 is equivalent to the maximal HGR correla- tion problem (7), when the density ratio is modeled in the form of DR y) = 1 1+ f(z)g(y), and g are restricted zero mean and unit variance functions.
17

---

Under review as a conference paper at ICLR 2023
Proof. We substitute (2, y) into Equation 17
JD (f.g) = // (1 +1+ (())a//)(( //
// (1 (())(((()))((l)e) y)dzdy + Const.
(B((())) var(f(X))var - 2-2B(((()g())) + Const.
Then it is not hard to see, finf-gg ln JDR (J, g). subject to B(J) (E(g) 0, var(s) var(g) 1, is equivalent to the maximal HGR correlation problem
Proposition The density ratio estimation problem 17 is equivalent to the Soft-HGR problem O when the density ratio is modeled in the form of DRo (T. = 1 f())g((), and J. 9 are restricted to zero mean functions.
We further note that the above density ratio estimation can be regard truncated singular value decomposition DR (2 y) = 1 + += f((z) (y), where m 1. When letting m > 1 and imposing zero mean and unit variance constraints on all f. and g.. the corresponding JDR minimization problem
is equivalent to solving the m-mode HGR maximal correlation (Wang al. 2019 Lee 2021).
Definition 3 (m-mode HGR maximal correlation). Given I m S min (X) DI), the m-mode maximal correlation problem for random variables E X, is,
(f" gf") A
arg max F- th gyy-RR
IE (X)(())(
(18)
B[ =
where f = J2, ...,.,.......................... mm] .g = [/1 92.... .mll are referred as the maximal correlation func- tions.
E CONNECTIONS AMONG MI, ML, AND MC
In this section, we build connections among minimizing MI, maximum likelihood, and maximal correlation. The equivalence between minimizing MI and maximizing likelihood was built in Zhang & Hyvärinen (2009). The following proposition shows the connection to maximal correlation.
1 is generated from PNL model = gg- (f() )), where f.g are both invertible functions, and the noise follows Gaussian density p(e, with zero mean and variance 0, then maximizing the log- lig-liieldoood log p( (Tx, y:)) y:)) is equivalent to solving the regression problem 8
Proposition 3. Suppose the dataset yr),,
Proof. Under proper assumptions in proposition 3 the log-likelihood can be written as follows,
(n.gg)
log logp(a pa.plly...... f,
log p(x:))
log bog(yy
Σ log p(x))
logrlg(v.)/(x:)))))) (f.g are invertible)
are
Σ
log p(rr)
log logr(g(v:) (xx))=0) (from PNL model)
(((()) f(x:))2
log p(x))
+ log
(Gaussianity) (19)
20
/270
is not hard to see, with fixed 6. maximizing the log-likelihook L, (f.gg g) is equivalent to minimizing [xf g(x) with invertible f and g. Without loss of generality, one make f, g zero mean. To avoid trivial solutions, one can further restrict to have unit-variance. Then the equivalence to the regression problem 8 is build.
n
18

---

Under review as a conference paper at ICLR 2023
Corollary 1. When - the ground truth transformations f .9", minimize the MI(x, i' to zero. achieve optimum of (8) and maximize the log-likelihood (nnff,
Proof. The proof is directly follows Theorem 3 in ZZhng et al. 2015).
The reformulating to or allows efficient BCD-like optimization algorithms to be exploited.
F RANDOM FEATURE GENERATION
We generate the random features as described López-Paz et al. (2013). The generation process has
the following two steps: copula transformation (optional) and random nonlinear projection.
Step Copula transformation. We first estimate the empirical cumulative distribution of both X and Y by,
PX (x) == 1
Σ (x, S x) (x) = Σ
R
1 5))
Then we can apply the empirical copula transformation to data samples ((z,, y:) ux P. (x:)
and uz P. (3)), where the marginals U- and U. follow uniform distribution (O0, 1).
Step 2. Random nonlinear projection. We design a k-dimensional random feature vector ((x) [sin(w;æ: + b.),... sin(wg + by)]", where wi, by NO0,, s²). The random feature matrix + gkxn is stacked as
sin (wir) + b1)
sin (Wien + bi)
(x;
sin (wy21 by)
sin (wgg + bk)
One replace the H; here by ux from the first step form the random feature matrix. Similar procedures can be applied to as well to generate V. The number of random Fourier features is user-defined, which typically chosen from few tens few thousands (Rahimi & Recht 2008 Theodoridis 2015). In our experiments, we set = 30 and s 2.
G
ON THE OPTIMIZATION OF PROBLEM (14)
G.1 SUBPROBLEM: EQUALITY CONSTRAINED QUADRATIC PROGRAMMING
To simplify the notation, rewrite the sub- sub-problem into the following form,
min f(x) = 1 Ax bTx,
ZER
(20)
s.t.
oJ Tz=
With the KKT conditions, one can find the unique optimal solution x" by solving the following linear system,
3 (-)).. A
I
(21)
0
=KKT
when the KKT matrix is nin-singular. our setting, we can choose and properly to make ΦΦT and VADT positive definite, or add small positive definite perturbation matrix e/, such that the unique optimum would be obtained. Besides, sub-problem of smaller size and easy to solve.
We that the optimal solution of is also solution of
19

---

Under review as a conference paper at ICLR 2023
G.2 LANDSCAPE STUDY WITH HESSIAN
To simplify the notation, we rewrite
J(a, A), a" AaB" B3 a²c3+ a Da B² EB,
(22)
where,
A DΦ"
2n
B V
C 1 DUUT
(n-12 Hyx
(23)
H K xx H4²
(n-17
11 VH VHKxxll u E 1n-13
And the corresponding Hessian is
2AB BB + 2D Ba Ba A - C
Aag B C 2Ba Aa + 2E
V))(a,3) =
(24)
Now are able to verify the property of the critical points via checking their Hessians numerically.
One obvious critical point is the all zero vector O. From our experiments, the Hessian at 0 is mostly indefinite, as long the convex regularization term is not too huge, which means 0 is saddle point. In practice, the algorithm rarely converges to O.
H FINE-TUNE WITH BANDED Loss / UNIVERSAL HSIC
In the PNL model, the injected noise are assumed to be independently and identically distributed. Thus, the residual plot should forms a "horizontal band". We design banded residual loss to fine-tune the models as follows. The data samples are separated into bins (x(*) .y( 16
1==1 according to the ordering of X, and we expect the residuals in those bins Res; f(x()) g (1) to have the
same distribution, see Figure 6 To this end, we adopt the empirical maximum mean discrepancy measure of distribution discrepancy. The banded residual loss
(MMD) (Gretton et al. 2012)
defined band (MM))
E (IMD(Res, Res..).. where Res f(x) g Then
append this y--penaiized banded loss to Problem (14) as.
min J(a, 3) + μ (MMD)
s.t. a a' = 3T v1 3
0.
(25)
1.0
Res:
Ros;
Ross
0.5
Residual
-0.5
-1.0
Res,11
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
X
Figure 6: The construction of banded residual loss.
The above banded residual loss involves MMD, which is highly non-convex and brings difficulties to the optimization. We used the projected gradient descent with momentum to optimize the loss function. The residual plot shows band shape, top row in Figure
20

---

Under review as a conference paper at ICLR 2023
1.3
ADC.0.11
RDC
Exiraaed NccAbbol
1.0
Nor Hornaktedd
Normullzed
-0.5
0
-1
-1.5
-10
Normalize...
Vermsizzex
Norma
Norallozd
0.5
1.5g
RDC.0.0
0.0
3.0
Detimsted (tandtul
1.5
1.0
-0.
-24
Norma
Normaltess
Namaa and
(a) Syn-1
(b) Syn-2
Figure Fine-tuning with the banded residual loss.
We also show the results of fine-tuning by enlarging the penalty (to À = 10000) HSIC term with universal Gaussian RBF kernel Figure
Definition 4 (Universal Kernel (Gretton al., 2005).. A continuous kernel k(-..) on compact metric space (X, d) is called universal and only ij the RKHS FF induced by the kernel is dense in (XX), the space of continuous functions on X, with respect the infinity norm [| g glloo.
RDC0011
1.5
RDC-0.11
1.0
Mamna bod
Exsimatec Rus
0.3
2
0.0
0.5
-1.0
Normolled.. X
Normatzed
Norma lood
Nomallze
RDC.0.0
Polmatedd SANGON
0.0
0.0
-0.
1.00
-15
-2-
Nomm
Nomal
(a) Syn-1
(b) Syn-2
Figure Fine-tuning with the HSIC-RBF loss.
BOOTSTRAP FOR TRUSTWORTHY CAUSAL DISCOVERY
Bootstrap is commonly used technique to estimate the confidence interval. In this section, we show a few examples of bootstrap with Tuebingen data (Mooij et al. 2016). We obtained 30 estimates of RDC from the data re-sampled with replacement, Figure The blue bars indicate the RDC distribution under the true causal direction; The orange bars indicate the RDC distribution under the false causal direction.
21

---

Under review as a conference paper at ICLR 2023
Pair
Pair
341
RDC
RDC
RDC
Pair
Pair
Pair
Pair
dentin
05
36
Pair
Pair 7
Pair 8
Figure 9: Bootstrap results of MC-PNL eight Tuebingen datasets. We plot the histogram of the RDC estimates and the estimated causal scores of 30 replications.
J ADDITIONAL CONVERGENCE RESULTS
In this section, we show the convergence results on Syn-1 as well.
-p
spoches
epach-10
epach-20
8.5
05
05
8920
6030
6.
03
0.2
02
62:
0.
0.0
00
-61
-61
-15-5 D-0.5 05 23 15
-1-110-3300 05 1.
135-0--300 05 15 15
-15-19--336 05 10 15 X
epoch o0
spoches
spoch-10
epoch-20
1.5
1.5
1. L.
800
0.0
2020
-35
-15
J.S. 65 13
0-0503 35 10
-1.-1000 05
-1.51.0-552 05
Figure 10: The Algorithm converges on Syn-1. We plot the snapshots the feature transformations training epochs 10, 20], under 15 random initializations (indicated by colors). Upper: 0, most initializations converge to local minimizers (symmetry: (ax + (ac, a -13). Lower:
most initializations converge to two local minimizers (symmetry: (a, B) + (a,,))
K DETAILED DATA DESCRIPTIONS
In this section, we describe the datasets in detail.
Gene Datasets:
22

---

Under review as a conference paper at ICLR 2023
For D4-SI, D4-S2A, D4-S2B, D4-S2C, used the preprocessed data in Duong & Nguyen 2022 D4-S1 contains variable pairs with 105 samples in each pair; D4-S2A, D4-S2B, D4-S2C contains 528, 747, and 579 variable pairs respectively, and each pair contains 210 samples.
The GSE57872 dataset built on Patel et al. (2014), in which the data has continuous values. Following Choi al. (2020). first screen 657 gene pairs that have corresponding labels in the TRRUST database (Han et al. (0177) The gene contains many repeated values. we examined each
gene pair and deleted those repeated expression values.
ON THE CHOICE OF
We tried seven different values for X, and report the AUC scores the PNL-A-unif dataset with different noise levels. We found that the MC-PNL is suitable to use the small noise regime. We also found that for the data with small noise, smaller preferred; and for the data with large injected noise, larger À is preferred.
0.2
0.7
0.4
AUC
0.6
0.6
0.5
0.8
1
0.4
0.3
Figure 11: The detailed AUC scores À under five noise levels PNL-A-unif data.
hhtt:://gtttuhccc//aasssCCCCC
23


Reviews:
[{"review_id": "_ktlvkaZP0", "reviewer": "Reviewer_E1Jp", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper proposes a new bivariate causal inference algorithm called MC-PNL (maximal correlation with independence regularisation). The authors propose to use the randomised dependence coefficient (RDC) instead of HSIC test often used by the causal inference community. \n", "summary_of_the_review": "Some points need to be clarified. \nI have some concerns about the numerical experiments, since the reported results outperform some existing methods but are much worse compared to some other state-of-the-art approaches. \nA number of assumptions made in the theoretical part are not mentioned explicitly. ", "strength_and_weaknesses": "Strengths. The paper is in general well-written. It does a good overview of the state-of-the art methods. The proposed method is novel and worth being studied. \n\nWeaknesses. It is stated that the proposed method outperforms the state-of-the-art independence test-based methods (bivariate causal inference), however, the method is much worse than some other state-of-the-art causal discovery methods. The results of the numerical experiments are not so convincing. I see that the idea is to show that the method outperforms the test-based methods, however, the IGCI method which is also based on the same asymmetry (independence of mechanism) hypothesis performs much better. Do you have some intuition? Are the independence tests used not efficient enough? And the future work should focus on a more efficient tests?  \n\nFigure 2 is unclear. What is the true causal direction? And how the plots should be interpreted? Taking, e.g., subplot a), above and below, we clearly see that the red line is well-fitted. What is the direction? And what can be deduced from the estimated residual values? \n\nAlgorithm 2 abstains from the decision only if C = 0. It is not too strict for real (noisy) applications? Were it be not better to introduce some \\epsilon, so that if abs(C) < \\epsilon the decision can not be made? As it is done, e.g. in F. Liu and L. Chan,“Causal inference on discrete data via estimating distance correlations,” Neural Computation, vol. 28, 2016.  \n\nFinally, in the problem formulation (equation 13), there are two parts: the HGR correlation and RDC (instead of HSIC which is, however, seems to be used in the paper). Could you study the efficiency of each of these parts?  \n\nSection 2.2: generating process and the causalities: I guess a number of assumptions are made behind the description of the idea, and these assumptions need to be mentioned. \n\nIn the numerical experiments (Section 5.2.) did you fix C = 0 (as in the Algorithm 2)? \n\nWhile reading, I have an impression that the HSIC takes too much room in the paper, finally the HSIC is criticised and is proposed to be replaced.\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "_iJJXNrldKb", "reviewer": "Reviewer_9qg4", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "The authors suggest using the randomized dependence coefficient (RDC) instead of the Hilbert-Schmidt independence criterion (HSIC) for the independent test for performing bivariate causal discovery using MLPs in the post-nonlinear setting. ", "summary_of_the_review": "Proposing using a different independence measure is not a sufficient contribution in my opinion. This is analogous to having a different algorithm around the PC algorithm by simply using different conditional independence testers. Please see below for some more minor comments.\n\nOn Justification of the Method:\n\"However, the state-of-the-art (SOTA) PNL-based algorithms involve highly non-convex objectives for neural network training, which are\ntime-consuming and unable to produce meaningful solutions with finite samples.\"\nCausal discovery literature does not generally rely on neural network training. So I find this justification a bit inadequate. How about fitting a function class then checking residual independence as is typically done with ANM models?\n\n\"moreover, the discovered DAG may not necessarily be causal.\"\nThis is not a fair criticism either. The experts are aware that GES would return a graph in the MEC. \n\n\"In this paper, we will focus on a more fundamental problem\"\nthe relative fundamental-ness of bivariate vs. full graph discovery is subjective. I would suggest authors refrain from such subjective comparative statements.\n\nOn Experiments:\n\nA typical benchmark used for bivariate causal discovery is the Tuebingen dataset. Did you test your algorithm on this real data? It would be nice to address why this was not added. \n\n", "strength_and_weaknesses": "The practical take-aways such as the strength of RCD with finite samples could be useful. \n\nAuthors have not tested their algorithm on Tuebingen dataset, which is the baseline for bivariate causal discovery.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "_Sf2OLxU8y", "reviewer": "Reviewer_7wPL", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper focuses on the challenges of estimating post-nonlinear models (PNLs) for causal discovery in the bivariate case. Indeed, solving the practical issues of estimating PNLs is an important topic.\n\nThe problems with existing methods are that they can produce local minima which are trivial and meaningless solutions and that the misspecification of noise distributions can lead to wrong causal discovery results.\nThe paper proposes an objective that avoids relying on noise distributions (similar property as HSIC minimization methods) and enforces the independence of residual and potential direction cause. \n\nIt combines the objective function in the alternating conditional expectation (ACE) algorithm (Breiman & Friedman, 1985) with a dependency measure as the penalty term for enforcing the residual to be independent of the potential direct cause. And it determines the causal direction by using the independent tests of the residual and potential direct cause. \n\nThe main experiments in Table 2 cannot show that the proposed method is significantly better than the other methods empirically.\n", "summary_of_the_review": "The paper works on the practical issues of estimating PNLs and shows the problems of existing methods. It proposes to use the objective,  maximum correlation with a dependency measure as the penalty term. My main concern is about the theoretical guarantee and the identifiability of the results given by the proposed objective. My minor concern is about the experiments and the causal score used in Algorithm 2.", "strength_and_weaknesses": "## Writing \n(+) The paper is well-written in general and makes the problem that it is solving and the contributions clear. \n\n(-) However, I would also suggest spending more effort on the proposed method (Page 6 - Page 7) than illustrating the problem (Page 2 - Page 5). Because I would expect the readers can know some about the problems but not really about the proposed method. It can be not the case but then perhaps readers care more about the concerns of the proposed method.\n\n## Experiments\n(+) The work provides transparent and nice-illustrated figures. \n\n(-) But in Table 2, it cannot show the proposed method stands out compared with the others empirically. A fairer claim is that it is better than the related works based similar framework. As for real-world data, the ANM-based method performs surprisingly well and is comparable with the proposed method. \n\n(-) Moreover, in this case, it would be good to show experimental results on the typical benchmark dataset, the 100 (or more) Tuebingen cause-effect pairs. Not necessary to show superiority, but to provide a comparison with the others.\n\n## Concerns about Algorithm 2\n(-) The causal score is a single number of the difference between the dependency measures in different directions. It would be more convincing to report the numbers for the experiments, especially to which degree (or within which threshold), we can believe the number is larger than zero, smaller than zero, or equal. Moreover, I am worried that in real-world data or finite samples, how much we can trust a single number that gives us the causal relationship. And it doesn't seem to show the experiments about the independent case.\n\n\n## Proposed method\n(+) The work clearly formulates the proposed method with a discussion of optimization process. \n\n(-) Nevertheless, __my main concern__ is that the proposed objective (13) lacks a justification and theoretical guarantee from the identifiability concern. ( Note that I am not asking for an identifiability proof, but the justification and theoretical guarantee for that given the PNLs under identifiability assumptions, can the method be used for determining the causal direction?)\n\nAs actionable feedback, can the authors show that  \n(i) for objective (8), is it in any way related to maximum likelihood or minimizing mutual information?\n\n(ii) as for the model in the causal direction and the one in the reverse causal direction, can objective (8) imply that the one in the causal direction has a smaller value of the optimal objective function value than the other one? \n\nFurthermore,\n\n(iii) as for objective (13), how can the property maintain by adding the penalty term? A concern is that the penalty term is added as a \"soft\" constraint which is not necessary to be exactly the case. Then, is it possible that the optimal solution by solving (13) can be the local minima which are taken as a trade-off between the objective (8) and the penalty term? Then, will this lead to a misspecified model ? will this lead to a problem for causal discovery?\n\n\nTo further elaborate on my point:\nThe authors introduced the PNL-MLP of Zhang & Hyvärinen (2009), which uses mutual information for estimating the model and later uses independent tests for causal discovery. And the authors point out the problem of using mutual information as the objective, which can be hard with large-scale datasets. But an important fact of using the maximum likelihood or minimizing mutual information is that they are well justified by the identifiability of PNLs as illustrated in [1], especially, the independent noise assumption. \nSimilarly, as for the regression by dependence minimization (Mooij et al., 2009), it directly minimizes the HSIC score to enforce the independent noise assumption and pick up the model with a smaller score as the causal direction. This is fine because it directly uses the assumption as objective from the perspective. But for the objective in (Uemura & Shimizu 2020) and this paper, they neither directly use the independent noise assumption nor have a theoretical guarantee of the identifiability as Thm.2 and Thm. 3 in [1]. Therefore, to fix my concerns, maybe the authors could consider my actionable feedback.\n\n[1] Zhang, K., Wang, Z., Zhang, J., & Schölkopf, B. (2015). On estimation of functional causal models: general results and application to the post-nonlinear causal model. ACM Transactions on Intelligent Systems and Technology (TIST), 7(2), 1-22.\n\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "NxRo9bVRCYz", "reviewer": "Reviewer_2uGF", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "This paper studies bivariate causal direction learning from observational data. The authors focus on the general bivariate model, i.e., the post-nonlinear model, and propose a new method to learn the model. The contributions are as follows, the authors first analyze the drawbacks of the existing estimation methods, e.g., the PNL-MLP algorithm, and the AbPNL algorithm. Then they propose a new optimized function to estimate the PNL model by using the Soft-HCR. Finally, the authors apply their methods in both synthetic and Gene datasets to verify the efficiency. ", "summary_of_the_review": "The authors make some progress on the hard problem of estimating the causal direction of interest only using the observed variables.  The proposed method is non-trivial and useful with finite samples. However, some results are based on simulation analysis rather than theoretical analysis.", "strength_and_weaknesses": "Strength\n\nThe authors focus on the challenge of causal direction inferring. This is an important but challenging problem. \n\nThe paper builds on prior work in the field of casual discovery in PNL model. The authors put some effort into mitigating the problem of time-consuming and unreliable results with finite samples. \n\nThe basic idea is interesting and useful.\n\n This paper is well-written and well-organized.\n\n\nWeakness\n\nSome results are based on simulation analysis rather than theoretical analysis.\n\n\nSome concerns or questions\n\n1. The last but one paragraph on Page 2: a biased HSIC? It may be unbiased?\n\n2. ﻿In Table 2, the ROC-AUCs of IGCI have better accuracy than the proposed method in the synthetic datasets. However, IGCI has worse accuracy in real-data sets. Can you explain this result?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject