Differentiable Structure Learning and Causal Discovery for General Binary Data
Authors: Chang Deng
Bryon Aragam
Venue: NeurIPS 2025 poster ()
Keywords: differentiable structure learning
causal discovery
discrete data
binary data
DAG
Primary area: probabilistic_methods
TL;DR: We propose a differentiable structure learning framework for general binary data that makes no parametric assumptions about the data‐generating process.

Abstract: Existing methods for differentiable structure learning in discrete data typically assume that the data are generated from specific structural equation models. However, these assumptions may not align with the true data-generating process, which limits the general applicability of such methods. Furthermore, current approaches often ignore the complex dependence structure inherent in discrete data and consider only linear effects. We propose a differentiable structure learning framework that is capable of capturing arbitrary dependencies among discrete variables. We show that although general discrete models are unidentifiable from purely observational data, it is possible to characterize the complete set of compatible parameters and structures. Additionally, we establish identifiability up to the Markov equivalence class (MEC) under mild assumptions. We formulate the learning problem as a single differentiable optimization task in the most general form, thereby avoiding the unrealistic simplifications adopted by previous methods. Empirical results demonstrate that our approach effectively captures complex relationships in discrete data.

Differentiable Structure Learning and Causal Discovery for General Binary Data
Bryon Aragam Booth School of Business, University of Chicago, Chicago, IL 60637
Chang Deng"*
Abstract
Existing methods for differentiable structure learning in discrete data typically as- sume that the data are generated from specific structural equation models. However, these assumptions may not align with the true data-generating process, which limits the general applicability of such methods. Furthermore, current approaches often ignore the complex dependence structure inherent in discrete data and consider only linear effects. We propose a differentiable structure learning framework that is capable of capturing arbitrary dependencies among discrete variables. We show that although general discrete models are unidentifiable from purely observational data, it is possible to characterize the complete set of compatible parameters and structures. Additionally, we establish identifiability up Markov equivalence under mild assumptions. We formulate the learning problem as single differen- tiable optimization task in the most general form, thereby avoiding the unrealistic simplifications adopted by previous methods. Empirical results demonstrate that our approach effectively captures complex relationships in discrete data.
1 Introduction
Causal relationships are often represented by directed acyclic graphs (DAGs), where nodes correspond to variables and directed edges indicate cause-effect links [33, 42, 35]. Learning a DAG from observational data (structure learning, also known as causal discovery) well-known NP-complete problem [8, 10]. Recent advances have recast this combinatorial search as differentiable constrained program, enabling the use gradient- baased methods for structure learning [54, 55]. In this approach, the adjacency matrix of the graph is treated as continuous variable and one optimizes a score (e.g. negative log-likelihood) subject to a differentiable acyclicity constraint that ensures the solution is valid DAG [54, 5]. Early work in this area has largely been limited to continuous data and relied on specific structural equation models (SEMs) example, linear or nonlinear functions with additive Gaussian noise [36]. While such assumptions facilitate tractable optimization, they may misrepresent the true generative process for non-Gaussian or discrete data, leading to biased or inconsistent structure learning.
Many real-world datasets involve binary or other discrete variables (e.g. presence/absence of condition, binary genetic markers, survey responses), whose complex dependency structures are not well handled by methods designed for continuous data. Traditional Gaussian [53, 17] or continuous assumptions often break down in these cases. However, only few works have attempted to extend differentiable structure learning to discrete data. These existing approaches typically impose specific parametric forms on the data-generating process. For instance, recent method [49] for discrete differentiable structure learning assumes generalized linear SEM. Moreover, the theoretical identifiability guarantees of these methods rely particular distributional constraints that may not hold for general discrete data. These gaps call for more flexible approach that can capture the rich dependence patterns in discrete data without relying untested generative assumptions.
"[changdeng,
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Outside the differentiable paradigm, decades of work have produced both score-based [e.g. 9] and constraint-bar -assed [e.g. 41] algorithms for causal discovery. Many discrete methods also impose additional simplifications (additive-noise models [34], hidden compact representations [7], linear effects [27], latent Gaussian variables [1]), which may not apply to díscrete data. Even the most recent differentiable approaches evaluated on discrete data lack formal justification: For instance, Bello et al. [5] apply continuous optimization methods to discrete data under a logistic regression model, although the underlying identifiability and theory for such a model has yet to be worked out. These observations underscore the need for a general, theoretically-soun framework for differentiable structure learning in discrete data.
In this work, we propose a general differentiable structure learning framework for discrete data, which captures arbitrary dependencies without assuming specific data data-generating process. We first show that, in this general setting, DAGs are non-identifiable from observational data alone. Nevertheless, we characterize the complete set of compatible structures and parameters. To select among these, we adopt a sparsity principle, aiming for the sparsest DAG consistent with the observed distribution. We formulate this single differentiable optimization problem. Lastly, we establish theoretical guarantees showing that our method recovers the correct Markov equivalence class under mild assumptions. Unlike prior work, our approach avoids assumptions like linearity or additivity, enabling it to model richer causal dependencies.
In summary, our key contributions are as follows:
1. We start with general binary data, without assuming a particular data generating process, and use the multivariate Bernoulli distribution (Definition 1) general representation for any such model. We prove that, in this fully general setting, the underlying DAG is con-identifiable and characterize which graph-parameter pairs are compatible with the observed distribution (Theorem 1).
2. We formulate learning of the sparsest graph and its parameters single differentiable program (13). We show that even without assuming faithfulness, any global minimizer yields the sparsest valid DAG (Theorem 3). Moreover, based our framework, we provide theoretical justification for prior works [13, 5].
3. We conduct experiments demonstrating that our method can capture general causal rela- tionships between discrete variables. Furthermore, we introduce a two-stage approach, BINOTEARS, which reduces computational complexity and outperforms existing base- lines.
2 Related work
In this section, survey existing approaches to structure learning and causal discovery, highlighting methods for discrete data, differentiable DAG learning in continuous settings, and recent extensions of continuous-data techniques discrete domains. We discuss their key ideas, assumptions, and limitations, which motivate the need for a general framework for multivariate discrete data.
Traditional structure learning for discrete data Classical causal-discovery methods for discrete data fall into two main categories. Constraint-basee algorithms (e.g. PC [41], MMHC [43], Copula- PC [11]) use conditionalinndependenee tests to infer edges. Score-based approaches assign each candidate DAG a goodness-of-fit score (e.g. BIC, BDeu [29, 20], generalized scores [22, 27]) and then search for the optimal structure. Greedy strategies such as GES [9] are known to be NP-hard [8, 10] and can get trapped in local optima [26, 2]. Moreover, score-based search traditionally relies on strong parametric assumptions- e.g. additive-noise models [34], hidden compact representations [7], linear effects [27], latent Gaussian variables [1]-which may not hold in practice. recent exception for general nonparametric models is [3]. A few methods target bivariate causal discovery [45, 28, 6]. but they do not to multivariate problems. As result, existing discrete causal-dis -iicooedddyy techniques either lack statistical robustness in finite samples or impose restrictive assumptions, motivating the need for scalable, assumption-lige methods that handle dimeediiensional binary data.
Continuous DAG learning for continuous data Since the introduction of smooth characteriza- tion of acyclicity by Zheng et al. [54], which allows treating the adjacency matrix continuous

---

optimization variable, a large body of work has emerged in differentiable structure learning. Subse- quent research has extended this framework in several directions: Nonlinear models [55, 48, 25, 56], optimization theory [44, 32, 14], designing more stable and computationally efficient acyclicity constraints [5. 30, 48, 50, 51], understanding score function properties [31, 15, 40], and incorporating extra side information [4]. Despite this progress, most existing methods are fundamentally designed for continuous variables and do not address the specific challenges of discrete data. The common reliance on Gaussianity or continuous-varian variable assumption--eeadiie to objectives like mean--quured error or Gaussian log- discrete data.
potential biases and inconsistencies when applied to
Continuous DAG learning for discrete data Only handful of recent studies have extended differentiable structure learning to discrete data by imposing specific parametric forms on the conditionals. For example, most works [49, 5, 13] assume a linear SEM in which each discrete node's probability is given by logistic link on a linear combination of its parents. Such restrictive assumptions can fail to capture the rich, higher-order dependencies present in general discrete distributions. Another line of work in differentiable structure learning studies nonlinear models via neural networks [55, 25, 56], which in principle could accommodate discrete inputs. However, these methods treat inputs as generic real-valued vectors and do not formally address the unique features of discrete data, neither in theory nor in empirical evaluation. A recent contribution by Deng et al. [15] proposes general differentiable framework for parametric families, but the necessary details to connect generic parametrizations of discrete data to structure learning (in particular, the adjacency matrix), are missing. By contrast, one the of the key contributions in this work is to develop tractable parametrization via the multivariate Bernoulli distribution that allows explicit representation of the adjacency matrix in a fully differentiable framework without restricting the class of discrete distributions.
3 Preliminaries
Let G (V, E) denote a directed graph on p nodes, with vertex set V = [p] == {1, ,., p} and edge set E C V x V, where (i.j) E E indicates the presence of directed edge from node to node j. We associate each node i V to a random variable X,, and let X (X1..... .,. .pp). In this work we focus on discrete random variables. Any variable with T categories can be represented by (T 1) binary indicators [46]. For instance, if Z E (0. 1, 2). then introduce z* = 1 1 k). k = 0, 1, indicating whether Z for & 0, 2, so that Z is replaced by three binary indicators. Hence, without loss of generality, we only consider binary data in the sequel. Extending the results to general discrete data only needs additional notation and does not require new technical ideas.
We review the multivariate Bernoulli distribution (MVB) [12], a flexible family that captures all joint dependencies among high-din sigh-dimensional binary variables. Because it assigns a probability to every possible configuration, allows us to study binary data without imposing additional assumptions.
Definition 1 (Multivariate Bernoulli distribution [12]. Let X (X1, (Xi.... Xp) correlated Bernoulli random variables with X; E {0, for i = =,...,, We say that X follows multivariate Bernoulli distribution with p E R²?", written X ~ MultiBernouli(p), if
be a vector of possibly
P(X x) II p(1s) my 1lzes(1-x)), SCIp
where 1s E f0, 13 for the vector with ones on and zeros elsewhere and each entry of p = (p(0, ........, 0). ,(...).,((........... 1)) E RR" is the probability mass of a distinct configuration. Such called the general parameter. Let 12x E gg2 denote the all-ones vector; normalization requires 1 1.p = 1.
pl.....,,11
Based on (1), is not clear how to determine the dependence between each variable. For this, define (x)x IIes Tyy for any S C [p] with convention B (x)) 1. Rewrite (1) in exponential-family
form:
P(X x) exp
Σ S B²(x) Sch2
(2)
3

---

where the coefficient { are the natural parameters. The corresponding natural parammeer vector is f= ( ...., fr ,...,f..... fl...)) F²a Throughout, the superseript set S in f is treated as unordered; e.g. for [1,2] we have f'2 f Because the MVB is an exponential family, there is one-to-one correspondence between the natural parameters f and the general parameter p: explicit conversion formulas are provided Appendix C.I. Expressing the model as in (2) will later allow us to identify variable dependencies directly from the coefficients f when reconstructing causal graphs. A simple but useful closure property follows immediately.
Corollary 1. If X MultiBernouli(p)) then every marginal distribution and every conditional distribution of X is again multivariate Bernoulli.
4 General discrete data: A nonidentifiable model
We first establish that, for fully general binary data, neither the causal structure nor the associated parameters can be identified from a single observational environment [19, 23]. The intuition is simple: For any fixed topological sort' of the variables there exists unique causal structure and associated parameters that reproduces the observed distribution exactly. Since there are p! distinct topological sorts, this renders the model fundamentally non-identifiable.
4.1 Conditional distributions and higher order interaction
First, we provide a useful expression for conditional distributions in this model. Given any permutation on V, the joint distribution admits the standard factorization
P(X) = II P(X=)) X (x()))..xG-1)))
Let us use the joint law P(X) given by the exponential form in (2). Without loss of generality let (P. and focus on j = p; other choices of π and j lead to the same algebra with heavier notation. Then
Pr(Xp = 1| X = x-p) =0
fr + f fl. t... ff. fn-1.P 1p-1ppp-' f 112.0 x1 x1 +... fl 11. I1...
Σ
SSP (x)x)
sC|p-1]
where [p P.....PP 1} and (=)) 1/(1=xxp(=2)) The full derivation appears in Appendix C.2. Equation (3) highlights that, for general binary data, all higher order interaction terms are present [12, 16]. Omitting these higher-order interactions-c.g., restricting attention first-order (additive) effects [47, can misrepresent real data.
4.2 Nonidentifiability and Equivalence
In the example above, Xxp is the last node in the ordering so any variable in Xp-1) may serve as a potential parent of Xp- A variable X, is deemed a parent precisely when some interaction coefficient involving and p is non-zero; equivalently,
X, + Xp + 3 C Ip- - 1] with S, such that SSP 00
Σ (sSw)²>
0. (4)
8cIe-11,j88
All coefficients in (3) can be estimated simultaneously via a single logistic regression [24]; deciding
whether an edge X, + Xp exists thus reduces to checking if the corresponding coefficient set is non-zero. Because (3) involves higher-order interaction terms, we introduce the following notation.
Let 21rl denote the power set of the index set [p] (1..... p). For any vector X (X1,.... Xp) RPP and subset S C (p), define the interaction feature B (X) Iljes X, with Bx( 1.
A topological sort permutation Vf such that Xa(s) Xa i<j.

---

Collecting all 2" such terms yields the extended feature map: (X) = [B (X)
Se2 E R²"
ordered according to the lexicorriico rule described in Appendix C.3. Explicitly,
((X) [1. XI..... Xp: X1X2,.... Xp-1Xp= X1 XI...X, Xp]
(5)
For example, when p- and X (x,, X3). then the extended feature vector (X)
[1. X1,X11X22 X3, ,, X2 X E R2. When & is applied to data matrix X E
RExp, operates row-wise.
For each j and order n, define the relevant natural parameters in the vector fx. E R22 stored in the same -rxioo-lexicoogaphic order, but relative to the (Xx(1)....; Xx6-1)) Concretely,
f fI
(x))
p (1),x())
(1)
(x-11)
1)
(1()))))
(1 .t x)) ()-))))))
(6)
n
constant
(2)
X (1--1)
X X
Each underbrace indicates the interaction term either constant, a single feature. or product of features from (Xx(I)...... (.ai))....,,..-
that the coefficient corresponds to. Thus, the elements of f n.
align exactly with the entries of the vector (((x())....xxt--)) The special case used in (3) is recovered by taking = (I, and = p.
With these definitions in place we can now recover the full causal structure and its param- eters. Fix an ordering π and consider any j E [p]. By Corollary 1, both the marginal P(Xx(j), Xx()) Xx()-)) and conditional P(Xx(j) xx())...,()-= remain multivariate Bernoulli. Let X (X....., Xp) E R2xp be i.idd samples from X ~ MultiBernouli(p). one sample per row. Hence we can apply the same logic as in Section 4.1:
Form extended feature vector O (((x(1)) x))) and run single logistic regression of X=()) on these features to recover natural parameters fx.j E g2--
Apply the edge criterion in (4) to F.j to determine the parent set PA(T(i)).
Repeat steps above for [p] to construct the DAG G, and parameter f = (
Iterate over all permutations to enumerate every graph-paramete pair.
Complete details on such procedure can be found in Appendix C.4.
Theorem 1. Let p and fix any topological sort æ. Consider the population case (n + 0o). Let (fz, Gx) be the output of procedure above. If Y E IP are generated by the following structural equations
Y=() Bernoulli (o (FT, (
Vj = 1 L.....P
where (2) 1/(1 + exp(-2)), then Y MultiBernoulli(p)
We work the population case to eliminate finite sample estimation error. The assumption p 0 is mild regularity condition routinely adopted in structure learning problems [49, 36, 21] and guarantees that the logistic regressions used in the procedure above always have unique solutions.
Remark 1. Equation (7) expresses a general discrete model as a structural equation model (SEM)), framework widely used for structure learning because it provides generative description of how each variable arises from its direct causes. A central appeal of SEMs is their universality: Every probability distribution can, in principle, be represented by a suitably chosen SEM (Proposition 7.1 in 35). practice, however, SEM-based methods almost always impose strong parametric formss--foomm-o-ee forms-linear additive, or low-order interactions [47, 49, 18, 52]-which risk ignoring important structure in real data [12, 16]. The problem is particularly acute for discrete data, where higher-order dependencies are easily overlooked, sharply limiting the model's 's expressive power. By contrast, Theorem I shows that general binary distributions necessarily involve higher-order interaction terms. Omitting these terms therefore precludes exact representation and may lead to incorrect causal conclusions.
Thus, the Theorem indicates that, for every topological order и we can construct a distinct SEM that reproduces the same distribution. The model is therefore non- non-ide non-identifiable: multiple parameter-grap
Formal details can be found in Appendix C.5
5

---

pairs (frr: G+) obtained by procedure before give rise to the identical law MultiBernoullipp) Collect all such pairs into the equivalence class
{p) = ((S GT) Gz)
The cardinality of {p) at most pl, corresponding to the number of permutations of p variables. Remark 2. It is well-known that causal DAGs are non-identifiable from observational data alone [20, 19, 23). However, Theorem I provides a novel and comprehensive characterization of all DAG- parameter pairs that exactly reproduce given distribution within the general mullivariate Bernoulli framework.
4.3 Minimal equivalence class
Because observational data cannot distinguish among members of {p)), we seek the simplest representation. Specifically, the DAG with the fewest edges. Let 8G denote the number of directed edges in graph G; our objective is to find the pair (f Gx) E(p) with min SC..
Definition (Minimality [15]). (fx: Gx) is the minimal element in equivalence class E(p) if
SG, S SG+, V(fx, G+) E Ep)). The set of all the minimal element in equivalence class E(p) is referred to as the minimal equivalence class Emiu(p)
(p)n(p) = : G+)= (fx: Ga) is the minimal element, (J Ga) E E(p))
(8)
For the procedure in Section 4.2, we retains those pairs (fx. with the fewest edges, thereby recovering Emin (p). Minimality is closely related the Sparsest Markov Representation (SMR)
assumption [38, 26), which posits that, among all DAGs compatible with a distribution, the sparsest one is unique up Markov equivalence, i.e., encoding the same conditional independence relationship. SMR is weaker than the well-known faithfulness assumption [35]: if distribution P is faithful to DAG G, then (G, P) satisfies SMR. Hence, targeting the sparsest graph both principled and practically appealing. Let M(G) denote the Markov Equivalence class of a DAG G, i.e. the set of all iaarriiona--depedddeee relations. Formal definitions of these concepts
DAGs encoding the same
appears in Appendix C.6.
Theorem For any p > O and consider the population case (n + c0). Under the SMR assumption (or faithfulness assumption), if (fxx, Gm, ), (f m. Gra) E (p)n(p), then M(G,,) M(G...)
5 Continuous structure learning for general binary data
The procedure in previous section is purely combinatorial and thus fails scale: It requires fitting p! logistic regressions, so the computational cost grows exponentially with the number of nodes To overcome this bottleneck, we leverage recent advances in differentiable structure learning [54, 55]. Specifically, we recast the search for the sparsest binary-dar data DAG as a single differentiable program, which can then be tackled by standard gradient- optimizers.
Parameterization Let X = (X,,.... Xp) be n i..dd samples from X MultiBernoulli(p), one sample per row. Detine the parameter matrix
H
=( h U.j
12. hP.. 1225
h (P-11)p.) 1.p.
128...p.j)
E 12
(9)
first order
second order
third
1)-th order
p-th order
Let H = (H,,.. Hp) E 2 Here H, plays the role of column in the weighted adjacency matrix. Unlike linear model, where single entry Wiy signals the edge X, - X, [54], the multiple coefficients in H, jointly determine whether such an edge exists.
Weighted adjacency matrix
Define the induced adjacency matrix
[W (H)) =
Σ (h&')²
(10)
scipliie
so [W (H)), > 0f and only if some interaction involving X, contributes to the equation for X,. that some coefficient S, with E non-zero. Self-loops are forbidden, so we impose
h = whenever j E S, Vje lp), S S [p].

---

Score function By Theorem the negative log- ligelikei of the multivariate Bernoulli model- i.e. the logistic or cross entropy loss provides suitable score. Applying the extended feature map (5) row-wise to the data matrix X yields the score function
{(;
I (log(1n + exp(*(X)H)) - X, (*(X)H)
where denotes the Hadamard product and exp is applied element-wise see Appendix C.7 for details.
Regularization While {HH; X) identifies the equivalence class Epp). our objective is a minimal equivalence class Emi (p). Simply adding & [20] penalty would break differentiability. and an & penalty is both biased and imprecise in edge counting. Instead, we adopt the smooth quasi minimax-concave penalty (MCP) [15]:
[(1
t² 28
quasi-MCP:
pt,(())
( < +)+ + (1d) > 5)
(11)
where I(-) is the indicator function. Px.5 is quadratic on [0. 6] and flat beyond 6, smoothly approxi- mating an to penalty and thereby encouraging sparse W/ (H) without sacrificing differentiability.
Let Px, (H) = EEiz (H))(()))) Define the regularized score
(H;2, = ((H; X) + Px.6 ( (H))
(12)
Recovering Emin (p)
We formulate this task as the single continuous optimization problem
min H
s(H; X, 6)
subject to
h W (H) =0 Su if VjE pl..
(13)
where Rexp + [0, oc) is the differentiable acyclicity constraint [54, 5, 30, 48], satisfying h(W) 0 iff W is a DAG. To study the optimal solutions of (13), define the set of global minimizers
O,,2,5 ( GG(((H") H" minimizer of (13))
(14)
where G(W) denotes the graph encoded by the adjacency matrix W via (10). Although (13) optimizes only over H, we include the induced graph to match the notation of & {p)) and Emin (p).
Ideally. we want the optimal solution set to coincide with the minimal equivalence class, i.c. O,,A.5 Emin (p). However, it is unclear whether suitable values of (A, exist. The next result shows that, in the population limit, such values can always be chosen. Let Oc,,A,5 denote the set of minimizers of (13) when the empirical loss s(H; X, ó, X.) is replaced by its population counterpart E[s(H; X, 6, x.
Theorem 3. For any p > there exist X, 0 sufficiently small such that Ox0 Emin (p). Moreover, under Sparsest Markov Representation (or faithfulness) assumption, for any (Zz, G.y), (fm,, GTz)) Oxc,2,6, s, holds M(Gm.) = M..
Here, M(G) is the Markov equivalence class of G. The proof adapts Theorem of Deng et al. [15], but requires additional technical work accommodate the multivariate Bernoulli setting: The theorem implies that, with appropriately small regularization, every global optimum (13) recovers sparsest DAG and preserves all conditional-indeeedeeee relations.
Connection to prior continuous methods on discrete data. Existing empirical studies [13. 5] often generate binary data by passing linear combination of parent variables through logistic link, then apply continuous-open cotruuuuuoootimzzztionn techniques without formal justification. Theorem provides that justification: Such simulations instantiate special case of our general model. To formalize this setting, we characterize this with the following assumption.
Assumption A (First-order logistic model). There exist DAG Go such that X (X1,.... X,) are generated according to the following structure equation model
Xj XPA()) Bernoulli (a (w] X + o;))
(15)
where XPAC) are the parents node of X, in GO, Here [w)lb))AG)] 0, (WWW)PAAG) #0,c, R.
7

---

Under Assumption A, only first-order interaction terms are present; all higher-order coefficients in (9) vanish. Consequently H shrinks from 2 p to (p + 1) x p. and the optimization problem becomes tractable even for moderate p. A detailed formulation and corresponding consistency theorem are provided in Appendix C.8.
Experiments
We solve (13) using either DAGMA [5] NOTEARS-MIP [55]. Both methods tackle a sequence of unconstrained problems by incorporating the acyclicity constraint as penalty term with an increasing weight. Our method-denoted BINOTEARS- compared against several baselines, including DAGMA [5], PC [41], and FGES [37]. Our primary empirical results appear in Figures 1 and 2. We evaluate accuracy using structural Hamming distance (SHD), which measures the discrepancy between the estimated and ground-truth graphs, with lower values indicating better performance. Because the multivariate Bernoulli model nonidentifiable, we compare completed partially directed acyclic graphs (CPDAGs) of the estimates and the ground truth. Full experimental details are provided in Appendix D.
In Figure 1 (a), we simulate data according to the SEM in (7) for small DAG (d < 10), including both first-order effects and second-order interactions among parents whenever node has multiple parents. We observe that BINOTEARS (ours) achieves competitive performance relative to the baselines. In contrast, DAGMA which models only first-order effects, exhibits markedly poor recovery when second-order interactions are present. These results underscore that a purely linear form is insufficient to capture the complex dependencies in general discrete data and can lead to substantial estimation eπтоrs whenever higher-order terms play role. In Figure (b), we repeat the simulation while ablating all first-order effects for nodes with multiple parents, retaining only the highest-order orddrr interaction term. Under these conditions, BINOTEARS (ours) consistently attains near-optimal recovery performance across all settings, showing the robustness of our method.
In Figure 2, we consider a larger DAG, where including all possible interaction terms would cause the feature map ((X) to grow exponentially and render the optimization in (13) intractable. To address this, our ANOTEARS works in two-stage fashion. First, we adapt the original NOTEARS-MLA framework of Zheng et al. [55] to handle discrete data, and use it learn an initial graph. From this graph we extract a topological ordering π. In the second (finetuning) stage, we follow the procedure of Section 4.2: we construct only those interaction features consistent with π,, and then logistic regressions with first order and second order to recover the final edge structure G. BINOTEARS attains performance on par with existing baselines. While all methods degrade as graph size and density increase, our approach exhibits robustness and consistently recovers meaningful structures.
We further evaluate NOTEARS (linear), NOTEARS-MLP. and BINOTEARS (ours) the real- world dataset of Sachs et al. [39]. The data contain n = 7,466 measurements of d = 11 proteins and phospholipids in human immune cells and are accompanied by widely used consensus network that serves as gold standard. Results are summarized in the Table. Further details on the experimental setup and additional results be found in the Appendices D and E.
Method
SHD
Num. of Edges
NOTEARS [54]
22
18
NOTEARS-MLP [55]
16
13
BINOTEARS (Ours)
13
15
7 Conclusion
We have introduced fully general differentiable framework for structure learning in binary data, based on the multivariate Bernoulli distribution. Our key contributions include (i) a complete characterization of non-identifiability in the general discrete setting and constructive procedure
'Concretely, we augment DAGMA with higher-order terms on small graphs to pursue exact recovery, and (ii) employ NOTEARS-MIL MML to obtain a topological order on large graphs and subsequently fit logistic regressions for edge selection. We use these two pathways because, in our experiments, they strike different ppeedaaccuran trade-offs. See Appendix D.2 for implementation details.
8

---

7.5
Hamning Desonce (S
Stucura Panning Desonce
3.0
0.8
10.0
15
Soucural
10.
7.5
3.0
2.5
0.0
Muumber maden)
strodes)
6luumber nondes)
DAEMA
BINOTEARS 10urs] FEES FC
DAGMA
BINOTEARS (Ours)
FGES - PC
(a) First order and second order interaction are involved
(b) Only highest order interaction is involved
Figure I: Results in terms of SHD between MECs of estimated graph and ground truth. Lower is better. Column: & 21. Row: random graph types. (ER,SF)-k [Scale-Free, Erdōs-Rényi] graphs with k expected edges. Here p (5, 6, 7, 8, 9). Error bars denote the standard error computed over 10 replications.
125-
30
Homning Detorce (SHD)
20-
27
75
10
60-
60
80
40
20
20-
010
(INumaertt (nddess)
iNuubbrrtt
[Wuumbrro of
Methel
DAGMA
FGES
BINOTEARS (Ours)
PC
Figure 2: Results in terms of SHD between MECs of estimated graph and ground truth. Lower is better. Column: kt 2, 4]. Row: random graph types. (ER.SF)-k (Scale-Free, Erdōs-Rényi) graphs with kd expected edges. Here {10, 20, 30, 40). BINOTEARS is our two-stage approach.
to recover the minimal equivalence class (Theorems 1-2), (ii) fully differentiable constrained programming that provably recovers sparsest DAG up to Markov equivalence (Theorem 3), and (iii) empirical validation on synthetic graphs of varying size, demonstrating that our method captures complex higher-order dependencies where existing approaches fail.
To scale beyond small graphs, we further propose a practical two-stage heuristic, BINOTEARS, which first learns coarse structure via NOTEARS-MLE adapted to discrete data and then refines edge estimates by fitting logistic regressions along the inferred topological order. While this approach yields strong numerical performance on larger networks, its theoretical properties remain unexplored. In future work, we aim to develop principled, scalable algorithms for discrete structure learning that retain both computational tractability and rigorous identifiability guarantees in regimes with thousands nodes.

---

References
[1] Andrews, B., Ramsey, and Cooper, G. F. [2019], Learning sigh-dimensiona directed acyclic graphs with mixed data-types, in 'The 2019 ACM SIGKDD Workshop Causal Discovery', PMLR, pp. 4-21.
[2] Andrews, B., Ramsey, J., Sanchez Romero, R., Camchong, J. and Kummerfeld, E. [2023], 'Fast scalable and accurate discovery of dags using the best order score search and grow shrink trees", Advances in neural information processing systems 36, 63945-63956.
[3] Aragam, B. [2024], 'Greedy equivalence search for nonparametric graphical models', arXiv preprint arXiv:2406. 17228
[4] Ban, T., Chen, L., Wang. X., Wang, X., Lyu, D. and Chen, H. [2024], 'Differentiable structure learning with partial orders'. Advances in Neural Information Processing Systems 37, 117 117426- 117455.
[5] Bello, K., Aragam, and Ravikumar, [2022], 'Dagma: Learning dags via m-matrices and a log-determinant acyclicity characterization" Advances in Neural Information Processing Systems 35, 8226-8239.
[6] Budhathoki, K. and Vreeken, J. [2018], Accurate causal inference on discrete data, in '2018 IEEE International Conference on Data Mining (ICDM),, IEEE, pp. 881-886.
[7] Cai, R., Qiao, J., Zhang, K., Zhang. Z. and Hao, Z. [2018]. 'Causal discovery from discrete data using hidden compact representation! Advvnccss in neural information processing systems 31.
[8] Chickering, D. M. [1996], 'Learning bayesian networks is np-complete', Learning from data: Artificial intelligence and statistics V pp. 121-130.
[9] Chickering, D. M. [2002], 'Optimal structure identification with greedy search', Journal of machine learning research 3(Nov), 507-554.
[10] Chickering, D. M., Heckerman, D. and Meek, C. [2004], 'Large-samplc learning of bayesian networks is np-hard', Journal of Machine Learning Research 5(Oct), 1287-1330.
[11] Cui, R., Groot, P. and Heskes, T. [2016], Copula pe algorithm for causal discovery from mixed data, in 'Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part 16', Springer, pp. 377-392.
[12] Dai, B., Ding, S. and Wahba, G. [2013], 'Multivariate bernuuli distribution'.
[13] Deng, C., Bello, K., Aragam, B. and Ravikumar, P. K. [2023]. Optimizing notears objectives via topological swaps, in 'International Conference on Machine Learning'. PMLR, pp. 7563-7595.
[14] Deng, C., Bello, K., Ravikumar, P. and Aragam, B. [2023], 'Global optimality in bivariate gradient-based dag learning', Advances in Neural Information Processing Systems 36, 17929- 17968.
[15] Deng, C., Bello, K., Ravikumar, and Aragam, B. [2024], 'Markov equivalence and consis- tency in differentiable structure learning', Advances in Neural Information Processing Systems 37, 91756-91797.
[16] Enouen, J. and Sugiyama, M. [2024], 'A complete decomposition of kl error using refined information and mode interaction selection', arXiv preprint arXiv:2410.11966
[17] Gao, M., Tai, W. M. and Aragam, B. [2022], Optimal estimation of gaussian dag models, in "International Conference on Artificial Intelligence and Statistics', PMLR, pp.
[18] Gu, J., Fu, F. and Zhou, Q. [2019], 'Penalized estimation of directed acyclic graphs from discrete data', Statistics and Computing 29(1), 161-176.
[19] Guo, S., Tóth, V., Schölkopf, B. and Huszár, F. [2022], 'Causal de finetti: On the identification of invariant causal structure in exchangeable data', arXiv preprint 2735522003.15756
10

---

[20] Heckerman, D., Geiger, D. and Chickering, D. M. [1995]. 'Learning bayesian networks: The combination of knowledge and statistical data', Machine learning 20, 197-243.
[21] Hoyer, P., Janzing, D., Mooij, J. M., Peters, and Schölkopf, B. [2008], 'Nonlinear causal discovery with additive noise models', Advances in neural information processing systems 21.
[22] Huang. B., Zhang. K., Lin, Y., Schölkopf,, B. and Glymour, C. [2018], Generalized score func- tions for causal discovery, in 'Proceedings the 24th ACM SIGKDD international conference on knowledge discovery & data mining', pp. 1551-1560.
[23] Huang, B., Zhang, K., Zhang, J.. Ramsey, J., Sanchez- R., Glymour, C. and Schölkopf, B. [2020], 'Causal discovery from helcrogeccooushoonstttinnny data', Journal of Machine Learning Research 21(89), 1-53.
[24] Kleinbaum, D. G., Dietz, K., Gail, M., Klein, M. and Klein, M. [2002], Logistic regression, Springer.
[25] Lachapelle, S., Brouillard, P., Deleu, T. and Lacoste-Julien, S. [2019], 'Gradient-based neural dag learning' arXiv preprint 1006:0:900.0022
[26] Lam, W.-Y., Andrews, B. and Ramsey, J. [2022], Greedy relaxations of the sparsest permutation algorithm, in 'Uncertainty in Artificial Intelligence', PMLR, pp. 1052-1062.
[27] Li, C. and Shimizu, S. [2018], "Combining linear non-gaussian acyclic model with logistic regression model for estimating causal structure from mixed continuous and discrete data', arXiv preprint
[28] Liu, F. and Chan, L. [2016], 'Causal inference on discrete data via estimating distance correla- tions', Neural computation 28(5), 801-814.
[29] Maxwell Chickering, D. and Heckerman, D. [1997], 'Efficient approximations for the marginal likelihood of bayesian networks with hidden variables', Machine learning 29, 181 181-212.
[30] Nazaret, A., Hong. J., Azizi, E. and Blei, D. [2023], 'Stable differentiable causal discovery', arXiv preprint ar 2xxv:23311.1066.
[31] Ng, I., Ghassami, A. and Zhang, K. [2020], 'On the role of sparsity and dag constraints for learning linear dags', Advances in Neural Information Processing Systems 33, 17943-17954.
[32] Ng, I., Lachapelle, S., Ke, N. R., Lacoste-Julien S. and Zhang, K. [2022], On the convergence of continuous constrained optimization for structure learning, in 'International Conference on Artificial Intelligence and Statistics', Pmlr, pp. 8176-8198.
[33] Pearl, J. [2009], Causality, Cambridge university press.
[34] Peters, J., Janzing, D. and Schölkopf, B. [2010], Identifying cause and effect discrete data using additive noise models, in 'Proceedings of the thirteenth international conference on artificial intelligence and statistics', JMLR Workshop and Conference Proceedings, pp. 597-696-44
[35] Peters, J., Janzing, D. and Schölkopf, B. [2017], Elements of causal inference: foundations and learning algorithms, The MIT Press.
[36] Peters, J., Mooij. J. M., Janzing, D. and Schölkopf, B. [2014], 'Causal discovery with continuous additive noise models', JMLR.
[37] Ramsey, J., Glymour, M., Sanchez-Romero, R. and Glymour, C. [2017], 'A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application functional magnetic resonance images', International journal of data science and analytics 3, 121-129.
[38] Raskutti, G. and Uhler, C. [2018), 'Learning directed acyclic graph models based on sparsest permutations', Stat 7(1), e183.
[39] Sachs, K., Perez, O., Pe'er, D., Lauffenburger, A. and Nolan, G. P. [2005]. 'Causal protein- signaling networks derived from multiparameter single-cell data", Science 308(5721), 523-529.
11

---

[40] Seng, J., Zečević, M., Dhami, D. S. and Kersting, K. [2024], Learning large dags is harder than you think: Many losses are minimal for the wrong dag, in "The Twelfth International Conference on Learning Representations.
[41] Spirtes, P. and Glymour, C. [1991], 'An algorithm for fast recovery of sparse causal graphs', Social science computer review 9(),, 62-72.
[42] Spirtes, P., Glymour, C. N. and Scheines, R. [2000], Causation, prediction, and search, MIT press.
[43] Tsamardinos, I., Aliferis, C. Statnikov, A. R. and Statnikov, E. [2003], Algorithms for large scale markov blanket discovery., in 'FLAIRS', Vol. 2, pp. 376-81.
[44] Wei, D., Gao, T. and Yu, Y. [2020], 'Dags with no fears: A closer look at continuous optimization for learning bayesian networks', Advances in Neural Information Processing Systems 33, 3895- 3906.
[45] Wei, Y., Li. X., Lin, L., Zhu, D. and Li, Q. [2022], 'Causal discovery discrete data via weighted normalized wasserstein distance', JEEE Transactions on Neural Networks and Learn- ing Systems 35(4), 4911-4923.
[46] Wenjuan, W., Lu, F. and Chunchen, L. [2018], Mixed causal structure discovery with application to prescriptive pricing, in PProceedings of the 27th International Joint Conference on Artificial Intelligence', pp. 51 5126-5134.
[47] Ye, Q.. Amini, networks', IEEE Transactions on Pattern Analysis and Machine Intelligence.
and Zhou, Q. [2024], 'Federated learning of generalized linear causal
[48] Yu, Y., Chen, J., Gao, T. and Yu, M. [2019], Dag-gnn: Dag structure learning with graph neural networks, in 'International conference on machine learning', PMLR, pp. 7154-7163.
[49] Zeng, Y., Shimizu, S., Matsui, H. and Sun, F. [2022], Causal discovery for linear mixed data, in "Conference on Causal Learning and Reasoning', PMLR, pp. 994-1009.
[50] Zhang, Z., Ng, I., Gong, D., Liu, Y., Abbasnejad, E., Gong, M., Zhang, K. and Shi, Q. [2022], 'Truncated matrix power iteration for differentiable dag learning Advances in Neural Information Processing Systems 35, 18390-18402.
[51] Zhang, Z., Ng. I., Gong, D., Liu, Y., Gong, M., Huang. B., Zhang, K., Hengel, v. and Shi, J. Q. [2025], 'Analytic dag constraints for differentiable dag learning". arXiv preprint arXiv:2503.19218.
[52] Zhao, B., Wang, W., Zhu, D., Liu, Z., Wang, D., Zhang, Z., Zhou, and Kolar, M. [2024], 'Personalized binomial dags learning with network structured covariates', arXiv preprint arXiv:2406.06829.
[53] Zhao, B., Zhai, P. S., Wang. Y. and Kolar, M. [2024], High-ddmensionnl functional graphical model structure learning via neighborhood selection approach', Electronic Journal of Statisties 18(1), 1042-1129.
[54] Zheng, X., Aragam, B., Ravikumar, P. K. and Xing, E. P. [2018], 'Dags with no tears: Continu- ous optimization for structure learning'. Advances neural information processing systems 31.
[55] Zheng, X., Dan, C., Aragam, B., Ravikumar, P. and Xing, E. [2020], Learning sparse non- parametric dags, in 'International Conference on Artificial Intelligence and Statisties", Pmlr, pp. 3414-3425.
[56] Zhu, S., Ng, I. and Chen, Z. [2019], 'Causal discovery with reinforcement learning', arXiv preprint 1xXiv:1906.04477
12

---

SUPPLEMENTARY MATERIAL
Differentiable Structure Learning for General Binary Data
A Preliminary Technical Results
In this appendix, we present several key technical results that are essential to our proofs. Lemma 1. Suppose X ~ MultiBernoull(p)) with > 0. Then for every subset S S [p], the marginal vector Xg satisfies
Xs ~ MultiBernoollipps)
(16)
for certain Ps E R where > 0.
The marginal vector Xs remains multivariate Bernoulli with natural parameter Ps. and positivity of p indicates the positivity of ps-
Lemma 2. Suppose X MultiBernoull((p) in the natural parameterization, or equivalently X MultiBernoulli(() f in the general parameterization. Then
p
>
If 3c.
(17)
The general parameter vector p is strictly positive if and only if its corresponding natural parameter f is strictly finite.
Lemma 3. Let p > and suppose X MultiBernoulli(p). Fix any topological order and index j E [p]. Define the population negative log-likelihood
((w) lo (1 + exp(ui (w
(x(-)))) Xx(2) w (x...)).........
(18)
where w E 2 Then E(w) is strictly convex and therefore admits a unique minimizer wa."
This optimization problem is equivalent to fitting, in the population limit, a logistic regression of Xm(i) on ((X=1))..... ())) Because f=.j. as defined in Section 4.2, is one solution to the optimization above, and we show optimal solution is unique, then w = f n. Consequently, logistic regression perfectly recovers the true natural parameter.
Corollary 2. Let f Rd + R be differentiable and m-strongly corvex, i.e.
f(v) f(x) + Vf(x) (y x)) + m |y x:12 Denote by x² the unique minimizer f. so f(x") -ii,in f(x). Then for any constant e 2 f(x*) the subset
Vx,y E TRd
(19)
2
Le = fx R. f(x) cc
(20)
is bounded.
Corollary asserts that any strongly convex function possesses bounded level sets, which is used in the proof of Theorem 3.
B Proofs
In this appendix, we provide detailed proofs of the main results.
B.1 Proof of Lemma
From the proof of Corollary we know that Xs ~ MultiBernoulli(ss) Moreover, because
0, then
P(Xs Iz))
Σ
P(Xs Zs, X -bp](s)
(21)
#bbhsE 1-11
Therefore, for any æs. P(Xs IZ) > 0. So, ps 0.

---

B.2
Proof of Lemma 2
Sufficiency By (67), each natural-paramenn probability satisfies exp(phrammsh)
(22)
p # zeros among J1, j2. Jr components and other components are all zero)
p(odd # zeros among f1, (2.J2..... .2....... Jr components and other components are all zero)
(23)
and since p 0, every term exp(fti-.r) ) is strictly positive and finite. Hence < exp (xp(f))m5r) 0c
| co.
(24)
Necessity
Note that
8irb2... Σ
f Σ
fūnh +...+ prrden
(25)
ISeSh
ISeciSh
If [f] co, then [ < 00,
0 < exp(Stiit...r)) < 00
(26)
The joint probability of observing ones in positions j1, ,,.....3r I is
P2..J....... Ir positions are one, others are zero) exp(Stntem))))
exp(b(f))
(27)
exp Sxp((Siitarr..))
E 1+ (Essnsse exp(Shiaandt))
K
which lies in the interval (0, 1):
0 < f2............. in positions are one, others are zero) < 1
(28)
B.3 Proof of Lemma 3
Fix topological order and index j E [p). To simplify notation, set
Y Xn(i {0,1), Z + (Xx,-())..,, E 10.
(29)
and let w E 2=2 be the parameter vector. The population objective is
((w) = E log(1 l + exp(w" T 2)) - Y (w T 2)
(30)
straightforward calculation shows
V2((w) = Elo(wT (0 o(w (1 2)) 2 _
o(c) = 1
(31)
Levery inner product w Z is finite and thus oww Z)(1 oww ((TT22) 2)) Since Z E 1)
0.
Define
m(w) min aww 2)(1 ow 2)) 0. =E(0.1)25
(32)
The marginal distribution of (X (Xx(1),... X ()) is still a multivariate Bernoulli distribution, with
natural parameters Pz.-- i.e., (X (Xx(I).... X......,,,-.. x))) ~ MultiBernoulli (Pr.-1)). By lemma 1,
Jmin..j-1 min 2/ [( li: O. Hence for any nonzero E 2
Σ
p(2)o(w] 2) - o(w? 2)) 2T
=E(0,1)2/-1
(33)
4min.--1
Σ
(8T2)² 2)² > 0.
=
14

---

because the sum is over all E {0.1) So (8T2)2 > O for at least one when 3 * 0.
Since the Hessian is everywhere positive definite, {(w) is strictly convex and therefore has unique minimizer.
B.4 Proof of Theorem 1
Let be any permutation of (1....,P)
p). Then the joint distribution of X can be written as
P(X) II
P(X+0) (=())...*--1)))
(34)
==1
In the population limit (n oo), each conditional law has the Bernoulli form
P(Xx() Xx(i)..... Xx(--1) = Nooi) (1 9)) (xxx())
(35)
where
g 9= a (S7, 2 ((X
(36)
Lemma guarantees that logistic regression uniquely recovers f x... Moreover, the structural equation model
Y=o) Bernoulli (o ( (
(37)
induces exactly the same conditionals as X. i.e.
Ya() Ya(I).-. Y=(-11) Xx()) Xat--))-
(38)
Hence
Mulcierrouull(p),
(39)
establishing the desired result.
B.5
Proof of Theorem 2
By Theorem for every G.) E (min the vector f, defines, via the structural equation model (7), distribution
X ~ MultiBernoulli())
(40)
that is Markov with respect to G.. By definition of Emin(p) each G3 is a sparsest graph in the equivalence class Ep)). Since all sparsest Markov representations lie in the same Markov equivalence class, it follows that for any two pairs
( Gr.), (f GTa) Emin(p)
(41)
their Markov classes coincide:
M(G = MMG
(42)
B.6 Proof of Theorem 3
The proof relies on Theorems and Deng et al. [15]. We verify Assumptions A and B from that work.
Assumption A (1) This requires the equivalence class to be finite. Since
Emin (p) < p!
(43)
the condition holds.
15

---

Assumption A (2) that in (10),
This requires the weighted adjacency matrix W (H) to be L-Lipschitz. Recall
(H)
Σ
(h
(44)
For simplicity we instead use the equivalent form
[ (W))
Σ
|AS
(45)
scipl. ES
The reason of equivalence is that they encode the same nonzero pattern of [W(H]].y Let H and H2 be two parameter values. We show there exists L such that
[[WH)) - W (H2))12 S L || H2ll2.
(46)
First,
[[ (H2))22
ΣΣ iCS
Σ
Σ
I
SJ
(47)
Scio iES
Meanwhile,
H Hell2
(h SJ 1
(48)
J sc
By Cauchy-Schree wccch-,
2
2
Σ
Σ
hg
- (1g (1550) ))
SCipl,iES
SciP),ES
scIp).iES
3|5 S |l
Σ
[
S.
11
1)
2
8cbliiS
(49)
52p-1
11 - 11
sEip).is
<2P-1
(1S 8.. (..5 )
SCIp
Hence
- (H2))
Σ
S.
h S.j
ScIp)iis
ScIplies
2p-1 1
Σ
h1
8.j
h
2
8C1P]
(50)
hg 1,8,5 2
$
2 p
Scipl
-/22(0111/2
Σ Σ
(ngs - 1990 2
scipl
-V (0 1)/2
H Hzllz
Thus one may take = 122(2-1)/2
16

---

Assumption B This requires XX) to have bounded level sets. From Lemma 3 we know that under p > 0, each population logistic loss is strongly convex. Since E[s(H: X)] is sum of p such strongly convex terms, it is itself strongly convex. By Corollary 2, any strongly convex function has bounded sublevel sets.
B.7 Proof of Theorem 4
Let x² be the topological ordering of G guaranteed by Assumption A. Then
P(X) II P(X=6) ((-))
(51)
/=1
Write Sa.j = (- (1), (" 11). Under the linear SEM assumption,
X Xs..; = US,
(ern Xs... + cr-(1))
(52)
where ()) and ()) are finite. Hence for any IS..''
< o ( 8+, TSns + Cr) <1,
(53)
so both
P( = Xs., ISra) and P(X+-() = 0 Xs,,, =SS.))
(54)
lie in (0, 1), implying p
Although every (f G=) E {p) generates X ~ MultiBernoullip)) via (7), the linear optimization (99) only admits first-order models. Thus any xf Gx) with higher-order terms misspecified. Define
Elinear (p) = ((fx Gx): (f Gx) where fa only has first order term, Vr}
(55)
By Assumption A, (Wo, GO) E glinear (p). Let
Emiu
clinnar (p) = [(f G=): : (fx: Gz) the minimal element, (Sx: Gx) glinear (p)
(56)
If (W', GO) is not already minimal, simply replace it by the sparsest element in glinear (p). so clinear
that (Wo, GO) E &min (p).
In the proof of Theorem [15]. Therefore, by Theorem 4 of Deng et al. [15], we have (W", O
verified that model meets Assumptions A and B of Deng al.
Elinear (p). and hence
linear xc.
(57)
B.8 Proof of Corollary 1
Since X MulttiBernoullip) and (p))..,(p)). Then, VS C [p]. the range of Xs
{0, 1si The corresponding density mass
P(Xs = =xa)=
Σ
P(Xs Is, X xlpss)
(58)
(0.)sE
In such way, could enumerate all the value for the v. {0, 1) and put them together to get the natural parameter for Xs. i.e., Ps- As a consequence, Xs ~ M ulti Bernouli(ps))
For any E [p], and any S E [p]Vi. V. since X, {0, so the conditional distribution P(X, Xxs) is Bernoulli distribution, and the probability
P(X, = Xs x8)
P(x, = 2j X = xz)
P(Xs = xs)
(59)
= = rs)
P(X, = 1 + P(X, = 0, X3 = F3)
17

---

B.9 Proof of Corollary 2
Since f is strongly convex,
f(v) + (y- - 2) +
m 1l x1
Vx,y E Rd
(60)
2
Take x = x² then
f(0) 2 f(x')((==)))))))>)))) + m 1ly - x²1
Vv,yER
(61)
2
Rearranging,
2(c (2 - f(x*))
m bly - x²12 se- f(x'{(('))))))))))) x-1x"1 s
f(y) c 2
(62)
{ Thus, Le S { y 11y -= v ((ce_/(=") [(x'))
a bounded ball.
m
C Supplementary Technical Details and Examples
In this appendix, we collect additionsl technical derivations, algorithms,, and definitions that support our main theorems.
In C.I, we give the explicit parameter transformation between the general parameter p and the natural parameter f of the multivariate Bernoulli model [12];
In C.2, we derive the logistic form of the conditional distributions in equation (3)
In C.3, we formalize the graded-lexicograppe ordering used to index interaction features.
In C.4, present the population-level recovery algorithms for each topological order
In C.5, we define the structural equation model framework underlying Theorem 1.
In C.6, we review faithfulness, Markov equivalence, and the Sparsest Markov Representation necessary for Theorems 2 and
In C.7, we provide the derivation of our score function X) (negative log-likelihood function).
In C.8, we provide theoretical justification for the previous work [13, 5] with general framework.
C.1 Parameter Transformation the Multivariate Bernoulli Distribution
All the material in this subsection can be found in [12]. We include here for the completeness.
The multivariate Bernoulli distribution has two different parameterization, one is using general parameter p and another one is using natural parameter f.
The density is expressed by general parameter p.
P(Y = y1, Y2 = 1/2.. YK = Уk) = p(11. ......., .9K)
p(0,0. 0, 0) ()-15-(((-15)
p(1 ......... 0)(w (5
(63)
(0, 0) (1- (1,s))
p(1, 1, 1.... 17
The density is expressed by natural parameter f.
P(Y 11, Y2 12. Yx =yK) exp 1
Σ
pirda Birdandr (v)
(64)
18

---

To simplify the notation, we could define the quantity S to be
ShiJand
Σ
p
Σ
fish + +...
(65)
1585
Isectsr
and also define the interaction function B
Biidzwrr (4)
M/ 1/2 Mj,
(66)
The following lemma shows the one-to-one mapping between general parameter p and natural parameter f.
Lemma (Parameter transformation). For the multivariate Bernoulli model, the general parameters and natural parameters have the following relationship.
(67)
p # zeros among J1. components and other components are all zero)
ptodd # zeros among j1. s. components and other components are all zero)
(68)
where # refers to the number of zeros among the superscript y;, 95 of f. In addition,
exp(Shher)))) p2i,, Ja...... ....... jrr positions one, others are zero)
(69)
p(O, 0....
and conversely the general parameters be represented by the natural parameters
..... positions are one, others are zero) exp(Siibami.))
(70)
exp(b(f))
where
exp(Sh ]
K
b() = log Σ
1
(71)
C.2 Conditional distribution of multivariate Bernoulli distribution
In this part, we derive the conditional distribution of multivariate Bernoulli distribution. Especially,
P(X, = 1 X1 = E1...., Xp-1 = 1p--1
Σ
fir 19I ... Tp
=0
(72)
r=1
1Shn <3 -.....---ppp
=a (f + fl + f fp- fp-1 + 122 x1 x2 +... + fl.P x xp-1)
f12
is known the multivariate Bernoulli distribution in exponential form can be written as
pirtand Birbant (x) (2)) )) (73)
P(X = 21
X, = xp) = exp 8
Σ
r=1
Then,
P(Xp - x-p)
=
P(Xp = X x
P(X x-p))
P(Xp = 1,
= =x-p)
P(Xp = 1, + P(X, = 0, x-p) (74)
1+ -1
19

---

where X = Xp-1).
P(X1
21,..., Xp =Tp)
phita Birtande(1))) )
8 Σ
Σ
cxp
15h 199
exx
f +
Σ
firde- Biide..Je (x) +
Σ
r=1
15h (5n
15h <SpE PE (11m-te)
PE pE(
(75)
Then,
P(X1 = 21..... [x,...,,p=111 11
exp
1
Σ
firta- Birkendr ((I-p)) 1)) +
E
finda- Birlande (10)
==1
15ji St pé (E(
ISji Sje Pt (E(eennt))
(76)
P(X1 = F1, C......,,--111.
=p = exp
f +
Σ
fūrde- Birdente
(77)
1Si Cii <i.S0 pE(Jimmee)
Finally, put them together,
P(Xp X- x-p)
1 - x=
X1Xp=1.X--=--p)
(78)
1 exp -E
Eishepee.. <Srsp firdande Borde.. pe(iimntt
((x f 1))
=0 (fp + f + f2 fp- fp-1 f12 +... + PI. ... Ipp1)
where (x) I+exp(-r)
C.3
Graded-lexicographh order
To index the 2" interaction features and corresponding parameters in a consistent way, we use the (the number of elements in S).
order subsets of [p). where for any finite set S, |S| denotes its cardinality
Definition 3 (Graded-lexicographic order). Let w be a permulation of [p]. For any subsets S,T C [p], we say S grlex T if either
1. S (T).
2. J (T and, when listing the elements of S and T in ascending order under æ,, the first index at which they differ belongs to S.
Under this rule, all subsets are grouped by increasing cardinality, and ties are broken by the usual lex order induced by x.
Example. Take p = 3 and the identity order (1, Then the graded-lexicographic sequence of subsets is
sequene
0. (1), (2), (3), (1,2), (1,3), {2,3). (1,2,3).
(79)
20

---

Algorithm 1: RECOVERPARENTSCO. T..)
PA(a(j)) O.. and X1,j f (Xx(1)... X x()))
2: Compute the general parameters P of X=, // Compute the marginal distribution of X
3: Convert Pa.j to natural parameters f. x.j utilizing Lemma 4
4: for 1.....j- do
(Fa (Ia0))0)))) > 0 then
if Σ (L) att (.NN)))))))))) PA(T()) & PA(T()
5: return (PA(T()).f=>))
Algorithm 2: DECOVERDAG(p..)
Input: Probability vector p (or empirical count) and topological sort
Output: DAG G.. and natural parameters fx
G -mmpy graph and f + 0
2 for j= 1,2
do
(PA(T()), J nn)
RECOVERPARENTSS F. j. PA(+(j)): the parents of node TT()
for PA(T(i) do Add edge X; + Xj to G fm fz (frr)
Accordingly, the extended feature map (XX) (X))ss) becomes ((X) = X1, X2, X3. X,X2, X1 X2X3, X, X1X2X3]
(80)
Similarly, for any j and order the parameter block .. R2 is arranged that its entries align one-to-one with &
in graded-lexicorraphi order.
C.4 Procedure for recovering causal graph and parameters
To formalize the recovery procedure from Section 4.2, we present Algorithms 1 and 2.
Algorithm 1: Parent and parameter recovery. For fixed topological order и and node index j E [p], this algorithm
1. computes the marginal probabilities Px.j of (Xx(I)).... xx))).
2. converts Pxjj to the natural- parameeer block f via Lemma 4, and
3. selects the parent set (z) using the nonzero-coefficien criterion in (4).
For estimating the natural-parameter block Jn. Section 4.2 uses a logistic regression approach [24]. In Algorithm I, we instead compute fx., by applying the mapping of Lemma to the marginal probabilities. Under the positivity assumption p 0, these two methods are equivalent and yield identical estimates for f..j.
Algorithm 2: Equivalence-class enumeration. This algorithm iterates over all p! permutations G, and, for each, calls Algorithm for every j l....,p. p. It assembles the corresponding DAG G₂ and parameter collection fr [f The output is the full equivalence class
E(p) (f G+) (f.. Gz) is returned for some n
Although we state the algorithm in terms of the population vector in practice the data matrix X E IR' xp, since empirical frequencies convert X into p.
simply input
Building on Algorithms and 2, Algorithm 3 then enumerates all graph-parameter pairs in E(p) and retains only those with the fewest edges, thereby recovering the minimal equivalence class Emin (p) as defined in (8).
21

---

Algorithm 3: RECOVVESSP.RSISSDACCG))
Input: Probability vector p (empirical count)
S 5+
for each π B do
3
(Gr, x) RRROOERRGRRPP(H..... T// Gy: set of all the permutation on p variables
SU S ((),)))
return ((G+,fz)) E EG: sG., V(G fx) S)
C.5 Structural equation model
An structural equation model (SEM) [35] (X, f, P(N)) over the random vector X (X1, Xp)
a collection structural equations the form:
= fx(,NN) Oxf, = Oif PA(J).
(81)
where f - (f) is a collection of functions fz= 1p+1 - I. here N = (Ni....,N) Np) of independent noises with distribution P(N), and PA(j) denotes the set of parents of node j. Here, OcLr f, denotes the partial derivative of [, w.r.f. Xx. which identically zero when fJ is independent of X4. i.e. (X, be a DAG, will be represented by the following p x weighted adjacency matrix B:
is a vector
(X, N2) (X(AA)) Nz). The graphical structure induced by the SEM, assumed to
B B
B 110,f,
(82)
and GBB) denote the corresponding binary adjacency matrix.
The structural equuttoo framework in (81) provides a fully generative foundation for causal discovery by specifying, for each variable, explicit functional dependence its direct causes. Its generality encompasses a wide range of models-lineat SEMs, additive-noise models, post-nonlinear models, generalized linear models, and more expressive non-linear SEMs. A major virtue of SEMs is their universality: any joint distribution can be represented in the form (81) ([35, Prop. 7.1]). practice, however, rescarchers impose strong parametric restrictions such linearity, additivity, or low-order interactions- that may fail to capture the full complexity of real-world data.
C.6 Faithfulness, Sparsest Markov representation, Markov equivalence class
We formally define the concepts mentioned in Section 5,
Definition 4 (Faithfulness [35]). A pair (G, P) is said to be faithful if:
X; X3 Xx =
X; and Xj are d-separated by Xx in G,
(83)
for all disjoint subsets {i.j}, V. That is, every conditional independence in P corresponds exactly d-separation in G. and vice versa.
Definition 5 (Markov Equivalence Class [42]). Two DAGs G1 and G2 on the same vertex set V are Markov equivalent if they encode the same set of conditional independence relations equivalently, they have the same skeleton (undirected edges) and the same set of v-structures (induced subgraphs of the form i > j- j with and j not aujjacent). The Markov equivalence class of a DAG
M(G) (G' G a DAG and G' is Markov equivalent to G}.
(84)
Definition 6 (Sparsest Markov Representation [38], A pair (GO,P) satisfies the Sparsest Markov Representation (SMR) assumption if:
(Go, P) satisfies the Markov property, i.e. every d-separation in G implies the correspond- ing conditional independence in P.
2. For any other DAG G &M(G)) satisfying the Markov property with respect to P, we have
|E(G) > [[
Equivalently, GO is the (unique up Markov equivalence) sparsest DAG compatible with P.
22

---

C.7
Derivation of the Logistic Loss
First, consider a single example (X,y) with feature vector e R binary label y E {0, 1}, and parameter vector w E R.. Let
1
q o=lww X)
(85)
1 1 exp( w x)
The (negative) log-likelihood is
f(w; X, y) = log 9) 1)
(86)
y log y) log(1 q)-
Noting that
log 2 w X, log(1 9) -loog(1 + exp(w X)),
(87)
q
we obtain the familiar logistic-loss form:
[(w; y) log(1 + (xp(wTx]) y (wTTX X).
(88)
In setting, each "feature" is replaced by the extended feature matrix (XX) R and each "label" is one column of the data matrix X E Raxp, Stacking over all p columns and averaging over samples yields
(H;
log (1 + exp (X))H)) X,J [*()H]} R
j=1 1=1
(89)
(log(1, + exp(((X)H)) - -X, (6(X))).
where denotes the Hadamard product and 1,, E R is the all-ones vector.
C.8 Theoretical justification for previous works
Under Assumption A, there exists topological ordering ग consistent with G such that each
Xx() is generated by a linear combination of (Xx(I)..... Xx(1-1)) E IR--1 via the logistic link,
(90)
rather than via the logistic link (((xa()...... (((())))....))))))) E R2--1 Importantly, Algorithms and 2 remain valid under this assumption. For any other topological sort ă * the output (G*, from Algorithm 2 may include higher-order interaction terms; nevertheless, by the structural equation model (7), it still recovers the exact distribution
X MultiBernoulli(p),
(91)
so Theorem I continues to hold. Moreover, under Assumption A reduce the dimensionality of the optimization (13) from H E R22 xp to H E ((++1)xp,, We formalize this reduction below.
Define the parameter matrix
H,
hP. T
RP+1
H (H I Hp) p(p+1)xp
(92)
constant
first order
where 1. is the intercept and h-. The induced adjacency matrix is
,..,, hP are the first-order coefficients.
[W |
(93)
Self-loops are forbidden, so we impose
his = 0
Vj E Ip]
(94)
23

---

Redefine the feature map row-wise
( [1, X1,
Xp],
(95)
so that for a data matrix X E R (X) applies + to each row.
The score (negative oog-likelihood) remains
{(H; X)
1n (log(In exp (xp(((X)HH) x, (((X)))
(96)
We continue to use the quasi-MCP penalty [15], defined by
Pt,s(t) = A [
quasi-MCP:
1 ([t) < 6) + (1d) > 6)
(97)
26
Our final score function is as below
sH; 6, X) = sHH; X) + PA.a ( (H))
(98)
We formulate this task the single continuous optimization problem
min H
s(H; X, 8,
subject
h(W(H))
(99)
hJ =0
Vj E (p
Define the global optimal solution of (99)
O ginner
(HH', G(W(H)))) H* a minimizer of (99)}
(100)
Let Obivaar denote the set of minimizers of (99) when the empirical loss s s(H; X, 6, X) is replaced by its population counterpart [H:(H: X, 6, X)] Let collect all the parameters in assumption A.
H C1
E Rp+iixp
(101)
W1
Wpp
Theorem 4. Suppose Assumption A holds, then ~ MulttiBernoullip) where p > O. Moreover, there exist X, 8 > O sufficiently small such that (H", , E Obinnar where G is the ground truth graph in Assumption A.
is important to note that under this assumption
Olincar * Emin(p))
(102)
because there may exist topological sort for which (f Gz)) E (min(p) involves higher-order terms, whereas every solution in Olincur contains only first-order terms.
D Experiments
In this section, present comprehensive experimental details, including the graph types evaluated, the data-generation process, the bascline methods for comparison, the steps required to reproduce our implementation, and the evaluation metrics employed.
D.1 Experimental Setting
In this section, we outline the process for generating graphs and data. For each model, random graph G is generated using one of two types of random graph models: Erdös-Rényi(END or Scale-Free (SF). The models are specified to have, on average, kp edges, where {1, 2, 4). These configurations are denoted as ER or SFk, respectively.
Random graphs whose edges are add independently with equal probability, We simulated models with p. 2p and 4p edges (in expectation) each, denoted by ERI, ER2, and ER4 respectively.
Scale-free netW Network simulated according to the preferential attachment process. We simulated scale-free network with p, 2p and 4y edges and / 1, where 3 the exponent used in the preferential attachment process.
24

---

Algorithm 4: Generate data matrix X
Input: DAG G. sample size n., interaction type 7 {1st+2nd, 1st+pth, 1st+2nd+pth 2nd, pth)
Output: X E {0, 1Jexp
Compute topological ordering x of G
for j < to do
Sample "PA()(()) = (w
E IR?
PA
iid wk ~ Unif([-2, -1] -I [1.
if = 1st 1st+2nd then
q ("PA(=G)) T 1sr+22d (PP(F)))))
else if r = 1st+pth then
9 - o ("TA(+(()) 1st+ psh (XPA(T())))
else if r = 2nd then
+ σ (PA(()())) 2nd (PA(#()))))
else
o ((P)))())) @uch (XP(((O))
11
12
XT() Bernoulli(g)
General binary data Since we wish to study structure learning for general binary data, Theorem implies that for any p > 0,
X MultiBerrouuli(p)
(103)
can be generated via the SEM (7). To allow different interaction orders, define the following extended feature maps for X- (X1...., Xp):
2st (X)
1
X1,
Xps x2x.....,
0_
) 1
constant
first order
second order
third order forth 1) order p-th onder
T Ist+2nd+pphh (X)
X X1X2, x X X,
R
constant
first order
second order
third order fourth (p- 1) th order
p-th order
Φ Ist (X)
X1 Xp
E R
constant
first order
second order third order
forth to (p-1)-th order
p-th order
2nd (x)
X2X Xp-
T E R
constant
first order
second order
third order forth (p- 1)-th order
p-th onder
BPth (X)
X1 X1 Xp)
E R22
constant first order
second order third order forth to (p- 1) order
p-th order
(104)
where each vector the nonzero blocks correspond respectively to the constant term, first-order terms, second-order terms, and highest-order term. By convention, if X = 0, then all four maps reduce the scalar (1) When applied to data matrix X E Raxp, each ((X) operates row-wise.
Finally, given random DAG B E {0, 1 sampled from one of our graph models, we generate X using Algorithm 4, choosing the desired interaction map according to whether we study first+second, first+highest, second, or highest-order interactions.
Simulation We generate random dataset X E 1xp by sampling i.i.d from the models described above. For each simulation, we produce datasets with n samples cross graphs with p nodes.
(Small Graph) p (5, 9) & (1. 10000 and graph types: (ER,SF)
(Large Graph) p 210,20,30 30, 40), A (1, 4). 1000 and graph types: (ER,SF)
25

---

D.2 Implementation
For each dataset, we applied several structural learning algorithms, including fast greedy equivalence search (FGES [37]), constraint-based methods (PC [42]), DAGMA [5], NOTEARS [55]. The implementation details are provided in the following paragraph. After running the algorithms, a post-processing threshold of 0.3 was applied to the estimated Best to prune small values, following the same procedure in [54, 5].
Fast Greedy Equivalence Search (FGES [37]) is based on greedy search and as- sumes linear dependency between variables. The implementation is based on the py-tetrad package, available at attps://gg //tg/gtttub... structuree_riio==)).
hatptupphhpppyttpeepe We use
search. use_ (sample...iorreeeeeeeee.eeeeeeee
PC [42] is constraint-based method and based on uses conditional independence induced by causal relationships to learn those causal relationships. The implementation based on the py-tetrad package, available at https: / / github..c.. com/ ha//cmmu-ppa comu-philppy-ttteeee We use search. use_chi teee_chi_suuurreeeeeeeeeeee square (alpha=0.)
NOTEARS-MLP [55] is continuous DAG-learning method that employs a least- squares loss with 4 regularization. Its Python implementation is available at https://githute com/xunzheng wwunzheng/ootears...
notears.
DAGMA [5] is a continuous DAG-leaming algorithm that achieves improved accuracy and faster computation, with barrier methods. Its implementation can be found at https://githut. cwm//eevinselll//aggeeoooo We use the default parameters from origi- nal implementation.
BINOTEARS builds on DAGMA and NOTEARS-MLP with the choice of solver tailored to problem size. For small graphs (d E {5, 6, 7, 7,,,9))), we optimize the full formulation (13) (13)-including all higher-order interaction terms using the DAGMA framework. For larger graphs (d E (10, 20, 30, 40)), we adopt the NOTEARS-MLP framework. The original implementation optimizes an & loss with an & penalty; in our variant, we insert a sigmoid activation a (z) = 1/(1 + exp(-x), on the final layer and replace the squared loss with the cross-entropy (logistic) loss to accommodate binary data. After estimating the weighted adjacency Best via NOTEARS-MLE. ordering of the resulting graph, and then apply Algorithm 2 with first- and second-orde -odder terms obtain the final structure. Finally, remove any remaining edges whose weight does not exceed 1.0 to eliminate spurious connections.
prune entries below 0.3, compute a topological
Hyperparameter tuning Theorem 3 indicates that one should ideally choose small values of À and for the quasi-MCP penalty. In practice, however, achieving the global optimum of (13) is infeasible, and if X and 8 are too small the penalty becomes ineffective and the algorithm may fail to recover the sparsest solution. To mitigate this, we adopt the continuation strategy of Deng et al. [15]: start with relatively large À and ò, solve (13) via BINOTEARS to obtain an initial estimate Hest, then iteratively shrink and by a factor y < 1, using the previous estimate as the warm start for the next run of BINOTEARS. We terminate when the negative log-likelihood s sHest, X) ceases to decrease.
Empirically, y = 0.5, X = 0.05, and = 0.2 perform well in our experiments.
Equipment
The experiments are conducted in the following CPU architectures
Intel Broadwell-28 cores @ 2.4 GHz with 64 GB memory per node
Intel Skylake cores @ 2.4 GHz with 96 GB memory per node
D.3
Metrics
Structural Hamming distance (SHD): A standard benchmark in the structure learning literature that counts the total number of edges additions, deletions, and reversals needed to convert the estimated graph into the true graph. Since our data specified in (1) is nonidentifiable, the Structural Hamming Distance (SHD) calculated with respect to the completed partially directed acyclic graph (CPDAG) of the ground truth and B
26

---

E
Additional Figures
E.1
Small graphs
15-
10.
Soutural Hamming Disbance (SHO)
Stautural Hamming Distance (5HD)
2-
0.0
7.5
7.5
5.0
0.0
#Numbor modes)
4 d'Nunborof recceal
ONumbbrr of modes)
Mitedd on DAGMA
DINOTEANS (Dan)
rocs sne PC
Miticd mns
DHOMA
DINOTEARS rOux: FGES s
(a) Extended feature map o Istr2ndeeth (104)
(b) Extended feature map 2nd in (104)
Figure 3: Results in terms of SHD between MECS of estimated graph and ground truth. Lower is better. Column: & = f1, Row: random graph types. {ER,SF)-k = [Scale-Free Rényi graphs with kd expected edges. Here p {5, 6, 7, 9}. Error bars denote the standard error computed over 10 replications.
E.2 Large graphs
60
125
30-
100
Structural Hamminng Distance (SHD)
20-
40
EN
75
10.
20
50
40-
60
30
100
20
10
20
10
20
30
40
10
20
30
40
10
20
30
40
(Number of modes)
d dWummer of nodes)
d (N ummber of nodes)
Method
DAGMA
FGES BINOTEARS (Ours)
PC
Figure 4: Results in terms of SHD between MECS of estimated graph and ground truth. Lower is better. Data are generated using extended feature map Isispih in (104). Column: k = {1, 2, 4). Row: random graph types. (ER,SF)-& (Scale-Free,, Rédōs-Rényi] graphs with kd expected edges. Here {10,20 20, 30, 40J. BINOTEARS is our two stage approach.
27

---

80-
50-
40
100
Structural Hamming Distance (SHD)
30
ER
40
20
50
10-
20
80
50-
40
100
30
15
45
40
20
50
10
20
10
20
30
40
10
20
30
40
30 (Number of nodes)
20
(Number of nodes)
(Number of nodes)
Method
FGES
PC
DAGMA
BINOTEARS (Ours)
Figure 5: Results in terms of SHD between MECs of estimated graph and ground truth. Lower is better. Data are generated using extended feature map 2 in (104). Column: k {1, 4}. Row: random graph types. (ER,SF)-k = (Scale-Free, Erdōs Erdōs-Rényi] graphs with kd expected edges. Here
p {10, 20, 30, 40). BINOTEARS is our two stage approach.
160
40
120
30
Structural Harming Distance (SD)
20
80
10-
40
60
160
40-
120
15
80
20
25
40
10
20
30
40
20
30
40
30 (Number of nodes)
20
(Number of nodes)
(Number nodes)
Method
DAGMA
FGES
BINOTEARS (Ours)
PC
Figure 6: Results in terms of SHD between MBCs of estimated graph and ground truth. Lower is better. Data are generated using extended feature map 4pth in (104). Column: & = {1, 4J. Row: random graph types. (ER,SF)-K = (Scale-Free, Erdös rrdōs-Rééni] graphs with kd expected edges. Here p {10, 20, 30, 40). BINOTEARS is our two stage approach.
28

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Each claim in the abstract can be founded in the main paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: the conclusion section, we have briefly discussed about future work, which address the limitation of the paper.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]

---

For each theorem statement, we include the condition/aassumtiion needed, see Theorem 1,2,3.
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification:
We can provide the experiments' details in the appendix D.
Guidelines:
The answer NA means that the paper does include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
30

---

Answer: [Yes]
Justification:
The code will available if the paper get accepted.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips..
oneiionnolcyy for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: /CodeSubmic / //nips. wttpccccpppccccceeececccccccccct / ssisnnoillllo onPolicy))) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
See the experiment section in the appendix for reproducible purpose.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: Error bars are provided in the Experiment in the main paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, at least for the experiments that support main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
31

---

It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Such details are included in the appendix D.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g.. preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: /ttp///eers///e//ee//ee httpppic/tthiscuuiiiiteepappeaee
Answer: [Yes]
Justification: Our work is more related to theory and methodology.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is a special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: We mainly focus the theoretical findings.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to
32

---

generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
there are negative socictal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
It is not relevant to research topic.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: is unclear to me.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's is creators.
13. New assets
Question: Are new assets introduced in the paper well documented and the documentation provided alongside the assets?
33

---

Answer:
[NA]
Justification: Unclear to me.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: Not relevant.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: No human involved.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Answer: [No]
34

---

Justification:
Our work is original.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy for what should or should not be described.
: http///rrpss.ccoooofeeeeeeeeee/e
35


Reviews:
[{"review_id": "Vou352P4Pu", "reviewer": "Reviewer_eR7S", "summary": "This paper presents a rigorous and well-executed study on differentiable structure learning for general binary data. Its main theoretical contribution is the proof that causal discovery is non-identifiable from purely observational data in the most general binary setting. By leveraging the multivariate Bernoulli distribution, the authors fully characterize the set of graph–parameter pairs consistent with the data, and introduce the notion of a minimal equivalence class to select among these.\n\nThe second major contribution is a novel differentiable optimization framework that recovers the causal structure up to Markov equivalence, without relying on restrictive assumptions such as linearity or additive noise. This formulation enables the recovery of the sparsest DAG consistent with the data, supported by strong theoretical guarantees.", "questions": "See weakness parts.", "limitations": "My main concern is about scalability. The size of $\\mathbf{H}$ matrix increase exponentially against the number of nodes in the graph (or at lease cubic if use approximation), which prevents the  algorithm to scale to even graph with near 100 nodes. Also it would be good if the authors can provide runtime comparison of different algorithms. algoriithms", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths\n### Theoretical Rigor:\nThe paper provides a strong theoretical foundation by proving the non-identifiability of causal discovery for general binary data. It rigorously characterizes all graph–parameter pairs compatible with the data using the multivariate Bernoulli distribution.\n\n### Novel Optimization Framework:\nA differentiable formulation is introduced to recover the structure up to the Markov equivalence class, enabling the use of gradient-based methods without restrictive parametric assumptions (e.g., linear or additive noise models).\n\n### General Applicability:\nThe approach handles arbitrary higher-order dependencies in binary data, overcoming limitations of existing methods that assume first-order or linear models.\n\n### Empirical Validation:\nExperiments on synthetic datasets demonstrate the method’s superiority in capturing complex causal structures, especially when higher-order interactions are present.\n\n\n## Weaknesses\n\n### Minor Writing/Typesetting Issues:\nThe paper appears to have been written in a rush. There are several typesetting and referencing errors—for example, broken or missing references to appendices (e.g., “Appendix ??”), which slightly detract from readability and polish.\n\n\n### Limited Empirical Analysis of Optimization Aspects (Minor):\nWhile the authors adopt the DAGMA framework for optimization, the analysis of the DAG constraints themselves is underdeveloped. The paper could be strengthened by a more detailed empirical or theoretical investigation of how the chosen DAG constraint affects optimization performance and recovery quality.\n\nIn particular, recent work suggests that the acyclicity constraint in DAGMA could be replaced with:\n\n- Matrix inversion-based constraints [1,2], if vanishing gradients are a concern.\n\n- NOTEARS-style constraints, if non-convexity of the objective is the dominant challenge.\n\nIncorporating or comparing such alternatives could improve the overall robustness and insight of the paper’s optimization component.\n\n[1] Zhang, Zhen, et al. \"Truncated matrix power iteration for differentiable DAG learning.\" Advances in Neural Information Processing Systems 35 (2022): 18390-18402.\n[2] Zhang, Zhen, et al. \"Analytic DAG Constraints for Differentiable DAG Learning.\" arXiv preprint arXiv:2503.19218 (2025).", "was_revised": true, "final_justification": "I believe my concern has been resolved.", "rebuttal": "We thank the reviewer for acknowledging that our paper offers a rigorous theoretical framework for general binary data—capable of modeling arbitrary higher‑order interactions—as well as a novel optimization framework. We appreciate this recognition and welcome the opportunity to address any remaining concerns.\n\n> Minor Writing/Typesetting Issues: .... readability and polish.\n> \n\nWe apologize for this oversight and any inconvenience it may have caused during the review. The issue arose because the manuscript was compiled without the appendix content. A clean version—including all appendices—is available in the supplementary material, and we will ensure this does not happen again. Thank you for bringing it to our attention.\n\n> Limited Empirical Analysis of Optimization Aspects\n> \nThank you for highlighting these foundational works. Our method’s success hinges on how we enforce the acyclicity constraint in optimization (13), both theoretically and numerically. In the revised manuscript, we will add a dedicated discussion of prior approaches to DAG constraints—particularly the theoretical contributions of [1,2]—and report new experiments in which we replace our DAGMA constraint with the acyclicity formulation from [1,2], presenting those comparative results. \n\nBelow, we report the runtime and SHD when the exponential acyclicity constraint is replaced by the method of [1], yielding our variant NOTEARS‑MLP‑REG‑TMPI. The source code for [2] is not yet publicly available; we will perform a similar comparison once it is released on GitHub. These results demonstrate that incorporating the enhanced acyclicity constraints from [1,2] can further improve our method’s performance. Thanks for pointing out this.\n\n\n| Times(second) | NOTEARS-MLP-REG(Ours)| NOTEARS-MLP-REG-TMPI |\n|---------------|:---------------:|:--------------------:|\n| ER2 (d=40)    |     295.333     |        272.389       |\n| ER4 (d=40)    |     337.205     |        301.232       |\n| SF2 (d=40)    |     293.903     |        273.130       |\n| SF4 (d=40)    |     313.349     |        292.179       |\n\n| SHD        | NOTEARS-MLP-REG(Ours)| NOTEARS-MLP-REG-TMPI |\n|------------|:---------------:|:--------------------:|\n| ER2 (d=40) |      42.323     |        43.172        |\n| ER4 (d=40) |      73.323     |        67.293        |\n| SF2 (d=40) |       48.2      |        43.890        |\n| SF4 (d=40) |      103.2      |        92.333        |\n\n> My main concern is about scalability..... Also it would be good if the authors can provide runtime comparison of different algorithms.\n> \n\nWe appreciate the reviewer’s concern about scalability. As noted in Lines 314–322, the size of our feature map $\\Phi(\\mathbf X)$ (and thus the cost of solving (13)) grows exponentially with the number of variables $p$, so a direct solution quickly becomes infeasible for even moderately large graphs. To address this, we introduce a two‑stage heuristic called NOTEARS‑MLP‑REG. In the first stage, we train an MLP (NOTEARS‑MLP[3]) to capture higher‑order interactions in discrete data and recover an initial graph structure. From this preliminary graph we extract a node ordering $\\hat{\\pi}$, which we then use in the second stage to fit a logistic‑regression model—including both first‑ and second‑order terms—to obtain the final graph. Because logistic regression under a fixed ordering is very efficient, the MLP‑based graph recovery remains the dominant computational cost. The runtime of our method (NOTEARS‑MLP-REG) is quite similar to NOTEARS‑MLP[3]. \n\nThe following table provides the runtime for different methods and will be included in the final version. Here,  DAGMA is renamed as “DAGMA-1st”, to emphasize only linear term is used.\n\n\n|Times(second)| NOTEARS-MLP-REG(Ours) | DAGMA-1st | FGES   | PC     |\n|------------ |:---------------:|:---------:|:------:|:------:|\n| ER2 (d=10)  | 43.223          | 11.291    | 2.2565 | 3.3855 |\n| ER4 (d=10)  | 38.972          | 12.887    | 2.3755 | 5.522  |\n| SF2 (d=10)  | 37.925          | 12.789    | 2.192  | 3.3545 |\n| SF4 (d=10)  | 41.259          | 11.567    | 2.658  | 4.4015 |\n\n|Times(second)| NOTEARS-MLP-REG(Ours) | DAGMA-1st | FGES   | PC      |\n|-------------|:---------------:|:---------:|:------:|:-------:|\n| ER2 (d=20)  | 84.392          | 26.291    | 3.7165 | 6.6575  |\n| ER4 (d=20)  | 87.267          | 28.987    | 6.2005 | 14.349  |\n| SF2 (d=20)  | 80.952          | 25.781    | 3.7605 | 6.5365  |\n| SF4 (d=20)  | 93.136          | 28.659    | 5.6295 | 11.7385 |\n\n|Times(second)| NOTEARS-MLP-REG(Ours) | DAGMA-1st |  FGES  |  PC    |\n|-------------|:---------------:|:---------:|:------:|:------:|\n| ER2 (d=40)  |     295.333     |   43.232  | 4.2745 |  14.49 |\n| ER4 (d=40)  |     337.205     |  58.2323  | 7.6775 | 37.806 |\n| SF2 (d=40)  |     293.903     |   42.232  |  4.215 |  23.11 |\n| SF4 (d=40)  |     313.349     |  48.6789  | 11.616 |  21.21 |\n\nFrom the reported runtimes, it is clear that our method does not exhibit exponential time complexity and can deliver meaningful discoveries within a reasonable timeframe.\n\n__Reference__\n\n[1] Zhang, Zhen, et al. \"Truncated matrix power iteration for differentiable DAG learning.\" Advances in Neural Information Processing Systems 35 (2022) \n\n[2] Zhang, Zhen, et al. \"Analytic DAG Constraints for Differentiable DAG Learning.\" arXiv preprint arXiv:2503.19218 (2025).\n\n[3] Zheng, Xun, et al. \"Learning sparse nonparametric dags.\" *International conference on artificial intelligence and statistics*. Pmlr, 2020."}, {"review_id": "Fw2NAm4z5j", "reviewer": "Reviewer_coPu", "summary": "This paper studies the problem of causal discovery in a set of discrete binary variables without any specific model. They propose a differentiable learning method that formulates the learning problem as a single differentiable optimization task. They show that general discrete models are unidentifiable from purely observational data and characterise the Markov equivalence class.", "questions": "Please see above.", "limitations": "Please see above comments.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper is well-written. The problem is a relevant and important one. \nThe logistic formulation of the conditional multivariate Bernoulli probabilities is quite interesting. This helps the authors to form the differentiable form and might be also useful for further researches. \n\nThere are couple of major weaknesses with this work. The contribution is limited given that the previous results on differentiable method for causal discovery such as notears. \n\nThe proposed method also suffers from the limitations of the notears like methods such as getting trapped into local minima, speed of convergence, good initialization, and the choice of parameters such as regularisation coefficient. There are no discussion about these points in the paper. \n\nNotation B^{j_1,…,j_r} was defined for simpler representation but in some place it has not been used while it should eg. Eq. (3).\n\nReferences are missing to the Appendix and in some places equations are out of page border eg. (1).\n\nBefore remark 2, “The cardinality of E(p) is at most p!”. Why not exactly p! ?\nThe summation in score function definition should be to n instead of p.\nFunction (11) is quadratic on [-delta, delta] not only on [0, delta].", "was_revised": true, "final_justification": "I increase my score as the authors' rebuttal could address some of my concerns.", "rebuttal": "We thank the reviewer for recognizing the clarity of our presentation, for finding our work both interesting and important, and for providing many constructive comments that will help us improve the manuscript. We appreciate the opportunity to address any points of confusion and clarify our contributions.\n\n> The contribution is limited given that the previous results on differentiable method for causal discovery such as notears.\n> \n\nThank you for raising this concern. While NOTEARS introduced an elegant differentiable acyclicity constraint that makes gradient‑based causal discovery possible, it offers **no** theoretical guarantees on its solutions and leaves unanswered what is actually learned when and how applied to discrete data. Our work fills this gap by developing a unified framework for general binary data that delivers clear theoretical properties for the estimated graph by adapting the NOTEARS formulation (see Theorem 3). Specifically, Theorem 3 shows that any optimal solution of Equation (13) recovers parameters that **exactly** reproduce the observed data distribution and yields a DAG whose conditional independencies coincide with those of the data. By modeling the **general** distribution of binary variables, we avoid the restrictive linearity assumptions of earlier methods and can capture arbitrary interactions in real‑world data. To address the exponential blow‑up in complexity inherent to general binary models, we further introduce a scalable two‑stage heuristic (NOTEARS‑MLP‑REG) that retains full modeling flexibility while remaining computationally practical on large datasets. Together, these contributions go well beyond a straightforward application of existing differentiable approaches, offering new theoretical insights, enhanced modeling power, and genuine scalability. We respectfully request the reviewer to reconsider our contributions in light of the explanations provided above.\n\n> The proposed method also suffers from the limitations of the notears like methods such as getting trapped into local minima, speed of convergence, good initialization, and the choice of parameters such as regularization coefficient. There are no discussion about these points in the paper.\n> \n\nThank you for this helpful suggestion! Practical implementation details are indeed critical to our method. In the main text (Line 96) we briefly review related optimization‐focused studies, and in Appendix D.1 we describe our hyperparameter‐tuning procedure for the regularization coefficients. Several works have rigorously examined NOTEARS‐style methods from these implementation perspectives: [1] characterizes the KKT and local‐optimality conditions under a linear Gaussian model; [2] extends these analyses to more general settings; [3] proves that an appropriate optimization scheme converges to the global minimum in the bivariate case; and [4, 5] investigate the convergence behavior and nonconvexity of NOTEARS. Although we only mention these studies in the current submission due to space constraints, we will expand both our implementation details and our discussion of these related work in the revised manuscript.\n\n> Notation B^{j_1,…,j_r} was defined for simpler representation but in some place it has not been used while it should eg. Eq. (3).\n> \n\nThanks! We will modify Eq. (3) in the updated version and remove the usage of $B^{j_1,…,j_r}$.\n\n> References to the Appendix are missing and in some places equations are out of page border eg. (1).\n> \n\nWe apologize for these oversights and any inconvenience they may have caused. In the revised manuscript, we will split Equation (1) across multiple lines so that it fits within the page borders. The missing references were due to compiling the paper without the appendix; a complete version without missing reference is available in the supplementary material. We will ensure that these issues do not recur in future submissions.\n\n> Before remark 2, “The cardinality of E(p) is at most p!”. Why not exactly p! ?\n> \n\nWe write **“at most  $p!$”** because different node orderings can induce the same DAG and parameter mapping, so not all $p!$ permutations yield distinct models. Formally, one can have $\\pi_1\\ne \\pi_2$ with $f_{\\pi_1} = f_{\\pi_2}$ and $G_{\\pi_1} = G_{\\pi_2}$. For example, if $X = (X_1,….,X_p)$ are mutally independent any topological sort $\\pi$ produces the same empty graph and identical parameters, giving $|\\mathcal{E}(\\mathbf{p})|=1$, rather than $p!$\n\n> The summation in score function definition should be to n instead of p.\n> \n\nWe appreciate your keen attention to detail. The summation over $p$ features is indeed correct, however. In our score function $\\ell(\\mathbf{H};\\mathbf{X}),$ the index $i$ runs over the $p$ features $X_1,\\dots,X_p$, where $\\mathbf{X}_i$ denotes the $i$th column of the data matrix $\\mathbf X$. The sum over the $n$ observations is implicitly carried out by the inner product with the vector $\\mathbf1_n^\\top$ in the term that follows. We recognize that this notation can be confusing, and we will clarify it in the revised manuscript.\n\n> Function (11) is quadratic on [-delta, delta] not only on [0, delta].\n> \n\nThank you for pointing this out. In fact, the quasi‑MCP penalty $p_{\\lambda,\\delta}(t)$ is an even function, symmetric about zero. On each of the subintervals $[0,\\delta]$ and $[-\\delta,0]$, it is given by a quadratic expression with negative curvature that attains its maximum at $|t| = \\delta$. To clarify this in the revised manuscript, we will include a plot of $p_{\\lambda,\\delta}(t)$ and explicitly state that it is quadratic separately on $[-\\delta,0]$ and $[0,\\delta]$, though no single quadratic formula covers the full range $[-\\delta,\\delta]$.\n\n\n\n[1] Wei, et al \"DAGs with No Fears: A closer look at continuous optimization for learning Bayesian networks.\" *Advances in Neural Information Processing Systems* 33 (2020)\n\n[2] Deng, et al. \"Optimizing notears objectives via topological swaps.\" *International Conference on Machine Learning*. PMLR, 2023.\n\n[3] Deng, et al. \"Global optimality in bivariate gradient-based DAG learning.\" *Advances in Neural Information Processing Systems* 36 (2023)\n\n[4] Ng, et al \"Structure learning with continuous optimization: A sober look and beyond.\" *Causal Learning and Reasoning*. PMLR, 2024.\n\n[5] Ng, et al. \"On the convergence of continuous constrained optimization for structure learning.\" *International Conference on Artificial Intelligence and Statistics*. Pmlr, 2022."}, {"review_id": "1eHjlKQJUv", "reviewer": "Reviewer_QLcb", "summary": "This paper addresses differentiable structure learning for discrete data. Arguing that prior work relies on overly restrictive parametric assumptions, the authors propose a general framework starting from the multivariate Bernoulli (MVB) distribution to capture arbitrary dependencies. This paper proposes finding the sparsest graph, relying on the Sparsest Markov Representation (SMR) assumption to guarantee uniqueness up to the Markov Equivalence Class, and it translates this search into a differentiable optimization problem and introduces a two-stage heuristic to address the primary method's exponential complexity.", "questions": "I'd appreciate the authors' feedback on the weakness section.", "limitations": "Please see the weakness section.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n\n1. The paper is exceptionally well-written. The argument flows logically from the limitations of prior work to the general theoretical characterization, the proposed solution, and the empirical validation. The core ideas, definitions, and theorems are presented with precision and clarity.\n\n2. The formulation of the problem as a continuous optimization program is elegant. The connection established in Theorem 3, which guarantees that the solution to the practical optimization problem corresponds to the theoretically desired minimal equivalence class, is a crucial and impressive result that bridges theory and practice.\n\n\nWeaknesses:\n\n1. The paper's solution to the non-identifiability problem relies heavily on invoking the Sparsest Markov Representation (SMR) assumption. Theorem 2, which states that the sparsest graphs belong to a single MEC, is essentially a direct corollary of this assumption rather than a novel derived result. The paper expertly characterizes a problem in Theorem 1, but the principle for solving the resulting identifiability issue is assumed rather than discovered.\n\n2. The experiments are performed only on synthetic data. While this demonstrates that DAGMA-HO works in specific cases where higher-order interactions are known to exist, the lack of any application to real-world data makes it difficult to assess the practical relevance of the problem and the utility of the proposed heuristic.\n\n\nMinor: Format issues -- all appendix references are not rendered.", "was_revised": true, "final_justification": "After reading the authors' response and the other reviewers' comments, I've decided to maintain my current rating.", "rebuttal": "We thank the reviewer for acknowledging that our paper is well‑written and that our problem formulation is elegant. We welcome the opportunity to address and clarify the reviewer’s concerns.\n\n> The paper's solution to the non-identifiability problem relies heavily on invoking the Sparsest Markov Representation (SMR) assumption. Theorem 2, which states that the sparsest graphs belong to a single MEC, is essentially a direct corollary of this assumption rather than a novel derived result. The paper expertly characterizes a problem in Theorem 1, but the principle for solving the resulting identifiability issue is assumed rather than discovered.\n> \n\nThank you for this insightful comment. The identifiability issue arises because, with purely observational data, multiple DAGs can encode the same joint distribution—a phenomenon we formalize in Theorem 1. Without additional assumptions or interventional data, the true graph cannot be uniquely determined. While the classical faithfulness assumption addresses this, it has been shown to be overly restrictive [1]. Instead, we adopt the Sparsest Markov Representation (SMR), which is strictly weaker than faithfulness [2] and guarantees that every DAG with the minimal number of edges lies in the same Markov equivalence class. This aligns with the principle of parsimony—selecting the simplest model among those that fit the data—which is known to improve interpretability, generalization, and robustness. Notably, even without SMR, solving (13) recovers a sparsest model that exactly reproduces the observed data (Theorems 1 and 3). Furthermore, SMR is not merely theoretical: in our simulations for $p\\le8$, Algorithm 3 (which exhaustively recovers all sparsest graphs, as detailed in Appendix C.4) empirically confirms that all sparsest DAGs occupy the same MEC. We will expand this discussion in the revised manuscript to clarify the role and justification of SMR.  Thanks for bring up this point to our attentions.\n\n> The experiments are performed only on synthetic data. While this demonstrates that DAGMA-HO works in specific cases where higher-order interactions are known to exist, the lack of any application to real-world data makes it difficult to assess the practical relevance of the problem and the utility of the proposed heuristic.\n> \n\nThank you for emphasizing the need for a real‑data evaluation. We applied DAGMA‑HO to the [3] benchmark dataset, which comprises $n=7466$ continuous expression measurements of $d=11$ proteins and phospholipids in human immune cells and includes a consensus network widely accepted by the biological community. We worked with discretized dataset, and ran DAGMA‑HO. Our method recovers the consensus structure with competitive accuracy compared to existing approaches, demonstrating its potential existence of higher order interaction.\n\n\n|                | SHD | Num of Edges |\n|----------------|:---:|:------------:|\n| NOTEARS-Linear [4] |  22 |      18      |\n| NOTEARS-MLP [5]    |  16 |      13      |\n| DAGMA-HO (Ours) |  13 |      15      |\n\n\n> Minor: Format issues -- all appendix references are not rendered.\n> \n\nWe apologize for these oversights and any inconvenience they may have caused. The missing references were due to compiling the paper without the appendix; a complete version without missing reference is available in the supplementary material. We will ensure that these issues do not recur in future submissions.\n\n### Reference\n\n[1] Uhler, Caroline, et al. \"Geometry of the faithfulness assumption in causal inference.\" *The Annals of Statistics* (2013): 436-463.\n\n[2] Raskutti, Garvesh, and Caroline Uhler. \"Learning directed acyclic graph models based on sparsest permutations.\" *Stat* 7.1 (2018): e183.\n\n[3] Sachs, Karen, et al. ”Causal protein-signaling networks derived from multiparameter single-cell data.” Science 2005.\n\n[4] Zheng, Xun, et al. \"Dags with no tears: Continuous optimization for structure learning.\" *Advances in neural information processing systems* 31 (2018).\n\n[5] Zheng, Xun, et al. \"Learning sparse nonparametric dags.\" *International conference on artificial intelligence and statistics*. Pmlr, 2020."}, {"review_id": "OQ2vwPxYxw", "reviewer": "Reviewer_1Kmy", "summary": "This paper provides a general differentiable structure learning framework for causal discovery on binary data from the multivariate Bernoulli distribution, embedding all higher-order interactions without restrictive parametric assumptions. \nIt first provides a non-identifiability result under matching observable distributions. \nThen it shows the proposed algorithm identifies the sparest DAG up to Markov equivalence class under faithfulness. \nExperiments are conducted on simulated synthetic data, compared with a set of causal discovery baselines.", "questions": "1. Following up on the first weakness: it would be helpful to include a real-world motivational example where modeling higher-order interactions is crucial for causal discovery. This would better justify the need for such a general formulation beyond synthetic settings.\n2. What is the computational complexity of the proposed DAGMA-HO method, particularly with respect to the number of nodes and the order of interactions?\n3. Is it always beneficial to model higher-order interaction terms? How should one choose the order of interaction terms in practice? Is there a principled way to select this, or is it treated as a hyperparameter?", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "### Strengths\n1. The paper tackles a broad nonparametric framework that better reflects real-world scenarios, and backs its claims with comprehensive theoretical proofs.  \n2. The narrative is overall well-structured.\n3. On synthetic Erdős–Rényi and scale-free graphs, the method achieves competitive performance compared to existing baselines.  \n\n### Weaknesses\n1. While modelling higher-order interactions is motivated by real-world use cases, a motivational example is missing. Since the experiments are also focusing on synthetic data, it is not clear to me why higher-order interactions should play an important role in real-world causal discovery tasks.  \n2. As noted by the authors, the full higher-order feature map scales exponentially with the number of nodes, limiting its applicability to larger networks. The NOTEARS-MLP-REG heuristic helps in practice but lacks theoretical guarantees, making its connection to the paper’s core identifiability results less clear.  \n3. The notation in Section 4 is difficult to follow. Simplifying or redesigning symbols would greatly improve readability, e.g., by reducing super- and subscripts.\n4. Visualization: Please keep the color-coding aligned for all figures in the paper! For example, FEGS was marked as purple in figure 1 but orange in figure 2.", "was_revised": true, "final_justification": "The rebuttal has addressed most of my concerns, and I am happy to increase my score to 4. I kindly encourage the authors to ensure these clarifications are incorporated into the revised manuscript to improve overall clarity.", "rebuttal": "We thank the reviewer for their constructive comments and are pleased that they find our paper well‑structured and theoretically comprehensive. We will address all concerns and incorporate their suggestions in the revised manuscript.\n\n> While modelling higher-order interactions is motivated by real-world use cases, a motivational example is missing. Since the experiments are also focusing on synthetic data, it is not clear to me why higher-order interactions should play an important role in real-world causal discovery tasks.\n>\n\n> Following up on the first weakness: it would be helpful to include a real-world motivational example where modeling higher-order interactions is crucial for causal discovery. This would better justify the need for such a general formulation beyond synthetic settings.\n>\n\nThank you for emphasizing the need for a real‑data evaluation. We applied DAGMA‑HO to the [1] benchmark dataset, which comprises $n=7466$ continuous expression measurements of $d=11$ proteins and phospholipids in human immune cells and includes a consensus network widely accepted by the biological community. We worked with discretized dataset, and ran DAGMA‑HO. Our method recovers the consensus structure with competitive accuracy compared to existing approaches, demonstrating its potential existence of higher order interaction in real data application.\n\n\n|                | SHD | Num of Edges |\n|----------------|:---:|:------------:|\n| NOTEARS-Linear[2] |  22 |      18      |\n| NOTEARS-MLP[3]    |  16 |      13      |\n| DAGMA-HO(Ours) |  13 |      15      |\n\n> As noted by the authors, the full higher-order feature map scales exponentially with the number of nodes, limiting its applicability to larger networks. The NOTEARS-MLP-REG heuristic helps in practice but lacks theoretical guarantees, making its connection to the paper’s core identifiability results less clear.\n> \n\nThank you for raising this important point. Our manuscript is primarily concerned with the **identifiability** of discrete‐valued DAGs—i.e., with characterizing when the global minimizers of (13) recover the true minimal‐equivalence class under a general binary model. We deliberately *abstract away* the algorithmic question of how to solve (13) at scale, since:\n\n1. **DAG learning is NP‑hard** in general, and most existing methods resort to heuristics without strong guarantees.\n2. Our focus is on the properties of the solution set, not on optimization algorithms.\n\nThat said, we agree that addressing computational tractability is a crucial—and distinct—research direction. **NoTEARS‑MLP‑REG** represents our first attempt to bridge this gap in practice: by utilizing the neural network to recover graph and extracting a node ordering and then fitting penalized regressions up to second order, we avoid an exponential search over DAG space. Although it does not yet come with a formal consistency proof, it demonstrates that moderate‐size problems are solvable in reasonable time.\n\nWe’ll add a paragraph in the discussion to: \n\n1. Emphasize that identifiability and efficient optimization are orthogonal challenges.\n2. Position NoTEARS‑MLP‑REG as a *heuristic prototype*, analogous to other NP‑hard DAG‑learning schemes (e.g.,GES, NOTEARS, DAGMA).\n3. Highlight the open question of designing **provably** efficient algorithms that trade off between (i) restricted model space (for identifiability) such as Assumption A (only first order interaction is allowed) and (ii) polynomial‑time solvability.\n\n> The notation in Section 4 is difficult to follow. Simplifying or redesigning symbols would greatly improve readability, e.g., by reducing super- and subscripts.\n> \n\nWe apologize for the heavy notation in Section 4. In the revised manuscript, we will streamline our symbols by reducing unnecessary superscripts and subscripts, especially in Eq (6) where the natural parameter is defined. Concretely, we introduce \n\n$$\n\\mathcal{T}_{\\pi,j} = \\\\{T\\subseteq\\\\{\\pi(1),\\ldots,\\pi(j-1)\\\\} \\\\}\n$$\n\nand then write the parameter vector as  \n\n$$\n\\boldsymbol{f}\\_{\\pi,j} =\\left(f\\_{T\\cup \\{\\pi(j)\\} }\\right)\\_{T\\in\\mathcal{T}_{\\pi,j}}\\in \\mathbb{R}^{2^{j-1}}\n$$\n\nor equivalently  \n\n$$\n\\boldsymbol{f}\\_{\\pi,j}  = \\left(f\\_{S}\\right)\\_{S\\subseteq \\\\{\\pi(1),\\ldots,\\pi(j)\\\\},\\pi(j)\\in S }\n$$\n\nLikewise, we define the feature vector  \n\n$$\n\\mathbf{\\Phi}\\_{\\pi,j}(X) = \\left( B_S(X)\\right)\\_{S\\in \\mathcal{T}\\_{\\pi,j}\\cup \\\\{\\pi(j)\\\\} }\n$$\n\nUnder this notation, the conditional probability becomes  \n\n$$\nP(X\\_{\\pi(j)}=1\\mid X_{\\pi(1:j-1)})  =\\mathrm{logistic}(\\boldsymbol{f}\\_{\\pi,j}^\\top\\mathbf{\\Phi}\\_{\\pi,j}(X))\n$$\n\nA similar simplification applies to each $\\mathbf{H}_j$. We will also include a simple illustrative example to ensure full rigor while greatly improving clarity.\n\n> Visualization: Please keep the color-coding aligned for all figures in the paper! For example, FEGS was marked as purple in figure 1 but orange in figure 2.\n> \n\nWe appreciate this suggestion and will ensure that all color mappings remain consistent.\n\n\n\n> What is the computational complexity of the proposed DAGMA-HO method, particularly with respect to the number of nodes and the order of interactions?\n> \n\nThis question is quite interesting. Suppose we have $n$ samples and $p$  variables, and include all interactions up to order $k$, so that the total feature dimension is $\\sum_{i=0}^k \\binom{p}{i}\\approx p^k$. Computing the gradient of the regularized score function then costs $O(n\\cdot p\\cdot p^k)$ for the score‑function terms plus $O(p^2\\cdot p^k)$ for the quasi‑MCP penalty, which together give $O\\bigl((n+p)\\cdot p^{k+1}\\bigr)$. The gradient of the acyclicity constraint $h\\bigl(W(H)\\bigr)$ requires $O\\bigl(p^3 + p^2\\cdot p^k\\bigr)=O\\bigl(p^{k+2}\\bigr)$, so a single gradient evaluation remains $O\\bigl((n+p)\\cdot p^{k+1}\\bigr)$. If the DAGMA uses $T_{\\rm out}$ outer updates (to decrease the penalty weight) and $T_{\\rm inner}$ inner iterations (e.g. Adam steps), the total complexity becomes $O\\bigl(T_{\\rm inner}\\cdot T_{\\rm out}\\cdot (n+p)\\cdot p^{k+1}\\bigr)$. In the extreme case where all higher‑order interactions are included (i.e. $k=p$), replacing $p^k$ with $2^p$ yields a total cost of $O\\bigl(T_{\\rm inner}T_{\\rm out}(n+p)p 2^p \\bigr).$\n\n> Is it always beneficial to model higher-order interaction terms? How should one choose the order of interaction terms in practice? Is there a principled way to select this, or is it treated as a hyperparameter?\n> \n\nIn theory, allowing every possible interaction can only increase the expressiveness of the model—and recovering the true DAG would benefit from including arbitrarily high orders. In practice, however, the number of interaction terms grows combinatorially in the number of nodes $p$, making full higher‐order expansions infeasible except for very small $p$.\n\nEmpirically, we and others have found that truncating at second order strikes a reliable balance: it captures most of the “nonlinear” dependencies you’d see in real data, yet keeps the total feature count—and hence both memory and runtime—within practical limits.\n\nWhen $p$ is very small (say under 10), exploring up to third or even fourth order can be done directly; once $p$ grows beyond a dozen or so, second order is almost always the sweet spot. \n\n[1] Sachs, Karen, et al. ”Causal protein-signaling networks derived from multiparameter single-cell data.” Science 2005.\n\n[2] Zheng, Xun, et al. \"Dags with no tears: Continuous optimization for structure learning.\" *Advances in neural information processing systems* 31 (2018).\n\n[3] Zheng, Xun, et al. \"Learning sparse nonparametric dags.\" *International conference on artificial intelligence and statistics*. Pmlr, 2020."}, {"review_id": "ur6CSmdTGE", "reviewer": "Reviewer_WmEn", "summary": "The paper proposes a general differential structure learning framework and continuous DAG learning algorithm for discrete data that could capture arbitrary dependencies with no distribution restrictions. To address the combinatorial search problems in discrete data, the paper proposes Notears-MLP-Reg method to reduce computational complexity.", "questions": "See above.", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper is very well-written with clear structured theoretical investigation and empirical evidence. The novelty is to spot the equivalence of modelling binary data as exponential families with interaction terms as potential candidates for causal parents. \n\nComments: \n- in the paper there are multiple places referring that general discrete models are unidentifiable from purely observational data alone (e.g. in abstract line 8, line 218, line 222) and there are also places restricting to the general multivariate Bernoulli distribution (see line 62, line 245) that indicates single environement /  distribution rather than observational data alone. However, it is well-known that multi-environment observational data alone could identify causal structure [1, 2]. To make the paper more precise, it would be good to clarify the differences. \n\n- Theorem 3 says there exists $\\lambda$, $\\delta$ sufficiently small such that the predicted set is the same as minimal equivalence set. Could the authors discuss in practice (when performing experiments) how to select $\\lambda$, $\\delta$. Is it similar as the theorem suggests, the smaller the better?\n\n- The natural problem in discrete data is combinatorial search. The authors address this problem via proposing NoTEARS-MLP-Reg - which first extracts topological ordering and then fit interaction features consistent with $\\pi$ and restrict to first and second order. Could the authors clarify, under this setup, how did the number of searches scale in terms of number of node increases?\n\n[1] S. Guo, V. Toth, B. Schölkopf, and F. Huszar. Causal de Finetti: On the Identification of Invariant Causal Structure in Exchangeable Data. In A. Oh, T. Neumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 36463–36475. Curran Associates, Inc., 2023a\n\n[2] Huang, B., Zhang, K., Zhang, J., Ramsey, J., Sanchez-Romero, R., Glymour, C., & Schölkopf, B. (2020). Causal discovery from heterogeneous/nonstationary data. Journal of Machine Learning Research, 21(89), 1-53.", "was_revised": true, "final_justification": "I believe my concerns are resolved especially around computational time.", "rebuttal": "We thank the reviewer for their time and effort, and for the positive feedback on our paper’s clarity, theoretical guarantees, and empirical evaluation. We will incorporate the following discussion into the revised manuscript.\n\n> in the paper there are multiple places referring that general discrete models are unidentifiable from purely observational data alone (e.g. in abstract line 8, line 218, line 222) and there are also places restricting to the general multivariate Bernoulli distribution (see line 62, line 245) that indicates single environment / distribution rather than observational data alone. However, it is well-known that multi-environment observational data alone could identify causal structure [1, 2]. To make the paper more precise, it would be good to clarify the differences.\n> \n\nThank you for raising this important point. We will revise the manuscript to make it clear that our identifiability results pertain strictly to the single‑environment, purely observational setting—under which general discrete models remain non‑identifiable. We will then contrast this with the multi‑environment case (e.g. data collected under different interventions), in which causal structure can be identified, as shown in [1,2]. In the updated draft we’ll explicitly distinguish these two scenarios and discuss the relevant literature. \n\n> Theorem 3 says there exists $\\lambda,\\delta$, sufficiently small such that the predicted set is the same as minimal equivalence set. Could the authors discuss in practice (when performing experiments) how to select $\\lambda,\\delta$, . Is it similar as the theorem suggests, the smaller the better?\n> \n\nThank you for this insightful question. The related discussion is included at Appendix D.1 (Lines 1088–1095) in supplementary file due to space limit. In theory, any$(\\lambda,\\delta)$ below the threshold specified in Theorem 3 will recover the minimal equivalence set. In practice, however, choosing overly small values can trigger numerical‐precision issues, causing the quasi‐MCP penalty to under‑regularize and produce an overly dense graph. To strike the right balance, we employ a simple “weight‑decay” method proposed in [3]: \n\n- **Initialization**: Start with relatively large $\\lambda,\\delta$ and solve (13) to obtain an initial estimate.\n- **Decay**: Multiply both $\\lambda$ and $\\delta$ by a factor $\\gamma<1$ and re‑solve, using the previous solution as a warm start.\n- **Stopping criterion**: Stop once the negative log‑likelihood ceases to decrease.\n\nEmpirically, this procedure yields sparse, stable, and accurate graph estimates.\n\n> The natural problem in discrete data is combinatorial search. The authors address this problem via proposing NoTEARS-MLP-Reg - which first extracts topological ordering and then fit interaction features consistent with and restrict to first and second order. Could the authors clarify, under this setup, how did the number of searches scale in terms of number of node increases?\n> \n\nThank you for raising this important question. NoTEARS‑MLP‑Reg is a two‑stage procedure. In the first stage, we train an MLP (NOTEARS‑MLP [4]) to capture higher‑order interactions in discrete data and recover an initial graph structure. We then extract a node ordering $\\hat{\\pi}$ from this preliminary graph and, in the second stage, fit a logistic regression model—including both first‑ and second‑order terms—to recover the final graph. Since logistic regression under a fixed ordering is very efficient, the MLP‑based stage remains the dominant computational cost. Empirically, the overall runtime of NoTEARS‑MLP‑Reg closely matches that of NOTEARS‑MLP [4]. The table below summarizes these runtimes and will be included in the final manuscript.\n\n\n|Times(second)| NOTEARS-MLP-REG | \n|:-----------:|:---------------:|\n| ER2 (d=10)  | 43.223          | \n| ER4 (d=10)  | 38.972          |\n| SF2 (d=10)  | 37.925          | \n| SF4 (d=10)  | 41.259          |\n| ER2 (d=20)  | 84.392          | \n| ER4 (d=20)  | 87.267          |\n| SF2 (d=20)  | 80.952          | \n| SF4 (d=20)  | 93.136          | \n| ER2 (d=40)  |     295.333     | \n| ER4 (d=40)  |     337.205     | \n| SF2 (d=40)  |     293.903     |  \n| SF4 (d=40)  |     313.349     |\n\n\nFrom the reported runtimes, it is clear that our method does not exhibit exponential time complexity and can deliver meaningful discoveries within a reasonable timeframe.\n\n[1] S. Guo, V. Toth, B. Schölkopf, and F. Huszar. Causal de Finetti: On the Identification of Invariant Causal Structure in Exchangeable Data. In A. Oh, T. Neumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 36463–36475. Curran Associates, Inc., 2023a\n\n[2] Huang, B., Zhang, K., Zhang, J., Ramsey, J., Sanchez-Romero, R., Glymour, C., & Schölkopf, B. (2020). Causal discovery from heterogeneous/nonstationary data. Journal of Machine Learning Research, 21(89), 1-53.\n\n[3] Deng, Chang, et al. \"Markov equivalence and consistency in differentiable structure learning.\" *Advances in Neural Information Processing Systems* 37 (2024)\n\n[4] Zheng, Xun, et al. \"Learning sparse nonparametric dags.\" *International conference on artificial intelligence and statistics*. Pmlr, 2020."}]

Author rebuttal:

Meta-review: This paper introduces a differentiable structure learning framework for
general binary data, moving beyond the restrictive parametric assumptions
of prior work to capture arbitrary higher-order interactions. The work is
grounded in a rigorous theoretical analysis, starting from the general
multivariate Bernoulli distribution to first prove non-identifiability and
then characterise the complete set of compatible causal structures.

The reviewers initially identified this as a borderline case, raising
several valid concerns. The primary issues included the method's
exponential complexity and questionable scalability, an evaluation limited
to synthetic data, and the reliance on the Sparsest Markov Representation
(SMR) assumption for identifiability.

However, the authors provided a comprehensive and convincing rebuttal that
successfully addressed these key points. They clarified the practical
performance of their proposed heuristic with detailed runtime comparisons,
demonstrating its feasibility. Crucially, they added a new experiment on a
well-known real-world biological dataset, which provided strong evidence of
the method's practical utility and directly addressed a major limitation.
The authors also effectively justified their use of the SMR assumption and
committed to improving the paper's clarity. The productive discussion and
the authors' substantive responses led multiple reviewers to increase their
scores or strengthen their support for the paper.

Given the strong theoretical foundation, the novelty of addressing
structure learning in a general binary setting, and the new empirical
evidence demonstrating its real-world applicability, the paper now stands
as a solid contribution. Hence the AC recommend acceptance.
Decision: Accept (poster)