On the Complexity of Identification in Linear Structural Causal Models

Answer to questions:

1. At the moment, we consider the main contribution of the paper on the theoretical side. We are currently working on implementations.

2. We would like to emphasize that since the formulation of the parameter identification problem in linear SEMs in the 1960s, this problem has been intensively studied by econometricians and sociologists, and more recently by the ML, AI and statistics communities. The key importance of parameter identification (from observed data) is that the inferred parameters can be used to estimate causal effects, especially in cases where randomized controlled trials are not possible.

However, despite decades of extensive research, no (complete) algorithm has been found that matches the general applicability of the Gröbner basis approach and thus runs faster than in doubly-exponential-time. Moreover, no non-trivial lower bound (hardness result) on the computational complexity of this problem has been established which would justify the high computational complexity of state-of-the-art algorithms. In our work, we give a novel algorithm that provides a constructive upper bound on the computational complexity of the generic identification problem, which significantly improves on the best-known doubly-exponential bound. Moreover, we give the first hardness results for numerical parameter identification in linear SEMs showing that it is hard for the complexity class ∀R.

The proof technique provided in our paper is novel. In hindsight, the proofs in Sections 5 and 6 look easy, but no one made this connection to the existential theory of the reals as well as the result by Koiran before. To the best of our knowledge, this proof technique has not been used before for linear causal model analysis, which would be of interest to the NeurIPS community. Furthermore, the reduction in Section 4 is quite astonishing, since already SCMs with only two layers turn out to be hard.

Regarding your recommendation to revise Sections 4, 5, and 6:

We think that our proofs are quite clear. However, when writing for diverse communities, it is always challenging to balance between intuition and mathematical rigor. We thank you for your advice and will follow it by moving more technical parts to the appendix and adding more intuition.

• Recursive (or, in graphical language, acyclic) models are fairly standard and most commonly used in causality. They are the basic structural causal models used in causal inference (e.g. in the famous Pearl's do-calculus, Markov-equivalence of causal models, etc.) and in causal structure learning (e.g., the well-known PC algorithm by Spirtes and Glymour or GES algorithm by Chickering learn from observed data acyclic models). Also (causal) Bayesian models are represented by such acyclic graphs. Furthermore, when restricting general causal models to linear ones (which are the subject of our work), recursive models are most commonly used. In particular the papers studying parameter identification in linear SEMs cited in our work [8, 10, 11, 14, 15, 24, 28, 29, 31, 36, 37, 38, 40, 41] assume recursive causal models. Since some recent works try to consider cyclic models as well, we have indicated in line 35 that our work assumes the standard model. Nevertheless, as we observed very recently, this restriction is not needed and our results can be generalized to non-recursive models. We will add an explanation in the final version.

• Identification in linear structural causal models naturally reduces via Equation (10) in the paper to solving polynomial systems of equations over the reals. Therefore, tools from commutative algebra and real algebraic geometry are inherently needed. We agree that we should have given short inline descriptions of the terms you mention. We will do so and also add a more detailed appendix on the basic algebraic concepts used.

• We hoped that it would be easier and more accessible if we develop our constructions in a textual manner. We find this quite hard to balance, since the audience of the paper is quite diverse. Other reviewers wanted more informal explanations and less proofs. Note that the proof of Theorem 2 for instance is still a rigorous proof, it is a so-called "gadget-based" proof, which are used commonly. The two SCMs in Figure 4 are not just examples. They are our basic building blocks. They can be instantiated in such a way that they implement the equations given in (4), the basic building blocks of a QUAD^{++} instance. The gadget on the right-hand side implements the equation ab - c = 0, and the other gadget all other four equations with different instantiations of the $\sigma$-values. This is described right after Figure 4. Then one simply takes these basic building blocks and combines them as described in the proof of Theorem 2 to encode the whole QUAD^{++} instance. We will expand the explanation of the proof structure accordingly.

• Large parts of our work are conceptual. The complexity of generic identification has been an open problem since the 1960s. In hindsight, the proof looks easy, but no one made this connection to the existential theory of the reals as well as the result by Koiran before. Moreover, no hardness results for identification were known before.

• Thank you for pointing out the grammatical errors/typos, we will correct them.

Regarding weaknesses:

• We agree that learning the graph structure of the causal model and/or the topological ordering of variables is a difficult and very challenging task. However, we would like to emphasize that learning causal structures and/or the topological orderings was not the subject of this work and, as it is commonly accepted in studies on the parameter identification in linear SEMs, we assume that such a structure is given as an instance of the problem. It is also common to assume that noises are normal. We would like to point out that while these assumptions make the upper bounds easier in some sense, they make hardness results more difficult to prove. So they only work partly in our favor.

• We would like to explain that our algorithm (in the proof of Theorem 11) does not need to have access to an oracle. To make the description of the algorithm modular and easy to understand, we first provide an algorithm having oracle access to ETR and next we explain that using Renegar’s algorithm the ETR queries can be computed in PSPACE (Corollary 3). Therefore, we only use them to make the statements more convenient. You can think of it as procedures that we call. The oracle calls can always replaced by actual algorithms. We will point this out more clearly.

• We thank the reviewer for commenting on the need to add some supporting explanations to the designations.

Regarding Questions:

• Being positive definite has some implications on the diagonal, but depending on the definition, it is not always stated explicitly. For instance, in the first definition used in the proof of Lemma 5, there is some explicit restriction stated, while the second definition in the proof does not have such an explicit restriction. Since both definitions are equivalent, there is an implicit restriction in the second definition. Being diagonally dominant implies positive definiteness but the other direction is not true in general.

• We apologize for the confusion and thank you for pointing this out. A few indices are not correct, see below for corrections. (We made some last minute changes and should not have done this.)

  -- On the left-hand side in Figure 3, the r's should be i and j, respectively.

  -- In the sum in the middle of (11), $\lambda_{a_\ell,i}$ should be $\lambda_{b_k,j}$. We are sorry for the typo. The equation is obtained from (10) by simply writing down the entries of the matrix product and looking at entry (i,j).

  -- Obseration 4 follows from (10) and the fact the our graph has only two layers. Since i and j are both in the top layer, $\sigma_{i,j}$ can only appear in the equation in (10) that corresponds to entry (i,j) in the matrix. It can never appear multiplied with any $\lambda_{k,\ell}$, because then a least one of i or j needs to be in the bottom layer. The value $\sigma_{a_\ell,j}$ can only appear in an equation corresponding to a missing edge $\{h,j\}$ that contains j and when there is a directed edge $(a_\ell,h)$. But there is only one leaving $a_\ell$ by construction and this edge ends in i, hence $h = i$. The same is true for $\sigma_{b_k,i}$. For $\sigma_{a_\ell,b_k}$ to appear in an equation given by a missing edge $\{g,h\}$, there have to be directed edges $(a_\ell,g)$ and $(b_k,h)$. Again by construction, $g = i$ and $h = j$.

• The algorithm itself starts in line 320. But it calls subroutines that are described earlier. In step 1, we obviously call Koiran's algorithm. In step 2, we call the algorithm from Lemma 5 to decide membership in PD(B). In step 3, we use the algorithm from Lemma 12 as well as Koiran's algorithm to compute $\dim\ S_G$. We will state this more explicitly.

• The algorithm decides whether the instance is generically identifiable. If it is, it will automatically also obtain the parameters, since the ETR solver computes them on the way.

Answers to main questions:

1. The Gaussian recursive mixed graph models imply causal assumptions. The value (distribution) of each node is only affected by its parents and the bidirected edges. That yields Reichenbach's assumption, that if two nodes are not independent, there is either a direct cause (through their parents) or a common cause (through a bidirected edge). Some variant of the Markov assumption also holds if one considers the bidirected edges as parents.

2. There are different kinds of DAG models, with unobserved variables and without unobserved variables.

The one without unobserved variables is trivially identifiable. It is equivalent to the Gaussian recursive mixed graph models without bidirected edges.

The one with unobserved variables is more general than mixed graph models, so one can expect the lower bound to hold. The classic reduction of a mixed graph model to the DAG model is to replace any bidirected edge X <--> Y by a path through an "unobserved" node U, i.e. X <-- U --> Y. The covariance of U influences the covariances of X and Y. If the identification problem in the original mixed graph means solving equation system (1), then the reduced identification problem is solving an equation system (1') which is obtained from (1) by replacing each $\omega_{XY}$ with the product $\lambda_{XU}\lambda_{YU}$ where the solution must not contain any $\sigma$ involving the unobserved variables. Now a solution to (1) is also a solution to (1') after replacing the omegas (because (1) has no variables involving the unobserved variables).

Answers to technical questions:

3. You are right about lines 33 and 129, the index is wrong. Thanks for pointing this out, we will correct it.

4. That is right, we assume (as do all the previous results that we cite) that $\Sigma$ is given. We will point this out more clearly in the mentioned lines.

5. That is true, thanks for pointing this out as well. We will correct this.

6. The classes $\forall \mathbb{R}$ and $\exists R$ are co-classes in the sense that when a problem $L \in \forall \mathbb{R}$, then its complement $\bar L$ is in $\exists R$ and vice versa. The complement of QUAD$^{++}$ is the following problem: Given a set of quadratic equations with at least one solution, is this solution the only solution. It is a $\forall \mathbb{R}$ complete problem. By Lemma 3, our constructed instance is identifiable iff the given satisfiable set of quadratic equations has exactly one solution. So we have a reduction from the complement of QUAD$^{++}$, which is $\forall \mathbb{R}$ hard, to Numerical Identifiability and therefore, Numerical Identifiability is also $\forall \mathbb{R}$ hard.

7. Standard simulation means that PSPACE is contained in EXPTIME. It is easy to see that an algorithm that is polynomially space bounded is also exponential time bounded (modulo some small technical details). It can be found in standard text books on complexity, like Papadimitriou, Computational Complexity, Pearson, 1993. We will add the reference.

We agree it is not guaranteed that theoretically faster algorithms will be better in practice. We are working on an implementation of our algorithms. However, experience shows that drastic improvements on the theoretical running time typically come with similar improvements on the practical side. (But we admit that there are sometimes exceptions.)

We would like to emphasize that since the formulation of the parameter identification problem in linear SEMs in the 1960s, this problem has been widely studied by econometricians and sociologists, and more recently by the ML, AI, and statistics communities, but no (complete) algorithm has been found that matches the general applicability of the Gröbner basis approach and that runs faster than in doubly-exponential-time. For the first time, our work provides a new algorithm that states a constructive upper bound on the computational complexity of the generic identification problem which significantly improves the doubly-exponential bound of Gröbner bases.

We decided to publish our results in its current theoretical form because (1) it shows (for the first time) that it is possible to avoid the method based on Gröbner bases and (2) it proposes a new approach (using semi-algebraic sets and Koiran's algorithm) to solve the identification problem, which would be of interest to the NeurIPS community. We believe that our results also provide an important first step towards achieving an implementable algorithm that can handle cases with sizes much larger than those possible using state-of-the-art methods (see, e.g., Garcia-Puente et al. [22]).

For the related problem of numerical identifiability, we get an almost tight complexity classification. A $\forall \mathbb{R}$ lower bound, which is to our knowledge the first lower bound for any such identification problem, and a(n almost) matching upper bound, namely promise $\forall \mathbb{R}$ or $\forall \exists \mathbb{R}$.