Faster Generic Identification in Tree-Shaped Structural Causal Models

Thank your for your thoughtful comments:

Regarding weaknesses:

• We are not aware of any existing implementation of the algorithm by Gupta and Bläser. With some additional effort, we could have implemented it within our code framework. However, as already pointed out in our response to Reviewer 8jdN, the algorithm will not be practical even for moderate values of $n$. We provide more details in the answer to your first question below.

• Thank you for your suggestions for improving the presentation. We will take them into account and incorporate them. Since the paper lies at the intersection of algebraic statistics and causality and employs methods from algebraic algorithms and computer algebra, many concepts from different areas need to be introduced. Depending on the background of the actual reader, concepts might be more or less well-known. So we are really grateful for your input.

Regarding questions (our answer items follow the same order as the question items):

• We did not compare our algorithm with that of Gupta and Bläser for two main reasons. First, we are not aware of any existing implementation of their algorithm. Their result is primarily of theoretical interest, identifying a graph class for which generic identification is provably feasible in the sense of algorithm theory. Second, their algorithm, in its current form, is not practical even for moderate $n$. With some additional effort, we could have adapted our framework to implement the Gupta and Bläser algorithm. However, the running time of their algorithm is $O(n^6)$, while ours is $O(n^3 \log^2 n)$. As shown in Figure 5, for $n = 100$, our algorithm already performs clearly below $O(n^4)$ (with the break-even point around $n = 50$). Given that both algorithms use similar techniques, the constants hidden in the $O$-notation are likely comparable. Thus, our implementation is faster by a factor of at least $n^2$ on moderately sized instances — for example, a factor of $10.000$ for $n = 100$.

• Both factors contribute. The treeID algorithm has exponential worst-case running time as it brute-forces over all potential identifying cycles, which is one reason why our method is faster. Additionally, the implementation of identity tests in treeID appears to be inefficient. From our understanding of the treeID code, identity tests are performed by expanding polynomials and passing them to a computer algebra system. This approach can take exponential time (and even space) in the worst case. We believe this, at least partly, explains the irregular behavior observed when $n$ reaches 20.

• $\omega_{0,0}$ denotes the variance of the error term $\epsilon_0$ (see line 31) at the root node $0$. $\omega_{0,0} = \sigma_{0,0}$, since the root has no incoming directed edges. One could have used the notation $\omega_0$, but using double indices makes the summation formulas derived from Wright’s trek rule more elegant. (Similarly, $\omega_{i,i}$ denotes the variance of the error term at node $i$.) We will clarify this in the revised version.

• $t$ always means the length of an identifying cycle, which we are currently looking for. This is stated in line 321. In the layered graph in the self reduction part, time steps go from $1$ to $t$, so $t$ is the "total time". We understand your confusion and maybe it is better to rename $t$ into $\ell$ for length.

• Of course, running times are affected by the choice of programming language and other implementation details. However, these factors only influence the leading constants in the running time. If we have only one input size as in Figure 4, it is very hard to estimate the influence of the different effects. But as shown in Figure 5, we are also faster when $n$ increases, which means that the main effect comes from our faster algorithm. The reason we tested our algorithm on the benchmark of 879 graphs in Figure 4 was that the expected output was known from the paper by van der Zander et al. (2022), and we used it to verify correctness.

• The runtime of our algorithm is very stable, and two main features influence it: whether there are missing edges to the root, and the length of the identifying cycles. Beyond that, the specific graph structure has little impact. If there are missing edges to the root, we can identify almost immediately (as shown in Figures 6 and 7). After that, the runtime depends only on the length of the missing cycle, since we use iterated matrix products to search for all potential missing cycles in parallel. However, this influence is also small, since we use binary search and we always need to check for length $n/2$. In the worst case of length $n$, we will always branch in the more costly half of the binary search, but the overall effect is minor. Similar observations apply to the runtime of the treeID algorithm, aside from the inefficiencies in identity testing, and HTC. Execution times were consistent across multiple runs. For simplicity, we used a single input instance to generate the plots, as noted under "7. Experiment statistical significance" in the paper checklist. The actual edge distribution can influence the running time somewhat, but only lower order terms: the tree structure slightly affects variance computation (an $O(n^2)$ term), and the density of bidirected edges influences the self-reduction step ($O(n^3)$), though the overall impact is small (as seen in Figures 8, 9, and 10).

Regarding suggestions:

• Thank you, we will correct this.

• We adopted this convention from previous papers on the topic. It is convenient because nodes can be used directly as indices. However, we understand the potential confusion and will add a clarifying remark.

• We are using equations (3) and (4) on page 4 in the proof of Lemma 1, which are exactly your calculations. We should have mentioned this explicitly — thank you for pointing it out. Regarding Proposition 2: Let $n$ be the number of different prime factors. Since the smallest possible prime factor is $2$ and $z$ is the product of all its prime factors, $2^n \le z$, which gives the bound. We will add this to the proof.

• Thank you, we will correct this.

• This behavior can occur because matrix entries may be negative. For example, consider two disjoint paths between nodes $i$ and $j$ of the same length, where all weights are equal except for the last edges, which are negatives of each other. With a bit more care, one can construct such examples that stem from tree-shaped SCMs. We will add an example to the appendix of our paper.

Thank you for your thoughtful comments.

Regarding the weaknesses:

(1) It is true that tree-shaped SCMs represent a restricted class. However, the complexity of generic identification has remained an open problem for more than 20 years. In light of the (related) hardness results by Dörfler et al. (2024), the problem may indeed be computationally difficult. Therefore, it is reasonable to search either for criteria that are efficient but not complete, or for classes of SCMs for which the identification problem can be solved. We here follow the second path, as has also been done by van der Zander et al. (2022) or Gupta and Bläser.

(2) While it is true that some of the underlying techniques—such as the Schwartz-Zippel lemma or self-reduction—are well-known in algorithm theory, they are rarely used in practical implementations. Translating the algorithm by Gupta and Bläser into practice was a highly nontrivial task. In particular, the idea of choosing random values only once at the beginning is novel, and care must be taken with dependencies, especially in the self-reduction part, which required new local update techniques, too. We improve upon the algorithm by Gupta and Bläser by a factor of more than $n^2$ (nearly $n^3$ for large $n$), which is crucial for making it practical. It saves a factor of about $10.000$ for SCMs of size $n = 100$. (More details on this can be found in the answer to reviewer 8jdN).

(3) While it is true that we only used one graph per input size to draw the figures for simplicity, the running time of our algorithm does only depend on two parameters of the graph: whether there are missing edges to the root and the length of an identifying cycle. In the first case identification is very fast. In the second case, the most costly part of the algorithm uses linear algebra techniques (iterated matrix products) to find the identifying cycles and is oblivious to the actual edges in the SCM, since the product matrices get dense very fast. The longer the cycle, the more time we need, since we have to compute longer matrix products. However, this influence is also small, since we use binary search and always need to test for length $n/2$. If the cycle is long ($n$ in the worst case), then we always branch in the more costly part in the search. But the overall effect is small. There are some parts in the algorithm where the actual edges play some role (computation of covariances, self reduction), but they are of lower complexity and their influence is minor (as can be seen in Figures 7,8,9).

Regarding the questions:

(1) Finding efficient criteria for determining whether a graph is generically identifiable has been a long-standing open problem—see, for example, Foygel et al. (2012, Problem 1), which remains unresolved. (Full references provided in our paper.) Consequently, investigating subclasses of SCMs for which identifiability can be determined is a natural approach. Moreover, in cases involving hierarchical data, tree-shaped SCMs arise naturally—see, for instance, Thorson et al. (2023). Of course, in practice, the tree structure may not be perfect. In such cases, one could attempt to remove a few edges to obtain a tree. If this resulting SCM is not generically identifiable, then the original one will not be either. The converse, however, does not hold, and in such cases, the identification result must be manually transferred to the original instance.

(2) Generic identifiability is a property of an SCM, not of a particular data set. It means that all parameters are uniquely determined for almost all sets of observed covariances. (This is why, for example, Brito and Pearl (2002) refer to it as "almost everywhere identifiability.") Typically, one also obtains circuits that compute these parameters. So the "easy" answer to your question is that our algorithm will determine whether the perturbed SCM is generically identifiable, and if so, it will produce expressions for the parameters. More interesting, of course, is the question you intended: whether we can transfer the results from the perturbed SCM to the original one. Since we still lack a clear understanding of the complexity of generic identification, we are far from having general methods to determine how modifying a few edges affects the number of solutions. However, depending on the perturbations, certain properties might be preserved. For instance, adding directed edges will preserve nonidentifiability.

Thank you for your thoughtful comments.

Regarding weaknesses:

• The complexity of generic identifiability is a long-standing open problem: This means one can either look for methods that are not complete or for restricted graph classes for which we can design an efficient complete algorithm. This work takes the second route as have done others before, like van der Zander et al. (2022) or Gupta and Bläser (2024), all references are in the paper.

• While the algorithm by Gupta and Bläser is also polynomial time, our algorithm is the first that is practical. We are not aware of an implementation of the algorithm by Gupta and Bläser. With some extra effort, we could have also implemented it using our code framework. However, as already pointed out in the answer to reviewer 8jdN, the algorithm is not practical even for moderate $n$: The running time of the Gupta and Bläser algorithm is $O(n^6)$, while ours is $O(n^3 \log^2 n)$. Since the Gupta and Bläser algorithm uses similar techniques to ours, the constants hidden in the $O$-notation will be comparable. Therefore, our implementation is faster by a factor of at least $n^2$ on moderately sized instances, which for instance means a factor of $10000$ for $n = 100$.

Regarding questions (our answer items follow the same order as the question items):

• Thank you for providing this nice idea. Page 5 contains an outline of the algorithm at the top, but we will add a more detailed description with complexities.

• Finding efficient criteria for deciding whether an SCM is generically identifiable has been a long-standing open problem, see e.g. Foygel et al. (2012, Problem 1). Therefore, looking for subclasses of SCMs for which this is possible is a natural approach to this problem. Moreover, whenever we have hierarchical data, tree-shaped SCMs naturally occur, see for instance Thorson et al. (2023).

• Regarding benchmarking: We are happy to extend the benchmarking. However, the running time of our algorithm is very stable, it essentially depends on only two parameters: whether there are missing edges to the root, and the length of the identifying cycles. Beyond that, the specific graph structure has little impact. If there are missing edges to the root, we can identify almost immediately (as shown in Figures 6 and 7). After that, the runtime depends only on the length of the missing cycle, since we use iterated matrix products to search for all potential missing cycles in parallel. However, since we use binary search to find the length of the missing cycle (which is the most costly step), this also has little influence, since we always need to check for length $n/2$. If the length is $n$, then we will always branch in the more costly part during binary search, so this is the worst case. Similar observations apply to the runtime of the treeID algorithm, aside from the inefficiencies in identity testing, and HTC. Regarding your concrete questions:

  • For exemple, how do your algorithm perform when the size of the cycles of missing bidirected edges vary: In Figure 5 we use the worst case for our algorithm, only one missing cycle which contains all nodes. In the random graphs in Figures 8,9,10, the identifying cycles are short, since this happens with high probability in a random graph (see also Section D.2). These are the two extreme cases. One could add cases in between, these would however be quite comparable.
  • You explain that treeID performs especially bad because of the implementation of the identity test. Could you show how treeID with a better identity test implementation perform: This would require recoding treeID, since we think, after we inspected the code, that the current implementation uses a computer algebra system to perform the identity test, which would be hard to remove. Even with more efficient identity tests, the treeID algorithm still brute-forces over all potential cycles, which can be exponentially many.
  • The plotted graphs should also include the theoretical complexities of the different algorithms: One can do this, however, this would mean that one has to set the constant of the O-notation to some "guessed" value to fit the observed running times. But we agree that this is helpful for comparison reasons and we will add them.
  • What is the goal of adding the black plot which shows a worst version of your algorithm (with O(n^4) detection of identifying cycles)? I do not understand the added value: While $O(n^3 \log^2n)$ is asymptotically faster than $O(n^4)$, it is not clear whether this is true for small $n$ (due to hidden constants). (See lines 965-966 in the paper.) The plot confirms that the break-even point is indeed small.
  • For figure 5: You plot the execution time of SEMID HTC while saying that it cannot identify any of the cases. I do not understand the added value: This still tells you something about the running time of HTC, even if it fails to identify any parameters. In this case, it will perform $n$ max-cut computations, which is a lower bound for its running time in general, since it performs between $n$ and $n^2$ many max-cut computations depending on the instance.
  • I believe it would be more relevant to compare with another algorithm that is complete (e.g. the ID algorithm) ideally the most efficient one: There are only two complete methods known. One is based on Gröbner bases. This only works for small graphs. For instance, van der Zander et al. (2022) report that the Gröbner basis algorithm took more than four months (sic!) on the benchmark of Figure 4 (where we need 8s). The second one is based on the existential theory of the reals, see Dörfler et al. (2024). We are not aware of any implementation, but it will, like the Gröbner basis approach, only work for small $n$. So there are no efficient complete algorithms known for general graphs. In the light of the related hardness results by Dörfler et al. (2024), the problem might be hard in general. We do not know what you mean by "the ID algorithm". The most general algorithm we know is the ACID algorithm by Kumor et al. (2020). However, on tree-shaped SCMs it has the same problems as HTC, as shown by van der Zander et al. (2022), see lines 1180 - 1184. As ACID is an extension of HTC, it is more involved and slower.

• Thank you for pointing this out, we will add the bibliographic details of the reference to the abstract.

Thank you for your thoughtful comments.

Regarding the weaknesses:

W1: While we generally follow the strategy of Gupta and Bläser, we would like to emphasize that our improvements are nontrivial and crucial for making the algorithm practical. Moreover, one can interpret the soundness proof by Gupta and Bläser as suggesting that identifying cycles is, in some sense, the only viable approach to solving the problem.

The running time of the Gupta and Bläser algorithm is $O(n^6)$, whereas ours is $O(n^3 \log^2 n)$. Of course, the additional polylogarithmic factor in our bound is non-negligible for small and moderate values of $n$. However, as shown in Figure 5, for $n = 100$, our algorithm is already faster than $O(n^4)$ (as mentioned in the paper, the break-even point is around $n = 50$). Since the Gupta and Bläser algorithm employs techniques similar to ours, the constants hidden in the $O$-notation of a potential implementation of the Gupta and Bläser algorithm will be comparable to ours. Therefore, our implementation is faster by a factor of at least $n^2$ on moderately sized instances — for example, a factor of $10.000$ (!) for $n = 100$.

W2: We appreciate your comments regarding the write-up. Since the paper lies at the intersection of algebraic statistics and causality and employs methods from algebraic algorithms and computer algebra, many concepts from different areas need to be introduced. Our intention was to first provide a gentle introduction to the topic, then explain our main ideas for the improvements in Sections 4.1 and 4.2 using the example of covariance computation, and finally present the further improvements in Section 4.3, with more details provided in the appendix. But we see your points and will incorporate your suggestions.

Regarding the questions:

Q1: Generic identifiability is a property of an SCM, not of a particular data set. Given a tree-shaped SCM, we produce circuits (or determine that none exist) that compute the parameters for almost all possible observed sets of covariances. This does not directly apply to the scenario outlined in your question. However, one could attempt to find a tree-shaped SCM in a first step that approximates the observed scenario, and then use our methods to estimate the parameters. This seems like an interesting direction for future research. Thank you for raising this question.