<figure><img src="image_1.png" /></figure>

HAL
open science

# Multi-View Causal Discovery without Non-Gaussianity: Identifiability and Algorithms

Ambroise Heurtebise, Omar Chehab, Pierre Ablin, Alexandre Gramfort, Aapo Hyvärinen

► To cite this version:

Ambroise Heurtebise, Omar Chehab, Pierre Ablin, Alexandre Gramfort, Aapo Hyvärinen. Multi-View Causal Discovery without Non-Gaussianity: Identifiability and Algorithms. 2025. hal-05285946

HAL Id: hal-05285946

https://hal.science/hal-05285946v1

Preprint submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multi-View Causal Discovery without Non-Gaussianity: Identifiability and Algorithms

Ambroise Heurtebise¹*; Omar Chehab², Pierre Ablin³, Alexandre Gramfort¹, Aapo Hyvärinen⁴

¹Inria, CEA, University of Paris-Saclay, France

²University of Carnegie Mellon, Pittsburgh PA, USA

³Apple

⁴University of Helsinki, Finland

## Abstract

Causal discovery is a difficult problem that typically relies on strong assumptions on the data-generating model, such as non-Gaussianity. In practice, many modern applications provide multiple related views of the same system, which has rarely been considered for causal discovery. Here, we leverage this multi-view structure to achieve causal discovery with weak assumptions. We propose a multi-view linear Structural Equation Model (SEM) that extends the well-known framework of non-Gaussian disturbances by alternatively leveraging correlation over views. We prove the identifiability of the model for acyclic SEMs. Subsequently, we propose several multi-view causal discovery algorithms, inspired by single-view algorithms (DirectLiNGAM, PairwiseLiNGAM, and ICA-LiNGAM). The new methods are validated through simulations and applications on neuroimaging data, where they enable the estimation of causal graphs between brain regions.

## 1 Introduction

Causal discovery is a fundamental problem in scientific data analysis as well as several technological applications [22]. The basic problem is that we are given a number of observed variables, and we need to infer which variables cause which, and with what “connection strengths”. In the typical case, we only observe the data passively without any possibility of performing interventions, and this purely observational quality of the data makes the problem difficult.

Often, the problem is formalized as a structural equation model (SEM) [3], also called a functional causal model [20]. Then the question is how to estimate the SEM parameters using statistical theory. In fact, before even considering estimation methods, we need to know if it is at all possible to solve the problem. This is the question of *identifiability* of the model: can the parameters that describe the causal relations and directions be (uniquely) estimated?

A well-known fact is that if all variables are Gaussian and the model is linear, identifiability of the SEM is very problematic. Some of the earliest work showed that the directions can sometimes be recovered by an analysis of the *conditional independencies* between variables [33]; however, this is only possible in some special cases and notably not possible when observing only two variables. Recent literature has, therefore, focused on different departures from the linear-Gaussian framework to achieve identifiability. A major advance was to consider a linear model together with *non-Gaussianity* [30], which leads to full identifiability of the model under weak conditions such as acyclicity; however, such non-Gaussianity may not be strong enough in many data sets. A related framework was proposed by assuming that the cause undergoes a

<footnote>*Corresponding author ambroise.heurtebise@inria.fr</footnote>
nonlinear transform, while the disturbance is still additive [13]; however the nonlinearity may not be strong enough in many applications, and such nonlinearity slightly contradicts the linear additivity of the disturbance. Further related frameworks were proposed by [21, 38, 18]. Each framework has lead to development of a number of algorithms, the discussion of which we defer to the Background Section below.

An alternative recent framework is given by the multi-view setting, where data is collected from different sources whose similarities can be leveraged for estimation and identifiability. However, current approaches using multi-view structure are scarce and limited; either they assume the views are statistically independent [28] or make strong assumptions on their dependencies [5]. Here, we consider the linear multi-view SEM, providing a general theory and algorithms that leverage that multi-view structure.

Our main contributions are:

1. We generalize some of the most well-known algorithms for single-view causal discovery to the multi-view setting. The generalizations of PairwiseLiNGAM and DirectLiNGAM result in completely new and fast causal discovery using second-order statistics only, while our adaptation of ICA-LiNGAM improves on Chen et al. [5].

2. We show that our algorithms are consistent under weak conditions: we can replace the non-Gaussianity assumption on the disturbances by conditions on correlations between the views and “diversity” of such second-order statistics (SOS). Still, if the disturbances are non-Gaussian, we generalize Shimizu [28], Chen et al. [5].

3. The proofs above directly lead to rigorous identifiability results on the model, including cases where all the disturbances are Gaussian, or, equivalently, identifiability results based on SOS only. Importantly, only weak conditions on SOS is necessary for identifiability, greatly generalizing previous work.

4. We make our algorithms usable by practitioners by providing them as open-source, and demonstrate their usefulness by benchmarking them on real neuroimaging data.

## 2 Background

**Causal ordering** In the following, we consider a random vector $x \in \mathbb{R}^p$ whose entries are nodes of a directed acyclic graph (DAG). The DAG is represented by an adjacency matrix $B \in \mathbb{R}^{p \times p}$, where non-zero values indicate the edges of the graph. The entries of $x$ can be reordered such that each “causes” the next. This is formalized by the adjacency matrix admitting a special decomposition $B = P^\top T P$, where $T$ is a strictly lower triangular matrix and $P$ is a permutation matrix referred to as the causal ordering. In general, the causal ordering is not unique.

**LiNGAM: A single-view model for causal discovery** We consider a structural equation model (SEM), also known as a functional causal model [3, 20], that is linear and where the data follows

$$x = Bx + e \qquad (1)$$

and the entries $e_1, \dots, e_p$ of the random vector $e \in \mathbb{R}^p$ are independent noise terms, or disturbances. Causal discovery consists in inferring the parameters $B$ of the model, from observations of $x$. Yet, the identifiability of the model and therefore the uniqueness of $B$ is not straightforward. In fact, it is well-known that with Gaussian disturbances, the model is unidentifiable in general [26, 8]. A major advance by Shimizu et al. [30] was assuming non-Gaussian disturbances, leading to their Linear Non-Gaussian Acyclic Model (LiNGAM).
**Identifiability of LiNGAM** Shimizu et al. [30] showed that the LiNGAM model is identifiable in terms of the matrix $B$. A rigorous re-statement of this result is given in Appendix A.2. A further question that has received less attention is under what conditions the model is identifiable in terms of the causal ordering $P$. In fact, it is not in general: specifically, there may exist many permutation matrices $P$ and strictly lower triangular matrices $T$ such that the generated data has the same distribution and the generating permutation cannot be identified. For instance, in the degenerate case $B = T = 0$, any permutation matrix $P$ is equally valid and gives the same data distribution; any very sparse $B$ leads to several equally valid orderings. As is well-known, a DAG in general defines only a partial ordering in the sense that for some pairs of variables, we cannot necessarily say which is “earlier” and which is “later”. Thus, if we want to make the causal ordering unique — i.e. to obtain a total ordering — we need further assumptions, as considered in later sections.

**Estimation of LiNGAM** Since the LiNGAM model was proposed, two important categories of algorithms for estimating this model have been developed. First, algorithms based on recursively computing the residuals of pairwise regressions include DirectLiNGAM [32] and PairwiseLiNGAM [15]. These compare pairs of variables from $x$ and deduce what is the causal direction between the two variables based on regressing the variables in the two directions. By aggregating this information, they finally determine a causal ordering $P$. Once an ordering is found, the matrix $T$ can easily be obtained with conventional methods such as least-squares (LS) regression, which then gives $B = P^TTP$.

The second category is the ICA-based algorithm that was proposed in the original LiNGAM paper [31]. It is based on the observation that LiNGAM model can be rewritten as a latent variable model, in particular an ICA model [14] as

$$x = As. \qquad (2)$$

Again, the entries in $s$ are independent and non-Gaussian, and the “mixing” matrix $A$ expresses how the data is generated from the latent variables $s$. Many methods developed for ICA can be used to estimate the matrix $A$. However, it is important to consider the special structure of the mixing, since $A = (I - B)^{-1}$, where $B$ is a DAG [30]. This method first recovers the adjacency matrix $B$, and then estimates a causal ordering $P$ from $B$.

## 3 A Multi-View Model for Causal Discovery

We now extend the linear acyclic model in Eq. 1 to the multi-view setting, leading to a Linear Multi-View Acyclic Model (LiMVAM). Suppose we observe $m$ different views of a $p$-dimensional vector $x^i$. The most general linear acyclic extension is

$$x^i = B^i x^i + e^i \qquad (3)$$

for views $i \in [1, m]$. The $B^i$ are adjacency matrices of DAGs, and $e^i$ are view-specific disturbance vectors (with entries of non-zero variance). Similarly to Eq. 2, the LiMVAM model can equivalently be written as a latent-variable model

$$x^i = A^i e^i \qquad (4)$$

where the mixing matrix has a special structure, $A^i = (I - B^i)^{-1}$, which we will exploit later in Section 5. Our model generalizes that by Chen et al. [5], who assumed a specific structure on the influences, as will also be explained in Section 5.

To exploit the multi-view structure, we follow Shimizu [28], Chen et al. [5] and assume that all adjacency matrices $B^i$ share the same causal ordering. Equivalently, they admit a decomposition $B^i = P^T T^i P$ where $T^i$ is strictly lower triangular and $P$ is a permutation matrix that is
shared among the views i. Moreover, as commonly assumed in causal discovery, we require the components of each e^i to be mutually independent, for any given i. (In contrast, correlations between views i and i' will be seen to be essential for identifiability.) Contrary to most existing approaches, we do not impose any restriction on the distribution of the disturbances, such as non-Gaussianity, which greatly broadens the applicability of the model.

The next two sections develop efficient methods for estimating the adjacency matrices $B^i$, together with identifiability theory. We begin in Section 4 with methods that estimate the causal ordering using recursive residuals and using only second-order statistics, based on the single-view algorithms DirectLiNGAM [32] and PairwiseLiNGAM [15]. Then, we propose a modified version of the ICA-based multi-view model from Chen et al. [5] in Section 5.

# 4 Estimation by SOS and recursive residuals

First, we show how to estimate the model using Second-Order Statistics (SOS) only, and how the SOS alone lead to identifiability. To do so, we will rely on a principle that we call recursive residuals. While it has been used by Shimizu et al. [32], Hyvärinen and Smith [15], Shimizu [28], we provide here a unified treatment.

## 4.1 Causal ordering using recursive residuals

Methods based on what we call recursive residuals are built on the following principle: 1) we can find a root variable (i.e. one with no parents) by analyzing the causal directions between all the pairs of variables, and 2) when all other variables are regressed on a root variable, the resulting residuals preserve the causal ordering of the original variables [32, Lemma 2 and Corollary 1]. In the multi-view case, the same principles apply since all views are assumed to share a common ordering. A formal proof is provided in Appendix D.3.1 and D.3.2.

Algorithms based on recursive residuals therefore follow the same recursive scheme to recover the causal ordering: (i) perform pairwise regressions to obtain residuals, (ii) identify a root variable according to a specific criterion of causal direction computed on the residuals, and (iii) remove the selected root and repeat the procedure on the residuals until all variables are ordered. The main difference between algorithms lies in step (ii), i.e. in how the root variable is selected.

Since a root variable has no parents, we can find one by testing the causal direction for every pair of variables and taking as a root a variable that is never inferred to be an effect. A practical way to find the root variable from such tests with a finite sample is described in Appendix E.4.

Thus, root selection reduces to designing a criterion of causal direction capable of reliably distinguishing the causal direction between two variables. In the following, we introduce two new such criteria in the multi-view setting. For simplicity of notation, we consider two random variables aggregated across views: $x_1 = (x_1^1, \dots, x_1^m)^\top$ and $x_2 = (x_2^1, \dots, x_2^m)^\top$; either $x_1$ causes $x_2$, or the converse. Covariance matrices are denoted by $\Sigma_{x_1} = \mathbb{E}[x_1x_1^\top]$, and likewise for $x_2$, $e_1$ and $e_2$.

## 4.2 Likelihood-based criterion: PairwiseLiMVAM

We first derive a criterion that extends the PairwiseLiNGAM criterion of Hyvärinen and Smith [15] to the multi-view setting. Assuming Gaussian disturbances, we compute the expected log-likelihoods of two models: $x_1 \to x_2$ and $x_2 \to x_1$ (where an arrow denotes a directed edge between two adjacent variables in the causal graph). This leads to the following criterion, which is the logarithm of the maximum likelihood ratio (LR) between the two models:

$$ LR_{12} = - \log \det \Sigma_{x_1} - \log \det \Sigma_{e_2} + \log \det \Sigma_{x_2} + \log \det \Sigma_{e_1} . \quad (5) $$

We prove in Appendix B that, under very weak constraints on the correlations of variables and views, this criterion is positive when $x_1 \to x_2$, and negative in the opposite direction,
no matter whether the disturbances are Gaussian or non-Gaussian. Importantly, this makes the causal direction — and thus the entire causal ordering — identifiable under very general circumstances. So, although we developed this criterion for Gaussian disturbances, we emphasize that the resulting method is consistent regardless of the Gaussianity of the disturbances.

The necessary and sufficient conditions for the consistency of this criterion are given in Appendix B.6; a sufficient condition, assuming for simplicity that $x_1$ is the cause, is that two views $i, i'$ exist such that $|\text{corr}(e_2^i, e_2^{i'})| \neq |\text{corr}(x_1^i, x_1^{i'})|$, while $\text{corr}(x_1^i, e_2^i) \neq 0$ or $\text{corr}(x_1^{i'}, e_2^{i'}) \neq 0$. This essentially says that the views need to be correlated and the correlations need to be diverse, and will be explained in more detail below. In practice, we replace the covariance matrices $\Sigma_{x_1}$, $\Sigma_{e_1}$, $\Sigma_{x_2}$, and $\Sigma_{e_2}$ with consistent estimators (see Appendix E.1); we use Ordinary LS to compute the residuals throughout this paper. Thus, we obtain a consistent estimator of the true direction.

## 4.3 Cross-covariance-based criterion: DirectLiMVAM

Our next criterion is a new multi-view extension of the DirectLiNGAM criterion [32]. In a single-view setting, Shimizu et al. [32] relied on the fact that $x_1$ and $e_2$ are independent if the causal direction is $x_1 \to x_2$, and measured the independence using a kernel-based estimator of mutual information called KGV [2]. A multi-view extension was actually proposed by Shimizu [28], where the correct direction was chosen by comparing the sums over views of these scores for the two directions. However, such summation over views does not fully take advantage of the multi-view structure. Thus, the disturbances still needed to be non-Gaussian.

Here, we efficiently exploit the multi-view framework while only looking at SOS. We choose to evaluate the cross-covariance between the entire vectors $x_1 = (x_1^1, \dots, x_1^m)^\top$ and $e_2 = (e_2^1, \dots, e_2^m)^\top$ for the direction $x_1 \to x_2$, and between $x_2$ and $e_1$ for the direction $x_2 \to x_1$. Specifically, the correct direction is chosen by comparing the Frobenius norms of the cross-covariances between regressors and residuals. Thus, we obtain the following criterion:

$$ \mathrm{FC}_{12} = \left\| \mathbb{E}[x_2 e_1^\top] \right\|_F - \left\| \mathbb{E}[x_1 e_2^\top] \right\|_F . \qquad (6) $$

We prove in Appendix C that this criterion consistently finds the correct direction under the same conditions as the likelihood-based criterion, again no matter whether the disturbances are Gaussian or non-Gaussian.

## 4.4 Estimating the causal coefficients

Applying the recursive scheme explained in Section 4.1, and using the criteria from Eq. 5 or Eq. 6, we obtain a causal ordering encoded by $P$. Importantly, since the ordering is now known, the causal matrices $B^i$ can be uniquely recovered with little effort. In fact, the model becomes

$$ x^i = T^i x^i + e^i, \quad i \in [1, m] \qquad (7) $$

where $x^i = P x^i$ are the reordered observations and $e^i = P e^i$ are the reordered disturbances. By construction, each variable $x_j^i$ only depends on its predecessors $x_1^i, \dots, x_{j-1}^i$. For each $j \in [2, p]$, we estimate the $j$-th row of all matrices $T^i$ jointly using one-step Feasible Generalized Least Squares (FGLS) [37], which is more efficient than Ordinary LS. A detailed description of one-step FGLS is provided in Appendix E.5. Under the assumption that disturbances are mutually independent within each view, FGLS yields a unique regression for every variable on its parents. Finally, the adjacency matrices are recovered as $B^i = P^\top T^i P$.

The overall procedure is summarized in Algorithm 3 in the Appendix; it notably requires no hyperparameter tuning. Its worst-case computational complexity is in $O(m^3 \cdot p^3 \cdot n)$, which we show in Appendix E.7. In practice, these algorithms are highly parallelizable as they rely on basic algebraic operations (e.g. multiplying pairs of matrices), and fast as they do not require numerical optimization of an objective function.
## 4.5 Identifiability using second-order statistics

We can build an identifiability theory based on the conditions required by the criteria LR and FC (given in Appendices B and C). We start by formulating the following central assumption, which will be seen to be a sufficient condition for identifiability. Even milder sufficient conditions are discussed in Appendix D.1 but they are harder to interpret.

**Assumption 1** (Correlation and diversity across views) For any two variables $x_j$ and $x_{j'}$ such that $x_j \to x_{j'}$ (in the sense that $B_{j'j}^i \neq 0$ for at least one $i$), there exist two distinct views $i \neq i'$ such that a) those variables are correlated in at least one of those views, $\text{corr}(x_j^i, x_{j'}^i) \neq 0$ or $\text{corr}(x_j^{i'}, x_{j'}^{i'}) \neq 0$, and b) their correlations fulfill the “diversity” condition:

$$ |\text{corr}(x_j^i, x_{j'}^{i'})| \neq |\text{corr}(e_{j'j}^i, e_{j'j}^{i'})| . \qquad (8) $$

To better understand this assumption, note first that it requires some correlation across views: it rules out the trivial case where all cross-view correlations are zero. Moreover, it demands some diversity in these correlations. Consider the extreme situation where, for standardized variables, $x_j^i = \alpha x_{j'}^i$ and $x_{j'}^i = \alpha x_{j}^i$ (and analogously for the residuals), for some scalar $\alpha$. In this case, all correlations equal 1, so condition Eq. 8 is violated. Appendix D.2 provides more realistic examples of violations. We next introduce an additional assumption that ensures identifiability of the total causal ordering.

**Assumption 2** (Dense connectivity of the graph when pooled across views) Let $\mathcal{G}$ denote the union graph of the $m$ views, obtained by taking the union of their edge sets, and assume that the ordering induced by $\mathcal{G}$ is total, i.e. there exists a directed path between any two variables.

With these assumptions, we prove in Appendix D.3 the following identifiability result

**Theorem 1** (Identifiability of the ordering and adjacency matrices $\mathbf{B}^i$) Assume the data follows the model in Eq. 3, and that Assumption 1 holds. Then, causal (partial) ordering and causal coefficients $\mathbf{B}^i$ are identifiable. If we further make Assumption 2, then $\mathbf{P}$ is identifiable as well.

The identifiability of the partial ordering does not necessarily mean that the permutation matrix $\mathbf{P}$ can be recovered uniquely. If the ordering is only partial (for example, in a three-nodes graph: $x_1 \to x_2$ and $x_1 \to x_3$, but there is no relation between $x_2$ and $x_3$, so both orderings $1 \to 2 \to 3$ and $1 \to 3 \to 2$ are valid), then $\mathbf{P}$ is only identifiable up to a class of permutations that induce the same partial ordering. The above theorem shows a case where $\mathbf{P}$ can be (uniquely) identified.

# 5 Estimation by Non-Gaussianity assuming shared disturbances

Next, we provide an alternative approach to the LiMVAM model definition and estimation. We utilize the non-Gaussianity of the disturbances to obtain identifiability conditions that are different from the previous section. However, we still do not impose non-Gaussianity as a necessary condition.

**Model definition** The approach here is to propose a particular model for the disturbances by supposing that they can be additively composed into two terms: one term is view-specific, while the other is shared across views and therefore makes the views dependent. Specifically:

**Assumption 3** (Shared disturbances) The disturbances can be decomposed as

$$ e^i = \mathbf{D}^i \mathbf{s} + \mathbf{n}^i \qquad (9) $$

where the $\mathbf{D}^i$ are diagonal matrices with positive entries on the diagonal, the $\mathbf{s}$ has mutually independent entries and second-order moment $\mathbb{E}[\mathbf{s}\mathbf{s}^\top] = \mathbf{I}_p$, and the view-specific noises are
Gaussian $\mathbf{n}^i \sim \mathcal{N}(\mathbf{0}, \mathbf{\Sigma}^i)$ and depend on the view $i$ via the second-order moments $\mathbf{\Sigma}^i$ that are diagonal matrices. Lastly, the vectors $\mathbf{s}$ and $\mathbf{n}^i, i = 1, \dots, m$, are assumed to be mutually independent.

This Assumption and Eq. 9 are related to Chen et al. [5] who where originally inspired by the multi-view Shared ICA of Richard et al. [25]. Specifically, Chen et al. [5] considered the special case where the matrices $\mathbf{D}^i = \mathbf{I}_p$, so that there are the same disturbances $\mathbf{s}$ over views. We argue here that their model was too restrictive and unrealistic as it imposes an arbitrary scaling for the data variables. Consider the root variable of the DAG, denoting it by $x_i$. Its disturbance is $s_i + n_i$ which has variance of at least one since the variance of $s_i$ is defined as one. Thus, the root variable has variance of at least one which is absurd: the model should be invariant to the scaling of the variables. Likewise, suppose we create a new data set by dividing just one view, say $\mathbf{x}^1$ by, say 10, giving new data $\mathbf{x}'^1$. This new data does not follow the model by Chen et al. [5] anymore since $\mathbf{s}$ should be rescaled to $\mathbf{s}/10$ for just this view; now the $\mathbf{s}/10$ here is different from the $\mathbf{s}$ in other views. Such problems motivated us to develop a model where the scaling terms $\mathbf{D}^i$ are new parameters, and thus to the definition in Eq. 9.

**Identifiability** Now we proceed to show the identifiability of the model based on Assumption 3. The starting point is that our model can be rewritten as a multi-view ICA model as in Eq. 4, just like in the case of LiNGAM. Thus, the methods developed for multi-view ICA can be used to analyze its identifiability. We first discuss the assumptions used for our theorems.

**Assumption 4** (Diversity of the views with Gaussian disturbances) If $j$ and $j', j \neq j'$, are the indices of two Gaussian common disturbances in $\mathbf{s}$, then the number of views is at least three, and there exists a view $i$ such that $\frac{\mathbf{\Sigma}_{jj}^i}{(\mathbf{D}_{jj}^i)^2} \neq \frac{\mathbf{\Sigma}_{j'j}^i}{(\mathbf{D}_{j'j}^i)^2}$.

This assumption, inspired by Richard et al. [25], is trivially fulfilled if all disturbances are non-Gaussian. For Gaussian disturbances, it essentially states that the views should be diverse in terms of disturbances' SOS. This follows from the theory of Shared ICA which does not require the common disturbances $\mathbf{s}$ to be non-Gaussian if there is enough diversity [25, 1]. The next assumption assumes that the union graph of the $\mathbf{B}^i$ is dense enough across views.

**Assumption 5** (Dense connectivity of the graph when pooled across views) There exists a permutation $\bar{\mathbf{P}}$ such that the matrices $\bar{\mathbf{T}}^i = \bar{\mathbf{P}}\bar{\mathbf{B}}^i\bar{\mathbf{P}}^\top$ are strictly lower triangular, and for each entry $(j, k)$ with $j > k$, there is at least one view $i$ such that $\bar{\mathbf{T}}_{j,k}^i \neq 0$.

Note that this assumption is quite weak as the $\mathbf{T}^i$ can be very sparse. This is especially the case when the number of views is large, since only one non-zero entry is required over all views. We can now state our main identifiability result for this variant of the model.

**Theorem 2 (Identifiability with shared disturbances)** Assume the LiMVAM model in Eq. (3) together with Assumption 3. Under Assumption 4, the model is identifiable in the sense that the parameters $\mathbf{B}^i$, $\mathbf{D}^i$, and $\mathbf{\Sigma}^i$ are identifiable. If we further make Assumption 5, $\mathbf{P}$ is identifiable as well.

The proof of the theorem is given in Appendix F. We thus see that for non-Gaussian disturbances, identifiability is achieved with weak conditions; in fact, no special conditions on the covariances are needed. Still, the theory of Shared ICA leads to identifiability even for Gaussian disturbances, but with stronger conditions than obtained in Section 4. The utility of the Shared ICA theory lies also in the fact that it allows us to construct an algorithm very easily.
<figure><img src="image_3.png" /><figcaption>Figure 1: Separation performance of three recursive residuals algorithms and three ICA-based algorithms. We assume shared disturbances as in Eq. 9. We used 50 different seeds. Lower is better.</figcaption></figure>

**Estimation** The model can be estimated by a combination of the ICA-based estimation procedure by Shimizu et al. [30], combined with the Shared ICA methods by Richard et al. [25]. We call it ICA-LiMVAM and explain it in detail in Algorithm 4. Surprisingly, the resulting algorithm gives estimates of $B_i$ that are identical to Chen et al. [5], even though our model is more general. In contrast, our algorithm estimates different noise variances and scalings $D_i$. The worst-case computational complexity is in $O(tnmp^3 + mp^5)$, for $t$ iterations of the Shared ICA (ML) algorithm, $n$ samples, $m$ views and $p$ components. We detail this in Appendix G.2.

# 6 Experiments

## 6.1 Simulations

We benchmark a total of six multi-view causal discovery algorithms on synthetic data. Three of these algorithms are based on recursive residuals: our PairwiseLiMVAM and DirectLiMVAM, as well as Multi-GroupDirectLiNGAM from Shimizu [28]. The other three algorithms are ICA-based: ICA-LiMVAM which is essentially the same as Chen et al. [5], as well as ICA-LiNGAM [30] which is naively applied to each view, separately. Note that ICA-LiMVAM comes in two variants, each using a different Shared ICA algorithm: either Shared ICA “ML”, a likelihood-based method that can handle both Gaussian and non-Gaussian disturbances, or Shared ICA “J”, which jointly diagonalizes covariance matrices without using non-Gaussianity.

The first synthetic experiment is inspired by Richard et al. [25]. We monitor the performance of causal discovery algorithms across varying sample sizes and disturbance distributions. More specifically, the data are generated according to the LiMVAM model with shared disturbances

<figure><img src="image_4.png" /><figcaption>Figure 2: Estimation error of the adjacency matrices $B^i$ versus total fitting time for each method, across 200 random runs with $m = 6$ views, $p = 5$ variables, and $n = 1000$ samples. Disturbances are drawn from a standard Gaussian. Note that the point clouds for PairwiseLiMVAM and DirectLiMVAM largely overlap.</figcaption></figure>
from Eq. 9. The performance of the causal discovery algorithms is measured by the $l_2$ distance between true and estimated adjacency matrices $B^i$, and the percentage of incorrectly estimated causal orders $P$ over 50 random runs.

Results of the first synthetic experiment are plotted in Figure 1. Our algorithms Pairwise-LiMVAM and DirectLiMVAM consistently outperform all others by a significant margin. ICA-LiMVAM yields reliable estimates but with larger errors. In contrast, MultiGroupDirectLiNGAM and ICA-LiNGAM, which are designed for non-Gaussian disturbances, fail entirely in the Gaussian setting, as expected. Note that the high error rate of ICA-LiNGAM in recovering the causal ordering $P$ is likely explained by its inability to exploit that the ordering is shared across views. For the second synthetic experiment, we report the trade-off between computational complexity and estimation error of the adjacency matrices. We use a different data generation model, described by Eq. 3 where the disturbances are Gaussian with equal variances. This breaks the assumptions for all methods but our PairwiseLiMVAM and DirectLiMVAM. As expected, we see in Figure 2 that our algorithms PairwiseLiMVAM and DirectLiMVAM provide a major reduction in estimation errors, at little or no computational cost. In contrast, the four other algorithms struggle with this data.

In Appendix H, we include further experiments that test Assumptions 4 and 5, as well as scalability with more views and variables. In particular, we find that DirectLiMVAM outperforms PairwiseLiMVAM in some cases, which is less obvious from Figures 1 and 2 alone.

## 6.2 Real data experiments with MEG

Next, we performed experiments on magnetoencephalography (MEG) data measuring human brain activity. MEG data has a high temporal resolution, so typically methods related to Granger causality are used. However, we are here interested in analyzing the *energies* of *oscillatory* signals. The energies change very slowly (on the time-scale of seconds), while the underlying brain activity being measured is likely to change much more quickly. Thus, it is appropriate to use a model of instantaneous causality as in our model.

We used the Cam-CAN dataset, which is the largest publicly available MEG dataset [27, 36], and considered the “sensorimotor task” during which each participant had to respond with a right index finger button press to auditory/visual stimuli. Further details are available in Appendix H.2.1. We applied PairwiseLiMVAM to these data, each of the 98 participants being one view. The experiment ran in 2 hours and 20 minutes using 5 CPUs. Although our model does not require or assume any similarity across the matrices $T^i$, one may reasonably hypothesize that the participants’ brains share structural patterns. So, to facilitate population-level interpretation, we computed the element-wise median of the $B^i$ matrices.

<figure><img src="image_5.png" /><figcaption>Figure 3: (a) Top ten strongest median causal effects (estimated by PairwiseLiMVAM) across 98 Cam-CAN (MEG) participants. (b) Pearson correlations between median causal matrices across 50 runs in MEG. (c) Top ten strongest median causal effects from PairwiseLiMVAM across 9 fMRI participants. Arrows indicate direction and strength of effects, with width proportional to magnitude.</figcaption></figure>
Fig. 3a displays the ten strongest (in absolute value) median causal effects. Notably, many causal connections are between homologous regions across hemispheres — particularly within the primary motor and somatosensory cortices — consistent with prior findings that symmetric brain regions exhibit strong inter-hemispheric correlations [17]. Additionally, numerous arrows involve the postcentral regions, which have been implicated in fine motor control and hand-related tasks [4]. To evaluate the robustness of the results, we repeated the experiment 50 times, each time randomly selecting 30 participants from the full cohort. We assessed consistency by measuring the Pearson correlation between the resulting median matrices across runs. As shown in Figure 3b, the median effects are highly correlated across different subsets (average correlation: 0.67), demonstrating the stability of our method. Additional experiments using ICA-LiMVAM-ML are reported in Appendix H.2.2.

## 6.3 Real data experiments with fMRI

We also evaluated PairwiseLiMVAM on an fMRI rhyming judgment task dataset [24]. The data consist of recordings from nine participants, each represented by nine variables: one task regressor (Input I), obtained by convolving the boxcar design of the rhyming task with a canonical hemodynamic response function, and eight regional time series from bilateral regions of interest. The dataset was acquired on a 3T scanner with a repetition time of 2 s, yielding 160 samples per subject, and is available via the OpenfMRI project (preprocessed version from <fnref n="1" />).

We report the ten strongest effects from the element-wise median of the estimated adjacency matrices in Figure 3c. The recovered structure is consistent with established findings for this task. In particular, we find a strong edge from the task regressor to the left occipital cortex (Input → LOCC), reflecting the expected visual drive of the stimuli, and prominent connections from left to right homologous regions (e.g. LOCC → ROCC, LIFG → RIFG, LIPL → RIPL), in line with the left-hemisphere dominance and inter-hemispheric influences reported by Ramsey et al. [24]. Edges such as LOCC → LACC and connections between LACC and RACC also match relationships highlighted in their analyses.

# 7 Discussion and Related work

A multi-view causal model based on LiNGAM was considered by Shimizu [28], but their model was very different. The fundamental difference is that we assume the disturbances have some shared information across views and can be Gaussian, whereas Shimizu [28] estimates disturbances that are independent across views and necessarily non-Gaussian. Thus, they only showed how to use the multiple views to improve the estimation while still assuming non-Gaussianity of the disturbances.

Chen et al. [5] considered multi-view SEM using a special case of our shared disturbances model, which is in its turn a special case of our general LiMVAM model. Algorithmically, there is little difference between their work and our ICA-LiMVAM. However, their identifiability proofs are not available nor did they consider the identifiability of the causal ordering, while we do both. Moreover, as we already argued in Section 5, their model formulation has the serious flaw that it forces an arbitrary scaling of the data variables. Crucially, we propose new algorithms, PairwiseLiMVAM and DirectLiMVAM, that use only SOS; they are shown to be empirically superior to the ICA-LiMVAM in most cases.

A related very general framework, allowing for various kinds of interventions, was proposed by Mooij et al. [19]. On the other hand, a related “multi-context” approach, where the disturbances are not necessarily related at all while the DAGs are, was proposed by Sturma et al. [34].

<footnote n="1">¹https://github.com/cabal-cmu/Feedback-Discovery</footnote>
# 8 Conclusion

We propose a new approach for causal discovery in the multi-view linear-DAG case, improving the earlier work by Shimizu [28], Chen et al. [5]. We show that the model is fully identifiable by second-order statistics (covariance structure) only, and under much weaker conditions than by Chen et al. [5]. We develop algorithms that either use the SOS only or combine them with non-Gaussian statistics. In particular, we adapt the Pairwise and Direct approaches to the multi-view case, resulting in highly efficient estimation. Results on brain imaging data are promising.

# References

[1] Anderson, M., Fu, G., Phlypo, R., and Adali, T. (2013). Independent vector analysis: Identification conditions and performance bounds. *IEEE Transactions on Signal Processing*, 62:4399-4410.

[2] Bach, F. R. and Jordan, M. I. (2002). Kernel independent component analysis. *Journal of machine learning research*, 3(Jul):1-48.

[3] Bollen, K. A. (1989). *Structural equations with latent variables*, volume 210. John Wiley & Sons.

[4] Braun, C., Haug, M., Wiech, K., Birbaumer, N., Elbert, T., and Roberts, L. E. (2002). Functional organization of primary somatosensory cortex depends on the focus of attention. *Neuroimage*, 17(3):1451-1458.

[5] Chen, W., Huang, X., Li, Z., Cai, R., Huang, Z., and Hao, Z. (2024). Individual causal structure learning from population data. In *Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence*, pages 7109-7117.

[6] Comon, P. (1994). Independent component analysis, a new concept? *Signal Processing*, 36(3):287-314. Higher Order Statistics.

[7] Dale, A. M., Fischl, B., and Sereno, M. I. (1999). Cortical surface-based analysis: I. segmentation and surface reconstruction. *Neuroimage*, 9(2):179-194.

[8] Genin, K. (2021). Statistical undecidability in linear, non-gaussian causal models in the presence of latent confounders. In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W., editors, *Advances in Neural Information Processing Systems*, volume 34, pages 13564-13574. Curran Associates, Inc.

[9] Gramfort, A., Luessi, M., Larson, E., Engemann, D. A., Strohmeier, D., Brodbeck, C., Goj, R., Jas, M., Brooks, T., Parkkonen, L., and Hämäläinen, M. (2013). Meg and eeg data analysis with mne-python. *Frontiers in Neuroscience*, 7.

[10] Gramfort, A., Luessi, M., Larson, E., Engemann, D. A., Strohmeier, D., Brodbeck, C., Parkkonen, L., and Hämäläinen, M. S. (2014). Mne software for processing meg and eeg data. *NeuroImage*, 86:446-460.

[11] Greene, W. H. (2003). *Econometric analysis*. Pearson education india.

[12] Hämäläinen, M. S. and Ilmoniemi, R. J. (1994). Interpreting magnetic fields of the brain: minimum norm estimates. *Medical & Biological Engineering & Computing*, 32(1):35-42.

[13] Hoyer, P., Janzing, D., Mooij, J. M., Peters, J., and Schölkopf, B. (2008). Nonlinear causal discovery with additive noise models. *Advances in neural information processing systems*, 21.
[14] Hyvärinen, A., Karhunen, J., and Oja, E. (2001). *Independent Component Analysis*. Wiley Interscience.

[15] Hyvärinen, A. and Smith, S. M. (2013). Pairwise likelihood ratios for estimation of non-Gaussian structural equation models. *Journal of Machine Learning Research*, 14:111–152.

[16] Khan, S., Hashmi, J. A., Mamashli, F., Michmizos, K., Kitzbichler, M. G., Bharadwaj, H., Bekhti, Y., Ganesan, S., Garel, K.-L. A., Whitfield-Gabrieli, S., Gollub, R. L., Kong, J., Vaina, L. M., Rana, K. D., Stufflebeam, S. M., Hämäläinen, M. S., and Kenet, T. (2018). Maturation trajectories of cortical resting-state networks depend on the mediating frequency band. *NeuroImage*, 174:57–68.

[17] Li, A., Yetkin, F. Z., Cox, R., and Haughton, V. M. (1996). Ipsilateral hemisphere activation during motor and sensory tasks. *American Journal of Neuroradiology*, 17(4):651–655.

[18] Monti, R. P., Zhang, K., and Hyvärinen, A. (2019). Causal discovery with general non-linear relationships using non-linear ICA. In *Proc. 35th Conf. on Uncertainty in Artificial Intelligence (UAI2019)*, Tel Aviv, Israel.

[19] Mooij, J. M., Magliacane, S., and Claassen, T. (2020). Joint causal inference from multiple contexts. *Journal of machine learning research*, 21(99):1–108.

[20] Pearl, J. (2009). *Causality: Models, Reasoning and Inference*. Cambridge University Press, USA, 2nd edition.

[21] Peters, J. and Bühlmann, P. (2014). Identifiability of gaussian structural equation models with equal error variances. *Biometrika*, 101(1):219–228.

[22] Peters, J., Janzing, D., and Schölkopf, B. (2017). *Elements of causal inference: foundations and learning algorithms*. The MIT Press.

[23] Power, L., Allain, C., Moreau, T., Gramfort, A., and Bardouille, T. (2023). Using convolutional dictionary learning to detect task-related neuromagnetic transients and ageing trends in a large open-access dataset. *NeuroImage*, 267:119809.

[24] Ramsey, J. D., Hanson, S. J., Hanson, C., Halchenko, Y. O., Poldrack, R. A., and Glymour, C. (2010). Six problems for causal inference from fmri. *neuroimage*, 49(2):1545–1558.

[25] Richard, H., Ablin, P., Thirion, B., Gramfort, A., and Hyvarinen, A. (2021). Shared independent component analysis for multi-subject neuroimaging. In *Advances in Neural Information Processing Systems*, volume 34, pages 29962–29971. Curran Associates, Inc.

[26] Richardson, T. and Spirtes, P. (2002). Ancestral graph Markov models. *The Annals of Statistics*, 30(4):962–1030.

[27] Shafto, M. A., Tyler, L. K., Dixon, M., Taylor, J. R., Rowe, J. B., Cusack, R., Calder, A. J., Marslen-Wilson, W. D., Duncan, J. S., Dalgleish, T., Henson, R. N. A., Brayne, C., and Matthews, F. E. (2014). The cambridge centre for ageing and neuroscience (cam-can) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing. *BMC Neurology*, 14.

[28] Shimizu, S. (2012). Joint estimation of linear non-gaussian acyclic models. *Neurocomputing*, 81:104–107.

[29] Shimizu, S. (2022). *Statistical Causal Discovery: LiNGAM Approach*. Springer.

[30] Shimizu, S., Hoyer, P. O., Hyvärinen, A., Kerminen, A., and Jordan, M. (2006). A linear non-gaussian acyclic model for causal discovery. *Journal of Machine Learning Research*, 7(10).
[31] Shimizu, S., Hyvärinen, A., Kano, Y., and Hoyer, P. O. (2005). Discovery of non-gaussian linear causal models using ica. In *Conference on Uncertainty in Artificial Intelligence*.

[32] Shimizu, S., Inazumi, T., Sogawa, Y., Hyvarinen, A., Kawahara, Y., Washio, T., Hoyer, P. O., Bollen, K., and Hoyer, P. (2011). Directlingam: A direct method for learning a linear non-gaussian structural equation model. *Journal of Machine Learning Research-JMLR*, 12(Apr):1225–1248.

[33] Spirtes, P., Glymour, C., and Scheines, R. (2001). *Causation, prediction, and search*. MIT press, second edition edition.

[34] Sturma, N., Squires, C., Drton, M., and Uhler, C. (2023). Unpaired multi-domain causal representation learning. In Oh, A., Naumann, T., Globerson, A., Saenko, K., Hardt, M., and Levine, S., editors, *Advances in Neural Information Processing Systems*, volume 36, pages 34465–34492. Curran Associates, Inc.

[35] Taulu, S. and Simola, J. (2006). Spatiotemporal signal space separation method for rejecting nearby interference in meg measurements. *Physics in Medicine & Biology*, 51(7):1759.

[36] Taylor, J. R., Williams, N., Cusack, R., Auer, T., Shafto, M. A., Dixon, M., Tyler, L. K., Cam-CAN, and Henson, R. N. (2017). The cambridge centre for ageing and neuroscience (cam-can) data repository: Structural and functional mri, meg, and cognitive data from a cross-sectional adult lifespan sample. *NeuroImage*, 144:262–269. Data Sharing Part II.

[37] Zellner, A. (1962). An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias. *Journal of the American Statistical Association*, 57(298):348–368.

[38] Zhang, K. and Hyvärinen, A. (2009). On the identifiability of the post-nonlinear causal model. In *Proc. 25th Conference on Uncertainty in Artificial Intelligence (UAI2009)*, pages 647–655, Montréal, Canada.
# Appendix

The appendix is organized as follows.

In Section A, we rigorously restate known identifiability results for LiNGAM.

In Section B, we prove the consistency of the likelihood-based criterion for a pair of variables.

In Section C, we prove the consistency of the cross-covariance-based criterion for a pair of variables.

In Section D, we extend the consistency of these criteria for a pair of variables to *many* pairs, and show that the LiMVAM model is identifiable.

In Section E, we show how to implement algorithms used to identify the causal graph.

In Section F, we prove the identifiability of a special case of LiMVAM with shared disturbances.

In Section G, we present an algorithm for estimating the causal graph in this special case.

In Section H, we detail our experiments on real and synthetic data.

<table><thead><tr><th></th><th></th></tr></thead><tbody><tr><td><strong>A LiNGAM</strong></td><td><strong>16</strong></td></tr><tr><td>A.1 Technical lemmas</td><td>16</td></tr><tr><td>A.2 Identifiability of LiNGAM</td><td>16</td></tr><tr><td>A.3 Proofs of technical lemmas</td><td>17</td></tr><tr><td><strong>B Likelihood-based criterion</strong></td><td><strong>19</strong></td></tr><tr><td>B.1 Basic setting</td><td>19</td></tr><tr><td>B.2 Derivation of the expected log-likelihood in one direction</td><td>19</td></tr><tr><td>B.3 Derivation of the likelihood-based criterion</td><td>19</td></tr><tr><td>B.4 Non-negativity of the likelihood-based criterion</td><td>20</td></tr><tr><td>B.5 Proof of a useful lemma</td><td>21</td></tr><tr><td>B.6 Conditions for strict positivity of the likelihood-based criterion</td><td>23</td></tr><tr><td><strong>C Cross-covariance-based criterion</strong></td><td><strong>25</strong></td></tr><tr><td>C.1 Derivation of the cross-covariance-based criterion</td><td>25</td></tr><tr><td>C.2 Conditions for strict positivity of the cross-covariance criterion</td><td>25</td></tr><tr><td><strong>D Identifiability of LiMVAM</strong></td><td><strong>26</strong></td></tr><tr><td>D.1 Multivariate conditions for strict positivity of a criterion</td><td>26</td></tr><tr><td>D.2 An example for interpreting the SOS-based identifiability condition in Assumption 1</td><td>26</td></tr><tr><td>D.3 Proof of Theorem 1</td><td>28</td></tr><tr><td><strong>E SOS-based Algorithms for LiMVAM</strong></td><td><strong>31</strong></td></tr><tr><td>E.1 Consistent estimator of the likelihood-based criterion</td><td>31</td></tr><tr><td>E.2 Consistent estimator of the cross-covariance-based criterion</td><td>31</td></tr><tr><td>E.3 Algorithms for estimating the two criteria</td><td>32</td></tr><tr><td>E.4 How to find the root variable</td><td>32</td></tr><tr><td>E.5 One-step Feasible Generalized Least Squares (FGLS)</td><td>33</td></tr><tr><td>E.6 Full algorithm</td><td>34</td></tr><tr><td>E.7 Computational complexity</td><td>35</td></tr><tr><td><strong>F Identifiability of LiMVAM with shared disturbances</strong></td><td><strong>36</strong></td></tr><tr><td>F.1 Proof of Theorem 2 (first claim)</td><td>36</td></tr><tr><td>F.2 Proof of Theorem 2 (second claim)</td><td>36</td></tr><tr><td>F.3 Results for view-specific causal orderings</td><td>37</td></tr><tr><td><strong>G ICA-based Algorithm for LiMVAM with shared disturbances</strong></td><td><strong>39</strong></td></tr><tr><td>G.1 ICA-based Algorithm</td><td>39</td></tr><tr><td>G.2 Computational complexity</td><td>40</td></tr></tbody></table>
<table><tr><td>H Experiments</td><td>41</td></tr><tr><td>H.1 Synthetic experiments</td><td>41</td></tr><tr><td>H.2 Real data experiments</td><td>43</td></tr></table>
# A LiNGAM

**Defining a domain connecting ICA and SEM parameters** In the following proofs, we define the set

$$ W = \{ \mathbf{W} \in \mathbb{R}^{p \times p} | \text{ there exist a diagonal matrix } \mathbf{D} \text{ with positive diagonal entries and a DAG matrix } \mathbf{B}, \text{ such that } \mathbf{W} = \mathbf{D}^{-1}(\mathbf{I} - \mathbf{B}) \} \quad (10) $$

as the domain of the “unmixing matrices” $\mathbf{W} = \mathbf{A}^{-1}$ in the ICA model from Eq. 2, that are compatible with the DAG structure in the structural equation model (SEM) from Eq. 1. In general, unmixing matrices are not constrained in the ICA model, except by invertibility (and possibly by normalizing their rows). However, in the context of SEM estimation, the requirement that $\mathbf{W} = \mathbf{D}^{-1}(\mathbf{I} - \mathbf{B})$ belongs to the set $W$ is a key structural constraint. Furthermore, it should be noted that finding $\mathbf{W}$ is equivalent to finding $\mathbf{D}$ and $\mathbf{B}$, since they are related by $\mathbf{B} = \mathbf{D}^{-1}(\mathbf{I} - \mathbf{W})$.

Finally, we will use the fact that any DAG matrix $\mathbf{B}$ can be decomposed into $\mathbf{B} = \mathbf{P}^{\top}\mathbf{T}\mathbf{P}$, where $\mathbf{P}$ is a permutation matrix and $\mathbf{T}$ a strictly lower triangular matrix.

## A.1 Technical lemmas

The following two lemmas are about matrices in the context of DAGs. They are used in the proofs of Theorems 5, 2, and 8, and are proven in Appendix A.3.

The basic identifiability results for the causal matrices $\mathbf{B}^i$ use the following:

**Lemma 3** Let $\mathbf{W} = \mathbf{D}^{-1}(\mathbf{I} - \mathbf{B})$ and $\mathbf{W}' = \mathbf{D}'^{-1}(\mathbf{I} - \mathbf{B}')$ be two matrices that belong to $W$ (defined in Eq. 10), and let $\mathbf{Q}$ be a sign-permutation matrix. Then

$$ \mathbf{W}' = \mathbf{Q}^{\top}\mathbf{W} \implies \mathbf{Q} = \mathbf{I}, \mathbf{D}' = \mathbf{D}, \text{ and } \mathbf{B}' = \mathbf{B}. \quad (11) $$

On the other hand, the identifiability results for the causal ordering $\mathbf{P}$ (or causal orderings $\mathbf{P}^i$) use the following:

**Lemma 4** Let $\mathbf{P}$ and $\mathbf{P}'$ be permutation matrices and $\mathbf{T}$ and $\mathbf{T}'$ be strictly lower triangular matrices, such that $\mathbf{P}^{\top}\mathbf{T}\mathbf{P} = \mathbf{P}'^{\top}\mathbf{T}'\mathbf{P}'$. Assume that $\mathbf{T}$ contains only non-zero elements in its strictly lower triangular part. Then, $\mathbf{P} = \mathbf{P}'$ and $\mathbf{T} = \mathbf{T}'$.

## A.2 Identifiability of LiNGAM

As a special case, we restate here the identifiability of the single-view LiNGAM model (in terms of the matrix $\mathbf{B}$).

**Theorem 5** (Identifiability of LiNGAM) In the statistical model defined by Eq. 1, the parameter $\mathbf{B}$ is identifiable, provided that the entries in $\mathbf{s}$ are mutually independent, that at most one of them is Gaussian, and that $\mathbf{B}$ is a DAG.

This proof is based on the one from Shimizu et al. [30], which we attempt to make a bit more rigorous. Identifiability was also shown in Shimizu [29, Section 2.3] in the case with 2 components.

Let us transform the SEM

$$ \mathbf{x} = \mathbf{B}\mathbf{x} + \mathbf{s} \quad (12) $$

into an ICA model

$$ \mathbf{x} = \mathbf{A}\mathbf{s} \quad (13) $$
where the mixing matrix $A$ is constrained as $A^{-1} = I - B$, and where $B$ is a DAG. The assumptions on $s$ are identical in SEM and ICA. Let us parameterize by $W = A^{-1} = I - B$, instead of $B$.

Consider two matrices $W = I - B$ and $W' = I - B'$, with $B$ and $B'$ being DAG matrices, that parameterize the same statistical model in Eq. 12. We have $det(W) = det(W') = 1$, so $W$ and $W'$ are invertible, which makes them valid unmixing matrices for the same ICA model in Eq. 13. So, from the identifiability theory of ICA [6], we know that there exist a sign-permutation matrix $Q$ and a scaling matrix $D$ (i.e. a diagonal matrix with positive entries on the diagonal) such that

$$W' = DQ^{\top}W \qquad (14)$$

hence

$$W'' = Q^{\top}W \qquad (15)$$

where we defined $W'' = D^{-1}W'$. We thus have that $W'' = D^{-1}(I - B')$ and $W = I - B$ both belong to the set $W$, so applying Lemma 3 to $W$ and $W''$ imposes $Q = I$, $D = I$, and $B' = B$. This makes $B$ fully identifiable (“fully” identifiable is in contrast to conventional ICA theory where some indeterminacies always remain), which concludes the proof.

## A.3 Proofs of technical lemmas

### A.3.1 Proof of Lemma 3

Consider $W, W' \in W$ and $Q$ a sign-permutation matrix. Suppose that

$$W' = Q^{\top}W. \qquad (16)$$

A permutation inequality Using the definition of $W$ in Eq. 10 and the decompositions of DAG matrices, we have

$$W = D^{-1}P^{\top}(I - T)P, \quad W' = D'^{-1}P'^{\top}(I - T')P' \qquad (17)$$

where $D, D'$ are diagonal matrices with positive entries on the diagonal, $P, P'$ are permutation matrices, and $T, T'$ are strictly lower triangular matrices. Denote $L = I - T$ and $L' = I - T'$; these are lower triangular matrices that have unit diagonals. We plug the decompositions into Eq. 16, and get

$$D'^{-1}P^{\top}L'P' = Q^{\top}D^{-1}P^{\top}LP. \qquad (18)$$

To exploit the structure of the lower triangular matrices $L$ and $L'$ and show how they constrain $Q$, we now switch notations from permutation — or sign-permutation — matrices ($P, P', Q$), to their corresponding permutation functions $(\phi, \phi', \psi)$. Eq. 18 thus yields

$$\frac{L'_{\phi'(i), \phi'(j)}}{D'_{ii}} = \pm \frac{L_{\phi(\psi(i)), \phi(j)}}{D_{\psi(i), \psi(i)}}, \quad \forall i, j \in [1, p] \qquad (19)$$

where the $\pm$ symbol comes from $Q$. In particular, $L'$ has unit diagonal, so

$$L_{\phi(\psi(i)), \phi(i)} = \pm \frac{D_{\psi(i), \psi(i)}}{D'_{ii}} \neq 0, \quad \forall i \in [1, p] \qquad (20)$$

and $L$ is lower triangular, so that its non-zero entries must satisfy

$$\phi(i) \le \phi(\psi(i)), \quad \forall i \in [1, p]. \qquad (21)$$

More generally, we can replace $i$ with $\psi(i)$, and so on, and obtain

$$\phi(i) \le \phi(\psi(i)) \le \phi(\psi^{(2)}(i)) \le \dots \quad \forall i \in [1, p] \qquad (22)$$

where the superscript denotes composition and where we can apply the permutation $\psi$ an arbitrary number of times.
**Q must be a sign matrix** Suppose that $\psi$ is not the identity. We can pick an index $k \in [1, p]$ where $k \neq \psi(k)$. Because $\phi$ and $\psi$ are injective, we can apply them to the inequality any number of times $n \in \mathbb{N}$ and get

$$ \phi(\psi^{(n)}(k)) \neq \phi(\psi^{(n+1)}(k)) \qquad (23) $$

which together with Eq. 22 implies

$$ \phi(k) < \phi(\psi(k)) < \phi(\psi^{(2)}(k)) < \dots \qquad (24) $$

which can be applied any number of times. However, each inequality increases the index by at least one, while the index cannot go above $p$. Thus, applying the chain at least $p$ times, we have a contradiction. So, the permutation $\psi$ in $Q$ must be the identity, and $Q$ boils down to a sign matrix, i.e. a diagonal matrix with 1 and -1 on the diagonal.

**Q must be the identity** Since $Q^\top = Q$, Eq. 18 can be reformulated as

$$ P'^{\top} L' P' = D' Q D^{-1} P'^{\top} L P \qquad (25) $$

where we now know that $D'QD^{-1}$ is a diagonal matrix. The matrix $P'^{\top}L'P'$ has a diagonal of ones and so too does $P'^{\top}LP$. It follows that $D'QD^{-1} = I$. Since $D$ and $D'$ have positive diagonal entries, it follows that the entries of $Q$ cannot equal -1. So, $Q = I$, and thus $D' = D$. This concludes the proof of the Lemma.

### A.3.2 Proof of Lemma 4

From the equality $P'^{\top}TP = P'^{\top}T'P'$, we deduce that

$$ T = \tilde{P}^{\top} T' \tilde{P} \qquad (26) $$

where $\tilde{P} = P'P'^{\top}$ is the permutation matrix represented by function $\phi$, which is defined such that: $\forall k \in [1, p]$, $\tilde{P}_{k,\phi(k)} = 1$ and $\forall l \neq \phi(k)$, $\tilde{P}_{k,l} = 0$. Using the notation $\phi$ instead of $\tilde{P}$, Eq. 26 can be rewritten as, $\forall k, l \in [1, p]$,

$$ T_{k,l} = T'_{\phi(k),\phi(l)} \qquad (27) $$

We proceed by a proof by contradiction: assume that $\tilde{P} \neq I$. As a consequence, there exists an index $k$ such that $\phi(k) \neq k$. Let us fix such an index as $k$. We can assume without loss of generality that

$$ \phi(k) > k \qquad (28) $$

otherwise we just have to invert signs and switch rows and columns in the following. By assumption, the strictly lower triangular part of $T$ only has non-zero elements, so we have

$$ T_{\phi(k),k} \neq 0 \qquad (29) $$

Yet, from Eq. 27, we know that $T_{\phi(k),k} = T'_{\phi(\phi(k)),\phi(k)}$ and all the non-zero elements of $T'$ are in its strictly lower triangular part, which implies

$$ \phi(\phi(k)) > \phi(k) \qquad (30) $$

This logic can be applied repeatedly as in the preceding lemma's proof, and we see that

$$ \phi(k) < \phi^{(2)}(k) < \phi^{(3)}(k) < \dots \qquad (31) $$

Now, this leads to an infinite strictly increasing sequence, which is contradictory since the index cannot grow greater than $p$. So, $\tilde{P} = I$, which means that $P'P'^{\top}P = P$, and thus, using the orthogonality of $P$,

$$ P' = P \qquad (32) $$

It also follows that $T' = T$, which concludes the proof of Lemma 4.
# B Likelihood-based criterion

## B.1 Basic setting

The fundamental question here is to choose the causal direction for two variables, having many views i. We denote the two variables by x and y to avoid too many indices. It could be that the causal direction is $x \to y$:

$$y^i = b_i x^i + e^i \quad \text{for all } i \qquad (33)$$

or $y \to x$:

$$x^i = c_i y^i + d^i \quad \text{for all } i. \qquad (34)$$

Recall that superscripts are for view in the case of the random variables; however, for the regression coefficients we use subscripts because we need to take squares of such quantities. For the direction $x \to y$, collect the regressor and residuals in the vectors $\mathbf{x} = (x^1, \dots, x^m)^\top$ and $\mathbf{e} = (e^1, \dots, e^m)^\top$, respectively; and likewise for $\mathbf{y}$ and $\mathbf{d}$ in the opposite direction. Second-order moment matrices are denoted by $\Sigma_x = \mathbb{E}[\mathbf{x}\mathbf{x}^\top]$, $\Sigma_y = \mathbb{E}[\mathbf{y}\mathbf{y}^\top]$, and likewise for $\mathbf{d}$ and $\mathbf{e}$. Note that the independence between $\mathbf{x}$ and $\mathbf{e}$ and between $\mathbf{y}$ and $\mathbf{d}$, and the fact that $\mathbf{x}$ and $\mathbf{y}$ are standardized, imply that

$$b_i = c_i = \operatorname{cov}(x_i, y_i) \quad \text{for all } i. \qquad (35)$$

## B.2 Derivation of the expected log-likelihood in one direction

We start by formulating the likelihood of the model in the case where disturbances are assumed to be Gaussian.

Consider the log-likelihood of $(\mathbf{x}, \mathbf{y})$ for the direction $x \to y$. Under the model, the prior (marginal) distribution of $\mathbf{x}$ is $\mathcal{N}(\mathbf{0}, \Sigma_x)$, and the prior distribution of $\mathbf{e} = \mathbf{y} - \mathbf{b} \odot \mathbf{x}$ is $\mathcal{N}(\mathbf{0}, \Sigma_e)$ — where $\odot$ denotes the elementwise product. Thus, we have

$$\begin{align}
\log p(\mathbf{x}, \mathbf{y}; \mathbf{b}, \Sigma_x, \Sigma_e, x \to y) &= \log p(\mathbf{y}|\mathbf{x}; \mathbf{b}, \Sigma_x, \Sigma_e, x \to y) + \log p(\mathbf{x}; \mathbf{b}, x \to y) \nonumber \\
&= -(\mathbf{y} - \mathbf{b} \odot \mathbf{x})^T \Sigma_e^{-1} (\mathbf{y} - \mathbf{b} \odot \mathbf{x}) - \log |\Sigma_e| - \mathbf{x}^T \Sigma_x^{-1} \mathbf{x} - \log |\Sigma_x| \nonumber \\
&= -\operatorname{tr}(\Sigma_e^{-1}(\mathbf{y} - \mathbf{b} \odot \mathbf{x})(\mathbf{y} - \mathbf{b} \odot \mathbf{x})^\top) - \log |\Sigma_e| - \operatorname{tr}(\Sigma_x^{-1}\mathbf{x}\mathbf{x}^\top) - \log |\Sigma_x| \tag{36}
\end{align}$$

where $\operatorname{tr}$ is the trace operator, and where we neglect the part of the normalization constant that does not depend on any parameters.

The expected log-likelihood is thus

$$\mathbb{E}[\log p(\mathbf{x}, \mathbf{y}; \mathbf{b}, \Sigma_x, \Sigma_e, x \to y)] = -\operatorname{tr}(\Sigma_e^{-1}\mathbb{E}[(\mathbf{y} - \mathbf{b} \odot \mathbf{x})(\mathbf{y} - \mathbf{b} \odot \mathbf{x})^\top]) - \log |\Sigma_e| \qquad (37)$$

$$-\operatorname{tr}(\Sigma_x^{-1}\mathbb{E}[\mathbf{x}\mathbf{x}^\top]) - \log |\Sigma_x| \qquad (38)$$

$$= -\log |\Sigma_x| - \log |\Sigma_e| - 2m. \qquad (39)$$

Denote by $\mathcal{L}$ this expected log-likelihood. We have (up to constants):

$$\mathcal{L}(\Sigma_x, \Sigma_e; x \to y) = -\log |\Sigma_x| - \log |\Sigma_e|. \qquad (40)$$

## B.3 Derivation of the likelihood-based criterion

We now derive a criterion that can be used to find the causal direction.

The expected log-likelihood of the model for the direction $y \to x$ can be obtained by switching the roles of $x$ and $y$ in Eq. 40. Thus, we have

$$\mathcal{L}(\Sigma_y, \Sigma_d; y \to x) = -\log |\Sigma_y| - \log |\Sigma_d| \qquad (41)$$

where we recall that $\Sigma_y = \mathbb{E}[\mathbf{y}\mathbf{y}^\top]$ and $\Sigma_d = \mathbb{E}[\mathbf{d}\mathbf{d}^\top]$.
The correct direction can then be found by comparing the expected log-likelihoods of both directions. In particular, we calculate the difference between the two log-likelihoods, which can be reformulated as a likelihood ratio (LR):

$$LR := \mathcal{L}(\Sigma_x, \Sigma_e; x \to y) - \mathcal{L}(\Sigma_y, \Sigma_d; y \to x) = \mathbb{E} \left[ \log \frac{p(x, y; b, \Sigma_x, \Sigma_e, x \to y)}{p(x, y; c, \Sigma_y, \Sigma_d, y \to x)} \right] . \quad (42)$$

Using Eq. 40 and Eq. 41, the criterion LR becomes

$$LR = -\log |\Sigma_x| - \log |\Sigma_e| + \log |\Sigma_y| + \log |\Sigma_d| . \quad (43)$$

Intuitively, the criterion in Eq. 43 will be large if the $x^i$ are highly correlated (more than the $y^i$) over the views, since then the first determinant will be small. This should make intuitive sense since it means the cause is highly structured. Likewise, if the residuals for $x \to y$ are highly structured, that direction is more likely. In other words, the direction that exhibits more structure wins.

Surprisingly, in this formulation, there is actually no particular model for the dependencies of the disturbances: the likelihood should simply prefer dependent disturbances. It is important to note that the criterion makes sense even if the distribution is non-Gaussian, as will be proven below. We thus obtain a method that works for data of any distribution, but using only second-order statistics.

## B.4 Non-negativity of the likelihood-based criterion

In this section, we prove that the criterion in Eq. 43 is non-negative for the direction $x \to y$. Importantly, no assumption of Gaussianity needs to be made here.

Consider the direction $x \to y$, and let us write $\mathbf{B} := \text{diag}(b_1, \dots, b_m)$. Using the matrix $\mathbf{B}$ rather than the vector $\mathbf{b}$ makes notations clearer in the proofs. In the following, we use interchangeably $\mathbf{b}$ and $\mathbf{B}$. The covariance matrix of $\mathbf{y} = \mathbf{B}\mathbf{x} + \mathbf{e}$ is

$$\Sigma_y = \mathbf{B}\Sigma_x\mathbf{B} + \Sigma_e . \quad (44)$$

Given that $b_i = c_i$ for all $i$, we have $\mathbf{d} = \mathbf{x} - \mathbf{B}\mathbf{y} = \mathbf{x} - \mathbf{B}(\mathbf{B}\mathbf{x} + \mathbf{e}) = (\mathbf{I} - \mathbf{B})^2\mathbf{x} - \mathbf{B}\mathbf{e}$. Thus, the covariance matrix of $\mathbf{d}$ is

$$\Sigma_d = (\mathbf{I} - \mathbf{B})^2 \Sigma_x (\mathbf{I} - \mathbf{B})^2 + \mathbf{B}\Sigma_e\mathbf{B} . \quad (45)$$

Using these expressions, the criterion in Eq. 43 becomes

$$LR = -\log |\Sigma_x| - \log |\Sigma_e| + \log |\mathbf{B}\Sigma_x\mathbf{B} + \Sigma_e| + \log |(\mathbf{I} - \mathbf{B})^2\Sigma_x(\mathbf{I} - \mathbf{B})^2 + \mathbf{B}\Sigma_e\mathbf{B}| . \quad (46)$$

Collect the variances of the residuals $\mathbf{e}$ in the diagonal matrix

$$\mathbf{L} := \text{diag} \left( \sqrt{\text{var}(e^1)}, \dots, \sqrt{\text{var}(e^m)} \right) \quad (47)$$

and define the correlation matrix of the residuals as $\tilde{\Sigma}_e$, so that

$$\Sigma_e = \mathbf{L}\tilde{\Sigma}_e\mathbf{L} . \quad (48)$$

Now, the objective becomes

$$\begin{aligned} LR = & -\log |\Sigma_x| - \log |\tilde{\Sigma}_e| - \log |\mathbf{L}^2| + \log |\mathbf{B}\Sigma_x\mathbf{B} + \mathbf{L}\tilde{\Sigma}_e\mathbf{L}| \\ & + \log |(\mathbf{I} - \mathbf{B})^2\Sigma_x(\mathbf{I} - \mathbf{B})^2 + \mathbf{B}\mathbf{L}\tilde{\Sigma}_e\mathbf{L}\mathbf{B}| . \end{aligned} \quad (49)$$
Importantly, the fact that $x$ and $y$ are standardized means that $\Sigma_x$ and $\Sigma_y$ have unit diagonal, so Eq. 44 implies

$$L^2 + B^2 = I \qquad (50)$$

and it follows that the diagonal entries of $L$ and $B$ are in $[-1, 1]$. Thus, we have

$$\text{LR} = -\log |\Sigma_x| - \log |\tilde{\Sigma}_e| - \log |L^2| + \log |B\Sigma_x B + L\tilde{\Sigma}_e L| + \log |L^2\Sigma_x L^2 + LB\tilde{\Sigma}_e BL| \qquad (51)$$

$$= -\log |\Sigma_x| - \log |\tilde{\Sigma}_e| + \log |B\Sigma_x B + L\tilde{\Sigma}_e L| + \log |L\Sigma_x L + B\tilde{\Sigma}_e B| \qquad (52)$$

where $L$ and $B$ commuted in the first equality because they are diagonal. The following lemma (proven in Appendix B.5) shows that this criterion is non-negative, and even positive under some assumptions.

**Lemma 6** Let $\Sigma_x$ and $\tilde{\Sigma}_e$ be $m \times m$ real symmetric positive definite matrices with unit diagonal entries. Let $B = \text{diag}(b_1, \dots, b_m)$ and $L = \text{diag}(l_1, \dots, l_m)$ be diagonal matrices with entries in $[-1, 1]$, satisfying $B^2 + L^2 = I$ (i.e. $b_i^2 + l_i^2 = 1$ for all $i$). Define the function

$$J(B, L) := -\log \det \Sigma_x - \log \det \tilde{\Sigma}_e + \log \det (B\Sigma_x B + L\tilde{\Sigma}_e L) + \log \det (L\Sigma_x L + B\tilde{\Sigma}_e B). \qquad (53)$$

Under the above assumptions, we have $J(B, L) \ge 0$. Moreover, $J(B, L) = 0$ if and only if

$$L \tilde{\Sigma}_e B = B \Sigma_x L. \qquad (54)$$

## B.5 Proof of a useful lemma

Let us prove Lemma 6.

**Step 1: define A and C — Define**

$$A := B\Sigma_x B + L\tilde{\Sigma}_e L \quad \text{and} \quad C := L\Sigma_x L + B\tilde{\Sigma}_e B. \qquad (55)$$

We have

$$J(B, L) := -\log \det \Sigma_x - \log \det \tilde{\Sigma}_e + \log \det A + \log \det C. \qquad (56)$$

**Step 2: define M — Define the $2m \times 2m$ block matrix**

$$M := \begin{pmatrix} B & L \\ -L & B \end{pmatrix}. \qquad (57)$$

Since $B$ and $L$ are diagonal, they commute, hence $BL = LB$. Using $B^2 + L^2 = I$, we compute

$$MM^\top = \begin{pmatrix} BB + LL & -BL + LB \\ -LB + BL & LL + BB \end{pmatrix} = \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix} = I_{2m}. \qquad (58)$$

Thus $M$ is orthonormal, so $\det M = \pm 1$.

**Step 3: define B — Consider the block-diagonal matrix $\text{diag}(\Sigma_x, \tilde{\Sigma}_e)$ and conjugate by $M$:**

$$B := M \begin{pmatrix} \Sigma_x & 0 \\ 0 & \tilde{\Sigma}_e \end{pmatrix} M^\top = \begin{pmatrix} A & S \\ S^\top & C \end{pmatrix}, \quad S := L \tilde{\Sigma}_e B - B \Sigma_x L. \qquad (59)$$

Because $M$ is orthonormal and $\Sigma_x, \tilde{\Sigma}_e > 0$, the matrix $B$ is symmetric positive definite (SPD). In particular, $A$ and $C$ are also SPD (as they are principal blocks of a SPD matrix) and thus $A^{-1}$ exists.
**Step 4:** Schur complement of A in B — For a symmetric block matrix $\begin{pmatrix} A & S \\ S^\top & C \end{pmatrix}$ with $A > 0$, the Schur complement of $A$ is $C - S^\top A^{-1} S$, and the Schur's formula yields

$$ \det\begin{pmatrix} A & S \\ S^\top & C \end{pmatrix} = \det(A) \det(C - S^\top A^{-1} S) . \qquad (60) $$

Applying this to $\mathcal{B}$ and using multiplicativity of determinant and the fact that $\det M = \pm 1$ gives

$$ \det(\Sigma_x) \det(\tilde{\Sigma}_e) = \det(\mathcal{B}) = \det(A) \det(C - S^\top A^{-1} S) . \qquad (61) $$

**Step 5:** derive a matrix inequality — Because $\mathcal{B} > 0$ and $A > 0$, it follows from the Schur complement characterization of positive definiteness that $C - S^\top A^{-1} S > 0$. In particular, we have $\det(C - S^\top A^{-1} S) > 0$. Moreover, since $A > 0$, we have $S^\top A^{-1} S \ge 0$, and thus $C - S^\top A^{-1} S \le C$.

**Step 6:** derive a determinant inequality — If $0 \preceq X \preceq Y$ and $Y > 0$, then conjugating by $Y^{-1/2}$ shows $0 \preceq Y^{-1/2} X Y^{-1/2} \preceq I$. Thus, all eigenvalues of $Y^{-1/2} X Y^{-1/2}$ lie in $[0, 1]$, so their product satisfies $\det(Y^{-1/2} X Y^{-1/2}) \le 1$. Hence $\det(X) \le \det(Y)$. Applying this to $X := C - S^\top A^{-1} S$ and $Y := C$ yields

$$ \det(C - S^\top A^{-1} S) \le \det(C) . \qquad (62) $$

**Step 7:** $J(B, L) \ge 0$ — From Steps 4 and 6, we get

$$ \det(\Sigma_x) \det(\tilde{\Sigma}_e) = \det(A) \det(C - S^\top A^{-1} S) \le \det(A) \det(C) . \qquad (63) $$

Taking natural logarithms and rearranging gives

$$ \log \det(A) + \log \det(C) \ge \log \det(\Sigma_x) + \log \det(\tilde{\Sigma}_e) , \qquad (64) $$

which is precisely $J(B, L) \ge 0$.

**Step 8:** equality condition — The equality $J(B, L) = 0$ holds iff

$$ \det(C - S^\top A^{-1} S) = \det(C) . \qquad (65) $$

Write $K := C^{-1/2}(S^\top A^{-1} S)C^{-1/2} \ge 0$. We have

$$ \det(C - S^\top A^{-1} S) = \det(C) \det(I - C^{-1} S^\top A^{-1} S) \qquad (66) $$

$$ = \det(C) \det(C^{-1/2} C^{1/2} - C^{-1/2} C^{-1/2} S^\top A^{-1} S C^{-1/2} C^{1/2}) \qquad (67) $$

$$ = \det(C) \det(C^{-1/2}) \det(I - K) \det(C^{1/2}) \qquad (68) $$

$$ = \det(C) \det(I - K) . \qquad (69) $$

Thus, the equality condition is $\det(I - K) = 1$. From Step 5, we know that $0 \preceq C - S^\top A^{-1} S \preceq C$. Conjugating by $C^{-1/2}$ gives $0 \preceq I - K \preceq I$, so each eigenvalue $\lambda_i$ of $K$ satisfies $0 \le \lambda_i \le 1$. The eigenvalues of $I - K$ are $1 - \lambda_i \in [0, 1]$, hence

$$ \prod_i (1 - \lambda_i) = \det(I - K) = 1 . \qquad (70) $$

The equality thus holds iff $\lambda_i = 0$ for all $i$, hence $K = 0$. Since $C^{-1/2}, A^{-1} > 0$, we get $S = 0$. Conversely, $S = 0$ clearly forces equality. Therefore, the equality occurs exactly when

$$ S = L \tilde{\Sigma}_e B - B \Sigma_x L = 0 , \qquad (71) $$
i.e.

$$ \forall i, j: \quad b_i l_j(\Sigma_x)_{ij} = l_i b_j(\tilde{\Sigma}_e)_{ij}. \qquad (72) $$

This completes the proof.

## B.6 Conditions for strict positivity of the likelihood-based criterion

From now on, we assume that the model is not degenerate in the sense that disturbances have a positive variance, that is $l_i = \text{var}(e_i) > 0$ for all i.

The criterion LR derived in Appendix B.4 is used to determine the causal direction from its sign, with LR = 0 indicating that the two directions cannot be distinguished. Lemma 6 shows that LR = 0 if and only if

$$ L \tilde{\Sigma}_e B = B \Sigma_x L . \qquad (73) $$

Entrywise, this is equivalent to

$$ \forall i, j: \quad l_i b_j(\tilde{\Sigma}_e)_{ij} = b_i l_j(\Sigma_x)_{ij}. \qquad (74) $$

Let us fix $i, j \in [1, m]$.

**Case 1:** $b_i = 0$ and $b_j = 0$ — The equality immediately holds.

**Case 2:** $b_i = 0$ and $b_j \neq 0$ — The right-hand side in Eq. 74 is equal to 0. Since $l_i = \text{var}(e^i) > 0$ and $b_j \neq 0$, we must have $(\tilde{\Sigma}_e)_{ij} = 0$.

**Case 3:** $b_i \neq 0$ and $b_j = 0$ — The left-hand side in Eq. 74 is equal to 0. Since $l_j = \text{var}(e^j) > 0$ and $b_i \neq 0$, we must have $(\Sigma_x)_{ij} = 0$.

**Case 4:** $b_i \neq 0$ and $b_j \neq 0$ — Multiplying Eq. 74 by $\frac{l_j}{b_i b_j^2}$, and using $l_j^2 = 1 - b_j^2$, we obtain

$$ \frac{l_i l_j}{b_i b_j} (\tilde{\Sigma}_e)_{ij} = \left( \frac{1}{b_j^2} - 1 \right) (\Sigma_x)_{ij}. \qquad (75) $$

Switching the roles of i and j, and using the symmetry of $(\tilde{\Sigma}_e)_{ij}$ and $(\Sigma_x)_{ij}$, yield

$$ \frac{l_i l_j}{b_i b_j} (\tilde{\Sigma}_e)_{ij} = \left( \frac{1}{b_i^2} - 1 \right) (\Sigma_x)_{ij}. \qquad (76) $$

hence

$$ b_i^2(\Sigma_x)_{ij} = b_j^2(\Sigma_x)_{ij}. \qquad (77) $$

Similarly, exchanging the roles of $(b_i, b_j, (\Sigma_x)_{ij})$ and $(l_i, l_j, (\tilde{\Sigma}_e)_{ij})$ yields

$$ l_i^2(\tilde{\Sigma}_e)_{ij} = l_j^2(\tilde{\Sigma}_e)_{ij}. \qquad (78) $$

Now suppose $(\Sigma_x)_{ij} \neq 0$ or $(\tilde{\Sigma}_e)_{ij} \neq 0$. From Eq. 77 and Eq. 78 it follows that $b_i^2 = b_j^2$ or $l_i^2 = l_j^2$. Since $b_i^2 + l_i^2 = 1$, both equalities hold simultaneously:

$$ |b_i| = |b_j| \quad \text{and} \quad l_i = l_j \qquad (79) $$

where the absolute value can be dropped for $l_i$ because $l_i = \text{var}(e_i) > 0$. Plugging these equalities into Eq. 75 yields

$$ l_j^2(\tilde{\Sigma}_e)_{ij} = \pm(1 - b_j^2)(\Sigma_x)_{ij} \qquad (80) $$
hence

$$ |(\tilde{\Sigma}_e)_{ij}| = |(\Sigma_x)_{ij}| . \qquad (81) $$

Moreover, for Eq. 74 to hold we must also have the sign equality

$$ \operatorname{sign}(b_i) \operatorname{sign}(b_j) = \operatorname{sign}((\Sigma_x)_{ij}) \operatorname{sign}((\tilde{\Sigma}_e)_{ij}) . \qquad (82) $$

Altogether, either both $(\Sigma_x)_{ij}$ and $(\tilde{\Sigma}_e)_{ij}$ vanish, or the following four conditions hold:

$$ |(\tilde{\Sigma}_e)_{ij}| = |(\Sigma_x)_{ij}| \neq 0 \qquad (83) $$

$$ |b_i| = |b_j| \qquad (84) $$

$$ l_i = l_j \qquad (85) $$

$$ \operatorname{sign}(b_i) \operatorname{sign}(b_j) = \operatorname{sign}((\Sigma_x)_{ij}) \operatorname{sign}((\tilde{\Sigma}_e)_{ij}) . \qquad (86) $$

Conversely, if $(\Sigma_x)_{ij} = (\tilde{\Sigma}_e)_{ij} = 0$, or if conditions Eq. 83-Eq. 86 are satisfied, then Eq. 74 holds.

**Conclusion** — By contraposition, Eq. 73 is violated whenever there exists $(i, j)$ such that one of the following cases holds:

* $b_i = 0, b_j \neq 0$, and $(\tilde{\Sigma}_e)_{ij} \neq 0$

* $b_i \neq 0, b_j = 0$, and $(\Sigma_x)_{ij} \neq 0$

* $b_i \neq 0, b_j \neq 0$, $(\Sigma_x)_{ij} \neq 0$ or $(\tilde{\Sigma}_e)_{ij} \neq 0$, and at least one of the conditions Eq. 83-Eq. 86 is not met.

In particular, a simple sufficient condition is

$$ \exists i, j : \quad (b_i \neq 0 \quad \text{or} \quad b_j \neq 0) \quad \text{and} \quad |(\tilde{\Sigma}_e)_{ij}| \neq |(\Sigma_x)_{ij}| . \qquad (87) $$
# C Cross-covariance-based criterion

Consider the setting described in Appendix B.1.

## C.1 Derivation of the cross-covariance-based criterion

A simple approach to infer the causal direction is to exploit the assumption that root variables and residuals are independent. Specifically, if the true direction is $x \to y$, then $x$ and $e$ are independent; if the true direction is $y \to x$, then $y$ and $d$ are independent. This implies that the cross-moment matrix $\mathbb{E}[ex^\top]$ should vanish when $x \to y$, while $\mathbb{E}[dy^\top]$ should vanish when $y \to x$. Comparing the Frobenius norms of these matrices therefore provides a natural criterion, based solely on cross-covariances:

$$ \mathrm{FC} := \left\| \mathbb{E}[dy^\top] \right\|_F - \left\| \mathbb{E}[ex^\top] \right\|_F . \qquad (88) $$

Intuitively, if this criterion is positive we deduce that $x \to y$, and if negative we deduce that $y \to x$.

## C.2 Conditions for strict positivity of the cross-covariance criterion

Assume without loss of generality that the true causal direction is $x \to y$, and let us analyze the sign of the criterion FC.

Given that $x \to y$, the true model is

$$ y = Bx + e \qquad (89) $$

where $B = \mathrm{diag}(\mathrm{cov}(x^1, y^1), \dots, \mathrm{cov}(x^m, y^m))$, and $e$ is independent from $x$. So, we have

$$ \mathbb{E}[xe^\top] = 0. \qquad (90) $$

In the wrong direction, the residuals take the form

$$ d = x - Cy, \quad C = B. \qquad (91) $$

The cross-covariance between observations $y$ and residuals $d$ is then

$$ \mathbb{E}[yd^\top] = \mathbb{E}[y(x - By)^\top] \qquad (92) $$

$$ = \mathbb{E}[(Bx + e)(x - B(Bx + e))^\top] \qquad (93) $$

$$ = \mathbb{E}[Bxx^\top(I - B^2)] + \mathbb{E}[ex^\top(I - B^2)] - \mathbb{E}[Bxe^\top B] - \mathbb{E}[ee^\top B] \qquad (94) $$

$$ = B \Sigma_x (I - B^2) - \Sigma_e B \qquad (95) $$

where $\Sigma_x = \mathbb{E}[xx^\top]$ and $\Sigma_e = \mathbb{E}[ee^\top]$.

From Eq. 90 and the non-negativity of the Frobenius norm, we see that FC is always non-negative. Moreover, using Eq. 95, we deduce that $FC = 0$ if and only if

$$ B \Sigma_x (I - B^2) - \Sigma_e B = 0. \qquad (96) $$

Note that the fact that $l_i > 0$, for all $i$, makes $L$ invertible. So, using the identity $I - B^2 = L^2$, and multiplying on the right by $L^{-1}$, we obtain

$$ B \Sigma_x L = \Sigma_e B L^{-1} \qquad (97) $$

hence

$$ B \Sigma_x L = L \tilde{\Sigma}_e B \qquad (98) $$

since $\Sigma_e = L \tilde{\Sigma}_e L$, and $B$ and $L$ commute.

This condition coincides exactly with Eq. 73. In particular, this shows that the likelihood-based criterion and the cross-covariance-based criterion share the same consistency assumptions.
# D Identifiability of LiMVAM

## D.1 Multivariate conditions for strict positivity of a criterion

The following assumption is a multivariate generalization of the assumption detailed in B.6 and C.2.

**Assumption 6** (Diversity of the views) Let $j \neq j'$ denote any two component indices such that $\mathbf{x}_j \to \mathbf{x}_{j'}$. Then, there exist two views i and i' such that the three following conditions hold:

1. $\text{corr}(\mathbf{x}_j^i, \mathbf{x}_{j'}^{i'}) \neq 0$ or $\text{corr}(\mathbf{x}_j^{i'}, \mathbf{x}_{j'}^{i'}) \neq 0$

2. $\text{corr}(\mathbf{x}_j^i, \mathbf{x}_j^{i'}) \neq 0$ or $\text{corr}(\mathbf{e}_{j',j}^i, \mathbf{e}_{j',j}^{i'}) \neq 0$

3. At least one of the four following inequalities is met:

$$|\text{corr}(\mathbf{x}_j^i, \mathbf{x}_j^{i'})| \neq |\text{corr}(\mathbf{e}_{j',j}^i, \mathbf{e}_{j',j}^{i'})|, |\text{corr}(\mathbf{x}_j^i, \mathbf{x}_{j'}^{i'})| \neq |\text{corr}(\mathbf{x}_j^{i'}, \mathbf{x}_{j'}^{i'})|, \text{var}(\mathbf{e}_{j',j}^i) \neq \text{var}(\mathbf{e}_{j',j}^{i'}) \\ \text{sign}(\text{corr}(\mathbf{x}_j^i, \mathbf{x}_{j'}^{i'})) \neq \text{sign}(\text{corr}(\mathbf{x}_j^{i'}, \mathbf{x}_{j'}^{i'})) \neq \text{sign}(\text{corr}(\mathbf{e}_{j',j}^i, \mathbf{e}_{j',j}^{i'}))$$

where corr denotes the correlation coefficient and sg denotes the sign function.

These conditions can only be met for indices $i \neq i'$, so they can be interpreted as a need for correlation and diversity across views. Because they may be hard to interpret at first glance, we next visualize in Figure 4 what these conditions mean for the causal graph.

In particular, Assumption 6 immediately implies its simpler version used in the main text, Assumption 1.

## D.2 An example for interpreting the SOS-based identifiability condition in Assumption 1

Consider one view and two variables denoted for notational simplicity $x$ and $y$. We have

$$y^1 = b^1 x^1 + e^1 \qquad (99)$$

and assume for simplicity that $x^1$ and $y^1$ have been standardized.

Now, suppose we observe two datasets from this one view, denote them by $(x^1, y^1)$ and $(x^{1'}, y^{1'})$, which have identical distributions and are independent of each other. Crucially, all this data is from one view. However, we could artificially generate a new “view” by adding those two datasets together:

$$x^2 = x^1 + x^{1'} \qquad (100)$$

and likewise for $y^2$ and $e^2$. Clearly the new data follows the SEM with the same $b$:

$$y^2 = b^1 x^2 + e^2 \qquad (101)$$

since we are simply adding each of the terms separately in the two datasets.

But this “multi-view” data is degenerate. In fact, the correlation coefficients are

$$\text{corr}(x^1, x^2) = \frac{\text{cov}(x^1, x^1) + \text{cov}(x^1, x^{1'})}{\text{std}(x^1)\text{std}(x^2)} = \frac{\text{var}(x^1)}{\sqrt{2}\text{std}(x^1)\text{std}(x^1)} = \frac{1}{\sqrt{2}}. \qquad (102)$$

An exactly identical calculation applies for $\text{corr}(e^1, e^2)$ which is also equal to $1/\sqrt{2}$.

Thus, we see that this case violates the simple identifiability condition in Assumption 1. This is understandable since no new information is created when computing the variable in the second view, and the two views are thus degenerate.
<figure><img src="image_6.png" /><figcaption>Figure 4: Assume the direction is $x_j \rightarrow x_{j'}$. Assumption 6 corresponds to being in one of the 11 situations showed in the figure. Note that no arrow means that no information is required.</figcaption></figure>
## D.3 Proof of Theorem 1

We proceed in three steps. (i) We show that the LiMVAM model always admits at least one root variable and that our LR and FC criteria recover one consistently. (ii) We prove a recursive residuals lemma: removing a recovered root and regressing on it preserves the model form and the causal ordering, which implies that the full (partial) ordering is consistently recoverable. (iii) Conditional on the recovered ordering, each view-specific strictly lower-triangular system can be estimated consistently, hence the adjacency matrices are identifiable.

**Preliminaries** For clarity, we unpack the “shared ordering” assumption. This assumption is equivalent to the existence of a main DAG $G$ that is the union of view-specific DAGs ${G^i}_{i=1}^m$: an edge appears in $G$ if it appears in at least one $G^i$. Root variables are therefore common candidates across views. Furthermore, in terms of vocabulary, we distinguish (a) directed paths (a sequence of edges with consistent direction) from (b) direct edges (one edge between adjacent variables). An indirect path is a directed path with at least two edges.

### D.3.1 Recovering a root variable

Let observations follow LiMVAM in Eq. 3. Since all views share an ordering, there exists a (possibly non-unique) permutation $P$ such that

$$T^i = P^\top B^i P \quad \text{is strictly lower triangular for all } i \in [1, m]. \qquad (103)$$

Reordering the observations and disturbances with $x'^i = P x^i$ and $e'^i = P e^i$, the model is

$$x'^i = T^i x'^i + e'^i, \quad i \in [1, m]. \qquad (104)$$

The first row of each $T^i$ is zero, so $x'^i_1$ is exogenous; hence a root variable always exists. To find a root, consider distinct indices $(j, k)$ and define $M_{jk}$ to be the criterion from Eq. 5 or Eq. 6 (with $M_{kj} = -M_{jk}$). Stack these into $M \in \mathbb{R}^{p \times p}$ with zero diagonal, and analyze signs under Assumption 6 (which concerns only direct relations):

1. **Direct edge:** $x_j \to x_k$ or $x_k \to x_j$. Under the bivariate conditions derived from Assumption 6, $M_{jk} > 0 > M_{kj}$ when $x_j \to x_k$, and $M_{jk} < 0 < M_{kj}$ when $x_k \to x_j$.

2. **Indirect path:** $x_j \to \dots \to x_k$ or $x_k \to \dots \to x_j$. Consider the reduced form of the model in Eq. 4, $x^i = A^i e^i$ with $A^i = (I - B^i)^{-1}$. The entry $A_{kj}^i$ is the total causal effect of $x_j^i$ on $x_k^i$, i.e. the sum over all directed paths from $j$ to $k$ of the products of edge coefficients. Hence for each view $i$,

$$x_k^i = A_{kj}^i x_j^i + \epsilon^i \qquad (105)$$

where $\epsilon^i$ collects all terms not caused by $x_j^i$ and is independent of $x_j^i$. Thus, the pair $(x_j^i, x_k^i)$ behaves as if there were a direct edge with coefficient $A_{kj}^i$, and the sign analysis of Item 1 applies but with non-strict inequalities: $M_{jk} \ge 0 \ge M_{kj}$ when $x_j \to \dots \to x_k$ (and conversely). If the bivariate conditions hold for $(j, k)$, the inequalities are strict; however, we do not require this here.

3. **No relation, no common ancestor:** then $x_j$ and $x_k$ depend on disjoint disturbance sets and are independent, hence $M_{jk} = M_{kj} = 0$.

4. **No relation but with a common ancestor.** If $x_j$ and $x_k$ are not connected by any directed path but share at least one ancestor, then their pairwise criteria $M_{jk}$ and $M_{kj}$ may take mixed signs. However, since each of $x_j$ and $x_k$ has at least one parent, there exist variables $x_{j'}$ and $x_{k'}$ such that $M_{jj'} < 0$ and $M_{kk'} < 0$. Thus, even if the signs of $M_{jk}$ and
$M_{kj}$ are not determined, the j-th and k-th rows of $\mathcal{M}$ necessarily contain negative entries. For instance, with three variables, $\mathbf{x}_1 \to \mathbf{x}_2$ and $\mathbf{x}_1 \to \mathbf{x}_3$ makes $\mathbf{x}_2$ and $\mathbf{x}_3$ siblings: there is no directed path between them, but both rows 2 and 3 of $\mathcal{M}$ contain negative entries due to their direct parent $\mathbf{x}_1$.

Consequently, $\mathbf{x}_j$ is a root (no incoming edges) if and only if the j-th row of $\mathcal{M}$ has only nonnegative entries. Hence a root exists and is consistently detectable. If multiple roots exist, we pick one at random.

### D.3.2 Recovering the causal ordering

We now show the recursive step.

**Lemma 7** (Residuals preserve LiMVAM and the ordering) Let $j$ be a root variable. Regress every other variable on $\mathbf{x}_j$ (within each view) and remove $\mathbf{x}_j$. The residual vector still follows a LiMVAM model with the same ordering over the remaining variables.

*Proof.* The proof is a direct adaptation of Shimizu et al. [32, Lemma 2 and Corollary 1] to the LiMVAM model: the original result is in the single-view setting with non-Gaussian disturbances. Choose a permutation $\mathcal{P}$ that places $\mathbf{x}_j$ first for simplicity (since $\mathbf{x}_j$ is a root variable, this permutation exists). In the permuted system we write

$$ \mathbf{x}^{i} = \mathcal{T}^{i} \mathbf{x}^{i} + \mathbf{e}^{i}, \quad \mathbf{x}^{i} = \mathcal{P} \mathbf{x}^{i}, \quad \mathbf{e}^{i} = \mathcal{P} \mathbf{e}^{i} \qquad (106) $$

where each $\mathcal{T}^i$ is strictly lower triangular and $\mathbf{x}'_1 = \mathbf{x}_j$. Equivalently, the model can be expressed in reduced form

$$ \mathbf{x}^{i} = \mathcal{A}^{i} \mathbf{e}^{i}, \quad \mathcal{A}^{i} := (\mathcal{I} - \mathcal{T}^{i})^{-1} \qquad (107) $$

where the matrix $\mathcal{A}^i$ is lower triangular with unit diagonal. The entry $A_{k1}^i$ ($k \ge 2$) quantifies the dependence of $x_k^{i}$ on $e_1^{i}$. Since $x_1^{i} = e_1^{i}$ is exogenous, we have

$$ \mathbb{E}[x_k^{i} x_1^{i}] = A_{k1}^i \mathbb{E}[(x_1^{i})^2] \qquad (108) $$

so the least-squares regression coefficient of $x_k^{i}$ on $x_1^{i}$ equals $A_{k1}^i$. Therefore, regressing each $x_k^{i}$ ($k \ge 2$) on $x_1^{i}$ yields residuals $\mathbf{r}_{(1)}^i \in \mathbb{R}^{p-1}$ such that

$$ \mathbf{r}_{(1)}^i = \mathbf{x}_{2:p}^i - \mathcal{A}_{2:p,1}^i \mathbf{x}_1^{i} = \mathcal{A}_{2:p,1}^i \mathbf{e}^{i} - \mathcal{A}_{2:p,1}^i \mathbf{e}_1^{i} = \mathcal{A}_{2:p,2:p}^i \mathbf{e}_{2:p}^{i} \qquad (109) $$

Here $\mathcal{A}_{2:p,2:p}^i$ is again lower triangular with unit diagonal, and $\mathbf{e}_{2:p}^{i}$ has mutually independent components. Hence $\mathbf{r}_{(1)}^i$ follows a LiMVAM model on $p-1$ variables. Moreover, since $\mathcal{A}_{2:p,2:p}^i$ is precisely the submatrix of $\mathcal{A}^i$ obtained by deleting the first row and column, the relative structure among the remaining variables is unchanged. In other words, the causal ordering of variables 2, ..., p is preserved. $\square$

Applying Lemma 7 recursively — find a root via $\mathcal{M}$ (under Assumption 6), regress out its effect, remove it, and repeat — recovers the full (partial) causal ordering consistently.

### D.3.3 Recovering the adjacency matrices

Sections D.3.1 and D.3.2 showed that, under Assumption 6, the algorithms PairwiseLiMVAM and DirectLiMVAM consistently recover the causal ordering. Equivalently, they recover a permutation matrix $\mathcal{P}$ that makes the $\mathcal{T}^i = \mathcal{P}^\top \mathcal{B}^i \mathcal{P}$ strictly lower triangular. So, we can reorder the variables as follows:

$$ \mathbf{x}^{i} = \mathcal{T}^{i} \mathbf{x}^{i} + \mathbf{e}^{i}, \quad i \in [1, m] \qquad (110) $$
Fix $j \in [1, p]$. We have

$$x_j^i = \sum_{k=1}^{j-1} T_{jk}^i x_k^i + e_j^i, \quad i \in [1, m]. \qquad (111)$$

Because $e_j^i$ is independent of $(x_1^i, \dots, x_{j-1}^i)$, the regressors are exogenous. This is a system of linear regressions across views with potentially cross-view-correlated errors. Feasible GLS is consistent in this setting: once a consistent estimate of the cross-view error covariance is obtained, the GLS estimator converges to the true coefficients [11]. Therefore, $\{T_{jk}^i\}_{k<j}$ are consistently estimated for every $(i, j)$, hence each $T^i$ is consistently recovered. Finally, $B^i = PT^iP^\top$, so the view-specific adjacency matrices are identifiable. This concludes the proof.
# E SOS-based Algorithms for LiMVAM

## E.1 Consistent estimator of the likelihood-based criterion

Assume we have observed $n$ independent samples, concatenated as $X = [x_1, ..., x_n]$ and $Y = [y_1, ..., y_n]$, with sample index $j$.

In practice, the four covariance matrices in Eq. 43 can be estimated by consistent estimators from the samples in $X$ and $Y$. Since $b_i = c_i = \text{cov}(x_i, y_i)$ for all $i$, we define

$$ \forall i: \quad \hat{b}_i := \hat{c}_i := \frac{1}{n} \sum_{j=1}^{n} x_{ij} y_{ij} \xrightarrow{p} \text{cov}(x_i, y_i) = b_i = c_i \quad (112) $$

so that $\hat{b} := (\hat{b}_1, ..., \hat{b}_m)^\top \xrightarrow{p} b$ and $\hat{c} := (\hat{c}_1, ..., \hat{c}_m)^\top \xrightarrow{p} c$. Using these vectors, we define the empirical covariances of the residuals

$$ \hat{\Sigma}_e := \frac{1}{n} \sum_{j=1}^{n} (y_j - \hat{b} \odot x_j)(y_j - \hat{b} \odot x_j)^\top \xrightarrow{p} \Sigma_e \quad (113) $$

$$ \hat{\Sigma}_d := \frac{1}{n} \sum_{j=1}^{n} (x_j - \hat{b} \odot y_j)(x_j - \hat{b} \odot y_j)^\top \xrightarrow{p} \Sigma_d \quad (114) $$

and we also define the empirical covariances of the observations

$$ \hat{\Sigma}_x := \frac{1}{n} \sum_{j=1}^{n} x_j x_j^\top \xrightarrow{p} \Sigma_x \quad (115) $$

$$ \hat{\Sigma}_y := \frac{1}{n} \sum_{j=1}^{n} y_j y_j^\top \xrightarrow{p} \Sigma_y . \quad (116) $$

We can then evaluate the criterion in these estimators.

$$ \widehat{\text{LR}} := \mathcal{L}(\hat{\Sigma}_x, \hat{\Sigma}_e; x \to y) - \mathcal{L}(\hat{\Sigma}_y, \hat{\Sigma}_d; y \to x) = -\log |\hat{\Sigma}_x| - \log |\hat{\Sigma}_e| + \log |\hat{\Sigma}_y| + \log |\hat{\Sigma}_d| \quad (117) $$

and it follows that this empirical criterion is a consistent estimator of the true criterion:

$$ \mathcal{L}(\hat{\Sigma}_x, \hat{\Sigma}_e; x \to y) - \mathcal{L}(\hat{\Sigma}_y, \hat{\Sigma}_d; y \to x) \xrightarrow{p} \mathcal{L}(\Sigma_x, \Sigma_e; x \to y) - \mathcal{L}(\Sigma_y, \Sigma_d; y \to x) =: \text{LR} . \quad (118) $$

## E.2 Consistent estimator of the cross-covariance-based criterion

Assume we have observed $n$ independent samples, concatenated as $X = [x_1, ..., x_n]$ and $Y = [y_1, ..., y_n]$, with sample index $j$.

In practice, the matrices $\mathbb{E}[ex^\top]$ and $\mathbb{E}[dy^\top]$ are estimated from the data. Performing univarite least-squares regression of $y$ on $x$ (resp. $x$ on $y$) yields the residual matrices $\hat{E} = [\hat{e}_1, ..., \hat{e}_n]$ (resp. $\hat{D} = [\hat{d}_1, ..., \hat{d}_n]$). The corresponding cross-covariance estimators are

$$ \frac{1}{n} \sum_{j=1}^{n} \hat{e}_j x_j^\top \xrightarrow{p} \mathbb{E}[ex^\top] \quad \text{and} \quad \frac{1}{n} \sum_{j=1}^{n} \hat{d}_j y_j^\top \xrightarrow{p} \mathbb{E}[dy^\top] . \quad (119) $$

Consequently, the empirical criterion

$$ \widehat{\text{FC}} := \left\| \frac{1}{n} \sum_{j=1}^{n} \hat{d}_j y_j^\top \right\|_F - \left\| \frac{1}{n} \sum_{j=1}^{n} \hat{e}_j x_j^\top \right\|_F \quad (120) $$

is a consistent estimator of the population criterion:

$$ \widehat{\text{FC}} \xrightarrow{p} \left\| \mathbb{E}[dy^\top] \right\|_F - \left\| \mathbb{E}[ex^\top] \right\|_F =: \text{FC} . \quad (121) $$
## E.3 Algorithms for estimating the two criteria

In this section, the observations $X, Y \in \mathbb{R}^{m \times n}$ lie in two dimensions: the $m$ views, and the $n$ independent samples. We use the index $i$ for the views, and $k$ for the samples.

**Algorithm 1** Likelihood-based criterion

**Input:** Observations $X, Y \in \mathbb{R}^{m \times n}$ corresponding to two variables, residuals $E, D \in \mathbb{R}^{m \times n}$ corresponding to the directions $x \to y$ and $y \to x$, respectively.

1. Compute the covariance matrices for the direction $x \to y$:

$$ \hat{\Sigma}_x = \frac{1}{n} \sum_{j=1}^{n} X_j X_j^\top \quad \text{and} \quad \hat{\Sigma}_e = \frac{1}{n} \sum_{j=1}^{n} E_j E_j^\top \quad (122) $$

and for the direction $y \to x$:

$$ \hat{\Sigma}_y = \frac{1}{n} \sum_{j=1}^{n} Y_j Y_j^\top \quad \text{and} \quad \hat{\Sigma}_d = \frac{1}{n} \sum_{j=1}^{n} D_j D_j^\top . \quad (123) $$

2. Compute the likelihood based-criterion:

$$ \widehat{LR} = - \log |\hat{\Sigma}_x| - \log |\hat{\Sigma}_e| + \log |\hat{\Sigma}_y| + \log |\hat{\Sigma}_d| . \quad (124) $$

**Output:** Likelihood-based criterion $\widehat{LR}$.

**Algorithm 2** Cross-covariance-based criterion

**Input:** Observations $X, Y \in \mathbb{R}^{m \times n}$ corresponding to two variables, residuals $E, D \in \mathbb{R}^{m \times n}$ corresponding to the directions $x \to y$ and $y \to x$, respectively.

1. Compute the second-moment matrices

$$ M_1 = \frac{1}{n} \sum_{j=1}^{n} E_j X_j^\top \quad \text{and} \quad M_2 = \frac{1}{n} \sum_{j=1}^{n} D_j Y_j^\top . \quad (125) $$

2. Compute the cross-covariance-based criterion:

$$ \widehat{FC} = \|M_2\|_F - \|M_1\|_F . \quad (126) $$

**Output:** Cross-covariance-based criterion $\widehat{FC}$.

## E.4 How to find the root variable

Assume we have a scalar criterion of the causal direction between two variables $x_j = (x_j^1, \dots, x_j^m)^\top$ and $x_k = (x_k^1, \dots, x_k^m)^\top$, such that the criterion is positive if $x_j \to x_k$ and negative in the opposite direction. Let $M_{jk}$ denote the criterion obtained for the pair $(j, k)$, and collect all pairwise criteria in a matrix $M$ with a zero diagonal. We now have a matrix $M$ whose entries determine the causal direction by their sign: $M_{ij} > 0$ indicates $x_i$ causes $x_j$, and $M_{ij} < 0$ indicates that $x_j$ causes $x_i$. In this setting, $x_j$ is a root variable if and only if the $j$-th row of $M$ contains only non-negative entries. We follow Hyvärinen and Smith [15, Section 3.2.2] to derive a principled criterion for finding the root variable given $M$. We next recall their argument.
Suppose that the entries in M follow a normal distribution, so that $M_{ij} \sim N(\mu_{ij}, \sigma^2)$ where $\sigma^2$ models the estimation error due to finite samples (and is supposed to be constant for simplicity). Then, the log-likelihood of $\boldsymbol{x}_i$ being the root variable is

$$ \log \prod_{j=1}^{p} P(\mu_{ij} > 0 | M_{ij}) = \sum_{j=1}^{p} \log P\left(\frac{\mu_{ij} - M_{ij}}{\sigma} > \frac{-M_{ij}}{\sigma} | M_{ij}\right) \quad (127) $$

$$ = \sum_{j=1}^{p} \log \Phi\left(\frac{M_{ij}}{\sigma}\right) \approx \frac{-1}{2\sigma^2} \sum_{j=1}^{p} \min(0, M_{ij})^2 \quad (128) $$

where $\Phi$ is the cumulative distribution function of a standardized Gaussian, whose log we then approximated. We therefore obtain the index of the root variable using

$$ \arg \min_{i} \sum_{j=1}^{p} \min(0, M_{ij})^2 . \quad (129) $$

## E.5 One-step Feasible Generalized Least Squares (FGLS)

The FGLS procedure begins with an Ordinary Least Squares (OLS) regression in each view to obtain residuals $(e'_j, \dots, e'^{m}_j)$, which are then used to estimate the cross-view covariance of the disturbances. In a second step, the OLS coefficients are adjusted via Generalized Least Squares using this estimated covariance. In the population limit, this estimator coincides with OLS but is asymptotically more efficient whenever disturbances are correlated across views, as is the case here.

The estimation proceeds as follows. Assume we have observed $n$ independent samples. Here, we consider the model

$$ \boldsymbol{x}'^i = \boldsymbol{T}^i \boldsymbol{x}'^i + \boldsymbol{e}'^i, \quad i \in [1, m] \quad (130) $$

where $\boldsymbol{x}'^i$ and $\boldsymbol{e}'^i$ are the reordered observations and disturbances, respectively, and $\boldsymbol{T}^i$ are strictly lower triangular matrices. By construction, each variable $x'_j$ only depends on its predecessors $x'_{<j} := (x'_{1}, \dots, x'_{j-1})$.

Fix a variable index $j \in [1, p]$. In view $i$, the structural equation reads

$$ x'_{j}^{i} = \sum_{k=1}^{j-1} T_{jk}^{i} x'_{k}^{i} + e'_{j}^{i} . \quad (131) $$

Stacking all views, we obtain the block regression model

$$ \boldsymbol{x}'_{j} = \boldsymbol{X}_{<j} \boldsymbol{T}_{<j} + \boldsymbol{e}'_{j} \quad (132) $$

where

* $\boldsymbol{x}'_j = (x'_j, \dots, x'^{m}_j)^{\top} \in \mathbb{R}^{mn}$ stacks all samples across the $m$ views,

* $\boldsymbol{X}_{<j} = \text{diag}(\boldsymbol{x}'_{<j}, \dots, \boldsymbol{x}'^{m}_{<j}) \in \mathbb{R}^{mn \times m(j-1)}$ is block-diagonal and contains the predecessors,

* $\boldsymbol{T}_{<j} = (\boldsymbol{T}_{j,<j}^{1}, \dots, \boldsymbol{T}_{j,<j}^{m})^{\top} \in \mathbb{R}^{m(j-1)}$ collects the coefficients,

* $\boldsymbol{e}'_j = (e'_j, \dots, e'^{m}_j)^{\top} \in \mathbb{R}^{mn}$ are disturbances.

The disturbance vector $\boldsymbol{e}'_j$ has covariance

$$ \Omega_j := \mathbb{E}[\boldsymbol{e}'_j(\boldsymbol{e}'_j)^{\top}] = \boldsymbol{\Sigma}_{e_j} \otimes \boldsymbol{I}_n \in \mathbb{R}^{mn \times mn} \quad (133) $$

where $\boldsymbol{\Sigma}_{e_j} \in \mathbb{R}^{m \times m}$ captures cross-view covariances for component $j$.

The one-step FGLS procedure is:
1. OLS step:
$$ \hat{T}_{<j}^{\text{OLS}} = (X_{<j}^{\top} X_{<j})^{-1} X_{<j}^{\top} x'_{j} \qquad (134) $$
yielding residuals $\hat{e}'_{j} = x'_{j} - X_{<j} \hat{T}_{<j}^{\text{OLS}}$.

2. Covariance estimation:
$$ \hat{\Sigma}_{e_j} = \frac{1}{n} \sum_{k=1}^{n} \hat{e}'_{j,k} (\hat{e}'_{j,k})^{\top} \qquad (135) $$

where k denotes the sample index.

3. FGLS update:
$$ \hat{T}_{<j}^{\text{FGLS}} = (X_{<j}^{\top} \hat{\Sigma}_{e_j}^{-1} X_{<j})^{-1} X_{<j}^{\top} \hat{\Sigma}_{e_j}^{-1} x'_{j} . \qquad (136) $$

This estimator reduces to OLS if $\Sigma_{e_j}$ is diagonal (no cross-view correlation). More generally, it achieves asymptotic efficiency by exploiting cross-view correlations.

## E.6 Full algorithm

**Algorithm 3 Pairwise-comparison algorithms for multi-view causal discovery**

**Input:** Observations $X \in \mathbb{R}^{m \times p \times n}$.

1. Preprocess the data $X$. Subtract the mean over the samples axis.

2. Estimate the causal ordering $P$.

(a) Perform pairwise regressions between each pair of variables, $x_j$ on $x_i$, such that $i \neq j$. Compute the residuals $e_{i \to j}$.

(b) Compute a skew-symmetric matrix $M \in \mathbb{R}^{p \times p}$ that determines the causal direction between pairs of variables.

The entries $M_{ij}$ are computed with Eq. 5 or Eq. 6 for all $i < j$, using the regression results $(x_i, x_j, e_{i \to j}, e_{j \to i})$. Then set $M_{ji} = -M_{ij}$. $M_{ij}$ positive means $x_i$ causes $x_j$, negative means $x_j$ causes $x_i$.

(c) Determine the root variable $k$: it causes all others, so $M_{kj}$ should be positive for all $j$.

$$ k = \underset{i}{\operatorname{argmin}} \sum_{j} \min(0, M_{ij})^2 . \qquad (137) $$

(d) Remove the root variable $k$ and replace observations with residuals obtained when regressing on $k$. Repeat.

(e) Store the estimated ordering in a permutation matrix $P$, and reorder the variables in $X$ according to this permutation. This yields a new model $x'^i = T^i x^i + e'^i$, where $x'^i$ and $e'^i$ are reordered versions of $x^i$ and $e^i$, respectively, and the $T^i$ are strictly lower triangular matrices.

3. Estimate the adjacency matrices $B^i$. For each variable $j \in [2, p]$, estimate the $j$-th row of all matrices $T^i$ with one-step Feasible Generalized Least Squares [37]. Recover matrices $B^i = P^\top T^i P$.

**Output:** Adjacency matrices $(B^1, \dots, B^m) \in \mathbb{R}^{m \times p \times p}$.
E.7 Computational complexity

The computational complexity of Algorithm 3 can be broken down across the different steps.

In Step 1, the preprocessing is in $O(m \cdot p \cdot n)$.

In Step 2.a., we solve $O(p^2)$ LS regression problems, each of which has complexity $O(m \cdot n)$.

In Step 2.b., we compute $O(p^2)$ entries of a matrix. Each entry is computed either using the cross-covariance-based criterion in Algorithm 2 which is in $O(m^2 \cdot n)$, or using the likelihood-based criterion in Algorithm 1 which is in $O(m^2 \cdot n + m^3)$ due to computing log determinants.

So the total complexity of this step is $O(m^2 \cdot p^2 \cdot n)$ or $O(m^2 \cdot p^2 \cdot n + m^3 \cdot p^2)$.

In Step 2.c., we make $O(p^2)$ operations.

In Step 2.d., we simplify redefine a matrix.

These steps 2.a-d., are repeated $O(p)$ times, until all variables are removed, so this yields a complexity in either $O(m^2 \cdot p^3 \cdot n)$ or $O(m^2 \cdot p^3 \cdot n + m^3 \cdot p^3)$.

In Step 2.e., we reorder a matrix, which costs $O(m \cdot p \cdot n)$ operations.

In Step 3, we run the one-step Feasible GLS procedure, which is in $O(m^3 \cdot p^3 \cdot n)$. It scales cubically in the $m$ views and $p$ dimensions as it involves matrix multiplications and inversions.

It then does $m$ matrix multiplications. Each of them involves permutation matrices which has a quadratic cost $O(p^2)$. So the cost of this step is dominated by $O(m^3 \cdot p^3 \cdot n)$.

The total computational complexity of the algorithm is thus $O(m^3 \cdot p^3 \cdot n)$.
# F Identifiability of LiMVAM with shared disturbances

In this section, we analyze the identifiability of the LiMVAM model with shared disturbances

$$x^i = B^i x^i + D^i s + n^i, \quad i \in [1, m] \qquad (138)$$

where the $D^i$ are diagonal matrices with positive entries on the diagonal, $s$ has mutually independent entries and second-order moment $E[ss^\top] = I_p$, the view-specific noises are Gaussian $n^i \sim N(0, \Sigma^i)$ with diagonal $\Sigma^i$, and the vectors $s$ and $n^i, i = 1, \dots, m$, are mutually independent.

## F.1 Proof of Theorem 2 (first claim)

Consider two sets $\Theta = (D^1, \dots, D^m, \Sigma^1, \dots, \Sigma^m, B^1, \dots, B^m)$ and $\Theta' = (D'^1, \dots, D'^m, \Sigma'^1, \dots, \Sigma'^m, B'^1, \dots, B'^m)$ that parameterize the same statistical model in Eq. 138, with $D^i$ and $D'^i$ being diagonal matrices with positive entries on the diagonal, $\Sigma^i$ and $\Sigma'^i$ being diagonal matrices with non-zero entries on the diagonal, and $B^i$ and $B'^i$ being DAG matrices. The model in Eq. 138 can be reformulated as

$$x^i = (I - B^i)^{-1} D^i (s + (D^i)^{-1} n^i) \qquad (139)$$

which corresponds to the Shared ICA model

$$x^i = A^i(s + \tilde{n}^i) \qquad (140)$$

for $A^i = (I - B^i)^{-1} D^i$ and $\tilde{n}^i = (D^i)^{-1} n^i$. We can observe that the unmixing matrices $W^i = (A^i)^{-1} = (D^i)^{-1}(I - B^i)$ and $W'^i = (A'^i)^{-1} = (D'^i)^{-1}(I - B'^i)$ belong to the domain $\mathcal{W}$, which is included in the space of invertible matrices. Thus, the two sets of $W^i$ and $W'^i$ matrices are valid sets of unmixing matrices for the same multi-view ICA model in Eq. 140. Moreover, Assumption 4 allows the rescaled noises $\tilde{n}^i$ to meet the noise diversity condition of Shared ICA. So, from the identifiability theory of multi-view ICA [25, Theorem 1], we know that there exist a sign-permutation matrix $Q$ such that for any view $i \in [1, m]$, we have

$$W'^i = Q^\top W^i \\ \Sigma'^i = Q^\top \Sigma^i Q. \qquad (141)$$

Note that, contrary to the single-view context, the fact that we obtain $W$ up to sign and permutation rather than scale and permutation is because $E[ss^\top] = I$. Then, we apply Lemma 3, which shows that being in the domain $\mathcal{W}$ imposes $Q = I$, $D'^i = D^i$, and $B'^i = B^i$.

## F.2 Proof of Theorem 2 (second claim)

In the context of shared causal ordering $P$, we prove that, under Assumptions 4 and 5, the decomposition of matrices $B^i$ into matrices $T^i$ and $P$ is unique.

Consider the two sets $\Theta = (\Sigma^1, \dots, \Sigma^m, P, T^1, \dots, T^m)$ and $\Theta' = (\bar{\Sigma}^1, \dots, \bar{\Sigma}^m, \bar{P}, \bar{T}^1, \dots, \bar{T}^m)$ that parameterize the same statistical model in Eq. 138, where $P$ is some permutation matrix, $T^i$ are strictly lower triangular matrices, and $\bar{P}$ and $\bar{T}^i$ are the particular matrices given by Assumption 5. Note that no assumption on the sparsity of the $T^i$ is made. From Theorem 2, we know that $\Sigma^i = \bar{\Sigma}^i$.

Assumption 5 amounts to saying that the “reunion” of the $\bar{T}^i$ represents a fully connected graph. In other words, if we define the reunion as $\bar{T}^U = \sum_{i=1}^m \text{abs}(\bar{T}^i)$, where abs denotes the element-wise absolute value function, then we are assuming that the strictly lower triangular part of $\bar{T}^U$ only has non-zero elements. In a similar way to the definition of $\bar{T}^U$, let us define

$$\bar{W}^U = \sum_{i=1}^{m} \text{abs}(\bar{W}^i) \qquad (142)$$
where $\bar{W}^i = I - \bar{P}^\top \bar{T}^i \bar{P}$. The non-zero elements of $\bar{P}^\top \bar{T}^i \bar{P}$ are outside of the diagonal, so matrices $I$ and $\bar{P}^\top \bar{T}^i \bar{P}$ contain non-zero elements at different locations. Thus, we have

$$ \text{abs}(I - \bar{P}^\top \bar{T}^i \bar{P}) = I + \text{abs}(\bar{P}^\top \bar{T}^i \bar{P}) . \quad (143) $$

Furthermore, applying $\bar{P}$ to the rows and columns of $\bar{T}^i$ only shuffles its entries, without modifying their values. So we have

$$ \text{abs}(\bar{P}^\top \bar{T}^i \bar{P}) = \bar{P}^\top \text{abs}(\bar{T}^i) \bar{P} . \quad (144) $$

Consequently,

$$ \bar{W}^U = \sum_{i=1}^{m} I + \bar{P}^\top \text{abs}(\bar{T}^i) \bar{P} = \sum_{i=1}^{m} \bar{P}^\top (I + \text{abs}(\bar{T}^i)) \bar{P} = \bar{P}^\top (mI + \bar{T}^U) \bar{P} . \quad (145) $$

So, dividing both sides by $m$ gives

$$ \frac{1}{m} \bar{W}^U = \bar{P}^\top \left( I - \left( -\frac{1}{m} \bar{T}^U \right) \right) \bar{P} \quad (146) $$

where $-\frac{1}{m}\bar{T}^U$ is strictly lower triangular, and its strictly lower triangular part only has non-zero elements.

Next, we apply the same reasoning to the alternative set of parameters, given by $W^i = I - P^\top T^i P$, and we consider $W^U = \sum_{i=1}^m \text{abs}(W^i)$. The proof of Theorem 2 already implied that $\bar{W}^i = W^i$ for all $i$. Thus, we have $\bar{W}^U = W^U$ and

$$ \bar{P}^\top \left( I - \left( -\frac{1}{m} \bar{T}^U \right) \right) \bar{P} = P^\top \left( I - \left( -\frac{1}{m} T^U \right) \right) P \quad (147) $$

where $T^U$ is defined in a similar way as $\bar{T}^U$, except that its strictly lower triangular part can be sparse. Using Lemma 4 on Eq. 147, we obtain that $P = \bar{P}$, and thus $T^i = \bar{T}^i$ for all $i$. In conclusion, all the sets of DAG decompositions that parameterize the same model are equal, as soon as for one of these decompositions, the reunion of the $T^i$ is dense. We conclude that matrices $P$ and $T^i$ are unique, and thus identifiable in our terminology. This proof also implies that there can be only one matrix $\bar{P}$ that fulfills Assumption 5.

## F.3 Results for view-specific causal orderings

Next, we consider a case which is outside of the theory of the main paper, although a simple extension: we allow the causal orderings to be different in different views. It turns out that in the view-specific $P^i$ case, identifiability is obtained by assuming that the directed acyclic graph $B^i$ is dense enough in each view, as formalized in the following assumption. However, since the $B^i$ can be permuted to strictly lower triangular matrices, they cannot be denser than having $\frac{p(p-1)}{2}$ non-zero entries.

**Assumption 7** (Dense connectivity in each view) For each view $i$, the matrix $B^i$ has exactly $\frac{p(p-1)}{2}$ non-zero entries.

Using Assumption 7, the following theorem states that, in addition to identifying $B^i$, $D^i$, and $\Sigma^i$, one can also identify $T^i$ and $P^i$.

**Theorem 8** (Identifiability of multiple causal orderings) Consider the LiMVAM model with multiple causal orderings. Consider the quantities $(D^1, \dots, D^m, \Sigma^1, \dots, \Sigma^m, P^1, \dots, P^m, T^1, \dots, T^m)$ as the set of parameters to be estimated. Under Assumptions 4 and 7, all these parameters are identifiable.
In the context of view-specific causal orderings $P^i$, we prove that, under Assumptions 4 and 7, the decomposition of matrices $B^i$ into matrices $T^i$ and $P^i$ is unique.

Consider two sets of parameters $\Theta = (\Sigma^1, \dots, \Sigma^m, P^1, \dots, P^m, T^1, \dots, T^m)$ and $\Theta' = (\Sigma'^1, \dots, \Sigma'^m, P'^1, \dots, P'^m, T'^1, \dots, T'^m)$ that parameterize the same statistical model in Eq. 138, where $P^i, P'^i$ are permutation matrices, and $T^i, T'^i$ are strictly lower triangular matrices. Note that here we parameterize by $P^i$ and $T^i$ rather than $B^i$ or $W^i$. From Theorem 2, we know that $\Sigma^i = \Sigma'^i$ and that the resulting causal matrices must be equal:

$$ (P^i)^\top T^i P^i = (P'^i)^\top T'^i P'^i. \qquad (148) $$

Assumption 7 states that, for each view $i$, the matrix $B^i = (P^i)^\top T^i P^i = (P'^i)^\top T'^i P'^i$ contains exactly $\frac{p(p-1)}{2}$ non-zero elements, so it is also the case for $T^i$ and $T'^i$ which thus represent fully connected graphs. So, from Lemma 4, we deduce that, for each view $i$, we have $P^i = P'^i$ and $T^i = T'^i$, which concludes the proof.
# G ICA-based Algorithm for LiMVAM with shared disturbances

## G.1 ICA-based Algorithm

**Algorithm 4 ICA-LiMVAM**

1. Estimate the adjacency matrices $B^i$.

(a) Estimate the unmixing matrices $W^i$ and the noise variance matrices $\Sigma^i$, by running the Shared ICA estimation algorithm on the data.

The algorithm returns an estimate of the true unmixing matrix $W^i = M(I - B^i)$ but up to sign-permutation matrix $M$ that is the same for all views [25].

(b) Determine the sign-permutation indeterminacy $M$, using the structure of the underlying $(A^i)^{-1}$ that has a diagonal of ones. Thus, we can adapt the simple two-step heuristic by Shimizu et al. [30].

First, find a permutation matrix $M$ such that $MW$ has a non-zero diagonal, by solving with the Hungarian algorithm

$$ \min_M \sum_{j=1}^{p} \frac{1}{|(MW)_{jj}|}, \quad W = \sum_{i=1}^{m} \text{abs}(W^i), \qquad (149) $$

and then rescale the rows of $M$ to ensure that $MW$ has a diagonal of ones

$$ M_{ij} \leftarrow \text{sign}((MW)_{ii})M_{ij}, \qquad (150) $$

where the sign function outputs 1 or -1.

(c) Determine the scale matrices $D^i = \text{diag}(\frac{1}{(MW^i)_{11}}, \dots, \frac{1}{(MW^i)_{pp}})$.

(d) Determine the causal matrices $B^i = I - D^i MW^i$ and update the noise variance matrices with $\Sigma^i \leftarrow (D^i)^2 \Sigma^i$.

2. Determine the causal ordering $P$.

The DAG decomposition states that $B^i := P^\top T^i P$ where $T^i$ is lower triangular, so we penalize for that

$$ \min_P \sum_{l \ge k} (PBP^\top)_{kl}^2, \quad B = \sum_{i=1}^{m} \text{abs}(B^i). \qquad (151) $$

We approximately minimize this objective using the heuristic algorithm presented in Shimizu et al. [30, Algorithm C].

Note that the algorithm described in Chen et al. [5] finds the same adjacency matrices as we do. This is because of the sign-permutation matrix $M$ from ICA: to determine that matrix, we ensure $MW$ has a diagonal of ones by dividing by the scalings $D^i$. For [5], these scalings are not part of the model: they are simply a part of the algorithm which determines $M$. For us, these scalings are part of our model.
## G.2 Computational complexity

We next detail the worst-case computational complexity of the ICA-LiMVAM algorithm which is in $O(tnmp^3 + mp^5)$, for $t$ iterations of the SharedICA-ML algorithm, $n$ samples, $m$ views and $p$ components. Note that the term in $p^5$ comes from Algorithm C of the original LiNGAM and is worst-case: in practice, it took a couple of iterations only in our experiments. Below are the detail of the derivations for each step of Algorithm 4.

**Computational complexity** The computational complexity of ICA-LiMVAM can be broken down across the four steps in Algorithm 4. Step 1 calls SharedICA, with its most costly variant (SharedICA-ML) running in $O(t \cdot n \cdot m \cdot p^3)$, where $t$ is the number of iterations, $n$ the number of samples, $m$ the number of views, and $p$ the number of components. Step 2 solves the permutation indeterminacy via a linear sum assignment in $O(m \cdot p^3)$, instead of trying $p!$ permutations. Step 3 involves $m$ multiplications with permutation matrices and additions, costing $O.mp^2)$. Step 4 optimizes over permutations. Instead of trying $p!$ permutations, we use the heuristic Algorithm C from LiNGAM. It is $O(m \cdot p^5)$ in the worst case (and much lower in practice). Further details are provided next:

1. Run the ICA algorithm SharedICA-ML, $O(tnmp^3)$

Each of the $t$ iterations of the optimization performs the so-called “E-step” and “M-step” of the E-M algorithm, detailed in Section 4 of Richard et al. [25]. Their complexity is in $mnp^3$: the intuition is that for each view $m$ and each sample $n$ require some matrix operations (e.g. addition, multiplication, computing gradients and Hessians) all of which are dominated by $p^3$. So the final complexity is in $O(tnmp^3)$.

2. Solve the permutation-sign indeterminacy, $O.mp^3)$

Computing the objective can be done in $O(p^2)$. The detail of this is: we begin by computing the matrix $W$ which is in $O.mp^2)$, then multiplying it with a permutation matrix which is in $O(p^2)$, and finally summing the diagonal elements which is in $O(p)$.

Finding the correct permutation matrix can be done by solving a linear-sum assignment problem; this can be done using the Hungarian algorithm which is in $O(p^3)$.

3. Update the entries of a matrix which is in $O(p^2)$.

4. Update the entries of a matrix which is in $O(p^2)$.

Adding the complexities thus far leads to a complexity of $O.mp^2 + p^3)$, which is dominated by $O.mp^3)$, for Step 2.

5. Compute the causal matrices, $O.mp^2)$

Each of the $m$ causal matrices requires multiplying a dense matrix with a permutation which is in $O(p^2)$ and then adding a matrix which is in $O(p^2)$. The final cost is $O.mp^2)$.

6. Find the causal ordering(s), $O.mp^5)$

Algorithm B in the original LiNGAM paper finds a row that has all zeros $O(p^2)$, removes it, does this $p$ times. Its complexity is in $O(p^3)$.

Algorithm C in the original LiNGAM paper calls Algorithm B at most $\frac{p(p-1)}{2} = O(p^2)$ times, so its final complexity is $O(p^2p^3) = O(p^5)$.

When the causal ordering is shared across views: first we compute sum $m$ matrices with $p^2$ entries which is in $O.mp^2)$, and then we use Algorithm C once which is in $O(p^5)$. The final complexity is in $O.mp^2 + p^5)$ which is dominated by $O.mp^5)$.

When the causal ordering is view-dependent: for each of the $m$ views, we use Algorithm C which is in $O(p^5)$, so the final complexity is in $O.mp^5)$.
# H Experiments

## H.1 Synthetic experiments

### H.1.1 Data generation in Figure 1

The causal effect matrices $B^i = P^\top T^i P$ are generated from a random permutation $P$ and strictly lower triangular $T^i$ obtained from a standard Gaussian; the diagonal matrices $D^i$ are drawn from a uniform density between 0.1 and 3; the common disturbances in $s$ can be either Gaussian or non-Gaussian: Gaussian disturbances $s_j$ are generated from a standardized Gaussian and their corresponding noises $n_j^i$ have standard deviation $\Sigma_{jj}^i$ obtained by sampling from a uniform density between 0 and 1, while non-Gaussian disturbances are generated from a Laplace distribution (with scale parameter equal to $\frac{1}{2}$) and their corresponding noises all have a std of $\frac{1}{2}$. We use $m = 5$ views, $p = 4$ variables, and vary the number of samples $n$ between $10^2$ and $10^4$.

### H.1.2 Simulation study in higher dimension

We evaluate the methods by measuring the estimation error on the matrices $B^i$ across a range of settings involving more views and components than those considered in the main text. Specifically, we vary the number of views among $\{3, 5, 8, 12, 16, 20\}$ and the number of components among $\{3, 6, 9, 12\}$, using 1000 samples and 50 random seeds for each configuration.

<figure><img src="image_7.png" /><figcaption>Figure 5: Separation performance of three pairwise-comparisons-based algorithms — Multi-GroupDirectLiNGAM from Shimizu [28], and DirectLiMVAM and PairwiseLiMVAM which use the criteria in Eq. 5 and Eq. 6, respectively — and three ICA-based algorithms — ICA-LiNGAM, which naively uses Shimizu et al. [30] separately in each view, and our implementation of the two variants of ICA-LiMVAM [5]. We varied the number of views and components. Disturbances are generated using the generalized normal distribution with shape parameter $\beta \in \{1.5, 2, 2.5\}$. The metric used is the $\ell_2$-distance between true and estimated causal effect matrices $B^i$ (lower is better).</figcaption></figure>

Data are generated from the model in Eq. 3, where disturbances $e_j$ are split into three types—sub-Gaussian, Gaussian, and super-Gaussian—which motivates the choice of disturbance counts as multiples of 3. Sub-Gaussian disturbances follow a generalized normal distribution with shape parameter $\beta = 2.5$; super-Gaussian disturbances use $\beta = 1.5$; Gaussian disturbances correspond to $\beta = 2$. Each $T^i$ is sampled as a strictly lower-triangular matrix with Gaussian entries, and a common permutation matrix $P$ is used to construct the $B^i = P^\top T^i P$. For each method, we
report the median error on the matrices $B^i$, with the 25th and 75th percentiles shown as error bars.

As shown in Fig. 5, DirectLiMVAM and PairwiseLiMVAM consistently achieve lower errors than the other approaches. DirectLiMVAM performs best when the number of views is limited (3, 5, or 8) while the number of components remains relatively high (9 or 12), thereby outperforming PairwiseLiMVAM in these settings. This slightly surprising observation that DirectLiNGAM is a bit better than PairwiseLiNGAM is perhaps due to the fact that there the particular non-Gaussianities used here are ill-suited for the Gaussian likelihood underlying PairwiseLiNGAM. In contrast, ICA-LiNGAM and MultiGroupDirectLiNGAM achieve slightly better-than-average performance when both the number of views and the number of components are small, but their performance quickly degrades to the average level as dimensionality increases. This behavior reflects their inability to properly handle Gaussian disturbances. Finally, both ICA-LiMVAM variants fail completely, as they are not designed to operate without shared disturbances. Overall, the estimation error tends to increase with the number of components.

### H.1.3 Testing Assumption 4

Assumption 4 provides a sufficient condition for the identifiability of the causal matrices $B^i$ in the LiMVAM model with shared, Gaussian disturbances. To evaluate the practical impact of this assumption, we simulated data from the model in Eq. 3, with shared disturbances as in Eq. 9. We used 5 views, 4 common disturbances in $s$ (2 Gaussian and 2 Laplacian), and 1000 samples.

The scale matrices $D^i$ were drawn uniformly from the interval $[0.5, 2]$, and the strictly lower triangular matrices $T^i$ were sampled from a Gaussian distribution and then permuted to form the causal matrices $B^i$. The noise variances $\Sigma_{jj}^i$ corresponding to the Laplacian disturbances were fixed to $\frac{1}{2}$, while the noise variances corresponding to the two Gaussian disturbances, indexed by $j$ and $j'$, were initially sampled uniformly in $[0, 1]$. To test the assumption, we selected an increasing number of views $i$ such that the scaled noise variances

$$ \frac{\Sigma_{jj}^i}{D_{jj}^i} = \frac{\Sigma_{j'j'}^i}{D_{j'j'}^i}, \text{ ranging from 0 to 5. Note that} $$

Assumption 4 only fails when this condition holds across all 5 views. We repeated this experiment over 50 random seeds.

Figure 6 shows that ICA-LiMVAM maintains low estimation error on the $B^i$ matrices as long as the noise diversity assumption holds, but its performance deteriorates sharply once the assumption is violated. This confirms the practical relevance of the assumption. Interestingly, PairwiseLiMVAM appears unaffected by the violation of this condition.

<figure><img src="image_8.png" /><figcaption>Figure 6: Effect of violating the noise diversity assumption (Assumption 4) on the estimation error of the causal matrices $B^i$. We report the average $\ell_2$ error over 50 repetitions for ICA-LiMVAM and PairwiseLiMVAM, as the number of views with identical scaled noise variances increases.</figcaption></figure>

### H.1.4 Testing Assumptions 5 and 7

Assumption 5 (resp. Assumption 7) gives a sufficient condition for the identifiability of the causal matrix $P$ (resp. causal orderings $P^i$) and matrices $T^i$, in the shared causal ordering (resp. multiple causal orderings) case. However, it is not straightforward how robust the algorithms are to these assumptions.
<figure><img src="image_9.png" /><figcaption>Figure 7: Spearman’s rank correlation between true and estimated causal orderings as a function of the number of sparse entries in each $T^i$. Left: In the multiple-ordering setting, performance of ICA-LiMVAM-ML and ICA-LiMVAM-J degrades linearly with increasing sparsity, as the assumption of full support is violated. PairwiseLiMVAM is excluded since it assumes a shared ordering. Right: In the shared-ordering setting, all methods maintain near-perfect recovery up to 11 sparse entries, after which performance drops sharply.</figcaption></figure>

To evaluate the practical impact of these two assumptions, we simulated data from the model in Eq. 3, with shared disturbances as in Eq. 9. We used 8 views, 6 common disturbances (two sub-Gaussians sampled from a generalized normal distribution with shape parameter $\beta = 2.5$; two super-Gaussians for $\beta = 1.5$; two Gaussians corresponding to $\beta = 2$), and 1000 samples. The scale matrices $D^i$ were drawn uniformly from the interval $[0.5, 2]$, the strictly lower triangular matrices $T^i$ were sampled from a Gaussian distribution and then permuted to form the causal matrices $B^i$, and the variances $\Sigma^i$ were sampled uniformly in $[0, 1]$. The experiment was repeated 30 times.

Fig. 7 reports the Spearman’s rank correlation between the true and estimated causal orderings as a function of the number of sparse entries in each $T^i$, randomly chosen among the $6 \times 5 = 30$ non-zero entries.

In the multiple causal orderings setting (left panel), both ICA-LiMVAM variants exhibit a roughly linear drop in performance as sparsity increases. This behavior is expected, as Assumption 7 no longer holds once any entry in the strictly lower triangular part of $T^i$ is set to zero. PairwiseLiMVAM is not included in this setting, as it is designed under the assumption of a shared causal ordering.

In the shared causal ordering setting (right panel), all methods benefit from the shared structure and recover the true ordering $P$ reliably up to 11 sparse entries. Beyond this point, performance degrades sharply for all methods, reflecting the increasing violation of Assumption 5. These results empirically support the practical relevance of Assumptions 5 and 7.

## H.2 Real data experiments

### H.2.1 Details on the preprocessing of MEG data

The MEG data measures participants’ responses to auditory and visual stimuli. The auditory stimuli were binaural pure tones. The visual stimuli were checkerboards presented both to the left and right of a central fixation for 34-ms duration. This task leads to strong signal power modulations during the motor preparation and motor execution.

The original MEG data were acquired with 306 sensors, recorded at 1000 Hz, and band-pass filtered between 0.03 and 330 Hz. All MEG processing was done using the MNE-Python library [9, 10] and we largely followed the pre-processing steps used in Power et al. [23]. We applied a Maxwell filter [35] to improve data quality and a band-pass filter between 8 and 27 Hz to focus on power effects spread over the alpha and beta bands of the brain. This range of waves is supposed to be particularly active in sensorimotor tasks, especially during movement
preparation and execution, and it typically shows a characteristic suppression (event-related desynchronization) during movement, followed by a rebound (event-related synchronization) after movement cessation, which is thought to reflect sensorimotor processing and inhibitory mechanisms. In particular, the suppression and rebound are reflected in the energies of the signals, not raw signals.

The data were then parsed into trials synchronized to each button press, with a duration of 4.5 s, including a 1.5 s pre-movement interval. The 4.5 s window length was selected to ensure a sufficient post-movement interval to capture the entire beta rebound response. Trials were excluded if the button press occurred more than 1 s after the audiovisual cue (indicating poor task performance) or if another button press occurred within the time window. Then, a baseline correction was applied using the pre-movement interval (–1.5 s, –1 s). The procedure led to about 60 trials per participant on average.

<figure><img src="image_10.png" /><figcaption>Figure 8: Example of time series obtained after the preprocessing. We averaged the data across subjects and trials, resulting in 10 time series, one for each brain region.</figcaption></figure>

Next, we performed a cortical projection. First, each participant's MRI (additionally given in the Cam-CAN database) was segmented using FreeSurfer [7]. The segmentation provided a digitization of the cortical surface for source estimation, a transformation to the average brain (i.e. fsaverage) for spatial normalization and group statistics, and a boundary element model of the head to provide more accurate calculation of the forward solution. The inverse solution, based on the MNE method [12], allowed to consider cortical region activations for further analysis.

We used the cortical parcellation from [16] to divide the cortical man-
tel (both hemispheres) into 448 distinct regions and summarized each
region by an averaged time course. Then, we selected 10 of the 448 re-
gions based on their known importance in sensorimotor tasks. Specifi-
cally, we picked for each hemisphere three regions in the motor cortex
(two parcels in the “postcentral” and one in the “precentral”; visible
in blue in Fig. 3a), one region in the auditory cortex (“superiortem-
poral” parcel; highlighted in pink), and one region in the visual cortex
(“pericalcarine” parcel; not visible in the figure). Note again that
the task for the participant is active as right index button presses are
triggered by audiovisual stimulations.

For each participant, we thus extracted one time series for each of
the 10 regions and fixed the number of trials to 40. Participants
with less than 40 trials were discarded and extra trials were averaged.
Participants for whom some of the parcels did not have any vertex in
the source space were also discarded, resulting in a total of 98 available
participants. We performed a Z-score normalization of the time series

to correct the depth bias and, importantly, computed the Hilbert envelope of the signals to study
modulations in cortical source power. Finally, the signals were centered and downsampled from
1000 Hz to 10 Hz due to the slowly changing nature of the envelopes (energies). The final dataset
consisted of 98 subjects, 10 brain regions, and 1760 time points.

<figure><img src="image_11.png" /><figcaption>Figure 9: Pearson correlations between median causal effect matrices across 50 runs of ICA-LiMVAM-ML, each using a different random subset of 30 participants from the Cam-CAN cohort.</figcaption></figure>
Figure 8 shows the typical time series obtained after preprocessing. For visualization, we averaged the time courses across participants and trials. As expected, we observe a prominent peak in the “superiortemporal” auditory regions shortly before the button press (t = 0 s), reflecting the stimulus onset. In the motor regions, a characteristic rebound pattern emerges: the signals decrease around the time of the button press, consistent with event-related desynchronization of beta rhythms, and subsequently increase after approximately 0.3 seconds, indicating beta resynchronization. This figure also labels the ten cortical parcels selected for analysis.

## H.2.2 MEG experiment using ICA-LiMVAM-ML

In the following, we present additional analyses that extend the experiments described in Section 6.2.

Figure 10 displays median causal effect maps estimated using ICA-LiMVAM-ML on the CamCAN dataset. For each of the six panels, ICA-LiMVAM-ML was applied to a randomly chosen subset of 50% of the participants. These results are consistent with the patterns observed in Figure 3a, notably the frequent presence of directed connections between motor areas in opposite hemispheres. This further supports the hypothesis of consistent inter-hemispheric causal influences in sensorimotor processing.

<figure><img src="image_12.png" /><figcaption>Figure 10: Top ten strongest median causal effects (estimated by ICA-LiMVAM-ML) of six different runs. Each run was performed on the data of 49 randomly chosen subjects. Red arrows represent positive effects and blue arrows negative effects.</figcaption></figure>

Finally, we conducted the same robustness analysis for ICA-LiMVAM-ML as previously done for PairwiseLiMVAM (Figure 3b). Specifically, we ran ICA-LiMVAM-ML 50 times, each time using a randomly selected subset of 30 participants from the full cohort, and computed the Pearson correlations between the element-wise median of the estimated individual matrices $B^i$ across runs. The resulting correlations are shown in Figure 9. While the correlations remain predominantly positive (with an average of 0.27), they are noticeably lower than those obtained with PairwiseLiMVAM. This suggests that, in terms of stability across subsets, PairwiseLiMVAM exhibits greater consistency than ICA-LiMVAM-ML.