Distributionally Robust Skeleton Learning of Discrete Bayesian Networks
Authors: Yeshu Li
Brian D Ziebart
Venue: NeurIPS 2023 spotlight (2023)
Keywords: structure learning
Bayesian network
robustness
TL;DR: We propose a distributionally robust method for skeleton learning of discrete Bayesian networks with tractable algorithms and out-of-sample exact recovery guarantees.

Abstract: We consider the problem of learning the exact skeleton of general discrete Bayesian networks from potentially corrupted data. Building on distributionally robust optimization and a regression approach, we propose to optimize the most adverse risk over a family of distributions within bounded Wasserstein distance or KL divergence to the empirical distribution. The worst-case risk accounts for the effect of outliers. The proposed approach applies for general categorical random variables without assuming faithfulness, an ordinal relationship or a specific form of conditional distribution. We present efficient algorithms and show the proposed methods are closely related to the standard regularized regression approach. Under mild assumptions, we derive non-asymptotic guarantees for successful structure learning with logarithmic sample complexities for bounded-degree graphs. Numerical study on synthetic and real datasets validates the effectiveness of our method.

Distributionally Robust Skeleton Learning of Discrete Bayesian Networks
Yeshu Li
Brian D. Ziebart
Alibaba Group
Department of Computer Science University of Illinois Chicago
liyeshu. byyeshu.lys@al..1 ibaba- inc. com
bziebartQuic.ddi
Abstract
We consider the problem of learning the exact skeleton of general discrete Bayesian networks from potentially corrupted data. Building on distributionally robust op- timization and a regression approach, we propose to optimize the most adverse risk over family of distributions within bounded Wasserstein distance KL divergence to the empirical distribution. The worst-case risk accounts for the effect of outliers. The proposed approach applies for general categorical random variables without assuming faithfulness, an ordinal relationship or specific form of condi- tional distribution. We present efficient algorithms and show the proposed methods are closely related to the standard regularized regression approach. Under mild assumptions, we derive con-asymptotic guarantees for successful structure learning with logarithmic sample complexities for bounded-degree graphs. Numerical study on synthetic and real datasets validates the effectiveness of method.
Introduction
A Bayesian network is prominent class of probabilistic graphical models that encodes the conditional dependencies among variables with a directed acyclic graph (DAG). It provides a mathematical framework for formally understanding the interaction among variables of interest, together with computationally attractive factorization for modeling multivariate distributions. If we impose causal relationships the edges between variables, the model becomes causal Bayesian network that encodes the more informative causation. Without such interpretation, a Bayesian network serves as dependency graph for factorization of multivariate distribution. We focus on discrete Bayesian networks with purely categorical random variables that are not ordinal, but will discuss related work on both discrete and continuous Bayesian networks for completeness.
The DAG structure a Bayesian network is typically unknown practice [Natori et al., 2017 7 Kitson et al. 2023 Structure learning is therefore an important task that infers the structure from data. The score-based approach defines a scoring function that measures the goodness-of-oo of each structure and aims to find optimal DAG that maximizes the score. Unfortunately. the resulting combinatorial optimization problem is known to be NP-hard [Chickering et al., 2004] without distributional assumptions. Representative approaches include those based heuristic search (Chickering. 2002), dynamic programming [Silander and Myllymäki, 2006), integer linear programming [Jaakkola al. 2010) or continuous optimization [Zheng
2018 which either
yields an approximate solution or an exact solution in worst-case exponential time. The construint- based approach [Spirtes and Glymour 1991 Spirtes al.) 1999 Colombo et al. 2014 performs conditional independence tests to determine the existence and directionality of edges. The time complexity is, however, exponential with the maximum in-degree. Furthermore, the independence test results may be unreliable or inconsistent with the true distribution because of finite samples or
* "ork done when Yeshu was a PhD student at UIC.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

even corrupted samples. In general, without interventional data or assumptions on the underlying distribution, we can only identify a Markov equivalence class (MEC) the true DAG belongs to from observational data where DAGs in the MEC are Markov equivalent, that is, encoding the same set of conditional independencies.
A super-structure is an undirected graph that contains the skeleton as a subgraph which removes directionality from the true DAG, has been shown that given super-structure possibly reduces the search space or the number of independence tests to be performed. For example, exact structure learning of Bayesian networks may be (xredpparameter tractable (Downey and Fellows [995] if the super-structure satisfies certain graph-theoretic properties such as bounded tree-width [Korhonen and Parviainen 2013 Loh and Bühlmann [2014). bounded maximum degree [Ordyniak and Szeider
2013] and the feedback edge number [Ganian and Korchemnal 202T]. An incomplete super-structure with missing edges also helps improve the learned DAG with post-processing hill-climbing method (Tsamardinos et al., 2006) Perrier et al. 2008). Furthermore, a combination of skeleton and variable ordering determines a unique DAG structure. Learning the exact skeleton rather than rough super-structure is desirable in Bayesian network structure learning.
Spirtes and Glymour [199T], Tsaaardnnos et al.
[006) make use of independence tests estimate the skeleton. [Loh and Bühlmann [014] learn a super-structure called moralized graph via graphical lasso [Friedman et al., 20088. Shojaie and Michailidis
(2010) learn the skeleton assuming an ordering of variables. Bank and Honorio [2020] leverage linear regression for skeleton recovery in polynomial time. These methods either rely on independence test results, which are unstable, or a regularized empirical risk minimization problem, where regularization is usually heuristically chosen to combat overfitting. In practice, the observational data is commonly contaminated by sensor failure, transmission error or adversarial perturbation (Lorch et al., 2022, Sannararaman et al.) [202]] Kitson (ct al. 2023]. Sometimes only a small amount of data is available for learning. As a result, the existing algorithms are vulnerable to such distributional uncertainty and may produce false or missing edges in the estimated skeleton.
In this paper, we propose a distributionally robust optimization (DRO) method [Rahimian and MMerotrra 2019] that solves node-wise multivariate regression problem (Bank and Honorio| [202] [0200 for skeleton learning of general discrete Bayesian networks to overcome the above limitations. We do not assume any specific form of conditional distributions. We take into account the settings with a small sample size and potential perturbations, which makes the true data generating distribution highly uncertain. Our method explicitly models the uncertainty by constructing an ambiguity set of distributions characterized by certain priori properties of the true distribution. The optimal parameter is learned by minimizing the worst-case expected loss over all the distributions within the ambiguity set so that it performs uniformly well on all the considered distributions. The ambiguity set usually defined in such a way that it includes all the distributions close to the empirical distribution in terms of some divergence. With an appropriately chosen divergence measure, the set contains the true distribution with high probability. Hence the worst-case risk can be interpreted as an upper confidence bound of the true risk. The fact that a discrete Bayesian network encompasses an exponential number of states may pose a challenge to solve the DRO problem. We develop efficient algorithms for problems with ambiguity sets defined by Wasserstein distances and Hellbackk-Leibblee (KL) divergences. We show that a group regularized regression method is a special case of our approach. We study statistical guarantees of the proposed estimators such as sample complexities. Experimental results on synthetic and real-world datasets contaminated by various perturbations validate the superior performance of the proposed methods.
1.1 Related Work
Bayesian networks have been widely adopted in number of applications such as gene regulatory networks [Werhli et al., 2006). medical decision making [Kyrimi et al., 2020 and spam filtering MMajushaa and Kumar 2010
In addition to the score- based structure learning methods and constraint- constraint-based methods discussed in the introduction section, there are a third class of hybrid algorithms leveraging constraint based methods to restrict the search space of a score-based method [Tsamardinos et al., 2006 Gasse et al. 2014| Nandy et al. 2018). There is also flurry of work on score-based methods based on neural
networks and continuous optimization (Zheng et al. 2018 Wei et al. 2020 Ng et al. 2020 Yu et al.,
2 Ng et al. 2022 Gao al. 2022 motivated by differentiable characterization of acyclicity

---

without rigorous theoretical guarantees. We refer the interested readers to survey papers (Drton and Maathuis 2017 Heinze-Deml et al. 2018 Constantinou et al. 202T for more thorough introduction
of DAG structure learning and causal discovery methods.
Recently, there is emerging line of work proposing polynomial-time algorithms for DAG learning
PPark and Raskutti 2017 Ghoshal and Honorio, 2017 2018 Chen et al. 2019. Bank and Honorio
2020 Gao et al. 2020 Rajendran et al., [021], among which Bank and Honorio [2020] particularly
Tocuses on general discrete Bayesian networks without resorting to independence tests.
Learning super-structure can be done by independence tests, graphical lasso regression, as discussed in introduction. Given a super-structurs how to determine the orientation has been studied
by Perrier et al. 20088, Ordyniak and Szeider [2013], Korhonen and Parviainen 2013 Loh and
Bühlmann
2014 Ng et al. 2021
Ganian and Korchemna 2021]
DRO is powerful framework emerging from operations research (Delage and Ye| 2010 Blanchet and Murthy 2019) Shafieezadeh-Abaaeeee Abadeh et al. 2019 Duchi and Namkoong 2019] and has seen effective
applications in many graph learning problems such inverse covariance estimation (Nguyen et al. 2020], graph Laplacian learning [Wang et al. 2021 Markov random field (MRF) parameter learning Fathony et al. 2018), MRF structure learning
2022], graphical lasso learning [Cisneros- Velarde et al.
Li et al.
2022] and causal inference [Bertsimas et al.
2022].
Preliminaries
We introduce necessary background and a baseline method for skeleton learning of Bayesian networks.
2.1 Notations
We refer to [n] as the index set (1. .... For vector æ R" we use 2, for its i-th element and Is for the subset of elements indexed by S = [m] with + [i). For a matrix A E R"xm we use Agj. A,. and A. denote its jit entry, i-th row and j-th column respectively. AST represents the submatrix of A with rows restricted to S and columns restricted to T S [m]] We define a row- partitioned block matrix as A [A ....... Ax]T E R ny xm where A; E R"*xm, The lp norm of a vector x is defined as ||~||p := E x 1)///p with |-1 being the absolute value function. The lp.g norm of matrix A is defined as || p. = (E //, When p q 2, it becomes the Frobenius norm ||-11 F. The operator norm written as ||| A plg :- suppol.l A The block matrix norm is defined || B.p.- :- |1,)14
The inner product of matrices is designated by (A, B- + Tr[ATB] where AT is the transpose of A. Denote by the tensor product operation. With a slight abuse of notation, [S] stands for the cardinality of a set S. We denote by 1 (0) a vector or matrix of all ones (zeros). Given distribution on E, we denote by Ep the expectation under P. The least c-Lipschitz constant of a function f f:E. - R with a metric e: E x E IR is written as (J) == inf Ac( where Af( == () 0: VE, E2 e =\t (S()) f(E2) E2))
2.2 Bayesian Network Skeleton Learning
Let IP be a discrete joint probability distribution on n categorical random variables V == Let G := (V, Emue) be DAG with edge set Eme.. We use X; to represent the i-th random variable or node interchangeably. We call (G,P) a Bayesian network it satisfies the Markov condition, i.e., each variable X, is independent of any subset of its non-descendants condi- tioned on its parents Pa,.. We denote the children of X,. by Ch,, its neighbors by Ne,- == Pa, u Ch, and the complement by Co, == [n] Ne, {r). The joint probability distribution thus be factorized in terms of local conditional distributions:
Xn). X2,.,
P(X) = (P X2.
I P(X. Pa.).
Let Gstel (V, Esee) be the undirected graph that removes directionality from g. Given m samples (az(() 1 drawn i.i.d. from P. the goal of skeleton learning is estimate Gseel from the samples.
We do not assume faithfulness (Spirtes et al. 2000 or any specific parametric form for the conditional distributions. The distribution is faithful to graph if all (conditional) independencies that hold true
3

---

in the distribution are entailed by the graph, which is commonly violated in practice [Uhler et al. 2013 Mabrouk et al., 2014). The unavailability of a true model entails a substitute model. Bank and Honorio 2020) propose such model based on encoding schemes and surrogate parameters.
Assume that each variable X, takes values from finite set C, with cardinality |Cr| 1. For an indexing set 8 S [n], define ps Eiss [c,] 1 and Ps. Eiss [Ci.. The maximum cardinality minus one is defined as Pmax = [a] [ - 1. Let S, == UseNe (P(--1] + 1,... P(g) be indices for Ne, in P|x] and its complement by S == [p|m]] S, - (Pr--1] + 1. .......rjj].. Let C. BP- be encoding mapping with a bounded and countable set B R.. We adopt encoding schemes with B (-1,0, 0, 1) such as dummy encoding and unweighted effects encoding which satisfy a linear independence condition. With little abuse of notation, encoding any X, and denote by {Xxs) E BPs the concatenation of the encoded vectors (es Consider linear structural equation model for each X,: E(x,) = W*T [W ... W , ... w e RP= with W E IRPiXPr is a surrogate parameter matrix and E IR is vector of errors not necessarily independent of other quantities. A natural choice of fixed W is the solution to the following problem given knowledge of the true Bayesian network:
reuse & for
e, where W P
W arg inf EE (l,,) W WTE(X;) W
W; ViE Co.-
Therefore W = O)) with W = E3
is the
optimal solution by the first-order optimality condition assuming that (x(E(X+)ss {())), is invertible. The expression of W... captures the intuitions that neighbor nodes should be highly related to the current node while the interaction among neighbor nodes should be weak for them to be distinguishable. We further assume that the errors are bounded:
Assumption (Bounded error). For the error vector, ||e|| s & and ||r||e]]l>> & μ.
Note that the true distribution does not have to follow a linear structural equation model. Equation (1) only serves as a surrogate model to find technical conditions for successful skeleton learning, which will be discussed in a moment,
The surrogate model under the true distribution indicates that W |2,2 > 0 X₃ Ne,, This suggests a regularized empirical risk minimization (ERM) problem to estimate W*;
1 F. [[((X)) - WT& M||W B,2, 2
W arg inf L(W)= =
(2)
where Ā 0 is a regularization coefficient, the block f2,1 norm is adopted to induce sparsity and P == 1 Σm (g))) stands for the empirical distribution with مς60 being the Dirac point measure (+)) This approach expected to succeed as long as only neighbor nodes have a non-trivial impact on the current node, namely, WV* l|2.2 > 0 X, Ne,.
-i=1
Define the risk of some W under a distribution 1P as
R ( = W (X) = Ex (17,11((,)) -
where (w)-) is the squared loss function. The Hessian of the empirical risk R.m (W) is block diagonal matrix V²R (W) * H X Ip- Re.PeXPoPr, where H = Eg [E,,[E((,))) [E(x)) ReexPr and I,, E RP- is the identity matrix of dimension pr.. Similarly under the true distribution, As a result, H is independent of the surrogate parameters W* thus
H = (x, conditions on the Hessian translate to conditions on a matrix of cross moments of encodings, which only depend on the encoding function E and P.
In order for this baseline method to work, we make the following assumptions.
Assumption 2 (Minimum weight). For each node r, the minimum norm of the true weight matrix 8 for neighbor nodes is lower bounded: mingeNe, | r > B > 0.
Assumption 3 (Positive definiteness of the Hessian). For each node r. H.,.s, > 0, or equivalently. Amin (...).) > A > where (-in denotes the minimum eigenvalue.
'If there are four variables, dummy encoding may adopt {(1,0, 0), (0, 1, 0) (0, 0, 1). (0. 0, 0)} whereas unweighted effects encoding may adopt ((1,0, (0. 1)0), O0, 1), --1,-1, -1)) encoding vectors.

---

Assumption 4 (Mutual incoherence). For each node |H Hs,,s |B,,000 & 1 - a for some 0 a 1.
Assumption guarantees that the influence of neighbor nodes is significant in terms of a non-zero value bounded away from zero, otherwise they will be indistnguuisable from those with zero weight. Assumption ensures that Equation (2) yields unique solution. Assumption 4lis is a widely adopted assumption that controls the impact non-neighbor nodes on r WWinwwiihh 2009 Ravikumar et al.) 2010 Daneshmand et al., 2014]. One interpretation is that the rows of Hsps, should be nearly orthogonal to the rows of S.s., Bank and Honorio (2200) show that these assumptions hold for common encoding schemes and umite-sample settings with high probability under mild conditions.
They also show that incoherence is more commonly satisfied for the neighbors than the Markov blanket, which justifies the significance of skeleton learning.
Finally, we take the union of all the learned neighbor nodes for each [n] by solving Equation to get the estimated skeleton 9 (V, Esta))
3 Method
As noted in Bank and Honorio [2020), due to model misspecification, even in the infinite sample setting, there is possible discrepancy between the ERM minimizer W and the true solution W* resulting false or missing edges. In the nigh-dimensional setting (m < n) or the adversarial setting. this issue becomes more serious due to limited knowledge about the data-generating mechanism P.
In this section, attempt to leverage a DRO framework incorporate distributional uncertainty into the estimation process. We present efficient algorithms and study the theoretical guarantees of methods. All technical proofs are deferred to the supplementary materials.
3.1 Basic Formulation
Let X a measurable space of all states of the Bayesian network (G. i.e. X X. Let P(X) be the space of all Borel probability measures on X. Denote by X² := (E( VX X3 the space of all the allowed encodings.
Instead of minimizing the empirical risk and relying on regularization, we seek a distributionally robust estimator that optimizes the worst-case risk over an ambiguity set of distributions:
W E arg inf sup 1 -oll W &
W QEA
where A P(X') is an ambiguity set typically defined by nominal probability measure IP equipped with discrepancy measure div(..-) -) for two distributions A = tQ E P(x) div P)) & e). where e is known as the ambiguity radius or size. This way of uncertainty quantification can be interpreted as adversary that captures out-of-sample effect by making perturbations on samples within some budget e. Some common statistical distances satisfy div(Q,
<> Q P.. In
this case, if E is set to zero, Equation 3 reduces to Equation (2) without regularization. We will show that the DRO estimator W can be found efficiently and encompasses attractive statistical properties with a judicious choice of A.
3.2 Wasserstein DRO
Wasserstein distances or Sectorovicc-Ruuiinsee metric in optimal transport theory can be inter- preted as the cost of the optimal transport plan to move the mass from P to Q with unit transport cost c X' x X Rt. Denote by (p')')) the space of all P e P(x') with finite p-th mo- ments for p > 1. Let M(x²2) be the set of probability measures on the product space X' x X. The p- D-Wasserstein distance between two distributions P. Q E Pp(X') is defined W, Q)
[Lre ad)
inf
:T((xx, X) P(dæ), ((x,dx))) Q(da))
We adopt the Wasserstein distance of order p as the discrepancy measure, the empirical distribution as the nominal distribution, and cost function c(x, x') [|E(z) E(x')] for some norm
5

---

||-Il- The primal DRO formulation becomes
1 BollE(X,) WT&((,))
W arg inf
sup Qe " (Pm) 2
According to Blanchet and Murthy 2019 the dual problem of Equation can be written as
inf TE + W0y00
sup 2E.
[((,,) - WTE(==)12 - - E(x
(5)
Strong duality holds according to Theorem in Gao and Kleywegt [022]. The inner supremum problems can be solved independently for each ,( Henceforth, we focus on solving it for some i [m]]
sup sEA 2 - >1(E(=) -(((=))((x())))
(6)
Equation 6 supremum of | convex functions W, thus convex. Since X is discrete set consisting of factorial number of points (m))) unlike the regression problem with continuous random variables in Chen and Paschalidis [2018]. we may not simplify Equation (6) into a regular- ization form by leveraging convex conjugate functions because X' is non-convex and not equal to RRP)) Moreover, since changing the value of x, for some E F is equivalent to changing WTE(()) by vector, unlike Li et al. [2022] where only set of discrete labels rather than encodings are dealt with, there may not be greedy algorithm based on sufficient statisties to find the optimal solution to Equation (6). In fact, let the norm be the & norm, we can rewrite Equation 6 by fixing the values of ||E(x) E(x(())]
|l((,,) WT&(x,)) -
sup mr x S.x.DeKRpp... {(=00)
If we fix k, Equation is generalization of the 0-1 quadratic programming problem, which can be
transformed into a maximizing quadratic programming (MAXQP) problem. result, Equation (6)
is an NP-hard problem with proof presented in Proposition 11 in appendix. [Charikar and Wirth [2004]
develop an algorithm to find an 12(1/ log n)) solution based on semi-detinite programming (SDP)
and sampling for the MAXQP problem. Instead of adopting similar SDP algorithm with quadratic constraints, we propose a random and greedy algorithm to approximate the optimal solution, which is illustrated in Algorithm in appendix, whose per- per-iteration time complexity is (n O(n2mpuax)) mmpmax.. It follows a simple idea that for a random node order π. we select a partial optimal solution sequentially from Ti to Tae We enumerate the possible states of the first node reduce uncertainty. In practice, we find that this algorithm always finds the exact solution that is NP-hard to find for random data with & 12 and Pmax in most cases.
Since X non-convex and not equal to RR using convex conjugate functions will not yield exact equivalence between Equation (5) and a regularized ERM problem. However, can draw such connection by imposing constraints on the dual variables as shown by the following proposition:
Proposition 5 (Regularization Equivalence). Ler W = [W; I 1 R K with Wr- - Ipr
If 2 P() [| the Wasserstein DRO problem in Equation (5) equivalent
1 inf P... Pm 1 -1(())) -TT(()+] [P(M) W
which subsumes linear regression approach regularized by the Frobenius norm special case.
This suggests that minimizing regularized empirical risk may not be enough to achieve distributional robustness. Note that exact equivalence between DRO and regularized ERM in Chen and Paschalidis (2188) requires x = d
Now we perform non oon-asymttotic analysis on the proposed DRO estimator W. First, we would like to show that the solution to the Wasserstein DRO estimator in Equation (4) is unique so that we refer to an estimator unambiguously. Note that Equation (4) is a convex optimization problem but not necessarily strictly convex. and actually never convex in the high nigh-dimensional setting, However, given

---

sufficient number of samples, the problem becomes strictly convex and yields unique solution with high probability. Second, we show that the correct skeleton Exel can be recovered with high probability. This is achieved by showing that, for each node X,. the estimator has 7ero weights for non-neighbor nodes Co, and has non-zero weights for its neighbors Ne, with high confidence. Before presenting the main results, we note that they are based on several important lemmas.
Lemma 6. Suppose E is separable Banach space and fix Po E P(=') for some E' =' Suppose e: c: Rgo is closed convex, k-positively homogeneous. Suppose f - y/ is mapping in the Lebesgue space of functions with finite first-order moment under Po and upper semi-continuous with finite Lipschitz constant lip (f). Then for all 0, the following inequality holds with probability 1:
Qe (Fo)))) 5 ( (Q(EE)) & (f) [{(E")PP((dE)
Lemma
2021 and allows us obtain upper bound between the worst-case risk and empirical risk. is crucial for the following shite-sample guarantees.
6 follows directly from Cranko et al.
Lemma 7. If Assumption holds, for any Q A Wr (PPm), with high probability, H is positive definite.
Lemma 8. Assumptionn 3 and Assumption f hold, for any Q A (Pmm) and E (0, 11. with high probability,
T SES, (H " < 61- 2
The above two lemmas illustrate that Assumption and Assumption 4 hold in the sample setting. Let the estimated skeleton, neighbor nodes and the complement be :- (V, Essa)) Ne, and Co respectively. We derive the following guarantees for the proposed Wasserstein DRO estimator.
Theorem 9. Given Bayesian network (G, P) of categorical random variables and its skeleton
Gskel = (V, Estal). Assume that the condition [ B,2,1 & B holds for some B >> O associated with optimal Lagrange multiplier XB 0 for W* defined in Equation (I). Suppose that W is
DRO risk minimizer of Equation (4) with a Wasserstein distance of order and an ambiguity radius
Eo/m where m is the number of samples drawn ii.dd from P. Under Assumptions
7 4 if the
number of samples satisfies
ClEa log (n/ō) log P(m)) ProoPaanc PancPiig O(
min(u²,1) 1)
where C only depends on u, A, and if the Lagrange multiplier satisfies
32/44maxx
3
AB
a
(a/(4- (2/(4) 2a)= + 2) Pmax / P(
then for any ò (0, 1], r e [n], with probability at least 1 - 6, the following properties hold:
(a) The optimal estimator W is unique.
(b) All the non-neighbor nodes are excluded: Co, Cor.
(c) All the neighbor nodes are identified: Ne, S Ne,,
(d) The true skeleton is successfully reconstructed: Gatel Gstel-
Proof sketch. The main idea in the proof follows that in the lasso estimator (Wainwright 2009). Based primal-dual witness construction method and Lemma 8 it can be shown that if we control X², a solution constrained to have zero weight for all the non-neighbor nodes is indeed optimal. Furthermore, Lemma 7 implies that there is unique solution given information about the true neighbors. The uniqueness of the aforementioned optimal solution without knowing the true skeleton then verified via convexity and a conjugate formulation of the block f,,1 norm. Hereby we have shown that the optimal solution to Equation (4) is unique and excluding all the non-neighbor nodes. Next, we derive conditions on X for the estimation bias I W | 8,2.00 < 3/22 to hold, which
allows us to recover all the neighbor nodes. In such manner, applying the union bound over all the nodes E [r]] leads to successful exact skeleton discovery with high probability.
7

---

3 follows directly from The results in Theorem encompass some intuitive interpretations. Compared to Theorem 1 in Bank and Honorio [2020] we make more explicit the relationship among m, X and 8. On one hand, the lower bound of B ensures that sparse solution excluding non-neighbor nodes is obtained. large error magnitude expectation j therefore elicits stronger regularization, On the other hand, the upper bound X is imposed to guarantee that all the neighbor nodes are identified with less restriction on W. There is naturally trade-off when choosing B in order to learn the exact skeleton. The sample complexity depends on cardinalities P[n]. confidence level ó, the number of nodes n, the ambiguity level čo and assumptions on errors. The dependence on indicates that higher uncertainty caused by larger error norms demands more samples whereas the dependence on μ results from the lower bound condition on X with respect to The ambiguity level is set to 6o/m based on the observation that obtaining more samples reduces ambiguity of the true distribution. In practice, we find that Eo is usually small thus negligible. Note that the sample complexity is polynomial in n. Furthermore, if we assume that the true graph has a bounded degree of d, we find that O Cleo+log (n/6))log n log hom)o" p is logarithmic with respect n, consistent with mini minlp"" (=.1)
the results in Wainwright 2009].
We introduce constants B and XB in order to find a condition for the statements in Theorem to hold. If there exists W incurring a finite loss, we always find solution W and let B + maxw [W] W1..,2,1 be the maximum norm of all solutions. Imposing || & B is equivalent to the original problem. By Lagrange duality and similar argument for the lasso estimator, there exists a X that finds all the solutions with [[ B. Therefore we have a mapping between and Tp-
3.3 Kullback-L DRO
In addition to optimal transport, d-divergence is also widely used to construct an ambiguity set for DRO problems. We consider the following definition of a special g-divergence called the KI. divergence: D(Q || P) S.x In Q(dw) Q(dx). where Q E P()) is absolutely continuous with respect
P(dz)
to PE P()) and ((da) denotes the Radon-Nikodym derivative. A noteworthy property of ambiguity sets based on the KI. divergence absolute continuity of all the candidate distributions with respect to the empirical distribution. It implies that if the true distribution is absolutely continuous probability distribution, the ambiguity set will never include it. In fact, any other point outside the support of the nominal distribution remains have zero probability. Unlike the Wasserstein metric, the KL divergence does not measure some closeness between two points, nor does it have some measure concentration results. However, we argue that adopting the KL divergence may bring advantages over the Wasserstein distance since the Bayesian network distribution we study is a discrete distribution over purely categorical random variables. Moreover, as illustrated below, adopting the KL divergence leads to better computational efficiency.
P(dar)
Let A A (Pm) be the ambiguity set, the dual formulation of Equation Theorem 4 in Hu and Hongl 2013):
Σ ) m ie(m
- w( + YE,
inf y In
W,y>0
which directly minimizes convex objective. In contrast to the approximate Wasserstein estimator, this KL DRO estimator finds the exact solution to the primal problem by strong duality.
The worst-case risk over a KI. divergence ball can be bounded by variance [Lam 2019), similar to Based on this observation, we derive the following results:
Lipschitz regularization in Lemma (
Theorem 10. Suppose that W a DRO risk minimizer of Equation (4) with the KI. divergence and an ambiguity radius s = En
cm/m. mm Given the same definitions of (G, PP), Gakes B, XBy m. in Theorem
Under Assumptions
4 if the number of samples satisfies
ClEa log (n/8) log P(m)) Pmaspfin) min(u2,11
m O
where C depends a, A while independent of n. and if the Lagrange multiplier satisfies the same condition as in Theorem 9 then for any 8 E (0, 1}, r [n], with probability least 8, the properties (a)-(d) in Theorem hold.
8

---

The sample complexities in Theorem and Theorem T0 differ in the constant C due to the difference between the two probability metrics. Note that C is independent of n in both methods. The dependency on 1/(x 2 is absorbed in the denominator because we require that X - 16/4mmmp The sample complexities provide a perspective of our confidence on upper bounding the true risk in terms of the ambiguity radius. En serves as our initial guess on distributional uncertainty and increases the sample complexity only slightly because it is usually dominated by other terms in practice: < log(n/6). Even though the samples are drawn from an adversarial distribution with proportion of noises, the proposed methods may still succeed as long as the true distribution can be made close to upper confidence bound.
0.
4 Experiments
We conduct experiments on benchmark datasets Scutari 2010 and real-world datasets
[Malone
et al.
2015 perturbed by the following contamination models:
Noisefree model. This is the baseline model without any noises.
Huber's contamination model. In this model, each sample has a fixed probability of Ç to be replaced by sample drawn from an arbitrary distribution.
Independent failure model. Each entry of sample is independently corrupted with probability Ç.
We conduct all experiments on laptop with an Intel Core i7 2.7 GHz processor. We adopt the proposed approaches based on Wasserstein DRO and DRO, the group norm regularization method
Bank and Honorio 2020]. the MMPC algorithm Tsamardinos et al. 2006) and the GRaSP algorithm [Lm et al. [202]] for skeleton learning. Based the learned skeletons, we infer DAG with the hill-climbing (HC) algorithm [samardinos et al. 2006]. For the Wasserstein- based method, we leverage Adam (Kingma and Bal [2144] to optimize the overall objective with 31 0.9, Bg 0.990, a learning rate of 1.0, a batch size of 500, a maximum of 200 iterations for optimization and 10 iterations for approximating the worst-case distribution. For the KL-based and standard regularization methods, we use the L-BFGS-B B [Byrd et al. [1955] optimization method with default parameters. We set the cardinality of the maximum conditional set in MMPC. The Bayesian information criterion (BIC) [Neath and Cavanaugh] 2012] score is adopted in the HC algorithm. A random mixture of 20 random Bayesian networks serves as the adversarial distribution for both contamination models. All hyper-parre parameters are chosen based on the best performance on random Bayesian networks with the same size as the input one. Each experimental result is taken as an average over 10 independent runs. When dealing with real-world datasets, we randomly split the data into two halves for training and testing.
We use the F1-score, or the Dice coefficient (regarding the label of each edge indicating its presence as binary random variable and considering all possible edges), to evaluate performance on benchmark datasets and BIC for real-world datasets. The results are reported in Table [ and more results can be found in Table 2 in appendix. We observe that in most cases the proposed DRO methods comparable to MMPC and MMHC, which are generally the best-performing methods in Bank and Honorio 2020)) We illustrate in Figure [ the results on earthquake by varying the number of samples, corruption level and ambiguity radius or regularization coefficient. Figure (a) suggests that all the methods perfectly recover the true skeleton given more than 2, 000 samples. The results in Figure (b-c) indicate that, in highly uncertain settings, Wasserstein DRO as well as KL DRO is superior to other approaches. Meanwhile, Table I and Figure [ (a) suggest that the DRO methods and the regularized ERM approach are comparable to MMPC and GRaSP when clean data is given. [ suggests a trade-off between robustness and target performance (Fl-score in our case). All the approaches have similar execution time except that Wasserstein DRO is several times slower due to the combinatorial sub-problem of computing the worst-case distribution.
The sensitivity analysis
5 Discussion and Conclusion
In this paper, we put forward a distributionally robust optimization method to recover the skeleton of general discrete Bayesian network. We discussed two specific probability metrics, developed
Our code is publicly available at htttts///ggttabb../////aaleeeel/

---

Table 1: Comparisons of F1 scores for benchmark datasets and BIC for real-world datasets (backache, voting), BIC is not applicable to skeletons. The best and runner-up up results are marked in bold. Significant differences are marked by f (paired t-test, 0.05).
Dessoce
Noise
Wiss
Reg
M0MPC
GRASP
WaveHC
KL+MIC
RegeHCC
MNDC-RC
CRASP-HC
HC
3200
Netafroz
0 po001
0 12221
0.10011
0.3007
0.5167
0.5125
010977
0.3743
0.0507
0.0383
015680
3:00
0.73131
0.72871
0:0/7
0.3915
0.3724
6-3-7
0.2917
113161
0.2383
3200
Nasafrez
1,00001
1.00301
1.00001
0L5333
00050
0.2800
0.2800
0.3000
0.2500
0.2500
0.2800
3000
Flooor
03.01
0.8321
01333
0.587
0-4333
0.363
6.41
0:394
3080
Natoiiis
0 01
82000
0.0776
1,0000
0.9778
0.200
0.2540
0.210
0.2300
0.2500
0.37791
B203
0.7
0155831
0-4018
0.0754)
0.63501
0.0000
0.000
0.0000
0.0000
suciaa
11
2004
Nandfoc
0.00371
0.5450
0.3037
1441001
0140351
0.081
0.0883
0-41121
0.3875
B200
Indeperalint
0.52581
0.52081
0.3121
1L5050
0.0000
0,0020
0.0010
0.0000
0.0000
0.2000
mansey
Nasifine
01635
0771
0.1748
6.178M
0.1799
0.178
wursey
$000
Iedeperdinn
0.5
0.62111
0.67051
018501
0.583
0:007
0.1071
6.1143
0.107)
0.000
0.307
2002
Natetroe
01-77411
0.8622
0.5330
1L38331
0.44TP
0.5633
0.4576
1616
slare
37
2004
Indeperdizi
02
0.14331
013641
0.54381
0.4309
112552
0.0000
0,0000
5.0000
0.0000
DJ0000
0,0000
vozing
17
S16
0.2
NGA
1887 1544
-156-240
1145-0775
-1116.8731
XON
44% 5208
4453 8296
465-007
MDA
NGA
-1710.7246
-1756
1713.7563
-1739 7961
baczeshe
12
10
Bulas
0.2
RGA
XON
3388
D001
5188 5001
$156 3001
5190
5001
5196 D401
3186. 5001
1.4
10-
0.4
Dgg
K2-
02:67
Siepls
Figure 1: Skeleton estimation results in Fc-score on the earthquake dataset for Wasserstein DRO, KL DRO, regularized linear regression, MMPC and GRaSP. From left to right: (a) noisefree, varying number of samples; (b) Huber with Ç = 0.5 noisy data, varying number of samples; (c) independent failure, varying S. the probability of noise; (d) independent failure with 0.2, varying the ambiguity radius, or À the regularization coefficient.
tractable algorithms to compute the estimators. We established the connection between the proposed method and regularization. We derived non-asymptotic bounds polynomial in the number of nodes for successful identification of the true skeleton. The sample complexities become logarithmic for bounded- graphs. Empirical results showcased the effectiveness of our methods.
The strength of making no distributional assumptions in our methods inherits from the regularized regression baseline, which is shown to be a special case of DRO. Besides the original benefits in Bank and Honorio (2020), our methods are explicitly robust and able to adjust to incorporate our uncertainty about the data generating mechanism.
Since we do not make any specific assumptions on the conditional probability distributions or on the corruption models, our methods may not be superior to the approaches proposed to tackle certain family of noises or parametric distributions. In addition, making assumptions such as ordinal relationship. continuous values or a Poisson structural equation model (SEM) may lead to more efficient algorithms and tighter bounds. In some cases, adopting an identity encoding mapping is sufficient to learn continuous Bayesian network [Aragam et al. 2015]. Furthermore, it would be interesting to incorporate prior distributional information into design of the ambiguity set for better performance. Another important topic is whether the underlying structure is identifiable in a robust manner [Sankararaman et al.) 2022]. Formulating the complete DAG learning problem as one optimization problem may lead non-convex problem. Nonetheless, leveraging adversarial training approach continuous optimization for Gaussian Bayesian networks and causal discovery is promising future direction pursue.
Acknowledgments and Disclosure of Funding
This material is based upon work supported by the National Science Foundation under Grant No. 1652530.
10

---

References
Bryon Aragam, Arash Amini, and Qing Zhou. Learning directed acyelic graphs with penalized neighbourhood regression. arXiv preprint arXiv 2015.
Adarsh Bank and Jean Honorio. Provable efficient skeleton learning of encodable discrete bayes nets in poly-time and sample complexity. In 2020 IEEE International Symposium on Information Theory (ISIT), pages
IEEE, 2020.
Dimitris Bertsimas, Kosuke Imai, and Michael Lingzhi Li. Distributionally robust causal inference with observational data. arXiv preprint arXiv:22 0832066 2022.
Jose Blanchet and Karthyek Murthy. Quantifying distributional model risk via optimal transport. Mathematics of Operations Research, 44(2):565-600, 2019.
Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A limited memory algorithm for 1995.
bound constrained optimization. SIAM Journal on scientific computing,
Moses Charikar and Anthony Wirth. Maximizing quadratic programs: Extending grothendieck's inequality. In 45th Annual IEEE Symposium on Foundations of Computer Science, pages 54-60. IEEE, 2004.
Ruidi Chen and Ioannis C Paschalidis. A robust learning approach for regression models based on distributionally robust optimization. Journal of Machine Learming Research, 19(13), 2018.
Wenyu Chen, Mathias Drton, and Y Samuel Wang. On causal discovery with equal-variance assumption. Biometrika, 106(4):973-880, 2019.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 30755):507-554, 2002.
Max Chickering, David Heckerman, and Chris Meek. Lampe learning of bayesian networks is np-hard. Journal of Machine Learning Research, 5:1287-1330, 2004.
Pedro Cisneros- Velarde, Alexander Petersen, and Sang- Yun Oh. Distributionally robust formulation and model selection for the graphical lasso. In International Conference on Artificial Intelligence and Statistics, pages 756-765, PMLR, 2020.
Diego Colombo, Marloes H Maathuis, et al. inderinndee constraint bonstraint-based causal structure learning. J. Mach. Learn. Res., 15(1):37413378, 2014.
Anthony Constantinou, Yang Liu, Kiattikun Chobtham, Zhigao Guo, and Neville K Kitson. Large- scale empirical validation of bayesian network structure learning algorithms with noisy data, International Journal of Approximate Reasoning, 131:151-188, 2021.
Zac Cranko, Zhan Shi, Xinhua Zhang, Richard Nock, and Simon Korbblith. Generalised lipschitz regularisation equals distributional robustness. In International Conference on Machine Learning, pages 2178-2188. PMLR, 2021.
Hadi Daneshmand, Manuel Gomez-Rodriguez, Le Song, and Bernhard Schoelkopf. Estimating diffu- sion network structures: Recovery conditions, sample complexity soft-thresholding algorithm. In International conference machine learning, pages 793-801. PMLR, 2014.
Erick Delage and Yinyu Ye. Distributionally robust optimization under moment uncertainty with application to data-driven problems. Operations research, 58(3):995-61,, 2010.
Rod G Downey and Michael R Fellows. Fixed-parameter tractability and completeness i: Basic results. SIAM Journal on computing, 24(4):873-911, 1995.
Mathias Drton and Marloes H Maathuis. Structure learning in graphical modeling. Annual Review of Statistics and Its Application, 4:365-393, 2017.
John Duchi and Hongseok Namkoong. Variance-based regularization with convex objectives. The Journal of Machine Learning Research, 20(1):2450 2019.

---

Rizal Fathony, Ashkan Rezaei, Mohammad Ali Bashiri, Xinhua Zhang, and Brian Ziebart. Dis- tributionally robust graphical models. Advances in Neural Information Processing Systems, 31, 2018.
Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432-44)) 2008.
Robert Ganian and Viktoriia Korchemna. The complexity of bayesian network learning: Revisiting the superstructure. Advances in Neural Information Processing Systems, 34:430-442, 2021.
Ming Gao, Yi Ding, and Bryon Aragam. A polynomial-tmme algorithm for learning nonparametric causal graphs. Advances in Neural Information Processing Systems, 33:11599-11611 2020.
Rui Gao and Anton Kleywegt. Distributionally robust stochastic optimization with wasserstein distance. Mathematics of Operations Research, 2022.
Tian Gao, Debarun Bhattachariya, Elliot Nelson, Miao Liu, and Yue Yu. Idyno: Learning nonpara- metric dags from interventional dynamic data. In International Conference on Machine Learning. pages 6988-7001. PMLR, 2022.
Maxime Gasse, Alex Aussem, and Haytham Elghazel. A hybrid algorithm for bayesian network structure learning with application to multi-label learning. Expert Systems with Applications, 41 (15):6755-6772, 2014.
Asish Ghoshal and Jean Honorio. Learning identifiable gaussian bayesian networks in polynomial time and sample complexity. Advances in Neural Information Processing Systems, 30, 2017.
Asish Ghoshal and Jean Honorio. Learning linear structural equation models in polynomial time and sample complexity. In International Conference on Artificial Intelligence and Statistics, pages 1475. PMLR, 2018.
Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Statistical learning with sparsity. Mono- graphs on statistics and applied probability, 143:143, 2015.
Christina Heinze-Deml, Marloes H Maathuis, and Nicolai Meinshausen. Causal structure learning. Annual Review of Statistics and Its Application, 5:371-394 2018
Zhaolin Hu and Jeff Hong. Kullback-leibler divergence constrained distributionally robust opti- mization. Available at Optimization Online, 1(2):9, 2013.
Tommi Jaakkola, David Sontag, Amir Globerson, and Marina Meila. Learning bayesian network structure using ip relaxations. In Proceedings of the Thirteenth International Conference Artificial Intelligence and Statistics, pages 358-365.11 JMLR Workshop and Conference Proceedings, 2010.
Leonid Vasilevich Kantorovich and SG Rubinshtein. On a space of totally additive functions. Vestnik of the St. Petersburg University: Mathematics, 13(7):52-59, 1958.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint 14Xiv:1412.6980 2014.
Neville Kenneth Kitson, Anthony Constantinou, Zhigao Guo, Yang Liu, and Kiattikun Chobtham. A survey of bayesian network structure learning. Artificial Intelligence Review, pages 1-94, 2023,
Janne Korhonen and Pekka Parviainen. Exact learning of bounded tree-width bayesian networks. In Artificial Intelligence and Statistics, pages 370-378. PMLR, 2013.
Evangelia Kyrimi, Mariana Raniere Neves, Scott McLachlan, Martin Neil, William Marsh, and Norman Fenton. Medical idioms for clinical bayesian network development. Journal of Biomedical Informatics, 108:103455. 2020.
Henry Lam. Recovering best statistical guarantees via the empirical divergence- distributionally robust optimization. Operations Research, 67(4):1090-1105 2019.
12

---

Wai- Yin Lam, Bryan Andrews, and Joseph Ramsey. Greedy relaxations of the sparsest permutation algorithm. In Uncertainty in Artificial Intelligence, pages 1052-1062. PMLR, 2022.
Yeshu Li, Zhan Shi, Xinhua Zhang, and Brian Ziebart. Distributionally robust structure learning for discrete pairwise markov networks. In International Conference on Artificial Intelligence and Statistics, pages 8997-9016 PMLR, 2022.
Po-Ling Loh and Peter Bühlmann. High-dimensie learning of linear causal networks via inverse covariance estimation. The Journal of Machine Learning Research, 15(1):3065-310 3105, 2014.
Lars Lorch, Scott Sussex, Jonas Rothfuss, Andreas Krause, and Bernhard Schölkopf. Amortized inference for causal structure learning. Advances Neural Information Processing Systems, 35: 13104-13118, 2022.
Ahmed Mabrouk, Christophe Gonzales, Karine Jabet- ChevvlChe,,,,er, and Eric Chojnacki. An efficient bayesian network structure learning algorithm in the presence of deterministic relations. In ECAI 2014, pages 567-572. IOS Press, 2014.
Brandon M Malone, Matti Järvisalo, and Petri Myllymäki. Impact of learning strategies on the quality of bayesian networks: An empirical evaluation. In UAI, pages 562-571. Citeseer, 2015.
Manjusha and Rakesh Kumar. Spam mail classification using combined approach of bayesian and neural network. In 2010 International Conference on Computational Intelligence and Communica- fion Networks, pages 149, IEEE, 2010.
Preetam Nandy, Alain Hauser, and Marloes Maathuis. High-dimensional consistency score-based and hybrid structure learning. The Annals of Statistics, 46(6A)):151-3111 2018.
Kazuki Natori, Masaki Uto, and Maomi Ueno. Consistent learning bayesian networks with thousands of variables. In Advanced Methodologies for Bayesian Networks, pages 57-68. PMLR, 2017.
Andrew Neath and Joseph E Cavanaugh. The bayesian information criterion: background, derivation, and applications. Wiley Interdisciplinary Reviews: Computational Statistics, 4(2): 199-203, 2012.
Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and dag constraints for learning linear dags. Advances in Neural Information Processing Systems, 33:1 13:17943-17954 2020.
Ignavier Ng,
Yujia Zheng, Jiji Zhang, and Kun Zhang. Reliable causal discovery with improved exact search and weaker assumptions. Advances in Neural Information Processing Systems, 34: 10308-20320, 2021.
Ignavier Ng. Shengyu Zhu, Zhuangyan Fang. Haoyang Li, Zhitang Chen, and Jun Wang. Masked gradient-based causal structure learning. In Proceedings of the 2022 SIAM International Conference on Data Mining (SDM), pages 424-432 SIAM, 2022.
Viet Anh Nguyen, Daniel Kuhn, and Peyman Mohajerin Esfahani. Distributionally robust inverse covariance estimation: The wasserstein shrinkage estimator. Operations Research, 70(1):490-515, 2022.
Sebastian Ordyniak and Stefan Szeider. Parameterized complexity results for exact bayesian network structure learning. Journal of Artificial Intelligence Research, 46:263-302, 2013.
Gunwoong Park and Garvesh Raskutti. Learning quadratic variance function (qvf) dag models via overdispersion scoring (ods). J. Mach. Learn. Res,, 18:224-1, 2017.
Eric Perrier, Seiya Imoto, and Satoru Miyano. Finding optimal bayesian network given a super- structure. Journal of Machine Learning Research, 9(10), 2008
Hamed Rahimian and Sanjay Mehrotra. Distributionally robust optimization: A review. arXiv preprint arXiv:1908.05659) 2019.
Goutham Rajendran, Bohdan Kivva, Ming Gao, and Bryon Aragam. Structure learning in polynomial time: Greedy algorithms, bregman information, and exponential families. Advances in Neural Information Processing Systems, 34: 186118666-17222 2021.
13

---

Pradeep Ravikumar, Martin J Wainwright, and John D Lafferty. High-dimensiioal ising model selection using -regullarized logistic regression. The Annals of Statistics, 38(3):1287-1319 2010.
Karthik A Sankararaman, Anand Louis, and Navin Goyal. Robust identifiability in linear structural equation models of causal inference. In Uncertainty in Artificial Intelligence, pages 1728-1737.
PMLR, 2022.
Marco Scutari. Learning bayesian networks with the bnlearn r package. Journal of Statistical Soffwaree 35:1-22, 2010.
Soroosh Baafieezadeh-bbadee Daniel Kuhn, and Peyman Mohajerin Esfahani. Regularization via mass transportation. Journal of Machine Learning Research, 20(103):1-68, 2019.
Ali Shojaie and George Michailidis. Penalized likelihood methods for estimation of sparse high- dimensional directed acyelic graphs. Biometrika, 97(3):51--533 2010.
Tomi Silander and Petri Myllymäki. simple approach for finding the globally optimal bayesian
network structure. In Proceedings of the Twenty-Second Conference Uncertainty in Artificial Intelligence, pages 445-452, 2006.
Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs. Social science computer review, 9(1):62-72, 1991.
Peter Spirtes, Christopher Meek, and Thomas Richardson. An algorithm for causal inference in the presence of latent variables and selection bias. Computation, causation, and discovery, 21: 211-252, 1999.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine learning, 65(1):31-78, 2006.
Caroline Uhler, Garvesh Raskutti, Peter Bühlmann, and Bin Yu. Geometry of the faithfulness assumption in causal inference. The Annals of Statistics, pages 436-463, 2013.
Martin J Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery using { ij-constrained quadratic programming (lasso). IEEE transactions on information theory, 55(5): 2183-2202. 2009.
Xiaolu Wang, Yuen-Man Pun, and Anthony Man-Cho So. Distributionally robust graph learning from smooth signals under moment uncertainty. arXiv e-prints, pages arXiv-2105, 2021.
Dennis Wei, Tian Gao, and Yue Yu. Dags with no fears: A closer look continuous optimization for learning bayesian networks. Advances in Neural Information Processing Systems, 33:3895-3906, 2020.
Adriano V Werhli, Marco Grzegorezyk, and Dirk Husmeier. Comparative evaluation of reverse engineering gene regulatory networks with relevance networks, graphical gaussian models and bayesian networks. Bioinformatics, 22(20):2523-2555
2006.
Yue Yu, Tian Gao, Naiyu Yin, and Qiang Ji. Dags with no curl: An efficient dag structure learning approach. In International Conference on Machine Learning, pages 12156-1221 12166. PMLR, 2021.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. Advances in Nearal Information Processing Systems, 31, 2018
14

---

end for A Algorithms
The pseudo-code of the greedy algorithm for solving Equation 66 in Wasserstein DRO illustrated in Algorithm
Algorithm Greedy Algorithm for the Wasserstein Worst-case Risk
Input: W. (+)
Output: solution di to Equation 6
Initialize ④ = (x)
for all (j,zf) e [n] C, do
Get a random permutation æ over [n] with mı =
for & = 2 to do
5, arg sup, (w - 7 >/(E(()))))) ((
end for if x yields a greater objective than ø then 2 x
end if
B Optimization Details
Define
1 - WT
(w(X) =
The Lagrangian dual problem of the Wasserstein DRO problem is
inf f(W, 7) 75
sup (w(x) y](E(x) )11
W W.000
xe
One of its sub- nub-gradients can be computed as
1 m
à 8(+())(((())))
E(a E(
TW
aw
m
à E(a())) {(x
:y f.
m
For the DRO problem based on the KL divergence:
Σ (w(wl00)/1 + +
inf f(W 7) ==
W
m
xe[m)
sub-gradient of which can be computed
Eiefm (a (0) (E(x ((x) JW - ((
à
Lie(m) (2((111))
(ww(e)l)///) - Exelml cl (w(
In( m v=[m]
+CE f.
Eve(ml elw (c(1))/9
C Technical Proofs
Proposition 11 (NP-hardness of Wasserstein DRO Supremum). The problem in Equation 6 is Na-hard.
15

---

Proof. Recall the MAXQP problem:
oyyty s.t. 24 E (-1,1] Vi.
In Equation 6 let y {x,,) 0, or correspond to n binary variables taking values in (-1, and ((x,) = wy. Let E [p=x=" For all i [n]. let k = (i - 1)m + j. The k-the column of W satisfies Wix 1, Wjx jk /2 and 0 for the other elements. We have obtained polynomial-time reduction from an NP-hard problem to Equation 66.
Proposition 5 (Regularization Equivalence). Ler w == - Ro with W I p..
y P(a]l W|],, the Wasserstein distributionally robust regression problem in Equation (5) is equivalent to
inf p. - eP(n] W WIE
which subsumes a linear regression approach regularized by the Frobenius norm as a special case.
Proof. Recapitulating on Equation c
[((x) - W WT -(11(())))1((())) E(x( me
sup
Observe that
[8(x,)) WTE + | W pm)
WT l2.2
W 12
P0gl
W
E7.
Therefore, for any * ((1)
1
[8(x,) WTE(x;)1B >1(E(z) {( - (1 -TTE(; (3 71(E(x(0)) - {(x(
2
(1E(x,) - WT&((;;)1 18(x!() WT WTE [13) >1(E()) - {(x(0))11
2
2 (27) 7118(z) E(2
$7-7
=0,
which implics that the supremum can always achieved at æ = ( Minimizing over y leads to
inf Fp.. 1 |1((x,) (CT(X))-TT((X)]] | W
Lemma Suppose E is separable Banach space and fix P. P(E') for some E' E Suppose c R is closed convex, k-positively homogeneous. Suppose f - is mapping in the Lebesgue space of functions with finite first-order moment under Po and upper semi-continuons with finite Lipschitz constant lip (f). Then for all 0, the following inequality holds with probability 1:
f(E')((dE)) < Elip.(f) + [(())o(dE))
sup C "2 (Po)..EPPE)))
16

---

Proof. The result follows directly from Theorem in Cranko et al. [002T]
f(E)Q(dE) & elip. + f(E))o(dEe)
sup CE (PA).QEP(E)
Since E' S E, observe
((5))((E)) &
sup ((S)(()))
sup QEA (EPP)
QEA (P..
Lemma 7. If Assumption holds, for any Q A (Pm)) with probability at least 1 S| exp (- 2[S
mt²
have
min H,s,)) > Amin Ami (Hs.s,) 4c|S,
Proof. The minimum eigenvalue of the true covariance matrix Hs,s, satisfies:
Amin (Hsss))
min ||
min w TH s,s, u u (H.., H S. v (H.... H.s.s)u |
Amin s.s + u (Hs,s, H u (Hs,s, Hs.s,) where |/2e|2 1 is an eigenvector of H with minimum eigenvalue.
Therefore, Amin H s.).)) can lower bounded as follows:
Amin (H) S. > Amin Hsss.) uT(Hs.s. H S. )u u uHH..s. Hu.s.Ju
Amin Hs.s.) - lu [HS.S... - H S.s,)u jul [Hs.s. Hs..s)lle,
due the fact that
u Hu & Amux H & (H), < H 2.2 2.2
We can obtain an upper bound on |uT (HH,,,
H S. )ul based on Lemma 6
Ju (Hs,s ) )ul < +
because for function g(E(x)) = u H.s.s u, it can be shown that for any ||E(x) E(x')ll and some |S| k
Ig(E(z)) ((E(('))) <
Σ Σ
H. uu + Has H uu & 4k|S,|
ieS.
Recall that we assume that the encoding schemes take values B = (-1,0,1}. Therefore lip.(g) = 1
We derive an upper bound of [HHsss - H...s)llF as follows. Consider a random variable and its expectation
Zig == (Hs.s,) = m
1
8(x (1 )(E()))) E [-1/m, 11/mm,/mm 1/m]
Ez (Hs.
By Hoeffding's inequality, we observe
Prob(I(Hs.s..e (HS.s.)ol > < exx (- mt?
2 )
17

---

for > 0. Setting t = [S.]
t for all ú,j E S, and applying the union bound,
ml²
Prob || (HHs,s, - (Hs.s,) F > < 2|S, exp (- S|
(8)
To conclude, with probability at least 2|Sr| exp (- mt²
2|8,12 ) we have
min H S. Amin (.. 4c|5,1 t.
Lemma 8. If Assumptiinn 3 and Assumption 4 hold, for any Q E A ((Pm) and a E (0, 1]. with
probability at least (exxp P S
Cm
+ log | + log (rrl)) and e & 87.37
H (HS 11B,1,,0 & a 2'
where C only depends min S.s..).
Proof. We would like to obtain an upper bound for H S:S, S.S.)
(H ) 11,1,000 We may write
H Sas, (H) s.,
Ms S..s, ) (Hs,s,)
[ s,
H S,s. Hs;,, ( S.s.) (Hs,s,)
+ Hs,s, (Hs.s.))
[ s;s, s.s.) (
< Hsps, [ S.S.) (Hs,s,) -1,,s,))1l((B,,,,,1
+ +|[ s.s. H5ys, 1 (Hs,,,) |1,,00
+ | s.s. -Hs S.s,)
+ Hs;s, (Hs,s,)
By Hoeffding's inequality,
Prob([(Hssss)) (Hsps.)isl 2exp (- mt²
22
)
for L > 0. Taking / = P.S.TT
f and applying the union bound over E Co,, we observe that
mt?
Prob Hss.. Hsps, IBI..x > t) <
Σ ioCo.
Sp exp (- |p?|S, S,
mt?
< S
Sr
exp
2p S,
Similarly, taking = [S
ml?2
Prob(]| Hs.s. Hs.s,lla, > 1) & Σ Σ
2exp
2|S
ieSr JeS.
mt?
-2|S, exp
2|5,1
In order to bound
|| |HS.
s.s. Hsss, |,,, for Q o consider
H s. Hsss, 1,,00 < H s.s. 11,0000 + H.s.s. 1...000
Ep))E(())s I + E, ((
sup F. Peed"
|| -
18

---

where Q and P.,, are probability measures on E E = (-1,+11 and identical marginals as Q
and P' respectively. We assume that Q * IP because otherwise || 5;5, Hs,s, 1l. B. 1,00 0 holds
trivially. In this way, the equality is always achieved by some Q', P m. i.e., setting ( - 1) - 1
and P ( -1) 1.
Define the transport cost function in the ambiguity set Ae " to be (( (x), [| 1 with zero cost for &. Let 9 (X {) = the Lipschitz constants of g:
B,1, Consider
1o(X x( :5) c(((X, (x'
(9) $
sup
x X
(( |1B, + (Ix(()))s |B,,,o
sup X,X'
[E(X)- - 1
2m S,
(9)
Therefore, by the Kantorovich-Rubinsee theorem Kantorovich and Rubinshtein, 1958
| ScS. -Hs,s, IB,,, <
sup Pmm T'EA P( m))
(xg(( L g ((x,)) E)]
sup P x
(g) Ba x, ()/ - Ep, ((x))))
lip. (QP
sup W. M'.. " (P.)
(g)e < 28 Pmax Srl.
Similarly,
H Hs.s, 2c|S,
Based on the above two inequalities, we find that
| S.S. Hsss, | | S.S.
H.s.s + |Hsss - Hsss,lBB,1
22 Pmax |
(10)
with probability at least 1 - 2| exp ( 2p S and
mf²
|| S.S. Hs..s 100,00 & 25|S,] + 4
(11)
with probability at least 1 - 2| exp (- 2)S.
Based on Equation (8) we also have
[.... - H sll & 2=|S,|+t,
(12)
with probability at least - 2|S exp
mf? 2|S
Next we look at the upper bound on the difference between the inverses of H H,S,, and Hs.s.. Observe that
( ( (HH,HH 10000000 = (HH,,))).) [HS.S H (HH 1,,000
<vIS, (Hs,s,) [Hs,s, HH S.S.) 2,2
<vIS SS (H
|||, Hs.s. - H S.S, 11ll (H S. -- 2,2
|8,1
[Hs,s, H 2. (H.)s.))) -
Amin (.,ss))
19

---

According to Lemma with probability at least 1 - -, exp
mf² 2|.
we have
Amin (HS, 2 Amin (Hs,s.) - 4e|S, -t. Ami (Hs,
Let 1 Amin s.sss) and
We get that, with probability at least
16|S,13
2|S,|2 exp (- mmAmml s.))
8|S.1
Amia H S.S.) > Amin (Hs.s.)
4
[l( S.
2,2 $ Aminl S.)
(13)
Set t (Ami(Hs, and &
tdmi S.. 8|8.1 Vs.
in Equation (12), get that, with probability least
1
1 - 2 exp (- mt? (Amr(Hs...))) (H)).,))) )
32|5,
tAminl H.S.s))
[|S.,, H H. 11122 [HsS. H H. 1.ll
2/1S,
2/8,11 exp (- me (Aain (Hs 3,)
Therefore,
probability at least 1 m(Amm (Hs s. 8|S.1"
with
32,5.11
2|S,| exp
and & min (H...H.....) 8/S-ly [8.]
Amn Hs..
16/5.1
S.S, )- (Hs. -1ll
o0,,, & t.
(14)
Now we are ready to obtain upper bounds for the four terms recapitulated here:
H (H (.Ss)) "11,,111,000 M.s:s[(H''. s.s.) (Hs.s.) '11l B,
+ +|[H (HS,.)(H.,,.) B,1,00
| Ses, Hss,, S.S.)
+ Hss (Hs,s,)
We derive the bounds separately.
For the first term, based on Assumption consider
Hsys, [H S.s.) "
Hsys, (Hs.s.) Hsss. HHH s. - 1,,,000
Hsys, (Hs.s,) ||,,, Hs.s. S.S ll 10,00 (H "
at)
Hs,s, H H.s. 1 v Sooovsss)) || S. 11112222
Taking
Amm(Hs and & 18[1-))|Ss1 Ann (H...,) in Equation (1) and adopting [.
24(1 xx)
[.
Equation
13
we conclude that, with probability at least 1 - 2|S exp
mo au (Hs.s.)
1152 1s
2|S,|² 2 exp
m (Hs. and & min [8]1-
Amia A (Hs.s.) [S+]
An
8[S.
16|8.
H Hs:s, (HH...)
"
6
For the second term, rewrite as
[[ s.s.
((s,, 1 11
Q Scs.
Hs. (Hs,,,,
Q
Hl.s.11lB...00/1 (Hs.s.)
S.s.
|
Hsss, Hss,s]llB,,oo Aminl Hs.S.)
scs.
20

---

Using Equation (10) by setting 121
/Nii S,) and & & 24pau [,
a
Amal (sss,))
we have, with
Srl
S.
probability at least 1 2/8 S,1 exp
mai Auia S...s..)
and & 24pmx
Auin Hs....
288 (S,
[.,
[ s,s.
(Hs,s)
a
6
For the third term, we obtain the upper bound
[ S;s. Hs:s,
[ s.s.) (Hs,s,) - Hsys, IIl.lllll s.S.) (Hs..,)) " 2 20
H C S,S.
Taking t in Equation 10
a 6
in Equation 14 Taking t
2 zV /e and 2m Sp & 3 \
& We establish The upper bound that, with probability at least S| exp ma(Am((Hs...)) (Hs.s.))
S, exp
48
ma
(5,T2)
1,2 Amia S.)
2|5,1 exp m(Amm (Hs, ))
8|S.1²
and
a (Hs, min 1pmn [.| 6
8|8.1
V 6|8.1'
16|S
[
Hs:s, I s.s,)
( &
For the fourth term, in accordance with Assumption 4
H Hsss, (Hs,s,)
& 11-6
In conclusion, we have shown that, with probability at least - 2| 2 exp
ma Am (Hs,
mo
1152(1 L-a 1s,
S| cxp (- mlAmm (Hamm (H..
cxp ( c 2|
Amm (Hs.s.)
-2S cxp 18 |8,
mi
8/8.
288pm mox [S,1"
2|S,|2 exp
maxAmii s.)) - -2,,, exp
m(Ania Hs. 11
and
192 S.1
8 8.S.1
& min(
a 48(1 a
Amin H H...s.)
Aminl (.S.s.)
a 24 Srl
Amin Hs,
S,|
16 1
SS
Amin (Hs.s.)
Amin Amin(Hs,... (Hs.s.)
4
S
8|S
6|S
16|5,1 b
the mutual incoherence condition holds for any worst-case distributions:
|| | Scs, (H 1 1,,,0 < 11- a
2
Simplifying Cyn
the above expressions, + log |S + log (S,1)) and
with C
probability at least
O(exp( Pan 1
Pres Pr.///////21
TH scS, (H..s))- 11 $1-2
2
where C' only depends on a, Aminl Hs.s.).
Lemma 12. If Assumptton
holds, then for any Q E A ((Pmm and a (0, 1], with probability at
least 1 - |S, pr exp
ma
20 C & and X;, 32 a/22 we have
a
B (X))s,
)s eT 2,00
8(1 a/2)
With probability at least 1 |Co, exp
mp² 20
E E 1 and X3x
/2jt/Paan Paeeee have
a²a
(x2E(X))sse)))
R,2,00
8
21

---

Proof. We start with [[ eT
2,,000 After some algeraic manipulation, we find that
[| eT
2,00 & max
(Xx
158,
max 158,
pp Pr max (eE
jEPr
max
Pr max
Bo E
158,
jep,
max v Pr max Rolesl- 188,
jepr
Since |eyl is bounded random variable according Assumption to get
apply Hoeffding's inequality
mt²
Prob( leyl > + t) & exp
20
Base on similar argument as Equation
9 we can derive
E E
Ep leyl & 280,
which leads to
mt²
Prob(Egles) > 2eo + 1 exp
2g²
Taking the union bound over all i S, and j j Pr. we find that
mt²
Prob [[E]B e1llo > (2r + p + 1)) & |r exp (-
2x
Setting
μ and while requiring X
32p/aal1-aa1
With probability at least 1
|S exp 2
my² we have
AB 8(1 - a/2)
12,00 &
(15)
Then we consider
1..2000
[[qE(X)) |,,,,
max ieCo
max PPPrr max Ex) jen; kE
ieCo,
mnax
Pipr max Eglex/-
kep.
Similarly, applying Hoeffding's inequality and the Knnbotovich-Rubinsteenn theorem gives us
Prob([EoE((Ap) 1,2,00 > Pr + p t)) < [Co,lp, exp (- mt?
20
Let t = & and
B
PPaPP hold, we have, with probability at least 1
|Co, exp (- mu²
20
A 8
[B2E(X)) 1,2,00 ||
Theorem 9. Given a Bayesian network (g. P) of n categorical random variables and its skeleton
Gstel (V, Esta)- Assume that the condition W || & B holds for some B 0 associated with an optimal Lagrange multiplier X for W* defined Equation 11. Suppose that W
DRO risk minimizer of Equation 4 with Wasserstein distance of order and ambiguity radius
22

---

cosm/m where m is the number of samples drawn i.i.d. from P. Under Assumptions 2 number of samples satisfies
3
if the
C( + log (n/o) log o2 P PmanPfin m O(
min(p²,1)
where C only depends on a, A, and if the Lagrange multiplier satisfies
32pmmaxx
3
A
a
(a/ 2a) + 2) Pmax v P(v)
then for any ó (0, 1], r [n], with probability at least - 6, the following properties hold:
(a) The optimal estimator W is unique.
(b) All the non-neighbor nodes are excluded: Co, = Cor.
(c) All the neighbor nodes are identified: Ne, Nee,
(d) The true skeleton is successfully reconstrucied: Gátet Gatet-
Proof. We prove the statements in this theorem in several steps. In order prove (a) and (b), will show that the DRO problem is strictly convex if true non-neighbors are known so that there an optimal solution. Next would like demonstrate that this solution with non-neighbor constraint is indeed unique for all the solutions without constraints. The proof for uniqueness comes with a conclusion that do not accidentally include any edge between the current node and its non-neighbors. Next, to prove (c), we present a generalization bound for the DRO estimator in terms of its true risk, which leads to a Ex bound of the difference between the estimator W and the true weight matrix W Combined with the assumption on the minimum weight, it implies that we include all the neighbor nodes successfully. Finally, by taking a union bound for all the nodes, we could conclude that the correct skeleton recovered with high probability, which proves (d).
(i) Given the true non-neighbors, there is a unique solution.
We start with the Wasserstein DRO problem, which recapitulate here for convenience:
1 Boll((X,) WT
W E arg inf W
sup Qe WP (?..) 2
The objective is convex because is supremum of convex functions.
For now, we assume that the non-neighbor nodes Co, are given. We can then explicitly restrict W. 0 for all i Co,. The Hessian of Ws,, is block diagonal matrix reads
H
S.S, 0
0
H S.S,
0
(2R.(Ws..)
PRP-P xpopser
:
0
H S.....
S.S.
where
H BylE(X,)8(xx))]] E Ree
is the covariance matrix of encodings of X₂ under some distribution Q A (PPm)
Since Wsy is fixed be zero and (2 (.s..) is a block diagonal matrix, we focus on showing that H > 0.
We apply Lemma to get the bound
Amin (H S. > Anin(H,s.) - 45|S,12 -t,
23

---

with probability at least 1 - -|, exp 2|8,1²) Amin (min......) 4c|S,| - />> 0 will guarantee
mf
that the DRO problem in Equation 4 has unique solution when the W, is satisfied for non-neighbor nodes.
(ii) Given the true non-neighbors, the solution is optimal.
We would like show that the solution Equation (4) with true non-neighbor constraints is optimal. In this way, we do not recover any non-neighbor nodes in the skeleton. We adopt the primal-dua
witness (PDW) [Wainwright 2009) method show optimality for the constrained unique solution.
Recall that we assume | Bl2,1 < B. To begin with, we write the dual problem as
We arg inf W
sup
BollE(X,)- w Zz
B)
2.2. al 2
(16)
s.t. ViE Co,
W = 0,
where Xy is the Lagrange multiplier for the norm constraint on W.
W is optimal if and only if there exists (Q*: Z*, X11)) that satisfies the KKT condition:
(x, T W Equ ((x,)
1= = 0
Q* A (Pm). z 1,2,00 1, 1 3 0, (z* W) W B. W|
W B,2, & B
B,2.1 B) 0.
Note that we assume that the constraint W B,2, & active such that XΦ > 0. This assumption only for convenience of theoretical analysis and restrictive. If it not active, we have which leads to an unconstrained problem similar to
[| Bl.,1 B < B for some B and 2 the ordinary least square problem, which is known to suffer from overfitting. Instead, we are usually interested in solutions that have finite norms so we can always find B c < for some small positive constant 0 to make the constraint active and thus X 0.
Substituting E(X,) WW*TE(X, + e into the first-order optimality condition yields
Ex,E(X)) W W*)
Ee E( + 1
H O* S.S. 5
Ws.
W
E (Xp)s,eT Bzg E (X (Xp)seeT
0 0
S,-
2 Z
(17)
H 0
Q+ s:S:
0
S.S.
Solving for z we find that
A Z X
H ( (X))(()))
SES.
scs,
which can be bounded such that
x 2 1B2.00
a s:s,
S.S.)
- -
H
scs. E + R |
B,2,
x S H S.S.
Zz ||,,, + H D*
(H ) E eT
B,2 [E
S.
11 25 2|,,,0000 + || (HS 11
1,00 (x,
S;S, S.
+
E(F))s |
Note that
Z 11.000 & 1 B,2, & 1.
24

---

Recall that 0 < a & 1 in Assumption 4 Based on Lemma 8 and Lemma 12
we may write
7 1,,,00
H Sps, (HS.s.) ) 1,
2g 11,,00 + || (.,,)) 111,,,,,,00
+ [ ( s |,2,00
a < )
a
1Ba
) ABa
+ (1
2
8(1 - ax
8
<18(11 a <x
7
with high probability and certain conditions on X and E.
Henceforth, - 8,2,00 satisfies strict dual feasibility and must have || s;. 11,211
according to complementary slackness:
Z** W) [ 1,2,1- In other words, have
Vi Co, W; O,
with high probability. This guarantees that we do not recover any node that is not a neighbor of with high probability.
(iil) Without information about the true skeleton, we have a unique and optimal solution.
We follow the proof of Lemma 11.2 in Hastie et al. [2015].
We have shown that W satisfying W, Vi E Co, is an optimal solution with optimal dual variables 5 1,2,00 1.
To avoid clutter of notations, define
(W)) == sup QEA (PP.))
1 Ee -TE(x)))
2
Let (W, be any other optimal solution infw sup, DRO (D) + X(| B,2, B). By defini- tion, DRO (W) X B,2.1 B) DRO W (z* W) B)
DRO (W) + À(I W B,2,1 - B) 2** W)
I, DRO ( + ( w W) - B)-
The first-order optimality condition for W says
VL DRO W + X h" 0,
which implies
X([ 8,2,1 B) (B ( W)) I DRO W + DRO (W), W - - 1-- DRO DR))
By definition, W ||,2,1 B and \ 0. Since T.DRO(.) is convex, the RHS of the above equation should be non-positive, or equivalently,
[ B22, & (Z* W)
On the other hand,
z W) < 1 0, W 1B2,1 & |
Therefore, the equality holds for the above inequalities, which leads
W W B,2,
- W)
Recall that || 1 < In order for W | Z W) hold, we must have
Ws= = 0.
25

---

In that wise, all the optimal solutions W have
W 0 Vi Co..
This implies that we have a unique solution that excludes all the non-neighbor nodes without information about the true skeleton. Until now, have proven properties (a) and (b).
(iv) The set of correct neighbors is recovered.
Consider again the first-order optimality condition in Equation [7].
Ss.- W = S.S, ) (Bq+ (xss x
Ws
W 8, 1B2,00 = -I( S.S.) (B (X+)s, 13 B,2,00
[1 S.S. |,,, R 8l X F)s.et
|| S.S.)
||,,1,00
e / 2,00 + | 11l,
< Pmaxx [H
S.S.
20,00
2 + \3)
Pmaax V S,
S.S.
2,2
B & (Xp)), e 1122,00 \
According to Equation 13 with probability at least
- expp
m(Ami HHs. 11² and & & 8|S.
(
16|8. 2
4 V Aminn (Hs,s,)
12222
(H
S.s.)
According to Equation /2p (1-a/2)
15 with probability at least
| exp (-
my² 20
) & and X >
we have
1Ba
(gE()))
2,00 & 8(1 - a/2)
On that account, with probability at least -2|S,12 exp m(Amin ( 8,.))
| exp ( my²
8|S
26
and & min (Am)(Hs, ()) while requiring XB 16[8,1
(1-0/(2))
W s... W w.. 11,2,0 & Pmnxx V | Amin H S.
ABl
+ 1).
8(1 a/2)
By Assumption 2 if the condition XØ
Annl (mn(H....)) is satisfied, the follow-
2( 8 3
+ 1) pm [.1]
ing inequality holds:
[Ws. W W. < 8/2.
In this way, we are able to recover all the neighbor nodes with threshold 8/2. This proves (c). (v) The true skeleton is recovered with high probability.
The above arguments tell us that with high probability and certain conditions for and A satisfied, for each node r, we do not recover any non-neighbor and do recover all the neighbor nodes. The correct Ne, and Co, are thus identified. Now we are ready to prove (d).
Putting everything together and taking the the union bound for all nodes r [m]. with prob- ability at least O(n exp
Cma² (- o2 Pin) + 2 log P(m))), E &
C#
and
32
X B
P
3/2
2 811 w/23 1oma V 4, where O only depends on a, A, we have
Gakel = Gskel-
26

---

Setting = and making the dependence on the sample size more explicit. We draw the conclusion that, if the number of samples satisfies
Clsa log (n/ō) log P
m O
min(p?,1)
where C' only depends on a, A, and if X satisfies
32/1pmax
3 AB < (a/ 2x) + 2) v [(n] 4
then with probability least 1 for (0. 1]:
Goucl Gskel-
Moreover, if we assume that the target graph has bounded degree of d, the sample complexity becomes logarithmic in n:
m = O C(so + log (n/8) + log n log Pmxx) o
Pinax d3
min(n².
Theorem 10. Suppose that W is a DRO risk minimizer of Equation (4) with the KL divergence and an ambiguity radius Eu /m. Given the same definitions of (g. P), Gshet. B, XB. in Theorem 9 Under Assumptions
Eu/m.
3 4 if the number of samples satisfies
Clen + log (n/0) + log P(n)) o Pmnstinn) m O
min(p?,1)
where C depends on a, A while independent of n, and if the Lagrange multiplier satisfies the same condition in Theorem 9 then for any 5 (0, 11, [n], with probability least 1 8, the properties (a)-(d) in Theorem hold.
Proof. Define
(x 1 -11((1))) WT
According to Theorem 7 in Lam 2019), the worst-case risk with KL divergence ambiguity set can be bounded as follows:
Σ (zw()) 22 Cs E 1e w 13 =E[m]
sup Ealw(X) <EE (X) + Ve, (P.
Eve(mng (e (2(()()) [ )
m
<B <Ep l (x() + Ve max ic(m)
(e lwl + Ce max is [m] (2(() -lwlp
where fw Lie(m] lw (2) and C' > 0 is constant independent of n.
27

---

Consider
max (w( - -lwl &
max w, W' 2,.
((((x) - (wr(
[m]
max I (x))
W
max mnax(1e((X,)l)) |
W.
max V/max + [l,,))
W. .m
max (BnP/x/Pmxx + [ W
max Max(V/mm (Vmmax /Pnnl W ||)
W. x
max Pmmax + VP[n] 1.2.1)
W.n
(v Pmaxx + B)²
= Bp
Define Emax wax(maaaeaaaa. e). Therefore, we find that
sup Eolw(X) & Ep l (X) + CEmax Bp.
Qe.A( (fm)
Similar to the Wasserstein robust risk, we observe that the following results hold for any Q AP(mm.)
With probability at least 1 - 2|5, exp 2/8, ), we have
mt²
Amin (H. > Amin Hs.ss)) - Cêmax |S, -6.
With probability at least 1 1 2|5, exp (- 2p 2p (8,12)
mt?
| S:S. Hsys, 118,111,,1,00 Ce [rr t.
With probability at least 1 - | exp (-21.
H S.S.
Hs.s, 10,000 C.max |S, t.
With probability at least 1 - 2|S, 2 exp
mtt (Amin (Hs... - 2/S,,12 exx (- m 8|S,1"
and
32(S.1"
Emax & C min
ANaii (H..s.))
Ami (HS.s.)). )
S|S vis.1
16|5.1
(HH,ss)
(Hs,s,) 11,,00 - & t.
With probability least 1 Olexp
Cm Pin P.
log |59 log ((()) and max Pass 3,7377
I S (H3.s.)
|l $ 11 2'
where C' only depends on a, Amin( Hs,ss).
Thanks to the boundedness of the error term e, we have similar conclusions to Lemma 12 if Emax holds.
In such wise, the properties in Theorem 9 hold with the same condition on A and the condition on
Emux that Emax
Cy 3/2 Since we set e and define Emax max( Vē, e), the condition on
(raes(()) (n)
Emax implies that
Eo 0 Pi [m]
3/2 co
max
12
#
28

---

Table 2: Comparisons of F1 scores for benchmark datasets and BIC for real-world datasets (backache, voting), BIC is not applicable to skeletons. The best and runner-up results are marked in bold. Significant differences marked by i (paired t-test, 0.05).
0LE0877
0.3143
014387
0.5553
04633
0.360
0.20m
0
0.3800
0.80885
09582
0103)
3.3300
01000
0.9933
0.0000
0L000
0.290
0.78001
0.75091
0.367
11:33
0100977
119907
1.973
0.0011
0.611
0.2133
14117
013677
11
1233
0.50081
0.5001
63075
198.0
0100
0.0000
06000
0.000
0.0933
0.880
0.000
(03)
0.72001
000754
031661
63800
0.111
0.2666
0.1653
0.0314
0.151
0.166
17
1130
Netectic
0.662)
01500
11500
0.3831
14801
64500
0.5035
0.4576
(40%)
0114100
6118.0
1.500
(4000)
018081
0.0010
barles
11:1
1115227
1.007)
14027
11347
0126
E331
0.808
1330
61531
11990
15
316
30A
304
304
104
304
-481390
-4655675
1435-3731
tecouce
Netectec
304
1179 134
-178-505
-1733330
-171975543
-175/2611
Longmabee
2486-1166
0314
2.990
Malee
304
.80006
-9472 (570
-1074012543
-10748.224
-107491.2543
-187 -187400.2543
The final sample complexity becomes
C(so + log (n/0) + log P Pomxx
m
O
min(u², 1)
D
More Empirical Results
Table 2 lists the complete experimental results.
29


Reviews:
[{"review_id": "n0z9PwqjhE", "reviewer": "Reviewer_Hyho", "summary": "The paper considers the problem of learning the exact skeleton of general discrete Bayesian networks from potentially corrupted data. The estimator optimizes the worst case risk over a family of distributions within bounded Wasserstein distance and KL divergence to the empirical distribution. Under mild conditions, the paper presents non-asymptotic guarantees for successful structure learning with logarithmic sample complexities for bounded degree graphs. Numerical studies on synthetic and real datasets validates the performance of the estimator.", "questions": "1) Is it possible to extend the DRO framework to learn the structure of continuous Bayesian networks from corrupted samples. A high level discussion regarding this would shed more light on the limitations of the present contributions.\n\n2) Since no assumptions on the underlying true distribution is made, it is only possible to recover the structure upto a Markov equivalence class. However the proposed approach overcomes this limitation and recovers the true skeleton with high probability. A discussion as to why DRO framework achieves this will be highly appreciated.", "limitations": " ", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": " ", "weaknesses": " ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the positive feedback. Please find our rebuttals below.\n\n**Q1: Is it possible to extend the DRO framework to learn the structure of continuous Bayesian networks from corrupted samples?**\n\nThe proposed DRO framework is based on a regression approach, which is naturally tailored to tackling continuous variables. For example, by assuming a linear Gaussian structural equation model, we can simply let the encoding function be the identity mapping. We have added the discussion on this point in our revised manuscript.\n\n**Q2: A discussion as to why DRO framework achieves this (recovers the true skeleton with high probability without distributional assumptions) will be highly appreciated.**\n\nOur methods are based on a series of work on DAG structure learning in polynomial time. Since the non-robust method our work is built upon is a special case of our methods, the desirable properties are naturally inherited. Moreover, the proposed DRO framework enjoys robustness that is modeled explicitly and converges to the true global optima as $m \\rightarrow \\infty$ given clean data samples as data distribution uncertainty vanishes. The Discussion section has been revised to address this issue."}, {"review_id": "zaMzAVjSMq", "reviewer": "Reviewer_AT8n", "summary": "The paper presents a method for learning the structure of discrete Bayesian networks from potentially corrupted data. It utilizes distributionally robust optimization and regression, optimizing the worst-case risk over a set of distributions. The approach works with general categorical variables without assuming specific conditions and provides efficient algorithms. Non-asymptotic guarantees for successful structure learning are derived, and numerical experiments validate its effectiveness.", "questions": "l 105: How is T defined in A_{ST}?", "limitations": "- The proposed method accounts for the effect of outliers by optimizing the worst-case risk over a family of distributions within bounded Wasserstein distance or KL divergence. However, the extent to which the method can effectively handle outliers and their impact on the learned network structure should be further explored and evaluated.\n- The proposed approach is designed for general categorical random variables without assuming faithfulness, an ordinal relationship, or a specific form of conditional distribution. While this makes the approach flexible, it may also limit its applicability to certain domains or types of data. Investigating the generalizability of the approach to a broader range of variable types or distributions could be a potential limitation.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Originality:\nThe paper addresses the challenge of missing data in distributionally robust skeleton learning. This aspect adds a unique and valuable contribution to the field.\n\nQuality:\nThe quality of writing in the paper is commendable. The ideas and concepts are presented effectively, and the paper demonstrates a strong grasp of the subject matter.\n\nClarity:\nThe paper's clarity is notable, as it is organized in a logical manner where each section builds upon the previous ones. The flow of information is well-structured, making it easy for readers to follow the paper's content.\n\nSignificance:\nThe task of building discrete Bayesian networks from potentially corrupted data holds significant importance, as it has wide-ranging applications in various fields. By addressing the challenge of distributionally robust skeleton learning, the paper contributes to the advancement of techniques that can handle data uncertainties and enable more accurate modeling and inference in real-world scenarios. The significance of this research lies in its potential to enhance decision-making processes, predictive modeling, and knowledge discovery across diverse domains.", "weaknesses": "- The author's assertion that they had limited space to relate their work to the existing literature is understandable. However, it is crucial to provide a clear positioning of their own work within the broader research landscape. Emphasizing the novelty and distinguishing features of their approach will strengthen the paper's contribution and demonstrate its unique value.\n\n- The paper lacks clarity on whether cross-validation or hyperparameter tuning was performed during the experiments. It would be beneficial to address this issue explicitly to ensure the robustness of the results. Additionally, to further improve the reliability of the findings, it is recommended to conduct a rigorous statistical test, such as the corrected Student's t-test proposed by Nadeau and Bengio [1], using 15 hold-out folds. Incorporating such statistical analyses will enhance the credibility and validity of the experimental results.\n\n[1] https://link.springer.com/article/10.1023/A:1024068626366", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your positive comments. Please see our point-by-point responses as follows.\n\n**Q1: The author's assertion that they had limited space to relate their work to the existing literature is understandable. However, it is crucial to provide a clear positioning of their own work within the broader research landscape.**\n\nThanks for the suggestion. We have added more related work to put our work in a more clear posistioning within the research landscape by assuming that extra space is available.\n\n**Q2: The paper lacks clarity on whether cross-validation or hyperparameter tuning was performed during the experiments. ... Additionally, to further improve the reliability of the findings, it is recommended to conduct a rigorous statistical test...**\n\nWe tune the hyper-parameters based on randomly generated Bayesian network with the same number of nodes as indicated by the input data. We will include this information as well as statistical tests in our revision.\n\n**Q3: l 105: How is T defined in A_{ST}?**\n\n$\\mathcal{T}$ is a subset of $\\\\{1, \\dots, m\\\\}$. Revised.\n\n**Q4: The proposed method accounts for the effect of outliers by optimizing the worst-case risk over a family of distributions within bounded Wasserstein distance or KL divergence. However, the extent to which the method can effectively handle outliers and their impact on the learned network structure should be further explored and evaluated.**\n\nThe class of distributionally robust optimization methods account for outliers in terms of uniform performance over a set of distributions and the outlier is assumed to be identifiable through divergence. The way of robustly modeling is general but lacking case-by-case analysis for a specific or given type of corruption. We aim to demonstrate a robust learning method that subsumes existing methods as special cases while enjoying certain robustness. There is also an important issue on robust identifiability of Bayesian network or causal structures. These are all interesting future directions to study further. Discussion on this point has been added in the revision.\n\n**Q5: The proposed approach is designed for general categorical random variables without assuming faithfulness, an ordinal relationship, or a specific form of conditional distribution. While this makes the approach flexible, it may also limit its applicability to certain domains or types of data. Investigating the generalizability of the approach to a broader range of variable types or distributions could be a potential limitation.**\n\nGood point. We agree that we sacrifice some efficiency and applicability to special domains for general formulations and bounds. For example, if the possible values of the random variables are ordinal or continuous, or the variables follow a Gaussian or Poisson distribution, a more efficient algorithm and a tighter bound can be established. We have added discussion on this point at the end of the manuscript in revision."}, {"review_id": "I6TvxTec4O", "reviewer": "Reviewer_YEGV", "summary": "In this paper, the authors revisit the problem of structure learning of Bayesian networks. Here, we are interested to learn the underlying DAG structure that encodes the conditional independence properties of a Bayes net from random samples. A number of algorithms have been proposed to solve this problem over the years which can be broadly classified into two groups: score-based and constraint-based. Despite a lot of research activity over the last several decades the general problem is not yet known to be solvable in polynomial time and several approaches are shown to be NP-hard. \n\nThe authors of this paper start with an approach first proposed by Bank and Honorio in 2020: to reduce this discrete structure learning problem to a continuous regression problem based on encoding schemes and surrogate parameters. The idea is that once cast as a regression problem, techniques from structure learning in continuous models such as regression can be employed. The main contribution of this paper is a distributionally robust optimization (DRO) method that builds upon the aforesaid framework to handle noisy data. In this proposed method, the best structure is searched among a ball of certain radius around the empirical distribution that minimizes the regression cost. Specifically, balls according to the Wasserstein distance and KL divergence are considered. The authors give two theorems for each of the above two distance functions under certain assumptions. They further verify the effectiveness of their methods by experimenting with benchmark real and synthetic datasets. \n\n ", "questions": "None", "limitations": "None", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The problem considered in extremely important from a theoretical and practical standpoint. The approach followed by the authors to reduce it to a regression problem as in the previous work could be promising for certain datasets.", "weaknesses": "A major weakness of the paper is that a number of assumptions are made in order to derive the results, which looks ad-hoc and too restrictive in my opinion. See page 4 of the paper for these assumptions. See also line number 271. \n\nI could not help but notice several vague statements throughout the paper, a few of which are listed below. \n\nLine number 224-233 the NP-hardness argument is very vague. Looks like you are using a brute-force approach, if so, the algorithms wont scale beyond 15 nodes. \n\nLine number 240-241: the citation of Chen and Paschalidis [2018] here is for the first time and sudden. What question are they solving? How is it related to your problem? \n\nUnfortunately, I am not convinced from the experimental results that the proposed algorithms are overall performing better than the existing methods. Many of the entries in table 1 show that other methods are performing superior to the proposed algorithms. \n\nFinally, I found the paper very hard to follow, maybe it reads more complicated than it actually is. The author may try to see whether writing it in a simpler manner with small examples makes the paper more readable or not. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your insightful comments.\n\n**Q1: A major weakness of the paper is that a number of assumptions are made in order to derive the results, which looks ad-hoc and too restrictive in my opinion. See page 4 of the paper for these assumptions. See also line number 271.**\n\nThanks for the good question. Assumption 1 simply finds the upper bounds of random errors but does not require them to be small. Assumption 2 makes sense such that if a neighbor node has weight close to zero, it will become indistinguishable from those having little impact on the current node. Our intuitive assumption is that neighbor nodes should have a larger impact than non-neighbor nodes do. Assumption 3 is made to make sure that the solution is unique. However, the original problem is convex without this assumption. Therefore, this assumption only serves as a tool for theoretical analysis but does not lessen its practicability. Assumption 4 has been shown to hold in practice by Bank and Honorio [2020]. We guess you mean the inequality right after line 271. This is not an assumption but the derived condition to be satisfied.\n\n**Q2: I could not help but notice several vague statements throughout the paper, a few of which are listed below. Line number 224-233 the NP-hardness argument is very vague. Looks like you are using a brute-force approach, if so, the algorithms wont scale beyond 15 nodes.**\n\nLine 224-233 argues that the supremum problem is a NP-hard problem where we omit the proof of reducing MAXQP to our problem. We are **not** using a brute-force approach but a random greedy algorithm whose execution time is $\\mathcal{O}(n^2m\\rho_{max}T)$ where $T$ is the specified number of iterations. The global optimal solution is obtained as $T \\rightarrow \\infty$. However, data in practice often elicits fewer than $10$ iterations to the find the solution. Given you question about it, we have revised this paragraph and added NP-hardness proof.\n\n**Q3: Line number 240-241: the citation of Chen and Paschalidis [2018] here is for the first time and sudden. What question are they solving? How is it related to your problem?**\n\nThanks for the good point. Chen and Paschalidis [2018] solved a distributionally robust regression problem with continuous random variables and showed that the distributionally robust optimization problem is equivalent to a regularized ERM problem under some conditions. However, this does not hold in our case. We have added clarifications in the revision.\n\n\n**Q4: Unfortunately, I am not convinced from the experimental results that the proposed algorithms are overall performing better than the existing methods. Many of the entries in table 1 show that other methods are performing superior to the proposed algorithms.**\n\nThe results in Table 1 show that our methods (Wass, KL) or the resulting DAG with the aid of HC achieve either the best or the runner-up performance. We will add experimental results on more datasets and make the advantages more clear in the revision.\n\n\n**Q5: Finally, I found the paper very hard to follow, maybe it reads more complicated than it actually is. The author may try to see whether writing it in a simpler manner with small examples makes the paper more readable or not.**\n\nThanks for the suggestion. We have added small examples, made more concrete explanations and reduced mathematical details in our revision.\n"}, {"review_id": "y1WvaWLLzK", "reviewer": "Reviewer_onGN", "summary": "The authors propose methods for learning Bayesian network skeletons based on discrete data when the data is corrupted by arbitrary noise.\n\nThe authors show guarantees on discovering the correct skeleton in terms of the required sample size. Some empirical experiments are also presented.", "questions": "To be honest, the level of mathematical detail in the paper is too much for me. It's not that I wouldn't be able to refresh my measure theory and calculus but I just can't afford to spend a week on a review. This is of course primarily my problem, but I have reviewed and published at NeurIPS before, and based on my experience, this paper is mathematically too dense for at least 90% of the NeurIPS community. But perhaps this is fine in case the 10% find it relevant and useful -- I just couldn't possibly tell.\n\nIn Sec. 4, explain what exactly you mean my the F1 score: the F1 score is used for binary classification, so I'm guessing you're perhaps looking at presence/absence of each of the edges in the DAG. Please clarify further which methods are proposed in this paper (Wass + KL?) and which are benchmarks (Reg, MMPC, GRASP?) so that it's easier to quickly grasp. \n\nAm I right in that the punchline is Figs 1.b and 1.c? If so, please make it easier for the reader to see this or otherwise they'll miss it. \n\nI take it the results indicate that the proposed methods are otherwise about as good as prior work except that they're better when the failure probability is close to 0.5 (Figs. 1.b and 1.c)? If so, doesn't this require an extremely high rate of failures, to the extent that the data is mostly random noise?\n\ndetailed comments:\n- p. 1 \"DAG structure of a Bayesian network is usually unknown\": citation needed (or maybe just say \"often\")\n- p. 2 \"contains as subgraphs the skeleton\": I'm confused by the plural/singular tenses here: if there's only one skeleton, how is it contained as subgraphs\n- p. 2: \"(fixed-parameter) tractable\": what does this mean?\n- p. 2: \"observational data is commonly contaminated [...]\": again, citation needed\n- p. 2: \"may produce false edges\": may it also miss edges?\n- p. 2: \"do not assume any specific form of conditional distributions\" vs \"set of distributions characterized by certain a priori properties\": aren't these statements contradictory? please elaborate\n- p. 2: \"small amount of samples (high-dimensional) and potential perturbations\": the word 'high-dimensional' is grammatically speaking placed in an odd spot; also, does \"small amount \" refer only to samples or also to perturbations -- please clarify the sentence", "limitations": "ok", "rating": 5, "confidence": 1, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "+ rigorous mathematical derivations", "weaknesses": "- very dense mathematical derivations\n- limited and inconclusive empirical results", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your very nice and helpful comments.\n\n**Q1: In Sec. 4, explain what exactly you mean my the F1 score: the F1 score is used for binary classification, so I'm guessing you're perhaps looking at presence/absence of each of the edges in the DAG. Please clarify further which methods are proposed in this paper (Wass + KL?) and which are benchmarks (Reg, MMPC, GRASP?) so that it's easier to quickly grasp.**\n\nThe F1 score can be leveraged as an evaluation criterion in structured prediction. Specifically, the label of each edge (zero for absence, one for presence) can be considered a binary random variable and the total number of predictions to be evaluated is the total number of edges. The methods proposed and evaluated in the experiments are Wass and KL corresponding to the two DRO methods with Wasserstein and KL ambiguity sets respectively. The baselines are Reg (Bank and Honorio [2020]), MMPC and GRASP. While the HC algorithm is adopted to orient the undirected edges in the predicted skeleton to obtain a valid DAG. We have elaborated on these points in our revision.\n\n**Q2: Am I right in that the punchline is Figs 1.b and 1.c? If so, please make it easier for the reader to see this or otherwise they'll miss it.**\n\nYes, you are right about this point. Figure 1.b and 1.c demonstrate our superior performance in challenging settings in addition to the comparable performance in noisefree settings in Table 1. We have clarified this punchline in our revision.\n\n**Q3: I take it the results indicate that the proposed methods are otherwise about as good as prior work except that they're better when the failure probability is close to 0.5 (Figs. 1.b and 1.c)? If so, doesn't this require an extremely high rate of failures, to the extent that the data is mostly random noise?**\n\nA failure probability close to 0.5 is a challenging setting that showcases the advantage of DRO methods. Due to space limitation, we did not illustrate results on other datasets with varying failure rate. Usually, a noise probability of about 0.2 is sufficient to distinguish our methods from non-robust methods. It is noteworthy that for constraint-based methods like MMPC, the selected contamination model and distribution is likely to introduce dependency if a variable is binary.\n\n**Q4: detailed comments: p. 1 \"DAG structure of a Bayesian network is usually unknown\": citation needed (or maybe just say \"often\")**\n\nThanks for pointing it out. We have added citations for this statement. In fact, there are often confounders in real world. In fact, we can only use data likelihood for evaluation on real-world datasets.\n\n**Q5: p. 2 \"contains as subgraphs the skeleton\": I'm confused by the plural/singular tenses here: if there's only one skeleton, how is it contained as subgraphs**\n\nWe agree that the usage here may be confusing. What we really mean is that the super-structure is a superset of the skeleton, or the skeleton is a subgraph of the super-structure. We have revised the statement accordingly in our revision.\n\n**Q6: p. 2: \"(fixed-parameter) tractable\": what does this mean?**\n\nFixed-parameter tractability is the parametrized analogue of polynomial time tractability, e.g., $f(k) n^{\\mathcal{O}(1)}$ for some computation function $f$. We have added a citation in the complexity literature for completeness.\n\n**Q7: p. 2: \"observational data is commonly contaminated [...]\": again, citation needed**\n\nCitations added.\n\n**Q8: p. 2: \"may produce false edges\": may it also miss edges?**\n\nYes, false edges, missing edges, or anything as long as the prediction label (present, absent) does not match the ground truth label. We have revised the sentence accordingly.\n\n**Q9: p. 2: \"do not assume any specific form of conditional distributions\" vs \"set of distributions characterized by certain a priori properties\": aren't these statements contradictory? please elaborate**\n\nThe former refers to the Bayesian network distribution of interest, e.g., a log-linear distribution, whereas the latter means that the ambiguity set is constructed based on certain a priori properties about our belief on how the uncertainty looks like.\n\n**Q10: p. 2: \"small amount of samples (high-dimensional) and potential perturbations\": the word 'high-dimensional' is grammatically speaking placed in an odd spot; also, does \"small amount \" refer only to samples or also to perturbations -- please clarify the sentence**\n\nThank you for pointing out the grammatical issue. Small amount only refers to samples and high-dimensional means $m < n$ (number of samples is less than number of features). We have revised this part in our revision.\n"}, {"review_id": "5L3klwQcV7", "reviewer": "Reviewer_Np2S", "summary": "The authors give novel algorithms to learn the structure of a Bayes net from noisy data.", "questions": "Line 21:\nWhat do you mean by \"ordinal?\"\n\nLine 69 to Line 70:\nUnclear.\n\nLine 127:\nCan you please elaborate on faithfulness?\n\nLine 135:\nPlease define \"unweighted effects.\"\n\nCan you please explain Equation (2)?\n\nCan you please justify Assumption 2?\n\nCan you please explain Equation (5)?\n\nProof of Theorem 9:\nI do not understand the union bound argument.\n\nLine 304:\nWhy is it always possible to find such a solution?\n\nLine 313 to Line 314:\nUnclear.\n\nFigure 1:\nCan you please add some short *further* explanation in the caption?\n\nExperiments look OK!", "limitations": "Yes.", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "Originality:\nWork seems original.\n\nQuality:\nQuality is high, writing-wise and content-wise; see questions.\n\nClarity:\nWriting is clear; see questions.\n\nSignificance:\nTopic is significant.", "weaknesses": "No significant issues detected; see questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your positive comments. Our rebuttals can be found below.\n\n**Q1: Line 21: What do you mean by \"ordinal?\"**\n\nIf there is an ordinal relationship among variables, the relation is important to learning and inference. For example, if we are to rate the performance of an employee out of very poor, poor, neutral, good, very good, and the employee is actually doing a very good job, assigning a rating of very poor to them is less desirable than a rating of neutral. In contrast, if discrete random variables are not ordinal, different values have no additional interpretation but difference.\n\n**Q2: Line 69 to Line 70: Unclear.**\n\nAn ambiguity set is a set distributions that reflects our uncertainty about the underly true distribution. It can be defined in various ways, for example, we can simply use all the distributions, distributions with finite moments or distributions whose support set is included in a box constraint. A typical and intuitive way is to make use of the available samples to form an empirical distribution and considers distributions close to it. The closeness can be measured by a divergence between two distributions.\n\n**Q3: Line 127: Can you please elaborate on faithfulness?**\n\nThe faithfulness assumption is a commonly adopted assupmtion for constraint-based causal structure learning methods. A distribution is faithful to a DAG if all (conditional) independences that hold true in the distribution are entailed by the DAG, which is commonly violated in practice.\n\n**Q4: Line 135: Please define \"unweighted effects.\"**\n\nIn the unweighted effect coding, the column of the coding matrix sums to zero, whereas the weighted effect coding considers weights defined by samples. For example, given four categories, the set of encoding vectors could be $\\\\{(1, 0, 0), (0, 1, 0), (0, 0, 1), (-1, -1, -1)\\\\}$.\n\n**Q5: Can you please explain Equation (2)?**\n\nEquation 2 optimizes the empirical version of Equation 1 with available data samples. The constraint of Equation 1 is abandoned since the neighborhood and graph information are not available in practice. Regularization is a widely adopted technique to combat overfitting in machine learning, which is also the one of the limitations in existing methods that we attempt to address in the paper.\n\n**Q6: Can you please justify Assumption 2?**\n\nThanks for pointing it out. Assumption 2 is made to guarantee that the neighbor nodes have a noteworthy impact on the node of interest such that the impact can be quantified as a non-zero value bounded away from zero. In this way, neighbor nodes can be identified and distinguished from nodes with an impact close to \"zero\". We have added this justification in the revised manuscript.\n\n**Q7: Can you please explain Equation (5)?**\n\nEquation 5 is a corollary of Theorem 1 in [1]. For the special case in our work, briefly speaking, the proof writes out the Largrangian dual problem and notes that the supremum over distributions can be simplified into a finite-dimensional problem based on each point of the empirical distribution. We have added the citation to [1] here.\n\n**Q8: Proof of Theorem 9: I do not understand the union bound argument.**\n\nBefore taking the union bound, we have attempted to show that for every node, there is a probabilistic bound to recover all of its true neighbor nodes. However, the events of neighborhood recovery for different nodes may or may not be independent. A union bound gives a strict condition for the event that all the sub-events for each node's neighborhood recovery hold true simultaneously, though the upper bound may be far from a tight bound for special cases. Taking the union bound thus allows us to combine the bounds for individual nodes and reach the conclusion that the whole skeleton is consistently recovered.\n\n**Q9: Line 304: Why is it always possible to find such a solution?**\n\nThanks for pointing it out. First, our objective is non-negative and thus always greater than negative infinity ($>-\\infty$). Second, by assuming that the loss function is finite for at least one $\\mathbf{W}$, we have a non-trivial optimization problem, otherwise the loss is always equal to $+\\infty$ thus trivial. For example, $\\mathbf{W} = 0$ yields a finite loss for any finite encodings. Third, there is a typo such that $\\bar{B}$ does not have to be finite in the unconstrained case with $\\lambda_B^* = 0$. Therefore, as long as it is not a trivial optimization problem, we can always find a solution incurring a finite loss value.\n\nWe have revised this paragraph in our revision to make the explanation more clear.\n\n**Q10: Line 313 to Line 314: Unclear.**\n\nPrecisely speaking, we refer to a continuous distribution as an absolutely continuous distribution. It means that a continuous distribution will **never** be accounted for by a KL ambiguity set defined by a discrete distribution.\n\nFor example, given that the nominal distribution is a univariate discrete distribution with $\\mathbb{P}(X = 1) = 0.2$, $\\mathbb{P}(X = 2) = 0.3$, $\\mathbb{P}(X = 3) = 0.5$, then a uniform distribution over $X \\in [1, 3]$ will never fall within the KL ambiguity set because any distribution in the set will have, for example, $\\text{Prob}(1.1 < X < 1.9) = 0$ (since only $X = 1, 2, 3$ can have non-zero mass) whereas the absolutely continuous uniform distribution induces $\\text{Prob}(1.1 < X < 1.9) = 0.4 > 0$.\n\nWe have rephrased the related sentences in our revision.\n\n**Q11: Figure 1: Can you please add some short further explanation in the caption?**\n\nThanks for your suggestion. The caption has been revised.\n\t\t\n\n[1] Jose Blanchet and Karthyek Murthy. Quantifying distributional model risk via optimal transport. Mathematics of Operations Research, 44(2):565–600, 2019.\n"}]

Author rebuttal:
Dear reviewers,

Thanks for your helpful reviews that would certainly improve our manuscript. Please see our rebuttals after individual reviews. We are more than glad to have further discussions or make clarifications in the discussion period.

We would like to restate that our main contribution is the first robust Bayesian network skeleton learning method with polynomial time and sample complexities without distributional assumptions.

Best,
Authors


Meta-review: This paper considers the problem of learning the structure of Discrete Bayesian Networks from noisy data. The review process has revealed that the paper is theoretically sound, well written, and makes an important contribution of interest to the NeurIPS audience.
Decision: Accept (spotlight)