Embracing Discrete Search: A Reasonable Approach to Causal Structure Learning
Authors: Marcel Wienöbst
Leonard Henckel
Sebastian Weichwald
Venue: ICLR 2026 Poster (2026)
Keywords: Causal Discovery
Bayesian Networks
DAGs
Structure Learning
Primary area: causal reasoning
TL;DR: Discrete graph search is a reasonable structure learning approach enabling accurate causal discovery; the size of the graph space was not the barrier.

Abstract: We present FLOP (Fast Learning of Order and Parents), a score-based causal discovery algorithm for linear models. It pairs fast parent selection with iterative Cholesky-based score updates, cutting run-times over prior algorithms. This makes it feasible to fully embrace discrete search, enabling iterated local search with principled order initialization to find graphs with scores at or close to the global optimum. The resulting structures are highly accurate across benchmarks, with near-perfect recovery in standard settings. This performance calls for revisiting discrete search over graphs as a reasonable approach to causal discovery.

Published as a conference paper at ICLR 2026
EMBRACING DISCRETE SEARCH: A REASONABLE APPROACH TO CAUSAL STRUCTURE LEARNING
Marcel Wienöbst', Leonard Henckel?, Sebastian Weichwald
University of Lübeck, ?University College Dublin, University of Copenhagen
ABSTRACT
We present FL.OP (Fast Learning of Order and Parents), score-based causal discovery algorithm for linear models. It pairs fast parent selection with itera- tive Cholesky-based score updates, cutting run-times over prior algorithms. This makes it feasible to fully embrace discrete search, enabling iterated local search with principled order initialization to find graphs with scores at or close the global optimum. The resulting structures are highly accurate across benchmarks, with near-perfect recovery in standard settings. This performance calls for revis- iting discrete search over graphs reasonable approach causal discovery.
INTRODUCTION
Learning about the directed acyclic graph (DAG) underlying system's eetta-genrrating process from observational data under causal sufficiency is a fundamental causal discovery task (Pearl, 2009). Score-based algorithms address this task by assigning penalized likelihood scores to each DAG and seeking graphs whose scores are optimal. Identifiability theory asks when such score- optimal graphs identify the target graph (or its equivalence class) in the infinite-sample limit, with various results under different assumptions and scores (Chickering, 2002; Nandy et al., 2018).
Exact algorithms, that are guaranteed to find a score-optimal graph, have exponential run-time and are feasible up to roughly 30 variables (Koivisto & Sood, 2004; Silander & Myllymäki, 2006). For larger graphs, local search must be employed, which evaluates neighbouring graphs find graphs with better scores; canonical moves for this hill climbing are single edge insertions, deletions, or reversals (Heckerman et al., 1995). the sample limit, greedy discrete search with neighbour- hood notion that respects score equivalence provably finds graph with optimal score (Chickering, 2002). In finite samples, scores are incxact and local search may get stuck in local optima or, as we demonstrate, even find graphs with better scores than the true graph. Finite-sample performance is a practical challenge, despite the mature identifiability theory and asymptotic guarantees.
Continuous optimization methods have emerged a popular alternative. For example, NOTEARS encodes acyclicity as a smooth constraint and optimizes a surrogate objective (Zheng et al., 2018), with many follow-ups (Bello et al., 2022; Rolland et al., 2022). Their supposed advantages have been questioned empirically and conceptually (Reisach et al., 2021; 2023; Ng et al., 2024). Further, NP-hardness results often cited to dismiss discrete search do not apply the commonly consid- ered discovery settings: The standard hardness constructions rely on data-generating processes that involve unobserved variables and cannot be represented by a DAG over only the observed vari- ables (Chickering, 1996; Chickering et al., 2004). When the distribution is representable by a sparse DAG, discrete procedures asymptotically recover the target graph with polynomially many indepen- dence tests or score evaluations (Claassen et al., 2013; Chickering & Meek, 2015).
One of the core issues of score-based methods in practice are sanite-sample induced local op- tima (Nielsen et al., 2003). Hence, the best-performing heuristics in benchmarks (Rios et al., 2025) are either able to escape local optima, for example through simulated annealing (Kuipers et al., 2022), or realize larger neighborhoods (Pisinger & Ropke, 2018), such as recent order-based meth- ods (Lam et al., 2022; Andrews et al., 2023) with effective reinsertion moves rather than only swap- ping neighboring nodes (Teyssier & Koller, 2005; Scanagatta et al., 2015). This helps explain the strong performance of the order-based BOSS algorithm (Andrews et al., 2023) and more recent
"Correspondence: swei chwal 1 dēmath. ku.

---

Published as a conference paper at ICLR 2026
ER, 50 nodes, average degree 8
ER, 50 nodes, average degree 8
SHD
AID
4
1500
DAGMA
400
300
GES
GES
1000
PC
200
PC
DAGMA
500
100
FLOP
FLOP
BOSS
BOSS
FLOP20 20
FLOPPo
FLOP
100
FLOP
0
0
2
3
4
5
2
3
4
5
Runtime in seconds
Runtime in seconds
Figure I: Run-time plotted against Structural Hamming Distance (left) and Ancestor Adjustment Identification Distance (Henckel al., 2024) (right) between the CPDAGs learned on linear ANM data and the target CPDAG corresponding to the underlying Erdōs-Renyi generated DAG with 50 nodes, average degree 8 and 1000 samples drawn. Every point corresponds one of 50 random instances; diamonds indicate averages. FI.OP variants differ in the number of ILS restarts to escape local optima. The fraction of instances with exact CPDAG recovery is 40% for BOSS and FLOPPo and 60% for FLOP20 and FLOP and zero for the remaining algorithms.
order-based local searches (Li et al., 2025) on common causal discovery benchmarks. Continu- ous relaxations also alter the search space traversal, yet they have not matched this performance and introduce additional challenges, for example optimization complexity, convergence issues, edge thresholding, and having resort to surrogate objectives.
Contributions. We introduce FL.OP (Fast Learning of Order and Parents), a score-based struc- ture learning algorithm for linear additive noise models that fully embraces discrete search, and ready-to-use from Python via pip install flopsearch. The FLOP algorithm adopts reinsertion- and order- based explo- ration of DAGs (Andrews et al., 2023) and adds four components that enable aggressive search for graphs with optimal BIC score. First, we simplify parent selection by re- re-initializing from the par- ent sets learned for the previous order, which reduces compute and memory cost without degrading performance (Section 3. Second, we accelerate score computations for the linear Gaussian BIC via efficient iterative updates of Cholesky factorizations, which amortize cost across local moves (Section 3.2). Third, develop principled initialization that, compared with a random initial order, reduces local parent selection failures on mncestor-descendane pairs that are far apart and only weakly dependent (Section 4.1). Fourth, the computational gains allow
offer Rust implementation at github.com//aaaa comm/CasalDisccofoossaaaa
employ an iterated local search (ILS) metaheuristic to escape local optima (Section 4.2). On standard benchmarks, the order-based methods BOSS and FLOP achieve strong accuracy at favorable run-time (see Figure 1). FLOP's run-time advantage can be translated into higher accuracy by extending the ILS budget, with FLOP; & denoting k iterations. By treating compute budget as a hyperparameter, our work highlights the link between run time and finite-sample accuracy in causal discovery.
PRELIMINARIES
We consider the problem of learning about the acyelic graph structure underlying linear additive noise models (ANMs) from observational data. For causal DAG G (V.E) with node set V = fl.....ph }} and edge set E t [i+jj i.j E V). a linear ANM is defined by a weight
matrix W E RPxP with Wa, # 0 = i-je E E and vector N [N1, [N., Nº1" of jointly in- dependent, real-valued zero-mean noise variables with finite fourth moment; the observed variables
x [X1,.., .., x,) are then defined by X = WT X + N.
2

---

Published as a conference paper at ICLR 2026
Algorithm 1: Reinsertion based local search proposed by Andrews et al. (2023).
input : Data set D over p variables.
output: A CPDAG G.
7 initial order of (I.....P) // Below, refers to the i-th element of
/* Per ly contain the parents and the local score of node e foreachie pp do (Pr,)n))) growShrink (7i Thi-1.
repeat
(ruldd gald) = (T for V E puli do
(T, P. 8) reinsert (T, P.P. D) // Optimal reinsertion for v.
end
until sum (f) > sum (lou)
return CPDAG of the DAG defined by parent sets Pl...., Pp
We denote the parents of v E V. that is all u E V such that + v E E, by Pa(u). Every at least one linear order т, also called topological order, of the nodes such that u + e E E' implies u coming before v in т. A DAG is called Markovian to a probability distribution if every variable is independent of its non-descendants (all nodes not reachable from it with a directed path) given its parents. We denote conditional independence by I.
DAG G = (V, E) can be associated
As score optimize. we choose the Bayesian Information Criterion (Schwarz, 1978). BIC for short, which is the common choice in score-based structure learning for ANMs and is asymp- totically consistent under the faithfulness assumption (Koller & Friedman, 2009), which we em- ploy throughout this paper. For DAG G (V, E)) and data set D containing n observa- tions, it is defined k In(n)= 2 (l((L with & being the number of parameters and É. the max- imized likelihood for the given DAG and data. The score can be decomposed into local scores BICD(G, X) Evev (X(,, Xpa(w)) and, for linear models with Gaussian noise, each local score is given by Xp((n) = n log( (ar rameter.' The BIC is score-equivalent Its value is the same for Markov equivalent DAGs (Verma & Pearl, 1990) that imply the same conditional independencies. In fact, with observational data under causal sufficiency and without additional assumptions, only such equivalence class of DAGs is identifiable. Throughout, our target object is therefore the equivalence class of the underlying DAG, represented by a completed partially directed acyclic graph (Andersson et al., 1997). Methods that internally optimize over DAGS, such as continuous relaxations, or rely on extra assumptions to iden- tify unique DAG, such as non- non- Gaussian linear models (Shimizu al., 2011). or noise-variance conditions (Park, 2020), need to be evaluated via the corresponding CPDAG for fair comparison; evaluating a single DAG as if identified is arbitrary under our assumptions and can be misleading.
+A(m))))Ann(())) [Pa(v)] with À being penalty pa-
We employ principles from the BOSS algorithm (Andrews et al., 2023) to optimize the BIC score over DAGs: We traverse the space of topological orders of DAGs, iteratively moving to orders that result in better scoring DAGs when selecting parents accordingly. Candidate orders are generated by taking a variable and reinserting procedure (Margaritis, 2003) to construct parent set for each variable from its prefix, the variables preceding it in т, and score the resulting DAG. Algorithm shows the BOSS reinsertion strategy, which we build on for the improved search in FLOP.
another position. Given an order т, we use the grow-shrink
SCALING UP ORDER -BASED SEARCH
This section presents two speedups for order-based local search with reinsertion moves, which yield significantly faster nu-timms than the grow-shrink trees used in BOSS. FLOP still uses grow-shrink to obtain DAGs from orders, but in a way that exploits the local, iterative search moves and scoring.
Here, one needs to assume that the empirical covariance matrix is non-degenerate and numerically well-conditioned, excluding cases such as vero-variance noise, high-dimensional settings with p. or vari- ance hlow-up severe enough cause numerical instability.
3

---

Published as a conference paper at ICLR 2026
Algorithm 2: Find the best-scoring reinsertion of node e in order given data D.
function re inser (T, P. (, v, D) // P stores parents, local scores.
i: position of u i
P,) = P.
foreach j E [1+ ,....,
1.. p) do // Test reinsertions at later positions. // Swap T/-1 one position to the right. // Compute parents for changed prefixes
(ī)-1, = (T1)
TJ-1 and Ty. (L,,, f=)==rowShrink growShrink (T,, 71j-1, D, Pc, b +Tz-1)
(6, 2) = growShrink (Tj-b T1j-2, D, Pr, - 1 -T))
if sum (E) sum (é) then (f,P,0)= (T, end foreach j E {i I..... .....,11 1 do... // Analogous for earlier positions.
10
11 end
12 return 1, P
3.1 STARTING GROW-SHRINK FROM THE PREVIOUS PARENT SET
During the scoring of node-reinsertioss each node's 's candidate parent set, that is, the nodes coming before it in the order, changes by at most one node being inserted or deleted from its prefix. Consider Algorithm that finds the best reinsertion for node v currently at position in order T. The possible reinsertions of can be efficiently evaluated by performing sweep from position to the right (and also to the left; analogous code omitted), moving it to position i 1, + 2, and so on, by swapping it rightward. At each step, the prefix of node increases by exactly one element, while the prefixes of nodes originally at positions i + 1, i + 2, and so on, lose exactly one node, namely v.
Instead of running grow-shrink from the empty set at every step as in BOSS, FLOP initializes grow- shrink with the previous parent set, that is, it continues from the result of grow-shrink for the previous prefix, now with additional one fewer node. The idea behind this strategy is that the parent set typically changes little when the prefix changes by just one node, and so this warm start makes parent selection far cheaper. Our implementation given in Algorithm 3. Moreover, our grow-shrink does not insist on inserting removing the single best parent with largest score improvement, making it non-greedy: it adds or removes any parent that improves the score, even if not maximally so. This eliminates the need for complicated grow-shrink tree caching as used in BOSS.
We show that the modified grow-shrink with warm start learns the restricted Markov boundary of a node with respect to a set Z (Lam et al., 2022). This yields theoretical guarantees that the DAG learned by FLOP the sparsest Markovian one for the considered order (Raskutti & Uhler, 2018).
Definition 3.1. Let P be a distribution over XI...., Xp. The restricted Markov boundary of X, relative to a set Z C (X1...., .Xp) ((X.). denoted by Mv, Z), is defined as set of nodes M C Z such that a) X,, 1 (Z M) M and b) there exists no M' C M such that X,, (Z \ M') M'.
Under mild assumptions, the Markov boundary is unique (Verma & Pearl, 1988). As in GRaSP (Lam et al., 2022) and BOSS (Andrews et al., 2023), we learn it using BIC score improvements in place of conditional independence tests (Margaritis, 2003). In the sample limit, the local BIC score (LX,,XPater U { ) is smaller than E(X,, Xpa(w) if, and only if, X, L X,, XPa(e) (Koller & Friedman, 2009). We show that this asymptotic guarantee also carries over to the modified grow- shrink algorithm that starts from an arbitrary initial parent set instead of the empty set.
Lemma 3.2. Let data set D consist of n i.i.d. observations of a probability distribution represented Then, the large sample limit of n. grow-shrink
by a Bayesian network over variables X1....,X; .p. finds the restricted Markov boundary of node v relative to a set Z C [x,,..., X.) (Xo} when sturted with any initial set P C Z.
Proof. This follows directly from the proof of correctness of the grow-shrink algorithm in (Margari- tis, 2003) and its generalization to restricted Markov boundaries in (Lam et al., 2022). Assume that at the end of the grow-phase, the current set of parents is P Parow Thus, it holds that X,, 1 X,, Paow for all X,, Z \ P Pgrw,, or rephrased X, 1. 7 Pgrow Pgrow. However, this would violate the
4

---

Published as a conference paper at ICLR 2026
Algorithm 3: Non-greedy grow-shrink with the option to start from previous parent set prev.
function growShrink (u, Z, D. Prev ppee 5)
i > then // 6 the node added to the candidate parents Z Prew == P prev U {8}
fnew = locall Soccalccores Score (u, Prew. D) // score with 8 added to parents
if Lnew feev then return pper. Ppee // return if no improvement
else if 5 < 0 tten/// [6] is the node removed from the candidates 7 if |6| Prrev then return (prevs Prev return if |6] was no parent
P new prev \ {[8}}
new = localScore (u, Phew. D)
else If no is provided, run from scratch. 10
(Pr (new) = (0, lol.llllloooo (u, 0, D)
12
end
grow (u, Pnews Cnews Z, D) shrink (u, Pnews knews Z, D)
14
15
end
16
function grow (u, P, D)
repeat
18
foreach v E P
if feew localScooe (u, PU PU (v}. D) < then (t (Loww PU {b})
20
end
21
until P is unchanged
22 end
23 function shrink (u, P. D) ... Analogous to grow, thus omitted.
uniqueness of the Markov boundary Mr,, 7) if it is not subset of P grow This argument does not depend on the initial set P. The correctness of the shrink-phase is unchanged, too.
In addition to the warm start, we implement another optimization. We pass the node v that we are either inserting to (coded +-) or removing from (coded as -v) the prefixes into grow-shrink. If has been removed and was not part of Ppev. we immediately return Ppevv node v has been inserted to the prefix and does not increase the score when added to Prev. we again immediately return P Ppev- We show that these modifications preserve the guarantees above. In the sample limit, FLOP returns a Markovian DAG, that is, one that induces additional conditional independencies.
Theorem 3.3. Let data set D consist of n i.i.d. observations of probability distribution represented by a Bayesian network over X1....,X, Xp. In the sample limit of n, the CPDAG returned by FLOP is Markovian P.
Proof. This statement holds assuming that the grow-shrink procedure FLOP finds Markovian graph for each scored order. As the grow-shrink routines depend on the previous runs, we prove this by induction. Initially, standard grow-shrink run for the starting order (line 2 of Algorithm 1), which yields parent sets corresponding to its sparsest Markovian DAG (Raskutti & Uhler, 2018; Lam al., 2022). Assume that the parent sets for the previous order have this property. By Lemma 3.2, the modified grow-shrink, if run fully, finds the restricted Markov boundary with respect to the prefix and thus yields parent sets of the sparsest Markovian DAG. It remains to show that the two early breaks in lines 5 and of Algorithm are correct, where the grow-shrink is not run.
removed node was not part of previous Markov boundary, clearly the Markov boundary remains unchanged for the reduced prefix, as neither the grow nor the shrink phase would add or remove node. added node 8 does not increase the score for the enlarged prefix, this is the case, too. Thus, by Lemma 3.2, the DAG learned by FLOP is Markovian and the statement follows by the fact that the CPDAG of such a DAG is returned.
We remark that further modifications to FLOP in the subsequent sections do not change this result. As with BOSS, can make FL.O asymptotically consistent, provably yielding the true graph in
5

---

Published as a conference paper at ICLR 2026
Time in
ER, average degree 16, 1000 samples
Pre-FLOP (grow-shrink from previous parents)
1000
Pre-FLOP (grow-shrink from empty set)
BOSS
FLOP
500
(grow-shrink from previous parents and Cholesky updates)
50
100
150
200
250
300
350
400
450
500
Number of nodes
Run-time in seconds, averaged over 50 repetitions with standard-deviation enror bars, for ER graphs with average degree 16, 1000 samples, and (50, 100, 150...., 500} nodes.
Figure 2:
the sample limit under the faithfulness assumption, by running the backwards phase of GES (Chick- ering, 2002) after termination of the local search. However, we refrain from this, since FLOP already reaches ofate-of-the-ar art finite-sample performance without it. More generally, any score-based dis- covery algorithm can be made asymptotically consistent by running another consistent algorithm, for example GES or PC, in parallel and returning whichever graph attains the better score.
3.2 DYNAMIC CHOLESKY UPDATES
We exploit the local structure of grow-shrink avoid recomputing the local score from scratch every step. Instead, we update the score from the previous parent set which at each step changes only by single insertion removal. This idea is generic and applies to any local search that adds or removes one edge at a time in a local search, not only to order-based methods.
In the multivariate normal setting, the local BIC
node u is f(X XPa())
log( (arr (Pa(w))) + In(n)(Pa(u)). Since the penalty depends only on the size of the parent set, the work computing or updating the local score is in the likelihood term, in the esti- mated conditional variance of X,. A direct way to compute this would be to invert submatrix of the sample covariance matrix, but numerically more stable and faster to avoid the matrix inver- sion in favor of using Cholesky factorizations. As shown in Appendix B, the bottom right entry of the Cholesky factor of the covariance submatrix corresponding to Xpa(w) and X., yields the square root of the conditional variance of X., given Xpa(u)-
Computing the Cholesky decomposition of l & matrix requires (1/3) floating- point operations. However, as discussed above, the submatrix which we -actesky- changes by adding or removing only single row and column (corresponding to the added removed parent node). We therefore update the Cholesky factor instead of recomputing it, using standard rank-one update and downdate routines (Gill et al., 1974; Golub & Van Loan, 2013). These updates require O(k2) floating-pont point operations, shaving off a factor & compared to a fresh Cholesky decomposition. This run-time improvement proves advantageous for larger and denser graphs. These Cholesky updates are applicable to other score-based causal discovery algorithms, for example GES or other hill- climbers. To our knowledge, this speedup has not been described in prior causal discovery work.
3.3 RUN-TIME COMPARISON
In Figure 2, we compare the run-time of FLOP, which includes the two run-time improvements de- scribed in this section, with two ablated versions, termed pre-FLOP in the plot, the first one using neither optimization (thus only differing from BOSS by using non-greedy grow-shrink) and the sec- ond one only using the grow-shrink started at the previous parent set. For reference, we provide the run-time of the BOSS implementation in the retrad software package. The comparison with Tetrad is not apples-to-appee since Tetrad is written in Java and multithreaded, while our code is Rust and tingle-threaded. The benchmark uses Erdös-Renyi graphs with average degree 16, ori-

---

Published as a conference paper at ICLR 2026
ented according to sefformml-randdon linear order (details the simulation setup are provided in Section 5). Each run has a 30 minute time limit.
Both the modificd grow-shrink and the Cholesky updates yield substantial run-time reductions. With both optimizations, FLOP is more than a factor 100 faster than BOSS for graphs with 100 nodes and scales to 500 nodes, whereas BOSS reaches the time limit for instances with 150 nodes. We note that accuracy of the discovered graphs is similar on these instances for both methods, both giving good, but not perfect results. In the following, we use the optimized FLOP and build these speedups to further improve the quality of the found graphs
4 IMPROVING THE ACCURACY OF ORDER-BASED SEARCH
This section presents two techniques to improve search accuracy. First, we replace random initial or- ders with a principled initial order construction putting strongly-correlated nodes next each other, which critical on directed paths in finite samples. Second, we use Iterated Local Search (ILS), which perturbs found solution and restarts the local search, trying to escape local optima through additional compute. With these techniques, FLOP attains state-of-the-ar accuracy in simulations.
4.1 INITIAL ORDER
Zp are challenging instances for order-based methods, which, to our knowledge, have not been previously discussed in this context before. On the left of Figure 3, compare different algorithms path graphs with 50 nodes. FL.OP with random initial order (FLOPE"d) and BOSS are in fact the worst-performing of all methods. A reason for this are far- apart ancestor-descendaal pairs with very weak marginal dependence, for which the grow-shrink procedure may fail to add edges, resulting in non-Markovian DAGs in finite samples. For example, if z; and 2j with i < appear first in the order, grow grow-shrink should, irrespective of the remaining order, make 2; parent of 2j for it to yield a Markovian DAG since z; and zj are marginally dependent. However, the dependence between I; and zj may be too small for grow-shrink to pick up on in finite samples.
Path graphs 21 - I2
As a remedy, we build the initial order so that strongly correlated nodes are adjacent, facilitating grow-shrink to find Markovian graph. To build the order, we start with the two most correlated nodes and append, at each step, the variable that can be best explained by variables already placed the order, that is, the one with the smallest residual variance when regressed onto the nodes in the order. We standardize the data beforehand efficiently, by iteratively constructing a Cholesky decomposition of the covariance matrix choosing the next node in the order according their residual variance (see Appendix B). On 50-node paths (Figure 3, left), FLOP with this initial order has an average SHD on-par with PC and GES, the best performing algorithms these instances.
avoid scale artefacts. We compute this order
4.2 ITERATED LOCAL SEARCH
Iterated Local Search (ILS) is classic metaheuristic in discrete optimization (Lourenço et al., 2018) that has been used in previous score-based search over DAGs and CPDAGs (Liu et al., 2023; Nazaret & Blei, 2024). It is a generic strategy that combines local search with perturbations to escape local optima: Run local search local optimum, perturb the best solution seen so far, then rerun local search starting from this perturbation; repeat. In principle, this procedure can be repeated indefinitely.
For FLOP, the first local search starts from the initial order constructed as described in the previous section. After that, the starting order for the next local search is obtained by perturbing the best- found order by k random swaps of two (not necessarily adjacent) elements. The idea being, that orders near local optima are better starting points than fully random ones. We set = Inp by default, which found yield robust results balancing moving far enough to escape while staying in promising basin.
On dense Erdös- Rendi graphs with 25 nodes and average degree of 16 (Figure 3, right), increasing the number of restarts of the local search (zero restarts amount to one local search, x restarts to perturbations and new local searches after that), consistently improves FLOP's accuracy. With

---

Published as a conference paper at ICLR 2026
Path, 50 nodes
ER, 25 nodes, average degree 16
SHD
SHD
FLOPPo
30
200
FL.OPO
BOSSa
150
20
BOSS
100
BOSS
DAGMA
FLOP
BOSS
FL.PPo
50
GES
FLOP
PC
0
FLOP
300
Exact
4
0
10
20
30
40
50
60
70
80
Runtime in seconds
Runtime in seconds
Figure 3: Run-time plotted against SHD on paths with 50 nodes for 1000 samples (left) and ER graphs with 25 nodes and average degree 16 for 50,000 samples (right). For the path graph. FLOP, finds the target graph in 72% of instances, PC in 32%, GES in 66% and the remaining algorithms in none; for the ER graphs, FLOP 20 does in 26% of cases, FLOP 100 in 50%, FLOP 500 in 56%, Exact in 58%, BOSS100 4% and the remaining algorithms none.
500 restarts, FLOP matches the exact score-based algorithm while having a faster run-time (the exact score-based algorithm implements the method by Silander & Myllymäki (2006) and multithreading). We compare against BOSS with full random restarts, that is x restarts mean x + independent runs of BOSS and returning the best-found solution. This is computationally heavy and yields substantially smaller gains than FLOP's ILS restarts.
ILS is an integral part of the FLOP algorithm. When calling FLOP, the user needs to specify either the number of restarts of the local search or time limit, and the solver runs ILS until the budget exhausted. This emphasizes the trade-off between run-time and accuracy inherent to scsee- causal discovery, but effectively ignored by the structure learning community with its focus on one- shot heuristic algorithms.
SIMULATIONS
We empirically compare FLOP to other causal discovery methods. For Figure 1, we generate Erdös- Renyi (ER) graphs with 50 nodes and average degree 8. We also consider scale-free (SF) graphs with density parameter i 4, generated by starting with a star graph of A: + nodes and adding further nodes by preferential attachment to & existing nodes, and DAGs from the bnlearn repository (Scutari, 2010), such as the Alarm network (Beinlich et al., 1989). We orient all graphs according to linear orders drawn uniformly random. For each graph, we generate 1000 samples from linear additive noise model with Gaussian noise (with mean 0 and variance uniformly drawn from [0.5, 2.0) and edge coefficients drawn uniformly from [-1, 0.255 U [0.25, 1]. Each setting is repeated for 50 random instances.
In addition FLOP and BOSS, we rum PC (Spirtes al., 2000) as classical constraint-based method, GES (Chickering, 2002) as traditional woorebbsed algorithm. and DAGMA (Bello et al., 2022) as gradient- continuous optimization method. For BOSS, PC, and GES, we rely on the implementation in Tetrad (Ramsey et al., 2018) through causal- -emd version 1.12.0, for DAGMA we use version 1.0 of the authors' implementation. The algorithms are run on a machine with 256GB of RAM and AMD Ryzen Threadripper 3970 CPU with 32 cores. We make no restrictions on the number of threads the implementations may use (FLOP only uses a single thread, whereas the other algorithms exploit multithreading) and report the wall-clock time of their execution. We use standard parameters in the literature, setting \BIc = 2 for the BIC-based algorithms (for a motivation of higher penalty parameter than prescribed by the standard BIC, Foygel & Drton, 2010),
8

---

Published as a conference paper at ICLR 2026
SF, 50 nodes, density parameter 4
Alarm network
SHD
SHD
300
50
GES
40
200
DAGMA
DAGMA
30
PC
FLOPo
100
20
FLOPP
BOSS
PC
BOSS
10
GES
FLOP2 20
FLOP
FLOP20
20
FL.OP
100
0
2
3
5
2
3
Runtime in seconds
Runtime in seconds
Figure 4: Run-time plotted against SHD on SF graphs (left) and the Alarm network, consisting of 37 nodes and 46 edges, (right), both for 1000 samples. For the SF graphs, FLOP20 finds the target CPDAG in 6% of cases, FLOP100 100 in 10%, the remaining algorithms in none: for the Alarm network, fLOP, does in 2% of cases, FLOP 20 in 74%, FLOP 100 in 82%, BOSS in 6%, GES in 16%, DAGMA and PC in none.
a = 0.01 for PC and ADAGMA = 0.02. As metric of accuracy, we report the Structural Hamming Distance (SHD) for CPDAGs, that is, the number of node pairs with differing edge relations in the compared graphs. If a method, such as DAGMA, returns a DAG, we first compute the corresponding CPDAG and compare this to the CPDAG of the true DAG (as generally consider assumptions where only the CPDAG is identifiable). For some settings, we also report the Ancestor Adjustment Identification Distance (AID), measuring the mistakes when using the learned instead of the true CPDAG for the downstream task of causal effect identification (Henckel et al., 2024). For the PC algorithm, which does not always return graph satisfying the invariants of CPDAGs, such as acyclicity, we report the AID only on runs that produced valid CPDAG.
Figure 1 shows run-time versus SHD (lower left is better). On SF graphs, the order-based algo- rithms clearly outperform PC, GES, and DAGMA; FLOP with ILS improves further. Even with 100 ILS restarts, FLOP's run time is comparable BOSS. On the Alarm network instances, the improvements through ILS are even more apparent, and with FLOP obtains near-perfect results.
Graphs returned by FLOOP achieve lower SHD than competing score-based methods due to better optimization of the BIC score. This is shown in Figure 5, where we report the BIC score difference to the ground-truth DAG. Generally, the results look qualitatively similar to the SHD plots for the presented settings. However, for the SF graphs, it can be seen that the BIC, c.g., for FLOP is close to zero, whereas the SHD for many instances lies clearly above zero. In fact, for majority of runs, the BIC score of the graph found by FLOP is even (slightly) better than the BIC score of the ground-trume graph showing that the global BIC optimum does not identify the ground truth in these cases.
We also evaluated the DAGMA loss function with MLE parameters fitted to the graph returned by FLOP and observed this to produce lower loss compared to the graph and parameters returned by DAGMA itself. This casts doubt on the idea that gradient-based methods relying on diffcrentiable DAG-constraints have an inherent advantage in optimizing their target score compared to discrete search. While these methods may offer other benefits, our results suggest that those likely come from aspects other than optimization quality.
DISCUSSION
In score-based causal discovery, two questions arise: (1) Is the true graph score-optimal? and (2) Can find score-optimal graph? Here introduce FLOOP, efficient discrete optimization

---

Published as a conference paper at ICLR 2026
ER, 50 nodes, average degree
SF, 50 nodes, density parameter 4
BIC true
BlCalgo
BIC irue
BlCalgo
5000
3000
4000
GES
GES
2000
3000
2000
1000
FLOPo
1000
FLOP&
FLOP20
FL.P 100 BOSS
FLOP
20
BOSS
0
0
FLOP 00
0
2
2
5
Runtime in seconds
Runtime in seconds
Figure 5: Run-time plotted against the BIC difference to the ground-truth graph for ER graphs on the left and SF graphs on the right. For the ER graphs, FLOP, finds a graph with better or equal BIC score than the true graph in 48% of cases, FLOP20 and FLOP 100 in 84% of cases, BOSS in 52% of cases and GES in 0% cases. For the SF graphs, FLOP, finds such graph in 6% of cases, FLOP in 76% of cases, FLOP100 in 94% of cases, BOSS in 6% of cases and GES again in 0% of cases.
algorithm to search for graphs minimizing the Gaussian BIC score. Indeed, when the true graph of a linear ANM has the globally optimal Gaussian BIC, FLOP typically recovers this graph or outputs one very close to it. This is further supported by simulations in Appendix C including uniform noise, unstandardized data, data based on an adaptation of the Onion method (Andrews & Kummerfeld, 2024), and real-world networks from bnlearn (Scutari, 2010). Across these settings, FLOP attains sere-oofthh-arr accuracy. typically achieving better BIC and lower SHD in a fraction of the run- oumtimme of competing methods. When assumptions are violated (Appendix C.7-C.9) or sample sizes are too small for asymptotic guarantees to hold (Appendix C.4 and Figure 11 in Appendix C.6), FLOP still optimizes the Gaussian BIC as intended and finds graphs with better BIC score than the ground truth, but graph recovery suffers because the scoring criterion does not identify the true graph.
These results highlight that it reasonable to revisit and embrace discrete search for causal struc- ture learning. FLOP often finds graphs that are score-optimal or score better than the target graph and makes the link between accuracy and speed explicit: A computationally efficient search permits more exploration through iterated local search and that leads to better-scoring graphs. Our findings also recalibrate what is considered hard. First, ER graphs with 50 nodes and about 200 edges are often presented as challenging, yet for linear ANMs order-based discrete search solves them reliably and quickly. Second, on widely used linear benchmarks, the discrete optimization of the Gaussian BIC feasible, not bottleneck, and often more efficient and reliable than other optimization meth- ods. Instead, key challenges in causal discovery lie in designing and selecting appropriate scoring criteria that identify the true graph as score-optimal not only asymptotically but with high probability also on finite samples.
At the same time, advancing causal discovery in practice remains difficult even on small graphs, since the ground truth is rarely known and assumptions are violated. It has been feasible for decades to find a global BIC optimum with exact exponential-time search up to roughly 30 vari- ables (Koivisto & Sood, 2004; Silander & Myllymäki, 2006). FLOP extends strong BIC optimiza- tion to substantially larger graphs, but that does not make the practical problems go away. Our work shifts the attention away from inflated combinatorial hardness rhetoric and from misattributed gap between asymptotic theory and observed finite-sample performance, toward the immense challenges causal discovery faces outside of synthetic benchmarks (Reisach et al., 2021; Göbler et al., 2024; Mogensen et al., 2024; Brouillard et al., 2025; Gamella et al., 2025; Gururaghavendran & Murray, 2025; Jørgensen al., 2025).
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
Marcel Wienöbst thanks Kenneth Langedal for fruitful discussions and introducing him iterated local search. Sebastian Weichwald was supported by a research grant (0069071) 1) from Novo Nordisk Fonden. This research was supported by the Pioneer Centre for AI, DNRF grant number PI.
REFERENCES
Steen A. Andersson, David Madigan, and Michael D. Perlman. A characterization of Markov equiv- alence classes for acyclic digraphs. The Annals of Statistics, 25(2):505-544 1997.
Bryan Andrews and Erich Kummerfeld. Better simulations for validating causal discovery with the DAG-adaptation of the onion method. arXiv preprint arXiv:2405.33/00 2024.
Bryan Andrews, Joseph Ramsey, Ruben Sanchez Romero, Jazmin Camchong, and Erich Kummer- feld. Fast scalable and accurate discovery of DAGs using the best order score search and grow shrink trees. Advances in Neural Information Processing Systems, pp. 63945-63956, 2023.
Ingo A. Beinlich, Henri Jacques Suermondt, R. Martin Chavez, and Gregory F. Cooper. The ALARM monitoring system: A case study with two probabilistic inference techniques for be- lief networks. In European Conference on Artificial Intelligence in Medicine, pp. 247-256, 1989.
Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. DAGMA: Learning DAGs via M-matrices and log-determinant acyclicity characterization. Advances in Neural Information Processing Systems, pp. 8226-8239, 2022.
Philippe Brouillard, Chandler Squires, Jonas Wahl, Karen Sachs, Alexandre Drouin, and Dhanya Sridhar. The landscape of causal discovery data: Grounding causal discovery in real-world appli- cations. In Causal Learning and Reasoning, pp. 834-873, 2025.
David Maxwell Chickering. Learning Bayesian networks is NP-complete. In Learning from Data: Artificial Intelligence and Statistics V. pp. 121-130. Springer, 1996.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3:507-554, 2002.
David Maxwell Chickering and Christopher Meek. Selective greedy equivalence search: Finding optimal Bayesian networks using a polynomial number score evaluations. In Uncertainty in Artificial Intelligence, pp. 211-219, 2015.
David Maxwell Chickering, David Heckerman, and Chris Meek. Large-sample learning of Bayesian networks is NP-hard. Journal of Machine Learning Research, 5: 3:128711300, 2004.
Tom Claassen, Joris M. Mooij, and Tom Heskes. Learning sparse causal models is not NP-hard. In Uncertainty in Artificial Intelligence, pp. 172-181, 2013.
Rina Foygel and Mathias Drton. Extended Bayesian information criteria for Gaussian graphical models. In Advances in Neural Information Processing Systems, pp. 604-612, 2010.
Juan L Gamella, Jonas Peters, and Peter Bühlmann. Causal chambers as a real-world physical testbed for Al methodology. Nature Machine Intelligence, 7(1):107-118 2025.
Philip Gill, Gene H. Golub, Walter Murray, and Michael A. Saunders. Methods for modifying matrix factorizations. Mathematics of Computation, 28( 56):555-333 1974.
Konstantin Göbler, Tobias Windisch, Mathias Drton, Tim Pychynski, Martin Roth, and Steffen Son- ntag. AassalAssembly: Generating realistic production data for benchmarking causal discovery. In Causal Learning and Reasoning, pp. 609-642, 2024.
Gene H. Golub and Charles F. Van Loan. Matrix Computations. JHU press, 2013.
Rajesh Gururaghavendran and Eleanor J. Murray. Can algorithms replace expert knowledge for causal inference? A case study on novice use of causal discovery, American Journal of Epidemi- ology, 13995): 2025.
11

---

Published as a conference paper at ICLR 2026
David Heckerman, Dan Geiger, and David M. Chickering. Learning Bayesian networks: The com- bination of knowledge and statistical data. Machine Learning, 2073):197-243, 1995.
David Earl Heckerman, Eric Horvitz, and Bharat N. Nathwani. Toward normative expert systems: Part i the pathfinder project. Methods of Information Medicine, 31(02)900-105, 1992.
Leonard Henckel, Theo Würtzen, and Sebastian Weichwald. Adjustment identification distance: A gadjid for causal structure learning. In Uncertainty in Artificial Intelligence, pp. 1569-1598, 2024.
Allan Leck Jensen and Finn Verner Jensen. MIDAS: An influence diagram for management of mildew in winter wheat. In Uncertainty in Artificial Intelligence, pp. 349-356, 1996.
Frederik Hytting Jørgensen, Luigi Gresele, and Sebastian Weichwald. What is causal about causal models and representations? arXiv preprint arXiv:2501.19333 2025.
Mikko Koivisto and Kismat Sood. Exact Bayesian structure discovery in Bayesian networks. Jour- nal of Machine Learning Research, 5:549-577 2004.
Daphne Koller and Nir Friedman. Probabilistic Graphical Models: Principles and Techniques. MIT Press, 2009,
Jack Kuipers, Polina Suter, and Giusi Moffa. Efficient sampling and structure learning of Bayesian networks. Journal of Computational and Graphical Statistics, 31(():639-650 2022.
Wai- Yin Lam, Bryan Andrews, and Joseph Ramsey. Greedy relaxations of the sparsest permutation algorithm. In Uncertuinly in Artificial Intelligence, pp. 1052-1062, 2022.
Mingjia Li, Hong Qian, Tian-Zuo Wang, Shujun Li, Min Zhang, and Aimin Zhou. Strong and weak identifiability of optimization-based causal discovery in non-linear additive noise models. 2025.
In International Conference on Machine Learning, pp.
Xiaohan Liu, Xiaoguang Gao, Zidong Wang, Xinxin Ru, and Qingfu Zhang. A metaheuristic causal discovery method in directed acyclic graphs space. Knowledge- -aased Systems, 276:1 10749, 2023.
Helena Ramalhinho Lourenço, Olivier C. Martin, and Thomas Stützle. Iterated local search: Frame- work and applications. In Handbook of Metaheuristics, pp. 129-168. Springer, 2018
Dimitris Margaritis. Learning Bayesian network model structure from data. PhD thesis, School of Computer Science, Carnegie Mellon University, 2003.
Søren Wengel Mogensen, Karin Rathsman, and Per Nilsson. Causal discovery in a complex in- dustrial system: time series benchmark. In Causal Learming and Reasoning, pp. 1218-1236, 2024.
Preetam Nandy, Alain Hauser, and Marloes H. Maathuis. High-dimensional consistency in score- based and hybrid structure learning. The Annals of Statistics, 46(6A):3151-3188 2018.
Achille Nazaret and David Blei. Extremely greedy equivalence search. In Uncertainty in Artificial Intelligence, pp. 2716-2745, 2745, 2024.
Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and DAG constraints for learning linear DAGs. Advances in Neural Information Processing Systems, 33:17943-17954 2020.
Ignavier Ng. Biwei Huang, and Kun Zhang. Structure learning with continuous optimization: A sober look and beyond. In Causal Learning and Reasoning, pp. 71-105, 2024.
Jens Dalgaard Nielsen, Tomás Kocka, and José M. Peña. On local optima in learning Bayesian networks. In Uncertainty in Artificial Intelligence, pp. 435-442, 2003.
Gunwoong Park. Identifiability of additive noise models using conditional variances. Journal of Machine Learning Research, 21(75):1-34, 2020.
Judea Pearl. Causality. Cambridge University Press, second edition, 2009.
12

---

Published as a conference paper at ICLR 2026
Jonas Peters and Peter Bühhmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101(1):219-25) 228, 2014.
David Pisinger and Stefan Ropke. Large neighborhood search. In Handbook of Metaheuristics, pp. 99-127. Springer, 2018.
Joseph D. Ramsey, Kun Zhang, Madelyn Glymour, Ruben Sanchez Romero, Biwei Huang, Imme Ebert-Uphoff, Savini Samarasinghe, Elizabeth A. Barnes, and Clark Glymour. TETRAD A toolbox for causal discovery. International Workshop on Climate Informaties, pp. 89 -2- 2018.
Garvesh Raskutti and Caroline Uhler. Learning directed acyclic graph models based on sparsest permutations. Stat, 7(1):e183, 2018.
Alexander Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the simulated DAG! Causal discovery benchmarks may be easy to game. Advances in Neural Information Processing Systems, pp. 27772-27784, 2021.
Alexander Reisach, Myriam Tami, Christof Seiler, Antoine Chambaz, and Sebastian Weichwald. A scale-invariant sorting criterion to find causal order in additive noise models. In Advances in Neural Information Processing Systems, pp. 785-807, 2023.
Felix L. Rios, Giusi Moffa, and Jack Kuipers. Benchpress: versatile platform for structure learning in causal and probabilisticc graphical models. Journal of Statistical Software, 114:1 2025.
Paul Rolland, Volkan Cevher, Matthäus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score matching enables causal discovery of nonlinear addi- tive noise models. In International Conference on Machine Learning, pp. 18741-18753, 2022.
Karen Sachs, Omar Perez, Dana Pe'er, Douglas A. Lauffenburger, and Garry P. Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721): 523-529, 2005.
Mauro Scanagatta, Cassio P. de Campos, Giorgio Corani, and Marco Zaffalon. Learning Bayesian networks with thousands of variables. Advances in Neural Information Processing Systems, pp- 1864-1872, 2015.
Gideon Schwarz. Estimating the dimension of a model. The Annals of Statistics, pp. 461-464, 1978.
Marco Scutari. Learning Bayesian networks with the bnlear R package. Journal of Statistical Software, 35:1-22, 2010.
Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvarinen, Yoshinobu Kawahara, Takashi Washio, Patrik O. Hoyer, and Kenneth Bollen. DirectLiNGAM: A direct method for learning a linear non-Gaussian structural equation model. Journal of Machine Learning Research, 12:1225-1248, 2011.
Tomi Silander and Petri Myllymäki. A simple approach for finding the globally optimal Bayesian network structure. In Uncertainty in Artificial Intelligence, pp. 145-452, 2006.
Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search. MIT Press, second edition, 2000.
Mare Teyssier and Daphne Koller. Ordering-basse search: a simple and effective algorithm for learning Bayesian networks. In Uncertainty in Artificial Intelligence, pp. 584-590, 2005.
Thomas Verma and Judea Pearl. Causal networks: Semantics and expressiveness. In Uncertainty in Artificial Intelligence, pp. 69-78, 1988.
Thomas Verma and Judea Pearl. Equivalence and synthesis of causal models. In Uncertainty in Artificial Intelligence, pp. 255-270, 1990.
Xun Zheng, Bryon Aragam, Pradeep K. Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Con- tinuous optimization for structure learning. Advances in Neural Information Processing Systems, pp. 9492-9503, 2018.
13

---

Published as a conference paper at ICLR 2026
Perturb Whett with random swaps
Thest
Repeat l times (ILS) Section 4.2
π
D
Find initial order cf. Section 4.1
π
Local search
n
Update bestt to x' if BIC(x') < BIC (best)
G (Test))
cf. Algorithms 1, 3,
Figure 6: Visualization of the general control flow in the FLOP algorithm.
A HIGH-LEVEL DESCRIPTION OF THE FLOP ALGORITHM
As an overview of the FLOP algorithm, we provide Figure 6, describing the high-level control flow of FLOP. The algorithm begins by computing an initial order for the given data set D. Afterwards, the ILS loop starts with a local search aiming to improve the order π through reinsertions. Reinser- tions done by moving a node possible (Algorithm 1). This optimal reinsertion is computed as described in Algorithm 2, which relies on the grow-shrink described Algorithm for updating the parents. This grow-shrink starts at the previous parent set and uses efficient Cholesky updates for scoring as described in Section 3.2.
its locally optimal position until no further improvements are
After the local search completes, the best found graph/order updated (Trest in Figure 6) if the score is lower than the previous best. This best- best-scoring order found thus far is then perturbed as described in Section 4.2 and the perturbed copy then used as the starting point for the next local search. We note that this procedure ensures that more ILS iterations can only improve the score of mest and thus of the returned graph G (Thea) because Thest is only updated if local search after the perturbation yields better scoring order.
We note that the reinsertion-basee local search follows the general principle of the BOSS algorithm, which consists of reinserting nodes their locally optimal position until no further (local) improve- ments are possible. However, FLOP and BOSS differ in the following key aspects:
The grow-shrink procedure of BOSS starts with the empty set. Moreover, it is a greedy grow-shrink that always inserts or removes the node with the largest local score improve- ment. In contrast, FLOP accepts any improving insertion in the grow-, and any improving removal in the shrink-phase For the implementation of grow-shrink BOSS relies at its core on an intricate data structure called grow-shrink trees, which FLOOP avoids. Overall, this allows FLOP to obtain a better run-time performance compared to BOSS.
FLOP uses Cholesky updates for efficient iterative scoring during the local search and, in particular, in the ghrwnk routine. This yields further run-time gains.
FLOOP makes use of iterated local search (ILS) that allows spending improved BIC optimization. As more ILS restarts can never yield worse scoring graphs, this effectively trades off compute with accuracy. Due FLOPs run-time improvements this yields free lunch with regard to accuracy gains.
compute for
BOSS starts the local search with a random order. This leads performance deteriorations on path instances as shown in Section 4. In contrast, FLOP explicitly constructs the initial order to avoid such problems.
B CHOLESKY DECOMPOSITION OF THE COVARIANCE MATRIX
Let X,...., Xp be real-valued centered random variables with finite second moments and with full-rank covariance matrix Σ
E RPxP, that is, for all j E (I,...,P) E(X;) and Σş Var(X,) < oo; further, Σ is symmetric positive definite and admits unique Cholesky factorization E = L.T with I. lower triangular and strictly positive diagonal.
14

---

Published as a conference paper at ICLR 2026
For E {1,.... p), let X, denote the best linear predictor of X, from its predecessors Xi, X5-1, that is, the ordinary least squares projection.
Then for all j E pl....,p.
12 Var(X, x))
that is, Lis is the standard deviation of the least squares residuals when linearly regressing X, onto its predecessors.
To obtain the statement, ji E (1,...,p) p). Block-partition the leading j principal submatrices of Σ and L:
Eigi- 1)---)
11:(i-1).5
L1. 1). 1:05-1)
1:1-11
Lig
From == LL we get the block identities
E' L'L',
L'T c- +l.
Since L' is full rank, the second identity gives T = 1'-18, hence
rr = (L' s=
s (L'L') == 1'-1s.
Substituting into rr 12 yields
= E'-ls.
On the other hand, the centered OLS problem
min aER
E [Xx a X1.0-1))
has normal equations E'a" 8, so a* E'-1s i the unique minimizer, X, -T--1 STE''1X1.G--11 and the minimal mean squared error is
Var(X; 82)= = Var (X,) + Var(X,) - Cov(
+ s E'-1 Var(X1.6 E'
2s TT' Cov(X X x1:(j-1))
= + 8 EE-EE'' - 28 Ty-
= c TT-1-
This equals e, that is, 12 Var(X, Xz), claimed.
We remark that if X1...., Xp are jointly Gaussian, then Var(X, - xX,) = Var(X, ( (1)
C FURTHER BENCHMARK SETTINGS
In this section, we consider further benchmark settings to investigate the stability of FLOP under dif- ferent graph and data generation procedures. If not specified otherwise, consider ER graphs with 50 nodes and average degree 8, and 1000 samples being drawn from the underlying linear additive noise model. We also consider settings where the assumptions of FLOP with linear Gaussian BIC are (potentially) violated, such as uniform noise (Subsection C.I) and non- non-linear relations (Subsec- tion C.7) as well as semi-synthetic (Subsection C.8) and real-world data (Subsection C.9). In some settings, such as the uniform-noise case this leads to no performance degradation, whereas in others leads to significantly higher SHDs compared to settings where the assumptions are satisfied. We observe that in these settings, FLOP still reliably optimizes the BIC, meaning the performance of FLOP large parts depends on how well the scoring criterion suited to the data. Hence, practi- tioners need to be careful that the assumptions hold when applying FLOP. Moreover, this shifts the focus research from designing optimization algorithms towards the development of efficient and practical scoring criteria (see also the discussion in Section 6).
15

---

Published as a conference paper at ICLR 2026
Uniform noise
Unstandardized
SHD
SHD
400
400
300
GES
300
GES
200
200
PC
DAGMA
PC
100
100
DAGMA
FLOPo
FLOP
FLOPAR
BOSS
LiNGAM
Fl.020
20
BOSS
0
FLOP 100
FLOP
0
5
10
15
0
5
10
Runtime in seconds
Runtime in seconds
Figure 7: Run-time against SHD for data sampled with uniform instead of Gaussian noise the left and for unstandardized data on the right (both settings are based on ER graphs with 50 nodes and average degree with 1000 samples drawn). In the uniform noise case, FLOP, and BOSS find the target CPDAG in 34% of cases, FLOP20 and FLOP in 54% of the cases, the remaining algorithms none. On unstandardized data, BOSS finds the target CPDAG 22% of cases, FLOP, in 34% and FL.O2o and 10.OP100 in 54% of cases, the remaining algorithms in none.
C.I UNIFORM NOISE
To check the performance of FLOP (using the Gaussian BIC to learn the CPDAG underlying a linear ANM) under non- Gaussian noise, we generate data with noise sampled uniformly from 1]. As the plot on the left Figure shows, there no performance degradation (of any algorithm). Moreover, compared the methods to DirectLiNGAM (Shimizu et al., 2011), which is based on identifiability theory for non-Gaussian noise. DirectLinGAN, gets low SHD on these instances, but in the 50 repetitions
recovered the ground truth.
C.2 Raw DATA
To avoid varsortability of the instances (Reisach et al., 2021), we typically standardize the data in the benchmarks as mention Section As an exception, we consider instances with unstandardized data the right of Figure 7. As expected, we find that DAGMA performs significantly better than the standardized settings. The performance of the other algorithms does not vary significantly. We also note that in FLOP choose to always standardize the data to obtain algorithm.
C.3 DAO DATA
We also consider the DAG-adaptation of the Onion method (Andrews & Kummerfeld, 2024) way to sample data from an ANM. This method has been proposed to avoid artefacts in the data, such as 122soottbyylity (Reisach al., 2023), which could be inadvertently or explicitly exploited game benchmarks. In line with the results by Andrews & Kummerfeld (2024), we find that this sampling method yields harder-to-identify instances, with FLOP nonetheless performing best; however, all methods produce SHDs greater than 50, as shown in the left panel of Figure 8. This may be caused by weak causal relationships or faithfulllu violations the data. is, however, not failure in the optimization, of the learned graph was better than the one of the ground truth, suggesting non identifiability the true CPDAG under the BIC for the provided number of samples.
observed that, for FLOP and other score-based algorithms, the BIC score
16

---

Published as a conference paper at ICLR 2026
Onion method
Dense ER graphs, 1000 samples
SHD
SHD
FLOPo
200
BOSS
200
BOSS
DAGMA
FLOH
BOSS20
150
FLOP 100
150
FLOPb
GES
PC
FLOP
500
Exact
100
BOSS
100
50
20.OPo FI.OP100 100
50
0
1
2
3
4
5
0
10
20
30
40
50
Runtime in seconds
Runtime seconds
Figure 8:
against SHD for data sampled with the DAG-adaptation the Onion method on the left (again on ER graphs with 50 nodes and average degree 8 with 1000 samples drawn) and for dense graphs (25 nodes, average degree 16) with 1000 samples generated in the standard way on the right. In both settings, none of the algorithms ever recover the target CPDAG.
ER, 250 nodes, average degree 8
ER, 500 nodes, average degree 8
SHD
SHD
1000
PC
PC
1000
500
GES
500
GES
FLOP
FLOP 100
BOSS
FLOP
FLOPa0
FLOP100
FL.OPp 20
0
100
200
300
500
1000
Runtime in seconds
Runtime in seconds
Figure 9: Run-time against SHD for ER graphs with 250 nodes on the left and with 500 nodes on the right (average degree 8 and 1000 samples drawn). BOSS times out the latter instances. In both settings, none of the algorithms ever recover the target CPDAG.
17

---

Published as a conference paper at ICLR 2026
Mildew network
Barley network
SHD
SHD
60
100
DAGMA
DAGMA
GES
40
F FLOPo
FLOP
30
50
GES BOSS
20
PC BOSS
10
20
PC
FLOPP
20
0
0
FLOP
0
Runtime in seconds
Runtime in seconds
Figure 10: Run-time against SHD for the Mildew network on the left (Jensen & Jensen, 1996), which consists of 35 nodes and 46 edges, and the Barley network on the right (Scutari, 2010), which consists of 48 nodes and 84 edges. For the Mildew network, BOSS finds the target CPDAG in 2% of cases, FLOP20 20 in 48% and FLOP in 52% of the cases, the remaining algorithms in none. For the Barley network, GES finds the ground truth in 4% of cases, FLOP, in 8%, BOSS finds the target CPDAG in 12%, FLOP 20 in 90% and FLOP 100 94% of cases, PC and DAGMA in none.
C.4 DENSE ER GRAPHS
In the main paper, we considered dense ER graphs (25 nodes and average degree 16) in a setting with 50,000 samples. Due to the denseness of the graph such large amount of samples is necessary to identify the target graph. Here, we show the performance of the algorithms for significantly fewer samples, namely 1000 samples, as in the other simulations. As can be seen on the right of Figure 8, FLOP still performs quite well, however, the algorithms are much closer with regard to the SHD.
We note that we again compared the BIC score of the graph returned by FLOP with the ground-truth graph as well as the other algorithms. We found that FLOP and the exact algorithm consistently found graphs with a better BIC score than the ground truth or the other discovery algorithms' output graphs. This illustrates that the sample size is too small for the BIC score to reliably identify the true CPDAG.
Another thing to note is that compared to the setting with 50000 samples in the main text, both BOSS and FLOP run faster on instances with 1000 samples, whereas there no noticeable difference for the exact algorithm. The reason for this increased run-time for larger samples sizes is that the BIC penalizes edges stronger for smaller sample sizes with the penalty term growing with In n and the likelihood term proportional with Thus, intermediate graphs in the search are typically denser for high-sample settings, which increases the computational effort.
C.5 LARGE ER GRAPHS
We also report the accuracy for large ER graphs with 250 and 500 nodes and average degree 8 in Figure 9. Here, DAGMA does not terminate within the time limit for either instances and BOSS does not for the graphs with 500 nodes. Overall, similar accuracy results as before can be observed though notably PC appears get worse with an increased number of variables in comparison with GES.
C.6 BNLEARN GRAPHS
In addition to the random graphs, also consider real-world networks from the bulearn repository, namely the Mildew (Jensen & Jensen, 1996), Barley and the Pathfinder network (Heckerman et al., 1992). All three large such that exact score-based algorithm based on dynamic programming
18

---

Published as a conference paper at ICLR 2026
Pathfinder network
Runtime in seconds
SHD
FLOPo
FLOP
100
FLOP20
BOSS
FLOP
BOSS
-100
50
-150
FL.PPo
FLOP 100
0
8
BIC
BIC
Runtime in seconds
Pathfinder network
Figure 11: Run-time against SHD for the Pathfinder network, which consists of 109 nodes and 195 edges, on the left. On the right, the BIC score difference compared to the ground-truth graph for the Pathfinder network. None of the algorithms ever recover the target CPDAG.
could be used, with Mildew consisting of 35 nodes and 46 edges, Barley of 48 nodes and 84 edges, and Pathfinder of 109 nodes and 195 arcs. same manner before. For Mildew and Barley on the left and right of Figure 10, FLOP performs significantly better than other methods and, in particular, that the ILS is needed to get close-to- perfect accuracy on these instances. For Pathfinder do not terminate within the time limit of 30 minutes. Here, FLOP, and BOSS yield roughly similar SHD. However, with an increasing number of ILS iterations, the SHD gets worse for FLOP and FLOP100 To analyze this behaviour further, we show the BIC score difference to the ground-truth DAG on the right of Figure 11. Indeed, all reported methods yield better BIC scores than the true DAG and ILS does find even better-scoring graphs, which, in this case, are further from the ground truth. Again, faithfulness violations promoted by the underlying graph structure may be the issue here, even though closer investigations are needed.
cases, we generate the data synthetically in the
the left of Figure 11, PC, GES and DAGMA
C.7 NON-LINEAR DATA
As settings where the linear Gaussian BIC is misspecified, we consider non-linear data generated from randomly initialized multi-layer perceptron (MLP) with single hidden layer of size 100 and sigmoid activation, as described in Appendix C.2.2 in (Bello et al., 2022) and from sampled Gaussian process regressions with unit bandwidth RBF kernel as proposed in (Rolland et al., 2022). In both settings, the ground-truth DAG is generated by orienting an ER graph with 25 nodes and average degree 4, thus containing on average 50 edges, according to linear order that drawn uniformly at random. We consider the same algorithms as before with the same parameter choice and score. They are hence not tuned towards the non-linear setting. Additionally, we include the non-linear DAGMA algorithm from Bello et al. (2022). In Figure 12, we plot the SHD of each method contrasted with the BIC difference to the optimal BIC score (for ABIC 2) for each of the algorithms (cases where PC does not return valid CPDAG are omitted). As can be seen, the BIC optimum does not correspond low SHD in both settings with GES, BOSS and FLOP 100 having similar performance, and the ground truth having suboptimal BIC scores. For the MLP setting, the non-linear version of DAGMA the best method for graph recovery, however, in the GP setting is not better than the other approaches. It is also the by far slowest method, taking over minutes per instance. The assumptions of the LiNGAM algorithm are also violated by the non-lineenn lincariiiies and clearly the worft-performing algorithm among the presented ones.
19

---

Published as a conference paper at ICLR 2026
Non-linear (MLP)
Non-linear (GP)
SHD
SHD
60
LINGAM
100
50
GES PC
DAGMA
LiNGAM
40
BOSS
DAGMA Ance lincur
FLOP 100
30
50
GES
PC
20
COSSS
DAGMA
FLOP
10
DAGMAppc--iiee
0
true
true m.0
1000
2000
3000
4000 5000 6000
100
200
300
500
BICopt BICalgo
BIC BICalgo
Figure 12: BIC score difference the BIC optimum plotted against SHD for non-linear data gener- ated based on MLPs with a single hidden layer on the left. On the right, non-linearities are generated from sampled Gaussian process regressions with a unit bandwidth RBF kernel. Both settings use ER graphs with 25 nodes and average degree 4, thus the ground truth contains 50 edges on average. In the MLP setting, FLOP 100 finds the BIC optimum in 44% of cases, BOSS finds it in 12% of cases and GES in 4% of cases. In the GP setting, FLOP100 finds the BIC optimum in 72% of cases, BOSS in 16% of cases and GES in 30% of cases.
C.8 CAUSALASSEMBLY DATA SET
We show the results on the causalAssembly dataset introduced by Göbler et al. (2024) in Figure 13. The goound-truth DAG consists of 98 nodes and 485 edges. We subsample 5000 observations with replacement from the data set 50 times and run the algorithms on this subsampled data. We exclude DAGMA and LiNGAM from the plots as they yield significantly larger SHD, which lies above 550, and take much longer than the competing algorithms, namely more than 30 seconds in the case of DAGMA and more than 200 seconds the case of LiNGAM. The remaining algorithms return results of similar quality, with notable improvements through the ILS restarts that FLOP uses. These small improvements stem from better BIC score optimization as shown in the right plot. Here, the BIC difference to the true graph is reported and it is clear that all methods return graphs with much better BIC scores than that of the ground truth, suggesting score misspecification.
C.9 SACHS DATA SET
We also evaluate the algorithms on the Sachs dataset (Sachs et al., 2005), which consists of 11 nodes, and compute the SHD with regard to the CPDAG of the ground truth consisting of 17 edges. We each algorithm using the same hyperparameters observations in the data set. As result, we observed FL.OPnn BOSS and GES performing on par, all yielding an average SHD of 12.58. The other algorithms yield similar results, with DAGMA having the best performance with an average SHD of 11.7. The PC algorithm obtains an average SHD of 12.56 and LiNGAM an average SHD of 14.18.
before 50 bootstrap samples of the 853
C.10 RESULTS FOR DIFFERENT PARAMETERS CHOICES
FLOP has two parameters that need to be chosen by the user. First, ABIC scales the penalty term of the BIC and, second, the number of ILS restarts control the amount of compute that is invested. For the latter parameter,
have typically shown the simulation results for multiple choices, such as FLOP,, when no restarts are performed, as well as FLOP2 and FLOP with 20 and a 100 restarts, respectively. We also note that more ILS iterations
only improve the BIC optimization.
For ABIC on the other hand, we choose the value 2, which is the standard setting from the literature. As Foygel & Drton (2010) have shown, larger value than should be chosen recover the structure
20

---

Published as a conference paper at ICLR 2026
causalAssemblyy 5000 samples
Runtime in seconds
SHD
3
5
400
GES
0
FI.OPo FL.OP20
BOSS
FLOP 100
-1000
300
-2000
200
-3000
100
GES
4000
FLOPp FLOP20
FLOP
0
BOSS
BICme BIC algo
Runtime in seconds
causalAssembly y,, 5000 samples
Run-time against SHD for the causalAssembly data on the left. On the right, the BIC score difference to the BIC of the ground truth reported. LiNGAM and DAGMA not in- cluded because both are considerably slower on these instances (with DAGMA needing more than 30 seconds per instance and LiNGAM more than 200 seconds) and obtain significantly worse SHD compared to the other methods, typically above 550 for both algorithms. As point-of-reference, the ground-truth graphs consists of 485 edges, thus these two methods give worse SHDs than the empty graph. All methods optimizing the BIC, shown in the right graph, yield BIC scores clearly lower than that of the true graph, indicating score misspecification for the linear Gaussian BIC.
Figure 13:
FLOP on ER, 50 nodes, average degree 8
FLOP on ER, 50 nodes, average degree 16
SHD
SHD
80
400
70
60
300
50
40
200
30
20
100
10
0
3
1
ABIC
In random swaps
Figure 14: SHD for FLOP with different choices of ABIC on the left (ER graphs with 50 nodes and average degree 8) and with different factors controlling the number of swaps a perturbation on the right (ER graphs with 50 nodes and average degree 16). For the ABIC. as explained by the derivation of the extended BIC (Foygel & Drton, 2010), values higher than samples with 2 being common choice. For the perturbations, it can be seen that many choices for number of random swaps are effective (the exception being no perturbations, and thus no ILS at all, which is shown at x 0, yielding an SHD in the hundreds for this setting), with outliers increasing for too few or too many swaps.
needed on finite
21

---

Published as a conference paper at ICLR 2026
of graphical models, while any constant value guarantees asymptotic consistency. The results on the left of Figure 14 confirm this, showing that for λBic larger or equal to 2, graphs close to the ground truth are recovered by FLOP while smaller choices of ABIC yields spurious edges and thus higher SHD.
Finally, for the ILS perturbations, FLOP defaults In many random swaps. The software interface of FLOP does not allow tuning this parameter, as we found it to be stable default choice. This is confirmed
the right of Figure 14, which runs FLOP with Inp many random swaps for & E {0, 1/, 1/2, 3/4, 1, 4/3, 2, 4). In the case that x is set to zero, which corresponds to not running ILS, this yields an SHD that often in the hundreds. Conversely, any of the positive choices of r lead to good performance. The best results are obtained for z between 3/4 and 2, while for the largest and smallest values of the number of outliers increases.
C.I1 ANCESTOR ADJUSTMENT DISTANCE
We show the Ancestor Adjustment Identification Distance (AID) as another metric for evaluating the learned graphs (Henckel et al., 2024). It effectively counts the number of mistakes one would make if one used the learned graph to select valid adjustment sets (using the ancestors of node) instead of the ground-truth graph. Figure 15 shows the AIDs for a selection of the previous simulation results. We note that we omit data points where PC does return CPDAG (as well-known to happen on finite samples). For example, on ER graphs with 500 nodes, the PC algorithm does not yield a single valid CPDAG.
C.12 OTHER GRADIENT-BR BASED ALGORITHMS
Due to the choice of the least-squares loss function, for DAGMA and other popular gradient-based methods, one implicitly assumes perfect varsortability (or even equal noise variances in the un- derlying linear ANM) to recover the underlying DAG as the unique score-optimal DAG (Peters & Bühhlann, 2014; Reisach et al., 2021; 2023). Because we compare the methods on standard- ized data in simulations (with the execption of Subsection C.2), this assumption is violated. In comparison, the GOLEM algorithm proposed by Ng et al. (2020) has a version that is developed specifically for the non-equal noise variance case. We ran this algorithm available in the gCastle package (version 1.0.4) with default parameters on the standard setting in this work (ER graphs with 50 nodes and average degree 8) and observed an average SHD of 176.6. This slightly better than DAGMA with an SHD of 195.08 (which may also be partially due to different hyperparameter choices), but clearly worse than FLOP 100 with an average SHD of 1.34. The GOLEM algorithm also had run-time of more than 200 seconds per instances, whereas running FLOP 100 took around seconds per instances.
22

---

Published as a conference paper at ICLR 2026
SF, nodes, density parameter
DAG-adaptation of the Onion
AID
AID
GES
GES
PC
1000
1000
DAGMA
FLOP20 220
DAGMA
PC
BOSS
FLOPo
F FLOP
FLOP
500
500
BOSS
FLOP20
FLOP100
0
0
1
2
3
4
5
7
0
1
3
4
5
Runtime in seconds
Runtime in seconds
ER, 250 nodes, average degree 8
ER, 500 nodes, average degree 8
AID
AID
20000
60000
PC
50000
15000
40000
GES
10000
30000
GES
20000
5000
FLOP
BOSS
10000
FLOP
FLOP 100.
FLOP 100
FLOP2u 20
FLOPzo
0
100
200
300
400
500
1000
Runtime in seconds
Runtime in seconds
Alarm network
Pathfinder network
AID
AID
9000
8000
FLOP 20
7000
FLOP 100
1000
DAGMA
6000
5000
FLOP
BOSS
FLOPo
GES
4000
500
3000
BOSS
2000
PC
FLOP 20
1000
FLOP
0
0
3
0
1 2
4
5
6
8
Runtime in seconds
Runtime in seconds
Figure 15: Run-time and AID for SF graphs (top left), data sampled by the DAG-adaptation of the Onion method (top right), ER graphs with 250 nodes (center left), ER graphs with 500 nodes (center right), the Alarm network (bottom left) and the Pathfinder network (bottom right).
23


Reviews:
[{"review_id": "WmKZh85fdo", "reviewer": "Reviewer_MbRD", "summary": "This paper introduces FLOP (Fast Learning of Order and Parents), a highly efficient score-based algorithm for causal structure learning in linear models. FLOP significantly accelerates the reinsertion-based local search over topological orders by \"warm-starting\" the grow-shrink parent selection and using $O(k^2)$ Cholesky rank-one updates for score computation, avoiding $O(k^3)$ re-computation . This speed advantage is leveraged to run an Iterated Local Search (ILS) with a principled order initialization, enabling the algorithm to escape local optima and achieve accurate graph recovery with dramatically lower runtime. The work compellingly challenges the conventional notion that discrete search cannot scale outside of small graphs.", "questions": "1. How does this approach compare to heuristic adaptions of GES such as XGES [1], where heuristics similar to ILS are used to escape local minima? Theoretical and experimental comparison would both be useful here.\n\n2. How does the proposed approach here compare empirically against recent score-based methods such as GRaSP [2], or against FCM methods designed for linear ANMs, such as DirectLiNGAM [3]? In general, it would be useful for the authors to undertake a more comprehensive empirical comparison to other recent methods, as many of the paper’s contributions relate to the promise of improved empirical performance, but without any formal guarantees of improvement.\n\n3. To what extent can the heuristic modifications made in FLOP be translated to the setting of discrete search in nonlinear ANMs? Some modifications seem like they might likely generalize such as the modified grow-shrink procedure, but other improvements such as dynamic Cholesky updates seem only possible due to the linearity assumptions.\n\n4. Is it possible to recalculate the ordering metric to use some more common notion of topological divergence, as seen unnormalized in [4] or normalized in [5,6]? This would provide additional insight into the performance of FLOP.\n\n5. Please give more details on this: “The standard hardness constructions rely on data-generating processes that cannot be represented by a DAG over the observed variables (Chickering, 1996; Chickering et al.,2004).” Not sure if I understand correctly what it means for a DGP to not be representable as a DAG over the observed variables (do you mean they require unobserved variables?).\n\n6. In Figure 2, can the authors clarify what the difference between FLOP naive and FLOP (reusing parents) is? My understanding is that the algorithm generated by modifying grow-shrink is exactly FLOP? If not, what exactly is the base FLOP algorithm? \n\nReferences:\n\n[1] Extremely Greedy Equivalence Search, Nazaret et al., UAI 2025.\n\n[2] Greedy Relaxations of the Sparsest Permutation Algorithm, Lam et al., UAI 2022.\n\n[3]DirectLiNGAM: A Direct Method for Learning a Linear Non-Gaussian Structural Equation Model, Shimizu et al., JMLR 2011.\n\n[4]Causal Discovery with Score Matching on Additive Models with Arbitrary Noise, Montagna et al., CleaR 2023.\n\n[5] Hybrid Top-Down Global Causal Discovery with Local Search for Linear and Nonlinear Additive Noise Models., Hiremath et al., NeurIPS 2024\n\n[6] LoSAM: Local Search in Additive Noise Models with Mixed Mechanisms and General Noise for Global Causal Discovery., Hiremath et al. UAI 2025", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "This paper builds upon a prior discrete search algorithm BOSS, by proposing 4 distinct modifications: 1) reusing parent sets for the grow-shrink procedure to improve runtime, 2) dynamically updating the choleksy decomposition to improve runtime, 3) initializing the first order to mitigate against known failure modes of order-based methods, and 4) running iterated local search to escape local minima. \n\nModifications 1), 2), and 3) seem entirely original to my knowledge (for 4) see the authors response to my first question in the Questions section). The quality of the work is high, as proofs provided make sense, all notation is clear, standard definitions and concepts in causal discovery are consistently used, and extensive coverage of the relevant literature is explored. The paper is generally clear, with a few caveats that I discuss in the Weaknesses and Questions section. Most importantly, the empirical results concerning the reduction in runtime are quite significant - in particular, Figure 2 shows how the addition of each modification leads to substantial decreases in the runtime needed for a fixed number of nodes, enabling scaling to graphs with over 500 variables. This is extremely significant given the fact that most discovery algorithms in the literature fail to scale past 30-50 nodes. In general, FLOP builds upon a growing trend of papers [1,2, 3] exploiting local structure previously ignored by methods aiming only for consistency. Additionally, this paper puts forward ideas that could be very helpful for other researchers in the field, even beyond score-based subfield.\n\nReferences:\n\n[1] Hybrid Top-Down Global Causal Discovery with Local Search for Linear and Nonlinear Additive Noise Models., Hiremath et al., NeurIPS 2024\n\n[2] LoSAM: Local Search in Additive Noise Models with Mixed Mechanisms and General Noise for Global Causal Discovery., Hiremath et al. UAI 2025\n\n[3] Strong and Weak Identifiability of Optimization-based Causal Discovery in Non-linear Additive Noise Models, Li et al., ICML 2025.", "weaknesses": "1. The empirical analysis is conducted largely on simulated linear settings, with only a few real-world dataset considered. Given that the main contribution of FLOP is strictly better empirical performance (no theoretical guarantees outlining provable improvements are given), a more substantial empirical study of the performance of FLOP in nonlinear ANM and real-world data would help contextualize the contribution of FLOP. This is true even if results are not as promising as in the linear setting, as it might highlight where modifications in future could be made to build on top of FLOP.\n\n 2. Although the authors differentiate FLOP from BOSS at some specific subparts of the algorithm, (for ex. at the end of Section 2 they highlight the BOSS reinsertion strategy to the reader to build upon), it is difficult for a reader who is not already familiar with the BOSS paper to parse. For example, at the beginning it’s not clear whether Algorithm 1 constitutes the entirety of the BOSS algorithm, or just a subpart (it is, which is discussed later). Would be nice for the full algorithm environments of BOSS and FLOP and corresponding walk-throughs to be present (or if necessary signaled to the reader early on and placed in the Appendix) to facilitate easy comparison for the reader. Or perhaps even adding a figure describing the general pipeline of these specific score-based methods could be helpful, for better clarity.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "baGHGowG1l", "reviewer": "Reviewer_cdXL", "summary": "Summary. The authors propose a modification of the BOSS algorithm (Andrews et al., 2023), a score-based approach for Markov Equivalence Class discovery from data generated according to a linear Gaussian model. The authors propose 4 modifications to improve the efficiency and accuracy of the algorithm. They show that their modifications preserve theoretical guarantees of convergence of the method to a solution in the Markov equivalence class. The empirical results support the claim of improved scalability and accuracy.", "questions": "What are all the hyperparameters of the method? In the weaknesses I ask clarifications about $\\lambda_{BIC}$, but are there others? If yes, can the authors clarify how the hyperparameters are set, and empirically analyse the sensitivity of FLOP to changes?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 4, "contribution": 3, "strengths": "The main contribution of this paper is an algorithm with really impressive empirical performance, both in terms of accuracy and scalability.  Theoretical guarantees of identifiability of the Markov equivalence class that are common in score-based approaches are maintained, as the authors demonstrate. Finally, the empirical results are convincingly and well presented. Figure 2 clearly shows improvements in the run time. Figure 1,3,4 clearly show the improvements in terms of accuracy due to the ingredients they add to their algorithm. Nice work!", "weaknesses": "1. **Real world  / beyond linear Gaussian experiments.** The proposed method achieves impressive performance on linear Gaussian data, which are covered by the theory of the algorithm. However, the scope of these models is quite limited. It would be relevant to experiment beyond linear Gaussian synthetic data (at least nonlinear ANM + vary the noise distribution — analytic or NN generated nonlinearities, e.g., see [Montagna et al. 2023](https://causally.readthedocs.io/en/latest/generated/causally.scm.causal_mechanism.NeuralNetMechanism.html#causally.scm.causal_mechanism.NeuralNetMechanism)), and also with real-world datasets (that’s where causal discovery is interesting at the end of the day). E.g. you can consider the following datasets that are easy to find: bnlearn networks considered in the paper (please use real data when available), Sachs, Lucas, Dream (relevant for high dimensions), and the Syntren synthetic data generator. \n2. **Hyperparameter tuning**. The authors set $\\lambda_{BIC} = 2$. It is, however, unclear whether this is optimized for the analysed datasets, and how sensitive FLOP is to this hyperparameter. Could the authors discuss that, and present experimental results at different values of $\\lambda_{BIC}$?\n3. **Incremental nature of the contribution**. It is worth noting that the new method consists of smart but incremental additions to existing methodologies — this doesn’t change my opinion on the positive value of the work. In this sense, I am curious about whether the *warm start* tweak of Section 3.1 and the use of ILS in Section 4.2 are novel in causal discovery, or whether the authors have a reference to past use of these ideas. \n\nI am happy to raise my score upon satisfactory execution of the required experiments!", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "aiiNAhVbk3", "reviewer": "Reviewer_ZiZP", "summary": "This paper discusses the acceleration of score-based causal discovery algorithms. Traditionally, estimating causal DAGs from observational data using score-based criteria (e.g., BIC) has been regarded as computationally intractable due to the NP-hardness of discrete search, raising concerns about scalability in large datasets. In recent years, continuous optimization methods such as NOTEARS and DAGMA have emerged as popular alternatives. However, these approaches suffer from challenges including non-convex optimization, approximate enforcement of acyclicity constraints, and the need for post-hoc thresholding.\n\nThe paper proposes a new score-based causal structure learning algorithm, FLOP (Fast Learning of Order and Parents), which dramatically improves the efficiency of discrete search and achieves better performance than existing continuous optimization approaches. FLOP targets causal DAG estimation under a linear additive noise model (linear ANM). It introduces parent set reuse, where the parent sets of each node are initialized from the results of previous searches, reducing both computational and memory costs. A modified grow-shrink procedure is proposed and theoretically proven to correctly identify the Markov boundary. In addition, dynamic Cholesky updates enable efficient computation of conditional variances in the BIC score by reusing Cholesky decompositions, reducing the cost from $O(k^3)$ to $O(k^2)$. The method also incorporates principled order initialization, which constructs an initial variable order by placing highly correlated nodes adjacent to each other, and Iterated Local Search (ILS) to substantially improve the accuracy of discrete optimization.\n\nExperimental results on synthetic data demonstrate that FLOP runs more than 100 times faster than conventional discrete search algorithms while achieving nearly the same accuracy as exact search methods. Moreover, it outperforms continuous optimization methods, highlighting the renewed importance of discrete optimization in score-based causal discovery.", "questions": "Could the authors discuss situations in which each of the proposed techniques, namely, parent set reuse, dynamic Cholesky updates, principled order initialization, and iterated local search, might fail to perform as expected?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "While discrete search–based score methods have often been overshadowed by the recent rise of continuous optimization approaches due to their NP-hardness, this paper makes a significant contribution by demonstrating that discrete search can be made dramatically more efficient through theoretically grounded and well-designed algorithmic improvements. The proposed method introduces principled modifications to key components of the discrete search process, resulting in substantial computational gains.\n\nThrough simulation experiments, the authors show that their method runs over 100 times faster than conventional discrete search algorithms, achieving nearly the same accuracy as exact search methods. Moreover, it also outperforms representative continuous optimization methods. Overall, this work makes an important contribution to the field by reestablishing the relevance and potential of discrete optimization in score-based causal discovery, a domain recently dominated by continuous optimization approaches.", "weaknesses": "While this paper makes a valuable contribution by demonstrating that score-based causal discovery with discrete search can be made computationally practical through well-engineered algorithmic improvements, its theoretical and methodological novelty is somewhat limited. Many of the proposed components, such as parent set reuse, dynamic Cholesky updates, and local search, are not fundamentally new, but rather effective adaptations of existing techniques.\n\nMoreover, although each of the algorithmic enhancements is well-motivated and complementary, their effectiveness may depend on specific data conditions. For example, the parent set reuse and dynamic Cholesky updates assume relatively stable local structures and well-conditioned covariance matrices; the principled order initialization relies on correlations being informative about causal directions; and the iterated local search may require tuning depending on the ruggedness of the score landscape. A more explicit discussion of these limitations and potential failure cases would further strengthen the paper’s credibility and generality.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "v1gGi5I8oX", "reviewer": "Reviewer_UFtS", "summary": "This paper introduces FLOP (Fast Learning of Order and Parents), a novel score-based, discrete search algorithm for causal structure learning for linear additive noise models (ANMs). FLOP builds upon an order-based search framework, similar to the BOSS algorithm, but introduces four key components that lead to significant gains in both speed and accuracy: warm-start parent selection, dynamic Cholesky updates, principled order initialization, iterated local search. Experiments across a wide range of benchmarks (ER, SF, Alarm, etc.) demonstrate that FLOP variants achieve state-of-the-art accuracy (measured by SHD and AID) while being orders of magnitude faster than competing methods like BOSS, GES, and DAGMA.", "questions": "N/A", "limitations": "", "rating": 2, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 1, "strengths": "1. This paper is well-written and easy to follow.\n\n2. By combining warm-start parent selection with dynamic Cholesky updates, FLOP achieves massive computational speed-up. Furthermore, this speed-up enables more exhaustive search, leading to higher accuracy.\n\n3. The performance of FLOP is validated across many benchmarks.", "weaknesses": "1. The authors do not provide theoretical analysis on finite-sample performance of FLOP. I think this is a major weakness of this paper. First, the authors repeatedly emphasize the importance of finite-sample performance, e.g., \"Finite-sample performance is a practical challenge, despite the mature identifiability theory and asymptotic guarantees\", \"one of the core issues of score-based methods in practice are finite-sample induced local optima\". Second, the authors consider the simplest structural causal model where structural functions are linear, exogenous noises are Gaussian, and there is neither latent confounder nor selection bias.\n\n2. The performance of FLOP is only validated on synthetic data. Although the authors consider real-world networks from the bnlearn repository, the data is still synthetic. This is not convincing for a paper whose primary contribution does not lie in theory.\n\n3. The ILS perturbation strategy is defined as $k = \\ln p$ random swaps of two elements. The paper simply states this is \"robust\" without an ablation study or sensitivity analysis.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)