arXiv CS open access

I-FLOP: Fast Learning of Order and Parents from Interventional Data

2026 · arxiv_cs

arXiv CS · Papers · 2026

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

Proceedings of Machine Learning Research TBD:1–28, 2026

Probabilistic Graphical Models (PGM)

I-FLOP: Fast Learning of Order and Parents from Interventional Data Liuting Chen†

jnk758@alumni.ku.dk

Alex Markham†

alex.markham@causal.dev

arXiv:2608.28245v1 [stat.ML] 28 Aug 2026

†

Department of Mathematical Sciences University of Copenhagen Denmark

Abstract We extend the FLOP (fast learning of order and parents) algorithm recently developed by Wienöbst et al. (2026) from observational to interventional data. In particular, we use the interventional BIC score of Hauser and Bühlmann (2012), adapting it to be used with the iterative Cholesky-based score updates that are partly responsible for FLOP’s speed. We show that, in the sample limit, I-FLOP recovers a DAG in the same interventional Markov equivalence class as the data-generating DAG. We compare I-FLOP to existing causal structure learning algorithms on simulated and real interventional data, where it performs favorably in terms of both performance and run time. Keywords: Causal Structure Learning, Interventional Data, Discrete Search

1. Introduction With observational data, we can typically only determine a Markov equivalence class (MEC) (Verma and Pearl, 1990)—a set of causal graphs that cannot be distinguished from such data alone. Interventions with known targets break some of these ambiguities, and with multiple interventional environments we obtain a finer equivalence class, the interventional Markov equivalence class (I-MEC) (Hauser and Bühlmann, 2012), which determines the direction of more edges than the MEC alone. To learn causal structures from interventional data, researchers have proposed various algorithms. These algorithms can be broadly categorized into three types based on how they utilize intervention information. The first type is score-based methods, such as Greedy Interventional Equivalence Search (GIES) (Hauser and Bühlmann, 2012): by defining a scoring function that incorporates intervention information, candidate graph structures are scored, and the optimal solution is searched within the graph space. The second type is permutation or ordering-based methods, such as Interventional Greedy Sparsest Permutation (I-GSP) (Wang et al., 2017): using intervention invariance constraints to guide the inference of node order, followed by selecting parent sets under the order constraints. The third type is invariance-based methods, such as Invariant Causal Prediction (ICP) (Peters et al., 2016), which leverage the principle of mechanism invariance where true causal mechanisms remain stable across multiple interventional environments while non-causal associations change. By identifying this invariance, causal relationships are pinpointed. © 2026 L. Chen & A. Markham.

Chen Markham

Among these approaches, ICP primarily tests invariance of predictive relationships for a specified response variable but is not designed as an interventional completed partially directed acyclic graph (I-CPDAG) search procedure. We therefore focus on GIES and I-GSP, which are more directly aligned with our goal of learning a global graph or an equivalence-class representative under known intervention targets. These two methods serve as the main reference points for our setting, but they face different challenges in practical use. GIES uses a local greedy strategy to search the space of I-MEC, gradually adjusting the graph structure through three phases: forward, backward, and turning, which respectively consider score-improving local modifications of the current equivalence-class representative. At each step within a phase, GIES selects the admissible local move that yields the largest improvement in the score. While this greedy approach is computationally efficient, it does not guarantee convergence to the globally optimal structure (Wang et al., 2017). More critically, the candidate search at each step may have exponential complexity (Hauser and Bühlmann, 2012) in the worst case, making it difficult to fully explore the search space in large-scale problems or scenarios with limited computational resources. I-GSP uses conditional independence (CI) tests under observational and interventional distributions to guide its permutation search. Its consistency result is asymptotic, whereas in finite samples the CI decisions can be sensitive to sample size, test choice, and significance threshold. Motivated by these challenges, we adopt an order-based score-and-search approach for known-target interventional data. Searching over node orderings provides an attractive alternative to direct structure search, as acyclicity is guaranteed by construction and large neighborhoods can be explored efficiently (Bouckaert, 1992; Singh and Valtorta, 1993; Larranaga et al., 1996; Friedman and Koller, 2003; Teyssier and Koller, 2005; Scanagatta et al., 2017). Building on this paradigm, FLOP (Wienöbst et al., 2026) combines decomposable scoring with efficient reinsertion-based search and local parent-set updates. We extend FLOP to known-target interventional data, yielding interventional FLOP (I-FLOP)1 . Our main contributions are threefold: 1) We introduce a target-filtered pooled local score that extends order-based search to multi-environment interventional data with known targets. 2) We establish both local and global theoretical consistency guarantees: we prove that the target-filtered grow–shrink procedure consistently recovers the optimal prefixrestricted parent sets from any initialization, and that exact global minimization of our target-filtered pooled score asymptotically identifies the true I-MEC. 3) We conduct synthetic and real-world experiments, showing that I-FLOP achieves competitive scalability and superior I-CPDAG accuracy compared with established interventional baselines, including GIES and I-GSP. 1. Code is available at https://github.com/L2-T2/I-FLOP

I-FLOP: Fast Learning of Order and Parents from Interventional Data

2. Preliminaries Let G = (V, EG ) be a directed acyclic graph (DAG) with node set V = {1, . . . , p} and edge set EG ⊊ {i → j | i, j ∈ V, i ̸= j}. We write PaG (v) for the parent set of node v. A linear Gaussian structural equation model (SEM) P (Peters and Bühlmann, 2014) associated with G is given by structural equations Xv = u∈PaG (v) Wuv Xu + εv where the noise variables εv are jointly independent Gaussian random variables with εv ∼ N (0, σv2 ) and σv2 > 0. We consider multiple environments E = {0, 1, . . . , m} associated with a known intervention target family I = {Ie : e ∈ E}, where I0 = ∅ denotes the observational environment. For a target set Ie ⊆ V , a known-target hard intervention replaces the structural equations of the targeted variables so that they are instead drawn jointly from an intervention (e) (e) (e) distribution QIe , i.e., XIe ∼ QIe , and leaves the structural equations of all non-targeted variables v ∈ V \ Ie unchanged. Every DAG can be associated to at least one topological order τ such that u → v ∈ EG implies that u precedes v in τ . Given τ , the candidate parent set Pv of node v consists only of nodes that precede v in the order, which can be represented as Pv ⊆ Predτ (v) = {u ∈ V : u precedes v in τ }. The ordered Markov property states that Xv ⊥⊥ XPredτ (v)\PaG (v) | XPaG (v) for every v ∈ V . Since the linear Gaussian SEM above has strictly positive noise variances, the induced joint Gaussian distributions are non-degenerate and hence positive. Therefore, the ordered, local, and global Markov properties are equivalent in our setting (Lauritzen, 2026). A joint distribution P over XV is Markov with respect to G if it satisfies these equivalent Markov properties. Observational conditional independence relations do not generally identify a unique DAG: different DAGs may encode the same set of conditional independence relations and are therefore observationally indistinguishable. Such DAGs form a Markov equivalence class (MEC), represented by a completed partially directed acyclic graph (CPDAG) (Andersson et al., 1997). Following Hauser and Bühlmann (2012), two DAGs G and H are interventionally Markov equivalent with respect to the target family I if they induce the same interventional Markov model MI (G) = MI (H). The corresponding equivalence class is the I-MEC. For the true DAG G⋆ , we denote the true I-MEC by GI⋆ = {G ∈ GDAG : MI (G) = MI (G⋆ )}, where GDAG denotes the set of all DAGs on V . The unique graphical representative of an I-MEC is the interventional essential graph, also called the I-CPDAG. Let D = {X (e) }e∈E denote the multi-environment dataset. We propose I-FLOP as an extension of FLOP P that retains FLOP’s search framework for optimizing a decomposable score S(G; D) = v∈V sv (PaG (v); D) over the space of topological orders (Algorithm 1). Given a topological order τ , each node v selects its parent set PaG (v) only from its prefix Predτ (v), thereby inducing a DAG G. Candidate orders are explored by reinsertion moves (Algorithm 2). During the evaluation of a reinsertion move, the affected parent sets and local scores are updated locally using a warm-start grow-shrink procedure (Algorithm 3), avoiding a full recomputation for all nodes. This reinsertion-move and warm-start growshrink machinery is unchanged from FLOP; the component that adapts it to interventional data is the target-filtered pooled local score sv , introduced next, which is the score these procedures optimize. After the reinsertion-based local search reaches a local optimum, Iterated Local Search (ILS) applies a random perturbation to the current best order and 3

Chen Markham

restarts local search from the perturbed order; the best-scoring order found across these restarts is retained (Lourenço et al., 2018; De Campos et al., 2003).

3. I-FLOP: order-based search for interventional data Algorithm 1: Target-filtered order search for known-target interventional data Input: Multi-environment data {X (e) }e∈E , target family I = {Ie }e∈E , initial order τ0 b and its induced I-CPDAG C bI Output: Estimated DAG G // Phase 1: Precompute environment statistics foreach e ∈ E do Center X (e) within environment e and compute scatter matrix C (e) end // Phase 2: Construct target-filtered pooled scores foreach v ∈ V do Ev ← {e / Ie } P∈ E : v ∈ Nv ← e∈Ev ne P Cvpool ← e∈Ev C (e) from Cvpool and Nv Define sv := stf-pool v end // Phase 3: Initialize parent sets and local scores under τ0 τ ← τ0 foreach v ∈ V do Zv ← Predτ (v) ℓ∅,v ← sv (∅) (Pv , ℓv ) ← GrowShrinktf (v, Zv , ∅, ℓ∅,v , 0) // Algorithm 3 end Let P := (Pv )v∈V and ℓ := (ℓv )v∈V // Phase 4: Reinsertion-based local order search improved ← true while improved do (τ old , ℓold ) ← (τ, ℓ) foreach v ∈ τ old do (τ, P, ℓ) ← Reinserttf (τ, P, ℓ, v) end P P old improved ← v∈V ℓv < v∈V ℓv end // Phase 5: Graph construction and output b ← V, {u → v : u ∈ Pv , v ∈ V } G bI ← I-CPDAG(G, b I) C b C bI ) return (G, 4

// Algorithm 2

I-FLOP: Fast Learning of Order and Parents from Interventional Data

Algorithm 1 presents the core reinsertion-based local order search. In our implementation, the optional ILS wrapper described above repeatedly perturbs the locally optimal order and reruns the reinsertion search, retaining the best-scoring solution across restarts. We first state two assumptions for the known-target hard-intervention setting: one on the data-generating process and one on the conservativeness of the intervention target family. Assumption 1 (Data-generating process and causal sufficiency) We assume causal sufficiency over V and a linear Gaussian SEM with true DAG G⋆ = (V, E). In each environment e ∈ E, the known hard-intervention targets Ie ⊆ V satisfy (e) (e) (e) XIe ∼ QIe , while each non-targeted node v ∈ / Ie retains its structural equation Xv = P (e) (e) (e) ⋆ ∼ N (0, (σv⋆ )2 ). The noise variables are mutually u∈paG⋆ (v) Wuv Xu + εv , where εv independent, have positive variances, and are invariant across environments. Assumption 2 (Conservative family of targets) Following Hauser and Bühlmann (2012), we assume the intervention target family I = {Ie : e ∈ E} is conservative, meaning that for every node v ∈ V , there exists at least one environment e ∈ E such that v ∈ / Ie . Intuitively, this assumption ensures that no variable is perpetually perturbed across all experimental settings, so that the unperturbed, natural causal mechanism P (v | pa(v)) of each variable can be observed and learned from at least one environment. 3.1. Valid environments under interventions In the observational setting, FLOP scores each candidate parent set by fitting a single local regression on samples assumed to come from one observational regime. This is no longer appropriate for interventional multi-environment data: different environments may induce different joint distributions, and samples from environments where node v is directly intervened upon do not follow the original mechanism of v. The invariant object is instead the local mechanism of v across environments in which v is not targeted. This leads to response-target filtering: Definition 1 For each node v ∈ V , we define the valid environment set Ev := {e ∈ E : v ∈ / Ie }, which consists of the environments in Pwhich v is not directly intervened upon, and the corresponding total sample size is Nv := e∈Ev ne . We use N to denote the total sample size over all environments. The valid environment set Ev is node-specific: different response nodes may use different subsets of environments for local scoring. Importantly, the validity of an environment for node v depends only on whether v itself is directly intervened upon. In particular, let P ⊆ V \ {v} be a candidate parent set for v. If u ∈ P is directly intervened upon in environment e, while v ∈ / Ie , then e ∈ Ev remains valid for scoring v. By Assumption 1, intervening on u may change the (e) distribution of Xu , but it does not replace the structural equation of v when v ∈ V \ Ie . Therefore, observations from all environments in Ev can be pooled to estimate the same structural equation for v. We next use these node-specific valid environments to define the target-filtered pooled local score. 5

Chen Markham

3.2. Target-filtered pooled local score Definition 2 For node v ∈ V and candidate parent set P ⊆ V \{v}, define the target-filtered pooled local score by stf-pool (P ) = v

1 Nv 2 1 + log σ bv|P + log N (|P | + 1). 2 2

2 Here, σ bv|P is estimated only from the valid environments Ev .

For each environment, the centered data and scatter matrix are e (e) , e (e) )⊤ X C (e) = (X

e (e) = X (e) − X̄ (e) , X

P and the pooled scatter matrix for node v over its valid environments is Cvpool = e∈Ev C (e) . We impose the following non-singularity assumption on the relevant principal submatrices of Cvpool . Assumption 3 (Local non-singularity) For every node v ∈ V and every candidate parent set P ⊆ V \ {v} considered in the local score evaluation, the principal submatrix (Cvpool )P P is non-singular. 2 Under Assumption 3, the pooled residual variance estimator σ bv|P for node v and candidate parent set P is defined as 2 σ bv|P =

i 1 h pool (Cv )vv − (Cvpool )vP ((Cvpool )P P )−1 (Cvpool )P v , Nv

where the second term is set to zero when P = ∅ (the Schur complement of the pooled 2 scatter matrix; see the Frisch–Waugh–Lovell identity used in Appendix A). Thus, σ bv|P corresponds to the residual variance obtained by fitting a shared local linear mechanism for v across its valid environments, rather than across all environments. By construction, this score preserves node-wise decomposability and is compatible with order-restricted parent sets. For a candidate DAG G′ , the total score is Stf-pool (G′ ; D, I) =

stf-pool (PaG′ (v)). v

v∈V

Given a fixed order τ , score-based I-FLOP selects parent sets by Pv (τ ) ∈ arg min stf-pool (P ), v P ⊆Predτ (v)

and induces the order-restricted DAG Gτ = {u → v : u ∈ Pv (τ )}. The corresponding order objective is Qscore (τ ) =

X v∈V

stf-pool (Pv (τ )). v

I-FLOP: Fast Learning of Order and Parents from Interventional Data

3.3. Warm-start parent-set selection For each node v and prefix Z, I-FLOP updates the candidate parent set using a grow–shrink procedure with the target-filtered pooled local score (Algorithm 3).2 This strategy is inherited from FLOP, which recovers the restricted Markov boundary under a single observational distribution. In our setting, however, valid environments preserve the response mechanism but may exhibit different joint distributions due to interventions on other nodes. Consequently, the restricted Markov boundary characterization does not directly apply to multiple data sources. To characterize the asymptotic objective optimized by the grow–shrink procedure, we introduce a pooled population residual risk over the valid environments. The following assumption specifies the asymptotic sample proportions contributed by each environment. Assumption 4 (Environment proportions) For every environment e ∈ E, the sample proportion P satisfies ne /N → πe , where πe > 0. Consequently, for every node v, Nv /N → αv := e∈Ev πe > 0. e (e) denote the covariance matrix For each environment e ∈ E, let Σ(e) = Cov X of the environment-centered variables. Under Assumptions 1, 2, and 4, the normalized valid-environment pooled scatter matrix Nv−1 Cvpool converges in probability to Σpool = v P 1 (e) (by the weak law of large numbers applied within each valid environment, e∈Ev πe Σ αv combined via Assumption 4’s sample-proportion limits). For a candidate parent set P ⊆ Z, the corresponding population pooled residual risk is −1 pool 2 σv|P,pool = (Σpool )vv − (Σpool )vP (Σpool )P P (Σv )P v . v v v 2 This quantity is the population counterpart of σ bv|P and characterizes the asymptotic fit of the shared response mechanism across the valid environments. It therefore induces the population target of the grow–shrink procedure:

Definition 3 The pooled local parent target of v relative to Z is defined as the sparsest candidate parent set in Z that gives the optimal valid-environment structural-mechanism fit: 2 Mvpool (Z) ∈ arg min σv|P,pool , P ⊆Z

where ties are resolved in favor of a sparsest minimizer. Although single-variable greedy search procedures could in principle be susceptible to local optima, the underlying linear Gaussian model and positive definiteness of the pooled covariance guarantee that the pooled residual-risk landscape has no local traps: Proposition 4 Fix v ∈ V and Z ⊆ V \ {v}. Under Assumptions 1, 2, and 4, the 2 pooled residual risk criterion R(P ) = σv|P,pool has a unique sparsest minimizer Mvpool (Z) = −1 pool supp(wZ ), where wZ = (Σpool ZZ ) ΣZv . Furthermore, for every candidate set P ⊆ Z with 2 2 Mvpool (Z) ⊈ P , there exists u ∈ Mvpool (Z) \ P satisfying σv|P ∪{u},pool < σv|P,pool .

2. Proofs of all theoretical results in this section are provided in Appendix A.

Chen Markham

Proposition 4 establishes that whenever the current candidate parent set P does not fully cover the optimal parent set Mvpool (Z), there always exists at least one missing parent whose addition strictly reduces the pooled residual risk. Consequently, target-filtered grow–shrink provably converges to the global parent target from any initialization: Theorem 5 Fix v ∈ V and Z ⊆ V \ {v}. Let Pbv,N (Z) be the parent set returned by target-filtered grow–shrink using stf-pool . Under Assumptions 1–4, starting from any initial v parent set P0 ⊆ Z, P Pbv,N (Z) = Mvpool (Z) −→ 1 as N → ∞. Since a fixed order τ induces a fixed-prefix problem for each node, applying Theorem 5 node-wise immediately yields the following fixed-order recovery guarantee: Corollary 6 Fix an order τ , and let Zv (τ ) = Predτ (v) for each v ∈ V . For each node v, let Pbv,N (τ ) be the parent set returned by target-filtered grow–shrink with prefix Zv (τ ), b N (τ ) = {u → v : u ∈ Pbv,N (τ ), v ∈ V }. Then, under and define the induced graph G Assumptions 1–4, P PaGbN (τ ) (v) = Mvpool (Zv (τ )) for all v ∈ V −→ 1 as N → ∞. Theorem 5 and Corollary 6 establish the large-sample behavior of the target-filtered grow–shrink updates for a fixed prefix and a fixed order, respectively. These results justify the local parent-set selection used in I-FLOP. We next show that the target-filtered pooled score admits a valid likelihood interpretation under known interventions and use this connection to state its global recovery property. 3.4. Global recovery from the target-filtered pooled score Under known hard interventions, the replacement mechanism of a directly targeted node does not depend on its candidate parent set in the original structural mechanism and therefore contributes only a graph-independent nuisance term. After profiling out environmentspecific intercepts, equivalently by working with the environment-centered observations, let XX ℓN (G, θG ) = ℓv,e PaG (v); θv , v∈V e∈Ev

denote the resulting log-likelihood contribution of the original structural mechanisms, evaluated only over the valid environments of each response node. Here, θG = {θv : v ∈ V }, where θv contains the regression coefficients associated with PaG (v) and the residual variance of the original mechanism of node v. Lemma 7 Under Assumptions 1–3, the target-filtered pooled graph score is the decomposable score associated with ℓN (G, θG ), up to terms independent of the candidate DAG. In particular, for any candidate DAG G′ , X 1 Stf-pool (G′ ; D, I) = − sup ℓN (G′ , θG′ ) + log N | PaG′ (v)| + 1 + CN . 2 θG′ v∈V 1P ′ where CN = − 2 v∈V Nv log(2π) does not depend on G . 8

I-FLOP: Fast Learning of Order and Parents from Interventional Data

We now characterize the global target of exact score minimization in terms of I-MEC. For a candidate DAG G, define the population counterpart of the graph-dependent fit component of the target-filtered pooled score as X αv 2 , G ∈ GDAG . L∞ (G) = log σv|Pa G (v),pool 2 v∈V

Let A = arg min L∞ (G),

Amin = arg min |EG |.

G∈GDAG

G∈A

The restriction to Amin corresponds to the graph-dependent complexity term X | PaG (v)| + 1 = |EG | + |V |. v∈V

By definition, A contains the DAGs attaining the optimal population fit of the validenvironment structural mechanisms, and Amin retains the sparsest DAGs among these representations. To conclude that the members of Amin belong to the true I-MEC, we require the conditional independence relations induced by the known intervention targets to coincide with those encoded by the true DAG G⋆ . We therefore impose the following interventional faithfulness assumption. Assumption 5 (Interventional faithfulness assumption) Let GI⋆ := {G ∈ GDAG : MI (G) = MI (G⋆ )} denote the true I-MEC. The interventional distribution family generated by G⋆ is faithful to G⋆ with respect to the known intervention target family I. Theorem 8 Under Assumptions 1, 2, 3, 4 and 5, the sparsest global residual-risk minimizers belong to the true I-MEC: Amin ⊆ GI⋆ = {G ∈ GDAG : MI (G) = MI (G⋆ )} . Theorem 8 is a global statement. For the empirical score, uniform convergence of the target-filtered residual risks and the BIC sparsity penalty imply that an exact global minimizer selects a DAG in Amin in the limit N → ∞. Combining this with Amin ⊆ GI⋆ gives the following I-CPDAG recovery result: b N be an exact global minimizer of Stf-pool (G; D, I) over DAG G. Under Corollary 9 Let G b N belongs to the true I-MEC with probability tending to one. Assumptions 1, 2, 3, 4, and 5, G b N = G⋆ ) → 1. If the target family I fully identifies the DAG, so that GI⋆ = {G⋆ }, then Pr(G

4. Experiments In this section, we evaluate I-FLOP for interventional causal discovery. We compare it with two external baselines, GIES and I-GSP, and two FLOP-based variants: FLOP-obs, the original observational method, and FLOP-envwise, which uses environment-wise variance estimation without target-based environment filtering. We generate synthetic data 9

Chen Markham

Figure 1: I-CPDAG SHD decreases as sample size increases from 102 to 106 per environment with 20 nodes and expected degree 2 for one observational and five known-target singleton hard intervention environments.

by Erdős–Rényi (ER) DAGs (Erdős and Rényi, 1959) under linear Gaussian SEMs. Each dataset contains one observational environment and multiple known-target hard interventional environments. The experiment can be divided into three main parts. First, we assess finite-sample convergence, specifically whether the algorithms recover the I-CPDAG (Section 4.1). Second, we evaluate the algorithms’ performances under well-specified synthetic settings (Section 4.2). Finally, we demonstrate performance on real interventional data from the Causal Chamber Light Tunnel system (Section 4.3). 4.1. Finite sample convergence toward the asymptotic guarantee We first investigate the finite-sample convergence behavior of I-FLOP as the sample size increases. Our theoretical results establish that, in the large-sample limit, optimizing the target-filtered pooled score recovers the true I-CPDAG. This experiment therefore examines whether the structural estimates produced by I-FLOP become increasingly accurate and stable with more data. Figure 1 shows that I-FLOP recovers the I-CPDAG in large sample limit, which supports our claim. 4.2. Main accuracy and runtime comparison Figure 2 compares accuracy and runtime in a representative well-specified synthetic setting. I-FLOP achieves the lowest I-CPDAG SHD in both ER graph regimes, with a clearer advantage on denser graphs, where GIES and I-GSP have substantially larger structural errors. I-FLOP is also markedly faster than GIES in the denser setting while remaining substantially more accurate than the faster I-GSP baseline, yielding a favorable accuracy– runtime trade-off. Table 1 shows the decomposition comparison of I-CPDAG SHD between our I-FLOP and the external baselines. 10

I-FLOP: Fast Learning of Order and Parents from Interventional Data

Figure 2: The I-CPDAG Structural Hamming Distance (SHD) is plotted against runtime for linear-Gaussian ER DAGs. Each setting uses one observational environment and five known-target singleton hard intervention environments, with 1,000 samples per environment.

Table 1: SHD decomposition of experiment in well-specified synthetic setting Expected degree d = 3 Algorithm

Mean

Std.

Expected degree d = 5

Skel. Orient. Mean

I-FLOP (ILS=0) 2.50 3.78 1.30 I-FLOP (ILS=20) 0.70 0.67 0.60 I-FLOP (ILS=50) 0.70 0.67 0.60 GIES 2.70 3.13 2.00 I-GSP 28.40 16.93 19.60

Std.

Skel. Orient.

1.20 1.60 1.58 0.90 0.10 1.50 1.27 1.10 0.10 1.50 1.27 1.10 0.70 47.40 49.02 39.10 8.80 94.70 28.79 73.10

0.70 0.40 0.40 8.30 21.60

4.3. Real-data evaluation on Causal Chamber We evaluate algorithms on real interventional data from the Causal Chamber Light Tunnel system (Gamella et al., 2025). The view we use contains 20 variables, including one observational environment (10,000 samples) and five singleton hard-intervention environments with known targets (red, green, blue, pol 1, pol 2, each containing 1,000 samples), all derived from lt interventions standard v1. Figure 3 shows that I-FLOP exhibits the best accuracy–runtime trade-off on this real interventional dataset.

5. Conclusion We introduce I-FLOP, an order-based score search method for causal structure learning from multi-environment interventional data with known targets. By replacing FLOP’s observational local score with a target-filtered pooled interventional score, I-FLOP retains the efficiency of warm-start grow–shrink updates and reinsertion-based order search while incorporating intervention information in a principled manner. Under Assumptions 1– 11

Chen Markham

Figure 3: I-FLOP (ILS = 0) attains the lowest mean I-CPDAG SHD while remaining faster than GIES. All I-FLOP configurations outperform FLOP-obs and I-GSP, and positive ILS budgets reduce variability (details in Appendix D).

5, exact global minimization of the target-filtered pooled score asymptotically recovers the true I-MEC, while under Assumptions 1–4, the fixed-prefix target-filtered grow–shrink update recovers the pooled local parent target in the large-sample limit. Empirically, IFLOP achieves a favorable accuracy–runtime trade-off on synthetic interventional data and competitive performance on the Causal Chamber benchmark. The main limitations of the current framework stem from its reliance on known hard interventions, linear-Gaussian mechanisms with strictly positive error variances, and conservative intervention designs. These assumptions restrict the model class to stochastic linear systems and exclude deterministic or functional relationships with degenerate noise, as well as nonlinear or uncertain intervention mechanisms. This mismatch is particularly relevant for physical systems such as the Causal Chamber, where variables may obey deterministic conservation laws or bounded nonlinear responses, and helps explain why more aggressive score optimization does not necessarily yield lower SHD when the data-generating process departs from the assumed model class. Future work could therefore extend target-filtered order search to deterministic or functional relationships, nonlinear additive noise models, and soft, uncertain, or partially known intervention targets (i.e., settings where the identity of the intervened variable(s) in a given environment is not fully or reliably known, as considered by Squires et al. (2020) and Gamella et al. (2022)), potentially through more flexible non-parametric or generalized scoring criteria.

Acknowledgments A.M. was partially supported by Novo Nordisk Foundation grant NNF20OC0062897 and also by the Pioneer Centre for Statistical and computational Methods for Advanced Re12

I-FLOP: Fast Learning of Order and Parents from Interventional Data

search to Transform Biomedicine (SMARTbiomed), DNRF grant number P4. We thank Ming Cai from the Graduate School of Informatics, Kyoto University, for helpful discussions and comments on an earlier version of the manuscript.

References Steen A. Andersson, David Madigan, and Michael D. Perlman. A characterization of Markov equivalence classes for acyclic digraphs. The Annals of Statistics, 25(2):505–541, 1997. doi: 10.1214/aos/1031833662. Remco R. Bouckaert. Optimizing causal orderings for generating DAGs from data. In Uncertainty in Artificial Intelligence, pages 9–16. Morgan Kaufmann, 1992. doi: 10. 1016/B978-1-4832-8287-9.50006-2. Luis M. De Campos, Juan M. Fernández-Luna, and J. Miguel Puerta. An iterated local search algorithm for learning Bayesian networks with restarts based on conditional independence tests. International Journal of Intelligent Systems, 18(2):221–235, 2003. doi: 10.1002/int.10085. Paul Erdős and Alfred Rényi. On random graphs I. Publicationes Mathematicae Debrecen, 6:290–297, 1959. doi: 10.5486/PMD.1959.6.3-4.12. Nir Friedman and Daphne Koller. Being Bayesian about network structure. a Bayesian approach to structure discovery in Bayesian networks. Machine Learning, 50:95–125, 2003. doi: 10.1023/A:1020249912095. Juan L. Gamella, Armeen Taeb, Christina Heinze-Deml, and Peter Bühlmann. Characterization and greedy learning of Gaussian structural causal models under unknown interventions. arXiv preprint arXiv:2211.14897, 2022. doi: 10.48550/arXiv.2211.14897. Juan L. Gamella, Jonas Peters, and Peter Bühlmann. Causal chambers as a real-world physical testbed for AI methodology. Nature Machine Intelligence, 7(1):107–118, 2025. doi: 10.1038/s42256-024-00964-x. Alain Hauser and Peter Bühlmann. Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs. Journal of Machine Learning Research, 13(79):2409–2464, 2012. URL http://jmlr.org/papers/v13/hauser12a.html. Alain Hauser and Peter Bühlmann. Jointly interventional and observational data: Estimation of interventional Markov equivalence classes of directed acyclic graphs. Journal of the Royal Statistical Society Series B: Statistical Methodology, 77(1):291–318, 01 2015. doi: 10.1111/rssb.12071. P. Larranaga, C.M.H. Kuijpers, R.H. Murga, and Y. Yurramendi. Learning Bayesian network structures by searching for the best ordering with genetic algorithms. IEEE Transactions on Systems, Man, and Cybernetics - Part A: Systems and Humans, 26(4):487–493, 1996. doi: 10.1109/3468.508827. 13

Chen Markham

Steffen L. Lauritzen. Graphical Models. Oxford University Press, 2nd edition, June 2026. doi: 10.1093/9780191775284.001.0001. Helena Ramalhinho Lourenço, Olivier C. Martin, and Thomas Stützle. Iterated local search: Framework and applications. In Handbook of metaheuristics, pages 129–168. Springer, 2018. doi: 10.1007/978-3-319-91086-4 5. Michael C. Lovell. A simple proof of the FWL theorem. Journal of Economic Education, 39(1):88–91, 2008. doi: 10.3200/JECE.39.1.88-91. Jonas Peters and Peter Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101(1):219–228, 2014. doi: 10.1093/biomet/ ast043. Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen. Causal inference by using invariant prediction: Identification and confidence intervals. Journal of the Royal Statistical Society Series B: Statistical Methodology, 78(5):947–1012, 2016. doi: 10.1111/rssb.12167. Mauro Scanagatta, Giorgio Corani, and Marco Zaffalon. Improved local search in bayesian networks structure learning. In Proceedings of The 3rd International Workshop on Advanced Methodologies for Bayesian Networks, volume 73 of Proceedings of Machine Learning Research, pages 45–56. PMLR, 20–22 Sep 2017. URL https://proceedings.mlr. press/v73/scanagatta17a.html. Moninder Singh and Marco Valtorta. An algorithm for the construction of Bayesian network structures from data. In Proceedings of the Ninth Conference on Uncertainty in Artificial Intelligence (UAI), pages 259–265, 1993. doi: 10.48550/arXiv.1303.1485. Chandler Squires, Yuhao Wang, and Caroline Uhler. Permutation-based causal structure learning with unknown intervention targets. In Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), volume 124 of Proceedings of Machine Learning Research, pages 1039–1048. PMLR, 2020. URL https://proceedings.mlr.press/ v124/squires20a.html. Marc Teyssier and Daphne Koller. Ordering-based search: A simple and effective algorithm for learning Bayesian networks. In Proceedings of the Twenty-First Conference on Uncertainty in Artificial Intelligence, UAI’05, pages 584–590, Arlington, Virginia, USA, 2005. AUAI Press. doi: 10.48550/arXiv.1207.1429. Thomas Verma and Judea Pearl. Equivalence and synthesis of causal models. In Proceedings of the Sixth Annual Conference on Uncertainty in Artificial Intelligence (UAI), pages 255–270. Elsevier Science Inc., 1990. doi: 10.48550/arXiv.1304.1108. Yuhao Wang, Liam Solus, Karren Yang, and Caroline Uhler. Permutationbased causal inference algorithms with interventions. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL https://proceedings.neurips.cc/paper_files/paper/2017/file/ 275d7fb2fd45098ad5c3ece2ed4a2824-Paper.pdf. 14

I-FLOP: Fast Learning of Order and Parents from Interventional Data

Marcel Wienöbst, Leonard Henckel, and Sebastian Weichwald. Embracing discrete search: A reasonable approach to causal structure learning. In International Conference on Learning Representations, volume 2026, pages 54630–54652, 2026. URL https://openreview. net/forum?id=lejOV6j3cj.

Chen Markham

Appendix A. Proofs of Theoretical Results A.1. Single-Move Score Comparisons We first establish the large-sample score comparisons which will be used in the proof of Theorem 5. Lemma 10 Fix v ∈ V and Z ⊆ V \ {v}. Under Assumptions 1–4, the following statements hold simultaneously over the finitely many candidate parent sets contained in Z with probability tending to one. 2 2 1. If P, Q ⊆ Z satisfy σv|Q,pool < σv|P,pool , then stf-pool (Q) < stf-pool (P ). v v tf-pool 2 2 2. If u ∈ P ⊆ Z and σv|P (P \ {u}) < stf-pool (P ). v \{u},pool = σv|P,pool , then sv

Proof By Assumptions 1 and 4, the environment-wise sample covariances converge to their population counterparts. Since ne /N → πe > 0 for every e ∈ E, it follows that, for every fixed P ⊆ Z, p 2 2 σ bv|P → − σv|P,pool . Since Z is fixed, there are only finitely many candidate parent sets, and this convergence holds simultaneously over all P ⊆ Z. 2 2 For the first statement, suppose that σv|Q,pool < σv|P,pool . Then 2 p 1 σv|Q,pool 1 tf-pool < 0. sv (Q) − stf-pool (P ) log → − v 2 Nv 2 σv|P,pool

The fit contribution is of order Nv , whereas the difference in complexity penalties is of order log N , and Nv /N → αv > 0. Therefore, any move that strictly improves the population fit decreases the empirical local score with probability tending to one. For the second statement, suppose that deleting u leaves the population residual risk unchanged. The two nested local models then attain the same population optimum for the centered valid-environment criterion. Under Assumptions 1–4, a standard fixed-dimensional quadratic expansion of the maximized criterion gives an Op (1) fit improvement from retaining the redundant parent u. Retaining u, however, incurs the additional penalty 12 log N , which diverges. Therefore, stf-pool (P \ {u}) < stf-pool (P ). v v

A.2. Proof of Proposition 4 Proof Fix v ∈ V and Z ⊆ V \ {v}. Under Assumptions 1, 2 and 4, the valid-environment P pooled covariance matrix Σpool = α1v e∈Ev πe Σ(e) is positive definite, since each Σ(e) > 0 v

and πe /αv > 0 for every e ∈ Ev . Hence every principal submatrix of Σpool , in particular v Σpool , is positive definite and therefore invertible. ZZ 16

I-FLOP: Fast Learning of Order and Parents from Interventional Data

Define the pooled population regression coefficient vector of Xv on XZ by −1 pool wZ := (Σpool ZZ ) ΣZv ,

S := supp(wZ ).

Since Σpool ZZ > 0, the least-squares coefficient vector wZ is unique. For each P ⊆ Z, let 2 R(P ) := σv|P,pool denote the pooled population residual variance from regressing Xv on XP . We first show that R(P ) = R(Z) if and only if S ⊆ P . If S ⊆ P , then wZ vanishes on Z \ P , so the optimal predictor based on XZ depends only on XP . Hence R(P ) ≤ R(Z). Since P ⊆ Z, enlarging the regressor set cannot increase the least-squares residual variance, so R(Z) ≤ R(P ). Therefore R(P ) = R(Z). Conversely, suppose that R(P ) = R(Z). Extend the least-squares coefficient vector for the regression on XP to a vector indexed by Z by assigning zero coefficients to Z \ P . This extended vector achieves residual variance R(Z), and is therefore also a minimizer of the least-squares problem over XZ . By uniqueness of wZ , it must coincide with wZ . Thus (wZ )u = 0 for every u ∈ Z \ P , which implies S ⊆ P . Therefore, the minimizers of R(P ) over P ⊆ Z are exactly the supersets of S. In particular, S is the unique minimizer of smallest cardinality. We have Mvpool (Z) = S = supp(wZ ). It remains to prove the absence of local traps. Let P ⊆ Z satisfy Mvpool (Z) ⊈ P , and set A := Mvpool (Z) \ P , so that A ̸= ∅. Suppose, toward a contradiction, that R(P ∪ {u}) = R(P ) for every u ∈ A. For u ∈ A, let rvP and ruP denote the population residuals obtained by projecting Xv and Xu , respectively, onto XP under the pooled distribution. By the Frisch–Waugh–Lovell theorem (Lovell, 2008), equivalently the Schur complement identity, R(P ) − R(P ∪ {u}) =

Covpool (rvP , ruP )2 . Varpool (ruP )

Since Σpool > 0, we have Varpool (ruP ) > 0. Hence the assumption R(P ∪ {u}) = R(P ) v implies Covpool (rvP , ruP ) = 0 for every u ∈ A. P denote the vector of residuals obtained by projecting X Let rA A onto XP . Then P P Covpool (rA , rv ) = 0. Applying the multivariate Frisch–Waugh–Lovell identity gives P P −1 P P R(P ) − R(P ∪ A) = Covpool (rvP , rA ) Varpool (rA ) Covpool (rA , rv ) = 0, P ) > 0 by positive definiteness of the corresponding Schur complement. where Varpool (rA Thus R(P ∪ A) = R(P ). On the other hand, P ∪ A = P ∪ Mvpool (Z) ⊇ Mvpool (Z) = S. By the characterization above, this implies R(P ∪ A) = R(Z). Therefore R(P ) = R(Z), which in turn implies S ⊆ P . Since S = Mvpool (Z), this contradicts the assumption Mvpool (Z) ⊈ P . 2 2 Hence there exists at least one u ∈ Mvpool (Z) \ P such that σv|P ∪{u},pool < σv|P,pool .

Chen Markham

A.3. Proof of Theorem 5 Proof Let Pgrow denote the parent set obtained at the end of the grow phase. We first show that M ⊆ Pgrow . Suppose otherwise that M ⊈ Pgrow . By Proposition 4, there exists u ∈ M \ Pgrow such that 2 2 σv|P < σv|P . grow ∪{u},pool grow ,pool While Lemma 10 gives stf-pool Pgrow ∪ {u} < stf-pool Pgrow , v v which contradicts the grow phase. Hence, M ⊆ Pgrow . We next show inductively that the shrink phase preserves M while removing all variables outside M . Initially, M ⊆ Pgrow . Suppose that the current parent set P during the shrink phase satisfies M ⊆ P . Since M is a global minimizer of the pooled residual-risk criterion and adding nodes cannot increase the optimal population residual risk, then we have 2 2 σv|P,pool = σv|M,pool .

If M ⊊ P , take any u ∈ P \ M . Since M ⊆ P \ {u} ⊆ P , the same monotonicity and optimality argument gives 2 2 2 σv|P \{u},pool = σv|P,pool = σv|M,pool .

By Lemma 10, deleting u decreases the empirical local score: stf-pool (P \ {u}) < stf-pool (P ). v v Thus, every strict superset of M admits a score-decreasing deletion of an extraneous parent. It remains to show that no score-decreasing deletion removes a member of M . Let u ∈ M . Since M ⊆ P , we have M \ (P \ {u}) = {u}. Suppose, for contradiction, that P \ {u} attains the same minimum population residual risk as M . Since it does not contain M , the no-local-trap assumption implies that adding u strictly reduces its population residual risk: 2 2 2 σv|P,pool = σv|(P \{u})∪{u},pool < σv|P \{u},pool .

which contradicts the assumed equality with the minimum risk attained by M . Therefore, 2 2 σv|P \{u},pool > σv|P,pool .

Applying the first part of Lemma 10, with Q = P and with P there replaced by P \ {u}, gives stf-pool (P ) < stf-pool (P \ {u}). v v Hence deleting a node in M does not improve the score. It follows inductively that every accepted shrink move removes a node in Pgrow \ M , and the shrink phase terminates at M . Therefore, one complete grow–shrink cycle returns Mvpool (Z) in the limit N → ∞. Since the argument does not depend on the initial parent set P0 ⊆ Z, target-filtered grow–shrink returns Mvpool (Z) with probability tending to one from any warm start.

I-FLOP: Fast Learning of Order and Parents from Interventional Data

A.4. Proof of Corollary 6 Proof Fix an order τ . For every node v ∈ V , let Zv (τ ) = Predτ (v) be the set of variables preceding v in τ . Then applying Theorem 5 to each node v gives Pr Pbv,N (τ ) = Mvpool (Zv (τ )) −→ 1. Since V is finite, these events hold simultaneously for all nodes with probability tending to one. Therefore, Pr Pbv,N (τ ) = Mvpool (Zv (τ )) for all v ∈ V −→ 1. b N (τ ) returned under the fixed order τ has parent sets By construction, the DAG G PaGbN (τ ) (v) = Pbv,N (τ ),

v ∈ V.

It follows that PaGbN (τ ) (v) = Mvpool (Zv (τ )),

v ∈ V,

simultaneously with probability tending to one. A.5. Proof of Lemma 7 Proof The local centered likelihood contribution of the original mechanism of v over its valid environments is X 2 Nv 1 X e (e) e (e) βv . ℓN (P ; βv , σv2 ) = ℓv,e (P ; βv , σv2 ) = − log(2πσv2 ) − 2 Xv − X P 2 2σv 2 e∈Ev

e∈Ev

For fixed P , maximizing over βv yields the pooled least-squares estimator computed from the valid environments of v. Maximizing further over σv2 gives 2 1 X e (e) 2 e (e) βbv|P . Xv − X σ bv|P = P Nv 2 e∈Ev

Therefore, − sup ℓN (P ; βv , σv2 ) = βv ,σv2

Nv 2 log(2π) + 1 + log σ bv|P . 2

Since ℓN (G, θG ) is node-wise decomposable, 1X X Nv 2 − sup ℓN (G, θG ) = Nv log(2π). 1 + log σ bv|Pa + G (v) 2 2 θG v∈V

v∈V

Adding the complexity penalty in the target-filtered pooled score yields X 1 Stf-pool (G; D, I) = − sup ℓN (G, θG ) + log N (| PaG (v)| + 1) + CN , 2 θG v∈V 1P where CN = − 2 v∈V Nv log(2π) does not depend on G. Thus, Stf-pool is a decomposable criterion for the valid-environment structural-mechanism likelihood after profiling out environment-specific intercepts. By Assumption 1, the likelihood contributions associated with directly targeted replacement mechanisms are independent of the candidate parent sets of the original mechanisms and therefore do not affect graph comparison.

Chen Markham

A.6. Proof of Theorem 8 Proof By Lemma 7, the target-filtered pooled score agrees, up to terms independent of the candidate DAG, with the negative maximized valid-environment structural-mechanism likelihood after profiling out environment-specific intercepts. Under known hard interventions, the omitted likelihood contributions of directly targeted replacement mechanisms are graph-independent nuisance terms. Hence, L∞ (G) is the population fit criterion induced by the graph-dependent part of the known-target interventional likelihood. The true DAG G⋆ generates the observational and interventional distribution family under the known target family I. Therefore, it attains the minimum value of L∞ , and G⋆ ∈ A. Now let H ∈ A. Since H attains the same minimum value of L∞ as G⋆ , the maximized population graph-dependent likelihood under H has no loss relative to the true known-target interventional distribution family. Under the regular finite-dimensional Gaussian model, a misspecified candidate model would incur a strictly positive population Kullback–Leibler loss. Hence, the true observational distribution is Markov with respect to H, and, for every intervention target Ie ∈ I, the corresponding interventional distribution is Markov with respect to the intervention graph H Ie . Thus, every H ∈ A is a causal independence map of the true known-target interventional distribution family. Fix any G ∈ Amin . Since G⋆ ∈ A and Amin retains the sparsest members of A, we have |EG | ≤ |EG⋆ |. By the interventional faithfulness assumption, Hauser and Bühlmann (2015, Remark 1) show that the minimum causal independence maps of the true known-target interventional distribution family are exactly the DAGs in the I-Markov equivalence class of G⋆ . In particular, G⋆ is a minimum causal independence map, so no causal independence map of this distribution family has fewer edges than G⋆ . Since G ∈ Amin ⊆ A, the preceding argument shows that G is itself a causal independence map. Therefore, |EG | ≥ |EG⋆ |. Combining the two inequalities yields |EG | = |EG⋆ |. Hence, G is a minimum causal independence map of the true known-target interventional distribution family and G ∈ GI⋆ . Since G ∈ Amin was arbitrary, we conclude that Amin ⊆ GI⋆ .

A.7. Proof of Corollary 9 Proof Let b N ∈ arg min Stf-pool (G; D, I) G G∈GDAG

be an exact global minimizer of the empirical target-filtered pooled score. Define its empirical fit component by X Nv 2 log σ bv|Pa . RN (G) = G (v) 2 v∈V

I-FLOP: Fast Learning of Order and Parents from Interventional Data

Up to terms independent of G, Stf-pool (G; D, I) = RN (G) +

X 1 log N (| PaG (v)| + 1) . 2 v∈V

For every fixed DAG G, consistency of the valid-environment pooled residual variances and the assumption on the asymptotic environment sample proportions imply 1 p RN (G) → − L∞ (G). N Since V is fixed, GDAG is finite, and this convergence holds simultaneously over all candidate DAGs. Consider any G ∈ / A and any G0 ∈ A. By the definition of A, we have L∞ (G) > L∞ (G0 ). Therefore, RN (G) − RN (G0 ) = N (L∞ (G) − L∞ (G0 )) + op (N ). The difference in complexity penalties is of order log N , whereas the fit difference is of order N . Hence, p Stf-pool (G; D, I) − Stf-pool (G0 ; D, I) → − +∞. b N ∈ A −→ 1. Since the DAG space is finite, it follows that Pr G It remains to distinguish among the members of A. Let G ∈ A \ Amin and G0 ∈ Amin . Since G, G0 ∈ A, both attain the same value of L∞ and, by the proof of Theorem 8, are causal independence maps of the true known-target interventional distribution family. Under Assumptions 1, 3, and 4, the difference between their maximized valid-environment structural-mechanism likelihood contributions is Op (1) for the fixed-dimensional Gaussian conditional-mechanism models considered here. Equivalently, RN (G) − RN (G0 ) = Op (1). Since G ∈ / Amin and G0 ∈ Amin , we have |EG | > |EG0 |. Moreover, X (| PaG (v)| + 1) = |EG | + |V |. v∈V

Hence, the complexity penalty difference between G and G0 is 21 log N (|EG | − |EG0 |), which diverges to +∞. Therefore, Pr(Stf-pool (G; D, I) > Stf-pool (G0 ; D, I)) −→ 1. b N ∈ Amin −→ 1. By Theorem 8, we also Again using finiteness of GDAG , we obtain Pr G b N ∈ G ⋆ −→ 1. have Amin ⊆ GI⋆ . Consequently, Pr G I If the known target family fully identifies the true DAG, so that GI⋆ = {G⋆ }, then b N = G⋆ −→ 1. Pr G

Chen Markham

Appendix B. Detailed Algorithmic Specifications Algorithms 2 and 3 detail the reinsertion and warm-start parent-set updates invoked in Algorithm 1. For brevity, let sv := stf-pool , v

P := (Px )x∈V ,

ℓ := (ℓx )x∈V ,

where Px is the current parent set of node x and ℓx = sx (Px ) is its current local score. We write X Q(ℓ) := ℓx x∈V

for the total score. For a prefix change of node x, +x means that x becomes admissible as a parent, whereas −x means that x ceases to be admissible as a parent. To describe the prefix changes induced by an adjacent swap, define ( +x, d = +1, ∆(d, x) := d ∈ {−1, +1}. −x, d = −1, Algorithm 2: Target-filtered reinsertion of node v Input: Current order τ , parent sets P , local scores ℓ, node v b Output: Best reinsertion state (b τ , Pb, ℓ) // Initialize the best state with the current state b ← (τ, P, ℓ) (b τ , Pb, ℓ) // Explore moving v to the left and to the right foreach d ∈ {−1, +1} do (τ ′ , P ′ , ℓ′ ) ← (τ, P, ℓ) j ← positionτ ′ (v) while 1 ≤ j + d ≤ p do ′ u ← τj+d // Swap v with the adjacent node u ′ ) ← (τ ′ , τ ′ ) (τj′ , τj+d j+d j // Update only the two parent sets affected by the swap (Pv′ , ℓ′v ) ← GrowShrinktf v, Predτ ′ (v), Pv′ , ℓ′v , ∆(d, u) (Pu′ , ℓ′u ) ← GrowShrinktf u, Predτ ′ (u), Pu′ , ℓ′u , ∆(−d, v) j ←j+d // Retain the best reinsertion position found so far b then if Q(ℓ′ ) < Q(ℓ) b ← (τ ′ , P ′ , ℓ′ ) (b τ , Pb, ℓ) end end end b return (b τ , Pb, ℓ) 22

I-FLOP: Fast Learning of Order and Parents from Interventional Data

Algorithm 3: Target-filtered warm-start grow–shrink Function GrowShrinktf (v, Z, P0 , ℓ0 , δ): if δ = +u then P ← P0 ∪ {u}; ℓ ← sv (P ) if ℓ ≥ ℓ0 then return P0 , ℓ0 end else if δ = −u then if u ∈ / P0 then return P0 , ℓ0 end P ← P0 \ {u}; ℓ ← sv (P ) else P ← ∅; ℓ ← sv (∅) end foreach d ∈ {+, −} do repeat P old ← P ( Z \ P old , d = +, U← P old , d = −; foreach u ∈ U do if d = + then P ′ ← P ∪ {u} end else if u ∈ / P then continue end P ′ ← P \ {u} end ℓ′ ← sv (P ′ ) if ℓ′ < ℓ then (P, ℓ) ← (P ′ , ℓ′ ) end end until P = P old end return P, ℓ The δ = +u short-circuit (returning P0 , ℓ0 immediately when adding u does not improve the score) is justified whenever P0 is itself already a grow–shrink fixed point over Z \ {u}, which holds at every call site in Algorithm 2: by Proposition 4, if no single admissible addition to such a P0 improves the score, no further local move can either, so the full grow–shrink procedure would return P0 unchanged regardless. 23

Chen Markham

Figure 4: Top: Scaling with the number of samples per environment, with p = 40, expected degree 3, and sample size varied from 102 to 106 in each environment. Bottom: Scaling with graph size, with expected degree 3, 1000 samples per environment, and p varied from 20 to 100. Each setting contains one observational environment and five known-target singleton perfect hard-intervention environments. Runtime is shown on a logarithmic scale, and accuracy is measured by I-CPDAG SHD.

Appendix C. Computational Scaling with Graph Size and Sample Size Figure 4 reports additional scaling results on sparse ER DAGs. We examine two complementary regimes: increasing the number of samples per environment while holding the graph regime fixed, and increasing the number of variables while fixing the sample size, expected degree, and intervention design. Scaling with sample size. As the number of samples per environment increases, I-FLOP substantially reduces its I-CPDAG SHD, approaching near-perfect recovery. GIES exhibits a comparable reduction in structural error while maintaining nearly constant runtime in this experiment. By contrast, the runtime of the FLOP-based methods increases markedly with sample size, particularly beyond 104 samples per environment. FLOP-envwise does not match the accuracy of I-FLOP at large sample sizes, indicating that environment24

I-FLOP: Fast Learning of Order and Parents from Interventional Data

wise estimation alone is insufficient without target-based filtering. I-GSP improves with additional samples but retains noticeably higher structural error than I-FLOP and GIES. Scaling with graph size. When the number of variables increases, I-FLOP maintains the lowest structural error among the evaluated methods while remaining computationally inexpensive relative to GIES. Its runtime grows moderately from p = 20 to p = 100, whereas GIES becomes substantially more expensive, reaching several hundred seconds at the largest graph size. The accuracy difference is also pronounced for larger graphs: I-FLOP retains comparatively low I-CPDAG SHD, while the errors of GIES, FLOP-obs, and especially I-GSP increase considerably. FLOP-envwise remains computationally close to I-FLOP but incurs consistently higher recovery error, again supporting the contribution of target-based environment filtering.

Figure 5: Runtime of GIES is decomposed into covariance/scatter precomputation, graph search, and end-to-end wrapper execution. All the settings are the same as Figure 4.

Explanation of the runtime of GIES. The nearly flat runtime curve of GIES in the sample-size experiment is explained by the timing decomposition in Figure 5. GIES performs its graph search using covariance computation. As the number of samples per environment increases from 102 to 106 , the cost of computing these statistics rises from the millisecond scale to below one second. However, the search stage consistently requires approximately 20 seconds and therefore dominates the end-to-end runtime throughout the evaluated range. End-to-end wrapper execution, covering input conversion into the GIES implementation’s required format and output conversion back into a standardized DAG representation for evaluation, contributes only a small, roughly constant overhead relative to the search stage. Consequently, increasing the sample size has little visible effect on total GIES runtime in Figure 4. The small non-monotonic variation at lower sample sizes is 25

Chen Markham

Figure 6: Left: Mean optimized total BIC score across 50 seeds as the ILS budget increases from 0 to 50. Right: Mean total SHD and its skeleton and orientation components across the same seeds. The shaded regions denote the interquartile range across seeds.

attributable to changes in the estimated scores, which will alter the greedy search trajectory and its number of steps, rather than to a systematic reduction in computational cost.

Appendix D. Diagnostic ILS Trace on the Real Data Experiment To investigate the observed mismatch between score improvement and graph recovery as the ILS budget increases, we conducted an additional diagnostic analysis on the Chamber dataset from the real-data experiment (Section 4.3). This analysis follows exactly the same experimental setting as the original real-data evaluation and examines how the optimized total BIC score and the corresponding SHD change as the ILS budget increases. As shown in Figure 6 (left), the mean total score decreases rapidly during the first several ILS steps and then gradually stabilizes. From ILS 0 to ILS 50, the mean total BIC score decreases by approximately 21.55 units, and 49 of the 50 seeds achieve a lower score at ILS 50 than at ILS 0. However, this improvement in the optimized total score does not translate into better graph recovery. Figure 6 (right) shows that the mean total SHD increases from 29.42 at ILS 0 to 30.12 at ILS 50. This increase is almost entirely attributable to the orientation component: the mean skeleton SHD remains essentially unchanged, decreasing slightly from 28.14 to 28.12, whereas the mean orientation SHD increases from 1.28 to 2.00. However, the aggregate results mask substantial heterogeneity across seeds. As shown in Figure 7, comparing ILS 0 with ILS 50, total SHD worsens for 27 seeds, improves for 14 seeds, and remains unchanged for 9 seeds. Skeleton SHD increases for 20 seeds and decreases for 20 seeds, whereas orientation SHD increases for 20 seeds and decreases for none. Although skeleton changes constitute the larger component-wise changes for 33 seeds, they occur in both directions and therefore largely cancel in the aggregate mean. By contrast, orientation changes are one-sided and consequently account for the net increase in mean total SHD. Thus, in the real data experiment, increasing the ILS budget generally yields lower-scoring solutions, but these solutions are not necessarily closer to the reference graph under SHD. This score–recovery mismatch is plausibly related to a mismatch between the assumptions of the I-FLOP algorithm and the data-generating process of the Causal Chamber 26

I-FLOP: Fast Learning of Order and Parents from Interventional Data

Figure 7: Seed-level changes in SHD on the real data dataset when increasing the ILS budget from 0 to 50. The horizontal axis reports the change in SHD, computed as the value at ILS 50 minus that at ILS 0. Seeds are grouped according to the change in total SHD: improved (n = 14), unchanged (n = 9), and worsened (n = 27). The black line denotes the change in total SHD, while the orange and purple bars denote the changes in skeleton SHD and orientation SHD, respectively.

dataset. In the synthetic experiments, the data are generated from a linear-Gaussian SEM, and interventions are implemented as Gaussian shifts on target variables while preserving the same linear mechanisms for the remaining variables. However, Causal Chamber is obtained from a real optical system rather than from this synthetic SEM. Its intervention 27

Chen Markham

environments are created by changing the operating ranges of physical actuators to bounded intervention-specific ranges. These interventions therefore alter the support and operating regime of the intervened variables, rather than merely applying Gaussian mean shifts within the same assumed model. These differences affect the interpretation of the optimized score. The I-FLOP score evaluates candidate graphs through linear-Gaussian regressions pooled across the relevant observational and interventional environments. Consequently, a lower total BIC score on the real data experiment represents improved fit under the assumed linear-Gaussian scoring approximation, rather than direct evidence of improved recovery of the physical reference graph. The ILS results suggest that increasing the search budget improves the optimized score, while the corresponding solutions retain essentially unchanged skeleton accuracy and accumulate additional orientation errors.

I-FLOP: Fast Learning of Order and Parents from Interventional Data

Analyze & explore related research in Conceptio

Related documents