Causal discovery of linear acyclic models with arbitrary distributions

Patrik O. Hoyer
Aapo Hyvärinen
Helsinki Institute forInformation Technology &Department of Computer Science
University of Helsinki
FinlandRichard Scheines
Peter Spirtes
Joseph Ramsey
Department of Philosophy
Carnegie Mellon University
Pittsburgh, PA, USA

Abstract
An important task in data analysis is thediscovery of causal relationships between observed variables. For continuous-valued data,linear acyclic causal models are commonlyused to model the data-generating process,and the inference of such models is a wellstudied problem. However, existing methods
have significant limitations. Methods basedon conditional independencies (Spirtes et al.
1993; Pearl 2000) cannot distinguish between
independence-equivalent models, whereas approaches purely based on Independent Component Analysis (Shimizu et al. 2006) areinapplicable to data which is partially Gaussian. In this paper, we generalize and combine the two approaches, to yield a methodable to learn the model structure in manycases for which the previous methods provide answers that are either incorrect or are
not as informative as possible. We give exact graphical conditions for when two distinct
models represent the same family of distributions, and empirically demonstrate the powerof our method through thorough simulations.

1

INTRODUCTION

In much of science, the primary focus is on the discovery of causal relationships between quantities of interest. The randomized controlled experiment is geared
specifically to inferring such relationships. Unfortunately, in many studies it is unethical, technically extremely difficult, or simply too expensive to conduct
such experiments. In such cases causal discovery must

Gustavo Lacerda
Machine Learning Department
Carnegie Mellon University
Pittsburgh, PA, USA
Shohei Shimizu
Osaka University
Japan

be based on uncontrolled, purely observational data
combined with prior information and reasonable assumptions.
In cases in which the observed data is continuousvalued, linear acyclic models (also known as recursive
Structural Equation Models) have been widely used
in a variety of fields such as econometrics, psychology,sociology, and biology; for some examples, see (Bollen
1989). In much of this work, the structure of the models has been assumed to be known or, at most, only
a few different models have been compared. During
the past 20 years, however, a number of methods havebeen developed to learn the model structure in an unsupervised way (Spirtes et al. 1993; Pearl 2000; Geiger
and Heckerman 1994; Shimizu et al. 2006). Nevertheless, all approaches so far presented have either required distributional assumptions or have been overly
restricted in the amount of structure they can inferfrom the data. In this contribution we show how to
combine the strenghts of existing approaches, yielding
a method capable of inferring the model structure inmany cases where previous methods give incorrect or
uninformative answers.
The paper is structured as follows: Section 2 precisely
defines the models under study, and Section 3 discusses existing methods for causal discovery of suchmodels. In Section 4 we formalize the discovery problem and give exact theoretical results on identifiability.
Then, in Section 5 we introduce and analyze a method
termed PClingam that combines the strenghts of existing methods and overcomes some of their weaknesses,
and is, in the limit, able to estimate all identifiable
aspects of the underlying model. Section 6 provides
empirical demonstrations of the power of our method.
Finally, Section 7 maps out future work and Section 8
provides a summary of the main points of the paper.

2

aLINEAR MODELSIn this paper, we assume that the observed data has
been generated by the following process:1. The observed variables xi , i = {1 . . . n} can bearranged in a causal order, such that no later variable causes any earlier variable. We denote sucha causal order by k(i). That is, the generatingprocess is recursive (Bollen 1989), meaning it canbe represented graphically by a directed acyclic
graph (DAG) (Pearl 2000; Spirtes et al. 1993).
2. The value assigned to each variable xi is a linear
function of the values already assigned to the earlier variables, plus a ‘disturbance’ (noise) term ei ,and plus an optional constant term ci , that isX
bij xj + ei + ci ,(1)
xi =k(j)<k(i)where we only include non-zero coefficients bij inthe equation.
3. The disturbances ei are all continuous randomvariables with arbitrary densities pi (ei ), and the ei
are
Q independent of each other, i.e. p(e1 , . . . , en ) =i pi (ei ).
This formulation neither requires the disturbances tobe normally distributed nor does it require them tohave non-Gaussian (non-normal) densities. In general,
some of the distributions can be Gaussian and some
not, and we do not a priori know which are which.
We assume that we are able to observe a large number of data vectors x (which contain the variables xi ),and each data vector is generated according to theabove described process, with the same causal orderk(i), same coefficients bij , same constants ci , and the
disturbances ei sampled independently from the same
distributions. Note that the independence of the disturbances implies that there are no unobserved confounders (Pearl 2000). Spirtes et al. (1993) call thisthe causally sufficient case.
Finally, we assume that the observed distribution isfaithful to the generating graph (Spirtes et al. 1993),
i.e. the model is stable in the terminology of Pearl(2000). If the model parameters are in some sense
randomly generated, this is not a strong assumption,
as violations of faithfulness have Lebesgue measure 0
in the space of the linear coefficients.
An example of such a model is given in Figure 1a. Notethat the full model consists of a directed acyclic graphover the variables, the connection strenghts bij , theconstants ci , and the densities pi (ei ). In this example
we have chosen ci = 0 for all i, so these are not shown.

b
x

c

d

x

x

x

x

x

y

y

y

y

y

z

z

z

z

z

3

y
-2

z

Figure 1: An example case used to illustrate the concepts described in Sections 2–4. (a) A linear, acycliccausal model for x, y and z. The data is generated as
x := ex , y := 3x + ey , and z := −2y + ez , with ex andey drawn from Gaussian distributions and ez from a
non-Gaussian distribution, and ex , ey and ez are all
mutually independent. Note that we show variables
with Gaussian disturbances using circles whereas variables with non-Gaussian disturbances are marked bysquares. (b) The three directed acyclic graphs overx, y and z which all entail the same conditional independence relationships as the generating model. (c)The three DAGs in (b) succintly represented as a dseparation-equivalence pattern. (d) The distributionequivalence pattern of the original model.

3

EXISTING METHODS

Given our set of data vectors x, to what extent can
we estimate the data generating process? Obviously,
if the number N of data vectors is small, estimation
may be quite unreliable. Therefore we will here mainly
focus on the theoretical question: To what extent (andwith what methods) can we identify the true model inthe limit as N → ∞?
The most well-known approach to inference of thistype of causal networks is based on (conditional) independencies between the variables (Spirtes et al. 1993;
Pearl 2000). When, as in our case, there are assumed
to be no hidden confounding variables and no selection
bias, one can in the large-sample limit identify the setof networks which represent the same independencies
as the true data generating model. To illustrate, in
Figure 1b we show all three DAGs which imply theset of independencies produced by the true model.
This set is known as the d-separation-equivalence class,and is often represented in the form of a d-separationequivalence pattern: a partly directed graph in which
undirected edges represent edges for which both directions are present in the equivalence class (Spirtes et al.
1993), as illustrated in Figure 1c. We want to emphasize that, using conditional independence information
alone, it is impossible to distinguish between members
inside a d-separation-equivalence class because these(by definition) represent the same set of conditional

independencies between the observed variables.
Fortunately, in many cases there is additional information available that can be used to further distinguishbetween different DAGs. In particular, it can be shown(Shimizu et al. 2006) that if all (or all but one) distributions of the error variables are non-Gaussian, it isin fact possible to identify the complete causal model,
including all the parameters. This is possible using
a method based on Independent Component Analysis (ICA) (Hyvärinen et al. 2001). Unfortunately,
however, when two or more disturbances are Gaussian the standard method based on ICA will fail. As
an extreme example, when all disturbances are Gaussian, standard ICA-based methods return nonsense
and are not even able to find the correct d-separationequivalence class.
These considerations raise the question of whether it
is possible to combine the methods so as to obtainrobustness with respect to Gaussian distributions but
not forgo the possibility of identifying the full model infavourable circumstances. Indeed, such a combinationis possible and is presented in Section 5. Here, we simply note that the naı̈ve solution of first running sometest and then selecting one of the two methods, will not
be optimal. Consider, for instance, our example modelin Figure 1a. Because there is more than one Gaussian error variable the standard ICA-based method(Shimizu et al. 2006) is not applicable, and hence one
would have to settle for the d-separation-equivalenceclass (Figure 1c) given by independence-based methods. However, as we show in the next section, in thisexample we can actually reject one of the DAGs in
the equivalence class and hence obtain a smaller set of
possible generating models.

In general, an ngDAG D is instantiated by many different models M which differ in their connection
strengths bij as well as in their distributions pi (ei ).
Next, we define the important concept of distributionequivalence between ngDAGs, which defines to what extent it is possible to infer the ngDAG which represents
the true data generating causal model, from observational data alone.
Definition 3 Two ngDAGs D1 and D2 are
distribution-equivalent if and only if for any linear acyclic causal model M1 which instantiates D1
there exists an instantiation M2 of D2 which yields
the same joint observed distribution as M1 , and viceversa.
Distribution-equivalence partitions the set of ngDAGs
into distribution-equivalence classes, and these may be
represented using simplified graphs:
Definition 4 An ngDAG pattern representing an
ngDAG D is a mixed graph (consisting of potentiallyboth directed und undirected edges), obtained in the following way:
1. Derive the d-separation-equivalence pattern corresponding to the DAG in D
2. Orient any unoriented edges which originate from,or terminate in, a node positively marked in ng of
D, in the orientation given by the DAG in D
3. Finally, orient any edges which follow from theorientations given in the previous step and dseparation-equivalence, according to the rules derived by Meek (1995).

DISTRIBUTION-EQUIVALENCE

We say that a mixed graph is an ngDAG pattern if it
represents some ngDAG.

First, we need to extend a DAG object to include information on the non-Gaussianity of associated disturbance variables.

An ngDAG pattern is similar in many respects to dseparation-equivalence patterns. For example, we havethe following result:Definition 1 An ngDAG is a pair (G, ng) where G is adirected acyclic graph over a set of variables V and ng
is a binary vector of length |V |, each element of whichis associated with one of the variables of V .

Lemma 1 An ngDAG pattern is a chain graph.

4

Definition 2 We say that a linear acyclic causal
model M instantiates an ngDAG D (alternatively, D
represents M ) if and only if the directed acyclic graphassociated with M is equal to that specified in D, andfurther if the set of variables with non-Gaussian disturbance variables in M is equal to the set of positive
entries in the binary vector specified in D.

The proof is given in the Appendix.
Our main result connects ngDAG patterns with
distribution-equivalence in mixed Gaussian and nonGaussian models in the same way that d-separationequivalence patterns are associated with distributionequivalence in purely Gaussian models:Theorem 1 Two ngDAGs are distribution-equivalent
if and only if they are represented by the same ngDAG
pattern.

The proof of this theorem is provided in the Appendix.
The important point is that we now know exactly
which models are indistinguishable from each other on
the basis of observational data alone.

(c) Calculate the corresponding ICA objectivefunctionX
2
Uf =(E{f (ei )} − k)(2)

As a simple illustration, in Figure 1d we show thengDAG pattern representing the ngDAG corresponding
to the generating model of Figure 1a. Note thatthe ngDAG pattern is more informative than the dseparation-equivalence pattern of Figure 1c. Nevertheless, there are still two ngDAGs (leftmost two in
Figure 1b) which cannot be distinguished based on
non-experimental data.

where k is the expected value of f applied toa zero-mean, unit variance Gaussian variable,
i.e. k = E{f (g)}, g ∼ N (0, 1). In the ICA
literature, many different choices of f havebeen utilized; here we suggest simply taking
the absolute value function f (ei ) = |ei |, giving
2
X
p
(3)
U=
E{|ei |} − 2/π

Henceforth in the paper we shall use the terms
ngDAG pattern and distribution-equivalence pattern interchangably.

5

PC-LINGAM

Although an important goal in this study was to lookat the theoretical aspects of identifying DAGs in mixed
Gaussian / non-Gaussian acyclic linear causal models, an equally significant objective is to give a practical method with which to infer models from a finite data set. Although there are a number of possible approaches, we here give a simple combinationof independence-based techniques and the ICA-based
method. The method, termed PClingam, consists ofthree steps:
1. Use methods based on conditional independence
tests to estimate the d-separation-equivalenceclass within which the generating model lies. In
particular, we advocate using the PC algorithm(Spirtes et al. 1993) which is computationally efficient even for a large number of variables. Note
that, for linear models, to obtain the d-separationequivalence class it is sufficient to identify the zeropartial correlations in the data, as these depend
only on the linear coefficients and the variancesof the disturbances (and not on non-Gaussianity
aspects of the distributions). However, since thedata may well be signficantly non-Gaussian, nonparametric tests should optimally be used to findthe zero partial correlations.
2. For each DAG G in the estimated d-separationequivalence class:(a) Estimate the coefficients bij using ordinary
least-squares regression. (Note that this provides consistent estimates regardless of nonGaussianity of the variables.)(b) Calculate the corresponding residuals ei andrescale them to zero mean and unit variancefor each iii

Of course, since we only have samples wehave to take the sample mean rather thanthe expectation.
3. Select the highest-scoring DAG Gopt from Step
2 and apply a statistical test for normality foreach of the corresponding residuals ei . Using Definition 4, compute and return the ngDAG patternrepresenting the ngDAG (Gopt , ng) where ng is thevector indicating those residuals whose normality
was rejected by the normality tests.
The objective function U is commonly used in ICA asa measure of the non-Gaussianity of a random variable, and it can be shown to give a consistent estimator for finding independent components under weak
conditions (Hyvärinen et al. 2001). ICA estimation
is closely related to choosing the right DAG because
statistical independence of the estimated residuals is anecessary condition for the correct model: Any DAG
for which the estimated residuals are not independent
violates the assumptions of the model (see Section 2)
and hence cannot be the data-generating DAG. On
the other hand, any DAG which results in statisticallyindependent residuals represents one valid model that
could have generated the data.
Note that if we could disregard sampling effects,
distribution-equivalent models would attain exactly
the same value of U . However, in the practical case ofa finite sample this is not the case, thus Step 3 in the
PClingam algorithm is required to identify the correct
distribution-equivalence class.
The method as presented above has at least a couple ofshortcomings. One is that, for any given function f (ei )used, there always exist distributions which are nonGaussian yet are not distinguished from the Gaussian
by this measure. This is a well-known issue in ICA
which fortunately tends to have little practical significance since few such distributions are encountered inpractice. If needed, non-parametric Gaussianity measures could be used to remedy this potential problem.

a

b

2
5

-1.3

3

1.1

2

0.8

2

5

-1.4
0.5 6

c

d

2

5

3

6

e

2

5

3

6

5

3

6

3

6

0.6

4

1

4

1

4

1

4

1

4

1

Figure 2: One of the networks used in the simulations. Variables with non-Gaussian disturbances are shown insquares, while those with Gaussian disturbances are plotted as circles. (a) True data-generating model. (b)True d-separation-equivalence pattern. (c) True distribution-equivalence pattern. (d) Estimated DAG Gopt .
(e) Estimated distribution-equivalence pattern. See main text for details.
Naturally, in some cases many of the disturbances may
be slightly non-Gaussian yet sufficiently close to Gaussian that the available samples may not be sufficient
to distinguish the two and utilize the information for
determining causal directions in the model. Of course,
this is not a shortcoming of this particular method but
is a more general phenomenon.
Another important limitation is that the ICA objective function given above will only provide a propercomparison of different DAGs for which the residuals
ei are linearly uncorrelated. This is guaranteed to be
the case when the search is in the correct d-separationequivalence class, but if in Step 1 of the procedure
we select a too simple model (i.e. containing too few
edges) then the estimated disturbances may be linearly correlated and the objective function misleading.
Thus, it might be wise to include a term penalizing linear correlations such as is used in maximum likelihood
estimation of ICA (Hyvärinen et al. 2001). However,to keep our method as simple as possible, we haveomitted such a penalty term in this paper.

6

SIMULATIONSIn this section we report on simulations used to test
the performance of the PClingam method. First,
we tested the ability of the non-Gaussianity objective function (3) of Step 2 and the normality testsof Step 3 of PClingam to identify the correct ngDAG
pattern (distribution-equivalence class) when the true
d-separation-equivalence pattern was known. In otherwords, we tested how well the algorithm would function if Step 1 of the method worked flawlessly. Subsequently, we experimented with the full method incorporating the necessary estimation of the d-separationequivalence class (Step 1).
Figure 2a displays one of the models used to test theprocedure. The disturbance distributions of variables
X1 and X5 were a standard Gaussian the values ofwhich were squared (but keeping the original sign)

while the disturbance of X4 was produced in a similar
way but instead raising the values to the third power.
The disturbances of X2 , X3 , and X6 were Gaussian.
The disturbance variables were scaled such that theirvariances ranged from 1.0 to 3.0. A sample of 1000
data vectors was generated from the model.
Figure 2b shows the true d-separation-equivalence pattern of the model in (a). The equivalence class consistsof 12 different DAGs. However, the non-Gaussianityof the disturbances of X1 , X4 , and X5 means that
there are actually only 2 DAGs which are distributionequivalent; these are represented by the distributionequivalence pattern of Figure 2c. Figure 2d shows the
DAG Gopt found by Step 2 of PClingam from the data,when the true d-separation-equivalence class was given
to the algorithm. An Anderson-Darling test for normality (Anderson and Darling 1954) gave the p-values
0.000, 0.3145, 0.2181, 0.000, 0.000, and 0.0197 for thecorresponding residuals e1 to e6 . Inferring a residual to be non-Gaussian when p < 0.01 in Step 3 ofthe method produced the ngDAG pattern of Figure 2e,which turns out identical to the true ngDAG pattern in(c).
This basic procedure was repeated 20 times, with theresults summarized in Table 1a. In each simulation,
we randomly generated a linear acyclic causal modelover 6 variables, with each variable randomly chosen to have either a Gaussian or a non-Gaussian disturbance. The non-Gaussian distributions used were
those mentioned above as well as a Student’s t (2 degrees of freedom), a bimodal Mixture of Gaussians(0.5 N (−2, 1) + 0.5 N (2, 1)), a log-normal distribution(exponentiated standard normal) and a uniform distribution. The true d-separation-equivalence pattern was
input to the algorithm, to test the functioning of thePClingam method when the correct pattern is selectedin Step 1. The panel shows how often a specific type oftrue edge (in the true distribution-equivalence pattern)
gave rise to a specific type of estimated edge (in the estimated distribution-equivalence pattern). Rows cor-

Table 1: Summary of the simulations employing various methods for inferring the d-separation-equivalenceclass in Step 1 of PClingam. Each table is a confusion
matrix of arcs in the true distribution-equivalence patterns vs arcs in the estimated distribution-equivalencepattern. See main text for details.

aUsing the true
d-sep-equiv pattern

b

PC

185

0

0

0

183

0

0

2

0

12

2

0

0

12

2

0

0

0

61

0

9

1

43

8

0

0

0

40

3

2

6

29

c

d

CPC

GES

183

0

1

1

173

0

8

4

0

12

2

0

0

10

2

2

9

0

50

2

2

0

51

8

3

1

2

34

1

0

1

38

7

While the theoretical aspects of identifiability aresolved, at least a couple of important issues regarding the estimation of the model from finite samples
remain.
First and foremost, non-parametric methods for identifying zero partial correlations in non-Gaussian settings should be used so as to obtain better estimates of
the appropriate d-separation-equivalence class within
which to search. Although the methods developed for
Gaussian variables seem to work relatively well in ourpartly non-Gaussian setting, it is likely they will beoutperformed by methods that take into account thepossibility of non-Gaussian distributions.
Another important question is how to make the procedure scalable to data involving many (tens or even
hundreds of) variables. Although the current approachrelies on a brute-force enumeration of all DAGs in thed-separation-equivalence class, it would not be difficult to adapt the method to do a local search among
DAGs in an equivalence class. The extent to whichsuch a method would be hampered by local maxima is
unknown.

8

respond to the true edges, columns to estimated ones.
Optimally all off-diagonal elements would be zero. It
can be seen that the results are close to perfect; themethod misclassifies two undirected edges as directed,
but correctly estimates all others.
These simulations confirm that the PClingam method
works well at least when the d-separation-equivalenceclass can reliably be estimated. But in practice, withfinite datasets, there may be significant errors in inferring the d-separation-equivalence class. The degree towhich this affects the algorithm is an important practical issue.
Thus, in further simulations, we applied several different methods for learning d-separation-equivalencepatterns from the simulated data, as Step 1 in the
PClingam method. The methods we compared were
the PC algorithm (Spirtes et al. 1993), the Conservative PC algorithm (Ramsey et al. 2006), and the GES
algorithm (Chickering 2002). Panels b-d of Table 1
summarize the results. Although all of the methods
assumed Gaussianity when learning the d-separationequivalence pattern, the results are still quite encouraging, and a clear majority of edges were correctly
estimated.

FUTURE WORKSUMMARY

The discovery of linear acyclic causal models is a topic
which has been thoroughly investigated in the lasttwo decades. Both the Gaussian and the fully nonGaussian special cases are well understood, but thegeneral mixed case has not been previously discussed.
In this paper we have provided a complete characterization of distribution-equivalence and a practical estimation method in this setting.
Acknowledgements
The authors wish to thank Clark Glymour for helpful
and stimulating discussions. P.O.H. was funded by a
postdoctoral researcher grant from the University of
Helsinki.
References
Anderson, T. W. and D. A. Darling (1954). A testof goodness of fit. Journal of the American Statistical Association 49 (268), 765–769.
Andersson, S. A., D. Madigan, and M. D. Perlman(1997). A characterization of markov equivalenceclasses for acyclic digraphs. Annals of Statistics 25, 505–541.
Bollen, K. A. (1989). Structural Equations with Latent Variables. John Wiley & Sons.

Chickering, D. M. (2002). Optimal structure identification with greedy search. Journal of Machine
Learning Research 3, 507–554.
Geiger, D. and D. Heckerman (1994). Learning gaussian networks. In Proc. of the 10th Annual Conference on Uncertainty in Artificial Intelligence,
pp. 235–243.
Hyvärinen, A., J. Karhunen, and E. Oja (2001). Independent Component Analysis. Wiley.
Meek, C. (1995). Causal inference and causal explanation with background knowledge. In Proc. of
the 11th Annual Conference on Uncertainty in
Artificial Intelligence, pp. 403–411.
Pearl, J. (2000). Causality: Models, Reasoning, andInference. Cambridge University Press.
Ramsey, J., J. Zhang, and P. Spirtes (2006).
Adjacency-faithfulness and conservative causal
inference. In Proc. of the 22nd Annual Conference on Uncertainty in Artificial Intelligence.
Shimizu, S., P. O. Hoyer, A. Hyvärinen, and A. J.
Kerminen (2006). A linear non-gaussian acyclic
model for causal discovery. Journal of Machine
Learning Research 7, 2003–2030.
Spirtes, P., C. Glymour, and R. Scheines (1993).
Causation, Prediction, and Search. SpringerVerlag. (2nd ed. MIT Press 2000).

APPENDIX
Proof of Lemma 1: Select any ngDAG D represented
by the ngDAG pattern P in question. Now, for eachvariable Xi in D with a non-Gaussian disturbance,
add two auxiliary variables both of which have no parents and both of which have Xi as their only child.
Now, consider what the d-separation-equivalence pattern looks like for this augmented graph. Any edge in
D into such an Xi may be oriented in the d-separationequivalence pattern on the basis of the resulting unshielded collider at Xi . Furthermore, any edge in D
out of such an Xi may similarly be oriented on the basis of the lack of an unshielded collider. Thus, all edges
either originating from, or terminating in, a variablewith a non-Gaussian disturbance will be oriented in
the d-separation-equivalence pattern of the augmented
graph. Additionally, any edges whose orientations canbe deduced as a result of knowing the newly oriented
edges are oriented as well. Hence, adding the auxiliary
variables has exactly the same effect as simply orienting any edges connected to a non-Gaussian variable.
Because d-separation-equivalence patterns are always chain graphs (Andersson et al. 1997), thed-separation-equivalence pattern for the augmented
graph is a chain graph. Since each of the auxiliary

variables is connected to the rest of the graph only bya single oriented edge, removing these variables cannot change the chain graph property of the graph.
Hence, the distribution-equivalence pattern obtainedby orienting any edges originating from, or terminating in, variables with non-Gaussian variables (and subsequently orienting any edges deduced from these) has
to be a chain graph. This completes the proof.
Proof of Theorem 1: We will first prove that if
two ngDAGs D1 and D2 are distribution-equivalentthen they must share the same ngDAG pattern.
Since distribution-equivalence implies d-separationequivalence it is clear they have to share the dseparation-equivalence pattern. Thus we need to showthat it is always possible to correctly orient any edges
directly connected to any variable with a non-Gaussian
disturbance.
Consider a linear acyclic causal model M1 which instantiates D1 . Let us for simplicity of notation assume
that the variables have been named X1 , . . . , Xn suchthat j < i whenever Xj is an ancestor of Xi . Then it ispossible to collect the linear coefficients bij which represent the direct causal effects into a lower-triangularmatrix B such that we have x = Bx + e, where x
collects the observed variables X1 , . . . , Xn and e denotes the disturbance variables. Note that we have forsimplicity of notation and without loss of generality
assumed the constants in (1) are equal to zero. Solving for x we obtain x = Ae where A = (I − B)−1 isthe reduced-form matrix (representing the total effects
between the variables), which here is lower-triangular
due to the above causal ordering of the variables. By
the assumptions of the model, the components of e
are mutually independent. This means that our generating model is an ICA model, although several (potentially even all) of the components may be Gaussian. Although it is impossible to completely estimate A when there are two or more Gaussian components, it is well known (Hyvärinen et al. 2001) thatthe basis vectors (columns of A) corresponding to the
non-Gaussian components are identifiable (except forthe standard indeterminancy of permutation and scaling), as is the covariance matrix of the (generally multidimensional) Gaussian component which groups all
Gaussian disturbances.
Thus each non-Gaussian disturbance in the generating model essentially gives us a vector ai containing
the total effects of its corresponding variable Xi on
the other variables. As no variable can affect its nondescendants, and as we have assumed faithfulness, theset of non-zero entries of each such vector exactly represents the union of the correct corresponding variable
and all its descendants. If we could somehow know towhich observed variable Xi each non-Gaussian basis

vector ai should be paired we could easily orient alledges connecting to that variable as follows: Any observed variable Xj which we know (by d-separation)
to be connected to Xi but which has a zero entry inai has to be a parent of Xi . Similarly, any observed
variable Xj which we know (by d-separation) to beconnected to Xi but which has a non-zero entry in ai
has to be a child of Xi .
But is it always possible to identify to which observed
variable a given basis vector should be connected? We
now show that this is indeed so. Consider what happens if we try to build a model in which a basis vector
ai is paired with a descendant Xj of Xi , rather thanwith the correct choice Xi . Because we are restricted
to acyclic models, the new model may be represented
by a reduced-form matrix A0 which is lower-triangular
for the new variable order (but not for the old order,
because in this new model Xi is necessarily a descendant of Xj ). Note also that, when the rows are identically ordered, the column ai must be equal to column
a0j of A0 . Now, in the new model, Xi has to be represented by some disturbance variable, so there must be
a column a0i of A0 which has a non-zero entry in therow corresponding to Xi . Furthermore, a0i must havezeros wherever ai has zeros. These properties imply
that a0i cannot be expressed as a linear combinationof the {ak }k6=i : Since ai cannot be used, one would
have to depend on ak corresponding to Xk which areancestors of Xi to properly represent the non-zero entry corresponding to Xi in a0i . But among these there
necessarily has to be an ak which corresponds to a(relative) source Xk whose effect cannot be cancelledout by the others to lead to a required zero in a0i .
Thus a0i is not in span({ak }k6=i ). This implies thatspan({ak }k6=i ) 6= span({a0k }k6=j ) and hence the covariance implied by A0 cannot equal that implied by A
and the two models cannot be distribution-equivalent.
Hence to represent the correct distribution one has topair basis vectors corresponding to non-Gaussian components to the correct observed variable, allowing theidentification of the orientation of all edges directly
connected to that observed variable.
Next we need to prove that two ngDAGs which share thesame ngDAG pattern are distribution-equivalent. Callthe two ngDAGs D1 and D2 . Since they are represented
by the same ngDAG pattern they may only differ interms of the orientations of arrows within each chain
component of the chain graph that is the ngDAG pattern. (A chain component is a maximal connected setof nodes such that there is no oriented edge between
any pair of nodes.) Let Vc stand for the set of variables(necessarily all with Gaussian disturbances) making
up any given chain component. First, note that we
can focus our attention on only the variables Vc andany immediate parents of any of these: If we can matchany conditional distribution P (Vc | pa(Vc )) of a modelM1 represented by D1 using a model M2 represented
by D2 then we can always match the full joint distribution over all variables. Next consider the parents ofthe chain component. It is clear that each such parent must be connected to all the variables of the chain
component, otherwise some edges in the component
could be oriented using d-separation-equivalence.
Select any parametrization M1 of D1 , and call the
resulting reduced-form matrix between the variables
Vc only A1 . Because of the old result that dseparation-equivalence among Gaussian variables implies distribution-equivalence we know that there is aparametrization M2 of D2 such that the corresponding reduced-form matrix is A2 , and we have A1 AT1 =A2 AT2 . Now consider the effect of a parent. Denoting the weigths of the parent onto the variables Vc , inmodel M1 , by the vector w1 , we have Vc = A1 (e+w1 p)
where e represents the independent Gaussian disturbances of the Vc and p the value of the parent. Thus
the conditional distribution of the Vc given the value pof the parent is a Gaussian with mean vector A1 w1 p
and covariance matrix A1 AT1 . For model M2 we similarly have a Gaussian with mean A2 w2 p and covariance matrix A2 AT2 . As noted above, because of dseparation-equivalence, these covariance matrices canbe made identical with a suitable choice of A2 . Then,
selecting w2 = A−1
2 A1 w1 yields identical means as
well. (Note that the reduced-form matrices are always
invertible.). Also note that additional parents always
add independently, and the weights can be selected in
the above fashion for each parent separately. Thus, insummary, for any chain component on which models
D1 and D2 have conflicting edge orientations, for any
parametrization M1 of D1 it is always possible to finda parametrization M2 of D2 such that the conditional
distribution of the chain component given its parents
is identical. Thus the two models represent the exact
same set of joint distributions over all the variables.
This concludes the proof of the theorem.