ExDBN: Learning Dynamic Bayesian Networks using Extended Mixed-Integer Programming Formulations
Authors: Pavel Rytíř
Aleš Wodecki
Georgios Korpas
Jakub Marecek
Venue: Accepted by TMLR ()

Abstract: Causal learning from data has received much attention recently. Bayesian networks can be used to capture causal relationships. There, one recovers a weighted directed acyclic graph in which random variables are represented by vertices, and the weights associated with each edge represent the strengths of the causal relationships between them. This concept is extended to capture dynamic effects by introducing a dependency on past data,  which may be captured by the structural equation model. This formalism is utilized in the present contribution to propose a score-based learning algorithm. A mixed-integer quadratic program is formulated and an algorithmic solution proposed, in which the pre-generation of exponentially many acyclicity constraints is avoided by utilizing the so-called branch-and-cut (``lazy constraint'') method. Comparing the novel approach to the state-of-the-art, we show that the proposed approach turns out to produce more accurate results when applied to small and medium-sized synthetic instances containing up to 80 time series. Lastly, two interesting applications in bioscience and finance, to which the method is directly applied, further stress the importance of developing highly accurate, globally convergent solvers that can handle instances of modest size.

Published in Transactions on Machine Learning Research (11/2025)
ExDBN: Learning Dynamic Bayesian Networks using Ex- tended Mixed-Integer Programming Formulations
Pavel RytfF
rytir. pptrr.pavel@fel.con.co
Caech Technical University in Progue
Ales Wodecki
wodecki.oess@eeL.com
Caech Technical Universiby in Prague
Georgios Korpas
cenngios.orpos@@hsse..mmm
HSBC Holdings Plc., Singopore
Czech Technicol University in Prague
Archimedes AI, Athena Research Center, Greece
Jakub Mareček
Caech Technicol University in Prague
Reviewed on OpenReview: https: // owtpr//veeervvee....m net/ forum? id= (d=46JJ29
Abstract
Causal learning from data has received much attention recently. Bayesian networks can be used to capture cansal relationships. There, one recovers weighted directed acyelic graph in which random variables are represented by vertices, and the weights associated with each edge represent the strengths of the causal relationships between them. This concept is extended to capture dynamic effects by introducing a dependency on past data, which may be captured by the structural equation model. This formalism is utilized in the present contribution to propose score-based learning algorithm. A mixed-integer quadratic program is formulated and an algorithmic solution proposed, in which the pre-generation of exponentially many acyclicity constraints avoided by utilizing the so-called branch-and- cut. (flazy constraint") ") method. Comparing the novel approach to the state-of-the-arl, we show that the proposed approach turns out to produce more accurate results when applied to small and medium-sized synthetic instances containing up to 80 time series. Lastly, two interesting applications in bioscience and finance, to which the method is directly applied, further stress the importance of developing highly accurate, globally convergent. solvers that can handle instances of modest size.
1 Introduction
Dynamic Bayesian Networks (DBNs) were conceived by Dagum et al. (1992) to unify and extend traditional linear state-space models, linear and normal forecasting models (such as ARMA), and simple dependency models (such as hidden Markov models) into a general probabilistic represontation and inference mechanism for arbitrary nonlinear and non-normal time-dependent domains. MMurphy| (2002) made important contribu- tions the inference. The network may be used to compute posterior probabilities for different situations, or the learned graph structure may be inspected, and dependencies of particular interest analyzed. This explainability is a key benefit of DBN. To the present day, there are no scalable methods for learning DBNs.
There exist a number methods for learning DBN, One can learn underlying contimuus dynamics rep- resented by stochastic differential equations (Bellot et al., (021), or related dynamical systems (Murphy) 2002). Another option assume priori knowledge about timo-lagged data. and incorporate this knowl- edge into the solver (Sun et al. 2021 Furthermore, deal with the general problem and propose local

---

Published in Transactions on Machine Learning Research (11/2025)
methods (Pamfil et al. 2020 Gao et al. 2022). which can scale further at the cost of some loss of accuracy. Under the assumptions of Markovianity and laithfulness, conditional independence-baset methods, such as PCMCI+ (Runge, 2020), , can be applied. Note that many of the previous works also combine several of these approaches to design solvers that are efficient and applicable to a wide range of applications. However, it should be noted that many methods may not identify DAG representations of causal dependencies under certain conditions (Kaiser & Sipos] 2022 Reisach et al. 2021). One of the possible causes is that many of them only converge to finst-order stationary point for the corresponding non-convex optimization problem, which is known (Daniely & Shalev-Shwartz 2016) to be hard to solve.
In our contribution, we revisit the score-based learning of DBNs utilizing a directed acyelic graph (DAG) structure augmented by additional time lagged dependencies (Murphy) 2002 Dean & Kanazawa 1989 Assaad et al., 2022). We utilize mixed-integer programming to learn the underlying dynamic Bayesian networks. While all of the previous methods focus mostly on scaling with adequate precision by utilizing a varicty of heuristics, we focus on leveraging quadratic mixod-integer programs to find near global optimum solutions to score-based DAG learning problem, which results in high-quality reconstruction of the DAG.
Notice that the number of constraints needed enforce the structure of directed acyclic graph ("cycle- exclusion constraints") is super-exponential in the number of time series on the input. This is sometimes referred to as the curse of dimensionality. We tackle the curse of dimensionality by avoiding the pre-generation of these cycle- exclusion constraints and subsequently adding these constraints only needed. is shown that, only a small number of these constraints are actually needed to ensure the acyclicity of the directed graph, which allows for a large speedup over the naive version of the algorithm that pre-generates all of the constraints. Additionally, this technique allows the method to solve problems much larger compared to the naive counterpart, since the memory consumption is greatly reduced. The formulation and its implementation are easily reproducible, making it accessible to a wide range of potential practitioners.
Learning DBNs has been successfully applied to a variety of problems, many of which are related to applica- tions in medicine (Zandoná et al. 2019 van Gerven et al. 2008 Michoel & Zhang 2023 Zhong et al. 2023 In addition to medical a applications, dynamic Bayesian net works are widely used in econometrics (Demiralp & Hoover 2003) and financial risk modeling (Ballester al. 2023). The applications of DBN are not limited to the ones mentioned and the curious reader may refer to (Kunurtsev et al.) 2024) for further details. After
we present our approach, we revisit these applications in Sections 4.5 and 4.4
2 Problem Formulation
Before formulating the problem of score-based Bayesian network learning as mixed-integer program, let us describe the problem as a structural vector autoregressive model DDmirall & Hoover 2003 Kilian 2011 Let d E N be the number of variables in an autorogressive model. Let T E N be the total number of time periods and let X,, be a set of endogenous variables at the given time f, where i E 21. d} and E (1. T}.
We assume that the variables interact with each other both simultancously and with time lags. These interactions can be captured by the following equation:
X. X,W + AI + Xr-2A2 +.... + X. X. + Zz,
(1)
where p E N is the autoregressive order, and
W E Rxdd A E dxd Z, E Rd, (1,
Simultaneous interactions are called intra-slice and are represented by the matrix W and time-lagged in- teractions are called inter-slice and are represented by the matrices A,, where i is the number of timesteps (length) of the time lag. The vector Z, is an error term. We will not assume any particular distribution of error terms.
We will assume that there is no cycle in intra-slice interactions. In other words, we view W as an adjacency matrix of a directed graph G, then a directed acyclic graph (DAG). The matrices A, could also be viewed as adjacency matrices of directed graphs, but we do not make assumptions about these graphs. Note that non-linear autoregressive models can also be formulated in an analogous way.
2

---

Published in Transactions on Machine Learning Research (11/2025)
Our goal is to learn the matrices W, which corresponds to DAG, and A, for i = 1,.... p. Assume that we have n data samples of each of d random variables organized into a data matrix X E IRaxd Then, the following equation holds for the solution W, A,, i P.
X XW V1A1 Ap + Z.
(2)
where Y are time lagged version of X for ,,..., and Z Raxd is a matrix with error terms for all samples.
To maximize the fit of the data over the model, we want to minimize the Frobenius norm of the error matrix
Z. We formulate it as the following cost function:
J
= | XW Y1A, - ................................ -...--........... GA + WIll + (1All + ... ||Ap)), (3)
where r denotes Frobenious matrix norm and ||-Il denotes an arbitrary matrix norm and X, > are sufficiently small regularization coefficients. The data fitting problem as follows
min W
J JW, . Ap)
W acyclic, W E Raxd
(4)
A, E Rdddd ,iE (1,2.....
Remark The identifiability of W and A, using 4 has been studied for Gaussian and non- Gaussian noise. Regardless of noise, the identifiability of A is consequence of the basic theory of autoregressive models Killian 2011). The identifiability of W is a bit more involved and must be separated into the Gaussian and non- Gaussian case. However, in either case, idenlifiabiily is possible under mild condilions Hyvürinen et ol. 2010: Peters & Bühhmann 2012).
Note that if the autoregressive order p is set to 0, ExDBN could be used to learn (static) Bayesian Networks. Dynamic Bayesian Networks differ from Bayesian Networks (BNs) primarily in the statistical structure of their data. While BNs are trained on sets of independent and identically distributed samples, DBNs rely on trajectories of time series where successive samples may be temporally dependent. The inference goal of DEN models typically forecasting, whereas for BNs it generative modeling. To evaluate learning algorithms, different types of benchmarks are required for BNs and DBNs. Therefore, in this work, we focus strictly on DBNs, i.c., autoregressive orders p 2
2.1 Brief Introduction To Mixed Integer Quadratic Programming
To better frame the content Section 2.2 we provide a short introduction to mixed-integer quadratic pro- gramming. An optimization problem, called mixed-inte mixed-integer quadratically constrained quadratic program (MIQCQP) if it is of the form
min ZER
x TTQx FT
(5)
s.t. x² Q2 + 95 S ai
(6)
Ax b,
(7)
:; E T for $ I
3, E 7, for E
where Q, Q: E Raxn, 9,4: E R", A Rmxn, a R R E. for some m, m,n, n, E.NENN Expression 55 is often
called the objective, cost, or loss function, Inequality are the linear constraints, and index set denotes the integral variables.
6 represents the quadratic constraints, Inequality
Mixed-intege quadratic programs have been shown be in NP DDll Pia et al. 2014). which often leads to exhaustive demand for computational resources. The algorithms used to solve MIQP typically
3

---

Published in Transactions on Machine Learning Research (11/2025)
branch-and-bound cutting plane Dakin 1965
Bonami et al. 2009 Westerlund & Pettersson 1995
Kronqvist et al. 2015). Both of these algorithmic treatments are often employed together, often with the addition of a presolving step, the use of heuristics and parallelism. The aforementioned allows many modern solvers to solve even large problems despite the NP hardness. Some of these solvers are open source (like SCIP (Achterberg 2009) and GLPK (GNU Project 2025 and others are commercial (GUROBI (Gurobi Optimization, LLC 2024) and CPLEX (TBM (025)). ). The powerful infrastructure present in these solvers can be made use of together with additional problem-specific modifications to deliver high-quality solutions.
Due to the exhaustive nature of the algorithms mentioned in the previous paragraph, global convergence is guaranteed (Belotti et al.] (013).. Furthermore, convergence to the global solution may be tracked and the error estimated by computing the dual problem of 5 19. computed the so called MIP gap as follows
The dual of the problem is then used to
MTP gap 1J (x*) - Janal (v*)
(10)
(xx))
where x² and y" are the current best solutions of the primal and dual problems respectively, and J and J are the cost functions of the primal and dual problems, respectively. The MIP GAP ensures that we can assess the quality of the minimization ion during solution time and lerminate the computation when the result is good enough (small enough MIP GAP). Furthermore, if the gap reaches 0 at any point, we are sure that the current solution is global optimum.
2.2 Mixed Integer Quadratic Programming Formulation
Formulating the learning problem allows us to use a globally convergent algorithm. Key formulations for the learning of directed acyelic graphs, closely related problem, have been proposed by Manzour al. (202]
Xu et al. b2024; Jaakkola et al. (2110).
First, we define the program variables. For i,j = l...., d, we define the binary variable eij whose value is interpreted as follows: Eij equals 1 if and only if there is an oriented edge from node to node j both in the intra-slice graph.
Similarly, we define binary variables for edges connecting nodes from inter-slices graphs to the nodes of intra-slice graph. For i,j 1, 1,...,d d and I = -, where is the autoregressive order. We define the binary variable e whose value is interpreted as follows: et equals if and only if there is oriented edge from node i in the inter-slice graph, corresponding to time lag to the node in the intra-slice graph
With each of the above variable ey and e/y, we associate a continuous variables wy and a'y- respectively, which encodes the weight of the corresponding edge.
Using these variables, the scoring function of problem becomes the following MTQP objective:
J (W, A AP)
Xij
Wik ΣΣ Y.naae + REG,
(11)
k=1
where REG is the regularization part of the objective function. We use two different regularization functions LO and L2 defined as follows. Regularization LO is defined as:
n
n
n
REG
1
(12)
Regularization L2 is defined as:
REG A
wi
(13)
where À > 0 and 1 >> 0 are regularization coefficients.
4

---

Published in Transactions on Machine Learning Research (11/2025)
Next, define the constraints of the optimization problem. In order to ensure that W encodes directed
acyclic intra-slice graph. We add the following constraints (Dantzig et al. 1954 that ensure that W is acyclic.
Let denote the set of all directed cycles in the complete directed graph on d vertices without loops. Then we add the following constraints:
Σ (3,)
ei S - for all C Ceec
(14)
Note that, the number above constraints exponential. We will discuss in the following section how to deal with such number of constraints.
In order to relate the binary variables egj and et with the weight variables wig and aly, we introduce the following constraints:
Wij Sceij:
(15)
Wij 2 - oij for all i,j E (1,2.....d) d).
di Scejj for all k,j {,, .,,,,.......... d},
(16)
d 2 ceis for all / {1, (1.
where e >> 0 is the maximal admissible magnitude of any weight. In our numerical experiments, we choose large enough to not affect the result of the learning.
3 Algorithmic Implementation Using Branch-and-Bound--ndd and-Cut
One of our main contributions is the development of a. brauchannd-houn--nd--or algorithm to solve the formulation mentioned above. The Bauuch-and-Bound-anddee mixed-integer optimization problems by combining three key steps: solving linear relaxations to obtain bounds, branching variables when solutions fractional, and tightening relaxations with valid inequal- ities. In addition, users supply so-called 'lazy" constraints, which constraints not included in the initial model but enforced only when violated by candidate integer solution during the search. This lets the solver keep the model smaller and only add these constraints when necessary, further improving efficiency while guarantceing correctness.
algorithm (Achterberg 2007 e.g.) solves
Since the acyclic coists 11 need to be imposed only for the edges of the graph representing the intra- slice level, all of what follows is only applied to the intra-slice graph. While we leverage the traditional branch- and-bound approach as described in (Achterberg) 2007] c.g.), we incorporate cycle exclusion constraints (14) using "lazy" constraints. These are only enforced once an integer-feasible solution candidate is found. If violation of lazy constraint occurs, the constraint. is added across all nodes in the branch-and-boung tree. At the root node, only O ([E)) constraints (15) and (16) are initially used. Cycle-exclusion constraints (14) are added later. Note that this method is not a heuristic and does not lead to a possibly harmful reduction (or extension) of the solution space leading to omitting possible solutions or returning solutions which are not DAGs. Furthermore, is shown that the number of constraints that are actually needed in computation is many orders of magnitude less than the number of all possible constraints.
Once mixed-integer feasible solution candidate is identified, detecting cycles becomes straightforward using depth-first search (Even & Even 2011 pp. 46-48). If a cycle is detected, the corresponding lazy constraint [14] is added to the problem. The depth-first search (DFS) algorithm solves the problem of cycle detection in a worst-case quadratic runtime relative to the number of vertices in the graph, which contrasts with algorithms that separate related inequalitics from a continuous relaxation (Borndörfer et al. 2020 Cook et al.) 2011), such as the quadratic program in our case. Three variants of adding lazy constraints for the problem tested.
5

---

Published in Transactions on Machine Learning Research (11/2025)
Adding lazy constraint only for the first cycle found.
Adding a lazy constraint only for the shortest cycle found.
Adding multiple lazy constraints for all cycles found in the current iteration in which an integer- feasible solution candidate is available.
The third mentioned variant was found to consistently deliver the best results, despite (ccherberg. 2007
Chapter 8.9). See Table 2 for numerical results. Therefore, is applied in all the numerical tests that follow.
4
Numerical Experiments
n=50
n=1000
ER-3-1-var
Structural Hamming distance SHD
ER-3-1
0
SF-3-1
0
ER-2-2
500
0
5
10
15
20
25
10
20
25
10
15
20
25
10
20
25
10
15
20
25
(number of nodes)
dynotears
exdbn
lingam
nts_notears
Figure 1: SHD for test cases using ER-3-1, SF-3-1, ER-2-2 random ensembles, where the first number is the edge-vertex ratio on intra graph, the second is the autoregressive order. The edge-vertex ratio of inter graph is always 1. For the problem denoted by var suffix, the variance of Gaussian noise is randomly sampled from uniform distribution on interval (0.6, 1.2) for each variable. For the other problems, the variance is set to 1. The mean, standard deviation, and maximum over 10 algorithm runs is depicted. Each run is performed on different randomly sampled dataset.
In recent years, many solvers have been developed to facilitate the graphical learning of Bayesian networks that represent causality (Pamfil et al. 2020 Hyvärinen et al. 2010 Malinsky & Spirtes 2018 Gao et al. 2022 Dallakyan 2023: Lorch et al. (021). Each of these solvers (including the one presented) faces the curse dimensionality, which somewhat restricts their applicability; thus, thorough testing is required. It is impossible lest the proposed solution against every solver that has been developed. However, there is significant branch of research that allows for direct comparison, and by the transitivity of results, this enables indirect comparison with many previous solvers.
Pamfil et al. 2020) have developed a locally convergent method, called DYNOTEARS, that learns causality as a Bayesian network that supersedes the solution methods previously developed (Hyvärinen et al.) 2010 Malinsky & Spirtes 2018 Zheng et al. 2018). Further developments based on previous publications include formulating the problem the frequency domain or defining differentiable Bayesian structures (Dallakyan)
2023 Lorch et al. 2021). In this section, provide comparison with DYNOTEARS, LiNGAM (Hyääinnen)

---

Published in Transactions on Machine Learning Research (11/2025)
al. 2010 and NTS-NOTEARS (Sun et al. 2021 and thus, by transitivity, with the methods documented
by Hyvärinen et al. 2010); Malinsky & Spirtes (0088).
n=
n=100
n=250
n=500
n=1000
1.0
ER-3-1-var
0.5
0.0
1.0
ER-3-1
0.5
0.0
F1 score
1.0
SF-3-1
0.5
0.0
1.0
ER-2-22
0.5
0.0
5
10
15
20
25
10
15
20
25
5
10
15
20
25
5
10
15
20
25
5
10
15
20
25
d (number of nodes)
dynotears
exdbn
lingam
nts_noteea
Figure 2: Fl score for test cases using ER-3-1, SF-3-1, ER-2-2 random ensembles.
4.1 Synthetic Data
One of the evaluations of ExDBN was performed the synthetic data generated by the following process. First, random intra-slice directed acyclic graph (DAG) generated using either the Erdös-Rény (ER) model or the scale-free Barabási Albert (SF) model. Then, the DAG weights were sampled uniformly from the union of intervals 2.0. -.55 U [0.5, 2.0].
Next, the inter-slice graphs were generated using the ER model. For each inter-slice graph, weights were sampled from the interval [-0.5a, -..2a] [0.2x, 0.5a], where a = 1/1²-1, ") 2 the decay parameter, and is the time lag of the slice. inter-slices.
corresponds the intra-slice, while E [l.....p] represents the
The data samples are then generated using the structural equation model equation and adding Gaussian noise with either variance 1 or different variance for each variable sampled uniformly from a given interval.
4.2 Benchmark Setup and Quantities of Interest
Let Wtrue denotes the adjacency matrix representing the intra-slice ground truth graph of given problem and let A denotes the adjacency matrices of inter-slice ground truth graphs of time-lagged interactions. In the case of synthetic data, Wtrue and At,ree known, since they were used for data generation.
Let Was and Ar,st be the solution of the learning problem. We will evaluate their quality by comparing to the ground truth and computing structural Hamming distance and Fj score.
We apply a small threshold 0.15 to West and Ar, We set all the elements that are smaller than 0. We do that in order to remove small rounding errors that would negatively affect computation of discrete metrics such structural Hamming distance.

---

Published in Transactions on Machine Learning Research (11/2025)
Structural Hamming distance (SHD) is defined follows:
SHD (Werue; At,true: West Ag,ett
distay (Werue: West) ) +
disti (Ageruer Acest), (17)
1=1 1,
where
if Cú + 0 and Dij + 0
0
if Ci = 0 and Dij = 0
distay (C,
(18)
2/2 if Ct # and Djs # 0 otherwise.
SHD used as a score that describes the structural similarity of two DAGs in terms of edge placement and is commonly used to assess the quality of solutions Zheng et al. 2018 Pamfil et al. 2020 Besides SHD,
true positive
precision true positive false positive
(19)
true positive
recall true positive false negalive"
(20)
are used (Andrews et al. 2024) to evaluate the quality of structural recovery. It is important to note that precision and recall isolate false positives and negatives, respectively, in contrast to SHD, where these quantities are both accounted for simultancously. The last metric that can be used to evaluate structural similarity is the Fi score and reads
2
F precision recall
(21)
Note that all of the quantities evaluated in over both inter-slice and intra-slice matrices between the solution and the ground truth.
19 and 21 are result of summing up all of the differences
4.3 Synthetic Benchmark Results
In the following benchmark, the generation methods described in Section 4.1 are used to compare ExDBN with DYNOTEARS, LiNGAM, and NTS-NOTEARS under the assumption of Gaussian noise. Even though the cost function is a maximum likelihood cstimator (see Section [ for non-Gaussian noise, we leave this evaluation for future publication. The performance was studied for different numbers of variables, samples, and graph generation methods, with the relevant metrics: SHD and F1 score, reported Figures 1 and 2
A statistical ensemble with 10 different seeds used for each experiment, and the mean, standard deviation, and worst-case values are reported in the plots. We used Gurobi MIQP solver. The running time for
SF-3-1
ER-2-2
200
100
150
SHD
100
50-
2000
4000
6000
2000
6000
Time limit
Time limit
Figure Comparison of ExDBN solution quality (SHD) and running time for the ER-2-2 and SF-3-1 ensembles with 25 variables and 250 samples. The mean, and standard deviation 10 algorithm runs is depicted. Each run is performed on different randomly sampled dataset.
8

---

Published in Transactions on Machine Learning Research (11/2025)
Number nodes
DYNOTEARS
ExDBN
LINGAM
NTS-NOTEARS
5
1.1
0.4
530.5 + 1504.2
0.7
74.7 = 52.6
7
2.0
0.8
2188.7 1
0.2
142.5
87.4
10
4.5
2.0
5163.1 3014.5
0.2
379.9 + 237.8
15
12.2
5.1
7214.6 66.4
1.0
875.6 = 490.1
20
11.8
7216.1 + 182.0
1.7 0.3
1571.4 953.4
25
4
7226.3 + 20.7
2.9 + 0.6
1975.7 + 1114.5
Table 1: Rumning times of the evaluated methods (mean + standard deviation) across all scenarios, sample sizes, and ten runs on different randomly generated datasets.
ExDBN was capped at 7200 seconds, and the memory was limited to 32 GB. The regularization applied in ExDBN needs to be scaled appropriately with the number of samples, as the optimal rogularization constant is assumed to be a decreasing function of sample size. We use this assumption as a non-strict guideline select the appropriate regularization for a given sample size. This follows from the requirement that regularization should remain proportionally small compared to the main objective expressed by equation [I1p. It was also found that switching from L1 to L2 regularization improves identification performance when the number of samples is large. When the ground-truth graph is unknown, the algorithm can be run for multiple values of X and 77, and the configuration yielding a better MIP GAP can be selected. For smaller sample sizes, L1 regularizal ion performs belter, whereas for larger sample sizes, L2 regularization lypically yields good results and is faster.
As noted in (Reisach al., 2021), the noise variances and data scale may be important for some algorithms perform well. We tested ExDBN on normalized data and noticed significant performance drop. Therefore, ExDBN is suitable for problems in which the data of the samples have true scale.
The results of the tests can be divided into two categories based on the variance of the applied noise. In the non-equal variance scenario, ExDBN performs better than DYNOTEARS (see Figures] 1 and 2 In the scenario where all variances are equal to 1, the results are more balanced between ExDBN and DYNOTEARS. LiNGAM and NTS-NOTEARS perform substantially worse than both ExDBN and DYNOTEARS. This can be attributed the fact that LiNGAM is optimized for non-Gaussian data, while NTS-NOTEARS is designed for non-linear models. Note that FxDBN shows a higher standard deviation in some scenarios with 25 variables. This deviation may be reduced by allowing more computation time for ExDBN, as suggested by Figure 3 which compares standard deviation and running time.
1.0
Algorithm
1200
nts its_notears
1000
dynotears
Algorithm nts, mos_nooears dynotears
800
milp
F1 score
0.6
lingam
SHD
600
0.4
milp
400
lingam
0.2
200
0-
30
30
40
50
60
70
80
Number of variables
Number variables
Figure 4: SHD and Fl score the SF-3-1 benchmark for larger problem instances. The number of variables d ranges from 30 to 80, and the number of samples n is set to 5d.
Focusing on the 50-sample case, while also considering the previous results, we observe that the performance gap between the solvers widens in favor of ExDBN the number of variables increases. In the lower-sample
9

---

Published in Transactions on Machine Learning Research (11/2025)
scenarios, ExDBN outperforms DYNOTEARS for many graph sizes average and consistently outperforms DYNOTEARS in the worst-case results (minimum or maximum, depending on the metric).
Note that the global convergence of the method, which is rooted in the fundamentals of mixed-integer quadratic programming, allows for increased computation time, leading to further improvements in the reported metrics. While some time-sensitive applications, such as short-term stock evaluation, may not be able to benefit from this, others such biomedical applications can, as computations lasting several days with measurably improved accuracy (monitored via the duality gap) often desirable. See Figure for comparison of rumning time and solution quality, and Table for comparison of the running times of the evaluated methods.
To evaluate the scalability of ExDBN, we conducted experiments on simplified benchmark based on SF-3-1. The number of variables ranged from 30 to 80, and the number of samples was set to 5d. The running time was limited to 22 hours, and the memory limit was increased 128 GB. See Figure f for the results. We can conclude that within these limits, ExDBN scales up 80 variables and performs better than the other compared methods this benchmark.
Lazy constraints added for
mean MIP GAP
all cycles found
0.335
the shortest cycle found
0.344
the first cycle found
0.346
Table 2: Evaluation of variants for adding lazy constraints on SF-3-1
We ran experiments on three variants adding lazy constraints on SF-3-1. The results are shown in They indicate that adding lazy constraints for all cycles found the current iteration performs best on average. See Figure for the number of added lazy constraints during the ExDBN run in each scenario.
Table 2
n=50
n=100
n=250
n=500
n=1000
40000
ER-3-1-var
20000
40000
Number lazy constraints
ER-3-1
20000
40000
5F-3-1
20000
40000
ER-2-2
20000
10
15
25
10
20
5
15
20
25
10
15
20
25
10
20
(number of nodes)
Figure 5: Number lazy constraints added by ExDBN using 'all cycles" variant on test cases using ER-3-1, SF-3-1, ER-2-2 random ensembles.
10

---

Published in Transactions on Machine Learning Research (11/2025)
4.4 Application in Finance
In financial services, there also several important applications. The original DYNOTEARS paper consid- ered model of diversification of investments in stocks based dynamic Bayesian networks. Independently, (Ballester et al. 2023) consider systemic credit risk, which is one of the most important concerns within the financial system, using dynamic Bayesian networks. They found that the transport and manufacturing sectors transmit risk to many other sectors, while the energy sector and banking receive risk from most other sectors. To a lesser extent, there is a risk transmission present between approximately 25% of the sectors pairs, and these network relationships explain between % and single systemic risks. Notice that these instances are much denser than the commonly used random ensembles.
0H9987
0H9987
1.0
0SABBF
0SABBF
GG6EBT
GG6EBT
В869AР
0.5
9869AР
0.5
06DABK
06DABK
DD359M
DD359M
0.0
0.0
EFAGG9
EFAGG9
48DGFE
8D8575
-05
48DGFE
8D8575
-0.5
NN7A8G
NN7A8G
6A516F
-1.0
6A516F
FH49GG
FH49GG
0H9987
2H6637
8A87AG
0SADBF
GG6EBT
BB09AP
06DABK
0D339M
FF667M
EFAGG9
48DGFE
8D8575
NNZA8G
6A516F
FH490G
2H6677
049987
2H6637
8A87AG
0SADOF
GO6EBT
BB6BAP
06DA8K
0D359M
FF667M
EFAGG9
48DGFE
8D8575
NNZ48G
6A518F
FH49GG
2H6677
(a) Intra-slice graph
(b) Lag inter-slice graph
0H9987
0.0008
2H66B7
8A87AG
0SABBF
GG6EBT
-0,0004
8869AР
0.0002
06DABK
DD359M
0.0000
FF667M
EFAGG9
48DGFE
0D8575
0.0004
NN2A8G
6A518
0.0006
PH49GG
2H6677
0159077
2H6687
BASTAG
OSABBF
GCREST
0dDABK
DD359M
FF667M
EFAGG9
48DGFE
6D8575
NNZABG
GAS16F
FMSGG
246677
(c) Lag 2 inter-slice graph
Figure 6: ExDBN Heatmaps of the adjacency matrices of the learned intra-slice and inter-slice graphs for the CDS datasct.
We elaborate on the example of (Ballester et. al. 2023 3, where 10 time series capture the spreads of 10 European credit default swaps (CDS). Considering the strict licensing terms Refinitiv, the data from (Ballester al., (202)) are not available from the authors, but we have downloaded 16 time-serics capturing the spreads of 16 European CDS with RED6 codes 05ABBE, 06DABK, 0H99B7, 2H6677, 2H66B7, 48DGFE, 6A516F, 8A87AG, 8B69AP, 8D8575, DD359M, EFAGG9, FF667M, FH49GG, GGGEBT, NN2A8G, from
11

---

Published in Transactions on Machine Learning Research (11/2025)
Jamuary 1st, 2007, September 25th, 2024. The codes 48DGFE, 05ABBF, 8B69AP, 06DABK, EFAGG9, 2H6677, FH49GG, and 8D8575 belong to the Banks sector. The codes GG6EBT, DD359M, and FF667M belong to the insurance sector. Finally, the codes 0H99B7, 2H66B7, 8A87AG, NN2A8G, and 6A516F belong to transportation and manufacturing.
0H9987
0119977
2H6887
2H66B7
0.75
BAB7AG
8A87AG
OSABBF
0.50
GGBERT
8869AР
В669АР
06DABK
DD359M
DD359M
FF667M
FF667M
0.00
EFAGG9
EFAGG9
48DGFE
48DGFE
0.25
8D8575
-0.5
8D8575
NN2A8G
NN2A8G
6A516F
6AS16F
FH49GG
-1.0
PH49GG
-0.75
2H6677
2H6677
0H9987
2H6687
BA87AG
OSABBF
GGGEBT
8B69AP
05DABK
DD359M
FF667M
EFAGG9
48DGFE
8D6575
NNZABG
6AS16F
FH496G
246677
0H9987
216837
BA87AG
0SABBF
G66E81
8869AP
06DAGK
DD359M
rT66714
EFAGG9
48DGFF
NNZA8G
6A516F
FH49GG
2H6677
800s
(a) Intra-slice graph
(b) Lag
inter-slice graph
0H9987
2H6697
8A87AG
05ABBF
0.1
06DABK
DD359M
7M
EFAGG9
48DGFE
-01
NN2A8G
6A516F
FH49GG
2H6677
049987
2H6687
8A87AG
0SABBF
6G6E8T
880802
06DABK
DD359M4
FF6677
EFAGG9
48DGFE
806575
NNZA8G
6A516F
FH49GG
2H6677
(c) Lag 2 inter-slice graph
Figure DYNOTEARS Heatmaps of the adjacency matrices of the Icarned intra-slice and inter-slice
graphs for the CDS dataset.
These data amount more than MB of time series data when stored as comma-delimited values plain text. Although the procedure for learning the dynamic Bayesian network in (Ballester et al, 2023) is rather heuristic,
solve the mixed integer programming (MIP) instance for the 16 European CDS in 30 min- utes. In the heuristic of (Ballester et al. 2023). they first perform unconditional independence tests each set of two time scries containing an original serics and lagged time scries, to reduce the subsequent number of unconditional independence lests performed. There are 45 unconditional and conditional independence tests performed first, to suggest another 200 conditional independence tests. We stress that the procedure of (Ballester et al. 2023) does not come with any guarantics, while our instance (11] is solved to a 20% MIP gap. The run-time of minutes to a solution with 20% MIP gap (using L. regularization) validates the scalability of mixed-integer programming solvers.
12

---

Published in Transactions on Machine Learning Research (11/2025)
For solution of auto-regressive order 2, see Figures 6a
6b and 6c We can see that each CDS mostly depends
on its own history. Furthermore, as expected, the identified dependencies are mostly within sectors, with two exceptions: FF667M depends on 05ABBF and a lag 2 dependency is observed, with 2H66B7 depending on 8B69AP. For a comparison, we also included a solution found by DYNOTEARS; see Figures ia 7b and 7c In this solution, we can observe that the Lag dependencies are very similar to those in the Ex DBN solution. There are many additional Lag 2 self-dependencies, which we found rather redundant given that the same Lag 1 dependencies already exist. When examining intra-slice dependencies, we can identify one dependency between different sectors [DD359M-05ABBR)
DYNOTEARS did not detect many strong intra- sector dependencies, unlike ExDBN. We would argue that the solution found by ExDBN is better, although this is somewhat subjective since the ground truth is unknown.
4.5 Application in Biomedical Sciences
In biomedical sciences, there is keen interest in learning dynamic Bayesian networks to estimate causal effects (Tennant et al. 2020) and identify confounding variables that require conditioning. A recent mela- analysis (Tennant et al. 2020) of 234 articles on learning DAGs in biomedical sciences found that the averaged DAG had 12 nodes (range: 3-28) and 29 arcs (range: 3-99). Interostingly, none of the DAGs were as sparse as the commonly considered random ensembles; median saturation was 46%, meaning that each of all possible ares appeared with a probability of 46% and did not converge to a global minimum of the problem.
As an example, consider recently proposed benchmark of (Ryšavý et al. 2024.. where the Krebs cycle is to be reconstructed from time series of reactant concentrations of varying lengths. There, DYNOTEARS cannot reach the (Ryšavý et al. 2024) Fl score of 0.5 even with a very long time series. In contrast, our method can solve instances (4) to global optimality. Using ExDBN, however, the global minimization is ensured given sufficient time and thus the maximun likelihood estimator is found. However, it should be noted that depending on the number of samples and noise, may be that even the maximun-likelihodd estimator is not sufficiently accurate. However, this does not reflect poorly on the method itself, but it is rather matter of the modification of data collection methods associated with the experiment. In one- one-hour time limit, ExDBN can find a solution with the 38% duality gap.
5 Conclusion
Dynamic Bayesian networks have wide-ranging applications, including those in biomedical sciences and computational finance, as illustrated above. Unfortunately, their use has been somewhat limited by the lack of well-performing methods to learn them. Our method, ExDBN, provides the best possible estimate of the DBN, in the sense of minimizing empirical risk 3 Significantly, our method does not suffer much from the curse of dimensionality, even for real-world dense instances, which are typically challenging for other solvers. This is demonstrated most clearly in the case of systemic risk transmission, detailed Section 4.4 where a solution with 20% MIP gap is found in 30 minutes. Additionally, the use of the guarantees on the distance to the global minimum (so-called MIP gap, available ahead of the convergence to the global minimum) provides significant tool for fine-tuning the parameters of the solver in the case of real-world application, where the ground truth is not available. Combined with global convergence guarantees of the maximum likelihood estimator, this provides a robust method with stale-of-the-art performance.
Acknowledgments
The authors acknowledge the support of National Recovery Plan funded project MPO (0273/24/21300) CEDMO 2.0 NPO. This work has been partially supported by project MIS 5154714 of the National Recovery and Resilience Plan Greece 2.0 funded by the European Union under the NextGeneration Program.
Disclaimer
This paper was prepared for information purposes and is not product of HSBC Bank Plc. or its affiliates Neither HSBC Bank Ple. nor any of its affiliates make any explicit or implied representation or warranty
13

---

Published in Transactions on Machine Learning Research (11/2025)
and none of them accept any liability in connection with this paper, including, but not limited to, the completeness, accuracy, reliability of information contained herein and the potential legal, compliance, tax or accounting effects thereof. Copyright HSBC Group 2025.
References
Tobias Achterberg. Constreint Integer Programming. PhD thesis, Technische Universitaet Berlin, 2007.
Tobias Achterberg. SCIP: Solving constraint integer programs. Mathematical Programming Computation, 1 (1):1-41, July 2009. ISSN 1867-2957. doi: 00.10077/122332-088000001-
Bryan Andrews, Joseph Ramsey, Ruben sanchez-romero, Jazmin Camchong, and Erich Kummerfeld. Fast Scalable and Accurate Discovery of DAGs Using the Best Order Score Search and Grow-Shrink Trees. Advances in neural information processing systems, 36:63945-63956, May 2024.
Charles Assaad, Emilic Devijver, and Eric Gaussicr. Survey and Evaluation of Causal Discovery Methods for Time Series. Journal of Artificial Intelligence Research, 73:767 819, February 2022. doi: 10. 10.1613/jair. 1.13428.
Laura Ballester, Jesúa López, and Jose M. Pavía. European systemic credit risk transmission using Bayesian networks. Research in International Business and Finance, 65: 65:101914, 2023. TSSN 0275-5319. doi: 10. 1016/j.ribaf.20233 wtt.:baff..22.101133 101914.
Alexis Bellot, Kim Branson, and Mihaela der Schaar. Neural graphical modelling in continuous-timee Consistency guarantees and algorithms. In International Conference on Learning Representations, 2021.
Pietro Belotti, Christian Kirches, Sven Leyffer, Jeff Linderoth, James Luedtke, and Ashutosh Mahajan. Mixed-integer nonlinear optimization. Acta Numerica, 22:1 131, 2013. doi: 10.1017/8006229991300022
Pierre Bonami, Mustafa Kilmç, and Jeff Linderoth. Algorithms and Software for Convex Mixed Integer Nonlinear Programs. In Mized Integer Nonlinear Programming, volume 154. October 2009. ISBN 978-1- 4614-1926-6. doi:
1.
Ralf Borndörfer, Heide Hoppmann, Marika Karbstein, and Niels Lindner. Separation of cycle inequalities in periodic timetabling. Discrete Optinization, 35:100552, February 2020. ISSN 15725286. doi: 10.1016/j. 2019 2.5000052...
William J Cook, David I, Applegate, Robert E. Bixby, and Vasek Chvátal. The Traveling Salesman Problem: A Computational Study. Princeton university press, 2011.
Paul Dagum, Adam Galper, and Eric Horvitz. Dynamic network models for forecasting. In Uncertainty in Artificial Intelligence, pp. 41-48. Elsevier, 1992.
J. Dakin. A tree-search algorithm for mixed integer programming problems. Comput. J., 1965.
Aramayis Dallakyan. On Learning Time Series Summary DAGs: A Frequency Domain Approach. 2023.
Amit Danicly and Shai Shalev-Shwarta Complexity theoretic limitations on learning dnf's. In Conference on Learning Theory, pp. 815 830. PMLR, 2016.
G. Dantzig, R. Fulkerson, and S. Johnson. Solution of Large-Scale of the Operations Research Society of America, 2(4):393-414 1954. doi: 10.1287 /0re887/prr.2.4.3333
Problem. Journal
Thomas Dean and Keiji Kanazawa. A model for reasoning about persistence and causation. Computational Intelligence, 5(2):142-1535 150, 1989. doi: 10.1111/j.1467-8661 01011/446660000099000000000004
Alberto Del Pia, Santanu Dey, and Marco Molinaro. Mixed-integer Quadratic Programming in NP.
Mathematical Programming, 162, July 2014. doi: 10.1007/s10107-000 1036--0
Selva Demiralp and Kevin Hoover. Searching for the Causal Structure of a Vector Autoregression* Ozford
Bulletin of Economics and Statistics, 65:745-761 December 2003. doi: 10.1046/1.0305-9049 2003..0087...
14

---

Published in Transactions on Machine Learning Research (11/2025)
S. Even and G. Even. Graph Algorithms. Cambridge University Press, 2011. ISBN 978-1-139-50415-
Tian Gao, Debarun Bhattacharjya, Elliot Nelson, Miaoyuan Liu, and Yue Yu. IDYNO: LearnMalinsky- ing Nonparametrie DAGs from Interventional Dynamic Data. In International Conference on Machine Learning, 2022.
GNU Project. GNU linear programming kit (GLPK), http://wwww...og/o/otoro//////pa 2025.
Gurobi Optimization, LLC. Gurobi optimizer reference manual, 2024.
Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O. Hoyer. Estimation of a Structural Vector Autoregression Model Using Non-Ganssianity. Journal of Machine Learning Research, 11(56):1709-1731, 2010.
IBM. IBM ILOG CPLEX optinization studio. studio, 2025.
Tommi Jaakkola, David Sontag, Amir Globerson, and Marina Meila. Learning Bayesian Network Structure using LP Relaxations. In Proceedings of the Thirtennth International Conference on Artificial Intelligence and Statistics, pp. 358-365. JMLR Workshop and Conference Procoedings, March 2010.
Marcus Kaiser and Maksim Sipos. Unsuitability of NOTEARS for Causal Graph Discovery when Dealing with Dimensional Quantities. Newral Processing Letters, 54(33:1587-1595) June 2022. ISSN 1370-4621, 1573-773X. doi:
Lutz Kilian. Structural Vector Autoregressions. CEPR Discussion Papers 8515, C.E.P.R. Discussion Papers, August 2011.
Jan Kronqvist, Andreas Lundell, and Tapio Westerlund. The extended supporting hyperplane algorithm for convex mixed-integer nonlinear programming. Journal of Global Optimization, 64, June 2015. doi: 10.1007/110890--155-02222
Vyacheslav Kungurtsev, Fadwa. Idlahcen, Petr Rysavy, Pavel Rytir, and Ales Wodecki. Learning Dynamic Bayesian Networks from Data: Foundations, First Principles and Numerical Comparisons, 2024.
Lars Lorch, Jonas Rothfuss, Bernhard Scholkopf, and Andreas Krause. DiBS: Differentiable Bayesian Struc- ture Learning. ArXiv, abs/2105.1188 2021.
Daniel Malinsky and Peter Spirtes. Causal Structure Learning from Multivariate Time Series in Settings with Unmeasured Confounding. In Thuc Duy Le, Kun Zhang, Emre Kıeman, Aapo Hyvärinen, and Lin Liu (eds.), Proceedings of 2018 ACM SIGKDD Workshop on Causal Disocvery, volume 92 of Proceedings of Machine Learning Research, pp. 23 47. PMLR, August 2018
Hasan Manzour, Simge Küggikyavuz, Hao-Hsiang Wu, and Ali Shojaie. Integer programming for learning directed acyclic graphs from continuous data. INFORMS Journal on Optimization, 3(1):46 73, 2021.
Tom Michoel and Jitao David Zhang, Causal inference in drug discovery and development. Drug Discovery Today, 28(1)):103737 2023. ISSN 1359-6446. doi:
2023.103337
Kevin Patrick Murphy. Dynamic Bayesian Networks: Representation, Inference and Learning. PhD thesis, University of California, Berkeley, 2002.
Roxana. Pamfil, Nisara Sriwattanaworachai, Shaan Desai, Philip Pilgerstorfer, Paul Beaumont, Konstantinos Georgatzis, and Bryon Aragam. DYNOTEARS: Structure Learning from Time-Series Data. In Interna- tional Conference on Artificial Intelligence and Statistics, 2020.
Jonas Peters and Peter Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101, May 2012. doi: 10.1093/biomee/
Alexander G. Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the Simulated DAG! Causal Discovery Benchmarks May Be Easy To Game. arXiv, November 2021. doi: 10.45550/arx1v.21122 13647.
15

---

Published in Transactions on Machine Learning Research (11/2025)
Jakob Runge. Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. March 2020. doi: 2000:5500000i.22000.0000
Petr Ryšavý, Xiaoyu He, and Jakub Mareček. Causal Learning in Biomedical Applications: A Benchmark, 2024.
Xiangyu Sun, Guiliang Liu, Pascal Poupart, and Oliver Schulte. NTS-NOTEARS: Learning Nonparametric
Temporal DAGs With Time-Series Data and Prior Knowledge. ArXio, abs/ abs/2100.04288,, 2021.
Peter W G Tennant, Eleanor J Murray, Kellyn F Arnold, Laurie Berrie, Matthew P Fox, Sarah C Gadd, Wendy J Harrison, Claire Keeble, Lynsie R Ranker, Johannes Textor, Georgia D Tomova, Mark S Gilthorpe, and George T H Ellison. Use of directed acyclic graphs (DAGs) to identify confounders in applied health research: Review and recommendations. International Journal of Epidemiology, 50(2); 020-632, December 2020. ISSN 0300-5771. doi: 10.1093 /100/ijc/dyaa213..
Marcel A. J. van Gerven, Babs G. Taal, and Peter J. F. Lucas. Dynamic Bayesian networks as prognostic models for clinical patient management. Journal of Biomedical Informatics, 41(4):515 529, 2008. ISSN 1532-0464. doi:
0..000
Tapio Westerlund and Frank Pettersson. An extended cutting plane method for solving convex MINLP prob- lems. Computers & Chemical Engineering, 19:131-136, 1995. ISSN 1394-1354. doi: 10.1016/0098 1354(95)
87027-X.
Tong Xu, Armeen Taeb, Simge Küççükyavuz, and Ali Shojaie. Integer Programming for Learning Directed Acyclic Graphs from Non-identifiable Gaussian Models. August 2024. doi: 10.48550/ar 20.555/42XX...244112533
Alessandro Zandonå, Rosario Vasta, Adriano Chió, and Barbara Di Camillo. A Dynamic Bayesian Network model for the simulation of Amyotrophic Lateral Sclerosis progression. BMC Bioinformatics, 20, April 2019. doi: 10.1186/s1 0158///11859---269999
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Continuous Optimization for Structure Learning. November 2018. doi: 48r50//aXivv... 180..0142..
Q Zhong, Y. Cheng, 7. Li, D. Wang, C. Rao, Y. Jiang, T.. Li, 7. Wang, P. Liu, Y. Zhao, P. Li, J. Suo, Q. Dai, and K. He. Ultra-efficient Causal Learning for Dynamic CSA-AKI Detection Using Minimal Variables. In medRriu, 2023.
16


Reviews:
[{"review_id": "k17QnOrXkq", "reviewer": "Reviewer_vDaK", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- The authors study a very important problem. The fact that the setting captures time dependencies has high practical potential in important real-world applications.\n- The mixed integer programming formulation is interesting, because it allows the use of globally convergent algorithms. As opposed to formulations in the prior literature, the setting under study involves a directed intra-dependency graph, which makes the problem more challenging.\n- The idea the authors propose of lazily adding constraints during branch-and-cut is interesting, as it enables a practical globally convergent algorithm (even though the actual convergence speed may be slow, depending on the setting).\n- The authors include synthetic experiments where the ground truth is known, so that the performance metrics can be accurately tracked and analyzed. The results seem to suggest that the proposed framework can have a positive impact for small- and medium-sized datasets (e..g, n=50 or n=100), especially as the number of nodes increases. On the other hand, for n>=250 the benefit seems to be negligible or non-existent. Furthermore, in the biomedical example, ExDBN can in principle find the global optimum, as opposed to DYNOTEARS that cannot even reach an F1 score of 0.5. \n\nWeaknesses\n- The authors do not explain why the third variant of lazily adding constraints is better than the other two, and do not provide any results at all.\n- I have concerns with regard to the competitors. The authors compare against DYNOTEARS, a method from 2020, and lingam, a method from 2010. I understand the argument that they do not need to include methods that DYNOTEARS outperforms. But what about more recent methods? The authors have cited more recent works in the field, but for some reason they did not include any comparison. I feel that this point needs clarifications.\n- Figure 1 seems to indeed suggest that ExDBN outperforms the studied competitors for small- and medium-sized synthetic datasets, especially as the number of nodes increases. On the other hand, for n=250 and n=500, we see that ExDBN can occasionally perform relatively poorly (e.g., in SF-3-1 and ER-2-2). But the biggest concern was probably with Figure 2, since the F1 score does not closely track the SHD score and the results are not fully consistent. For instance, for n=50 and ER-2-2, ExDBN is generally outperformed by DYNOTEARS in Figure 2, whereas the opposite is True in Figure 1. The two scores seems not to be so well aligned. And in Figure 2 the performance of ExDBN is not always better (look e.g., at AR-2-2 across all values of n).\n- ExDBN seems to suffer from high variance, compared to both lingam and DYNOTEARS. In some cases, the variance is especially pronounced, e.g., ER-2-2 for n=50. \n- The bioscience example provided an data point (DYNOTEARS cannot even reach an F1 score of 0.5, as opposed to ExDBN that can find a solution with 38% duality gap). I'm not sure there was such a comparison for the finance dataset. For instance, what solution did the procedure of Ballester et al. (2023) converge to?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work proposes a new mixed-integer quadratic program formulation for causal learning from data. To capture dynamic effects, the authors assume a Dynamic Bayesian network setting, where a structural equation model contains both intra-dependencies to the present as well as inter-dependencies on past data. The standard way of solving such a problem would incur exponentially many acyclicity constraints, making the method impractical. Instead, the authors avoid pre-generating all these constraints by lazily adding constraints in a branch-and-cut manner. The authors compare their approach to DYNOTEARS and lingam under various settings, and show that the proposed method produces more accurate results with small- and medium-sized synthetic instances containing up to 25 time series. Furthermore, they demonstrate the benefits of their method, particularly in terms of global optimum convergence, in two applications in bioscience and finance.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The proposed method is interesting, but I feel the paper would benefit from various revisions:\n- The authors could explain why the third variant of lazily adding constraints is better by conducting an ablation study. Currently, it looks as if they did experiments, but do not provide any details at all. Furthermore, it seems that their finding may contradict earlier findings, e.g., (Achterberg 2007, Chapter 8.9). If so, more details are necessary.\n- I'm not so sure that having only two competitors from 2010 and 2020 is sufficient proof about the superiority of the proposed method. The authors should either include more recent methods, or at least explain why they have not included more recent methods (e.g., due to scope, restrictions of more recent methods, etc.).\n- Why are Figures 1 and 2 not better aligned? As I argued above, the results on the two metrics are not totally consistent across all datasets. For example, in Figure 2 it is not the case that as the number of nodes in the DAG increases, the performance gap between ExDBN and DYNOTEARS increases, too (for n=20 or n=100). Furthermore, Figure 2 does not convincingly show that ExDBN is better than DYNOTEARS. The authors could elaborate on all these aspects in the revised manuscript. The results do not necessarily look very convincing, especially in Figure 2, which is further exacerbated by the fact that the authors have only included two old competitors.\n- ExDBN can often suffer from significant variance. Why is that the case? Can the authors elaborate on that, and argue why this is not a problem in practice?\n- I was not really clear about the purpose of the finance experiment. The authors do not provide any result on the solution quality of other methods. Is the point to simply argue that ExDBN can converge to the global optimum?", "broader_impact_concerns": "I do not have broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "7RelYwkd7h", "reviewer": "Reviewer_Tk3c", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n1. The paper propose a mixed integer progamming and avoid exponentially many acyclicity constraints by utilizing the so\"lazy constraint\".\n\n2. The paper includes real‐data case studies and reports improvements in both runtime and structural accuracy over baselines.\n\n**Weakness**\n\n1. There are numerous presentation/professionalism issues (e.g., typos, imprecise notation); see Requested Changes for a detailed list.\n\n2. Constraint (14) is unclear: as written, it seems to require selecting $|C|-1$ edges from a DAG $C$, which is unintuitive and, to me, likely incorrect—please clarify the intention and provide a justification or correction.\n\n3. The motivation for focusing specifically on Dynamic Bayesian Networks is not well articulated. The proposed MIQP/constraint trick appears applicable to static Bayesian network (DAG) structure learning as well; please clarify what is genuinely unique to the dynamic setting (or provide static-BN results to delineate the contribution).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a score-based framework—ExDBN—for learning Dynamic Bayesian Networks by jointly estimating intra-slice DAGs and inter-slice (lagged) dependencies via a mixed-integer quadratic program with sparsity or shrinkage regularization, enforcing acyclicity only within slices and certifying solution quality through the MIP gap. Algorithmically, it advances branch-and-cut for dynamic graphs by adding lazy cycle constraints on demand (with an empirical recommendation to add all detected cycles per iteration), which avoids exponential pre-enumeration while retaining global optimality. Empirically, across synthetic ER and scale-free benchmarks and multiple AR orders, ExDBN achieves lower SHD and higher F1 than strong baselines such as DYNOTEARS and LiNGAM, particularly in small-sample or heteroskedastic settings and at larger graph sizes. Two case studies—a biochemical Krebs-cycle network and a systemic-risk application using European CDS spreads—demonstrate practical viability on real data, including instances solved to global optimality.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Professionalism & presentation issues (requested fixes).**\n\n1. Abstract format. Make the abstract a single, self-contained paragraph.\n2. References. Several entries appear manually typed and contain errors; please generate the bibliography from a .bib file to ensure consistent author names, titles, venues, capitalization, and dates. For example, the following contain obvious issues and should be corrected using canonical metadata:\n– Tian Gao, Debarun Bhattacharjya, Elliot Nelson, Miaoyuan Liu, and Yue Yu. “Idyno: Learnmalinskying nonparametric dags from interventional dynamic data.” In ICML, 2022.\n– Jakob Runge. “Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets,” 03/2020.\n\n3. Citations on first mention. Add a citation immediately after the first mention of each solver/algorithm (e.g., SCIP, GLPK, Gurobi, CPLEX) and for standard procedures like DFS.\n4. Equation referencing style. Refer to numbered displays with parentheses: write Equation (7) or Inequality (7), not “Equation 7.” Use a single, consistent style throughout.\n5. Equation typesetting. Avoid splitting equations into multiple lines with repeated “==” or “++”. \n6. Redundant dating in prose. Replace formulations like “In 2020, Pamfil et al. (2020) have developed…” with a concise form, e.g., “Pamfil et al. (2020) developed…” or “DYNOTEARS (Pamfil et al., 2020) is…”. Avoid leading with the year and do not repeat it in parentheses.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "MjaDAJwEkh", "reviewer": "Reviewer_ghVQ", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n1. The proposed algorithm addresses the curse of dimensionality, which is a significant and relevant challenge both theoretically and practically.\n2. The algorithm has been evaluated through a comprehensive set of experiments, including both synthetic and real-world datasets.\n3. The method appears to be consistent.\n\n**Weaknesses**\n\n1. Section 3, which covers the methodology, should be a core part of the paper but currently lacks sufficient explanation and detail.\n2. Although the algorithm is intended to tackle the curse of dimensionality, it is only applied to small- and medium-sized causal graphs (up to 25 nodes), raising concerns about its scalability.\n3. The experimental evaluation is limited, as only DYNOTEARS is used as a baseline. In the biomedical application, comparing the F1 score of DYNOTEARS with the duality gap of the proposed method is not straightforward. Similarly, in the finance application, there is neither a clear performance comparison metric, nor an explanation of the real-world implications if the ground truth is unavailable.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a score-based algorithm for learning directed acyclic graphs with time lags using a mixed-integer quadratic programming formulation. By employing a Branch-and-Bound-and-Cut strategy, the algorithm effectively mitigates the curse of dimensionality by avoiding the need to regenerate exponentially many acyclicity constraints.\n\nEmpirical results on synthetic data show that the proposed method outperforms DYNOTEARS in certain scenarios, particularly with small sample sizes and a moderate number of nodes, while remaining competitive with DYNOTEARS on larger datasets.  The algorithm has also been applied to two real-world datasets.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. It would be helpful to provide a clearer explanation of the Branch-and-Bound-and-Cut approach, similar to the level of detail given in other parts of the methodology section.\n\n2. Including additional baselines would strengthen the empirical evaluation. While the practical limitation of benchmarking against all existing solvers is understandable, the use of small sample sizes and graph scales may help mitigate this constraint, allowing for the inclusion of multiple representative baselines. Additionally, there are some baseline methods in the same line of work as DYNOTEARS that are not computationally intensive. Please consider reviewing and incorporating such methods as additional baselines.\n\n3. Would it be feasible to report the number of constraints generated or applied by the Branch-and-Bound-and-Cut approach? Doing so could provide a more concrete illustration of how the proposed algorithm addresses the curse of dimensionality.\n\n4. For real-world datasets, please include comparisons using the same evaluation metric if ground truth is available. If not, please clarify whether the estimated DAG is consistent with domain knowledge and explain its practical implications. Figure 4 is not informative without contextual analysis.\n\n5. Please clarify why the proposed algorithm is only applied to graphs with up to 25 nodes. For reference, DYNOTEARS has been applied to graphs with up to 100 nodes.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)