On the Identifiability of the Post-Nonlinear Causal Model

Kun Zhang
Dept. of Computer Science and HIIT
University of Helsinki
Finland

Abstract
By taking into account the nonlinear effectof the cause, the inner noise effect, and themeasurement distortion effect in the observed
variables, the post-nonlinear (PNL) causal
model has demonstrated its excellent performance in distinguishing the cause from effect. However, its identifiability has not been
properly addressed, and how to apply it in
the case of more than two variables is also a
problem. In this paper, we conduct a systematic investigation on its identifiability in the
two-variable case. We show that this model is
identifiable in most cases; by enumerating all
possible situations in which the model is not
identifiable, we provide sufficient conditions
for its identifiability. Simulations are given
to support the theoretical results. Moreover,in the case of more than two variables, weshow that the whole causal structure can befound by applying the PNL causal model toeach structure in the Markov equivalent classand testing if the disturbance is independentof the direct causes for each variable. In thisway the exhaustive search over all possible
causal structures is avoided.

1

INTRODUCTION

Traditionally, causal discovery algorithms, which may
be constraint-based or score-based, produce a Markov
equivalent class of the causal models, in which somecausal directions may be undetermined (Spirtes et al.,
2001; Pearl, 2000). On the other hand, a functionalcausal model, which expresses each variable as a function of its direct causes and the independent disturbance (Pearl, 2000), if well specified, can explain thedata generating process and help find all causal relations among the variables. For example, under the

Aapo Hyvärinen
Dept. of Computer Science, HIIT
and Dept. of Mathematics and Statistics
University of Helsinki
Finlandcondition that causal relations are linear and acyclic,
the non-Gaussianity of the disturbances could help
find the whole causal model uniquely (Shimizu et al.,2006), by resorting to the independent component
analysis (ICA) technique (Hyvärinen et al., 2001).
However, if the assumed functional causal model is notcapable of approximating the true data generating process, the results may be misleading. Therefore, if thespecific knowledge about the data generating mechanism is not available, to make it useful in practice,
the assumed causal model should be general enough,such that it can reveal the data generating processes
approximately; at the same time, the model should
be identifiable, i.e., it is asymmetrical in causes andeffects and is capable of distinguishing between them.
Although the linearity assumption greatly simplifies
causal analysis, in some situations nonlinear effects
in the system are not negligible. In particular, therecently proposed post-nonlinear (PNL) causal modeltakes into account the nonlinear effect of the causes,
the noise effect, and sensor or measurement nonlinear distortion in the observed variables (Zhang &Hyvärinen, 2008). Mathematically, with the causal
structure represented by a directed acyclic graph(DAG), it expresses each variable xi as
xi = fi,2 (fi,1 (pai ) + ei ), i = 1, ..., n,(1)where pai contains the direct causes of xi , fi,1 denotes
the nonlinear effect of the causes, ei is the independent
disturbance, and fi,2 denotes invertible post-nonlinear
distortion in variable xi . It includes the so-called additive noise model (Hoyer et al., 2009) as a special casein which the nonlinear distortion fi,2 does not exist. It
is also related to the idea of causal reasoning based on
evaluating the complexity of conditional densities (Sunet al., 2008), since it intrinsically admits a simple expression for the conditional density of the effect given
the causes. This model has been used to distinguishbetween causes and effects for the “Cause-effect pairs”
task (Mooij et al., 2008) in the second causality challenge, and gave clearly the best results (the identifiedcausal directions are correct for all eight data sets).
Its good performance is partially due to the allowanceof the measurement distortion fi,2 , which is frequently
encountered in practice.

We assume the causal model holds in both directions
x1 → x2 and x2 → x1 , and show that this implies
some very strong conditions on the distributions andfunctions involved in the model.

Despite its success in solving some real-world problems, there are two unsolved problems related to the
PNL causal model, and they are addressed in this paper. One is the identifiability, a crucial issue, of thismodel. Although it was supported by empirical results and was touched in Zhang and Hyvärinen (2008),
it is far from complete. Here we give a systematicalinvestigation of its identifiability in the two-variablecase, under the assumption that the density of the disturbance has an unbounded support. We show thatthe model is generally identifiable, and give all the
non-identifiable cases, some of which are illustrated
by simulation studies. Our results also have some byproducts. Previously, Zhang and Hyvärinen (2008) investigated this issue, by relating this model to the PNL
mixing ICA problem (Achard & Jutten, 2005), under
the constraint that the nonlinear effect of the cause,
fi,1 , is invertible. Our findings reveal that their results
are not precise, and further provide counterexamples
to the separability theorem of the PNL mixing ICA
problem reported in Achard and Jutten (2005), which
their results are based on. This finding may also beof interest to the ICA community, since PNL mixing
ICA is an important and widely used nonlinear ICA
model. In addition, our results on the identifiability of
the PNL causal model also apply to the additive noisemodel (Hoyer et al., 2009), since the latter model is a
special case of the former one.

Assume that the data (x1 , x2 ) are generated by the
post-nonlinear (PNL) causal model with the the causal
relation x1 → x2 . This data generating process can be
described as

The other problem is how to find the causal relations among more than two variables implied by thismodel. One may search all possible causal structures
and test if they are consistent with the data in a bruteforce way, but it involves high computational load
and becomes impractical as the variable number increases. We provide some fundamental results for thisissue. It is shown that causal discovery based on the
PNL causal model for more than two variables can beachieved in two steps: after obtaining the equivalentclass, one can identify the undetermined causal relations by applying this model to the causal structuresin the equivalent class and testing if each disturbanceis independent of the parents associated with the same
variable. Consequently, the search space is greatly reduced, and statistical tests of mutual independence between more than two variables are avoided.

2

IDENTIFIABILITY

In this section we focus on the two-variable case. We
investigate the identifiability of the PNL causal model,
in particular, its direction, by a proof by contradiction.

x2 = f2 (f1 (x1 ) + e2 ),(2)where x1 and e2 are independent, function f1 is nonconstant, and f2 is invertible. If the other causal direction, x2 → x1 is true, the data generating process
given by the PNL causal model is
x1 = g2 (g1 (x2 ) + e1 ),(3)

where x2 and e1 are independent, g1 is non-constant,and g2 is invertible.
Notation. The following notations are used hereafter.
Suppose that both (2) and (3) hold. Random variables
t1 and z2 and functions h and h1 defined as follows:
t1 , g2−1 (x1 ),h , f1 ◦ g2 ,z2 , f2−1 (x2 ),h1 , g1 ◦ f2 .

That is, h(t1 ) = f1 (g2 (t1 )) = f1 (x1 ), and similarly,
h1 is a function of z2 . Moreover, we let η1 (t1 ) ,log pt1 (t1 ), and η2 (e2 ) , log pe2 (e2 ).1 The following
theorem gives the constraint that pt1 , pe2 , and h mustsatisfy to make both (2) and (3) hold.
Theorem 1 Assume that (x1 , x2 ) can be described byboth of the causal relations given in (2) and in (3).
Further suppose that involved densities and nonlinear
functions pt1 , pe2 , f1 , f2 , g1 , and g2 are third-orderdifferentiable, and that pe2 is positive on (−∞, +∞).
We then have the following equation for every (x1 , x2 )satisfying η2′′ h′ 6= 0:
η1′′′ −

η1′′ h′′
h′

=


η ′′′
− 2η2′′ · h′ h′′ − 2′′ · h′ η1′′
η2

′′2 
h
(4)
+ η2′ · h′′′ − ′ ,h η ′ η ′′′
2 2
η2′′

and h1 depends on η1 , η2 , and h in the following way:1
η ′′ + η2′′ h′2 − η2′ h′′
.
= 1
′
h1
η2′′ h′(5)

See Appendix for its proof. In the special case with f2
and g2 being the identity mapping, the above theorem
becomes Theorem 1 in Hoyer et al. (2009), which was
1
For the sake of conciseness, sometimes we drop thearguments of the functions in the presentation.

given for the additive noise model. In general, it is not
obvious if Theorem 1 holds for a practical problem.
Therefore, below we provide conditions which are easy
to verify, by investigating (4) in Theorem 1. To thisend, we first give some definitions and lemmas.
Definition 2 The density of the continuous variablev, denoted by pv , is said to be log-mixed-linearand-exponential (log-mix-lin-exp) if it is of the form
log pv = c1 ec2 v + c3 v + c4 , where c1 , c2 , c3 , and c4
are constants. Clearly, to make pv a valid density, we
have c1 < 0 and c2 c3 > 0. This type of distributions
includes the Type-1 Gumbel distribution as a special
case when c2 = c3 .
Definition 3 The density pv is said to be a generalized mixture of two exponentials if pv ∝ (c1 ec2 v +c3 ec4 v )c5 , with constants ci , or equivalently, if log pv =d1 v + d2 log(d3 + d4 ed5 v ) + d6 , with constant di .
Definition 4 The density pv is said to be one-sided
asymptotically exponential, if (log pv )′ → c, as v →−∞ or as v → +∞, where c is a non-zero constant.
It is said to be two-sided asymptotically exponential
if (log pv )′ → c1 , as v → −∞, and (log pv )′ → c2 , asv → +∞, where c1 and c2 are non-zero constants.
Clearly, log-mix-lin-exp densities are special cases ofthe one-sided asymptotically exponential densities,
and two-sided asymptotically exponential densities include those generalized mixtures of two exponentials.
Now we give a lemma to characterize some propertiesof a density function or its logarithm.
Lemma 5 Suppose the density pv is differentiable on(−∞, +∞) and monotonic for sufficiently large v. We
have the following. (i) If pv is positive on (−∞, +∞),
then log pv → −∞ as |v| → +∞. (ii) If there
only exists one point, denoted by v = C, satisfying(log pv )′ = 0, then on the support of pv , (log pv )′ > 0
for v < C, and (log pv )′ < 0 for v > C.
This lemma is obvious and proof is skipped. The following lemmas discuss some special and simple solutions to (4) in Theorem 1.
Lemma 6 Under the assumptions made in Theorem 1, if e2 is Gaussian, then h is linear and t1 isalso Gaussian.
Lemma 7 Under the assumptions made in Theorem 1, if function h = f1 ◦ g2 is linear, one of thefollowing must be true:i. t1 and e2 are both Gaussian, and h1 is also linear.

ii. The densities of t1 and e2 are log-mix-lin-exp,
h′1 (z2 ) is strictly monotonic, and h′1 (z2 ) → 0, as
z2 → +∞ or as z2 → −∞.
For proofs, see Appendix. We now give Theorem 8
to consider all the situations in which the PNL causal
model is not identifiable.
Theorem 8 Suppose that η2′′ h′ 6= 0 at every point orthat it is zero at only some discrete points. Under theassumptions made in Theorem 1, we have that pe2 and
pt1 , as well as h, must satisfy one of the five conditions
listed in Table 1.2
The proof is rather long and complicated, so we justgive the basic idea and its outline; see Appendix.
When the nonlinear distortions f2 and g2 are constrained to be the identity mapping, this result also applies to the identifiability of the additive noise model.
We give the following remarks on the situations listed
in Table 1. Situation I in Table 1 is the linear Gaussian case, which is well know to be not identifiable.
Situations II∼V are novel. As discussed in Zhang andHyvärinen (2008), if f1 is constrained to be invertible,
the PNL causal model (2) can be transformed to the
PNL mixing ICA model (Achard & Jutten, 2005), withf1 (x1 ) and e2 considered as sources. The identifiabilityof the former model is then implied by the separabilityof the latter one. Previously, under weak conditions,
it was shown that the PNL mixing ICA model is separable if at most one of the sources is Gaussian (Achard& Jutten, 2005). In Situations II∼V, pe2 is not Gaussian, but the causal model is not identifiable, and consequently the corresponding PNL mixing ICA model
is not separable. This means that the established separability results of the PNL mixing ICA problem havesome flaws and require further investigation.
Below we give some corollaries which follow from Theorem 8. They can be easily exploited to examine if thecausal relation between two variables is unique. The
following assumptions are made in the corollaries.
A1. The data (x1 , x2 ) are generated by the PNL causal
model (2), with f1 and f2 being third-order differentiable.
A2. Densities pe2 and px1 are third-order differentiable, pe2 is positive on (−∞, +∞), and η2′′ =(log pe2 )′′ is zero at most at some discrete points.
The following corollary immediately follows Theorem 8, so its proof is skipped.
2
Note that the identifiability of the PNL causal model(2) depends directly on the distribution of t1 and e2 , instead of the distribution of the observed variables x1 andx2 .

Table 1: All situations in which the PNL causal model is not identifiable.
I
II

pe2
Gaussianlog-mix-lin-exp

pt1 (t1 = g2−1 (x1 ))Gaussianlog-mix-lin-exp

h = f1 ◦ g2
linearlinearIII

log-mix-lin-exp

IV

log-mix-lin-exp

h strictly monotonic,
and h′ → 0, as t1 →+∞ or as t1 → −∞Same as aboveV

generalized mixtureof two exponentials

one-sided asymptotically exponential (butnot log-mix-lin-exp)
generalized mixture oftwo exponentials
two-sided asymptotically exponential

Corollary 9 Suppose that assumptions A1 & A2
hold. If pe2 , the density of the disturbance, is not
Gaussian, nor log-mix-lin-exp, nor a generalized mixture of two exponentials, then the PNL causal model(2) is identifiable.
Corollary 10 considers the identifiability of the PNL
causal model with f1 , the nonlinear effect of the cause,
being non-invertible.
Corollary 10 Suppose that assumptions A1 & A2
hold. If function f1 is not invertible, then the PNL
causal model (2) is identifiable.
For proof, see Appendix. This result is intuitively
appealing, and confirms the finding in Friedman andNachman (2000).

3NONLINEAR ICA-BASEDIDENTIFICATION METHOD

If data (x1 , x2 ) follow the PNL causal model withcausal direction x1 → x2 , from (2), we can see thatthe disturbance e2 , which is independent from x1 , canbe expressed in terms of x1 and x2 :e2 = f2−1 (x2 ) − f1 (x1 ).
This provides a way to verify if x1 → x2 holds according to the PNL causal model, as given below.
Under the hypothesis x1 → x2 , we can estimate thedisturbance e2 by finding functions l1 and l2 such thatê2 = l2 (x2 ) − l1 (x1 ) is independent from x1 . This isthen a constrained nonlinear ICA problem, and can beachieved by minimizing I(x1 , ê2 ), the mutual information between x1 and ê2 (Zhang & Hyvärinen, 2008).
After some simplifications, one can see I(x1 , ê2 ) =−E log pê2 (ê2 ) − E log |l2′ (x2 )| + H(x1 ) − H(x1 , x2 ).
As the last two terms do not depends on l1 (x1 ) andl2 (x2 ), minimizing I(x1 , ê2 ) is equivalent to maximizing E log pê2 (ê2 ) + E log |l2′ (x2 )|. Following Zhang
and Hyvärinen (2008), we use multi-layer perceptrons

Same as above

Remarkh1 also linear
h1 strictly monotonic, and h′1 →0, as z2 → +∞ or as z2 → −∞———(MLP’s) to represent l1 and l2 . The involved parameters can then be learned by gradient-based methods.
Finally, if statistical independence tests, such as thekernel-based test (Gretton et al., 2008), confirm thatê2 is independent from x1 , x1 → x2 is supported by the
PNL causal model, and the learned l1 and l2 provide
an estimate of f1 and f2−1 , respectively.
To find the causal relation between x1 and x2 impliedby the PNL causal model, one needs to test both hypotheses x1 → x2 and x2 → x1 , using the method justdescribed. If exactly one of them holds, the causal relation between x1 and x2 implied by the PNL model has
been successfully found. If neither of them holds, there
is no PNL causal relation between x1 and x2 . If bothhold, the cause and effect could not be distinguishedby the PNL causal model; additional information, suchas the smoothness of the involved nonlinearities, may
help find the causal model with a lower complexity.

4

WITH MORE THAN TWO
VARIABLES

The PNL acyclic causal model (1) is applicable in thecase of more than two variables. When there are only
very few variables, one may use a brute-force search tofind the causal relations, like the estimation of the additive noise model in Hoyer et al. (2009); for each possible acyclic causal structure, represented by a DAG,
we use the nonlinear ICA-based approach to estimate
the corresponding disturbances, and then verify if they
are mutually independent by performing independence
tests. The simplest causal model which gives independent disturbances is preferred. Clearly this approach
may encounter two difficulties. One is that the test of
mutual independence is difficult to do when we have
many variables. The other is that the search spaceof all possible DAG’s increases too rapidly with thevariable number. Consequently, this approach involves
high computational load, and does not scale well with
the number of variables.
A more practical approach to finding the causal re-lations implied by the PNL causal model consistsof the following two steps. We first use conditional
independence-based methods to find the d-separationequivalent class. Next, the PNL causal model is used
to identify the causal directions that cannot be determined in the first step: for each causal structure
contained in the equivalent class, we estimate the disturbances, and determine if this causal structure isplausible, by examining if the disturbance in a variable xi is independent from the parents of xi . Consequently, one avoids the exhaustive search over all
possible causal structures and high-dimensional statistical tests of mutual independence. The validity ofthis approach is supported by the following theorem.
Theorem 11 When fitting variables x1 , ..., xn to the
PNL acyclic causal model (1) with the causal structure represented by the DAG G, the disturbances ei aremutually independent if and only if the causal Markov
condition holds (i.e., each variable xi is independentof its non-descendants conditional on its parents paiin G), and the disturbance ei in xi is independent of
the parents of xi .
Proof is given in Appendix. An important issue inthis approach is how to perform conditional independence tests for variables with nonlinear causal relations. Generally speaking, the results of nonparametric conditional independence tests may be unreliable when the conditional set contains many variables,
due to the curse of dimensionality. Traditionally, in
the implementation of most conditional independencebased causal discovery algorithms, such as PC (Spirtes
et al., 2001) and IC (Verma & Pearl, 1990), it is assumed that the variables are either discrete or Gaussian with linear causal relations. Although this assumption greatly simplifies the difficulty in conditional
independent tests, it usually does not hold in ourcase which involves nonlinear causal effects and nonGaussian variables. Alternatively, in our case, one may
simplify the conditional independence test procedureby making use of the particular structure of the PNL
causal model. This is out of the scope of this paper
and left for future work.

5

SIMULATIONSThe PNL causal model has been applied for causal discovery of two variables with real-world data in Zhang
and Hyvärinen (2008). It successfully identified thecausal directions for all eight data sets in the “Causeeffect pairs” task (Mooij et al., 2008) included in
Causality Challenge #2, without any backgroundknowledge. This demonstrated the practical usefulness of this model and the identification method forsome real-world problems. Here we conduct simula-tions to verify the non-identifiable conditions of thismodel with two variables, given in Theorem 8. In particular, Situation I in Table 1 is well known to be not
identifiable, but the others are completely new and interesting. Here we use illustrative examples to show
the non-identifiability in Situations II and V.
5.1

ON SITUATION II IN TABLE 1

In Situation II in Table 1, function h = f1 ◦ g2 is linear and the densities of t1 and e2 are log-mix-lin-exp(given by (10) and (9) in Appendix, respectively), and
the PNL causal model is not identifiable. To illustrate that, one just needs to confirm that t1 and z2
are independent of e2 and e1 , respectively; independence between t1 and e2 means that the causal relation x1 → x2 holds, as described by (2); in addition, if
z2 and e1 given by (6) and (7) (see Appendix) are also
independent, one can see that x2 → x1 also holds, andthat the PNL model is not identifiable.
Given c1 , c2 , c3 , and c5 in (10) and (9), one can find theconstants c4 and c7 such that pt1 = eη1 and pe2 = eη2
are valid densities. We then generate random numbers
as realizations of the independent variables t1 and e2 ,by applying the inverse of their cumulative distribution functions (CDF’s) to independent and uniformly
distributed variables. According to (7), variable z2 involved in the causal relation x2 → x1 is calculated asz2 = e2 + h(t1 ). The other variable, e1 , can be foundaccording to (6), if h1 is known. According to (11), we
have

. 1
′ c1 z2 −c5 +c2.
+heh′1 = ec1 z2 −c5 +c2h′This yields h1 = h′1c1 log | h1′ +h′ ec1 z2 −c5 +c2 |+c8 , where
c8 is an arbitrary constant. Substituting h1 into (6),
we can find e1 .
In our simulation study, we set h(t1 ) = −t1 , c1 = 0.3,
c2 = −1, c3 = 1, c5 = −1 c8 = 0, and variables t1 ande2 were made zero-mean. 2000 samples were drawn.
Fig. 1 (a) gives the scatter plot of t1 and e2 , as wellas their marginal histograms. Finally, the scatter plotof z2 and e1 is shown in Fig. 1 (b). To verify if z2
and e1 are independent, we conducted the kernel-based
statistical independence tests (Gretton et al., 2008).
The left part of Table 2 presents the independence
test results for the pairs (t1 , e2 ) and (z2 , e1 ). For bothpairs, The independence hypothesis is accepted at thesignificance level α = 0.05. That is, in this situation,
the PNL causal model cannot distinguish the causefrom effect, if we do not have further knowledge about
the data generating process.
We further generated observed data (x1 , x2 ) from variables t1 and e2 , aiming to verify if the nonlinear ICAbased identification method can detect both causal

5

6
4
2

e2

0
−2

e1

2

1
0

0

−4

−2

−6

−5

−4

−8
−5

t

0

−10

5

0z

10

2

1

−10
−5

(a)

.

(b)

Figure 1: (a) Scatter plot of t1 = g2−1 (x1 ) and e2 inSimulation 1. (b) That of z2 = f2−1 (x2 ) and e1 .
Table 2: Result of kernel-based independence tests inSimulation 1, at the significance level α = 0.05. The
independence hypothesis is accepted in all the fourcases. See main text for further explanations.
Threshold
Statistic

t1 &e20.56690.3497z2 &e10.56650.1226x1 &ê2
0.43470.1430x2 &ê1
0.57780.1842

directions. Note that the identifiability of the PNL
causal model only depends the distributions of t1 ande2 and function h(t1 ) = f1 (g2 (t1 )). In situation II in
Table 1, f1 , f2 , and g2 can be arbitrary, given thath = f1 ◦ g2 satisfies the condition.3 In this simula1/3
tion, we let g2 (t1 ) = t1 /2 + t1 ; since h = f1 ◦ g2 ,f1 was then constructed as f1 = h ◦ g2−1 . we let
f2 (h(t1 ) + e2 ) = tanh((h(t1 ) + e2 )/8). Finally, (x1 , x2 )were constructed from t1 and e2 by x1 = g2 (t1 ) and(2). We applied the nonlinear ICA-based method
given in Section 3, and obtained ê2 and ê1 , which are
the estimate of the disturbance under the hypotheses
x1 → x2 and x2 → x1 , respectively. The right part ofTable 2 reports the results of independence tests forthe pairs (x1 , ê2 ) and (x2 , ê2 ). Both pairs accept theindependence hypothesis, meaning that both the estimated PNL causal model with x1 → x2 and that withx2 → x1 can explain the data.
5.2

pt1
pe2
h(t1)
h1(z2)

1.5

4

ON SITUATION V IN TABLE 1

We next verify Situation V in Table 1. In the proof of
Theorem 8, Solution 1 to the functional equation (12),given by (14), with the constraints A3 = −A2 andA4 6= 0, results in this situation. Generally speaking,in this situation, the solutions of η1 and h cannot be
expressed in terms of elementary functions. We thenfound numerical solutions of η1′ , h′ , and η2′ to the overdetermined system of ordinary differential equations
3
Note that if f1 is not invertible, the causal model isalways identifiable; see Corollary 10.

η1(t1)’0.5

η2(e2)’

0h(t1)’−0.5

h1(z2)’
0

5

(a)

10

−1
−5

0

5

10

(b)

Figure 2: An example of Situation V in Table 1.
(a) The curves of η1 (t1 )′ , η2 (e2 )′ , h(t1 )′ , and h1 (z2 )′ .
η1 (t1 )′ , η2 (e2 )′ , and h(t1 )′ were obtained by solving thesystem of ODE’s in (14) numerically. h1 (t1 )′ was givenby (5). (b) The curves of pt1 = eη1 , pe2 = eη2 , h(t1 ),and h1 (z2 ). See main text for detailed explanations.
(ODE’s) given in (14) using MATLAB. We let A1 =2, A2 = 0.4, A4 = 1, c6 = −10, c7 = 1.6, and theinitial conditions are η1′ (t1 = 0) = 0, η2′ (e2 = 0) = 0,h′ (t1 = 0) = 3. The solutions of η1′ , h′ , and η2′ aredepicted in Fig. 2(a). One can see that h′ is clearlyclose to zero for large t1 , consistent with Situation V.
Furthermore, we set h(t1 = 0) = 0, and obtained thesolutions of h; also noting that the total integral of
pt1 = eη1 and pe2 = eη2 is one, finally we obtained thesolutions of η1 , h, and η2 , as plotted in Fig. 2(b).
We used radial basis networks (RBF’s) to learn theprobability density functions (PDF’s) of t1 and e2 fromtheir numerical values, and then drew random samplesof t1 and e2 by applying the inverse of their CDF’s
to independent and uniformly distributed variables.
Fig. 3(a) gives the scatter plot of t1 and e2 used inthis simulation and their marginal histograms. After
that, we aimed to find variables z2 and e1 , which are
involved in the causal relation x2 → x1 , by making useof the transformation given in (6) and (7). Apparently,
z2 = h(t1 ) + e2 . In order to find variable e1 , one needs
to find h1 (z2 ). We first calculated h1 (z2 )′ according
to (5). It is also given in Fig. 2(a), and the solutionof h1 with the initial condition h1 (0) = 0 is shown in
Fig. 2(b). Variable e1 is then found as e1 = t1 +h1 (z2 ).
The scatter plot of z2 and e1 is shown in Fig. 3(b).
Results of statistical independence tests, given in theleft part of Table 3, confirm statistical independence
between t1 and e2 and that between z2 and e1 .
Analogously to Simulation I, we constructed observable data (x1 , x2 ) from t1 , e2 , and h given above. To dothat, we let g2 (t1 ) = t1 +t31 and f2 (h(t1 )+e2 ) = log(6+h(t1 ) + e2 ), and f1 was constructed as f1 = h ◦ g2−1 .
The nonlinear ICA-based method was then applied tothe data (x1 , x2 ). As shown in the right part of Table 3, under both hypotheses x1 → x2 and x2 → x1 ,

References
1

6

Achard, S., & Jutten, C. (2005). Identifiability of postnonlinear mixtures. IEEE Signal Processing Letters, 12,423–426.

0.5

2

1

4

0

−0.5

e

e2

0

−2

−1

−4
−1

t1

0

−4 −2

1

0

z2

4

6

8

2

(a)

Friedman, N., & Nachman, I. (2000). Gaussian process
networks. Proc. 16th Conference on Uncertainty in Artificial Intelligence (UAI) (pp. 211–219). Morgan Kaufmann.

(b)

Figure 3: (a) Scatter plot of t1 and e2 , which are independent variables with the densities given in Fig. 2(b).
(b) Scatter plot of z2 and e1 , which were obtained bymaking use of the transformation in (6) and (7).
Table 3: Result of independence tests in Simulation 2
at the significance level α = 0.05. The independence
hypothesis is accepted in all the four cases.
Threshold
Statistic

t1 &e2
0.54830.1641z2 &e10.1561
0.0669x1 &ê2
0.50810.3011x2 &ê1
0.69120.4142the assumed cause is independent from the estimate ofthe disturbance. Consequently, both causal directions
can explain the data, and in this example the PNL
causal model is not identifiable.

6

Braun, M. (1993). Differential equations and their applications: An introduction to applied mathematics. New
York: Springer-Verlag. 4th edition.

CONCLUSION

We have investigated the identifiability of the postnonlinear (PNL) causal model in the two-variable case,and a practical identification method based on nonlinear function approximations and independence tests
was discussed. Our results show that this model isgenerally identifiable, and consequently can be used
to distinguish the cause from effect. All the particular situations in which this model is not identifiable
were reported in Theorem 8, and some of them wereverified and illustrated by simulations. For the situation with more than two variables, we showed that
it is not necessary to apply the PNL causal model toall of the given variables directly, which becomes intractable as the variable number increase; instead, tofind the whole causal structure, one can apply thiscausal model to the Markov equivalent class and test
if the estimated disturbance is independent from theparents associated with the same variable.
Acknowledgements
The authors would like to thank Patrik Hoyer for helpful comments and suggestions.

Gretton, A., Fukumizu, K., Teo, C., Song, L., Schölkopf,
B., & Smola, A. (2008). A kernel statistical test of independence. NIPS 20 (pp. 585–592). Cambridge, MA:
MIT Press.
Hoyer, P., Janzing, D., Mooji, J., Peters, J., & Schölkopf,
B. (2009). Nonlinear causal discovery with additive noisemodels. NIPS 21. Vancouver, B.C., Canada.
Hyvärinen, A., Karhunen, J., & Oja, E. (2001). Independent component analysis. John Wiley & Sons, Inc.
Lin, J. (1998). Factorizing multivariate function classes.
NIPS 10 (pp. 563–569). Cambridge, MA: MIT Press.
Mooij, J., Janzing, D., & Schölkopf, B. (2008).
Distinguishing between cause and effect.
URL
http://www.kyb.tuebingen.mpg.de/bs/people/
jorism/causality-data/.
Pearl, J. (2000). Causality: Models, reasoning, and inference. Cambridge: Cambridge University Press.
Polyanin, A., & Zaitsev, V. (2003). Handbook of exact solutions for ordinary differential equations. Boca Raton:
Chapman & Hall/CRC Press. 2nd edition.
Polyanin, A., & Zaitsev, V. (2004). Handbook of nonlinear
partial differential equations (supplement s.4.4). Boca
Raton: Chapman & Hall/CRC Press.
Shimizu, S., Hoyer, P., Hyvärinen, A., & Kerminen, A.
(2006). A linear non-Gaussian acyclic model for causal
discovery. Journal of Machine Learning Research, 7,
2003–2030.
Spirtes, P., Glymour, C., & Scheines, R. (2001). Causation,
prediction, and search. Cambridge, MA: MIT Press. 2ndedition.
Sun, X., Janzing, D., & Schölkopf, B. (2008). Causal reasoning by evaluating the complexity of conditional densities with kernel methods. Neurocomputing, 71, 1248–
1256.
Verma, T., & Pearl, J. (1990). Equivalence and synthesisof causal models. Proceedings of 6th Conference on Uncertainty in Artificial Intelligence (UAI) (pp. 220–227).
Boston, MA: Morgan Kaufmann.
Zhang, K., & Hyvärinen, A. (2008).
Distinguishingcauses from effects using nonlinear acyclic causal models. Presented at NIPS 2008 Workshop on Causality. URL http://www.cs.helsinki.fi/u/ahyvarin/
papers/Zhang09NIPSworkshop.pdf.

APPENDIX: SOME PROOFS
Proof of Theorem 1: We prove this theorem using the linear separability of the logarithm of the joint
density of independent variables, which states the factthat for a set of independent random variables whose
joint density is twice differentiable, the Hessian of thelogarithm of their density is diagonal everywhere (Lin,1998). Since g2 is invertible, the independence between x1 and e2 is equivalent to that between t1 ande2 . Similarly, the independence between x2 and e1 is
equivalent to that between z2 and e1 . Combining thetwo causal models (2) and (3), one can see that thetransformation from (z2 , e1 ) to (t1 , e2 ) ist1
e2

= h1 (z2 ) + e1 ,= z2 − h(t1 ).

(6)
(7)

Denote by J the Jacobian matrix of this transformation. One can see that |J| = 1. Denote by p(z2 ,e1 )the joint density of (z2 , e1 ). We then have pt1 · pe2 =p(z2 ,e1 ) /|J| = p(z2 ,e1 ) , so, log p(z2 ,e1 ) = η1 (t1 ) + η2 (e2 ).
One can find the (1,2)-th entry of the Hessian matrix∂ 2 log p(z2 ,e1 )∂t1of log p(z2 ,e1 ) w.r.t. (z2 , e1 ):= η1′′ ∂z−∂e1 ∂z2
2′ ′′ ′′′ ′2 ′′′ ′′ ′′ ∂t1′′ ′′′ ′ ∂e2
η2 h ∂z2 − η2 h ∂z2 = η1 h1 − η2 h + η2 h h1 − η2 h h1 .

The independence between z2

and e1 implies∂ 2 log p(z2 ,e1 )= 0 for every possible (z2 , e1 ). That∂e1 ∂z2′′ ′′′ ′is, η1 h1 − η2 h + η2′′ h′2 h′1 − η2′ h′′ h′1 = 0. From thisequation one can see that h′1 = 0 implies η2′′ h′ = 0.
Consequently, the points which satisfy η2′′ h′ 6= 0 also
make h′1 6= 0. For such points, dividing both sidesof this equation by h′1 η2′′ h′ finally leads to (5). Furthermore, since h1 is a functions
 . of z2 and does not∂e1 = 0. According
depend on e1 , we have ∂ h1′
 ′′ ′′ ′21 ′ ′′ η +η h −η h/∂e1 = 0, whichto (5), we have ∂ 1 2η′′ h′ 22gives 2η2′′2 h′2 h′′ − η2′ η2′′ h′ h′′′ + η2′′ η1′′′ h′ − η2′ η2′′′ h′2 h′′ +η2′ η2′′ h′′2 + η2′′′ η1′′ h′2 − η2′′ η1′′ h′′ = 0. For the points satisfying η2′′ h′ 6= 0, we divide both sides of the above
equation by η2′′ h′ . After some simplifications, (4) is
obtained. 
Proof of Lemma 6: Gaussianity of e2 implies that
η2′′′ ≡ 0 and that η2′′ is constant. (4) then reduces toη ′′ h′′h′′2 
η1′′′ − 1 ′ + 2η2′′ · h′ h′′ = η2′ · h′′′ − ′ .
(8)hhSince the left-hand side does not depend on e2 and η2′′′2
is a function of e2 , we have h′′′ − hh′ = 0, which gives′′2h′′′ /h′ − hh′2 = 0. That is, (h′′ /h′ )′ = 0, so h′′ /h′ = c1 .
If c1 = 0, h is linear, and h′′ = h′′′ = 0. (8) then yields
η1′′′ = 0, meaning that t1 is Gaussian.
Otherwise, h′ = ±ec1 t2 +c2 . Moreover, the left-hand
η ′′ h′′η ′′′side of (8) must also be zero, which means h1′ − h1 ′2 +R  η1′′′
η1′′ h′′2η2′′ h′′ = 0. By integration, it givesh′ − h′2 +
2η2′′ h′′ dt1 = η1′′ /h′ + 2η2′′ h′ = c3 , so η1′′ = −2η2′′ h′2 +c3 h′ = −2η2′′ e2c1 t1 +2c2 ±c3 ec1 t1 +c2 . Consequently η1 =−η2′′ 2c1 t1 +2c2
± cc23 ec1 t1 +c2 +c4 t1 +c5 . Noting that η2′′ ise2c21
1a negative constant, we can see that η1 → +∞ whenc1 t1 → +∞. This contradicts Lemma 5. Thus, c1
must be zero, h is linear, and t1 is Gaussian. 
Proof of Lemma 7: Note that η2′′ (as well as η1′′ ) is
not constantly zero, as a density function of e2 could
not be proportional to ec1 e2 +c2 . When h is linear, h′is constant and h′′ = h′′′ = 0, and (4) becomes
η ′′′
η ′′′
η ′′′
η1′′′ = − 2′′ · h′ η1′′ , i.e., 1′′ = −h′ · 2′′
η2
η1
η2
η ′′′

η ′′′

1

2

Since η1′′ and η2′′ depends only on t1 and e2 , respecη ′′′tively, the equation above implies that both η1′′ and
η2′′′
η2′′

η ′′′

η ′′′

1

are constants. Let η2′′ = c1 , we have η1′′ = −h′ c1 .
2
1
If c1 = 0, clearly both e2 and t1 are Gaussian, andproposition (i) holds.
Otherwise, we have log |η2′′ | = c1 e2 + c2 . Consequently
η2′′ = ±ec1 e2 +c2 and η2 = ± c12 ec1 e2 +c2 + c3 e2 + c4 . If
1

η2 = c12 ec1 e2 +c2 + c3 e2 + c4 , clearly η2 → +∞ when1sgn(c1 ) · e2 → +∞, which contradicts Lemma 5. One
can verify that when c1 c3 > 0, the solution
1
η2′′ = −ec1 e2 +c2 , or η2 = − 2 ec1 e2 +c2 + c3 e2 + c4 (9)c1corresponds to valid densities. Analogously we have′′1
η1′′ = −e−h c1 t1 +c5 , or η1 = − ′2 2 e−h c1 t1 +c5 +c6 t1 +c7 .
h c1
(10)
Thus the densities of t1 and e2 are both log-mix-linexp.

Combining (5), (9) and (10), and recalling that h islinear, we have
h′1

=
=

1
η2′′ h′= ′′ ′′ ′′′′′η1 + η2 h′2 − η2′ h′′η1 /(η2 h ) + h′1
e−c1 z2 +c5 −c2 /h′ + h′

.

(11)

One can see that h′1 → 0, as sgn(c1 ) · z2 → +∞, andas sgn(c1 ) · z2 → −∞, h′1 → h′ , which is a constant.
Therefore Lemma 7 is true. Outline of Proof of Theorem 8: Eq. (4) can be
re-written as a bilinear functional equation (Polyanin& Zaitsev, 2004) of the form:Φ1 (t1 )Ψ1 (e2 ) + Φ2 (t1 )Ψ2 (e2 )+ Φ3 (t1 )Ψ3 (e2 ) + Φ4 (t1 )Ψ4 (e2 ) = 0, (12)

whereη1′′ h′′h′,
Φ1 (t1 ) = η1′′′ −Φ3 (t1 ) = h′ h′′ ,Ψ1 (e2 ) = −1,
η ′ η ′′′
Ψ3 (e2 ) = 2η′′2 − 2η2′′ ,2′′2

Φ2 (t1 ) = h′′′ − hh′ ,Φ4 (t1 ) = h′ η1′′ ,Ψ2 (e2 ) = η2′ ,η ′′′Ψ4 (e2 ) = − η2′′ .

(13)

2

Note that functionals Φi (t1 ) and Ψi (e2 ) depend only
on t1 and e2 , respectively. We then find all possible
situations in which (12) holds.
Clearly Φ4 (t1 ), Ψ1 (e2 ), and Ψ2 (e2 ) in (13) are not constantly zero. We first consider some simple cases ofthe solutions to (12), with Φ3 (t1 ) ≡ 0, Ψ4 (e2 ) ≡ 0, orΦ2 (t1 ) ≡ 0. These cases either have no valid solutionsfor η1 , η2 , and h, or the solutions are covered by thesituations in Table 1.
When none of the functions mentioned above is constantly zero, it can be shown that the functional
equation of the form (12) has three solutions, whichare (Polyanin & Zaitsev, 2004)
Solution 1:
Φ1 = A1 Φ3 + A2 Φ4 ,
Ψ3 = −A1 Ψ1 − A3 Ψ2 ,

Φ2 = A3 Φ3 + A4 Φ4 ,
Ψ4 = −A2 Ψ1 − A4 Ψ2 ,(14)

Solution 2:
Φ1 = B1 Φ3 , Φ2 = B2 Φ3 , Φ4 = B3 Φ3 ,
Ψ3 = −B1 Ψ1 − B2 Ψ2 − B3 Ψ4 ,Dn on which f ′ = 0 holds at most at some discrete
points. Without loss of generality, we assume that Dn
is on the right side of D0 . As h = f1 ◦ g2 , we haveh(t1 = C − )′ = 0. Suppose that both x1 → x2 andx2 → x1 can explain the data. In all the situations
listed in Table 1, h′ 6= 0. Hence, h(t1 = C + )′ 6= 0. As
h(t1 = C + )′ 6= h(t1 = C − ), h is not differentiable at
t1 = C, which causes a contradiction. So in this casethe causal direction between x1 and x2 given by the
PNL causal model is also unique. Proof of Theorem 11: The necessity part is obvious,and below we prove the sufficiency part, which statesthat disturbances e1 , ..., en are mutually independent,
or equivalently, x1 , ..., xn follow the PNL causal model
represented by G, if the causal Markov condition holds
and the disturbance ei in the variable xi is independent−1of pai . Let zi , fi,2(xi ). As the causal relations areacyclic, we can arrange xi in an order such that nolater variable causes any earlier one. Without loss ofgenerality, we assume that (x1 , x2 , ..., xn ) is one of suchorders. For any i = 1, ..., n, since xi = fi,2 (zi ), we have′p(xi |pai ) = p(zi |pai )/|fi,2(zi )|. Therefore,
H(ei ) ≥ H(ei |pai )= H(zi |pai ) = −E{log p(zi |pai )}′= −E{log p(xi |pai )} − E{log |fi,2(zi )|}′= H(xi |pai ) − E{log |fi,2 (zi )|}′≥ H(xi |x1 , ...xi−1 ) − E{log |fi,2(zi )|},(18)(15)

Solution 3:
Ψ2 = C1 Ψ1 , Ψ3 = C2 Ψ1 , Ψ4 = C3 Ψ1 ,
Φ1 = −C1 Φ2 − C2 Φ3 − C3 Φ4 ,(16)

where Ai , Bi , and Ci are arbitrary constants. Eachpossible solution given above is an over-determined
system of ordinary differential equations (ODE’s). In
most cases, solutions of η1 , η2 , and h to the system,
which also satisfy the properties given in Lemma 5,
can be found in closed form. In the remaining cases,
we use the way of phase portrait (Braun, 1993) to analyze the behavior of the solutions. We consider all
possible cases, and the solutions are always covered by
the five situations enumerated in Table 1. 
Proof of Corollary 10: First, consider the casewhere f1′ 6= 0 at every point or it is zero at only somediscrete points. As h = f1 ◦ g2 and g2 is invertible,h′ = 0 holds at most at some discrete points. All thesituations in which the PNL causal model is not identifiable are given in in Table 1. One can see that in allthese situations, h = f1 ◦ g2 is invertible. When f1 is
not invertible, no matter what g2 is, it is impossible to
make h = f1 ◦ g2 invertible, i.e., none of the conditionsin Table 1 holds. Consequently, in this case the causal
direction between x1 and x2 is uniquely identified.
Next, consider the case where f ′ ≡ 0 on the domain
D0 ∈ R. Let C be a point between D0 and the domain(17)

where H(·) denotes the entropy, the equality in (17)holds if and only if ei is independent from pai , and theequality in (18) holds if and only if the causal Markov
condition holds, i.e., elements of {xk |xk 6∈ pai , 1 ≤k ≤ i − 1} are independent of xi given pai . Taking thesummation of (18) over i gives
X
H(ei )i≥

X

H(xi |x1 , ..., xi−1 ) −i= H(x1 , ..., xn ) −X′
E{log |fi,2(zi )|}iX′
E{log |fi,2(zi )|}.

(19)i

In addition, since ei does not depend on xj (j > i), theJacobian matrix of the transformation from (x1 , ...xn )to (e1 , ...en ) is lower-triangular, with the (i, i)th en′
try being 1/fi,2(zi ). Consequently, the absolute valueof the determinant of this Jacobian matrix is |J| =′[Πi fi,2(zi )]−1 . Recalling (19), when the causal Markov
condition holds and ei is independent of pai , themutualinformation of e1 , ...eP
P n is then I(e1 , ..., en ) =
H(e
,
...,
e
)
=
− [H(x1 , ..., xn ) +1ni H(ei ) −i H(ei ) P
P′
log |J|] = i H(ei ) − H(x1 , ..., xn ) + i log |fi,2(vi )| =0. That is, e1 , ..., en are mutually independent. 