Inexact Column Generation for Bayesian Network Structure Learning via Difference-of-Submodular Optimization
Authors: Yiran Yang
Rui Chen
Venue: NeurIPS 2025 poster ()
Keywords: Causal Discovery
Integer Programming
Column Generation
Primary area: probabilistic_methods

Abstract: In this paper, we consider a score-based Integer Programming (IP) approach for solving the Bayesian Network Structure Learning (BNSL) problem. State-of-the-art BNSL IP formulations suffer from the exponentially large number of variables and constraints. A standard approach in IP to address such challenges is to employ row and column generation techniques, which dynamically generate rows and columns, while the complex pricing problem remains a computational bottleneck for BNSL. For the general class of $\ell_0$-penalized likelihood scores, we show how the pricing problem can be reformulated as a difference of submodular optimization problem, and how the Difference of Convex Algorithm (DCA) can be applied as an inexact method to efficiently solve the pricing problems. Empirically, we show that, for continuous Gaussian data, our row and column generation approach yields solutions with higher quality than state-of-the-art score-based approaches, especially when the graph density increases, and achieves comparable performance against benchmark constraint-based and hybrid approaches, even when the graph size increases.

Inexact Column Generation for Bayesian Network Structure Learning via Difference-of-Summoduaa Optimization
Yiran Yang
School Data Science
The Chinese University of Hong Kong, Shenzhen yiranyangol ink. cuhk edu. en
Rui Chen
School of Data Science
The Chinese University of Hong Kong, Shenzhen rchen@cuhk.. het.eccukkkuu..cc..
Abstract
In this paper, we consider score-based Integer Programming (IP) approach for solving the Bayesian Network Structure Learning (BNSL.) problem. State-of-the- art BNSI. IP formulations suffer from the exponentially large number of variables and constraints. standard approach ploy and column generation techniques, which dynamically generate rows and columns, while the complex pricing problem remains a computational bottle- neck for BNSL. For the general class of fu-penalized likelihood scores, show how the pricing problem can be reformulated as a difference submodular op- timization problem, and how the Difference of Convex Algorithm (DCA) can be applied as an inexact method to efficiently solve the pricing problems. Empir- ically, show that, for continuous Gaussian data, tion approach yields solutions with higher quality than state-of-the-arr score-based approaches, especially when the graph density increases, and achieves compara- ble performance against benchmark constraint-based and hybrid approaches, even when the graph size increases.
to address such challenges to em-
and column genera-
Introduction
Bayesian networks (BN) have wide applications in machine learning 27 Its structure learning problem, i.e., identifying directed acyelic graph (DAG) that represents the causal relationships between variables (X))ev from data, is central yet challenging topic causal inference (29]. Ex- isting methods for learning optimal Bayesian network structure from observational data fall into three main categories (45, 40 21]:: (I) constraint-based approaches (e.g., PC [55. 14] and FCI 145 50) that learn the structure through a constraint satisfaction problem and involve conditional independence tests to formulate the constraints; (2) score-based approaches (e.g., hill climbing (28). dynamic programming [35 44]], integer programming [17] and machine learning approaches (49]) that evaluate candidate DAGs using scoring functions and select the DAG structure with the highest score; (3) hybrid approaches that combine both score-based and constraint-bast techniques (e.g..
MMHC [47..
Score-based BNSL is known to be NP-hard, even under the assumption that each node has re- stricted number of parents [13]. This inherent hardness of BNSL imposes fundamental trade-off
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

between solution quality and computational efficiency. Consequently, there are two lines of work that focus on the either side: (I) exact methods that guarantee optimality but often suffer from higher computational cost and (2) heuristic approaches that offer computational efficiency but with limited optimality guarantees. Within the BNSI. literature, one of state-of-the
exact approaches to iden- tify the highest-scorine DAG is to formulate the problem as an Integer Program (IP) [31, 15. 12), as represented by the GOBNILP method [17]. While there are different ways to formulate the IP for this problem, standard formulation involves binary variables that indicate whether to choose the specific set as the parent set for each node, and constraints that enforce the acyclicity requirement for the DAG. Given the observational data, the objective of the IP is to maximize the score, which is a function of the potential DAG structure measuring how well the data fits the model. Solving this IP is challenging since both the number of variables and the number of constraints grow exponentially with the node size of the DAG.
Column Generation (CG) 26 199 is a standard approach for addressing huge IPs as it effectively hal- ances computational efficieccy and solution quality. By dynamically generating variables as needed through pricing problems, CG can be used to derive efficient heuristic algorithms (32]] or exact so- lution methods (often known as branch and price [S]) when combined with branch and cut. CG has been widely applied in operations research to solve real-world problems such as vehicle routing [20] and crew scheduling [41]. In [ 18), Cussens introduces CG to BNSL to construct an exact branch and price algorithm. However, due to the complex pricing problems (which are formulated as mixed integer nonlinear programs in (18)), the proposed approach has limited scalability and can only be applied to solve small BNSL instances.
Indeed, to accelerate CG algorithms, it is advised that one should not try to solve all pricing prob- lems exactly and an approximate solution often suffices until the last few iterations (for proving optimality) [37]. In this paper, built on top [8). we advance CG solution techniques for BNSL. We develop efficient pricing algorithm based on reformulation of the pricing problem as differ- ence of submodular (DS) optimization problem, enabling efficient solution of the pricing problem via the Difference of Convex Algorithm (DCA) (22 30]. Even though DCA does not necessarily pro- vide global optimal solutions, it guarantees convergence and often generates high-quality candidate columns for BNSL IP. Empirically, we demonstrate that it leads to effective primal heuristics, offering potential for developing more efficient exact solution approaches.
2 Preliminaries
2.1 Scored-based BNSL IP
In BNSL, a scoring function measures how well the observed data fits the DAG structure, where nodes represent individual random variables and directed edges represent direct causal dependencies. Given observational data D, score-based approaches aim to find a DAG with the highest score, i.e., solve the following optimization problem:
max score D) s.t. G being acyclic,
G
where score (G; D) denotes the score of the graph G = (V, E) under data D (see Appendix A.I. Given a partieular data type (continuous or discrete), there exist several different choices for the scoring function. In this paper, we focus on lo do-peallieed likelihood scores, i.c., scoring functions of the form sCore(G; D) = log(L(G; higher maximum likelihood L(G; D) and a lower number of free parameters to be estimated k k(G) in the graphical model. By setting A differently, one can recover some commonly used scoring functions, such as the AIC score 1 and the BIC score [42].
A k(G) for some A 2 O, which favor graphs with a
A crucial property of fo-penalized likelihood scores is their decomposability into node-specific components. Specifically, assuming that we have variables in our BNSL problem indexed by i V (I,...,n). the score score(G; D) of DAG G defined over nodes V satisfies
(G; = Σ score; (pa,
where (care;(pa;(G)) denotes the local score of node (where we omit its dependence on for simplicity), whose value only depends on the parent set pa,(G) of node i in graph G, i.e., the set

---

of nodes pointing towards node i in the directed graph G. Specific definitions of the Eo fo-penalized likelihood score and its local version, in both discrete and continuous data cases, can be found in Appendix A.I
The additive nature of the scoring function enables us reformulate ① as an integer (linear) pro-
gram. A classic IP formulation based on cluster constraints [15. 31 is as follows
max
(core;()) Ii+J
(2a)
JEP.
s.t.
Σ æitJ= 1,
i= 1=
(2b)
JEP:
Σ E
TieJ S -1 CeC,
(2c)
I JEP
JnC#0
LieJ E {0,1),
i= in JEP,
(2d)
Here Pi: == 2V((+) denotes the set of candidate parent sets of node i, the binary decision variable Jie-J indicates whether parent set J is selected for node i, with associated local score score:(J).
The linear objective reflects the score of the graph composed of the selected parent sets. Constraint
(2b) enforces exactly one parent set be chosen per node, while cluster constraints (2c) [31 guar- antec acyclicity by preventing directed cycles among any nonempty node subset C' E = 2 {0}.
Note that formulation 22 involves @(n2*-1) binary variables and ( constraints. Due to its
exponential size, even directly writing down 2 as an IP in a computer infeasible for reasonably large n, not say solving it. Although the exponential number of constraints can be addressed through row generation (i.e., replacing C by small subset C in (2) and iteratively adding violated ones into č, see Appendix B [15], the exponential number of variables remains significant chal- lenge. To deal with exponentially many variables, existing BNSL IP approaches often rely on a (compromising) assumption that the size of the parent set for each node is no larger than some con- stant k2 in which case one may replace P. by {J E P, [|] S k} in (2 for all to reduce the number of variables [6]. Instead of putting such a restriction on the parent sets, we consider CG for BNSL as proposed in 18)
2.2 Column Generation and the Pricing Problem for BNSL
Rather than including all the variables, CG handles IPs with huge number of variables by dynami- cally generating only the necessary ones through the so-called pricing problems. In the context of the BNSL IP (2), as demonstrated in (18), the CG procedure begins with restricted set of parent sets P. S P, for each node 1,..., n, and considers restricted version of 2 only over vari- ables TieJ; for i 1,..., n, E P,,, i.e., ( with small subset of columns. We refer to the linear programming (LP) relaxation of this restricted problem as the Restricted Master LP (RMLP). Based on the current RMLP solution, CG iteratively searches for additional parent choices that may improve the objective value through the solution of pricing problem for each node i, and add them to P.. This pricing problem aims to optimize the so-called reduced costs, to identify high-quality parent set choices to gradually enlarge the variable space and solve the LP relaxation of 2
As pointed out in 118, the pricing problem for generating candidate parent sets for node can be formulated as the following optimization problem:
min =(J;²*) score, (J) + 12+ + JEP.
Σ
X
CEC: aEC, JnC
where Xj and NO 2 O are the optimal dual values associated with constraints 2b and 2c in the RMLP, respectively. The value of (J; X*) is known as the reduced cost { of the variable JæeJ- Its negative measures how much the RMLP objective value might increase by adding the column associated with _i.J to RMLP. The column generation process terminates when the optimal reduced cost 33) becomes nonnegative for all the nodes. In this paper, instead of using CG to construct an exact solution approach like branch and price, we consider straightforward CG-based heuristic named the restricted master heuristic (32), where solve the IP 2 with the candidate set P. generated by CG, namely the Restricted Master IP (RMIP), after CG (approximately) terminates.
3

---

Algorithm The Difference of Convex Algorithm
1: Input: Convex functions g", fr. (and 2 = - initial point x threshold
2: Initialize t 4-0
3: repeat
4: y¹ E E5f((²))
5:
gd+1 - arg min (x) x))}
6:
++++1
7:
until 1x"(x") 2 (2--1) <E
8:
Return Solution x², objective value 24 (x²)
2.3 Difference of Submodular Optimization via DCA
The pricing problem 3 is a set function optimization problem. Although problem (3) in the con- tinuous Gaussian case, can be formulated it with an MINLP solver is highly computationally demanding as shown in [18]. Taking a direct set function optimization perspective, it is known that any set function can written difference of submodular (DS) functions although finding such DS decomposition has exponential complex- ity [30]. However, when such decomposition known, one can develop efficient algorithms to find local solutions [30, [22]. Here we briefly review a DS optimization algorithm [22 based on the well-known DCA (2)
mixed-integer nonlinear program (MINLP), solving
A DS function = 2 +
is a set function that can be expressed as:
z(J) = gJ) fJ).
where both g 2V -> I and f are submodular set functions (see Appendix A.2
25
For set
function h 2 + IR, its minimization can be equivalently formulated through its Lovász extension
ht [0. 1]" - R (see Appendix A (36). The Lovász extensions for submodular functions
convex [36), enabling the reformulation of min, x(J) = g(J) - fJ) as a difference of convex
(DC) program
min 2 2²(x) gz fx))..
To solve this DC optimization problem, DCA (Algorithm I iterates between two key steps. First, it computes linear approximation (z(x) (2)) of f (z) at the current iteration point x² defined by
(z) fz(²)) (V', =x'),
where yf is a subgradient (*x at x² Subsequently, it minimizes a convex approximation of the objective (x),, with f*)) replaced by fz)) This subproblem is solved to obtain the next iteration point zf+1
As the subgradient of the Lovász extension of a submodular function can be efficiently evaluated through its function evaluation, the convex program (Line[5) in Algorithm[[can be efficiently solved via bundle methods [4]. In our case, we adopt the classic Kelley's Algorithm (33] for numerical experiments.
One important advantage of DCA is that guarantees a non-increasing sequence objective values, ie., for t 2 0, we have
z'(x++1 < 2z (x").
is also easy to recover solution of the original set function optimization problem from DCA, as the function value of agrees with - integer solutions. Since the convex function g (v':
the Lovász extension (which characterizes the convex envelope) of a submodular function (i.e., plus modular function), all extreme point solutions of the convex program (Line 5 in Algorithm are integer solutions.
Solution of the Pricing Problem
3.1 The Pricing Problem DS Optimization
key step CG is the solution of the pricing problem 3 that iteratively selects "promising" columns to add into the restricted problem formulation. Existing MINLP approaches are known

---

to have very limited scalability to solve the pricing problem for BNSL [8] We show in this section
how the pricing problem can be explicitly rewritten as DS optimization problem, enabling us to take advantage of the DCA algorithm we describe in Section 2.3
Recall that the pricing problem for node is formulated as minimizing the reduced cost
(²):X*) score, + JEP,
min
Σ
X
GEc: SEC, JnCØE
For both continuous and discrete cases, we show how the to fo-penalized likelihood score 2x();²* can be expressed as a DS function.
Proposition 1 (Continuous Case). For continuous data with the to-penalized Gaussian likelihood score, 2,(J; X*) satisfies
N log(2n 1);
= -
2
where
N log det (e.(,.uu,,J))) + A|J1+
g(J)
Σ
No + x;
Cec. ÇEC, JnCØ0
and
N log det(E...)
[(J)
are both submodular functions, with Σg and Jos denoting the empirical covariance matrices associated with variables (X3 j J} and {X3 1 J U {i}}, respectively.
Proposition 2 (Discrete Case). For discrete data with the lo-penalized multinomial likelihood score, x*) satisfies
=xJ=A²) = 9( - f()),
where
= N- ((JU E
2 +
OEC: iEC, JNCYd
and
F(J) = N. N-H - A(a; - 1) II aj
jEJ
are both submodular functions, functioon,,, with HJJ and (JUU)) denoting the joint entropy associated with variables {X, j E J} and (X, j E J U }i)}, respectively, and 0j denoting the arity (i.e., the number of possible values it can take) of the variable X, for 1,
Proofs of Propositions D and can be found in Appendix C
3.2 Implementation of DCA for Pricing
We have reformulated the pricing objective DS function in Propositions and 2
As established in Section 2 minimizing the pricing objective is then equivalent to minimizing the difference of convex functions ==xx) g g (z) fz (z), where g0(2) and f*(x) represent the Lovász extension of submodular functions and f. respectively.
To optimize this objective, we apply DCA. i.e., Algorithm M to iteratively minimize (xr) through its successive convex approximations. The DCA procedure for solving the pricing problem for node i proceeds follows. The algorithm begins by initializing vector (iy))evv (i]- The required inputs include the dataset D R xn. node index é, optimal dual solutions X of RMLP, regularization parameter A, and convergence thresholds e.
(n -iddieensional solution
During iteration I, the algorithm first computes the subgradient y/ of fh at the current solution x.. Without loss of generality, assume for simplicity that Let a be permutation of (I...., n-1,
5

---

such that (1) 2 (2() ....
x (n--)--
For the continuous case, subgradient y² of fr. at x² is
given by
log det( (o(()()(1)))
k =
Vo(u) =
log det( (et(E((()))...( log det (e()) k 2 2,
where ((k)) denotes the empirical covariance matrix for the corresponding columns (o((),.... (k)) of the data matrix D E RgNx*.
To improve computational efficiency, we compute subgradients using Cholesky decomposition. By the lemma below, for any fixed permutation o, we can efficiently calculate log (e).(_())..... (m < n 1) directly from the Cholesky factor L of the permuted matrix (I( = LL
This approach requires only a single Cholesky decomposition of the permuted covariance matrix to obtain all necessary subgradients for the Lovász extension, rather than explicitly computing log det (S)) for (n) times.
Lemma Let M be a p x p positive definite symmetric matrix with lower-triangular Cholesky factor L. In block form, these matrices can be written as
M11
M12
L1 0 L2] --
M
L
M21
M22
where 11 and Li are blocks (m < p) and O is zero matrix of uppropriate size. Then L the Cholesky factor of M1,, and consequently,
m log Lk,k-
log (e1)
k=1
For the discrete case, a subgradient yt fr at xt is given by
(((((1))) A(a. ()(0.(1)) 1) (k)) (((a(1).... (((
k k=
1)}) A(a. 1)(a(()) 1)11 Go()
k 2.
After evaluating log ()) (-) during the computation of y is crucial to save the evaluated values to avoid repeating calculations in the future. This small trick greatly improves the efficiency of DCA in pricing problem solutions in practice.
4 Row and Column Generation Scheme
Algorithm demonstrates how we incorporate our DCA pricing method into a simultaneous row and column generation framework. The algorithm begins by initializing the candidate parent sets Pi.. and the cluster set C with basic candidates (Line 2). It then iterates among three main phases: (1) CG phase (Lines[4[44) where the RMLP is solved to obtain optimal dual solutions X", followed by the solution of the pricing problem for each node via DCA The pricing solutions that have negative reduced costs added to P. The CG phase repeats until all reduced costs become non- negative. In Line g the pricing problem on node is solved only when the current reduced cost rc, is negative, to avoid spending too much time on pricing; (2) a row generation phase (Lines 15 where the RMLP optimal solution is used to identify the most violated cluster constraint (see Appendix B for details) which is added to cluster set C; (3) the integer solution phase (Line 22). where the algorithm finds the optimal DAG G with current candidate Pi.n - This is achieved by solving RMIP using branch and cut with callbacks (a feature that allows users to interrupt the solution process in modern IP solvers such Gurobi), where check if any cluster constraint is violated when we meet an integer incumbent and add violated cluster constraints as lazy constraints. The three phases iterate until convergence to a valid DAG solution, progressively refining both the columns in P and cluster constraints in Ć.
When we implement the CG phase in Algorithm in addition to the final solution generated by DCA, we add to P; all candidate parent sets associated with intermediate x' solutions that have negative reduced costs. Empirically, we find that it helps accelerate convergence and improves solution quality. In addition, we initialize the RMIP with the optimal IP solution obtained from the last round, as it remains feasible and empirically reduces the solution time.

---

Algorithm 2 Row and Column Generation for BNSL
1: Input: Data matrix D IRN*", regularization parameter
2:
Initialize the cluster set é t- 0, and candidate parent sets P, + {0} for i 1,...,7
3: repeat
4: Initialize rc; <- -0 for i = 1,...,7 7
> initialize reduced costs
5: repeat
6: Solve RMLP to obtain the optimal dual values A
7: for i ,,...,77 do
8:
if rc; < 0 then
D selectively solve pricing problems
9:
Choose initial point =d for DCA (pa, rc,) t OCA-Pricing(D), i, x, X ,AA,))
10:
> DCA pricing
11:
P. + P. U (pa)
column generation
12:
end end for
13:
14:
until rc; 2 for 1...n
15:
repeat
16:
Solve RMLP with updated P and Identify the most violated cluster constraint C* if the cluster constraint associated with C* is violated then
17:
> row generation for RMLP
18:
19:
C
20:
end if
21:
until no new violated constraint can be found
22:
G optimal DAG chosen from Pian: simultaneously update ć D row generation for RMIP
23:
until G converges
return G
4.1 Choice of the Initial Point
step that remains ambiguous in Algorithm 2is the choice of the initial point x0 for DCA (Line which we find having significant impact the overall solution time and solution quality. Below we first describe two initialization strategies.
The first strategy is a warm-start initialization strategy, Note that by the complementarity property of optimal solutions in LP (], the parent sets selected (i.e., parent sets E P; associated with the primal variables J taking positive values) in the RMLP solution have reduced costs equal to zero. Therefore, we consider using the parent sets selected in the RMLP solution to construct the initial solution N0 for DCA. Due to the monotonicity property 4 of DCA, it often quickly finds columns with negative reduced costs, leading to rapid convergence of DCA. However, drawback the warm-start initialization strategy is that it may restrict exploration of parent sets as the evaluation remains localized around columns already contained in the candidate parent set choice set. Preliminary experiments show that a pure warm-start initialization strategy can quickly find DAG with good score while the solution often deviates lot from the ground truth due to its greedy nature.
The second strategy employs random initialization that samples point 2: uniformly from [0, ( This random initialization strategy facilitates global exploration on the parent set space and may discover patterns distant from the current ones in P,, However, we find that exclusive reliance on random initialization can significantly slow down DCA convergence, as computational resources may be expended on evaluating solutions far from the optimum.
hybrid initialization strategy that combines efficiency and diversity of the previous two strategies is often more effective in practice. In our implementation, we initially employ random initialization to broadly explore the solution space. Once the candidate set P. contains a reasonably large num- ber (set 50 our implementation) of parent sets, we switch to warm-start initialization for pricing at node i. This hybrid approach focuses local refinement around the current best pattern while building upon the foundation of global exploration. A comparison of the three initialization
strategies can be found in Appendix D.
7

---

Table 1: Score and Time Comparison between CG-DCA and GOBNILP on Gaussian Datasets
(n, Nr d)
Score Copp (%)
Time (seconds)
CG-DCA
GOBNILP
CG-DCA
GOBNILP
(15. 0.5)
(0.00. 0.00) [0.00]
(0.00, 1.46) [0.15]
[
1.12.
14,71)
4,06]
5.47.
24.28) | 8,33
(15. 5000. 1.0)
(0.00, 0.24) [0.03]
(0.00. 1.06) [0.22]
93.58)
25.76]
5.72.
56.11) [
19.74
(15, 5000. 1.5)
(0.00, 0.62) [0.16]
(0.00. 2.18) [0.62]
17.95.
196.01)
76.72]
11,51.
137.52)
60.09)
(15. 5000, 2.0)
(0.00, 0.84) [0.39)
(0:48, 4.88) [2.06]
41.13,
340.22)
170.39]
22.99.
641.90)
167.41]
(15, 20000, 0.5)
(0.00, 0.00) [0.00)
(0.00, 0.64) [0.06]
3.53.
35.11)
10.171
5.53.
12.62)
7.35]
(15, 20000, 1.0)
(0.00, 0.27) [0.03]
(0.00, 1.19) [0.23]
[
7.52,
100.20)
36.71]
1,31,
93.02) [
19.67]
(15, 20000, 1.5)
(0.00,
[0.02]
3.10) [0.86]
25.17.
265,33)
104.25]
12.34.
285.39)
76.29]
(15. 20000, 2.0)
(0.00, 0.94) [0.23]
(0.98, 3.97) 12.011
55.11,
335,96)
229
11.66.
190.19) 1 89.92
(20. 5000, 0,5)
(0.00, 0.01) [0.00]
(0.00, 0.18) [0.02]
20.20)
9.53]
1.36.
6.99)
3.38
(20. 5000, 1.0)
(0.00, 0.08) [0.01]
10.00, 1.38) [0:35]
5.97,
150,
84.34]
1.59,
1714,57)
[ 287.87]
(20. 5000, 1.5)
(0.01, 0.51) [0.22]
(0.00, 4.30) [1.63]
(144,94,
339.72)
312.67]
121.01.
3483,96) 96) 962.50]
(20. 5000, 2.0)
(0.22. 1.15) [0.70]
(0.07, 10.42) 14.771
(403,98,
1904.24)
11039.531
(1063.89,
8703.57)
[4603.81
(20. 20000, 0.5)
(0.00, 0.00) [0.00]
10.00, 0.35) [0.04]
[ 6.57,
36,37)
18.42]
1.50,
11.25) [ 4.92]
(20, 20000, 1,0)
(0.00, 0.22) [0.03]
1.49) [0.22]
13.72,
330.05)
133,69)
1.54,
1371.66) | 261.89]
(20. 20000, 1.5)
(0.00) 0.43) [0.15]
(0.00, 3.17) [1445]
(216.26.
821 [0]]
425.17]
135.59,
7269
7269.60)244 60)2457.261
(20, 20000, 2.0)
(0.09, 1,08) [0.54)
10:21, 7.05) [3,63]
(457.33,
1739.24)
997.99)
517.67,
10801 59) [4394.45]
(25. 5000, 0.5)
(0.00, 0.00) [0.00]
10.00, 0.82) [0.
[ 2.73, 67.05)
21.97]
1.61,
151.24)
21.95]
(25. 5000, 1.0)
(0.00, 0.25) [0.03)
(0.00, 1.83) [0.56]
[ 28.04, 953.40) [ 261.84]
11.14.
7451.77
(1964.11]
(25. 5000, 1.5)
(0.00. 1.11) [0.24]
(0.09, 5.27) [1.92]
(144.19, 2880,93)
949,90)
84.46,
10801
61) [5805.01]
(25. 5000, 2.0)
(0.00, 2.18) [0.86]
(1.84, 18. [5.15]
(448 93, (2665.38)
[3054.63]
723.38,
10801 63) [8500.77]
(25, 20000, 0.5)
(0.00) 0.01) [0.00]
(0.00, 0.78) (0.09]
9.50,
164.57)
50.221
1.62. 169.94) |
(25. 20000, 1.0)
(0.00, 0.22) [0.06]
(0.00, 2.72) [0.40]
70.34, 1287.17)
404.43)
11.86,
10801 48) [1948.05]
(25, 20000, 1.5)
(0.00, 1.19) [0.23]
10.66, 5.47) [2:21]
(237 86 4955.58)
[1372.59]
115.63, 10801.73 73) [5853.91]
(25, 20000, 2.0)
(0.00. 2.71) [0.67]
(1550 8.30) [4.51]
(733.04, 7371.52) (3189.98)
(2302.81.
10801.77 71) [7967.28]
# *The values are the (min, max) ] average] over len independent instances.
5 Numerical Experiments
To evaluate the efficiency of our method, we conducted numerical experiments comparing the per- formance of our Column Generation with DCA-pricing (CG-DCA) method against five baseline methods: GOBNILP 117, CG-MINLP [18] (in Section D. HC (28] (in Appendix (D.2), stable-PC [14] and MMHC [47]. For soore-based approaches, we use the BIC score as our scoring function. While CG-MINLP is evaluated on Gaussian. test dataset in the R package bnlearn [43], other methods are tested on both simulated and real-world datasets. We implement CG-DCA in python, and utilize the python implementation of GOBNILP. The last three baseline methods are tested using bnlearn in R. All IP-related experiments are conducted on Linux machine with two Intel XEON Platinum 8575C processors, with the number of threads used limited to 1. All IPs are solved using Gurobi 11.0.3.
5.1 Comparing DCA with MINLP for Pricing
We first compare the performance of CG with different pricing algorithms: our proposed CG- DCA versus the exact MINLP solver for pricing (CG-MINLP). The evaluation is conducted on the Gaussian. test dataset containing a ground truth graph with 7 nodes and 7 edges. For this small instance, CG-DCA successfully learns the exact ground truth structure in 1.2 seconds, while CG-MINLP requires 19.2 seconds to achieve the same result.
test
For larger instances, CG-MINLP usually fails to produce a valid solution DAG within a reasonable time limit as solving the pricing problem exactly as an MINLP is quite cumbersome in practice. For a graph with node size n = 20 and average in-degree d = 1 (N 5,000, simulated as in Section[5.2) solving the first pricing problem using MINLP takes an average of 35.1 seconds per node, while the DCA method takes 0.3 second per node for the same instance. The difference is even larger for later pricing problems. In addition, MINLP approaches suffer from numerical issues as they require an approximation for the nonlinear log function (at least in Gurobi), producing inaccurrate estimates for the pricing objective.
5.2 Results on Simulated Gaussian Data
Now we test on larger instances with Gaussian data. Following the experimental setup in 18. - we randomly generate Bayesian networks with node size n E (15, 20, 25). and simulate Gaussian data from those networks. To generate the DAGs, we first fix a topological order of the nodes
8

---

Table 2: Performance Comparison of CG-DCA, GOBNILP (GI), (Stable-) PC and MMHC on Gaus-
sian Datasets
(n,
Precision
Recall
SEID
CG-DCA
Gl
PC
MMHH.
CG-DCA
Gl
PCC
MMHC.
CG-DCA
Gl
PC
MMHHC
(15. 5000, 0.5)
0.78
0.66
0.59
0.50
0.80
0.73
0.62
0.51
1.70
2.90
3,00
3.50
(15. 5000. 1.0)
0.70
0.67
0,72
0.59
0.78
0.72
0.67
0.54
6.20
5.00
4.50
6.70
(15.
0.61
11:50
(15. 5000, 2.0)
0.55
0.58
0.54
0.52
0.72
0.65
0.39
0.36
21.90
16.90
20.90
20:60
(15. 20000, 0.51
0.71
0.69
0.68
0.53
0.74
0.73
0.68
0.53
2:20
2.30
1:90
3,00
20000,
0.78
0.71
0.09
4,00
20000,
0.66
(15. 200000 2.0)
0.61
0.55
0.53
0.55
0.75
0.60
0.43
0.39
20.00
19.00
20.01
19.80
(20. 5000, 0.5)
0.74
0.71
0.76
0.29
0.90
0.75
0.79
0.60
3.40
3.50
2.60
4,50
(20. 5000,
0.69
0.75
0.65
0.76
0.77
0.71
0.59
3.70
7.50
5.90
7,90
120, S0NK,
0.61
0.71
28:20
120. 5000, 2.0)
0.64
0.46
0.54
0.55
0.66
0.53
0.34
0.31
36.40
33,00
31 31.30
31.10
(20, 20000, 0.51
0.82
0.80
0.73
0.58
0.83
0.82
0.76
0.59
2.30
2,60
3,00
4,50
20000,
0.75
0.79
(20, 20000,
0.58
0.56
15.70
16.00
(20,
20000, 2.0)
0.52
0.48
0.51
0.56
0.73
0.55
1.57
0.33
35.00
32.40
51. 30
29.30
[25. 5000,
0.82
0.61
0.84
0.09
0.84
0.67
0.84
0.69
2.70
5,50
2,70
4,50
5000,
0.69
0.80
0.76
0.75
0.69
11.10
8.70
(25. 500K,
0.57
0.70
0.71
0.58
13,60
18:00
[25. 5000, 2.0)
0.54
0.47
0.55
0.59
0.73
0.54
1.35
0.35
45.70
40.30
38.70
36.50
(25,
20000, 0.5)
0.76
0,80
0.81
0.70
0.80
0.84
0.83
0.70
3,70
2,90
2,80
4.30
20000,
0.81
0.78
(25, 20000,
0.57
0.53
0.73
0.83
0.51
14.10
(25, 20000, 2.0)
0.54
0.48
0.61
0.59
0.66
0.56
0146
0.40
23.90
39.10
53.50
34.10
Fach value is the overgge over len indepondent instances.
and then simulate edges according to the specified probabilities where edges must be directed from lower-order nodes to higher-order nodes (there are n(n - 1)/2 such possible edges). While we do not enforce maximum in-degree of each node, we randomly simulate edges using their existence probability determining the overall density of the graph, i.e., for each possible edge we include in the network with probability p. The average in-degree d of nodes in the ground truth DAG is controlled by the edge existence probability since = 2d/ (n 1). We vary such that d is varied among {0.5, 1.5.2].
After generating the random graph structure, we simulate Gaussian data based on the network topol- ogy. The linear coefficients for the features randomly drawn from the interval = 2]. with additive noise following the Normal distribution o2). The true value of variance parameter o² uniformly drawn from the interval [0.7, 1.2].
Tables and present a comparison between CG-DCA and baseline methods across various node sizes (n), and average in-degrees (d) with N = 5,000 and N = 20,000. For each graph structure, we generate 10 independent data instances using different random seeds and record the average performance metrics.
Table reports the BIC score gap (difference between the achieved BIC score and BIC score of the true graph) and the runtime (with a 3-hour time limit per instance) of CG-DCA and GOBNILP. As demonstrated in Table the CG-DCA method consistently outperforms GOBNILP terms of average scores across simulated Gaussian datasets. Notably, CG-DCA achieves optimal BIC scores for instances with sparse graph structures (d 0.5). As the number of nodes increases, CG-DCA exhibits slower growth in runtime compared to GOBNILP while maintaining superior scoring performance. The robustness both computational efficiency and solution quality with respect to node count (ri) and graph density (d) makes CG-DCA more suitable for larger problems.
In Table we compare four methods with the following metrics on the quality of the solution graphs: (1) precision and recall [46]. which measure the proportion of true edges among predicted edges and predicted edges among true edges, respectively; (2) the structural Hamming distance (SHD) [47). which quantifies the dissimilarity between two graphs by counting the required edge additions, dele- tions, or reversals to make one graph identical to another (3). Prior to computing these metries, we convert both predicted and true graphs into their corresponding essential graphs [3], which represent their Markov equivalence classes. The conversion of essential graphs is performed following the methodology described in [38]. We also have results for HC but only present them in Appendix D.2 due to its inferior performance. Notably, CG-DCA achieves better average recall across most instances while maintaining competitive precision and SHD compared to other constraint-basee and

---

Table 3: Performance Comparison of CG-DCA, GOBNILP (GI), (Stable-)PC, and MMHC on Dis-
crete Datasets
Precision
Recall
SHD
Daraset CG-DCA
Gl
PC
MMHC
CG-DCA
Gl
PC
MMHC
CG-DCA
Gl
PC
MMHC
LUCAS
0.77
1,00
0.83
0.33
0.83
1.00
0.83
0.33
2
INSURANCE
0.53
0.90
0.91
0.44
0.83
0.58
0.35
37
34
ALARM
0.46
0.87
0.50
0.89
0.76
0.33
35
32
hybrid approaches. Better performance recall than on precision and SHD indicates that CG-DCA tends to select denser graphs than others.
5.3 Results on Discrete Data
Table 3 compares four methods discrete datasets LUCAS (3 with (n, d. (12, 2000, 1). INSURANCE (IO] with (n, N, d) = (27, 20000, 3.85) and ALARM [7] with (n, NN, d) (37, 20000, 2.49). GOBNILP demonstrates superior performance in all three discrete datasets.
Although CG-DCA an option for small discrete dataset, performance degrades with larger net- works. This is because the discrete BIC score contains highly supermodular penalization term (a; 1) aj, making DCA less effective for DS optimization. Thus, development of an im- proved decomposition strategy for discrete scores (potentially other scores like BDeu [1] is essen- tial to broaden the applicability of DCA to larger discrete Bayesian networks.
6 Conclusion
In this paper, we propose CG-DCA, method that leverages difference-off -otprrencce-o--uuuuuuuuuul minimization to solve the pricing problem within CG framework for BNSL. Empirical results demonstrate that CG- DCA outperforms statc-of-the-an score- score-based methods on simulated Gaussian datasets with varying node sizes, sample sizes, and graph densities, yielding solutions of high quality.
While CG-DCA is viable approach for small graphs with discrete data, its scalability larger graphs remains limited due to the high supermodularity induced by the penalty term in the discrete BIC function. To mitigate this limitation, one could explore alternative decomposition strategies for (potentially different) discrete scoring functions develop new DS optimization techniques tailored to this specific computational challenge. Furthermore, building upon our efficient pricing heuristic, future research could be focused developing exact pricing algorithms to enhance both the con- vergence guarantees and solution optimality of CG, by leveraging exact approaches for submodular [36] and supermodular (39] optimization.
10

---

Acknowledgements
This work was supported by the National Natural Science Foundation of China under Grant 72501249.
References
[1] H. Akaike. A new look at the statistical model identification problem. IEEE Transactions Automatic Control, 19:716, 1974.
[2] L. H. An and P. D. Tao. The DC (difference of convex functions) programming and DCA re- visited with DC models of real world nonconvex optimization problems. Annals of Operations Research, 133:23-46, 2005.
[3] S. Andersson, D. Madigan, and M. D. Perlman. characterization of markov equivalence classes for acyclic digraphs. The Annals of Statistics, 25(2):505- 5511, 1997.
[4] A. Bagirov, N. Karmitsa, and M. M. Mäkelä. Introduction to Nonsmooth Optimization: Theory, Practice and Software, volume 12. Springer, 2014.
[5] C. Barnhart, E. L. Johnson, G. L. Nemhauser, M. W. Savelsbergh, and P. H. Vance. Branch-and- price: Column generation for solving huge integer programs. Operations Research, 46(3):316 329, 1998.
[6] M. Bartlett and J. Cussens. Integer linear programming for the bayesian network structure learning problem. Artificial Intelligence, 144:258-271, 2017.
[7] Beinlich, J. Suermondt, M. Chavez, and G. Cooper. The alarm monitoring system:
A case study with two probabilistic inference techniques for belief networks. In AIME 89: Second European Conference on Artificial Intelligence Medicine, London, August 29th-31st 1989. Proceedings, pages 247-256. Springer, 1989.
[8] D. Bertsimas and J. N. Tsitsiklis. scientific Belmont, MA, 1997.
Introduction linear optimization, volume 6. Athena
[9] R. Bhattacharya, T. Nagarajan, D. Malinsky, and I. Shpitser. Differentiable causal discovery under unmeasured confounding. In International Conference on Artificial Intelligence and Statistics, pages 2314-2322, PMLR, 2021.
[10] Binder, D. Koller, S. Russell, and K. Kanazawa. Adaptive probabilistic networks with hidden variables. Machine Learning, 29:213-244 1997.
[11] W. Buntine. Theory refinement on bayesian networks. Uncertainty in artificial intelligence, pages 52-60. Elsevier, 1991.
[12] R. Chen, S. Dash, and Gao. Integer programming for causal structure learning in the pres- ence of latent variables. In International Conference on Machine Learning, pages 1550-1560. PMLR, 2021.
[13] D. M. Chickering, D. Heckerman, and C. Meek. Large-sample learning of bayesian networks NP-hard. Journal of Machine Learning Research, 5(0ct): 2004.
[14] D. Colombo and M. H. Maathuis. Order-independden constraint-based causal structure learn- ing. Journal of Machine Learning Research, 15(1):3741-3782, 2014,
[15] Cussens. Bayesian network learning with cutting planes. In Proceedings of the Twenty- Seventh Conference Uncertainty in Artificial Intelligence, UAI' 11, page 153-160, Arling- ton, Virginia, USA, 2011. AUAI Press.
[16] J. Cussens. Column generation for exact BN learning: Work in progress. In Proc. ECAI-2012 workshop on COmbining COnstraint solving with Mining and L.Earning (CoCoMile 2012). 2012.
11

---

[17] J. Cussens. GOBNILP: Learning Bayesian network structure with integer programming. In International Conference on Probabilistic Graphical Models, pages 605-608. PMLR. 2020.
[18] J. Cussens. franch-price-andee for causal discovery. Conference on Causal Learning and Reasoning, pages 642-661. PMLR, 2023.
[19] G. Desaulniers, Desrosiers, and M. M. Solomon. Column Generation, volume Springer Science & Business Media, 2006.
[20] M. Desrochers, J. Desrosiers, and M. Solomon. A new optimization algorithm for the vehicle routing problem with time windows. Operations Research, 40(2):342-354, 1992.
[21] M. Drton and M. H. Maathuis. Structure learning in graphical modeling, Annual Review of Statistics and Its Application, 4(1):365-393, 2017.
[22] M. El Halabi, G. Orfanides, and T. Hoheisel. Difference of submodular minimization via DC programming. In International Conference on Machine Learning, pages 9172-9201. PMLR,
2023.
[23] ETH Zurich.
LUCAS (LUng
CAncer
Simple
set).
https: http://wwww....tttttt..t inf. ethz. dadaaaaaaaaaa html
[24] S. Fujishige. Polymatroidal dependence structure of a set of random variables. and Control, 39(1):55-72, 1978.
Information
[25] S. Fujishige. Submodular Functions and Optimization, volume 58. Elsevier, 2005.
[26] P. C. Gilmore and R. E. Gomory. A linear programming approach to the cutting-stock problem. Operations Research, (66):449-555 1961.
[27] C. N. Glymour and G. F. Cooper. Computation, Causation, and Discovery. AAAI Press, 1999.
(28) D. Heckerman, D. Geiger, and D. M. Chickering. Learning bayesian networks: The combina- tion of knowledge and statistical data. Machine Learning, 20:197-243, 1995.
[29] C. Heinze-Deml, M. H. Maathuis, and N. Meinshausen. Causal structure learning. Annual Review of Statistics and Its Application, 5(1)3371-39,, 2018.
[30] R. Iyer and Bilmes. Algorithms for approximate minimization of the difference between submodular functions, with applications. Uncertainty in Artificial Intelligence (UAI), 2012.
[31] T. Jaakkola, D. Sontag, A. Globerson, and M. Meila. Learning Bayesian network structure using LP relaxations. In Proceedings of the Thirteenth International Conference on Artifi cial Intelligence and Statistics, pages 358-365. JMLR Workshop and Conference Proceedings, 2010.
[32] C. Joncour, S. Michel, R. Sadykov, D. Sverdlov, and F. Vanderbeck. Column generation based primal heuristics. Electronic Notes in Discrete Mathematics, 36:695-702, 2010.
[33] J. E. Kelley, Jr. The cutting-plant method for solving convex programs. Journal of the society for Industrial and Applied Mathematics, 8(4):703- 1960.
[34] A. K. Kelmans and B. N. Kimelfeld. Multiplicative submodularity a matrix's principal minor function of the set of its rows and some combinatorial applications. Discrete Mathematics, 44(1):113-1) 1983.
[35] M. Koivisto and K. Sood. Exact bayesian structure discovery in bayesian networks. Journal of Machine Learning Research, 5(May):549 573,995-33277 2004.
[36] I. Lovász. Submodular functions and convexity. Mathematical Programming The State of the Art: Bonn 1982, pages 235-257, 1983.
[37] M. F. Lübbecke and J. Desrosiers. Selected topics in column generation. Operations Research, 1307 2005.
12

---

[38] C. Meek. Causal inference and causal explanation with background knowledge. arXiv preprint 102.4:130..4977, 2013.
[39] G. L. Nemhauser and A. Wolsey. Maximizing submodular set functions: formulations and analysis of algorithms. In North-HHolland Mathemalics Studies, volume 59, pages 279-301. Elsevier, 1981.
[40] Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, 2000.
[41] D. Potthoff, D. Huisman, and G. Desaulniers. Column generation with dynamic duty selection for railway crew rescheduling. Transportation Science, 44(4))493-505, 2010.
[42] G. Schwarz. Estimating the dimension of model. The Annals of Statistics, pages 461-464,
1978.
[43] M. Scutari. Learning Bayesian networks with the bnlearn R package. Journal of Statistical Software, 35:1-22, 2010.
[44] T. Silander and P. Myllymaki. A simple approach for finding the globally optimal bayesian network structure. arXiv preprint
2012.
[45]
Spirtes, C. N. Glymour, and R. Scheines. Causation, Prediction, and Search. MIT Press, 2000.
[46] R. Tillman and P. Spirtes. Learning equivalence classes of acyclic models with latent and selec- tion variables from multiple datasets with overlapping variables. In Proceedings of the Four- teenth International Conference on Artificial Intelligence and Statistics, pages 3-15. JMLR Workshop and Conference Proceedings, 2011.
[47] Tsamardinos, L. E. Brown, and C. F. Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine Learning, 65:31-78, 2006.
[48] S. Van de Geer and P. Bühlmann. L0-penalized maximum likelihood for sparse directed acyelic graphs. The Annals of Statistics, 43(2):536-567, 2013.
[49] Y. Yu, J. Chen, T. Gao, and M. Yu. DAG-GNN: DAG structure learning with graph neural networks. In International Conference on Machine Learning, pages 7154-7163. PMLR, 2019.
[50] Zhang. On the completeness of orientation rules for causal discovery in the presence of 2008.
latent confounders and selection bias. Artificial Intelligence,
13

---

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: In the abstract and introduction (Section D. we clearly state the contributions and scope of this 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 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.
It is fine to include aspirational goals motivation as long as it 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 limitations of our work is mainly discussed in Section 5.3 and in the
conclusion (Section
Guidelines:
The answer NA means that the paper has 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 separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The au- thors should reflect on how these assumptions might be violated practice and what the implications would be.
The authors should reflect 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 on the factors that influence the performance of the ap- proach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-tex system might not be used reliably provide closed captions for online lectures because 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 ad- dress problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
14

---

Answer: [Yes]
Justification: Our theoretical results are in Proposition The proofs are in Appendix
C
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 theo- rems.
The proof' can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be comple- mented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data provided not)?
Answer: [Yes]
Justification: The implementation details of our algorithm and the experiment settings are fully discussed, which can be found in Section 3.2 and 5
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer 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 a dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, if the contribution is specific model and empirical evaluation, may be necessary to either make it 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 to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, other means that are appropriate the research performed.
While NeurIPS does not require releasing code, the conference does require all sub- missions to provide some reasonable avenue for reproducibility, which may depend on the nature the contribution. For example
(a) the contribution primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) 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 way to re- produce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case au- thors are welcome to describe 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., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
15

---

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, described in supplemental material?
Answer: [Yes]
We provide the access to the code of method and the code for generating data.
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please
data submission guidelines httpp:///iips.c/puulicccoo http:///eee/deeeeeasaaaiiooooooo
the NeurlPS code and
for more de-
tails.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to to reproduce the results. See the NeurIPS code and data submission guidelines (uttps://nips.. htttp:::/iips.cc/ppbi///////////
for more details. The authors should provide instructions data access and preparation, including how to 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. 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 as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
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?
Question:
Answer: [Yes]
Justification: We specify all the implementation details that is necessary for understanding the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper 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 appropri- ate information about the statistical significance of the experiments?
Answer: [No]
Justification: Our reported metrics for numerical experiments are discussed in Section 5 which does not include error bars or confidence intervals.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
16

---

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, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the enror bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should prefer- ably report 2-sigma error bar than state that they have a 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 out of 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 execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: See Section
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU 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 as 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 failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer:
[Yes]
Justification: Our research conform with the NeurlPS Code of Ethics in every respect.
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 circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: Our work is primarily methodological in nature without direct societal impact.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
17

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g,, deployment of technologies that could make decisions that unfairly impact spe- cific 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 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 to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used 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 unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitiga- tion strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
Safeguards
Question: Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification:
The paper poses no such risks.
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
allow for controlled use of the model, for example by re- quiring that users adhere to usage guidelines restrictions 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 is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or 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: [Yes]
Justification: Properly credited.
Guidelines:
The answer NA means that the paper does 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 used and, 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 (e.g., website), the copyright and terms of service of that source should be provided.
license, copyright information, and terms of use in the pack- age should be provided. For popular datasets, haperswithcoda com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the li- cense of dataset.
If assets are released,
18

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documenta- tion provided alongside the assets?
Answer:
[Yes]
Justification:
We provide details of the code as part of our submissions.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of 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 pa- per include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Our paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contri- bution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, cura- tion, or other labor should be paid at least the minimum wage in the country of the data collector.
15. 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 the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our paper does not involve crowdsoureing nor research with human subjects.
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 equiva- lent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS 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.
19

---

16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the 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.
Question:
Answer: [NA]
Justification: The core method development in this research does not involve LLMs as any important, original or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
/ Please refer to our L.LM policy (https://net) for what should or should be described.
20

---

A Definitions
A.1 The Data Matrix and Scoring Functions
We denote D as the dataset, which is an x matrix where the sample size and n is the number of random variables (nodes Bayesian Network). The entry D.j is the value of j-th variable Xj in the i-th observation. For discrete labeled data, each label is mapped to a distinct numerical value.
We use the eo-penalized liglliiel function as the scoring function of a DAG, which is
score D) :- (og(L(G; - k(G)
for some A 2 0. Here, L(G; D) is the likelihood of the graph structure G under data D, and k(G) the number free parameters to estimated in the graphical model. crucial property of this scoring function is that it can be decomposed into node-specific local scores
acore(G:D)
score,
where each local score score;(pa,()) only depends on a node i and its parent set pa,(G).
For continuous data, assume that variables (Xi mean-normalized such that [[] 0. We also assume that if J is the parent set of node i, then X; X,, [88] where X, = jjee The parameters Oij and o²ا_ are unknown. The likelihood function defining score; has included these parameters by their maximum likelihood estimator, thus its value only depends on PG,(G). The resulting local score score;(J) this context is
N log(d?_)) A|J,,
score; = (1 (og((T)
2
2
where
0 ,,, min Ep [(a x, x.) aER
denotes the empirical residual variance (under the empirical distribution P)) of the linear regression predicting X, (with B[X,] = 0) from the features X.,. For simplicity, we will ignore the constant
2 (1 + log(2r)) in our pricing problem optimization framework established based on Proposi- tion
For discrete data, we assume that if J is the parent set of node i, then X, follows a multinomial dis- tribution with parameters depending the configuration of the parent values X, {Xy]ee. The local score for discrete data using multinomial likelihood
#(m;,x)) #(x))
score:(J) Σ Σ #(x,, x)) log
A(a; 1) II aj.
2,ES, XJES,
JEJ
where S, and S, represent the sets of possible values that X, and X, take, respectively. The arity (i.e., the number of possible values it can take) of variable X, is denoted by a, for = 1,... The count function # (3r,, x)) Count Ti, X, x,) counts joint occurrences in the dataset D, and # (xx,) = Count(X. = x)) provides the corresponding marginal counts.
A.2 Submodular Function and Lovász Extension Function
We adopt standard definitions of the submodular set function and the Lovász extension 36). Definition (Submodular Set Function). Let V be a finite ground set. A set function f f :: + is submodular if it satisfies the diminishing return property, ie., for all A B V and V B,
[(AU(u)) - S(A) 2 ((BU J(B) -
A function is supermodular if -f is submodular, and is modular if it is both submodular and super- modular.
Definition 2 (Lovász Extension). Let f 2V R be function defined ground set V = (l,...,dl d}. Given a point z E 1], let a be a permutation of such that
Ez(1) To(2) Iz(d).
21

---

Define the nested subsets S = {o( for & E {l....d)} and So = 0. The Lovász
extension f 10, 11 7 R of the set function f at point x is defined as
('(x)
(Jo(H) To(k+1)) (S%))
where Io(0) 1 and (o(4+1)
==
B Row Generation and the Separation Problem
The exponential number of cluster constraints can be handled through generation 15 This cutting plane approach starts with restricted set Ć and sequentially adds constraints violated by the current solution x". For an integer solution x", the separation problem identifies violated cluster constraints via cycle detection in the decoded graph from xx (using, e-g., depth-first search). For fractional solution of R.M.P,, the separation problem identifies maximally violated cluster constraints through
max Σ ΣΣ
رirs 1
IEC JnC
which can be formulated as an Separation IP follows:
max M.
Σ Σ
ViLs Vid Σ 2i
(5a)
JEP.
V.J S 2,
i
n1, J E P,,
(5b)
- Σ
i 1,..,, n2 J E P.,
(5c)
Σ 2 1,
(5d)
1=1
V/J (0, 1), E (0,1) i= 1, ,, J E Pi.
(5e)
The optimal solution (y", 2*) of Problem 5 defines the new cluster C* = (i: 27 1 to be added
to C.
The complete row and column generation algorithm iteratively alternates between column genera-
tion for variable selection and row generation for constraint enforcement, dynamically refining both the solution space and constraint set.
C Proofs
C.1 Proof of Proposition I
Proof. Recall that for Gaussian data,
N
N
score; = 2
log(2r +1) log(o? AJJ.
2
The conditional variance o satisfies:
min Epl(a T x, x aERi
min a E. a 2x T + Ex,t
aER
= E,, $ 1 det((Eiui, /det (Ei
22

---

where the last equality is due to the Schur complement, IP denotes the empirical distribution, with E, Ei,, Jus and Σi,; representing the empirical covariance matrices of variable set X,, cross-covariance between X; and x, joint covariance of X, UX and variance of X,, respectively.
The logarithmic transformation yields
))) = log det (uur, log det (E.,),
where both log det({S.,)) and log det det(S..ui,..u) A|J + X is modular in J.
submodular functions of set J
34
Function
The remaining term involving J.E
No. in the pricing objective is also submodular in J.
Ced:
iEC, JnC+0
This follows from the property that for any J1 C J2 and j $ J2.
Σ
X
Σ
X
CEC: iEC, Junc=
CEÈ: iec, JanC=0
Given that X 2 0, the above inequality holds since J1 S J2 implies that C' n J2 0 necessitates CnJ = 0.
Consequently, the pricing objective to be minimized can be expressed as the following DS function:
N
N
NV.
z(J;A*) log Set(( + A|J +
Σ
log det E..J)
log(2x 1
2
CEC: iE,, JrCØ8
submodular
submodular
C.2 Proof of Proposition 2
Proof. Recall that for multinomial data,
# x.)
sJ) Σ Σ
#(x,,x, x,) log
A(a; - 1) II
aj.
#(x))
T,E&; xyE&,
jCJ
We reformulate the first term follows:
#(1,,,)) #(x))
Σ Σ x,),,,))l log
8,E8, xyES,
Σ
Σ
#a,,, (x,,,,) log( (I,,,
Σ
((xx))lo((*(X))
x, E,ES; ES
x,ES,
=N Σ L.ES XJES,
Σ
P(x x; x,, = xx log P(x; = Ri, X, x==x,)
Σ
pxxX, = x1)
x,ES,
=N(-H(JU + H(J)),
where N is the sample size. The functions
H(J)
= - Σ x log log = x,)
XJESS
and
X, =x, Fi, X, x,)
= Σ Σ -
X,ES,
denote the entropy functions, which are known to be submodular (44]
By the proof Proposition D Σ cEd: iEC, JnC40
N is submodular J. is also easy to verify that
(a: - 1) Jje aj is supermodular in J as aj 2 1 for j E J (1).
23

---

Therefore, the pricing objective can be expressed as the following DS function:
=(J;X)= N. +(J Σ +
N H()H()) A(a; - -1)
1) 4g
C.CEC iEC ChJ
JEJ
submodular
submedular
D Supplementary Numerical Results
D.1
Comparison DCA initialization methods
To evaluate the sensitivity of DCA initialization, we compare the BIC scores and time costs across three initialization methods (warm-start random, and hybrid). The results are summarized in Table
and S
Table 4: BIC Score Comparison of the Three Initialization Approaches for DCA
(m, N; d)
hybrid
rundken
wantssart
(15. 5000,0.5)
-101818.64
-10181864
101819.55
(15.
5000.1.0
-101954.63
-1019399775
1019975.70
2000.1.5)
101411.14
-101376.34
1012099599
(15,
101164622
-101556.54
(15, 20000.0.5)
-06645.22
-06345.225
2063346.2
(15, 20000.10)
410415 68
410415.68
410082.11
2.500.1.5)
401534424
-05487.11
4015999
(15, 20000,2.0)
501710:97
(20L 2000,0.5)
-136605.25
136605.25
136734.74
(20) 5000.1.0
-137222.15 13
1372112.7
-157363 27
-136624.8
-36677..8
(20, 2.00
366449
1336610.4
196750-04 44
(20, 20000.0.55
535504792
5350147.92
555001 201
(201 100
541859016
-41887.02
54226N 27
2.000.1.5)
S3869399
(20) 20000.200)
53411587
(25, 2000,0.5)
.170751.53
1770751.53
170611-74
(25, 5.00.100
-1702799.9
1100002.0 01
-171280.
(25, 5000,2.0)
-1713333 89
-171105
-170760.00
(25, 20000.0.5)
-677811 87
677815.87
-677812.65
(25, 20000.1.0)
-683475 44
6831899
-688866 26
1.5)
676267.45
-676649.
(25, 30000,2.00
-689041.3
408690.44
-687147.17
'Bach value the average ever indapendent instances.
Table 5: Runtime Comparison of the Three Initialization Approaches for DCA
(n, .N. d)
hybrid (s)
random (s)
wamnssart (s)
(15. 5000, 0.5)
4,06
4.04
1.17
(15. 5000, 1.0)
25.76
29: 40
4.49
$7.45
10.00
(15. 20000 0.5)
10.17
10.94
3,51
(15. 2000, 1.0)
36.71
36.57
10.63
2000,
106.25
109.73
23.95
2000,
229.80
(20. 5000, 0.5)
9.53
9.60
2.42
(20. 5000, 1.0)
84.34
90.87
12.51
1770.71
(20, 20000. 0.5)
18:42
18.54
8.71
(20, 20000, 1.0))
133.09
134.79
34,84
20000,
425.17
504.29
441.37
20000
1710.91
(25. 5000. 0.5)
21.97
22.35
3.87
(25. 5000, 1.0)
261,34
301.12
33.59
2434.05
131.69
(25.
6037.58
(25. 20000, 0.5)
50.22
53.15
12.48
(25, 20000. 1.4))
404..33
642.87
78.21
1372.59
3088.33
227.15
(25, 20000. 2.00
366.46
Tach value Uhe average over tum independent instomces.
From the tables, we observe that

---

Random initialization yields highest (best among three) BIC scores in 15 out of 24 instances among three methods but requires the most time in 22 out 24 instances.
Warm-start initialization achieves the lowest (worst among three) BIC scores in 16 out of 24 instances and is the most time-efficient in 23 out of 24 instances.
The hybrid approach strikes a balance between optimality and computational efficiency, since it focuses on local refinement around the current best pattern (through warm-start initialization) while building upon the foundation of global exploration (through random initialization) at the early stage.
25

---

D.2 Comparison of CG-DCA with Baselines
Table presents the comprehensive experimental results including the HC method, with N (5000, 20000} for Gaussian datasets. For the three discrete datasets, Table / provides the score and time comparison for CG-DCA and GOBNILP, while Provids the graph comparisons of all baselines. The runtime performance of all baseline methods (implemented in R) on all datasets summarized in Table 9 and
10
Table 6: Performance of CG-DCA, GOBNILP (GI), HC, (Stable-)PC, and MMHC on Gaussian Datasets
a.
Procision
Recall
CG-DCA
115. 5000, 0.5)
0.78
0.66
0.36
0.50
0.50
0.80
0.73
0.47
0.62
0.51
1.70
2.90
6.20
3.00
115. 5000. 1.0)
0.70
0.67
0.41
0.72
0.59
0.78
0.72
0.52
0.67
0.54
6.20
5.00
11.90
4.50
6.70
(15.
5000 1.5)
0.65
0.73
0.42
0.61
0.56
0.77
0.80
0.57
01.44
0.48
11.80
7.40
19.50
11.50
12:30
5000.
0.55
0.34
0.54
0.52
0.39
21.90
16.90
35,30
20.90
20.40
115. 20000 1.0)
0.76
0.78
0.49
0.71
0.61
0.78
0.81
0.61
0.69
0.59
4.10
4.00
9.50
4.70
6.20
(15, 20000 1.5)
0.76
0.73
0.47
0.72
0.66
0.85
0.79
0.61
0.66
0.57
6.90
9.40
18.30
5.60
10.00
115, 2000. 2.0)
0.61
0.53
0.38
0.5.
0.55
0.75
0.60
0.62
0.43
0.39
20.00
19.90
33.20
20.01
19,50
120. 5000,
0,74
0.71
0.45
0.76
0.59
0.80
0.75
0.56
0.79
3.40
3,50
6.80
2.60
4.40
(20) 5000, 1.5)
0.50
0.61
0.41
0.62
0.65
0.50
28.
17.90
16.80
120. 5000, 2.00
0.44
0.46
0.41
D:54
0.55
0.66
0.53
0.67
0.54
0.31
36 40
33.60
49,80
31.50
31.10
120, 20000.
0.82
0.80
0.54
0.73
0.58
0.83
0.82
0.62
0.76
0.59
2.20
2.60
5.70
3.00
4.60
120. 20000.
0.61
0.73
0.56
0.79
0.68
0.70
0.79
0.67
0.77
0.65
8.30
6,90
10,80
4,90
7.10
16.50
0.48
29.30
(25. 5000, 0.5)
0.82
0.61
0.60
0.84
0.69
0.67
0.67
0.84
0.69
2.70
3.90
6.70
2.70
125. 5000, 1.00
0.81
0.69
0.57
0.80
0.74
0.86
0.76
0.71
0.75
0.69
6.70
11.10
16.30
7.40
8.70
125. 5000)
0.66
0.57
0.44
0.70
0.71
0.75
0.68
0.65
0.58
0.56
13,60
22,10
59.50
18.50
18.00
5000,
0.54
0.38
0.53
0.35
45.70
38.70
36.50
7.30
0.76
0.71
9.90
125, 20000, 1.5)
0.57
0.53
0.40
0.73
0.63
0.67
0.66
0.62
0.63
0.51
14.10
25.00
41.90
16.50
19.50
(25,
20000 2.0)
0.54
0.48
0.35
0.61
0.59
0.66
0.56
0.64
0.46
0.40
23.90
39,10
77.00
33.80
34.10
Eaah volua is the average over 160 indppndentt instances.
Table 7: Score and Time Comparison of CG-DCA and GOBNILP on Discrete Datasets
BIC
Gap
(scconds)
CG-DCA
GOBNILP
CG-DCA
GOBNILP
LUCAS ALARM
0.03
0.00
11.59
7,99
2,85
0.00
4696.19
153,49
Table 8: Performance Comparison of CG-DCA, GOBNILP (GD) HC, PC, and MMHC Discrete Datasets
Datnee
CG-DCA
MMTIC
CG-DCA
MMIIC
CG-DCA
PC
MMII
LUCAS
0.77
1.00
0.83
0.33
0.83
1.00
0.33
0.83
0.33
10
INSURANCE
37
22.
ALARM
0.46
0.87
0.42
0.83
0.47
0.50
0.89
0.48
0.76
0.53
35
7
34
11
32
26

---

Table 9: Time Comparison of Stable PC, HC, and MMHC Algorithms on Gaussian Datasets
(n. N; (4)
Stalle (C
HC
MMHC
115.
0.148
115. 5000, 1.00
0.035
0.054
0.047
[15. 5000, 1.5)
0.084
0.101
0.062
0.092
[15. 2000, 1.0)
0.121
0.196
0.152
115. 20000 1.5)
0.188
0.203
0.326
115.
0.306
0.889
(20.
0.090
120. 5000, 1.00)
0.105
0.159
0.112
(20. 5000, 1.5)
0,180
0.309
0.291
0.266
120,
0.185
(20, 2000,, 1.0)
0.279
0,481
0.438
(20, 20000, 1.5)
0.660
0,950
1,366
120,
1.327
3.634
(25.
0.083
[25. 5000, 1.0))
0,955
0.133
0.09%
(25. 5000, 1.5)
0.227
0,348
0.218
(25,
[25. 2000, 1.0)
0.288
0.637
0,583
(25, 20000
0.725
1.568
2.612
2.940
"Lach value 35 the average over ten independent instances
Table 10: Time Comparison of Stable PC, HC, and MMHC Algorithms on Discrete Datasets
Datnect
Stale PC
HC
MMHC (4)
LLCAS
0.08
0.088
0.0N
ALARM
0.61
0.25
0.53
27


Reviews:
[{"review_id": "f5uA8NAsUZ", "reviewer": "Reviewer_oPou", "summary": "This paper addresses the computational bottleneck in Bayesian Network Structure Learning (BNSL) via Integer Programming (IP), where traditional column generation (CG) suffers from expensive pricing problems solved via Mixed-Integer Nonlinear Programming (MINLP). The authors propose CG-DCA: a CG framework that reformulates pricing as a Difference-of-Submodular (DS) optimization problem, solved approximately via the Difference-of-Convex Algorithm (DCA).", "questions": "1. Could you propose/explore alternative DS decompositions or scoring functions to handle the super-modular penalty in discrete BNSL? Would DCA with tailored sub-gradients help? \n\n2. How sensitive is DCA to initialization? Does the hybrid strategy consistently outperform pure warm-start/random?", "limitations": "yes", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths: \n1. Novelty: First to reformulate BNSL pricing as DS optimization, creatively leveraging DCA’s efficiency. Theoretically grounded (Propositions 1–2).  \n2. Practical Impact: Achieves 10–100× speedup vs. MINLP pricing (0.3s vs. 35.1s per node) and higher solution quality than GOBNILP on Gaussian data (Table 1). Handles denser graphs d=2 better than exact IP. \n3. Clarity:  Strong exposition of BNSL IP formulation (Section 2.1), CG pitfalls, and DS optimization. \n\nWeaknesses: \n1. Limited Scope:  Effectiveness is restricted to Gaussian data. Discrete-data performance degrades severely (Table 3), acknowledged but not mitigated.  \n2. Scalability Gaps: Runtime remains high for n=25,d=2 (4,850s), and no results for n>37. DCA’s convergence in pricing is empirical but untheorized.  \n3. Baseline Comparison:** Lacks runtime comparison against non-IP baselines (PC, MMHC). Some metrics (SHD) show PC/MMHC competitive or better (Table 2).", "was_revised": true, "final_justification": "All my concerns have been addressed.", "rebuttal": "We thank the reviewer for their time and insightful feedback. Our responses to the comments are provided below.\n\n**1. Limitations in Discrete Data Implementation**\\\nWe acknowledge the limitations of our current implementation for discrete data. Regarding the decomposition of the discrete BIC score, the approach we employ represents the most effective strategy we have identified so far, though we recognize that alternative decompositions may exist.\\\nTo address the supermodular penalty issue in discrete BNSL, we are currently investigating a strategy that involves convexifying the cardinality penalty term in the BIC score. To be specific, the pricing problem objective for node $i$ can be decomposed as $F^i(J) = G^i(J) - H^i(J) + L^i(J)$, where both $G^i(J)$ and $H^i(J)$ are submodular, and the supermodular penalty term $L^i(J)=\\Lambda(a_i-1) \\prod_{j\\in J} a_j$ depends on the number of possible values $a_j$ for each discrete variable $X_j$ , $j=1,\\ldots,n$. In our difference-of-submodular formulation, $G^i(J)$ represents the first submodular function, while $H^i(J) - L^i(J)$ constitutes the second. Unlike previous implementations that iteratively linearize the whole supermodular component $(-H^i(J)+ L^i(J))$ in DCA, leading to higher approximation errors, the new approach currently under investigation convexifies the penalty term $L^i(J)$ through cutting planes while maintaining $G^i(J) - H^i(J)$ in the DC program. Assuming binary data (such that $L^i(J)=\\Lambda 2^{|J|}$), we enhance the optimization program by adding cutting planes that represent the convex envelope of the supermodular penalty term $L^i(J)$. We have a preliminary experiment result on LUCAS data: \\\nPrecision: 0.85, Recall: 0.92, SHD: 2,\\\nand may improve and implement the approach to more general settings.\\\nAdditionally, we have developed a decomposition strategy for the BDeu score as a difference of submodular set functions. This approach may have potential to help mitigate the supermodularity-related difficulties, due to the structural differences between the BDeu and discrete BIC scores. While implementation is still underway, we plan to explore this in future work and will leave some relevant ideas in the conclusion section of this manuscript.\\\nRegarding the tailored subgradient method, we are not sure what you are referring to and have not yet incorporated it into our framework. However, we would greatly appreciate any references or additional information you could provide to guide its potential implementation.\n\n**2. Runtime Comparison with Baseline Methods**\\\nThe runtime performance of all baseline methods (implemented in R) is summarized in the table below. This table will be included in Appendix.\n| n, N, d      | stable PC | HC     | MMHC   |\n|---------|-----------|--------|--------|\n| 15, 5000, 0.5 | 0.038s | 0.048s | 0.048s |\n| 15, 5000, 1 | 0.035s | 0.054s | 0.047s |\n| 15, 5000, 1.5 | 0.054s | 0.101s | 0.062s |\n| 15, 5000, 2 | 0.095s | 0.098s | 0.092s |\n| 15, 20000, 0.5 | 0.097s | 0.154s | 0.126s |\n| 15, 20000, 1 | 0.121s | 0.196s | 0.152s |\n| 15, 20000, 1.5 | 0.188s | 0.263s | 0.326s |\n| 15, 20000, 2 | 0.306s | 0.420s | 0.889s |\n| 20, 5000, 0.5 | 0.050s | 0.066s | 0.090s |\n| 20, 5000, 1 | 0.105s | 0.159s | 0.112s |\n| 20, 5000, 1.5 | 0.180s | 0.309s | 0.291s |\n| 20, 5000, 2 | 0.266s | 0.342s | 0.425s |\n| 20, 20000, 0.5 | 0.146s | 0.255s | 0.185s |\n| 20, 20000, 1 | 0.279s | 0.481s | 0.438s |\n| 20, 20000, 1.5 | 0.660s | 0.950s | 1.366s |\n| 20, 20000, 2 | 1.327s | 2.331s | 3.634s |\n| 25, 5000, 0.5 | 0.083s | 0.099s | 0.071s |\n| 25, 5000, 1 | 0.095s | 0.133s | 0.098s |\n| 25, 5000, 1.5 | 0.227s | 0.348s | 0.218s |\n| 25, 5000, 2 | 0.273s | 0.403s | 0.351s |\n| 25, 20000, 0.5 | 0.200s | 0.407s | 0.235s |\n| 25, 20000, 1 | 0.288s | 0.637s | 0.583s |\n| 25, 20000, 1.5 | 0.725s | 1.568s | 2.612s |\n| 25, 20000, 2 | 1.158s | 2.087s | 2.940s |\n\nWe acknowledge that certain constraint-based methods exhibit significantly faster execution times. Compared to constraint-based method that operate sequentially and are susceptible to error propagation [1], exact score-based approaches consider the problem globally, though at potentially greater computational cost. Compared with the Hill Climbing (HC) method that aims at heuristically solving the problem at a low cost, the score-based IP methods transform the original problem complexity into that of solving the optimization problem, which our method specifically aim to addresses. Furthermore, our ultimate goal is to develop a truly exact score-based method, where our current framework serves as a subroutine to efficiently generate high-quality columns. Please refer to our response to Reviewer krir for some additional discussions.\n\n**3. Sensitivity to DCA Initialization**\\\nTo evaluate the sensitivity of DCA initialization, we compare the BIC scores and time costs across three initialization methods (warm-start, random, and hybrid). The results are summarized in the following tables:\n| n, N, d      | hybrid       | random       | warmstart    |\n|--------------|--------------|--------------|--------------|\n| 15, 5000, 0.5 | -101818.64   | -101818.64   | -101819.58   |\n| 15, 5000, 1   | -101954.63   | -101939.75   | -101975.70   |\n| 15, 5000, 1.5 | -101411.14   | -101376.34   | -101509.59   |\n| 15, 5000, 2   | -101646.26   | -101556.54   | -101608.25   |\n| 15, 20000, 0.5 | -406345.25   | -406345.25   | -406346.23   |\n| 15, 20000, 1  | -410415.68   | -410415.68   | -410382.11   |\n| 15, 20000, 1.5 | -405344.24   | -405487.19   | -405599.33   |\n| 15, 20000, 2  | -400909.42   | -400861.51   | -401710.97   |\n| 20, 5000, 0.5 | -136605.28   | -136605.28   | -136734.74   |\n| 20, 5000, 1   | -137222.13   | -137212.97   | -137363.27   |\n| 20, 5000, 1.5 | -136624.81   | -136577.03   | -136957.87   |\n| 20, 5000, 2   | -136849.69   | -136610.04   | -136750.44   |\n| 20, 20000, 0.5 | -535047.92   | -535047.92   | -535183.20   |\n| 20, 20000, 1  | -541850.16   | -541887.02   | -542168.27   |\n| 20, 20000, 1.5 | -538682.19   | -538693.99   | -539863.90   |\n| 20, 20000, 2  | -535301.38   | -534115.87   | -534939.03   |\n| 25, 5000, 0.5 | -170751.53   | -170751.53   | -170751.74   |\n| 25, 5000, 1   | -170279.79   | -170296.43   | -170302.01   |\n| 25, 5000, 1.5 | -171308.77   | -171157.64   | -171289.38   |\n| 25, 5000, 2   | -171333.89   | -171105.48   | -170760.08  |\n| 25, 20000, 0.5 | -677813.87   | -677813.87   | -677812.65   |\n| 25, 20000, 1  | -683475.44   | -683189.30   | -683826.26   |\n| 25, 20000, 1.5 | -676190.15   | -676267.45   | -676649.76   |\n| 25, 20000, 2  | -689041.33   | -688690.44   | -687147.17  |\n| **Average**   | **-312,582.26** | **-312,553.87** | **-312,712.62** |\n\n| n, N, d      | hybrid    | random    | warmstart |\n|--------------|-----------|-----------|-----------|\n| 15, 5000, 0.5 | 4.06s     | 4.04s     | 1.17s     |\n| 15, 5000, 1   | 25.76s    | 29.40s    | 4.49s     |\n| 15, 5000, 1.5 | 76.72s    | 87.43s    | 10.00s    |\n| 15, 5000, 2   | 170.39s   | 295.25s   | 23.74s    |\n| 15, 20000, 0.5 | 10.17s    | 10.94s    | 3.51s     |\n| 15, 20000, 1  | 36.71s    | 36.57s    | 10.63s    |\n| 15, 20000, 1.5 | 104.25s   | 109.73s   | 23.95s    |\n| 15, 20000, 2  | 229.80s   | 245.11s   | 41.94s    |\n| 20, 5000, 0.5 | 9.53s     | 9.60s     | 2.42s     |\n| 20, 5000, 1   | 84.34s    | 90.87s    | 12.51s    |\n| 20, 5000, 1.5 | 312.67s   | 430.37s   | 48.78s    |\n| 20, 5000, 2   | 1039.53s  | 1770.71s  | 95.36s    |\n| 20, 20000, 0.5 | 18.42s    | 18.54s    | 8.71s     |\n| 20, 20000, 1  | 133.69s   | 134.79s   | 34.84s    |\n| 20, 20000, 1.5 | 425.17s   | 504.29s   | 444.37s   |\n| 20, 20000, 2  | 997.99s   | 1710.91s  | 211.44s   |\n| 25, 5000, 0.5 | 21.97s    | 22.35s    | 3.87s     |\n| 25, 5000, 1   | 261.84s   | 304.12s   | 33.59s    |\n| 25, 5000, 1.5 | 949.90s   | 2434.03s  | 131.69s   |\n| 25, 5000, 2   | 3054.63s  | 6037.58s  | 224.73s  |\n| 25, 20000, 0.5 | 50.22s    | 53.13s    | 12.48s    |\n| 25, 20000, 1  | 404.43s   | 642.87s   | 78.21s    |\n| 25, 20000, 1.5 | 1372.59s  | 3688.33s  | 227.15s   |\n| 25, 20000, 2  | 3189.98s  | 6674.62s  | 366.46s  |\n| **Average**   | **495.62s** | **897.93s** | **88.71s** |\n\nFrom the table, we can see that:\n\n- Random initialization yields highest (best among three)\nBIC scores in 11 out of 24 instances among three methods but requires the most time in 22 out of 24 instances.\n\n- Warm-start initialization achieves the lowest (worst among three) BIC scores in 16 out of 24 instances and is the most time-efficient in 23 out of 24 instances.\n\n- The hybrid approach strikes a balance between optimality and computational efficiency, since it focuses on local refinement around the current best pattern (through warm-start initialization) while building upon the foundation of global exploration (through random initialization) at the early stage.\n\nWe appreciate the opportunity to clarify these points and remain open to further discussion or suggestions.\n\n## References\n[1] Federica Cugnata, Ron S Kenett, and Silvia Salini. Bayesian networks in survey data: Robustness and sensitivity issues. *Journal of Quality Technology*, 48(3):253–264, 2016."}, {"review_id": "3pZ8J2RwLw", "reviewer": "Reviewer_krir", "summary": "The paper describes a column generation integer-programming approach for Bayesian network structure (graph) learning. The novelty is an inexact pricing method based on algorithms for difference-of-convex-function minimization; the rest of the overall framework is a recent exact method for the same problem, with exact but very costly pricing subproblems. Embedding the inexact pricing routine proposed in the paper into the column generation framework is demonstrated to provide high-quality results in typically notably shorter time than the exact approach, and the overall procedure is competitive with several other (inexact) algorithms that have been proposed for this graph learning problem over the past decade or so.", "questions": "- Tables 2 and 3 (and 4 and 5 in the appendix) should be complemented by the associated running times of the different methods. The current presentation does not give any information about this, but it constitutes an important aspect when discussing solution algorithms (even, maybe especially, if talking about inexact methods). \n- Have the authors tested a version of their method that resorts to (costly) exact pricing if the novel heuristic pricing fails to provide a new column? That would turn the overall approach into an exact algorithm, which would still exploit the faster heuristic pricing whenever possible. It would be interesting to see how this compares to the inexact and the older exact method.\n- A few points could be made a bit clearer: Could you clarify what exactly is being \"observed\", i.e., what the \"given observational data $D$\" entails? Could you state clearly how you implemented your code (is it written in R?), and did you implement the other methods yourself or use codes provided by the respective authors? (For fairness of runtime comparisons: Are they all coded in the same language, or is it a mix?) What memory and CPU speed does the employed computer have? Finally, regarding reproducibility of the experiments, will a link to the code and data be provided upon acceptance? (Otherwise, the results are not reproducible exactly, due to missing information regarding programming language, random seeds used, etc.)\n\n* edit: I just saw the supplementary material file; somehow missed that initially. It contains the data and exact and inexact column generation approach coded in Python (and will presumably be made publicly available?), so that partly answers some of the questions in my last bullet point.", "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 and clearly explains the column generation approach, which is technically not the main contribution of the paper, but still commendable, especially given that most of the machine learning community is not overly familiar with the technical details of integer programming solution techniques and algorithms. The novel inexact pricing scheme is well-motived, described clearly, and then demonstrated in numerical experiments to yield an overall algorithm that is at least competitive with several others. Thus, the paper is a step forward to improved exact algorithms rooted in integer programming / column generation for the Bayesian network structure learning problem. The technical proofs also appear to be correct, as far as I could discern. The only thing I was missing is a runtime comparison -- this is only done in Table 1, but not Tables 2 and 3, i.e., only for the new method and the related exact approach, but not when comparing against other algorithms (GI, PC, MMHC).", "was_revised": true, "final_justification": "See comments below; overall rating retained at 5.", "rebuttal": "We sincerely appreciate your valuable feedback and insightful suggestions. Below, we address your comments in detail.\n\n**1. Runtime Comparison with Baseline Methods**\\\nThe runtime performance of all baseline methods is summarized in the table below. This table will be included in Appendix.\n| n, N, d      | stable PC | HC     | MMHC   |\n|---------|-----------|--------|--------|\n| 15, 5000, 0.5 | 0.038s | 0.048s | 0.048s |\n| 15, 5000, 1 | 0.035s | 0.054s | 0.047s |\n| 15, 5000, 1.5 | 0.054s | 0.101s | 0.062s |\n| 15, 5000, 2 | 0.095s | 0.098s | 0.092s |\n| 15, 20000, 0.5 | 0.097s | 0.154s | 0.126s |\n| 15, 20000, 1 | 0.121s | 0.196s | 0.152s |\n| 15, 20000, 1.5 | 0.188s | 0.263s | 0.326s |\n| 15, 20000, 2 | 0.306s | 0.420s | 0.889s |\n| 20, 5000, 0.5 | 0.050s | 0.066s | 0.090s |\n| 20, 5000, 1 | 0.105s | 0.159s | 0.112s |\n| 20, 5000, 1.5 | 0.180s | 0.309s | 0.291s |\n| 20, 5000, 2 | 0.266s | 0.342s | 0.425s |\n| 20, 20000, 0.5 | 0.146s | 0.255s | 0.185s |\n| 20, 20000, 1 | 0.279s | 0.481s | 0.438s |\n| 20, 20000, 1.5 | 0.660s | 0.950s | 1.366s |\n| 20, 20000, 2 | 1.327s | 2.331s | 3.634s |\n| 25, 5000, 0.5 | 0.083s | 0.099s | 0.071s |\n| 25, 5000, 1 | 0.095s | 0.133s | 0.098s |\n| 25, 5000, 1.5 | 0.227s | 0.348s | 0.218s |\n| 25, 5000, 2 | 0.273s | 0.403s | 0.351s |\n| 25, 20000, 0.5 | 0.200s | 0.407s | 0.235s |\n| 25, 20000, 1 | 0.288s | 0.637s | 0.583s |\n| 25, 20000, 1.5 | 0.725s | 1.568s | 2.612s |\n| 25, 20000, 2 | 1.158s | 2.087s | 2.940s |\n\nWe acknowledge that certain constraint-based methods exhibit significantly faster execution times. Compared to constraint-based method that operate sequentially and are susceptible to error propagation [1], our score-based approach considers the problem globally, though at potentially greater computational cost. Compared with the Hill Climbing (HC) method that aims at heuristically solving the problem at a low cost, the score-based IP methods transform the original problem complexity into that of solving the optimization problem, which our method specifically aim to addresses. Furthermore, our ultimate goal is to develop a truly exact score-based method, where our current framework serves as a subroutine to efficiently generate high-quality columns.\n\n**2. Exact Pricing Method**\\\nWe greatly appreciate your suggestion regarding the use of exact pricing. In response, we are currently developing an enhanced solution approach by incorporating the cutting plane method proposed in [2]. Unlike solving the pricing MINLP directly with Gurobi—which introduces numerical errors due to logarithmic term approximations—our exact solution avoids this issue.  This exact pricing method is activated when the heuristic algorithm fails to identify new columns. Below, we present some preliminary results obtained by integrating exact pricing into our framework (still using the restricted master heuristic though):\n| n, N, d   | Average BIC score (exact) | Average BIC score (heuristic) | Average Time (exact) | Average Time (heuristic) |\n|---------------------------------|------------------------|----------------------------|---------------------------|-------------------------------|\n| 15, 5000, 0.5                |       -101818.64             | -101818.64                 | 1738.09s                   | 4.06s                          |\n| 15, 5000, 1                  |           -101939.28             | -101954.63                 | 6923.13s                   | 25.76s                         |\n| 15, 5000, 1.5                  |            -101350.15             | -101411.14                 | 11597.61s                  | 76.72s                         |\n| 15, 5000, 2                   |         -101414.39             | -101646.26                 | 11903.48s                  | 170.39s                        |\n\nWhile these initial results demonstrate the effectiveness of the exact method in improving the objective value, we observe that the solution time scales with problem complexity, which motivates ongoing research into more efficient implementations. Our current efforts focus on optimizing the exact pricing framework to maintain its advantages while improving scalability for larger instances.\n\n**3. Clarification on Observational Data**\\\nWe apologize for any confusion regarding the observational data notation. To clarify, the dataset $D$ \nis  $N \\times n$\n data matrix where:\\\n    - $N$ : number of data samples\\\n    - $n$ : number of random variables (nodes in the Bayesian network)\\\n    - Each row $x_i = (x_{i,1},...,x_{i,n})$ represents a complete data sample\\\n    - Each column $D_{:,j} = (x_{1,j},...,x_{N,j})^\\top$ represents all observations for the $j$-th variable\\\n    - $x_{i,j}$ = value of the $j$-th variable in the $i$-th observation. For discrete labeled data, each label is mapped to a distinct numerical value.\\\n    We hope this clarification resolves any ambiguity, and we will clarify it in the paper. \n\n Thank you again for your constructive feedback, which helps us improve our work. We remain open to any further questions or suggestions.\n\n\n\n\n## References\n[1] Federica Cugnata, Ron S Kenett, and Silvia Salini. Bayesian networks in survey data: Robustness and sensitivity issues. *Journal of Quality Technology*, 48(3):253–264, 2016.\\\n[2] George L Nemhauser and Laurence A Wolsey. Maximizing submodular set functions: formulations and analysis of algorithms. In *North-Holland Mathematics Studies*, volume 59, pages 279–301. Elsevier, 1981."}, {"review_id": "ZU4MBoKBtp", "reviewer": "Reviewer_6rfE", "summary": "The paper considers the Bayesian network structure learning (BNSL) problem and adopts an existing score-based integer programming (IP) formulation with exponentially many variables and constraints. Row and column generation are standard techniques for handling large IPs. Existing exact score-based IP methods [16, 18] employ these techniques, but solving the pricing problem for column generation (CG) exactly is expensive, which limits their scalability.  This paper builds on [18], but instead of solving the pricing problem exactly, the authors reformulate it as a difference of submodular (DS) function minimization problem for $\\ell_0$-penalized likelihood scores, then apply the Difference of Convex Algorithm (DCA) which  was shown to be an effective heuristic in prior work. Empirical results comparing their method to various baselines on continuous and discrete data are provided.", "questions": "- Why are baselines not all tested on the same datasets and with all the performance metrics considered (BIC score gap, time, Precision, Recall and SHD)?\n- In Section 5.1, it's not clear how are larger instances generated? \n- Which implementation of GOBNILP do you use, the python one? It's good to specify that.\n- How do you break ties when computing the subgradient  of the Lovász extension $f^L$ on line 4 of Algorithm 1? Prior work showed that performance of DCA [22] / SubSup [30] is sensitive to this choice.\n- In [22], a variant of DCA which checks for local minimality at convergence and restarts the algorithm from the best neighboring point if  not satisfied (Algorithm 2 therein) is shown to converge to a local minimum. This could improve performance. Did you try that?\n- Can you clarify why the supermodular penalization term in the discrete BIC score makes DCA less effective?\n- Include a discussion of the computational complexity of the overall proposed algorithm (Algorithm 2)\n\n\n\nOther suggestions:\n- I suggest trying to initialize DCA from the solution of the pricing problem from the previous iteration in the for loop on line 5. Also initializing with zero at the first iteration might help avoid the dense solutions observed empirically.\n- Provide error bars in results tables \n- Discuss pros & cons of using score-based IP approach over other approaches to better motivate the paper\n- Provide better review of DS minimization in the main text, e.g., provide references to stated results (some are missing such as the form of the subgradient), explain how to obtain an integer solution.\n- In Section 5, it's helpful to clarify what type of methods are the baselines compared with (score-based, constrained-base, hybrid, exact or heuristics, IP, etc).", "limitations": "The authors discuss clearly the limitation of their method on discrete data. \nBut they do not discuss the societal impact of their work, nor do they explain why their work has no societal impact (as requested in the checklist).", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n- The approach of using DS minimization for solving the pricing problem is novel\n- The theoretical results formulating the pricing problem as a DS minimization problem for $\\ell_0$-penalized Gaussian and multinomial likelihood scores are correct\n- The paper is written clearly\n- The proposed method seems to outperform score-based baselines (GOBNILP, CG-MINLP, HC) on synthetic Gaussian data and somewhat comparable to constraint-based (stable-PC) and hybrid (MMHC) baselines.\n- Code is provided in supplementary material\n\nWeaknesses:\n- Experimental results are not well supported: different baselines are compared on different datasets and with different metrics, and some experimental details are missing (see questions for details).\n- The proposed method is computationally expensive, and does not work well on discrete datasets\n- The proposed method performs worse than the constraint-based (stable-PC) and hybrid (MMHC) baselines in some settings", "was_revised": true, "final_justification": "The authors addressed most of my concerns, and ones raised by other reviewers. The proposed method still have important limitations (high computational complexity, degraded performance on discrete datasets) but is still a good contribution towards better scoring-based methods for the BNSL problem.", "rebuttal": "We sincerely appreciate the reviewer's detailed feedback and thoughtful suggestions. Below we provide the responses.\n\n**1. Experimental Results Presentation**\\\nRegarding the experimental results presentation, we acknowledge the confusion and apologize for any lack of clarity. All baseline methods were evaluated on identical datasets using consistent performance metrics, with the exception of runtime comparisons for non-IP-based methods due to their R implementation versus our Python framework. Tables 1 and 2 in our paper share the same datasets, while Tables 2 and 3 employ identical metrics. For non-score-based methods (stable-PC and MMHC), we omitted score reporting as it wouldn't provide meaningful comparison. Results for $N=20,000$ samples are provided in the appendix, and we will additionally include the previously omitted table comparing scores and runtime for $N=20,000$ continuous and discrete data to ensure comprehensive reporting.\n\n**2. Computational Cost**\\\nWe fully recognize that our proposed approach incurs higher computational costs compared to non-IP methods. However, our primary objective is to develop a truly exact score-based framework where the current implementation serves as an efficient subroutine for column generation. In the meanwhile, we have achieved some efficiency improvement though some implementation enhancements, including:\n- Subgradient evaluation through Cholesky decomposition (for the Gaussian case). Indeed we only need to apply Cholesky decomposition once to a permuted sample covariance matrix to obtain each subgradient of the Lovász extension (rather than evaluating the function value O(n) times).\n- Warmstart initialization for the IP master problem. We initialize the master IP with the solution obtained from the last iteration to quickly drive up the primal bound.\n\nThe updated results demonstrating these enhancements are shown in the following table and will be incorporated in the revised manuscript.\n| n, N, d       | BIC score (previous) | BIC score (updated) | CG time (previous) | CG times (updated) |\n|-----------------|--------------------------|-------------------------|--------------------------|-------------------------|\n| 15, 5000, 1   | -101955.22               | -101954.63              | 24.92s                   | 25.76s                  |\n| 15, 5000, 2   | -101516.72               | -101646.26              | 217.87s                  | 170.39s                 |\n| 15, 20000, 1  | -410315.40               | -410415.68              | 46.32s                   | 36.71s                  |\n| 15, 20000, 2  | -400946.36               | -400909.42              | 284.91s                  | 229.80s                 |\n| 25, 5000, 1  | -170335.63               | -170279.79              | 272.19s                  | 261.84s                 |\n| 25, 5000, 2  | -171048.00               | -171333.89              | 4850.47s                 | 3054.63s                |\n| 25, 20000, 1  | -683322.64               | -683475.44              | 650.32s                  | 404.43s                 |\n| 25, 20000, 2  | -688321.14               | -689041.33              | 4162.49s                 | 3189.98s                |\n|Average  |-312,540.19\t\t   |-312,582.26\t\t    |674.10s\t\t     |495.62s\t\t |\n\n(Any n=20, d=0.5, 1.5 instance are omitted due to character limits)\n\nThe reason for different BIC score: for the updated version, we employ warm-start initialization for the Master IP solution. The BIC score differences arise because the warm-start initialization causes the IP callback to find different cycles than before, which in turn selects different clusters and ultimately yields varied results.\\\nConcerning the computational complexity of our row and column generation scheme for solving the Master IP, while IPs are generally not polynomial-time solvable, our framework demonstrates strong practical performance through a careful problem formulation. Although we cannot provide theoretical complexity guarantees for the pricing problem (like many practical problems solved using IP), empirical evidence suggests our method represents one of the most promising score-based approaches.\n\n**3. Tie-breaking Strategy**\\\nFor the tie-breaking strategy in subgradient computation, our manuscript currently implements tie-breaking up to the ranking function implemented in numpy *argsort* function. In response to the reviewer's suggestion, we've conducted additional experiments using the strategy proposed in [1]. This investigation confirms the robustness of our approach to different tie-breaking methods. In the table, 'Previous' represents our original implementation, and 'Strategic' represents that strategy in [1] is used for breaking ties.\n| n, N, d | BIC Score (Previous) | BIC Score (Strategic)|\n|----------------------|--------------------|---------------------|\n| 15, 5000, 1     | -101954.63         | -101935.35          |\n| 15, 5000, 2      | -101646.26         | -101680.72          |\n| 15, 20000, 1      | -410415.68         | -410364.40          |\n| 15, 20000, 2      | -400909.42         | -401300.69          |\n| 25, 5000, 1       | -170279.79         | -170367.37          |\n| 25, 5000, 2       | -171333.90         | -171165.26          |\n| 25, 20000, 1      | -683475.44         | -682028.55          |\n| 25, 20000, 2      | -689041.34         | -689041.34          |\n|Average      |-312,582.26\t\t|-312,621.99\t\t|\n\n(All n=20, d=0.5, 1.5 instances are omitted due to character limits)\n\n\n| n, N, d | Average Pricing Time (Previous) | Average Pricing Time (Strategic) |\n|----------------------|----------------|----------------|\n| 15, 5000, 1       | 0.04s          | 0.41s          |\n| 15, 5000, 2       | 0.05s          | 0.40s          |\n| 15, 20000, 1     | 0.13s          | 1.85s          |\n| 15, 20000, 2      | 0.15s          | 1.74s          |\n| 25, 5000, 1      | 0.13s          | 3.62s          |\n| 25, 5000, 2      | 0.20s          | 3.24s          |\n| 25, 20000, 1      | 0.44s          | 17.10s         |\n| 25, 20000, 2     | 0.68s          | 0.68s          |\n|Average\t     |0.19s\t\t|3.08s\t\t|\n\n**4. Score-based IP Approach v.s. Constraint-based Approach**\\\nTo better motivate our score-based IP approach, we emphasize that constraint-based methods, while often more scalable in practice, operate sequentially and are susceptible to error propagation [2]. In contrast, the holistic score-based formulation avoids this issue by considering the problem globally, though at potentially greater computational cost. The IP framework transforms the original problem complexity into that of solving the optimization problem, which our method specifically aims to address. We acknowledge that constraint-based methods like PC can be more scalable for sparse data, though they may still exhibit exponential worst-case complexity [3].\n\n**5. Effect of the Supermodular Term in Scoring Function**\\\nWe clarify the relationship between supermodular term and DCA performance: In our current implementation, each DCA iteration optimizes a convex program formed by linearizing the concave (supermodular) component of the pricing objective. The presence of additional supermodular penalty term (for $i\\gets J$) is of the form:\n$$L^i(J)=\\Lambda(a_i-1)\\prod_{j\\in J} a_j$$\nwhere $\\Lambda$ is a hyperparameter and $a_j,j=1\\ldots,n$ is the number of values can that variable $x_j$ can take. This exponetially growing penalty term is hard to approximate using a single linear approximation.\n\n**6. Error Bars of the Results**\\\nRegarding error bars, we initially omitted them due to space constraints but now include them in the table below to demonstrate the robustness of our results. We will include the complete results with error bars in the appendix of the manuscript.\nIn the table below, each row represents the range for the resulting BIC score across 10 instances (generated with different modeling coefficients). \n| n, N, d      | BIC score range              |\n|--------------|-------------------------------|\n| 15, 5000, 1  | [-108960.02, -97244.45]       |\n| 15, 5000, 2  | [-105363.73, -98345.70]       |\n| 15, 20000, 1 | [-416161.99, -400230.34]      |\n| 15, 20000, 2 | [-422703.77, -389511.92]      |\n| 25, 5000, 1  | [-175506.33, -164339.46]      |\n| 25, 5000, 2  | [-176386.77, -166138.00]      |\n| 25, 20000, 1 | [-719570.68, -655633.66]      |\n| 25, 20000, 2 | [-705351.05, -664208.20]      |\n\n(GOBNILP results are omitted here due to character limits)\n\n**7. Other**\\\n-We appreciate your suggestion regarding the verification of local minimality and the restart of DCA. Given the current time limitations, we will conduct further tests but are unable to include the results in the current version.\\\n-The larger instances referenced in Section 5.1 were generated following the methodology in Section 5.2 with $N$=5,000. \\\n-We use the Python implementation of GOBNILP for comparison in our numerical studies.\\\n-Finally, we confirm that our work is primarily methodological in nature without direct social impact, and sincerely apologize for any oversight in not explicitly declaring this in the checklist.\n\nWe thanks for the reviewer's insightful comments, and would be happy to provide any additional clarification or information as needed.\n\n\n## References\n[1] Marwa El Halabi, George Orfanides, and Tim Hoheisel. Difference of sub-modular minimization via DC programming. In *International Conference on Machine Learning*, pages 9172–9201. PMLR, 2023.\\\n[2] Federica Cugnata, Ron S Kenett, and Silvia Salini. Bayesian networks in survey data: Robustness and sensitivity issues. *Journal of Quality Technology*, 48(3):253–264, 2016.\\\n[3] Thuc Duy Le, Tao Hoang, Jiuyong Li, Lin Liu, Huawen Liu, and Shu Hu. A fast PC algorithm for high dimensional causal discovery with multi-core PCs. *IEEE/ACM transactions on computational biology and bioinformatics*, 16(5):1483–1495, 2016."}, {"review_id": "oXtBQ8oavb", "reviewer": "Reviewer_VSAx", "summary": "The authors propose a new algorithm for Bayesian Network Structure Learning (BNSL). The proposed method is based on the score-based Integer Programming (IP) approach and solves the formulated IP using row and column generation methods. The column generation method requires solving a set optimization problem called the pricing problem, which can be expressed as the difference of two submodular functions, and the problem can be solved efficiently by applying the Discrete Convex Algorithm (DCA) with Lovász extension. Experiments on discrete and continuous synthetic and real data confirm that the proposed method outperforms existing methods.", "questions": "- What part of the title refers to “Inexact”? Is the information “Inexact” important enough to include in the title? This part should be explained more emphatically in the text.\n- I did not understand the explanation in lines 145-148, where the output of Algorithm 1 is always an integer solution. I would appreciate a more detailed explanation. In particular, the part after “all extreme point solutions ...” should be explained.", "limitations": "- As the authors mention in the Conclusion, scalability is limited.\n- Theoretical guarantees regarding how well DCA can solve column generation problems are difficult to obtain.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 2, "strengths_and_weaknesses": "# Strengths\n- The proposed method, although heuristic, is practical and sound for the important problem of BNSL.\n- The superiority of the proposed method is experimentally confirmed on various data sets.\n- The paper is well-structured and very easy to read.\n\n# Weaknesses \n- In the machine learning field, the procedure of solving NP-hard problems by formulating them as Difference of Convex or Difference of Submodular is often used. Therefore, I feel that the proposed method, while effective for the problem, is not that novel.\n- DCA is not necessarily a globally optimal algorithm, and there are no theoretical guarantees regarding the quality of the solution.\n- Section 4.1 describes the importance of the initial value setting, but does not provide experimental results for this part, so there is a lack of information on what guidelines users of this algorithm should use for their initial value strategy.", "was_revised": true, "final_justification": "My concerns regarding\n* the global optimality of the DCA,\n* the choice of initial values, and\n* the other questions I raised\n\nhave all been satisfactorily resolved. Consequently, I have increased my assigned score.", "rebuttal": "We sincerely appreciate the reviewers' insightful comments. Below we provide our response:\n\n**1. Novelty of the Difference-of-Submodular Formulation**\\\nOur work introduces a key methodological innovation through a novel decomposition of the pricing objective. To the best of our knowledge, this represents the first successful application of difference-of-submodular optimization techniques to the Bayesian Network Structure Learning (BNSL) problem.\n\n**2. Global Performance of Our Algorithm**\\\nThe Difference-of-Convex Algorithm (DCA) serves as an efficient heuristic for the pricing subproblem, with its practical effectiveness empirically validated through our numerical study. Our implementation strategically employs DCA in column generation subproblems to rapidly produce high-quality columns. As is common in vehicle routing [1,2] (and pointed out by Reviewer krir), to turn our approach into an exact solution approach for BNSL, one may apply inexact column generation first to quickly generate good columns and then apply exact column generation to conclude that there is no need for additional columns. We believe our results empirically show that DCA serves this purpose, and our framework can be potentially applied to develop such an exact approach.\n\n**3. Initialization Strategies for DCA**\\\nWe apologize for the previous omission of the comparison. We conducted some addtional experiments. The tables below present a comparison of BIC scores and computational time obtained from our CG-DCA framework by applying the three initialization methods we mentioned in our manuscript. These additional results will be added into Appendix.\n\n| n, N, d      | hybrid       | random       | warmstart    |\n|--------------|--------------|--------------|--------------|\n| 15, 5000, 0.5 | -101818.64   | -101818.64   | -101819.58   |\n| 15, 5000, 1   | -101954.63   | -101939.75   | -101975.70   |\n| 15, 5000, 1.5 | -101411.14   | -101376.34   | -101509.59   |\n| 15, 5000, 2   | -101646.26   | -101556.54   | -101608.25   |\n| 15, 20000, 0.5 | -406345.25   | -406345.25   | -406346.23   |\n| 15, 20000, 1  | -410415.68   | -410415.68   | -410382.11   |\n| 15, 20000, 1.5 | -405344.24   | -405487.19   | -405599.33   |\n| 15, 20000, 2  | -400909.42   | -400861.51   | -401710.97   |\n| 20, 5000, 0.5 | -136605.28   | -136605.28   | -136734.74   |\n| 20, 5000, 1   | -137222.13   | -137212.97   | -137363.27   |\n| 20, 5000, 1.5 | -136624.81   | -136577.03   | -136957.87   |\n| 20, 5000, 2   | -136849.69   | -136610.04   | -136750.44   |\n| 20, 20000, 0.5 | -535047.92   | -535047.92   | -535183.20   |\n| 20, 20000, 1  | -541850.16   | -541887.02   | -542168.27   |\n| 20, 20000, 1.5 | -538682.19   | -538693.99   | -539863.90   |\n| 20, 20000, 2  | -535301.38   | -534115.87   | -534939.03   |\n| 25, 5000, 0.5 | -170751.53   | -170751.53   | -170751.74   |\n| 25, 5000, 1   | -170279.79   | -170296.43   | -170302.01   |\n| 25, 5000, 1.5 | -171308.77   | -171157.64   | -171289.38   |\n| 25, 5000, 2   | -171333.89   | -171105.48   | -170760.08  |\n| 25, 20000, 0.5 | -677813.87   | -677813.87   | -677812.65   |\n| 25, 20000, 1  | -683475.44   | -683189.30   | -683826.26   |\n| 25, 20000, 1.5 | -676190.15   | -676267.45   | -676649.76   |\n| 25, 20000, 2  | -689041.33   | -688690.44   | -687147.17  |\n| **Average**   | **-312,582.26** | **-312,553.87** | **-312,712.62** |\n\n| n, N, d      | hybrid    | random    | warmstart |\n|--------------|--------------|--------------|--------------|\n| 15, 5000, 0.5 | 4.06s     | 4.04s     | 1.17s     |\n| 15, 5000, 1   | 25.76s    | 29.40s    | 4.49s     |\n| 15, 5000, 1.5 | 76.72s    | 87.43s    | 10.00s    |\n| 15, 5000, 2   | 170.39s   | 295.25s   | 23.74s    |\n| 15, 20000, 0.5 | 10.17s    | 10.94s    | 3.51s     |\n| 15, 20000, 1  | 36.71s    | 36.57s    | 10.63s    |\n| 15, 20000, 1.5 | 104.25s   | 109.73s   | 23.95s    |\n| 15, 20000, 2  | 229.80s   | 245.11s   | 41.94s    |\n| 20, 5000, 0.5 | 9.53s     | 9.60s     | 2.42s     |\n| 20, 5000, 1   | 84.34s    | 90.87s    | 12.51s    |\n| 20, 5000, 1.5 | 312.67s   | 430.37s   | 48.78s    |\n| 20, 5000, 2   | 1039.53s  | 1770.71s  | 95.36s    |\n| 20, 20000, 0.5 | 18.42s    | 18.54s    | 8.71s     |\n| 20, 20000, 1  | 133.69s   | 134.79s   | 34.84s    |\n| 20, 20000, 1.5 | 425.17s   | 504.29s   | 444.37s   |\n| 20, 20000, 2  | 997.99s   | 1710.91s  | 211.44s   |\n| 25, 5000, 0.5 | 21.97s    | 22.35s    | 3.87s     |\n| 25, 5000, 1   | 261.84s   | 304.12s   | 33.59s    |\n| 25, 5000, 1.5 | 949.90s   | 2434.03s  | 131.69s   |\n| 25, 5000, 2   | 3054.63s  | 6037.58s  | 224.73s  |\n| 25, 20000, 0.5 | 50.22s    | 53.13s    | 12.48s    |\n| 25, 20000, 1  | 404.43s   | 642.87s   | 78.21s    |\n| 25, 20000, 1.5 | 1372.59s  | 3688.33s  | 227.15s   |\n| 25, 20000, 2  | 3189.98s  | 6674.62s  | 366.46s  |\n| **Average**   | **495.62s** | **897.93s** | **88.71s** |\n\nFrom the tables, we can observe that:\n- Random initialization yields highest (best amony three)\nBIC scores in 11 out of 24 instances among three methods but requires the most time in 22 out of 24 instances.\n- Warm-start initialization achieves the lowest (worst among three) BIC scores in 16 out of 24 instances and is the most time-efficient in 23 out of 24 instances.\n- The hybrid approach strikes a balance between optimality and computational efficiency, since it focuses on local refinement around the current best pattern (through warm-start initialization) while building upon the foundation of global exploration (through random initialization) at the early stage.\n\n**4. Clarification of \"Inexact\" in Title**\\\nThe term \"inexact\" in our title specifically refers to our heuristic approach for solving the column generation subproblems (pricing problems) approximately during the optimization process. We will make sure that this is explicitly clarified in the published version if accepted.\n\n**5. Integral Extreme Point Solutions**\\\nNote that the epigraph of the Lovász extension of a submodular function $h$ charaterizes the convex envelope of (the standard binarization of) the set function $h$ [3]. Since $g^L(\\cdot)$ is the Lovász extension of the submodular function $g$,  the expression $g^L(\\cdot) - \\langle y^t, \\cdot \\rangle$ is itself the Lovász extension of a submodular function (specifically, $g$ plus a modular function). Consequently, all extreme points of the epigraph of $g^L(\\cdot)-\\langle y^t,\\cdot\\rangle$ has binary $x$ components. Therefore, if we optimize $g^L(\\cdot)-\\langle y^t,\\cdot\\rangle$ exactly using the cutting plane method in Line 5 of Algorithm 1 and always retrieve an extreme point optimal solution, the solution $x^{t+1}$ would be binary/integral (subject to numerical errors). We will clarify this in an updated version of our manuscript.\n\nWe appreciate the opportunity to clarify these points and remain open to further discussion or suggestions.\n\n## References\n[1] Pessoa, Artur, et al. \"A generic exact solver for vehicle routing and related problems.\" *Mathematical Programming* 183.1 (2020): 483-523.\\\n[2] You, Zhengzhong, and Yu Yang. \"RouteOpt: An Open-Source Modular Exact Solver for Vehicle Routing Problems.\" *Available at SSRN* 5314242 (2025).\\\n[3] Lovász, László. \"Submodular functions and convexity.\" *Mathematical Programming The State of the Art: Bonn 1982*. Berlin, Heidelberg: Springer Berlin Heidelberg, 1983. 235-257."}]

Author rebuttal:

Meta-review: This paper introduces an inexact column generation framework for Bayesian Network Structure Learning (BNSL). In particular, the pricing problem is reformulated as a difference-of-submodular optimization problem and solved approximately with the Difference-of-Convex Algorithm (DCA). Empirical results demonstrate that the proposed method provides significant speedups over exact column generation and yields higher-quality solutions on Gaussian data compared with state-of-the-art score-based approaches. 

Overall, this paper reads as an excellent fit to NeurIPS: it is novel within the BNSL literature and well motivated. The presentation is also rigorous and makes the paper particularly relevant both to ML and optimization audiences. The review team also appreciated the theoretical grounding of the DS reformulation and the clarity of exposition. 

There were, however, important points raised during the discussion. Specifically, performance degrades substantially on discrete data, scalability is still restricted to relatively small graph sizes, and the heuristic nature of DCA and its structure make it difficult to provide strong theoretical guarantees. Initial reviewer concerns centered on runtime comparisons, initialization strategies, and experimental clarity. 

Nonetheless, the author-reviewer discussions were very productive, and all questions were addressed nicely in the rebuttal with additional results and clarifications. I must note that all are minor and can be easily incorporated into the current draft. I appreciate the author-reviewer engagement and the thorough authors' responses.

Given the combination of novelty, technical soundness, and demonstrated benefits, the review team supports the acceptance of the paper at NeurIPS. We all believe it is an impactful work, specifically by providing a valuable contribution at the interface of optimization and causal discovery.
Decision: Accept (poster)